Phylogenetic analysis of maximal oxygen consumption during exercise (V̇_O2,max) and ecological correlates among lizard species

Maximal aerobic capacity, measured as the maximum rate of oxygen consumption (V̇_O2,max) during forced exercise (Seeherman et al., 1981), sets the upper limit to the rate of work that can be sustained (Jones and Lindstedt, 1993; Levine, 2008; Seeherman et al., 1981; Spurway et al., 2012). In lizards, several studies have shown that animals exercised below the maximal aerobic speed can sustain activity for long periods (e.g. 30 min or more), whereas activity above this speed causes a rapid decline in time to exhaustion (Bennett and John-Alder, 1984; Huey et al., 1984; John-Alder and Bennett, 1981; John-Alder et al., 1983, 1986; see also Beck et al., 1995; Perry et al., 2004). As V̇_O2,max directly affects the maximum sustainable speed and hence endurance capacity in relation to running speed, it has the potential to limit activity levels during ecologically relevant tasks, such as patrolling a territory, fighting, courting, foraging or escaping from predators (Biro et al., 2018; Garland, 1993; Husak and Fox, 2008; Sinervo et al., 2000). However, relatively few studies have tested for evolutionary co-adaptation between V̇_O2,max and aspects of behavior, ecology or life history that intuitively might depend on high aerobic capacity (Albuquerque et al., 2015; Autumn et al., 1999; Bennett et al., 1984; Boratyński, 2020; Clemente et al., 2009; Killen et al., 2016).

Bennett et al. (1984) found that the widely foraging Eremias lugubris has higher V̇_O2,max and activity levels in the wild, as compared with the closely related Eremias lineoocellata, but two-species comparisons cannot provide strong evidence regarding adaptation or co-adaptation (Garland and Adolph, 1994). Autumn et al. (1999) compared eight species of nocturnal geckos with 24 species of diurnal lizards. Nocturnal geckos are active at lower body temperatures, and when measured at their field-relevant body temperatures, they had a lower V̇_O2,max than diurnal lizard species. However, assuming a Q₁₀ of 2.5 (from previous studies of lizards in general; Bennett, 1982; Bennett, 1983), correcting to the same body temperature yielded similar values for nocturnal and diurnal species. Finally, in a study of 18 varanid species, Clemente et al. (2009) showed that V̇_O2,max was positively correlated with endurance, measured as maximum distance and time running when chased around a circular racetrack until exhaustion. Endurance was higher in widely foraging species and in species from xeric as compared with tropical areas, but these comparisons were not reported for V̇_O2,max.

Variation in V̇_O2,max has also been studied in other taxa. Among mammals, Albuquerque et al. (2015) found that V̇_O2,max was marginally positively correlated with home range area. All else being equal, larger home ranges should be associated with higher activity levels. However, the correlation between home range and V̇_O2,max was weak (r=0.278) and statistically non-significant in phylogenetically based analyses. The authors argue that measurement error (see Ives et al., 2007), including a mismatch between the populations for which home range and V̇_O2,max were measured, and the fact that activity levels should be higher inside the defended portion of the home range (territory), might account for the weak correlation between home range and V̇_O2,max. Killen et al. (2016) showed that size-corrected V̇_O2,max is higher for pelagic teleost fishes and for species from higher trophic levels.

In the present study, we tested whether V̇_O2,max in lizards is associated with diet (herbivore, insectivore, insectivore/carnivore, carnivore and omnivore; which presumably relates to foraging mode), climate (tropical, temperate and arid), nocturnality and viviparity, while controlling for the effect of body mass and temperature using conventional ordinary least squares (OLS) and phylogenetic regressions to account for species relatedness. We included viviparity because retaining embryos for longer gestation periods should have important impacts on activity levels and perhaps the intensity of effort during daily tasks. Pregnancy can reduce locomotor capacity in lizards not only through the effects of added body mass but also through unidentified changes in physiology (Garland and Else, 1987; Olsson et al., 2000; Zani et al., 2008). Therefore, all else being equal (e.g. in the absence of compensatory behavior: Bauwens and Thoen, 1981; Brodie, 1989), a higher V̇_O2,max may be required to maintain similar activity levels in viviparous versus oviparous species. Finally, we also present allometric equations for V̇_O2,max.

We obtained V̇_O2,max data from scientific publications. We started with the dataset provided by Garland and Albuquerque (2017) and added data for species outside of the 35–40°C range used by them. We searched for ‘lizard V̇_O2,max’, ‘lizard aerobic capacity’ and ‘lizard maximal metabolic rate’ in Google Scholar and Web of Science. We considered only studies that measured oxygen consumption while gradually increasing speed on a treadmill. When the studies reported mass-specific V̇_O2,max, we multiplied the reported mass-specific V̇_O2,max by the average body mass.

To increase the number of data points available for analysis, we measured V̇_O2,max for 18 adult male Sceloporus occidentalis Baird and Girard 1852, from Riverside, CA, USA (California Department of Fish and Wildlife permit SC – 013534). All procedures were approved by the UCR Institutional Animal Care and Use Committee. We obtained each lizard's V̇_O2,max while it was running on a speed-controlled treadmill (e.g. Garland, 1984; Garland and Else, 1987). Prior to each trial, animals were maintained for at least 2 h in an environmental chamber at 35°C. We placed each lizard at the rear of the treadmill belt between adjustable Plexiglas walls. Lizards wore a light-weight translucent plastic mask through which room air was pulled and conducted along tubes to an oxygen analyzer (S-3A, Applied Electrochemistry, Sunnyvale, CA, USA). Flow rate through the mask (500 standard cubic centimeters per minute, controlled by a mass-flow controller) was sufficient to capture all exhaled gas. A set of heat lamps above the treadmill was adjusted such that body temperature remained near 35°C during the trial. We stimulated lizards to run by gentle taps on their tail and hindlimbs as we slowly increased the treadmill speed while monitoring oxygen concentration from the exhaled air. When we observed no further decline in oxygen concentration despite increases in speed for at least 60 s, we stopped the trial. V̇_O2,max was taken as the highest 1 min of O₂ consumption during the trial (e.g. Garland, 1984; Garland and Else, 1987). We ran each lizard twice, once on each of two consecutive days, and used the higher trial as the V̇_O2,max. A paired t-test indicated no significant difference between V̇_O2,max measured on the two trials (N=18, P=0.26). The two measures were also highly correlated (r=0.82, P<0.0001).

A small tube with calcium sulfate, soda lime and calcium sulfate in series removed carbon dioxide and water vapor from the air before it entered the oxygen analyzer. A data acquisition module (ADAM-4019) converted the analog signal from the oxygen analyzer into a digital signal sent to a Macintosh computer. We used LabHelper software (WarthogSystems, www.warthog.ucr.edu) to record and visualize the data in real time. We used Warthog LabAnalyst to calculate oxygen consumption from records of oxygen concentration using the Mask Mode 1 conversion.

We gathered literature data on climate, diet, nocturnality and viviparity for each species for which we had V̇_O2,max data. Climate was defined according to the first letter of the Köppen–Geiger classification: tropical (A), temperate (C) and arid (B). We used a Google Earth layer based on the updated Köppen–Geiger climate classification provided by Peel et al. (2007) to determine the climate category for the locality where each species was collected. If a species was collected in multiple localities and separate measures of V̇_O2,max were not reported, we used the climate category that encompassed the majority of collection sites for the individuals collected (that was the case only for Varanus gouldii). When captive species were used or when collection sites were not described in the original papers, we used the climate category that encompassed most of the distribution of the species as reported on the IUCN website (https://www.iucnredlist.org/).

Categorizing an animal's diet can be complicated (Peters, 1977). We classified species' diet as herbivore, insectivore, insectivore/carnivore, carnivore or omnivore following the 90% rule (Harestad and Bunnell, 1979; Perry and Garland, 2002; Peters, 1977) when studies that reported stomach contents in detail were available. For those studies, if under 10% of the volume of stomach contents reported was composed of plant material, then we considered the lizard to be primarily insectivorous. If more than 10% of the stomach contents volume was composed of plant material, the species was considered omnivorous. Stomach contents of exclusively carnivorous and herbivorous species might not be as well preserved as the contents of species that eat mostly arthropods, so we also considered descriptions in general ecology papers or field guides, which usually involve direct observations in the wild. Finally, for one species, Teratoscincus przewalskii, none of the sources of information described above were available, so we classified the species as insectivorous based on descriptions of what caretakers fed the animals in captivity, as taken from descriptions on pet websites or forums (e.g. Good Life Herps, http://goodlifeherps.weebly.com). We classified species as viviparous or oviparous and as diurnal or nocturnal according to the information provided in the Reptile Database (http://www.reptile-database.org, accessed in 2019) and Animal Diversity (https://animaldiversity.org/, accessed in 2019), or from data available in Mesquita et al. (2015). The final dataset used for analysis is available from figshare (https://doi.org/10.6084/m9.figshare.13150487.v2).

We used the time-calibrated phylogeny provided as online supplemental material in Tonini et al. (2016) and removed species not present in our study. We then added seven new branches to represent the two populations that we had for Sceloporus occidentalis, Tiliqua rugosa, Varanus gilleni, Varanus gouldii, Varanus panoptes, Varanus rosenbergi and Varanus tristis for which we found V̇_O2,max data. We arbitrarily set the branch length between the two populations of these seven species as equal to half the shortest branch length between any two species in our tree. The shortest bifurcation in the phylogeny for our species was between Varanus panoptes and Varanus gouldii (5.47 million years), so we set the seven bifurcations to 2.74 million years (Fig. 1). This tree is available from figshare (https://doi.org/10.6084/m9.figshare.13150487.v2).

We computed multiple regressions in two ways (reviewed in Garland et al., 2005; Lavin et al., 2008) using the Matlab REGRESSIONv2.m program (Lavin et al., 2008). First, we used conventional, non-phylogenetic, OLS and then regressions in which the residuals are modeled as having evolved via an Ornstein–Uhlenbeck process (RegOU), which is intended to mimic stabilizing selection on the specified phylogenetic tree. The OLS models assume a ‘star phylogeny’ with no hierarchical structure, whereas the RegOU models alter the branch lengths of the tree to increase or decrease its hierarchy (respectively pulling nodes towards the tips of the tree or closer to the root). As compared with an OLS regression, the RegOU regression model contains one additional parameter, d, that estimates the transformation of the phylogenetic tree (Blomberg et al., 2003; Lavin et al., 2008) to better fit the residuals obtained with a given set of independent variables. A d-value greater than 1 means that the analysis altered the tree to be more hierarchical; conversely, values less than 1 mean that the nodes were pulled closer to the root, as compared with the original tree. The minimum value for this parameter is zero, which would indicate that the tree was altered to a star phylogeny.

Our dataset included V̇_O2,max for 58 species distributed among 11 currently recognized monophyletic taxa (all recognized as taxonomic families, except for Gekkota, which is an infraorder) and the number of species per taxon ranged from one (Crotaphytidae) to 24 (Varanidae). These 11 taxa (henceforth referred to as ‘families’ for simplicity) correspond to major branches in our phylogenetic tree (Fig. 1). Several previous studies have suggested that some of these families have relatively high V̇_O2,max and/or endurance capacity (Beck et al., 1995; Clemente et al., 2009; see Garland, 1994). Therefore, we initially tested for differences among the 11 families (Foster et al., 2018; see Gartner et al., 2010). Preliminary analysis revealed that most of the effect of family was caused by Helodermatidae, Varanidae and Scincidae, all with relatively high V̇_O2,max. In addition, models with the 11-category family variable had relatively high AICc (Akaike information criterion corrected for sample size, see below), suggesting that this variable was too fine grained. We therefore decided to instead use a set of three dummy variables coding for helodermatids, varanids and skinks as compared with all other lizards. Only models that included (or excluded) all three of these dummy variables were considered. We made the a priori decision to include log₁₀ body mass and body temperature during V̇_O2,max trials in all models, because these are known to be highly influential predictors of V̇_O2,max (see Introduction, and also Garland, 1994, regarding endurance capacity).

We fitted a total of 48 multiple regression models, including all possible combinations of the independent variables climate (3 categories), diet (5 categories), nocturnality, viviparity, and either family (11 categories) or the set of three dummy variables for putatively aerobic taxa (helodermatids, varanids and skinks). Note that models never included both family and the 3-dummy set. All 48 models were fitted by both OLS and RegOU procedures, yielding a total of 96 regression models.

For each model, we report the partial regression coefficient and P-value for each independent variable, the d parameter (see above), the ln maximum likelihood (lnML) for the model, and the AICc {AICc=(−2×lnML)+[2×p×2/(n−p−1)]}, where p is the number of parameters and n is the sample size (lower AICc values indicate better model fit to the data) (Burnham and Anderson, 2002).

To assess the importance of each independent variable, we used the following approach, consistent with a previous analysis of mammalian V̇_O2,max (Dlugosz et al., 2013). First, we calculated the Akaike weight (w_i) for each model, which is the probability that model i would be the best-fitting model, if the data were collected again under identical circumstances. The sum of w_i for all 96 models is 1. Then, for each independent variable, we calculated the sum of Akaike weights (Σw_i) for all the models that included that variable. If an independent variable is mostly present in the best-fitting models (the ones with lowest AICc), it would have a high w_i and would therefore be a relatively important variable for model fit. We calculated a weighted average of the partial regression slope for body mass and the standard error (s.e.) of the coefficient using the relative importance (measured as the Akaike weight) of each regression model as the weighting factor (Burnham and Anderson, 2002). This weighted averaging procedure allowed us to obtain the average allometric scaling exponent of V̇_O2,max across all the models within 10 AICc from the best (these top eight models accounted for 99% of the cumulative w_i), while accounting for the proportional contribution of each model based on their AICc values. We used the average standard error of the body mass partial regression coefficient to calculate the confidence intervals of the allometric slope of V̇_O2,max. The average s.e. was calculated from models that ranged in degrees of freedom (d.f.) from 55 to 58 (some models included more variables). Therefore, we considered 55 d.f. to calculate the boundaries of the 95% confidence interval of the partial regression coefficient for body mass.

Preliminary analyses tested for violation of statistical assumptions, including the possibility that log₁₀ V̇_O2,max might be related to log₁₀ body mass or to temperature in a non-linear fashion (Bueno and López-Urrutia, 2014; see Garland, 1994; Kolokotrones et al., 2010; Verberk et al., 2020). Inspection of Fig. 2 suggests a possible non-linear relationship with log₁₀ body mass; however, this visual impression is driven mainly by the fact that many of the small lizards were measured at low temperatures (blue dots in Fig. 2A). Plots of residuals from the OLS models listed in Table 1 indicated no major outliers and no clear evidence of non-linearities. As a formal test for non-linearities, we added z-transformed log₁₀ body mass squared and z-transformed temperature squared (orthogonal polynomials) to the best regression model. Neither of these partial regression coefficients was statistically significant (results not shown).

The highest correlation between independent variables occurred between the dummy variables coding for viviparity and for skinks (r=0.70), which is probably not large enough to be a concern for multicollinearity (Slinker and Glantz, 1985). However, our dataset only included six data points for viviparous species and five data points for skinks, four of which were viviparous.

We obtained eight models within 10 AICc from the best (Table 1). The difference in AICc scores from the 8th to 9th model was 5.14, which can be considered a large difference in model fit (Burnham and Anderson, 2002). Additionally, the eight best models accounted for 99% of the cumulative evidence (cumulative w_i), so we are confident that these models encompass the most accurate explanation of the (co)variation in our dataset.

Considering these eight models, the average allometric slope for V̇_O2,max was 0.803 (95% confidence interval, CI 0.747–0.859) and the average slope for temperature using the same models was 0.051 (95% CI 0.041–0.061). All these models included the three dummy variables coding for helodermatids, varanids and skinks. The best model (based on AICc values) that included ‘family’ (with its 11 categories) had an AICc of −26.94 and was placed at the 23rd position. The best model that did not include family or the three dummy variables had an AICc of −31.80 and was placed at the 10th position. The allometric slope including only body mass and temperature was 0.861 (95% CI 0.795–0.927) from OLS regression and 0.857 (95% CI 0.782–0.931) from RegOU (Table 2).

The most influential independent variables were viviparity (Σw_i=0.996), the combined set of dummy variables coding for helodermatids, varanids and skinks (Σw_i=0.996), and climate (Σw_i=0.734). Nocturnality was the next most influential variable, but its cumulative evidence was only 0.291. Diet and family with 11 categories had very low cumulative evidence (Σw_i<0.01 each). The variables coding for viviparity, Varanidae and Helodermatidae were always statistically significant (P<0.05) in the top regression models (Table 1). Viviparity had a negative effect on V̇_O2,max, whereas being part of helodermatids, varanids or skinks positively affected V̇_O2,max, including in the phylogenetic regressions (RegOU).

The best eight models included a mix of OLS and RegOU (Table 1). The d-values of the RegOU models were all below 0.013. This value is much closer to 0 than it is to 1, meaning that the nodes of the tree were considerably moved towards the root. Consequently, the tree that best explains the residual variation in this dataset is much more similar to a star phylogeny than to its original configuration. For completeness, we also present full models including all independent variables (mass, temperature, family, diet, climate, nocturnality and viviparity) in Table 2.

In the present study, we estimated that V̇_O2,max scales as body mass^0.803 among species of lizards. This allometric slope is the average of eight models that included measurement temperature (range: 20–40°C) as a covariate and also different combinations of viviparity, climate and nocturnality as independent variables (Table 1). The average 95% confidence interval of the allometric slope of V̇_O2,max from these top eight models (0.747–0.859) includes the 0.779 allometric exponent reported in our previous study restricted to lizard species measured within the 35–40°C range and estimated without regard to phylogeny (Garland and Albuquerque, 2017). It also includes the value of 0.77 reported for 24 species of diurnal lizards (Autumn et al., 1999), but not the value of 0.96 for eight nocturnal species (Autumn et al., 1999) or the 0.74 slope reported for varanids (Clemente et al., 2009). Additionally, the allometric slope reported here is very similar to the value of 0.839 previously reported for a large size range of mammals analyzed with similar phylogenetic statistical methods (Dlugosz et al., 2013).

The best eight models included four OLS and four RegOU regressions with identical combinations of independent variables, as well as log₁₀ body mass and measurement temperature. These variables represented, in order of importance: viviparity, varanids, helodermatids, skinks, climate and nocturnality (discussed below). Moreover, the estimated values for the Ornstein–Uhlenbeck parameter (d) are close to zero (Table 1), indicating little hierarchical structure to the tree. Thus, once those independent variables are included as predictors, very little phylogenetic signal (Blomberg et al., 2003) remains in the residual variation in log₁₀ V̇_O2,max.

Partial regression coefficients showed that, on average, helodermatids and varanids had greater V̇_O2,max than other lizards. These lineages are predators (Beck et al., 1995; Gienger et al., 2014; King and Green, 1993) and their high V̇_O2,max has been noted previously (Beck et al., 1995; Clemente et al., 2009). Most varanid species use mobile foraging strategies that routinely cover long distances (Clemente et al., 2009) and so they probably benefit from relatively high aerobic capacity to sustain effort. Also, Clemente et al. (2009) showed that the coadaptation of high activity levels, endurance and V̇_O2,max is probably plesiomorphic for Australian varanids (Clemente et al., 2009), so even though some of the species included here are sit-and-wait foragers, their high V̇_O2,max levels may be explained by evolutionary history.

The two living species of helodermatids (Heloderma horridum and Heloderma suspectum, both included in this study) are sporadic, binge-feeders that prey on ephemeral and sparsely distributed food resources (Beck, 2005; Beck and Lowe, 1991; Stahnke, 1950). Gienger et al. (2014) showed that individual H. suspectum have a high proportional cost of activity relative to total energy expenditure during the active season, although they are largely inactive for most of the year. Previous studies also showed that helodermatids have relatively high endurance when accounting for body mass and temperature simultaneously (see fig. 11.7 in Garland, 1994). Perhaps the high V̇_O2,max observed in helodermatids is related to relatively higher activity levels during the short period when they are active above ground in search of food or mates (John-Alder et al., 1983; Jones, 1983). However, Beck et al. (1995) also suggested that the high V̇_O2,max (and aerobic scope) of helodermatids can be attributed to strenuous male–male disputes, rather than foraging or activity levels, a suggestion that may also apply to varanids (Bickler and Anderson, 1986). Future studies could test that hypothesis by relating V̇_O2,max to foraging success or the outcome of male–male disputes at the level of individual variation.

We also found that Scincidae tend to have higher V̇_O2,max than other lizard groups (although the difference was not statistically significant in some models: Table 1). The vast majority of skinks are active and widely foraging species (Cooper, 1994a,b; Huey and Pianka, 1981). Cooper et al. (2000) showed that omnivorous and insectivorous skinks quickly respond to chemical cues from active prey, and that this behavioral trait seems to be plesiomorphic for the family. The four skink species included in the present study (Tiliqua rugosa, Oligosoma macgregori, Oligosoma nigriplantare, Plestiodon skiltonianus) are either omnivores or insectivores that would probably show the same behavior. However, data on movement patterns or foraging strategy are not available for these species (but see Kerr and Bull, 2006), so we cannot make direct comparisons of activity levels and V̇_O2,max.

Previous studies suggested that active foraging lizards have higher V̇_O2,max (or other metabolic measures) than sit-and-wait ones (Bennett et al., 1984; Clemente et al., 2009; Nagy et al., 1984). The idea that active foragers have higher V̇_O2,max is intuitive and straightforward. If an animal moves actively in search for food, then it would benefit if it could move faster without fatiguing and thus cover a larger distance in a given amount of time (see Garland and Albuquerque, 2017). Previous studies within families of lizards have reported that active foragers do indeed have higher V̇_O2,max (Bennett et al., 1984; Clemente et al., 2009; Nagy et al., 1984). Here, we show that three families that are mostly composed of active foragers (Varanidae, Helodermatidae and Scincidae) have relatively high V̇_O2,max, but we did not observe such a pattern for Teiidae (also mostly composed of active foragers). However, foraging strategies and their associated movement speeds and patterns may be better represented by a continuum, rather than a simple dichotomous categorization (see Huey and Pianka, 1981; Huey et al., 2001). Unfortunately, detailed information on movement ecology in the field requires intensive field work, which limits our ability to associate interspecific variation in natural behaviors with locomotor capacity (Boratyński et al., 2020; Garland, 1999). Hopefully, future studies will expand the databases to allow more extensive comparisons.

Viviparity is strongly related to V̇_O2,max, with an overall importance of 0.996 and is present in all of the top eight models (Table 2). However, contrary to our prediction, it had a negative effect on V̇_O2,max. Viviparity is estimated to have evolved independently ∼115 times in squamates (three of which are included in the dataset: Fig. 1) and it seems to originate most commonly in cold climates and in lizards that typically have low body temperatures (Gibson et al., 2015; Meiri et al., 2013; Mesquita et al., 2016; Pyron and Burbrink, 2014; Shine, 2005). However, the low V̇_O2,max that we observed in viviparous lizards is not related to low body temperature per se, as all models included measurement temperature as a covariate. Additional ecological or behavioral factors (e.g. activity levels) may explain this negative correlation. Perhaps the direction of the causal relationship is the reverse of what we outline in the Introduction. Pregnancy is known to reduce locomotor capacity in lizards (Garland and Else, 1987; Olsson et al., 2000; Zani et al., 2008). If viviparity has a relatively strong negative effect on locomotor capacity, then perhaps only species that do not rely on high locomotor performance could afford that extra cost.

Climate was significantly related to V̇_O2,max in four of the eight top models, but its importance was 0.734, which is below the 0.8 threshold suggested to be considered an influential variable in model averaging (Calcagno and de Mazancourt, 2010). Climate was significantly related to endurance (and possibly V̇_O2,max) among 18 species of varanids (Clemente et al., 2009). Inspection of the partial regression coefficient of each climate category showed that ‘temperate’ had a significant negative effect on V̇_O2,max. It is possible that some of the relationship between climate and V̇_O2,max in lizards is driven by temperature variation. The structural complexity of the main habitats in each climate category (e.g. forests and deserts) might be related to the remaining correlation (e.g. Cloyed and Dell, 2020; Garland et al., 1988; Killen et al., 2016; but see Scales and Butler, 2016). Arid habitats should, on average, be more open and, perhaps, scarcer in food resources, which could lead to higher activity levels associated with foraging. Tropical habitats, such as ombrophilous forests, can also impose relative physical constraints on movements, which might increase the effort required to travel a given distance. Forests are also more shaded and offer fewer basking sites, creating the need for the lizard to move around more during thermoregulatory behavior.

None of the top eight models included diet (importance=0.003), and the best model that did had an AICc of −31.4303, which is considerably higher than any of the top eight models. These results indicate that diet, at least in the categories used here, is not associated with V̇_O2,max among species of lizards. Carnivorous and insectivorous lizards are higher in trophic chains and, all else being equal, animals that have to hunt for prey should have to move around more and possibly have larger home ranges (Garland et al., 1993; Perry and Garland, 2002; Peters, 1977). These broad diet categories should also be related to ecological and behavioral traits, such as home range (Harestad and Bunnell, 1979; Perry and Garland, 2002), and activity levels (all else being equal, herbivores need to move less to obtain food). However, expecting a simple correlation between diet and V̇_O2,max may be unrealistic because, among insectivorous and carnivorous lizards, a wide range of foraging strategies exist, including extreme sit-and-wait predators (e.g. Andrews, 1979; and to a smaller extent, Gambelia wislizenii in Rose, 2004) that might move even less than herbivores and have a reduced need for aerobic capacity. Conceivably, activity levels are more closely related to how active the prey is (Huey and Pianka, 1981). Diet is also a considerably variable trait throughout seasons and among populations within species (e.g. Albuquerque et al., 2018; Kolodiuk et al., 2009), so a mismatch between the populations used for V̇_O2,max and diet data would reduce the correlation.

Even though nocturnality occurred in four of the top eight models, its overall importance was only 0.291 and it was not a significant predictor of V̇_O2,max in any of them. Because we accounted for measurement temperature in all models, this result corroborates the previous finding that the lower V̇_O2,max in nocturnal lizards is mostly explained by their lower body temperature (Autumn et al., 1999).

One limitation of the present study is that it relies on data compiled from multiple sources. In comparative vertebrate physiology, this is usually unavoidable because no single investigator or even research group is likely to be able to measure dozens or hundreds of vertebrate species. Accordingly, the field is replete with data mining and synthetic analyses. Nonetheless, such datasets will often contain errors, and these can have important effects on conclusions (Genoud et al., 2018; McKechnie and Wolf, 2004).

Several factors unmeasured in the present study may affect V̇_O2,max and add noise to the potential correlations with traits evaluated. For example, sex differences may have added noise to the data (see Beck et al., 1995). Additionally, seasonality, especially in species exhibiting territoriality or seasonal reproduction, can have a significant effect on V̇_O2,max (John-Alder, 1984; John-Alder et al., 2009). The V̇_O2,max data used here came from different studies conducted during different seasons and using individuals with different sexes (some studies report values for males and females separately, but others pool them together or do not report the sex).

In future studies, it would also be important to consider variation in other aspects of biology that may co-adapt with aerobic capacity, such as mating system and sexual selection (Bickler and Anderson, 1986; Galliard and Ferrière, 2008; Garland, 1993; Husak and Fox, 2008) or foraging mode or activity levels (e.g. as indexed by home range size; Perry and Garland, 2002), daily movement distances (Garland and Albuquerque, 2017) or movement rates (Perry, 1999; Perry et al., 1990). Additional viviparous species and measures of their activity in nature would help us to understand why this reproductive mode is associated with relatively low V̇_O2,max. Measurement of species spanning a wider range of ecological traits and behavioral traits would also be key.

Laboratory measures and experiments can also help understanding of the evolution of V̇_O2,max. Studies of social dominance among male lizards have reported a positive correlation between dominance and sprint speed (Garland et al., 1990), endurance (Perry et al., 2004), and both sprint speed and endurance (Robson and Miles, 2000). The positive correlation with endurance is particularly interesting because during endurance trials energy is provided mostly by aerobic metabolism and, therefore, endurance capacity is expected to be influenced by V̇_O2,max (see Garland, 1984, on individual variation in lizards). Additionally, in the dominance studies mentioned, lizards were placed in relatively small arenas and had to compete for a single basking site. That is a helpful setup if we want to understand what immediate factors can influence the outcome of agonistic interactions. However, natural encounters may differ in various ways. For example, in their natural habitat, lizards can retreat for short periods and then re-engage, which could increase the overall length of an interaction. When fighting for a single basking site, whichever individual wins the first round has the benefit of keeping optimal temperature, which could confer an advantage if the adversary tries to re-engage after a few minutes. In nature, both lizards have the opportunity to maintain optimal temperature and keep fighting for a longer period.

We would like to thank Christopher Clemente and Graham Thompson for providing the collection location for 9 varanid species, and Marcell Cadney and Jessica Tingle for providing comments to earlier versions of the manuscript before submission. Most of the results and discussion in this paper are reproduced from the PhD dissertation of R.L.A. (Albuquerque, 2019).