No Evidence for Pace of Life Evolution Along Elevational Gradients in Squamate Reptiles

TL;DR

The study finds no clear evidence that life-history traits in squamate reptiles evolve consistently along elevation gradients as predicted by the pace of life syndrome hypothesis.

Contribution

The research challenges the universality of the pace of life syndrome hypothesis by showing inconsistent trait patterns in squamates across elevation.

Findings

Most traits in squamate reptiles showed no significant elevational patterns.

Snakes exhibited significant elevational differences in body size and fecundity, while lizards did not.

Elevation-latitude interactions suggest life-history evolution is taxonomically constrained.

Abstract

Ecological conditions can significantly influence the trade‐off between survival, growth and reproduction, driving life‐history divergence among populations in different environments. The pace of life syndrome hypothesis (POLS) proposes that additional suites of traits, from physiology to behaviour, adaptively co‐evolve with life‐history traits along a slow‐fast continuum. Using data from 192 studies representing 104 squamate species, we performed phylogenetic meta‐analyses testing whether traits vary predictably across elevation within species in accordance with POLS. Results show there is no clear evidence for an overall intraspecific elevational pace‐of‐life syndrome in squamate reptiles. While high elevation populations had significantly lower body temperatures and larger egg sizes, most traits—including body size, longevity, fecundity and thermal tolerance—exhibited non‐significant elevational patterns. Critically, lizards and snakes responded fundamentally differently: snakes showed significant differences in adult female and neonate body size and fecundity across elevation, while lizards showed no significant divergence. Elevation‐latitude interactions provide further evidence against a single pace‐of‐life solution to elevation. Our findings challenge the applicability of POLS theory within species, revealing that life‐history evolution is more context‐dependent and taxonomically constrained than syndrome‐based approaches suggest. These results contribute to growing evidence that universal trait syndromes rarely emerge across environmental gradients in diverse taxonomic groups.

Ecological conditions can significantly influence the trade‐off between survival, growth, and reproduction, leading to evolutionary divergence in life‐history traits among populations experiencing different environments. The pace of life syndrome hypothesis (POLS) proposes that additional suites of traits, from physiology to behaviour, adaptively co‐evolve with life‐history traits along this continuum. Using data from 192 studies representing 104 species, we find no clear evidence for an overall intraspecific elevational POLS in squamates, with most traits showing no significant patterns and a fundamental divergence between lizards (no response) and snakes (significant differences in adult female and neonate body size and fecundity). These findings challenge the broad applicability of POLS theory, suggesting that life‐history evolution is more context‐dependent and taxonomically constrained than universal syndrome‐based approaches predict.

Figures

Summary of mean responses (log response ratio, lnRR) and 95% confidence intervals for life‐history and physiological traits of squamates in response to elevational differences. Individual effect sizes are displayed as points scaled by precision (1/SE), where larger points indicate more precise estimates. The mean effect is shown as the central, outlined point and the black whiskers provide the 95% CI For each trait, k indicates the number of effect sizes, with the number of studies in parentheses.

Summary of mean variation (coefficient of variation, lnCVR) and 95% confidence intervals for life‐history and physiological traits of squamates in response to elevational differences. Individual effect sizes are displayed as points scaled by precision (1/SE), where larger points indicate more precise estimates. The mean effect is shown as the central, outlined point and the black whiskers provide the 95% CI For each trait, k indicates the number of effect sizes, with the number of studies in parentheses.

Tables

TABLE 1. Mean response ratios (β 0) modelled for each trait using multi‐level meta‐analyses without moderators. For each trait, the number of studies, number of effect sizes (k), standard errors (SE), 95% confidence intervals, significance (p) and heterogeneity (I 2 ) are summarised. In cells lacking I 2 values, the corresponding random factor was not included in the model based on the percentage of heterogeneity explained (see Methods section).

	Trait
	Adult body size	Adult female body size	Neonate body size	Egg size	Fecundity	Mean age	T b	CTmin	CTmax
N studies	129	93	26	11	63	18	40	14	20
k	250	106	45	16	130	18	46	26	34
Initial prediction	↑ with elevation	↑ with elevation	↑ with elevation	↑ with elevation	↓ with elevation	↑ with elevation	↓ with elevation	↓ with elevation	↓ with elevation
β 0	0.030	0.006	0.033	0.106	−0.053	0.077	−0.052	−0.115	−0.006
SE	0.029	0.013	0.024	0.047	0.076	0.098	0.009	0.093	0.008
p	0.306	0.671	0.188	0.038	0.489	0.441	< 0.001	0.228	0.462
Lower CI	−0.028	−0.021	−0.017	0.007	−0.206	−0.129	−0.070	−0.307	−0.022
Upper CI	0.089	0.032	0.083	0.206	0.100	0.284	−0.034	0.077	0.010
I 2 effect (%)	35.48	—	57.99	98.15	36.8	11.2	—	—	—
I 2 study (%)	—	5.76	20.52	—	—	11.2	66.1	—	—
I 2 species (%)	11.76	5.3	—	—	—	11.2	3.49	36.1	80.7
I 2 phylo (%)	—	2.31	—	—	26.04	46.6	—	57.8	14.9
I 2 location (%)	47.41	84.37	20.52	—	35.17	11.2	27.81	—	—
I 2 total (%)	99.78	97.73	99.04	98.15	98.02	91.2	97.41	94.18	95.6

TABLE 2. Summary of the moderators explaining the influence of elevation range (elvR) and latitude difference (latD) on the selected mean life‐history and physiological traits.

Moderator	Adult body size			Adult female body size			Neonate body size
	Est	SE	p	Est	SE	p	Est	SE	p
elvR	−0.013	0.028	0.659	−0.006	0.017	0.727	0.002	0.043	0.959
latD	0.017	0.011	0.122	0.021	0.007	0.005	−0.002	0.017	0.904
elvR × latD	−0.011	0.011	0.357	−0.002	0.008	0.808	−0.012	0.022	0.590

TABLE 3. The coefficient of variation (β 0) modelled for each trait using multi‐level meta‐analyses without moderators. For each trait, the number of studies, number of effect sizes (k), standard errors (SE), 95% confidence intervals and significance (p) and heterogeneity (I 2 ) are summarised. In cells lacking I 2 values, the corresponding random factor was not included in the model based on the percentage of heterogeneity explained (see Methods section).

	Trait
	Adult body size	Adult female body size	Fecundity	T b
N studies	129	85	63	40
k	250	89	130	46
Initial prediction	↑ variation at higher elevations	↑ variation at higher elevations	↑ variation at higher elevations	↓ variation at high elevations
β0	0.034	0.039	0.036	0.209
SE	0.036	0.041	0.041	0.084
p	0.343	0.346	0.382	0.017
Lower CI	−0.037	−0.044	−0.046	0.039
Upper CI	0.106	0.122	0.118	0.379
I 2 effect (%)	40.5	—	55.31	—
I 2 study (%)	16.0	31.9	—	91.59
I 2 species (%)	—	10.3	—	—
I 2 location (%)	30.5	29.5	23.35	6.02
I 2 total (%)	87.1	71.8	78.67	97.62

TABLE 4. Summary of the moderators explaining the influence of elevation range (elvR) and latitude difference (latD) on the variation of life‐history and physiological traits.

Moderator	Adult body size			Adult female body size
	Est	SE	p	Est	SE	p
elvR	−0.006	0.058	0.915	−0.023	0.060	0.705
latD	−0.007	0.022	0.744	−0.007	0.028	0.795
elvR × latD	0.002	0.024	0.919	−0.023	0.028	0.414

Full text

Introduction

1

Life‐history strategies can be arranged along a fast‐slow continuum (Stearns 1977; Promislow and Harvey 1990; Alice Boyle et al. 2016). This continuum emerges from the selection pressure to maximise the number of surviving offspring over an individual's lifetime, while simultaneously optimising growth and survival given limited resources (Stearns 1992; Hille and Cooper 2015). As these conditions cannot all be perfectly satisfied concurrently, this constraint is hypothesized to result in a co‐variation of life‐history traits across taxa called the pace‐of‐life (POL), which defines the fast‐slow continuum (Ricklefs 2000). Recent developments have extended this concept, initially to incorporate physiological (Ricklefs and Wikelski 2002; Wikelski et al. 2003; Martin II et al. 2006) and, later, behavioural traits (Réale et al. 2010), which co‐evolve with life‐history traits in a framework called the pace of life syndrome (POLS) hypothesis. According to POLS, suites of traits co‐vary across fast‐slow gradients, with a slow pace of life being characterised by small clutches, slow growth, long lifespan, and slow metabolism, while animals with a fast pace of life would have traits such as large clutches, fast growth, short lifespan, and fast metabolism.

To date, the POLS hypothesis has been tested for a variety of taxa, most comprehensively for birds (e.g., Hille and Cooper (2015), Wiersma et al. (2007), Wikelski et al. (2003)), mammals (e.g., Bielby et al. (2007), Promislow and Harvey (1990), Oli (2004)) and fish (e.g., Binder et al. (2016), Polverino et al. (2018), Turko et al. (2019)). Healy et al. (2019) concluded that a large proportion of variation in animal life history strategies can be explained by the pace of life, mortality risk and the number of breeding episodes throughout an individual's lifetime (i.e., degree of iteroparity). However, Dammhahn et al. (2018), based on a topical collection of meta‐analyses and reviews on POLS syndromes, as well as other integrative studies (Royauté et al. 2018; Mathot et al. 2019), emphasise that the empirical evidence for the POLS framework is mixed, especially at the within‐species level.

Variation in ecological conditions can mediate the trade‐off between survival, growth and reproduction, resulting in the evolution and divergence of traits across populations experiencing distinct environments (Mathot and Frankenhuis 2018), and, in the end, to an optimal pace of life solution. This hypothesis has been most prominently investigated in relation to latitude and, to a lesser extent, elevation (Hodkinson 2005; Hille and Cooper 2015). Spatial variation in latitude and elevation produces similar conditions (i.e., lower temperatures, greater seasonality, higher stochasticity, lower productivity, shorter breeding seasons) (Hodkinson 2005; Camfield et al. 2010; Hille and Cooper 2015; Alice Boyle et al. 2016; Balasubramaniam and Rotenberry 2016), but they also differ in several important aspects: high elevation environments are subject to increased solar insolation and lower atmospheric pressure (including partial pressure of gases), while maintaining the same photoperiod as the lowlands below (Körner 2007; Alice Boyle et al. 2016; Balasubramaniam and Rotenberry 2016). Moreover, elevational gradients are characterised by steep environmental changes over short horizontal distances (usually kilometres), which also leads to a decrease in land area (Körner 2007), while latitudinal gradients result in more gradual change over longer horizontal distances (usually hundreds of kilometres) (Hodkinson 2005) and do not necessarily imply a decrease in environmental space.

Understanding how predictably traits associated with pace of life respond to environmental gradients such as elevation can lead to more accurate predictions of organismal responses to current and future environmental change (Martin 1996; Buckley 2010; Alice Boyle et al. 2016; Martin et al. 2023). Indeed, elevational gradients are commonly used as a space‐for‐time substitution to better predict ecological and evolutionary consequences of climate change (Verheyen et al. 2019; Lovell et al. 2023). While organisms show widespread adaptive variation in phenotypic traits (especially life‐history traits) corresponding to changes in elevation (Keller et al. 2013), whether elevation consistently drives predictable and parallel pace‐of‐life evolution within species is largely unknown and has not been quantitatively tested for any taxa outside of birds (Alice Boyle et al. 2016; Balasubramaniam and Rotenberry 2016).

With a legacy of hundreds of millions of years of evolution, extant squamates are the most diverse clade of terrestrial vertebrates (Uetz et al. 2024). Squamates occupy most areas of the Earth, except polar regions and display a vast diversity of morphological and ecological forms (Vitt and Caldwell 2014). Thus, not surprisingly, they also exhibit a wide range of life histories and physiological and behavioural adaptations, shaped and evolved by the environments they live in (Shine 2005). Within squamates, the pace of life was investigated most comprehensively on the garter snake ( Thamnophis elegans ) (see Bronikowski and Arnold (1999), Bronikowski (2000), Bronikowski and Vleck (2010), Gangloff et al. (2017), Palacios et al. (2013)). Healy et al. (2019) also included reptile species from most living groups (chelonians, squamates and crocodilians) in their analysis of pace‐of‐life evolution. The adaptive impact of environmental factors on the phenotypic traits of reptiles has been a major focus of many studies, reviews and meta‐analyses, notably in relation to body size (Bergmann's rule) (Ashton and Feldman 2003; Angilletta et al. 2004; Ashton 2004; Sears and Angilletta Jr. 2004; Olalla‐Tárraga et al. 2006; Pincheira‐Donoso and Meiri 2013), life‐history traits (Sinervo 1990a; Rohr 1997; Sears and Angilletta Jr. 2003; Mesquita et al. 2016), and thermal physiology and behaviour (Sinervo 1990b; Garcia‐Porta et al. 2019; Muñoz and Bodensteiner 2019; Rivera‐Rea et al. 2023). Since single‐species studies lack generality, while single‐trait studies might overlook important factors driving variation (Alice Boyle et al. 2016), currently missing are integrative studies investigating the variation of traits within the POLS framework in squamates and the role environmental factors such as elevation play in shaping the pace of life (but also see Scharf et al. (2015)).

The goal of our study is to provide a quantitative synthesis testing intraspecific divergence and variation in squamate life history and physiology along the pace of life continuum in response to elevation. More precisely, we compare trait means and trait variation between conspecific populations occurring at different elevations, rather than comparing traits across species. We note that our approach tests whether individual POLS‐related traits diverge in the predicted direction along elevational gradients, rather than testing for the existence of POLS, which would require examining covariation among multiple traits (Réale et al. 2010; Dammhahn et al. 2018). Further, we also examine how these responses differ among squamate taxa (suborders, families). We broadly predict that, overall, populations of squamates from low elevations will exhibit traits corresponding to a fast pace of life, while populations from high elevations will conform to a slow pace of life, given the associated changes in environmental conditions and selection pressures with increasing elevation, specifically colder temperatures, shorter available activity times, and decreases in food availability and predation leading to convergent evolutionary (and/or plastic) responses (Adolph and Porter 1993; Shine 2005; Laiolo and Obeso 2017).

We also investigated how trait variability varied across elevation. Here, we predicted that populations living at higher elevations will show increased trait variation overall compared to lower elevations, because of the reduced availability and greater spatio‐temporal variation in resources (e.g., in the thermal environment or in food resources), resulting in greater variation in physiology and growth among individuals (Shine 2005). For the traits with data available to test our predictions, we expected populations from higher elevation to show higher heterogeneity in body size and reproductive output. In contrast however, we predicted higher elevation populations to be more precise thermoregulators due to lower thermal quality of their environments, leading to a decreased variation in body temperature at higher elevations (Blouin‐Demers and Nadeau 2005; Lourdais et al. 2013; Lymburner and Blouin‐Demers 2020).

Methods

2

Data Collection

2.1

We conducted an extensive literature search up to the end of 2024 for empirical field studies reporting trait data (life‐history or physiology) for two or more populations of the same species located at different elevations. One person (TCS) was responsible for the entire process, from performing the search to constructing the final database of articles used for the meta‐analyses. The search was conducted on Google Scholar and Scopus using multiple combinations of search terms (see Appendix S3, Table S1). First we looked at the title and abstract to see if it matched our area of interest in any way, then quickly continued to the methodology section to check if at least two populations located at different elevations were used, and then we looked for any type of summary data or Supporting Information; material from which we could extract information. Google Scholar does not allow the bulk download of references so in this case the entire process was performed online. We avoided the inconsistent search results displayed in Google Scholar by having one person conduct the entire search procedure and using browser bookmarks to continue with the results from the same search. We collected peer‐reviewed articles published in English which contained relevant data, based on a rigorous protocol detailed in Appendix S3—Supporting Information; methods. We started with 1916 articles and in the end we extracted data from 192. When faced with multiple populations or summary statistics divided into different categories we used statistical techniques for pooling data (N = 131 effect sizes) (see Appendix S1—dataPooling column, Appendix S3—Supporting Information; methods). We also extracted data from figures and, in some cases, specifically when we lacked one study to include a family in a meta‐regression, we e‐mailed authors of the primary papers (we contacted the authors of four papers and received one reply with three effect sizes for Vipera ammodytes ).

Upon finalising the database, we assigned traits to several categories: (1) adult body size, composed of snout‐vent length (SVL) and mass of adults, (2) adult female body size, composed of SVL of adult females, irrespective of reproductive status, (3) neonate body size, composed of SVL and mass of neonate individuals (not juveniles), from either oviparous or viviparous species, born in nature, or from clutches incubated to mimic the natural environment of the focal population, (4) body temperature (T b), comprising measurements taken in the field, (5) lower critical temperature (CT_min_), measuring the critical thermal minimum, (6) upper critical temperature (CT_max_), measuring the critical thermal maximum, (7) egg size, comprising measurements of egg volume, egg mass, egg length and egg width, (8) mean age, comprising the estimated mean age of individuals obtained from growth rings and (9) fecundity, comprising measurements of clutch size, clutch mass, clutch weight, mean number of live born, mean litter mass and relative clutch mass (RCM) (Appendix S2).

Effect Size Calculation

2.2

Trait categories collected varied in their unit of measure, so we chose the response ratio (or ratio of means) (lnRR) as our effect size metric (Hedges et al. 1999) for the analysis of mean differences between low and high altitude populations. The response ratio has been widely used in ecological studies since its introduction (Lajeunesse 2015) and, similar to the standardised mean difference (SMD), it is unitless, thus allowing for combining different outcomes (Friedrich et al. 2008). Unlike SMD, though, the response ratio has no indication of the magnitude of the effect, because it is calculated as the natural logarithm of the ratios between the two means (Nakagawa et al. 2015), which can be beneficial in cases when there is no knowledge regarding the expected magnitude of change. However, a sense of magnitude can be obtained by converting the result to a percentage change using the formula, where Ψ is the effect size (Pustejovsky 2018). Since the result is symmetric around 0, in our analyses negative values indicate a change towards low elevation populations, while positive values indicate a greater response in high elevation populations.

For the analysis of variation of traits between low and high elevation populations, we used the coefficient of variation (lnCVR) as our effect size statistic (Nakagawa et al. 2015), which is calculated as the natural logarithm of the ratio between the coefficients of variation from two groups and does not necessitate any additional data beyond what is used for the analysis of response ratios. Moreover, the coefficient of variation is most useful in systems where a mean–variance relationship (such as Taylor's law) exists (Senior et al. 2020) such as populations in higher elevation environments where total area decreases. Since the results are also symmetric around zero, we changed the sign of the results to match the meta‐analysis of response ratios: negative values signify greater variation of traits in low elevation populations, while positive values signify greater variation in high elevation populations.

Effect sizes were computed using the escalc function in the metafor package (Viechtbauer 2010). We calculated a separate effect size for each summary statistic when multiple summaries were derived for the same population (e.g., body mass and snout‐vent length) to maximise the amount of data available (Noble et al. 2017).

Meta‐Analytical Methods

2.3

Meta‐analyses and meta‐regressions were performed in R v4.4.1 (R Core Team 2022) using the metafor package with a restricted maximum likelihood estimator (REML). Datasets for meta‐analyses and meta‐regressions, both for the mean effect size and the coefficient of variation, were selected based on a set of rules detailed in the Supporting Information; methods. We used three distinct model types to address our research questions regarding mean effect and variation of traits: (Adolph and Porter 1993) overall meta‐analyses, which we used to test for intraspecific differences in mean and variability of squamate traits between low and high elevation populations (Agrawal 2020) categorical moderator meta‐regressions, where we compare intraspecific differences in trait response to elevation in higher taxonomic units, using suborder or family as categorical moderators, and (Alice Boyle et al. 2016) continuous moderator meta‐regressions, which tests the influence of elevational and latitudinal factors in determining the observed intraspecific differences.

We used multi‐level models to deal with shared measurements and phylogenetic relationships between species, and random factors were included based on their ability to explain the total amount of heterogeneity (I ^ 2 ^: representing the total relative within‐and‐between study variance) (factors explaining < 1% of heterogeneity were dropped), calculated for multi‐level models using the orchard 2.0 package (Nakagawa, Lagisz, et al. 2023). Phylogeny was controlled using the most recent phylogenetic subset downloaded from VertLife (https://vertlife.org) and the procedure detailed in Supporting Information; methods, as described by Noble et al. (2017), Nakagawa et al. (2022), and Nakagawa, Yang, et al. (2023).

Models used the following structure for overall meta‐analyses zi=β0+ueffecti+ustudyi+uspeciesi+uphylogenyi+ulocationi+mi, where z _ i _ is the effect size for the observation i, β _ 0 _ is the overall mean effect size (or intercept), u effect to u location are the random effects, and m _ i _ is the sampling error.

For the meta‐regressions with categorical moderators, the model structure was the following for suborders and families: zi=β1X1i+…+βkXki+ueffecti+ustudyi+uspeciesi+ uphylogenyi+ulocationi+mi, where β _ 1 _ to β _ x _ are the mean effects corresponding to dummy variables X _ 1 _ to X _ k _ for each family or suborder; the categorical models did not have an intercept as we were interested in differences in higher taxonomic levels not comparisons to a reference level.

For meta‐regressions with continuous moderators model assumed the following structure: zi=β0+β1X1+…+βkXk+ +ueffecti+ustudyi+uspeciesi+uphylogenyi+ulocationi+mi, where β _ 1 _ to β _ k _ are the regression coefficients for the variables X _ 1 _ to X _ k _. We note that the inclusion of the random factors in the final models was based on the amount of heterogeneity explained (see above and Tables 1 and 3).

**TABLE 1: Mean response ratios (β 0) modelled for each trait using multi‐level meta‐analyses without moderators. For each trait, the number of studies, number of effect sizes (k), standard errors (SE), 95% confidence intervals, significance (p) and heterogeneity (I

2 ) are summarised. In cells lacking I

2 values, the corresponding random factor was not included in the model based on the percentage of heterogeneity explained (see Methods section).**

We tested for the impact of shared measurements using a variance–covariance matrix as suggested by Nakagawa, Yang, et al. 2023 and continued the analyses with the original values based on the results (see methods).

Moderating Variables

2.4

In our meta‐regression models, we used elevation range (continuous), latitudinal difference (continuous), mean latitude (continuous), suborder (categorical) and family (categorical) as moderators (analogous to predictor variables in linear models) to test different hypotheses. Meta‐regressions were produced when there were at least four studies for each family and a minimum of 20 total effect sizes (Fu et al. 2011), and to keep results consistent, we applied the same criteria for all meta‐regressions, independent of moderator type. The trait categories that fit our criteria were adult body size, body temperature, adult female body size, neonate body size, and fecundity for mean effects, and adult body size, body temperature, adult female body size, and fecundity for the coefficient of variation.

Data regarding the elevation of the populations was extracted directly from the study (if available); otherwise, we used geographic coordinates or georeferencing to extract the data (see Supporting Information; Methods). Elevation of the population pairs was used to calculate elevation range, which we then log‐transformed to improve normality and included as a moderator to test if an increased separation corresponds to an increase in trait differentiation.

Because our database included all multi‐population studies which contained relevant data, it also meant that some populations had large latitudinal spacing which could create a bias in the results. We accounted for the effect of latitudinal differences between populations and between studies using two variables: (1) the log‐transformed absolute value of the latitude difference between low and high elevation populations, and (2) the absolute value of the mean latitude of the population pairs (see Supporting Information; Methods).

Publication Bias

2.5

We performed publication bias tests in the form of funnel plots and meta‐regressions aimed at detecting (a) a small study effect, which refers to cases where effect sizes based on small sample sizes tend to be larger (using Egger's regression adapted for multilevel models) or (b) a time‐lag bias or decline effect, where positive results are published earlier than negative results (Nakagawa et al. 2022; Nakagawa, Yang, et al. 2023). We performed publication bias tests for traits where meta‐regressions were also performed (see Supporting Information; methods). We detected a time‐lag bias for response ratios of neonate body size in both single moderator meta‐regression (t = −2.379, p = 0.0265) and the all‐in multilevel meta‐regression (t = −2.2342, p = 0.0365). Therefore, the results presented and discussed for neonate body size are the bias‐corrected values of the intercept (β 0) (Nakagawa, Yang, et al. 2023).

Results

3

Our final database used for the meta‐analyses consisted of 671 effect sizes extracted from 192 studies, comprising 104 species of squamates from 10 families (Appendix S2). Families Lacertidae and Phrynosomatidae had the greatest number of species (33 and 18 species respectively), while for Anguidae and Phyllodactylidae we only had one representative each. We collected data for six life‐history traits—adult body size, 250 effect sizes extracted from 140 studies, adult female body size, 106 effect sizes extracted from 93 studies, neonate body size, 45 effect sizes extracted from 26 studies, egg size, 16 effect sizes from 11 studies, mean age, 18 effect sizes from 18 studies, fecundity, 130 effect sizes from 63 studies and three physiological traits—body temperature (T b), 46 effect sizes from 40 studies, lower critical temperature (CT_min_), 26 effect sizes from 14 studies and upper critical temperature (CT_max_), 35 effect sizes from 21 studies. Coordinate pairs were available from around the world, including North, Central and South America, Eurasia, Australia and Africa, but from the latter we recovered the least amount of data (Appendix S3, Figure S1). The median difference in elevation between population pairs was 918 m (range 100–3500 m) and the highest elevation population was located at 4583 m.

Trait Divergence Across the Elevational Gradient

3.1

Overall Effect Size

3.1.1

Overall, our results do not show patterns of intraspecific trait divergence consistent with a POLS within squamate species in response to elevation. With few exceptions, effect sizes were small and non‐significant (Figure 1). Egg size increased significantly towards higher elevations by an average of 11% [CI: 0.7, 22.87] as predicted, while body temperature declined with elevation by an average of 5.3% [CI: 3.45, 7.25], also matching our expectations. On the other hand, we did not find significant differences between low and high elevation populations in adult body size (mean: 3.04% [CI: −2.83, 9.3]), adult female body size (mean: 0.6% [CI: −2.12, 3.25]), neonate body size (mean: 3.35% [CI: −1.71, 8.65]), fecundity (mean: −5.44% [CI: −22.87, 10.51]), mean age (mean: 8% [CI: −13.76, 32.84]), lower critical temperature (mean: −12.18% [CI: −35.93, 7.64]) or upper critical temperature (mean: −0.6% [CI: −2.22, 1]).

Summary of mean responses (log response ratio, lnRR) and 95% confidence intervals for life‐history and physiological traits of squamates in response to elevational differences. Individual effect sizes are displayed as points scaled by precision (1/SE), where larger points indicate more precise estimates. The mean effect is shown as the central, outlined point and the black whiskers provide the 95% CI For each trait, k indicates the number of effect sizes, with the number of studies in parentheses.

Heterogeneity (i.e., relative within‐and‐between study variance) was very high in all intercept models with I ^ 2 ^ values above 90%, indicative of a high degree of unexplained variation. Heterogeneity was a result of differences in measurements (I ^2^ effect), methodological differences (I ^2^ study), geographical location of the studies (I ^2^ location), and in some cases variability as a result of phylogenetic relatedness (I ^2^ phylo) or differences between species (I ^2^ species) (Table 1).

Divergence Among Suborders and Families

3.1.2

We found some divergence in trait responses between the two suborders (Lacertilia and Serpentes) and for some families, but differences were generally small and non‐significant. At the suborder level, only snakes showed significant differences between low and high elevation populations with increased adult female body sizes (mean: −6.82% [CI: −13.31, −0.7]) and fecundity (mean: −67.02% [CI: −110.22, −32.57]) at lower elevations, and bigger neonate body sizes at higher elevations (mean: 14.91% [CI: 2.12, 29.3]) (Appendix S3, Table S3, Figures S5 and S6, S8). However, we point out the sample sizes for snakes, which are considerably smaller for snakes than for lizards (Appendix S3, Figures S5 and S6, S8), therefore increasing the probability of false positives. The pattern identified for suborders was mostly mirrored at the family level. Notably, neonate body size increased significantly in high elevation populations of colubrids (mean: 14.56% [CI: 1.81, 28.91]), while fecundity was significantly reduced (mean: −66.36% [CI: −113.4, −29.69]) (Appendix S3, Table S2, Figures S6 and S7). Additionally, body temperature was significantly lower in high elevation populations of Anolidae (mean: −9.85% [CI: −16.41, −3.66]) and Phrynosomatidae (mean: −5.75% [CI: −10.18, −1.51]) (Appendix S3, Table S2, Figure S9).

Effect of Moderators

3.1.3

Elevation range was not a significant moderator of body size differences in any of the life stages analysed (adults, neonates), but differences in adult female body size were positively and significantly associated with increased latitude differences (Table 2). Neonate body size was also significantly moderated by mean latitude (increasing divergence at higher mean latitude), but only in analyses with a non‐significant interaction with elevation range, and not in the single moderator model (Appendix S3, Table S8).

We found a significant interaction between elevational range and latitudinal difference on differences in fecundity (Table 2). Over much of the parameter space, low elevation populations tended to have greater fecundity than high elevation populations, and this effect increased with increasing elevational range. In contrast, when populations were separated by large latitudinal differences, higher elevation populations tended to have greater fecundity (Appendix S3, Figure S10). However, there were few studies where populations were separated by large differences in latitude (Appendix S3, Figure S11), and so the biological interpretation of this interaction should be treated with some caution.

We then dropped the interaction to examine a reduced model with latitude difference and elevation range as main effects. In this model, elevation was marginally non‐significant (β = −0.082, SE = 0.041, t = −1.991, p = 0.052) with a trend of increasingly greater fecundity for low elevations when controlling for latitudinal differences. Latitudinal difference, in contrast, was significant, with higher elevation populations showing increased fecundity with increasing latitudinal differences (β = 0.037, SE = 0.017, t = 2.181, p = 0.034). Overall, and keeping in mind the smaller number of studies with large latitudinal extents, this suggests that latitude might moderate divergence in fecundity across elevation.

The body temperature difference between low and high‐altitude populations was significantly and negatively influenced by the difference in elevation (Table 2) in analyses with latitude difference but interacted with mean latitude of the population pair when we included it as a moderator (Appendix S3, Table S8). Body temperatures are greater for low elevation populations across almost all combinations of elevational range and mean latitude, whereas higher body temperatures for high elevation populations are only found at high values of mean latitude combined with small elevational differences between low and high populations (Appendix S3, Figure S12).

Trait Variation Along the Elevational Gradient

3.2

Overall Effect Size

3.2.1

The results on variation of life‐history and physiology traits completely contradicted our predictions. Body temperature variation was significantly higher in high elevation populations of squamates, on average by 23.24% [CI: 3.97, 46.08] (Figure 2, Table 3), contrary to our expectation of more precise thermoregulators at higher elevations. Moreover, also opposed to our expectations, there was no significant variation across elevational gradients in adult body size (mean: 3.45% [CI: −3.77, 11.18]), adult female body size (mean: 3.97% [CI: −4.5, 12.97]), or fecundity (mean: 3.66% [CI: −4.7, 12.52]) (Figure 2, Table 3).

Summary of mean variation (coefficient of variation, lnCVR) and 95% confidence intervals for life‐history and physiological traits of squamates in response to elevational differences. Individual effect sizes are displayed as points scaled by precision (1/SE), where larger points indicate more precise estimates. The mean effect is shown as the central, outlined point and the black whiskers provide the 95% CI For each trait, k indicates the number of effect sizes, with the number of studies in parentheses.

**TABLE 3: The coefficient of variation (β 0) modelled for each trait using multi‐level meta‐analyses without moderators. For each trait, the number of studies, number of effect sizes (k), standard errors (SE), 95% confidence intervals and significance (p) and heterogeneity (I

2 ) are summarised. In cells lacking I

2 values, the corresponding random factor was not included in the model based on the percentage of heterogeneity explained (see Methods section).**

Heterogeneity was high in the intercept models for adult body size (I ^2^ total = 87.1%), fecundity (I ^2^ _total _ = 78.67%) and body temperature (I ^2^ total = 97.62%) and moderate for adult female body size (I ^2^ total = 71.8%) (Table 3). Variability came mainly from differences in measurements (I ^2^ effect), the methodology used for the studies (I ^2^ study), the location of the studies (I ^2^ location), or was the result of differences between species (I ^2^ species) (Table 3).

Divergence Among Suborders and Families

3.2.2

There were no differences between low and high elevation populations in the two suborders (Appendix S3, Table S5). There was also no significant variation between low and high elevation populations in families for most of the traits analysed, with the exception of adult female body size, where viper populations from higher elevations showed significantly more variation compared to low elevations (mean: 63.88% [CI: 3.97, 158.05]) (Appendix S3, Table S4, Figure S4).

Effect of Moderators

3.2.3

Elevation difference played no significant role in predicting the variation of life‐history or physiology traits (Table 4), nor did latitude difference or mean latitude of the population pairs (Table 4, Appendix S3, Table S9).

Discussion

4

Our comprehensive synthesis found no clear evidence for intraspecific elevational trait divergence in life history or physiology consistent with the existence of an overall pace of life syndrome in squamate reptiles. Despite analysing data from 192 studies encompassing 104 species and 9 life‐history and physiological trait categories, we found consistently small and non‐significant effect sizes. This preponderance of negative results challenges the overall applicability of POLS theory within species in the context of environmental gradients.

This is the first comprehensive test of whether individual traits within squamate species diverge along elevational gradients in directions predicted by POLS, providing the first such data for squamates. These negative findings align with evidence that the detectability of POLS‐related trait shifts in response to environmental gradients varies across contexts. While our study tested for trait divergence in the direction predicted by POLS theory (e.g., slower pace of life at higher elevation) and not trait covariation per se, the lack of consistent directional divergence in individual traits suggests that elevation does not reliably generate the selection pressures assumed to drive POLS‐related trait shifts within species. Recent meta‐analyses across various taxa and environmental gradients have similarly shown that the detectability of POLS patterns is highly context dependent, varying with factors such as environmental conditions, sex or species studied (Santostefano et al. 2017; Royauté et al. 2018; Tarka et al. 2018). Our results thus contribute to understanding the conditions under which environmental gradients drive life‐history divergence and highlight the need for more nuanced approaches to understanding life‐history evolution. Indeed, there are scant formal models for predicting trait covariation under POLS (Mathot and Frankenhuis 2018; Laskowski et al. 2021) and there is ambiguity in which traits co‐vary as part of a syndrome or if composite variables should be considered (Araya‐Ajoy et al. 2018; Mathot and Frankenhuis 2018; Royauté et al. 2018; Gopal et al. 2023). Our own synthesis tested predictions derived from a verbal formulation of a POLS hypothesis based on the trade‐offs expected due to varying environmental conditions generated by the elevational gradient.

There are several potential reasons for the apparent lack of a unified elevational POLS within squamate reptile populations that likely reflect the extraordinary ecological and evolutionary diversity within this group. Although squamates are monophyletic and thus share a common evolutionary history, they exhibit a tremendous diversity of shape, function, and ecology that could yield disparate solutions to the shared challenge of adapting to elevation. Indeed, our suborder‐level analyses reveal that lizards and snakes respond fundamentally differently to elevation: snakes showed significant differences in adult female body size, neonate body size, and fecundity between low and high elevation populations, while lizards showed minimal and non‐significant responses across all life‐history traits. This divergence likely reflects the greater ecological and morphological disparity among lizards compared to snakes. Lizards have a great diversity of body shapes and diets, ranging from carnivory to herbivory, even within families, and some groups show increased thermal plasticity through behavioural thermoregulation (Labra et al. 2009; Muñoz and Losos 2018; Refsnider et al. 2018; Plasman et al. 2025). In contrast, snakes are limited in body shape, feeding ecology (all known species are predators) and rely on adjusting the timing of activity for coping with thermal environments (Weatherhead et al. 2012). Shine (1980) suggests that in lizards, reducing current fecundity for future gains in body size offers little benefit, as survivorship is likely the main reproductive cost. In snakes, on the other hand, adult female body size may decrease at higher elevations through a combination of reduced productivity (e.g., Alice Boyle et al. (2016)), restricted trophic breadth and reduced female mobility during pregnancy (Seigel et al. 1987).

Second, environmental complexity might moderate trait covariation. Multivariate environmental variation is often cited as a potential moderator of trait co‐variation, as variation in ecological conditions can induce or mask changes in morphology, life‐history, or physiology, with distinct effects across taxa and life stages (Hämäläinen et al. 2021). While elevation itself does not influence the life‐history and physiology of animals per se, it acts as a proxy for an aggregate of factors (such as temperature, seasonality, precipitation quantity and type, primary productivity, habitat size, predation risk) which affect ectotherms and endotherms differently (reviewed extensively in Alice Boyle et al. (2016), Hille and Cooper (2015), Hodkinson (2005) and Laiolo and Obeso (2017)). This environmental complexity means that elevation creates multiple, potentially conflicting selection pressures rather than a single coherent selective environment that would favour directional trait evolution.

Our results demonstrate this complexity empirically. The elevation‐latitude interaction we detected for fecundity exemplifies how multiple environmental gradients can interact. Our prediction that populations from lower elevations would exhibit greater fecundity, characteristic of a fast pace of life, did not match our results—with the exception of snakes. Although there was a trend towards increased fecundity at lower elevations, the effect was minimal and non‐significant. Instead, fecundity appears to be shaped by an interaction between elevational and latitudinal distance between populations. Similarly, divergence in body temperature was influenced by an interaction between elevational difference and mean latitude. That patterns of trait divergence associated with elevation may interact with latitude is perhaps not surprising as seasonality and annual temperature variation are much greater at higher compared to lower latitudes (Janzen 1967; Ghalambor et al. 2006). Moreover, these non‐linear patterns suggest regional adaptation, consistent with the interpretation that large‐scale ecological variation promotes the evolution of locally adapted traits rather than universal elevational response (Wang and Bradburd 2014).

Gene flow, especially along steep elevational gradients, could also constrain local adaptation (Bachmann et al. (2020); reviewed in Lenormand (2002)). Indeed, Bachmann et al. (2020) found that adaptive divergence in developmental rates in anurans was much weaker—and conversely local maladaptation was stronger—across steeper elevational clines compared to shallower latitudinal clines of comparable magnitude. Our results also show a lack of intraspecific divergence in thermal tolerance traits (CT_min_ and CT_max_), which is consistent with the pattern of gene flow constraining local adaptation. However, recent research (Lange et al. 2025) found that lower and upper critical thermal temperatures predict the cold and warm distribution limits of lizards in the genus Anolis along elevational gradients in the Caribbean. This contrast suggests that while there is a thermal physiology‐environment relationship at the species level, which is not subject to the homogenising effect of gene flow, this pattern will not be detected in intraspecific comparisons due to ongoing gene flow between populations at different elevations which prevents local adaptation.

Behaviour could also act to shield morphological and life history traits from selection in response to elevation (i.e., the Bogert effect Huey et al. (2003); Muñoz (2022)). For example, active thermoregulatory behaviour—especially widespread in lizards—could forestall evolutionary divergence of traits such as thermal tolerance across elevations (Huey et al. 2003; Muñoz and Bodensteiner 2019). Parental care, such as nest‐site selection (Shine and Harlow 1993; Shine and Harlow 1996; Telemeco et al. 2009) or maternal basking behaviour (Shine and Harlow 1993) could similarly influence their offspring's early environment and partially shield their phenotypes from selection (Diamond and Martin 2021).

Additionally, variation in resource acquisition, an important driver of life‐history evolution, but one which we did not account for in our analysis due to lack of data, could further obscure POLS patterns across elevation. Resource availability and competition for resources differ dramatically between elevations because lower elevations support higher productivity but also exhibit higher levels of competition (Morin and Chuine 2006; McCain 2010; Alice Boyle et al. 2016). If increased resource acquisition at low elevations is offset by intensified competition and, inversely, reduced competition at high elevations compensates for lower resource availability, this could result in shifts of life‐history strategies, masking or counteracting POLS predictions. Indeed, as variation of life‐history traits is fundamentally linked to resource acquisition rate (van Noordwijk and de Jong 1986), the lack of information regarding resource dynamics across the elevation gradient may contribute to the absence of a detectable syndrome in our analysis.

There were a few traits with responses that significantly differed across elevation both in support of and in contrast to our prior predictions. First, body temperature responses matched our predictions, with all families for which we had data (all lizards) exhibiting lower body temperatures at higher elevations and divergence increasing with greater elevational differences. This pattern recalls the Thermal Behaviour Syndrome (TBS) hypothesis proposed by Goulet et al. (2017), which argues that, given the strong and positive relationship between metabolism and temperature, ectotherms would be situated along a cold‐hot continuum, where individuals would select and perform better at low body temperatures (cold types) or at high body temperatures (hot types) (Goulet et al. 2017). Recent empirical work by Michelangeli et al. (2018) found evidence for a thermal‐behavioural syndrome in delicate skinks ( Lampropholis delicata ) and comprehensively demonstrated that thermal type also predicted habitat use, with ‘hot’ types spending more time in warmer, more exposed microhabitats, suggesting that thermal physiology could drive ecological niche partitioning within populations. Unfortunately, our dataset lacks behavioural or metabolic data needed to test this hypothesis. However, the potential for within‐population thermal divergence warrants further investigation, as individuals with different thermal types may respond distinctively to elevational selection pressures. The emergence of thermal patterns despite the absence of life‐history POLS reinforces the argument that the lack of syndrome‐level covariation reflects genuine biological differences rather than methodological limitations.

Our expectations regarding body temperature variation were that squamates would be precise thermoregulators (Díaz and Cabezas‐Díaz 2004; Blouin‐Demers and Nadeau 2005), therefore showing a reduced level of variability at higher elevations, owing to the strong dependence of ectotherm performance with body temperature (Huey and Kingsolver 1993). Instead, our results show that the mean variation of body temperature increases at higher elevations, which would suggest there are increased costs to thermoregulation at higher elevations (Huey and Slatkin 1976), probably related to reduced basking opportunities and the heterogenous and cooler nature of high‐altitude habitats. Alternatively, some squamate species can shift between thermoregulation and thermoconforming depending on ecological conditions (Hertz et al. 1993; Dubiner et al. 2024) and this behaviour flexibility could contribute to increased body temperature variation at high elevations. However, since many factors that could lead to increased variation were not controlled for in this body temperature data (e.g., season, time of day, sex, age, degree of behavioural thermoregulation), we caution that this result should not be overgeneralized. Moreover, increased variation in body temperature could instead be due to the evolution of wider thermal performance curves at higher elevations (Tüzün and Stoks 2018).

Lastly, as predicted, we did find overall support for increased egg size at higher elevations, and similarly an increase in neonate size within snakes at higher elevations. These results, together with the significant declines of fecundity at high versus low elevation in snakes, suggest a trade‐off between offspring investment and fecundity. Several factors may explain this pattern. First, elevation may alter age specific mortality—environments with low adult mortality relative to juvenile mortality are expected to favour larger eggs and smaller clutches and vice versa (Stearns 1992). The colder, shorter growing seasons and lower food availability of high elevations along with higher predation at lower elevations could drive patterns of higher juvenile vs. adult mortality at high elevations and higher adult vs. juvenile mortality at low elevations (Hille and Cooper 2015; Laiolo and Obeso 2017). Supporting this, longevity tended to increase with elevation, although not significantly so. Second, selection for viviparity could also lead to reduced clutch sizes, larger eggs and larger neonates at higher elevations. Viviparous species and populations are more common at higher elevations than oviparous ones, and high elevation oviparous populations tend to retain eggs longer before laying than low elevation populations (Mathies and Andrews 1995; Rodríguez‐Díaz and Braña 2012). The advantages of viviparity and longer egg retention are likely due to the benefits of maternal thermoregulation for embryo development in thermally challenging environments (Packard et al. 1977; Tinkle and Gibbons 1977; Shine 1985). In turn, longer egg retention increases clutch mass via water retention (Vleck 1991), which may favour smaller clutches, reducing negative effects on female mobility (Sinervo et al. 1991). Of course, despite these significant patterns, empirical counter examples exist, with the reverse trade‐off occurring across elevation (e.g., Sinervo (1990a)).

The preponderance of negative results has important theoretical implications that extend beyond squamate biology. The absence of intraspecific elevational POLS challenges fundamental assumptions about how predictable environmental gradients should generate predictable trait syndromes within species. Some authors argue that consistent pace of life is easier to detect across species than within species, as evolution has had more time to produce clear trade‐offs (Agrawal 2020). However, the more immediate constraint on intraspecific divergence may be gene flow, as ongoing genetic exchange between populations at different elevations can homogenise trait variation even under divergent selection (Lenormand 2002). Our results are consistent with this interpretation—we compared populations across a broad taxonomic sample encompassing 104 species and 10 families, yet found no consistent elevational patterns. Interspecific comparisons, which are generally not affected by gene flow and also benefit from the passage of time, may still reveal POLS signals that intraspecific tests cannot detect. Still other authors note that pace of life syndrome differences can be observed between individuals of the same species and between sexes (Dammhahn et al. 2018; Hämäläinen et al. 2018; Royauté et al. 2018; Arnqvist et al. 2022), suggesting that the relevant scale for detecting POLS may be narrower than environmental gradients. The formal development of POLS theory has lagged behind empirical applications (Mathot and Frankenhuis 2018; Laskowski et al. 2021), and our results underscore the need for more nuanced theoretical frameworks that account for phylogenetic constraints, ecological diversity, and context‐dependency in trait evolution.

Limitations

5

As part of our literature survey, we tried to assemble as broad a dataset as possible from various sources which were not always related to our own goal; therefore, some aspects may limit the generality of our findings. First, and the most obvious one, is the fact that we did not find enough data on some key traits such as growth rate, survival, metabolic rate, and our dataset does not include behavioural traits. Second, we were limited in our scope by the amount of data we could accrue, since we were looking for low‐high population pairs; as such we could only include 104 species belonging to 10 families, while as of January 2025 there are 11,991 named species of squamates in 76 families. Consequently, our results might not accurately capture patterns of trait divergence if our sample of studies overall or within taxonomic units is unrepresentative. Third, several geographic regions are also underrepresented in our study, and this includes most of the continents of South America, Africa and Oceania. Finally, our results were limited by selecting only studies published in English, which can lead to biases in overall results (Konno et al. 2020) as previous studies have shown that significant portions of literature related to biodiversity were not published in English (Amano et al. 2016).

Conclusions

6

Our meta‐analysis suggests that pace of life syndromes do not provide a useful framework for understanding intraspecific squamate responses to elevational gradients, although an elevational POLS signal may still be present in interspecific comparisons, which were beyond the scope of the current study. This negative result is both novel and important, contributing to growing evidence questioning the broad applicability of POLS theory while revealing that life‐history evolution is more complex and context‐dependent than syndrome‐based approaches suggest. Rather than representing a failure to detect existing patterns, our findings reveal fundamental limitations in applying universal syndrome concepts across diverse ecological and evolutionary contexts. Future research should focus on understanding the taxon‐specific and context‐dependent mechanisms that govern trait evolution, embracing the complexity and diversity of life history strategies rather than seeking universal syndromes.

Author Contributions

Tiberiu C. Sahlean: conceptualization, data collection, formal analysis, methodology, validation, visualisation, writing – original draft, writing – review and editing, funding acquisition. Ryan A. Martin: conceptualization, formal analysis, methodology, validation, visualisation, writing – original draft, writing – review and editing.

Funding

This work was supported by Academia Româna, RO1567‐IBB09/2026.

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Data S1: ele70343‐sup‐0001‐Supinfo.zip.