Exceptional parallelisms characterize the evolutionary transition to live birth in phrynosomatid lizards

Ethics statement

The data collection and experiments were conducted in accordance with the collecting permits (SGPA/DGVS/07946/08, 03369/12, 00228/13, 07587/13, 01629/16, 01205/17, 02490/17, 06768/17, 000998/18, 002463/18, 002490/18, 002491/18, 003209/18, and 02523/19) approved by Dirección General de Vida Silvestre, México.

Phylogeny and divergence time estimation

To estimate the phylogeny and divergence time among phrynosomatid species we used sequences of five mitochondrial and eight nuclear genes available in GenBank for 149 taxa (Supplementary Data 2). Accession numbers were the same as those used in Martínez-Méndez et al.⁵⁸ for the Sceloporus torquatus, S. poinsettii and S. megalepidurus groups and the same as those in Wiens et al.⁵⁹ for other phrynosomatid species. For taxa not included in the previous references, we searched GenBank for available sequences. We then performed alignments for each gene using MAFFT (ver. 7)⁶⁰ and concatenation and manual refinement using Mesquite (ver. 3.6);⁶¹ obtaining a concatenated matrix of 9837 bp for 149 taxa (Supplementary Data 3). For the relaxed clock analyses, three nodes were calibrated using lognormal distributions based on two previous studies^59,62. The first calibration was set for the Sceloporus clade (offset 15.97 million years ago (MYA)) based on a fossil Sceloporus specimen⁶³). The second calibration point was set for the Phrynosoma clade (offset 33.3 MYA) based on the fossil Paraphrynosoma greeni⁶⁴, and the last calibration point was for the Holbrookia–Cophosaurus stem group (offset 15.97 MYA) given the fossil Holbrookia antiqua⁶³. We conducted dating analysis with the concatenated sequences matrix, partitioned the mitochondrial and nuclear information, each gene under GTR + I + Γ model, and allowed independent parameter estimation. We performed Bayesian age estimation with the uncorrelated lognormal relaxed clock (UCLN) model in BEAST (ver. 2.5.2)^65,66 and run on CIPRES⁶⁷. Tree prior (evolutionary model) was under the Birth-Death model, and we ran two MCMC analyses for 100 million generations each and stored every 20,000 generations. We assessed the convergence and stationarity of chains from the posterior distribution using Tracer (ver. 1.7)⁶⁸. We combined independent runs using LogCombiner (ver. 2.5.2; BEAST distribution)⁶⁹ and discarded 30% of samples as burn-in, obtaining values of effective sample size (ESS) greater than 200. We estimated the maximum clade credibility tree from all post-burnin trees using TreeAnnotator (ver. 1.8.4)⁶⁹. The ultrametric tree is available as Supplementary Data 4. As we describe below, we accounted for phylogenetic uncertainty in our models by reperforming analyses using 500 trees that we randomly sampled from our posterior distribution. The 500 sampled trees are available as Supplementary Data 5.

Data collection

Parity mode

We categorized each species as either oviparous or viviparous based on previously published databases^21,37,51,70, published references, and unpublished data (Supplementary Data 1). Our assignations align with other studies, except for one species, Sceloporus goldmani, which has been previously considered a viviparous species^21,71,72,73. The only available sequence in GenBank (U88290) for that species is from a male (MZFC-05458) collected in Coahuila, Mexico⁷². However, in that same locality, one of us (F. R. Méndez-de la Cruz; unpubl. data) collected two females of the same species, and both laid eggs. Thus, the population of S. goldmani herein included is considered oviparous. Considering S. goldmani viviparous increases the number of originations of viviparity to 6 (from 5) in this lineage (Supplementary Fig. 4), but does not alter the outcome of our model-fitting analyses of trait evolution (Supplementary Table 7).

Thermal physiology

We compiled a database of four thermal physiological traits that influence the performance and fitness of ectotherms⁷⁴ for 104 phrynosomatid species. These data were gathered from both published sources and from our own field and laboratory work (Supplementary Data 1). The thermal physiological traits we examined were the field body temperature (T_b) of active lizards, the preferred body temperature (T_pref) in a laboratory thermal gradient⁷⁵, cold tolerance (critical thermal minimum, CT_min), and heat tolerance (critical thermal maximum, CT_max). These latter two traits (CT_min and CT_max) describe the thermal limits of locomotion; specifically, they describe the lower and upper temperatures, respectively, at which lizards fail to right themselves when flipped onto their backs^55,76. To minimize the confounding effects of experimental design, we limited our data selection to species that were measured with similar methods. Correspondingly, our new data collection approach mirrored that of the published studies from which we extracted data. To obtain mean values for each thermal physiological trait (CT_min, T_b, T_pref, and CT_max) we did not mix data measured from different locations (instead, we used data from the population with the highest sample size).

For species that we newly measured thermal physiological traits, we obtained the data as we describe below, and we based our methodology on the previous work^55,56,75,76. We captured active (perching) adult lizards by lasso or by hand, and immediately (<10 s) we measured their body (cloacal) temperature using a thermocouple (type K) connected to a digital quick-reading thermometer. We transported lizards to a field laboratory (which was at an ambient temperature of ~20 °C), measured the SVL and body mass of individuals, and the next day we measured their preferred body temperature by placing them into a laboratory thermal gradient from 8:00 to 17:00 h. The laboratory thermal gradient consisted of a wooden box (100 cm wide, 100 cm long, and 30 cm tall) divided into ten tracks. At one extreme of the laboratory thermal gradient, we placed ten 75 W bulbs (one per track) at the height of 25 cm above to generate a thermal gradient ranging from ~50 °C in the hot extreme to ~20 °C in the cold extreme. Then, we placed each lizard on a track, and we measured their cloacal temperature (using the same thermocouple and thermometer as in the fieldwork) every hour, during the length of the experiment. When we finished the experiment on thermal preferences, we performed the experiment on heat tolerance. For that, we placed individually each lizard into a plastic container (25 cm diameter and 30 cm height) and we increased their body temperature (1 °C/m) using a 90 W bulb suspended ~40 cm above the container. When lizards initiated panting behavior, we began flipping them onto their back every 20 s, and we recorded CT_max as the body temperature at which a lizard lost the ability to right itself. On the next day, we each placed individually each lizard into a plastic container (23 cm wide, 16 cm long, and 8 cm tall), and the plastic container on a bed of ice. Every 20 s we began flipping lizards onto their back, and we recorded CT_min as the body temperature at which a lizard lost the ability to right itself. For CT_max, and CT_min experiments, we did not include pregnant/gravid females, and after laboratory experiments, we hydrated lizards ad libitum, and released them at their capture sites.

In Supplementary data 1 we indicate coordinates where thermal physiological traits were measured and the environmental variables associated with each coordinate. In cases where locality details, but not coordinates, were available, we georeferenced sampling sites using Google Earth Pro (Version 7.3.3). All physiological data correspond only to adult lizards. Some studies have found that pregnant females reduce their core temperature to better match the optimal incubation temperature for their offspring^41,43,48. When we detected effects of reproductive condition on T_b or T_pref, we excluded data from pregnant (or gravid) females. To test whether behavioral and physiological traits differed between sexes, we performed t-tests for a sub-set of 25 species (Supplementary Table 1). We did not find significant behavioral and thermal physiological differences between (non-gravid/non-pregnant) females and males in T_b (t = 0.172, df = 48, P = 0.86), T_pref (t = −0.482, df = 48, P = 0.63), CT_min (t = 0.742, df = 45, P = 0.46), or CT_max (t = −0.407, df = 42, P = 0.69), so we combined data for both sexes. Ideally, we would rerun all analyses using thermal trait data from gravid/pregnant females, but such data are still lacking. Given that, in the few cases where robust data do exist, preferred temperatures in pregnant females tend to be even lower than in non-pregnant females^46,47,77, we suspect that our analyses provide a relatively conservative estimate of physiological differences among parity modes.

Operative temperatures

As we describe below, we were interested in estimating thermoregulatory patterns among phrynosomatid species. Doing so requires knowledge of the environmental operative temperatures (T_e) available to lizards. T_e represents the equilibrium temperature of an animal in the absence of behavioral thermoregulation⁷⁸. We recorded T_e using previously-calibrated pipe models (made of polyvinylchloride), which were similar in shape, size, and heat gain/loss with respect to lizards of each species (for examples of calibration, see refs. ^56,78,79). Into each pipe model, we inserted one temperature data logger (Thermochron iButton; model DS1921G), which recorded temperature (±0.1 °C) every ten minutes during the same periods (and days) during which we were measuring field-active body temperatures (T_b) in lizards. These models were placed randomly in microsites occupied by lizards in their activity period^56,79. Operative temperatures were typically measured during a sampling period of 1–5 days for each locality, which always occurred during times of the year when lizards exhibit surface activity. The pipe models also recorded temperature during the night (T_e-night), which we used as an approximation of core temperatures of individuals during their inactivity period to model mass-corrected metabolic rates during inactivity (described below). As a caveat, the T_e-night measured in microsites where individuals were active are likely to be somewhat cooler than microenvironmental temperatures experienced by inactive lizards in their retreats.

Thermoregulatory effectiveness

Several studies have found that viviparous species have lower field body temperatures than their oviparous counterparts^11,80. Less well known, however, is whether lower field body temperatures reflect a behaviorally passive property of viviparous lizards, perhaps because of their distributions in relatively cooler habitats, or whether those low field body temperatures reflect a more behaviorally active decision to thermoregulate. Therefore, we were particularly interested in the thermoregulatory patterns of oviparous and viviparous species. We calculated the effectiveness of temperature regulation (E), a ratio that describes how well lizards maintain their T_b within their T_set range (central 50% of data of T_pref; T_set25, and T_set75), given the operative temperatures (T_e) available in their habitat⁷⁵. We estimated E for each species following the equation proposed by Hertz et al.⁷⁵:

where (overline{{db}}) is the average of the accuracy of body temperature, and indicates the deviation of T_b from T_set range. If each T_b < T_set25, then each db = T_set25 − T_b, if each T_b > T_set75, then each db = T_b − T_set75, and if each T_b is within T_set range, then each db = 0. (overline{{de}}) is the average of thermal quality of the habitat, and indicates the deviation of T_e from T_set range. If each T_e < T_set25, then each de = T_set25 − T_e, if each T_e > T_set75, then each de = T_e − T_set75, and if each T_e is within T_set range, then each de = 0. Values of (overline{{db}}) close to 0 indicate that lizards accurately maintain their body temperature within their preferred range, and values of (overline{{de}}) close to 0 indicate that the habitat temperatures approximate (and/or fall within) the preferred range of lizards. As both (overline{{db}}) and (overline{{de}}) increase, body temperatures and operative temperatures, respectively, exceed species’ preferred thermal ranges. As such, E values close to 1 indicate that lizards are highly effective thermoregulators, and E values close to 0 indicate that individuals are more behaviorally passive with respect to the thermal environment. E was only estimated in cases where T_e and T_b were sampled during the same period, and if T_set was measured from the same population of lizards from which T_b was measured. In total, we were able to gather estimates of E from 64 species (37 oviparous and 27 viviparous) of phrynosomatid lizards (Supplementary Data 1).

Environmental temperature

In addition to the operative temperatures, which provided a detailed (but temporally limited) snapshot of the thermal environment, we gathered data on general air temperature trends for each species’ habitat. Specifically, we also gathered climatic measurements for each locality (Supplementary Data 1) from which any lizard trait data were gathered by extracting thermal variables from the environmental layers available in the WorldClim dataset (resolved to ~1 km²)⁸¹. These variables were mean annual temperature (bio1), mean temperature of the warmest quarter (bio10), and mean temperature of the coldest quarter (bio11). We did not use these data to calculate T_e for estimates of thermoregulatory effectiveness (as E should be calculated from T_e measured during the same time period as T_b). Instead, we used these bioclimatic variables as predictors of phenotypic trait variation using evolutionary regressions as described below.

Morphology and life-history traits

We gathered published and unpublished information for mean snout-vent length (SVL; mm), a common measure of body size in squamates, and body mass (g) of adult females and neonates. We also recorded clutch or litter size (i.e., the number of offspring produced per reproductive bout), and the number of clutches or litters produced during 1 year (Supplementary Data 1). We multiplied these two last traits to quantify annual fecundity, which reflects the total predicted annual reproductive output of a given species. We used annual fecundity for three reasons. First, in phrynosomatids (with the exception of some populations of three species^82,83,84), females have annual (seasonal) patterns of reproduction⁵¹. Second, oviparous species tend to produce eggs in multiple clutches per year⁸⁵, whereas viviparous species are typically able to produce only one litter in the same unit time⁵¹. Indeed, viviparous species tend to produce only one litter per year regardless of reproductive window length. For example, both Phrynosoma hernandesi, a species with shorter gestation (3 months)⁸⁶ and Sceloporus bicanthalis, a species with continuous reproduction^83,87, produce a single litter per year. Third, the maximum lifespan for phrynosomatid lizards varies considerably but does not differ between parity modes¹². For some species, the maximum lifespan in natural conditions is ~1 year (documented for the oviparous species, Sceloporus aeneus, and the viviparous species, Sceloporus bicanthalis^83,88), whereas for other species the maximal lifespan can approach ~10 years (documented for the oviparous species Phrynosoma asio⁸⁹ and for the viviparous species Sceloporus macdougalli⁹⁰). Thus, consistent with other studies^15,70, we consider that by standardizing production to one year, we have an estimate of reproductive output that can be readily compared among parity modes.

Metabolic rate

We modeled individual metabolic rate (I) for female lizards following the equation proposed by ref. ¹⁸:

where i₀ = is a normalization constant, M is the mean body mass (g) of females, e = Euler’s number, E = activation energy, k = Boltzmann’s constant, and T = mean-field body temperature (in Kelvin)¹⁸. Previous work has shown that the slope and intercept of the body size ~ metabolic rate relationship vary among vertebrate lineages, but vary much less within lineages (with the conspicuous exception of some amphibian lineages like salamanders)⁹¹. To tailor these equations to lizards, we used i₀ = ln(20.3) (normalization constant for reptiles) and an E value = 0.63 (activation energy excluding endotherms in hibernation and torpor)¹⁸. Then, mass-specific metabolic rate (B) can be modeled simply as I/M, or following the equation proposed by Gillooly et al.^18,92:

The universal validity of this equation has been (rightfully) debated^{93,94,95,96,97,98}, but is nonetheless a useful (and very widely applied) approximation of instantaneous energetic demand across a wide variety of physiological studies^{17,99,100,101,102}. Studies that provide empirical estimates of mass-specific metabolic rate typically do so by keeping the experimental temperature constant. This approach allows researchers to compare the mass-specific metabolic rate among different individuals, populations, and species across shared temperature regimes. Yet, this approach would tell us little about the energetic demands of organisms as expressed in their environments. Here, our goal was to understand how observed body sizes and field activity body temperatures impact the energetic demands of lizards based on their reproductive mode. This approach allowed us to predict how mass-specific metabolic rate should vary based on the observed morphology and thermal physiology of the species (during both activity and inactivity periods; see below), as opposed to comparing the metabolic demands across shared thermal conditions.

Mean body length (SVL) of phrynosomatids lizards is more frequently reported than mean body mass. We built a database of mean body mass and mean SVL of adult females for 30 phrynosomatid species (none were gravid or pregnant) via a combination of unpublished and published information (Supplementary Data 1). Using these data, we built a non-phylogenetic equation to predict log₁₀mean body mass from log₁₀mean SVL. Our equation is log₁₀mean body mass = 3.355log₁₀mean SVL −5.065 (R² = 0.88, P < 0.001). Then, we transformed the log₁₀mean body mass value into an integer value (mean body mass = 10^{log10mean body mass}). With our equation, we predicted the mean body mass of females for species for which mean SVL and mean-field body temperature were available. Based on this approach, we obtained a total database of the mass-specific metabolic rates of females for 95 phrynosomatid species (Supplementary Data 1).

Certainly, using mean body mass and mean-field body temperature of each species could under- or overestimate mass-specific metabolic rate (“the fallacy of averages”)^101,103. Likewise, it is also relevant to know the instantaneous energetic metabolic demands of species given the body mass or the field body temperature of individuals. Therefore, we also modeled temperature-corrected metabolic rate for each body mass measured in each adult female from 38 phrynosomatid species, and we modeled mass-corrected metabolic rate for each T_b measured in each individual of 65 species, following the equations proposed by ref. ¹⁸: (4) Temperature-corrected metabolic rate = 0.71*ln (mass) + 18.02, and (5) mass-corrected metabolic rate = −0.69*temperature (1/kT) + 20.3. Because T_b can change throughout the day, this approach includes a diel variation on energetic demand within our analyses. Phrynosomatids are diurnal lizards, and metabolic rates from field-estimated body temperatures reflect energetic demand during lizards’ period of activity. To test hypotheses about the energetic demands related to parity mode, it is also relevant to know metabolic rates during times of inactivity (like at night). We did not measure the field body temperature of inactive lizards. However, diurnal lizards exhibit a limited ability to thermoregulate at night (a time when thermal environments also tend to homogenize)¹⁰⁴, and their body temperatures, correspondingly, tend to correlate with equilibrium (operative) temperatures (i.e., T_e-night)^55,56.

So, we also modeled inactivity mass-corrected metabolic rate for each T_e-night recorded for each null model during the inactivity time of each species (typically from 18:00 to 08:00 h) in the same localities where field-active body temperatures were measured for 42 species, following the same equation mentioned above. As noted above, the surface perches from which we gathered operative temperatures are likely a bit cooler than nighttime retreats utilized by phrynosomatids and, correspondingly, should slightly underestimate energetic demand during inactivity. Lastly, we estimated the mean temperature-corrected MR, and mean mass-corrected MR during activity and inactivity for each species (Supplementary Data 1). Thus, our approach considers how mass-corrected metabolic rate varies among parity modes during both active and inactive periods.

Mass-specific production

We estimated mass-specific production (Pr) as the product of neonate mass and annual fecundity/female body mass¹⁶. Therefore, Pr describes the amount of energy converted into reproductive effort per year, normalized by maternal body mass.

Evolutionary analyses

All evolutionary analyses were conducted using the R environment ver. 3.6.0¹⁰⁵, with the exception of the phylogenetic logistic regressions, which were performed using ver. 4.1.1.

Stochastic character mapping of parity mode

To estimate the number of transitions between parity modes we performed stochastic character mapping¹⁰⁶ onto the ultrametric tree of Phrynosomatidae using the make.simmap function with 500 simulations and a transition model of equal rates (ER) in phytools (ver. 0.6.99) R package¹⁰⁷. We selected the ER model of character evolution because was it the least-complex, best-supported model (∆AICC = 1.3, weight = 0.26) in comparison to a symmetrical model (SYM; ∆AICC = 1.3, weight = 0.26) and with an all-rates-different model (ARD; ∆AICC = 0, weight = 0.48).

Ancestral state reconstruction

To fit the mean annual temperature through the Phrynosomatidae tree and graphically show the thermal environment where each population of each species used in this study inhabits, we performed ancestral state reconstruction using contMap function in phytools (ver. 0.6.99) R package¹⁰⁷.

Phylogenetic analyses of variance (ANOVA)

To test for differences in the effectiveness of temperature regulation, among parity modes we performed phylogenetic ANOVAs using the phylANOVA function with 500 simulations in phytools (ver. 0.6.99) R package¹⁰⁷.

Comparing trait evolution between viviparous and oviparous species

We were interested in whether transitions to viviparity are associated with predictable phenotypic shifts. To this end, we tested if parity mode (“oviparous” or “viviparous”) was associated with different evolutionary patterns of mass-specific metabolic rate, mass-corrected metabolic rate (both during periods of activity and inactivity), temperature-corrected metabolic rate, mass-specific production, body mass and size, thermal physiological traits, and life-history traits by fitting Brownian motion (BM) and Ornstein-Uhlenbeck (OU) models. To do so, we used the R package OUwie (ver. 1.57)¹⁰⁸ and the 500 stochastic character maps of parity mode built with the make.simmap function in the R package phytools (ver. 0.6.99)¹⁰⁷. We fitted three different models. The simplest (BM1) is a single-rate BM model in which a single rate of stochastic trait evolution (σ²) was estimated for all Phrynosomatidae, and phenotypic differences among species are proportionate to branch length (or time). The other two models were all adaptive OU models that varied in whether the estimated phenotypic optimum (θ) was either constrained to be equal among parity modes (OU1) or allowed to vary between oviparous and viviparous species (OUM). We fitted these three models separately for each physiological trait (T_b, T_pref, CT_min, CT_max, mass-specific metabolic rate, mass-corrected metabolic rate, and temperature-corrected metabolic rate), each morphological variable (adult body mass, and adult body size), and each life-history trait (offspring mass, offspring size, annual fecundity, and mass-specific production) (Supplementary Table 2). For these (and all) analyses, body mass, body size, offspring mass, offspring size, and annual fecundity were log₁₀ transformed. We assessed model fit using a modified Akaike information criterion (AIC_C) that incorporates a correction for a small sample size¹⁰⁹. Our approach, which was based on 500 stochastic character maps derived from the MCC tree, allowed us to account for uncertainty in reconstruction across the preferred tree, but could not account for uncertainty in the phylogeny itself. Therefore, we also repeated our stochastic character mapping across 500 individually sampled trees from the posterior distribution to account for this additional source of phylogenetic uncertainty and repeated all of our OUwie analyses using these 500 sampled trees. Our results in this latter approach are comparable to those using the MCC tree (Supplementary Table 3). Therefore, we present our results from the analyses based on the MCC tree in the main document.

More complex models, like OU models, can be incorrectly favored over simpler models if the statistical power of the analysis is weak^110,111. To assess the adequacy of our phylogenetic data for model-fitting procedures, we performed simulations to assess the probability of type-I error in model fit. We simulated trait evolution in two ways to be consistent with our OU analyses. First, we first simulated trait evolution 500 times on the MCC tree, and then simulated trait data once for each of the 500 trees from the posterior distribution. We then fitted BM and OU models to the simulated data to determine what percentage of analyses would incorrectly favor OU over BM. False positives were not an issue (Supplementary Table 2 and 3), supporting that we could reasonably fit OU models to our data.

Testing the strength of convergent evolution

Given that viviparity repeatedly evolved in phrynosomatids (Fig. 1a), we were interested in estimating the strength of convergence in morphology, thermal and metabolic physiology, and life-history traits of viviparous species. To this end, we estimated the Wheatsheaf index (w), which quantifies the strength of convergence (or lack thereof) for continuous traits in focal species¹¹². We assigned viviparous species as focal groups (recognized a priori with the ancestral state reconstruction), and we estimated w using the test.windex function and 500 bootstrap replicates (to estimate a P-value) in the R package windex (ver. 2.0.2)¹¹³. A high w value indicates stronger convergence, and a P-value <0.05 represents that convergence is significantly stronger after accounting for relatedness among species^112,113.

Phylogenetic Generalized Least Squares (PGLS)

We performed PGLS regressions using the gls function in the R package nlme (ver. 3.1.139)¹¹⁴ to know the evolutionary relationships between adult body size and adult body mass, between reproductive response traits (clutch size (or litter size) and offspring size), and adult body size, and between thermal physiological response variables (CT_min, T_b, T_pref, and CT_max) and environmental predictors (bio1, bio10, bio11, and T_e). For the PGLS analysis (and for the SLOUCH analysis, described below) between T_pref and environmental predictors we did not include Sceloporus graciosus, which represents one outlier point (as this species inhabits an extremely cold habitat). If we include S. graciosus, the PGLS regression is significant for oviparous species, but the strength of the relationship is weak (y = 0.094x + 32.78, P = 0.03).

Assessing environmental predictors of parity mode evolution

Given the strong conceptual framework linking the evolution of viviparity to cold environments, we tested whether changes in the thermal environment were strong predictors of parity mode shifts using two approaches. We tested for the evolutionary covariation between the thermal environment and reproductive mode (oviparous vs. viviparous) using the threshold model^115,116 using threshBayes function in the phytools (ver. 0.6.99) R package¹⁰⁷ and with phylogenetic logistic regression using phyloglm function in the phylolm (ver. 2.6.2) R package¹¹⁷. The threshold model is used to test for evolutionary covariation between continuous and discrete traits¹¹⁶. Under the threshold model, a discrete character (i.e., oviparity or viviparity) evolves as a function of a continuously varying feature (termed “liability”). When the value of “liability” crosses a certain threshold, the state of the discrete character evolves (i.e., a transition from oviparity to viviparity occurs)^115,116. We ran threshBayes for 1.0 × 10⁶ generations, sampling every 100 generations, and discarding the first 200 K generations as burn-in. We ran separate analyses for mean annual temperature (bio1), mean temperature of the coldest quarter (bio10), and mean temperature of the warmest quarter (bio11). The phylogenetic logistic regression is used to test if dependent variables (binary traits that switch between 0 and 1) are predicted by independent variables (continuous or discrete)⁵⁴. So, we coded oviparity and viviparity using 0 and 1, respectively. As above, our predictor variables in these analyses were bio1, bio10, and bio11.

Stochastic linear Ornstein-Uhlenbeck models

Our OUwie analyses revealed reductions in the phenotypic optimum (θ parameter) for thermal traits in viviparous lizards (see Results and discussion). Yet, it is unclear whether reductions in thermal physiology reflect adaptation to cool environments (given, for example, the greater representation of viviparous lineages at high elevation^20,21) or, instead, reflect energetic adjustments for life-history evolution (hypothesis i in Table 1a), which could be readily co-opted for life in cold environments. If cool-adapted physiology reflects adaptation to cool environments, there should be a strong evolutionary association between the local thermal environment and thermal physiology. However, if cool-adapted physiology more reflects a number of potential trade-offs like life-history energetics, then we expect viviparous species to exhibit more cool-adapted physiology than oviparous species regardless of ambient conditions, which should weaken the evolutionary relationship between the local thermal environment and thermal physiology.

To test these ideas, we used the SLOUCH model of ref. ¹¹⁸, which simultaneously estimates an “evolutionary regression” and an “optimal regression” in an OU framework. The evolutionary regression describes the observed relationship between climatic predictors (mean annual temperature (bio1), mean temperature of the warmest quarter (bio10), and mean temperature of the coldest quarter (bio11)) and physiological response variables (CT_min, T_b, T_pref, and CT_max), while accounting for the relatedness among species. The estimated “optimal regression”, by contrast, describes the relationship between these variables predicted under an OU model, and assumes adaptation of the response variables to the predictor variables. In addition to the regressions, the model permits the estimation of phylogenetic half-life (t_1/2) and rate of adaptation (α). Phylogenetic half-life represents the amount of time required for viviparous or oviparous lineages to get halfway to their thermal physiological optimum. So, a short t_1/2 (relative to the length of tree) indicates the phylogenetic signal degrades at a rapid pace. By contrast, a t_1/2 approaching (or exceeding) the length of the tree, indicates a strong phylogenetic signal. A rate of adaptation close to 0 represents a very slow physiological adaptation to thermal predictors, whereas α values higher than 100 (or approaching ∞ ) indicate instantaneous physiological adaptation to thermal predictors¹¹⁹.

The similarity between the evolutionary and optimal regressions is supported when t_1/2 is close to 0, which would indicate that transitions in the thermal environment are rapidly coupled with changes in thermal physiology. Differences in the slopes of these relationships, by contrast, are supported when the phylogenetic half-life (t_1/2) of the model is bounded away from zero, implying phylogenetic inertia, or a lag in physiological adaptation to the thermal environment. Such lags are consistent with the Bogert effect, in which behavioral preferences disrupt physiological adaptation to prevailing environmental conditions^120,121,122. Under this scenario, shifts in the thermal environment are predicted to be weakly associated with shifts in thermal physiology. To run the analyses, we simultaneously estimated the evolutionary regression, optimal regression, and t_1/2 for each thermal physiological trait (CT_min, T_b, T_pref, and CT_max) of phrynosomatid lizards, with respect to their thermal environment (bio1, bio10, and bio11) in an OU modeling framework using the slouch.fit function the R package SLOUCH (ver. 2.1.2)¹¹⁸.

Graphics

Figures 1b, 2, and 3 were generated using the R package ggplot2 (ver. 3.2.1)¹²³, and edited using Adobe Illustrator.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com