Alternative reproductive adaptations predict asymmetric responses to climate change in lizards

Anthropogenic climate change ranks among the major global-scale threats to modern biodiversity. Extinction risks are known to increase via the interactions between rapid climatic alterations and environmentally-sensitive species traits that fail to adapt to those changes. Accumulating evidence reveals the influence of ecophysiological, ecological and phenological factors as drivers underlying demographic collapses that lead to population extinctions. However, the extent to which life-history traits influence population responses to climate change remains largely unexplored. The emerging ‘cul-de-sac hypothesis’ predicts that reptilian viviparity (‘live-bearing’ reproduction), a ‘key innovation’ facilitating historical invasions of cold climates, increases extinction risks under progressively warming climates compared to oviparous reproduction – as warming advances polewards/mountainwards, historically cold-climates shrink, leading viviparous species to face demographic collapses. We present the first large-scale test of this prediction based on multiple lizard radiations and on future projections of climate-based ecological niche models. Viviparous species were found to experience stronger elevational range shifts (and potentially increased extinctions) in coming decades, compared to oviparous lizards. Therefore, our analyses support the hypothesis’s fundamental prediction that elevational shifts are more severe in viviparous species, and highlight the role that life-history adaptations play in the responses of biodiversity to ongoing climate change.

tolerance as a function of intrinsic factors such as genetic variation, population size, and generation time [23][24][25] . In recent years, research efforts have focused on identifying mechanisms behind such spatial asymmetries in responses to contemporary climate change, leading to reinforce the focus on a range of phenological, ecological, and physiological factors as primary factors underlying extinction risks [26][27][28] .
Surprisingly, the influence that life-history adaptations exert on species responses to climate change has only occasionally been addressed. In fact, although parity-mode transitions (from live-bearing to egg-laying) have been suggested to play a key role in accelerating extinction risks under rapid climate warming 19,26 , only a handful of studies have investigated the ecological basis of this phenomenon 19,26,29 . Given that life-history traits underlie reproductive success ("fitness"), and that climate change operates as rapid natural selection on fitness via its impact on environmentally-sensitive traits, the conceptual disjoint between both dimensions creates a major gap in our understanding of biodiversity declines via anthropogenic effects. In fact, selection arising from climatic pressures have been shown to generate predictable patterns of evolution and distribution of reproductive traits (e.g., parity mode 19 , transient fecundity 30 ), which therefore creates macroecological patterns of species interactions with the climate via life-history adaptations 19,[30][31][32][33] .
An emerging idea, the 'cul-de-sac' hypothesis, predicts that the evolution of viviparity -a key life-history adaptation believed to have facilitated the radiation of reptiles into cold climates 19,33,34 -is driving reptiles to extinctions due to climate change, relative to oviparous reproducers 19 . Thanks in part to evolving viviparity, reptiles have successfully proliferated across extreme high elevations and latitudes, where low temperatures create strong natural selection against multiple components of fitness via thermal demands on ecological (e.g., embryo survival and development 19,26 and life-history (e.g., limitations on reproductive output for viviparous species) functions 26,32,35 . In contrast, these pressures are relaxed in warm climates, where oviparous reptiles dominate. In viviparous squamates, embryo retention in the maternal body acts similar to an incubator to counteract the lethal effects that low and fluctuating environmental temperatures impose on egg development 19,[36][37][38] . However, cold-climate reptiles may be approaching an evolutionary dead-end in the face of climate change 19 . Given that most cold-climate reptiles are viviparous 30,32,39 , that the oviparity-to-viviparity transitions are largely irreversible (thus, fundamentally unidirectional) 19,32,39,40 , and that live-bearing parity (relative to oviparity) entails high fitness costs (e.g., reduced reproductive frequency, increased pregnancy-burden, prolonged basking time) 32,35,41 , the cul-de-sac hypothesis suggests that, as climate warming advances towards higher elevations and latitudes, viviparous reptiles will remain trapped in rapidly retracting cold environments 19 , while oviparous species are expected to colonize previously cold sites 19,26 . As viviparous species progressively run out of suitable environments as they approach mountaintops and continental edges, populations are expected to experience demographic collapse that will lead to extinctions 10,14,19,42 . So far, evidence consistent with this hypothesis remains limited to a few studies that have approached the question in a preliminary fashion 19,26,43 . However, a replicated, large-scale macroecological test assessing the responses of species to climate change as a function of life-history adaptations proposed by this hypothesis remains lacking.
Here, we present a comprehensive test of the cul-de-sac hypothesis, based on fine-scale distributional data spanning three prolific and evolutionarily contrasting lizard radiations that have diversified in the Andes and adjacent areas in South America. Their natural histories span all scenarios of parity mode variation: Liolaemus, with 270+ species is the world's second most diverse genus of living amniotes 44 , in which viviparity has evolved in several independent phylogenetic events 19,45 ; Phymaturus, a strictly cold-climate, viviparous genus of 60+ species 46 ; finally, Stenocercus consists of 60+ oviparous species, despite their distributions up to above 4000 m in the Andes 47 . Using an ecological niche modeling framework and future climate scenarios (IPPC Fifth Assessment 48 ) we tested whether viviparous (relative to oviparous) species (i) are predicted to experience greater range contractions, (ii) higher overall spatial displacements, and (iii) greater displacement towards higher elevations in relation with oviparous species.

Results
Our ecological niche model (ENM) analyses reveal that the geographic ranges and their changes are predominantly influenced by temperature-related variables across all three lineages (in 93.7% of the species), while precipitation-related variables tend to be considerably less relevant (only dominant in 6.3% of the species; Supplementary Table S1). Also, we found that potential impacts of climate change on species ranges were not evenly distributed across space. Range expansions and contractions were more likely in certain regions (e.g., the high Andes), and absent from other regions (e.g., Patagonian steppe, Temperate broadleaf and mixed forests in Southern Chile; Moran's spatial autocorrelation index >0.7, P < 0.001; Supplementary Table S8). Such range shifts patterns were observed consistently in viviparous species (Figs 1a,b and S1a,b). Although MIROC5-based models supported the hypothesis of a geographic directionality in the spatial shifts (Rayleigh's test range, P = 0.05-0.01), GISS-ER-based models showed centroid shifts that were not distinguishable from uniform ( Fig. 1c Geographic range shifts. Our analyses revealed that shifts in the percentage of range expansions or contractions driven by estimated magnitudes of climate change do not significantly differ between viviparous species combined (regardless of taxonomic group) and oviparous species combined (Table 1; Fig. 1a,b, Supplementary Tables S2-S5, Supplementary Fig. S1a,b). Likewise, the effects of climate change on the magnitude of spatial displacements of species' geographic range centroids (Fig. 2a,b) and their spatial direction (Fig. 1c, Supplementary Table S6, Supplementary Fig. S1c) do not significantly differ between viviparous species combined and oviparous species combined, in concordance with observations on expansions and contractions.
In contrast, as predicted by theory, we found significant differences in the magnitude of elevational shifts between parity modes (Liolaemus-viviparous + Phymaturus versus Liolaemus-oviparous + Stenocercus; Fig. 2c,d), expressed as higher elevational shifts for viviparous species relative to oviparous species (Phylogenetic t-test www.nature.com/scientificreports www.nature.com/scientificreports/ range, t 141 = 20.1-36.9, P = 0.02-0.004, Table 1; Fig. 2c). However, these differences in elevational shifts were not significant when compared among 'parity-by-taxonomy' groups (across Stenocercus, Liolaemus-oviparous, Liolaemus-viviparous and Phymaturus; Phylogenetic ANOVA range, F 3,139 = 8.3-12.4, P = 0.27-0.41, Table 2, Fig. 2d). More specifically, under RCP 2.6 scenarios, Liolaemus-viviparous and Phymaturus lizards were found to have a projected increase in their elevational ranges of around 12% (relative to the present), and up to 21% for RCP 8.5 scenarios. In contrast, oviparous groups (Stenocercus and Liolaemus-oviparous) showed a decrease of ~2% in both RCPs scenarios ( Fig. 1d     www.nature.com/scientificreports www.nature.com/scientificreports/ Future range overlaps between oviparous and viviparous species. Compared with the current range overlap between the geographic areas covered by oviparous and viviparous species combined regardless of taxonomy (currently 15.3%), our analyses suggest that climate change-driven range shifts are likely to cause an increase in the current levels of overlap between species with both parity modes. More specifically, range overlaps are projected to increase with the magnitude of greenhouse gas emissions (RCPs): while under the minimum climate change scenario (RCP 2.6), range overlap was observed to remain as today (~15% of overlap, relative to the modern 15.3%), under the maximum scenario (RCP 8.5) the predicted overlap is projected to increase up to 3%, leading to ~18% of overlap between both parity modes (Supplementary Table S9 and Fig. S2).

Discussion
Our study presents the first comprehensive and broad-scale quantitative test of the 'cul-de-sac' hypothesis 19 , which predicts that viviparous species are more vulnerable to extinctions as climate warming progressively shrinks the cold regions they are adapted to 19 . Our findings across three highly-diverse lizard radiations support the effects that parity modes are expected to exert on the patterns of elevational shifts projected for these organisms as climate warms up (whereas, no such differences were observed for direction of shifts in range centroids). Specifically, our findings revealed that the prediction strongly holds for alterations in elevational distributions, where viviparous species combined experience greater elevational shifts relative to oviparous species (see Fig. 1d). These elevational displacements observed in viviparous species were stronger in Phymaturus (strictly viviparous), suggesting the potential for greater magnitudes of demographic impacts, and potentially extinction risks faced by these lizards under the ongoing climate trends 49 . Importantly, our findings also suggest that as the degrees of emissions of greenhouse gases increase, the extent of spatial overlap between oviparous and viviparous reproducers will increase, thus leading to the predicted scenarios of novel forms of competition between oviparous and viviparous species 19,26 . Collectively, therefore, our study reinforces the need to add the life-history dimension to the search for factors that trigger mechanisms underlying biodiversity alterations caused by increased anthropogenic climate change 19 . It also adds a further component to the sustained efforts to establish the combination of factors that define the 'profile' (i.e., shared combinations of traits) of species that have entered a current phase of extinction risk (or, on the other hand, of species which have seen their populations remain stable or even increase) as a result of environmental alterations.
At a more general level, the pathway towards species declines is a function of alterations to the ecological and demographic stability of populations, which can occur due to different mechanisms. For example, as a consequence of loss of their suitable distributional area, which can erode the genetic diversity to respond to rapidly changing selection regimes 42,[50][51][52] . However, other species may not be experiencing shrinking of their geographic ranges, and yet, may still face demographic alterations via other mechanisms such as novel forms of competition with species in newly encountered assemblages 53 , or exposure to climatically suitable, but not necessarily structurally suitable (e.g., topographically) environments. Therefore, range shifts in diverse directions rather than towards higher latitudes (Fig. 1) are expected to impact on the ecological and demographic stability of those species by exposing them to novel environments and networks of species where selection regimes will be different. In some cases, these impacts may lead to population declines and thus, to increased risk of extinction via factors not necessarily connected directly with changing temperatures. For instance, a recent study on Liolaemus lizards showed that, in contrast with the traditionally established 'cold-climate hypothesis' that declines in environmental temperatures operate as drivers of natural selection for viviparity 32,54 , reductions in atmospheric oxygen levels play a central role as agents of selection, promoting embryo retention and evolution of viviparity 34 . Whether oxygen levels influence the magnitude and direction of alterations in range shifts of viviparous species via their interactions with thermal gradients remains an open question. However, this 'hypoxia' hypothesis contributes to the more general view that species' range stability can be broken down into interactions between a range of external factors (e.g., climate, topography) and population features (e.g., abundance, morphology, life-history, physiology). Future studies could develop predictions about the effects of environmental factors on parity modes taking into account this strong role of oxygen, in combination with classic ecophysiological factors. Such predictions could expand the findings of, for example, Medina et al. Medina,et al. 29 , who observed that the preferred body temperature in the laboratory (T pref ) does not differ among Liolaemus species with different parity models (while, however, field body temperature in the laboratory (T pref ) does not differ among Liolaemus species with different parity models (while, however, field body temperature, T b , does).
This study adds a further conceptual and empirical layer to the urgent need to develop 'profiles' of species facing risks of extinction. These profiles aim to summarize our current understanding of factors that contribute to altering the historical stability of populations, and thus, to reinforce the efficiency of actions developed to mitigate ongoing declines of biodiversity. However, one of the main limitations of this study is the bias generated by not considering the extremely endemic species in our analysis, which is related to the lack of enough geographic information to predict the potential distribution of these species.
Overall, this study shows that life-history adaptations, favoured by natural selection, can turn into determinant factors pushing species to decline under contemporary climate change. Life-history traits, given their direct effects on fitness, are therefore primary candidate traits to disentangle factors that contribute to extinction, and thus must be incorporated into empirical and theoretical studies aiming to develop estimations of climate change effects on life on earth.

Species geographic distribution.
We gathered a large-scale dataset encompassing three lizard radiations differing in their patterns of geographic distribution and in parity modes, which spans large species-samples of Liolaemus (40 oviparous and 52 viviparous species), Phymaturus (11 species, exclusively viviparous -although our dataset contained data for most species within this genus, we excluded most of them given their extremely www.nature.com/scientificreports www.nature.com/scientificreports/ limited geographic range sizes, which makes it inviable to perform the analyses employed to address our core predictions), and Stenocercus (40 species, exclusively oviparous) (Supplementary Fig. S1). The distributional and life-history data comprise 4,532 geographic occurrence records (after all duplicated points collected from different individuals at same localities were removed). This includes all known records of presence of these lizards following 20+ years of field and museum work, and museum-validated published occurrences 45,47,[55][56][57][58][59][60][61] . All the occurrences were individually checked and confirmed by experts (i.e., Daniel Pincheira-Donoso: Liolaemus + Phymaturus and Omar Torres-Carvajal: Stenocercus, to assure accuracy).

Environmental predictors.
To analyse the environmental space occupied by lizard lineages, we used the bioclimatic variables characterizing climate during the 1970-2000 period, obtained from the WorldClim 2 data repository 62 (available at: http://www.worldclim.org/version2) at a spatial resolution of 30 seconds (~1 km). To reduce collinearity between the environmental variables, we used VIF (Variance Inflation Factors) implemented in the "usdm" R-package 63 . Using this approach, we excluded all the highly correlated variables from the model (VIF greater than 10), which is associated with a signal that the model has a collinearity problem 64 . This method is based on the square of the multiple correlation coefficient (R 2 ) resulting from regressing the predictor variable against all other predictor variables. We used the remaining uncorrelated variables to calibrate models: mean annual temperature (bio1), mean diurnal range (bio2), isothermality (bio3), temperature seasonality (bio4), max temperature of warmest month (bio5), min temperature of coldest month (bio6), temperature annual range (bio7), mean temperature of warmest quarter (bio10), mean temperature of coldest quarter (bio11), annual precipitation (bio12) precipitation of wettest month (bio13), precipitation of driest month (bio14), and precipitation seasonality (bio15), precipitation of wettest quarter (bio16), and precipitation of driest quarter (bio17).
To explore the occupation of environmental space by Liolaemus, Phymaturus, and Stenocercus and their reproductive modes, we performed a principal component analysis (PCA) of present-day climate conditions in South America using NicheA v3.0 software 65 (available at: http://nichea.sourceforge.net), an open-source application that analyses ecological niches (or 'climatic spaces') in both environmental and geographic space ( Supplementary  Fig. S3a,b).
To estimate impacts of climate change on distributions of species, we employed two models of future climate conditions, the GISS-ER and the MIROC5 general circulation models for the period around 2070 (2061-2080). Each model consists of two Representative Concentration Pathways (RCPs) versions representing a stringent mitigation scenario (RCP2.6, which predicts magnitude of climatic variation for the parameters of the model) and one scenario assuming high anthropogenic greenhouse gas (GHG) emissions (RCP8.5, which predicts the highest magnitude of climatic change for the parameters of the model) as alternative scenarios of climate change 66 (IPCC, Fifth Assessment), because it is extremely likely that human activities caused more than half of the observed increase in GMST (Global Mean Surface Temperature) from 1951 to 2010 67 . We selected these models based on their high resolution and incorporation of covariates 68,69 . Species distribution modeling. Based on the occurrence data ( Supplementary Fig. 1), we calibrated Ecological Niche Models (ENMs) to estimate current and future potential distribution for each lizard species. To mitigate oversampling effects in our model, occurrences were re-sampled to one point per pixel with respect to the environmental grids. Coordinates were divided into two groups for calibration and other evaluation, based on four quadrants with similar numbers of points, using two off-diagonal quadrants for calibration and two for evaluation. ENMs were developed using Maxent 3.3.3k 70 , a presence-background software that estimates environmental suitability via an index of similarity that resembles a heterogeneous occurrence process or logistic regression function 70,71 . We used Maxent with clamping and extrapolation turned off (i.e., no prediction outside the range of environmental conditions in the calibration areas). To facilitate interpretations, we used the relative probability of presence as a proxy of environmental suitability 71 . Suitable area for each species was estimated as a Boolean (presence/absence) map that was thresholded based on the minimum training presence 72 . To determine the model parametrization with the best fit to the data available, we assessed six models for each species under different regularization multiplier values (i.e., 0.5, 1, 5, 10, 15, and 20) 73 . Then, we assessed AIC (Akaike Information Criterion; implemented in ENMTools 74 ) values to choose the model with the best performance (i.e., descriptive-model evaluation; Supplementary Table S10). Once the models with the best fit were selected, we assessed their performance using the software PartialROC 74 , considering an α = 0.05 to determine whether models anticipate independent occurrences better than by random expectations.
To determine changes between the current and future potential distributions, we identified the expansion, contraction, and stability ('no change') of areas predicted climatically suitable by using the SDM toolbox in ArcGIS v.10.2 75 . Direction of the potential distributional changes was calculated based on the centroid of the species' potential ranges. We applied Rayleigh's Test of Uniformity 76 to compare direction of centroid shifts vs a uniform circular distribution reflecting a null hypothesis of random centroid shifts. Rayleigh's statistic quantifies the angular dispersion among the vectors from 0 (representing uniform dispersion) to 1 (indicating complete concentration in a single direction).
Finally, to determine the spatial relationships between oviparous and viviparous species, we quantified the range overlap among the predicted distributions of oviparous and viviparous species. These calculations were made using the predicted present-day and future SDMs under GISS-ER and the MIROC5 models and each greenhouse emission scenarios (RCP 2.6 and 8.5).

Climate change impacts on life-histories.
To determine effects of climate change on parity mode, we performed an independent-sample t-test using the variables: (i) range extents of current and future ENMs, (ii) magnitude and direction of distributions based on range centroids, and (iii) average elevation of current and future ranges. Analyses were corrected for phylogenetic non-independence to test whether the predicted impacts www.nature.com/scientificreports www.nature.com/scientificreports/ of climate change (i.e., as geographic range, centroid or elevational shifts) were statistically significant between species with different parity mode (i.e., oviparous and viviparous), as well as among 'parity-by-taxonomy' groups (i.e., Stenocercus, Liolaemus-oviparous, Liolaemus-viviparous and Phymaturus) using the fully-sampled phylogeny of Squamates recently published by Tonini, et al. 77 which spans all 143 species in our dataset. We performed phylogenetic analyses (t-test and ANOVA) implemented in the R-package "phytools" 78 . Additionally, to calculate the differences between each independent pair of groups we applied ANOVA test with "holm" post-hoc pairwise method, owing to it is one of the most powerful multitests when the sample sizes of the groups are different, to perform this analysis, we used the R-package "DescTools". These analyses were performed for the two models of climate change (GISS-ER and MIROC5) and for each scenario of greenhouse gases emissions indicated above (RCP 2.6 and 8.5). Finally, we performed a spatial autocorrelation analysis (Global Moran's I) using the Spatial Statistics toolbox in ArcGIS to evaluate whether geographic rage shifts were clustered, dispersed, or randomly distributed.