Macroecology and macroevolution of body size in Anolis lizards

812 –––––––––––––––––––––––––––––––––––––––– © 2020 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Ken Kozak Editor-in-Chief: Miguel Araújo Accepted 15 January 2020 43: 812–822, 2020 doi: 10.1111/ecog.04583 43 812–822


Introduction
Body size is one of the most influential biological traits affecting ecological and physiological characteristics across all plants and animals (Peters 1983). Understanding the tempo and mode of body size evolution and the resulting present-day geographical gradients in body size are goals of the fields of macroecology and macroevolution (Smith and Lyons 2013). Researchers in these fields have postulated general trends in body size such as increase (Bergman's rule) or decrease (inverse Bergman's rule) across latitudes (Ashton and Feldman 2003) and increases through time (Cope's rule, Jablonski 1997, Clapham andKarr 2012). None of these trends has been shown to apply generally across taxa, but each has found some support in some lineages (Jablonski 1997).
Macroecological studies have usually asked whether geographical body size gradients (e.g. mean body sizes across assemblages) are associated with spatial environmental 813 variation (Terribile et al. 2009, Slavenko et al. 2019). These studies typically use an assemblage-based approach wherein body size is averaged for each set of species occurring in a single grid cell or latitudinal band . The resulting set of body size averages for each species assemblage is tested for correlation with geography (e.g. latitude) and/or climate across space. These approaches often uses phylogenetic comparative methods (PCMs) to evaluate the contribution of phylogenetic effects from the resulting geographical body size gradients but they do not explicitly address questions about macroevolutionary processes. That is, PCMs in macroecology have been commonly used for 'correction for phylogeny' in statistical analyses and not necessarily to infer evolutionary processes through deep temporal scales (Slavenko et al. 2019). Conversely, macroevolutionary approaches based on PCMs allow identifying changes in tempo (i.e. rates) and mode (e.g. process) of body size evolution over time (Harmon et al. 2010, Slater 2013, Clavel and Morlon 2017. These approaches, unlike those based on the fossil record (Marx andUhen 2010, Lehtonen et al. 2017), usually ignore spatial and/or climatic effects on trait diversification. Some macroevolutionary studies evaluates differences in trait evolution among regions/environments but fails to include spatially-explicit approaches as usually done in macroecological studies. However, some recent studies have begun to incorporate both spatial and climatic effects as predictors of trait macroevolution (Clavel andMorlon 2017, Caetano et al. 2018).
The integration of macroecological and macroevolutionary perspectives has been a challenging for researchers attempting to evaluate how biodiversity patterns emerge across spatial and temporal scales (Bromham andCardillo 2019, McGill et al. 2019). Most studies have usually explored each dimension, spatial or temporal, separately (Smith and Lyons 2013, Amado et al. 2019, Slavenko et al. 2019. Few studies have attempted to integrate both spatial and temporal dimensions simultaneously to evaluate geographical and evolutionary body size patterns (Carotenuto et al. 2015, Villalobos et al. 2017). Although these perspectives explore different dimensions of phenotypic patterns (e.g. rates of body size evolution versus spatially-averaged body size gradients), the simultaneous integration of both approaches allow us to understand how traits have evolved across time and space and how past and current climate have shaped the current patterns (McGill et al. 2019). The macroecological component describes how spatial gradients in phenotypic traits are related with past and current climatic variation. The macroevolutionary component describes the historical processes that produced the macroecological patterns we view today (e.g. rapid trait evolution can generate strong body size gradients across geography). Because some potential factors influencing body size patterns (e.g. temperature) vary across time and space these two approaches are complementary. For example, if both spatial patterns of body size and its evolutionary rates are correlated with climate (e.g. temperature), then we can conclude that climate have simultaenously influenced both spatial and evolutionary dimensions of body size. By contrast, if historical, but not spatial variation, in temperature explain patterns of body size, then phylogenetic effects can be more important than geography in shaping these trait patterns. Currently, a single joint statistical framework integrating spatial and temporal dimensions is lacking, but conceptual advances have been made in this direction (Diniz-Filho et al. 2013, Fritz et al. 2013, Lehtonen et al. 2017. A joint examination of geographical and phylogenetic components of body size, and their relationships with climatic factors, might illuminate how current body size patterns have emerged.
Anolis lizards are suitable for examining how body size patterns have emerged under both spatial and temporal dimensions. This clade is a classic example of adaptive radiation, diversifying extensively across multiple insular and mainland habitats (Losos 2009, Poe et al. 2017. Although several studies have evaluated the role of ecological interactions on body size on Caribbean anoles (Losos 2009, Mahler et al. 2010, the role of climate driving spatial body size gradients and body size evolution has been little explored (Muñoz et al. 2014a, Poe et al. 2018). In addition, spatial and phylogenetic body size variation in anoles might be related to species-or lineagespecific thermoregulatory strategies (e.g. thermoconformer [i.e. those whose body temperature fluctuates with environmental temperatures] versus thermoregulator [i.e. those with the ability to regulate their body temperature]; Huey andSlatkin 1976, Muñoz andLosos 2018). As a consequence, contrasting thermoregulatory strategies can explain both direct Bergmann's clines (i.e. increases in body sizes toward cold areas) as inverse clines (i.e. decreases in body size toward cold areas; Ashton and Feldman 2003, Sears and Angilleta 2004, Olalla-Tárraga 2011. Olalla-Tárraga and Rodríguez (2007; see also Olalla-Tárraga 2011) proposed the heat balance hypothesis to explain both direct Bergmann's clines (i.e. increases in body sizes toward cold areas) as inverse clines (i.e. decreases in body size toward cold areas) based on thermoregulatory strategies. Thermoconformer species (i.e. those whose body temperature fluctuates with environmental temperatures) occurring in colder regions are predicted to have small body sizes (i.e. a converse Bergmann's rule) as an adaptation to reduce heating times via their higher surface-volume ratio, thus being more efficient in energy expenditure (Ashton and Feldman 2003, Olalla-Tárraga and Rodríguez 2007, Terribile et al. 2009). By contrast, thermoregulator species (i.e. those with the ability to regulate their body temperature) occurring in colder regions can exhibit both body size clines. On the one hand, it is expected that large body sizes occur in colder regions being facilitated by heat conservation mechanisms given via their lower surface-volume ratio (Olalla-Tárraga and Rodríguez 2007, Olalla-Tárraga 2011). On the other hand, small body sizes are also expected in colder regions facilitated by rapid heat gain due to a higher surface-volume ratio (Olalla-Tárraga and Rodríguez 2007, Olalla-Tárraga 2011).
Montane lizards from the A. cybotes clade in Hispaniola are an example of this last case as they have small body size in mountain areas and heat themselves actively during the day (Muñoz et al. 2014a, b). It is unknown if this is a particular case or a general trend in Anolis lizards. We might expect that whether spatial and phylogenetic patterns of body size in Anolis lizards are driven by temperature, these patterns can be explained by underlying thermoregulatory strategies. However, the role of environmental factors on anole body sizes at large spatial and temporal scales remains unknown.
Here using macroecological and macroevolutionary approaches we evaluate the role of past and current environmental conditions on geographical and evolutionary patterns of body size in the entire Anolis lizard radiation. First, we evaluate the role of temperature, primary productivity and climatic seasonality on geographical patterns of anole body size. We test three classical ecological-physiological hypothesis driving geographical body size gradients (Ashton andFeldman 2003, Terribile et al. 2009): heat balance, primary productivity and seasonality. For heat balance, we predict that anole assemblages can exhibit either increase or decrease of body sizes towards colder regions depending of the thermoregulatory strategy (see above). For primary productivity, we predict that anole assemblages are expected to harbor species with smaller body sizes as a result of decreased food availability (i.e. in sites with low primary productivity) limiting growth rate and thus body size at sexual maturity (Terribile et al. 2009). An equally plausible hypothesis is that higher primary productivity is also associated to sites with complex vegetation structure where habitat specialists such as twig and crown-giant anoles can occur (Losos 2009). For seasonality, we predict that anole assemblages under high environmental seasonality maintain larger species relative to assemblages exposed to low environmental seasonality facilitated by weight-specific metabolic rates which engage long periods of starvation (Olalla-Tárraga et al. 2006, Terribile et al. 2009). Second, we test whether rates of body size evolution were affected by past temperature changes through time. We used a set of phylogenetic models including temperature-dependent, time-dependent, geographical-dependent and constant-rate models (Condamine et al. 2013, Clavel and Morlon 2017, Lewitus and Morlon 2018 to evaluate how body size evolution occurs through time and geography. In addition, we used a set of paleo-inspired models that test for specific shifts in tempo and mode of body size evolution after a specific time horizon (Slater 2013). These complementary (macroecological and macroevolutionary) approaches allow us to discern the role of climate in body size across spatial and phylogenetic scales.

Phylogenetic and body size data
For our comparative analyses, we used the most recent and complete phylogenetic hypothesis of Anolis lizards (379 species) from Poe et al. (2017;hereafter, MRT-tree). We also used the absolute time-calibrated phylogeny for 322 species based only on molecular data from Poe et al. (2017; hereafter MCC-tree). To account for phylogenetic uncertainty, we selected 100 trees from the Bayesian posterior distribution of the complete (hereafter, set-MRT-trees) and the absolute calibrated (hereafter, set-MCC-trees) trees. Body size data were taken from measurements of specimens in collections, and in a few cases from the literature. We compiled maximum snout to vent length (SVL) for all known species of Anolis, including males and females (n = 1-15 specimens per species). A recent study conducted by Armstead and Poe (2015) found that low sample sizes, even n = 1, are adequate for comparative interspecies analysis in Anolis lizards. We selected maximum SVL because it is an indication of the maximum size that a species can reach (Andrews 1982). SVL was log10transformed before analyses. All measurements were taken by the same person (SP) to avoid potential measurer effects.

Assemblage-based approach
We generated range maps for each of the 379 species based on occurrence records using minimum convex polygons and removing non-land areas afterwards (see Supplementary material Appendix 1 for details). We counted the number of species overlapping each grid cell by overlaying a grid of 1 × 1° (lat-long) resolution covering the region (Fig. 1). We also used points and species' range maps from Roll et al. (2017) to explore sensitivity of conclusions to alternative geographical data. Then, using an assemblage-based approach (Olalla-Tárraga et al. 2010) focusing on the grid-cell assemblages as observational units, we computed the mean and median log10-transformed SVL values for all species present within each grid cell. Results from mean and median body size values were very similar; therefore, we only discuss results from median values. We tested three environmental hypotheses suggested in the literature as potential explanations for geographical patterns of body size variation in non-avian reptiles (Ashton andFeldman 2003, Terribile et al. 2009) as follows: 1) heat balance (HEATBAL); 2) primary productivity (PRIMPROD); and 3) seasonality (SEAS). We used the following environmental variables grouped for each hypothesis: 1) HEATBAL: annual mean temperature (Hijmans et al. 2005); 2) PRIMPROD: actual evapotranspiration (< http:// csi.cgiar.org >) and 3) SEAS: seasonality of temperature and precipitation (Hijmans et al. 2005, Supplementary material Appendix 1 Fig. A1). These coarse-grain variables have been used extensively in macroecological studies examining their roles on geographical functional traits gradients, including body size (Olson et al. 2009, Terribile et al. 2009, Amado et al. 2019, Ochoa-Ochoa et al. 2019, Slavenko et al. 2019). We evaluated a series of spatial autoregressive models (SAR) to test each hypothesis on the relationship between environment and geographical body size gradients in Anolis lizards using the spdep and ncf packages with a script from Kissling and Carl (2007). The spatial weights were estimated from a neighborhood matrix based on the presence-absence matrix for anoles and a maximum distance of 1000 km. SAR models controlled spatial autocorrelation effects in residuals much better than OLS models (Supplementary material Appendix 1 Fig. A2). We selected the model with the highest explanatory power using the Akaike information criteria (AIC) and calculated the Nagelkerke pseudo-R 2 for each model. Finally, using a null modeling approach we evaluated whether our results were sensitive to pseudo-replication effects due to the problem of repeated co-occurrences in the assemblage approach (Hawkins et al 2017). For this, we randomized body size values (median log10-transformed SVL values) across species and generate 100 random spatial body size gradients. We performed SAR models as above on these randomized body size gradients and evaluated whether the observed Nagelkerke pseudo-R 2 values departures from values from 100 random Nagelkerke pseudo-R 2 .

Cross-species approach
Complementary to the assemblage-based approach, we also performed a cross-species analysis focusing on species as observational units. We used species' occurrences to calculate the average of each environmental variable across the range of each species. Then, we used the average environmental value of each species to perform regressions against species' maximum SVL. To account for the phylogenetic non-independence among species, we implemented a PGLS regression (Grafen 1989) between body size and environmental variables, averaged by species, for each hypothesis (i.e. HEATBAL, PRIMPRO and SEAS). Based on single values per species (across their individual ranges), this cross-species approach also avoids the repeated co-occurrence problem of the assemblage approach that has been shown to inflate type I error rates when relating aggregated species attributes (Hawkins et al. 2017). We performed the PGLS models under a Brownian motion model using the pgls function in the caper R package (Orme et al. 2012).

Macroevolutionary models of body size evolution
We fitted a series of macroevolutionary models to evaluate whether tempo (i.e. evolutionary rates) and/or mode (i.e. Brownian motion or Ornstein-Uhlenbeck processes) of body size evolution changed through time. These models include climatic-, geographical-and time-dependent models as well as models with constant evolutionary rates through time. For temperature-dependent models, we used recent models developed by Clavel and Morlon (2017) that are an extended BM model with time-varying evolutionary rates influenced by an environmental variable. This model has been used to test for an association between rates of body size evolution (i.e. σ 2 ) and temperatures throughout the Cenozoic in mammals and birds (Clavel and Morlon 2017). The method was implemented in the RPANDA R package (Morlon et al. 2016) and we used two versions of default functional forms of this model using the 'fit-t-env' function, one with an exponential association between evolutionary rates and temperature (the 'temperature-exponential' model where evolutionary rates are affected exponentially by temperature) and the other with a linear association (the 'temperature-linear' model where evolutionary rates are affected additively by temperature; see Clavel and Morlon 2017 for further details). A single parameter, beta, measures the strength and direction of temperature dependence. Positive beta values indicate that rates of body size evolution have decreased as a response to temperature. Negative beta values indicate that rates of body size evolution have increased as a response to temperature. We also tested whether the way the temperature curve was fitted affected parameter estimates. We tested three degrees of increased data smoothing using cubic splines. Temperature data were obtained from Cramer et al. (2011) available in the 'RPANDA' R package.
For geographical-dependent models, we fitted a series of models incorporating geographical information to establish whether evolutionary rates vary among selective regimes (i.e. geographical regions). In this case, we assumed that clades colonizing these regions were subject to strong selective regimes and that promoted shifts in rates of body size evolution according to the ecological opportunity hypothesis (Stroud and Losos 2016). Accordingly, if dispersal events influenced body size evolution then different evolutionary rates for body size in each one of the geographical locations are expected. We assigned species to the following geographical locations: mainland including endemics to islands close to mainland (e.g. Cocos, Malpelo and Gorgona island in the Pacific Ocean); Southern Lesser Antilles, Northern Lesser Antilles and Greater Antilles (including Cayman and Bahamas islands). We estimated geographical ancestral states across the anole lizard phylogeny using 100 histories of stochastic character mapping with the make.simmap function in the phytools R package (Revell 2012). We fitted seven different models of body size evolution for each geographical regime with the OUwie R package (Beaulieu and O'Meara 2019).
In addition, we fitted three paleo-informed models developed by Slater (2013) to evaluate whether body size evolution changed after a given temporal event. Specifically, we asked whether after a particular time horizon, indicating an abrupt climate change event, both tempo (i.e. rates) and mode (i.e. process; e.g. BM or OU) of body size evolution changed. The first model, 'shift', tests whether there is a significant shift in rates of body size evolution (i.e, σ 2 ; only tempo) before or after a temporal event. Accordingly, this model allowed us to detect if there was an increase or decrease in evolutionary rates (i.e. body size diversification, as opposed to a directional trend) after a given time horizon. The second model, 'ecological release', tests whether evolutionary rates change in mode after a given temporal event (i.e. only mode). Under this model, body size is modelled as a single-peak Ornstein-Uhlenbeck -OU -process before a given temporal event constraining its evolutionary rate; after that event, the trait is modeled as an unconstrained Brownian process with increases of evolutionary rate over time. In this model, the alpha parameter of the OU process is removed after the temporal event (Slater 2013). This model allowed us to capture an evolutionary dynamic where a constraining process is governing body size evolution (i.e. an OU process constraining body size variance) before a given temporal event and after that event body size evolved with increasing variance (i.e. a BM process increasing body size variance). This model was developed by Slater (2013) to test whether mammalian body size evolution was released after the Cretaceous-Paleogene extinction event. The third model, 'release and radiate', tests whether body size evolution varies in tempo (evolutionary rate; σ 2 ) and mode (process; e.g. a BM or OU model) simultaneously after a temporal event. In this model, body size evolution is modeled as a single-peak OU process but with a different evolutionary rate before and after a temporal event, with evolutionary rates increasing following an unconstrained Brownian motion process. In this model, the σ 2 parameter may change after the temporal event. Note that body size variance in these models refer to the variability in mean (or median) body size across species through time and not to body size variance within assemblages. We adjusted these models to evaluate whether tempo and mode in rates of body size evolution in Anolis shifted after three global climate change events during the Cenozoic: 1) Eocene Thermal Maximum (hereafter -ETM -52 mya; Payros et al. 2015); 2) Eocene-Oligocene Glacial Maximum (hereafter -EOGM -33 mya; Lear et al. 2008); and 3) Middle Miocene Climatic Optimum (hereafter -MMCO -14 mya; You et al. 2009).
Finally, we fitted a series of standard models with singlerates or single-optimum models as follows: 1) a single-rate BM model which assumes that rates of body size evolution are invariant through time (i.e. constant); 2) a single-optimum OU model which assumes that rates of body size are invariant through time (i.e. constant); 3) a Early Burst -EB -model which assumes that rates of body size evolution increases or decreases exponentially through time (Harmon et al. 2010); 4) a trend model which is a diffusion model with linear trends in rates through time (i.e. toward larger or smaller body sizes); 5) a lambda model which assumes that phylogenetic trees predicts covariance among trait values for species (Pagel 1999); 6) a kappa model which assumes that traits evolves according to a punctuational speciation model (Pagel 1999); 7) a delta model which assumes that trait evolution is timedependent (Pagel 1999). Delta values greater than 1 imply that recent body size evolution has been fast; by contrast, delta values less than 1 imply body size evolution has been slow; 8) a drift model which assumes that traits evolves with a directional drift component either towards smaller or larger values; 9) a white-noise model which assumes that trait evolution is decoupled from the phylogeny (i.e. trait dataset has no covariance structure among species). These models were fitted using the 'fitContinuous' function from the 'geiger' R package (Harmon et al. 2015). We compared the statistical support of all these 26 macroevolutionary models using Akaike information criteria. The best selected model was the one with the lowest AICc. We also compared the statistical support for all models using Akaike weights.

Assemblage and cross-species approaches
We did not find a clear geographical gradient of body size in Anolis lizards (Fig. 1). However, we observed that assemblages with smaller body size values were concentrated in Middle America, mainly in the Pacific coast. Using the assemblage-based and cross-species approach, we did not find any support for any of the three hypothesis tested. Although the pseudo-R 2 Nagelkerke values tend to be relatively high (Table 1), the fitted regression lines were weak (Fig. 2). In addition, the null model approach implemented here revealed that the observed R 2 Nagelkerke values are not different from random R 2 Nagelkerke values (Fig. 3). Furthermore, results from the assemblage-based approach were robust to different methodological approaches to generate species' ranges (Supplementary material Appendix 1 Fig.  A3). All these results suggest a very weak or absent signature of current environmental variables on the observed body size values across coarse-grain geographical scales (Fig. 2, 3).

Macrevolutionary models of body size evolution
We found a strong support for a climate-dependent model of body size evolution ( Table 2). The best selected models (i.e. with the lowest AICc value and highest Akaike weight) were those which explicitly considered an association between paleo-temperature changes throughout the Cenozoic and rates of body size evolution (Fig. 4a-b). In particular, the best model selected was one with an exponential temperature dependence on rates of body size evolution ( Table 2). Constant-rate models (i.e. those with evolutionary rates invariant through time) and time-dependent models (e.g. a delta model) had a relatively poor statistical support. In addition, paleo-inspired models testing whether tempo and mode of body size evolution changed after a particular abrupt climate change also had a relatively poor statistical support in comparison with climate-dependent models (Table 2). Finally, the geographical-dependent models also had a relatively poor statistical support suggesting that rates of body size evolution are not related with ancient dispersal events ( Table 2). All these findings were consistent regardless of whether relative or absolute-time calibrated trees were used (Supplementary material Appendix 1 Table A1). The strength and direction of temperature-dependence captured in the estimated beta parameter was consistently negative regardless of phylogenetic uncertainty (Supplementary material Appendix 1 Fig. A4) and male and female body size (Supplementary material Appendix 1 Fig. A5). These results suggest that rates of body size evolution increased as the climate became colder instead of after colonization of new areas.

Discussion
The evolution of body size has been one of the most studied topics in evolutionary biology (Smith and Lyons 2013). Few studies have attempted an explicit consideration of both spatial and temporal dimensions on body size evolution (Carotenuto et al. 2015, Villalobos et al. 2017, Rozzi 2018).
Here we found evidence that macroecological and macroevolutionary patterns are influenced differently by temperature. Table 1. Model summary for spatial autoregressive models (assemblage-based approach) and phylogenetic generalized squares models (cross-species approach) linking environmental variables with body size in Anolis lizards across geography and species. R 2 values for spatial autoregressive models corresponds to pseudo-R 2 Nagelkerke values. Bold text correspond to models with the lowest AICc values (Delta AICc = 0).  On one hand, spatial body size gradients are not related with current environmental variables. In addition, median body size assemblages do not to display a geographical gradient at a coarse-grain scale. On the other hand, we found that temperature had a strong effect on the rates of body size evolution through time. Evolutionary increases and/or decreases in body size of anole lizards are modulated by climatic cooling across the Cenozoic. We did not find a strong support for an ecological and/or physiological mechanism driving spatial body size gradients in Anolis lizards. By using an assemblage-based approach and a cross-species approach, we did not find a support for any ecological-physiological mechanism driving geographical body size gradients in Anolis lizards. A similar finding for squamate reptiles at a global scale was recently found by Slavenko et al. (2019), who reported a lack of an effect of current climate on spatial body size gradients across phylogenetic and spatial scales. It is probable that these lack of environmental signatures on body size are a direct reflection of two related phenomena. First, the macroecological analyses conducted here adopted a coarse-grained scale that likely does not capture body size differentiation associated with environmental variation at finer spatial scales (e.g. intraecomorph categories). For instance, Muñoz et al. (2014a) found inverse body size clines across elevation in A. cybotes and A. sagrei species complexes driven by interspecific and intraspecific variation, respectively. Given that each clade exhibits a different manifestation of trait-environment relationships, Muñoz et al. (2014a) suggested that Bergmannian patterns can emerge from disparate ecological and evolutionary processes. A similar pattern might be recovered here for mainland anoles but with some notable differences between clades. There are species in the Dactyloa and Draconura clades of the mainland Anolis radiation that exhibit large (e.g. species with SVL > 100 mm: A. apollinaris, A. danieli, A. eulaemus, A. fraseri, A. ginaelisae, A. microtus, A. petersi, A. vanzolinii) as well as small (e.g. species with SVL < 50 mm: A. antonii, A. forbesi, A. lamari, A. milleri, A. omiltemanus, A. pseudopachypus, A. quercorum) body sizes in higher elevations (e.g. elevations above 1000 m). Further analyses with these species/clades occurring in wide elevational gradients might help to highlight the role of environmental/physiological mechanisms driving body size variation under an explicitly phylogenetic comparative framework. Second, the low spatial variation of the environmental factors analyzed here across the clade's geographic distribution likely influences the weak environmental signal on geographical body size gradients.

Heat balance
Random r2 values The role of climatic events on phenotypic diversification has been little explored in the otherwise well studied Anolis radiation. We found a strong link between cooling temperatures and acceleration of rates of body size evolution over macroevolutionary time scales (Table 2). Similar findings linking cold climates with increases in rates of body size evolution have been found for extant mammals and birds (Clavel and Morlon 2017). Clavel and Morlon (2017) suggest that cold climates accelerating evolutionary rates are consistent with previous findings where rates of phenotypic evolution are higher in temperate latitudes (Cooper andPurvis 2010, Lawson andWeir 2014). However, for ectotherms, as in this case for Anolis lizards, physiological mechanisms underlying these evolutionary patterns are very different from endotherms (Rolland et al. 2018). In particular, we suggest that the link between cold temperatures and body size evolutionary Table 2. Model selection for macroevolutionary models of body size evolution in Anolis lizard through time and space. The models included were constant-rate, punctuated, adaptive radiation, time-, geographical-and temperature-dependent and paleo-inspired models (see main text for a full description of models). AICc: Akaike information criterion corrected.  Salazar et al. (2019) found that heat tolerance and field-active body temperatures differs between mainland and insular Anolis species suggesting notable differences in behavioral thermoregulation among regions. In addition, mainland species evolved faster in heat tolerance and field-active body temperature than insular species (Salazar et al. 2019), which was likely driven by differences in available climate space between regions (Velasco et al. 2016). In addition, Muñoz and Losos (2018) found a complex pattern of behavioral thermoregulation across an elevational gradient for the A. cybotes complex (the Cybotoid clade). Cybotoid anoles in these high elevation habitats are also active thermoregulators regardless of the limited availability of thermal space. By contrast, Cybotoid anoles in lowlands are larger and thermoconformers (Muñoz and Losos 2018). Muñoz and Losos (2018) found also that cybotoid anoles have shifted their behaviors to select habitat types (e.g. rocky substrates) with the aim to maintain their preferred body temperatures across different elevational ranges in Hispaniola. Cold climates may drive body size evolution in Anolis via specific thermoregulatory adaptations in each species/clade according to the heat balance hypothesis. We suggest that thermoregulating lineages in Anolis have evolved toward large sizes due to the physiological benefits of a reduced surfaceto-volume ratio. The heat balance hypothesis predicts that a reduced surface-to-volume ratio facilitates heat conservation (thermal inertia) in cold climates (Meiri 2011). This relationship might explain, for instance, why most large species from the Dactyloa clade tend to occupy high elevations in mainland areas (e.g. in the Andes region of Ecuador and Colombia). In contrast to the effect for thermoregulators, we predict that thermoconformer lineages in Anolis have evolved towards smaller body sizes due to the physiological benefits of a greater surface-to-volume ratio for ectotherms in warm climates. It is well-known that greater surface-to-volume ratios enhances heat gain in cold climates (Meiri 2011, Olalla-Tárraga 2011. This effect could also explain why most small species from the Draconura clade and some small species from the Dactyloa clade of Anolis occupy highlands in Middle and South America. However, to test these hypotheses, it is necessary to adopt an inter-and intra-specific approach with multiple species from both clades across a geographical gradient (Muñoz et al. 2014a). These non-mutually exclusive predictions might be tested using explicit phylogenetic comparative methods. For instance, whether a given thermoregulatory strategy (thermoconformers versus thermoregulators) promotes body size/species diversification in these two mainland clades could be evaluated. Additional research focus evaluating how these contrasting thermoregulatory mechanisms extend across spatial and phylogenetic scales in this iconic lizard radiation.
A recent study conducted by Poe et al. (2018) found evidence that mainland and insular areas exhibited similar rates of body size evolution. We fond similar results here, with no evidence of differences in climatic signatures between mainland and insular anole lineages. Accordingly, these results suggest a major role of past temperatures driving body size evolution in Anolis lizards either toward larger or smaller body sizes independently of geography. The complex dispersal history of anoles lizards across Caribbean islands, recently addressed by Poe et al. (2017), likely play a role on body size diversification driven by ecological opportunity (Mahler et al. 2010, Algar andMahler 2016). Although it remains an open question whether body size evolution was affected by multiple dispersal events inferred for Caribbean and mainland landmasses, our results showed that rates of body size are independently of geographical location. Further studies including geological information (e.g. sea level rise, mountain uplifts) should explicitly incorporate these biogeographical scenarios to test these macroecovolutionary hypotheses. In a similar vein, competitive interactions have been found to influence body size diversification in Anolis, particularly in Caribbean islands (Losos 2009, Mahler et al. 2010). However, the relative influence of such biotic factors in comparison with abiotic factors such as environmental conditions remain to be tested in order to discern if both factors have shaped trait diversification in Anolis lizards similarly (or not) across phylogenetic and spatial scales (Benton 2009).
Finally, we have shown that the adoption of macroecological and macroevolutionary approaches may be highly informative for explaining clade-wide patterns of body size variation across space and time. Although recent conceptual studies have emphasized the need to adopt a consolidated statistical framework (Pearse et al. 2018), such a methodological merger remains a work in progress. For instance, it should be clear that macroecological approaches tend to explore how environments drives the spatial variation in observed body size (e.g. mean or median) and macroevolutionary approaches are focused to explore how evolutionary rates (i.e. variance) changes through temporal scales (McGill et al. 2019). Although we recognize this limitation in our own study, we believe that the combination of both approaches allowed us to discern how current and past climate affect different body size dimensions. However, we suggest that new phylogenetic comparative methods integrating time and space explicitly as explanatory variables need to be developed, particularly when considering their joint effect on a particular trait. We suggest that the availability of paleo-environmental data covering the entire Cenozoic across regional scales (e.g. Neotropical realm) will help to evaluate macroecological mechanisms at deep time scales in Anolis lizards. The current lack of data describing the spatio-temporal variation in climate precludes a full integration of macroecological and macroevolutionary perspectives in a single joint framework. Compiling these data will be a significant challenge, but they are necessary for a comprehensive understanding of how environment drives spatial and temporal patterns of body size.