Evolutionary origins of viviparity consistent with palaeoclimate and lineage diversification

It is of fundamental importance for the field of evolutionary biology to understand when and why major evolutionary transitions occur. Live‐bearing young (viviparity) is a major evolutionary change and has evolved from egg‐laying (oviparity) independently in many vertebrate lineages and most abundantly in lizards and snakes. Although contemporary viviparous squamate species generally occupy cold climatic regions across the globe, it is not known whether viviparity evolved as a response to cold climate in the first place. Here, we used available published time‐calibrated squamate phylogenies and parity data on 3,498 taxa. We compared the accumulation of transitions from oviparity to viviparity relative to background diversification and a simulated binary trait. Extracting the date of each transition in the phylogenies and informed by 65 my of global palaeoclimatic data, we tested the nonexclusive hypotheses that viviparity evolved under the following: (a) cold, (b) long‐term stable climatic conditions and (c) with background diversification rate. We show that stable and long‐lasting cold climatic conditions are correlated with transitions to viviparity across squamates. This correlation of parity mode and palaeoclimate is mirrored by background diversification in squamates, and simulations of a binary trait also showed a similar association with palaeoclimate, meaning that trait evolution cannot be separated from squamate lineage diversification. We suggest that parity mode transitions depend on environmental and intrinsic effects and that background diversification rate may be a factor in trait diversification more generally.

and environmental contexts. Therefore, advantages and disadvantages associated with viviparity also differ across animal groups.
Disadvantages are reduced fecundity, maternal death resulting in death of all offspring, and a higher investment for the mother (Pollux et al., 2009;Recknagel & Elmer, 2019;Shine et al., 1978;Wourms & Lombardi, 1992). Squamates, the reptilian order including lizards and snakes, are the lineage in which the most transitions from oviparity to viviparity have occurred (Sites et al., 2011).
Climate plays a major role in modulating those advantages and disadvantages in squamate reptiles (Ma et al., 2018;Shine, 2005;Tinkle & Gibbons, 1977), yet previous studies did not assess the relationship between historical temperatures on the origins of viviparity on a global scale.
Their numerous origins of viviparity makes squamates an ideal model system to understand which environmental factors are associated with the evolution of viviparity. The frequency of viviparous relative to oviparous squamates increases with their distance from the equator. This led to the cold-climate hypothesis, which predicts that egg-layers are disfavoured in colder environments because their eggs take longer to develop in the cold soil or might die from exposure (Tinkle & Gibbons, 1977). Past studies have attempted to infer when viviparity originated in some squamate species groups and they generally supported the cold-climate hypothesis (Lambert & Wiens, 2013;Lynch, 2009). These studies found that the majority of transitions occurred during cold periods in the Cenozoic (the last 65 million years) (Lambert & Wiens, 2013;Lynch, 2009;Schulte & Moreno-Roark, 2010). In contrast, another study found no effect of an epoch-based increase in transition frequency to viviparity in reptiles and suggested that cold climate did not play a substantial role (King & Lee, 2015). Testing the cold-climate hypothesis poses difficulties because present and past conditions differ. A correlation between contemporary conditions in the present distribution and the proportion of viviparous species does not imply these conditions caused the origin of viviparity (Blackburn, 2006). The role of palaeoclimate has not been tested at a broad phylogenetic scale in a statistical framework. This is needed because only by statistically testing a link between high-resolution palaeoclimate data and phylogenetically robust transitions can the promoters of viviparity on geological and global scales be identified.
As a component of climate, the role of environmental stability for the evolution of viviparity is unclear. The maternal-manipulation hypothesis proposes that viviparity is favourable in variable environments as the female can control optimal temperatures for the development of the embryo (Shine, 1995;Webb et al., 2006).
One broad-scale analysis found no support for the maternalmanipulation hypothesis (Watson et al., 2014), and two other studies in snakes  and lizards (Pincheira-Donoso et al., 2017) concluded that viviparity is not necessarily adaptive in variable environments. However, in these contexts environmental stability usually relates to seasonal, short temporal scale variability.
Stearns proposed that the life-history tactic of a few well-developed offspring have a higher chance of survival than numerous but much less developed offspring in a stable environment (K-selection versus r-selection; Stearns, 1976) and such micro-scale ecological tradeoffs may ultimately result in differential phylogenetic patterns (Sibly et al., 2012), though the links from microevolutionary life-history strategies to macroevolutionary patterns remain difficult to draw.
Thus, the impact of environmental variation on the evolution of viviparity has not been tested with data from a geological timescale.
In contrast to extrinsic climate and environmental influences, another factor that may influence the number of trait transitions is the 'background' diversification rate. The contribution of lineage diversification to character transition frequency can be two-fold.
First, there can be purely stochastic effects, with diversification leading to more lineages that can change states. For example, previous research has shown that lineage diversification and phenotypic diversity correlate on macroevolutionary scales (Rabosky et al., 2013(Rabosky et al., , 2014, although this pattern seems to also vary considerably across lineages (Adams et al., 2009;Reaney et al., 2018).
Mechanistically, more genetic variation may result in trait diversification and repeated evolution (Blount et al., 2012;Martin & Richards, 2019), leading to an increase in phenotypic diversity.
Second, diversification itself may be favoured when transition events lead to speciation. Examples for such cases include evolution of an innovation followed by strong reproductive isolation between the lineages with the ancestral state and the innovation.
Subsequently, the lineage with the innovation may undergo further adaptive diversification, again resulting in an increase in lineages. For example, it has been suggested that viviparous species have a higher speciation rate than oviparous species in squamates (Lambert & Wiens, 2013;Pyron & Burbrink, 2014). However, although this would explain higher diversification rates in squamates, it would not explain why transitions to viviparity occur.
Distinguishing these alternatives is challenging as they may be biologically entwined.
Here, we assessed the nonexclusive hypotheses that evolutionary origins of viviparous reproductive mode from oviparous related to the following: (a) cold climate (as low temperature), (b) environmental stability (as long-term stability in temperature) and (c) background diversification rate. We did so by performing an ancestral state reconstruction using 3,498 taxa with parity mode data and published time-calibrated phylogenies. Identified evolutionary transitions from oviparity to viviparity were tested for association with 65 million years of high-resolution global palaeoclimate data and with lineage diversification. Although previous research has revealed a correlation between contemporary climate and transitions to viviparity, here we aim to elucidate the relationship between past climatic conditions at the time points of evolutionary transitions to viviparity across squamate reptiles on a global scale.

| Identifying transitions from oviparity to viviparity
Ancestral state reconstruction analyses were performed based on previously published time-calibrated phylogenies (Zheng & Wiens, 2016) and parity data (Pyron & Burbrink, 2014) and carried out in R using the package corHMM (Beaulieu et al., 2012). The final data set consisted of 3,498 taxa with parity mode information mapped to the time-calibrated phylogeny. We tested alternative transition models for discrete characters using the asymmetric rates ('rayDisc') and hidden rates model ('corHMM) functions to infer ancestral states. In addition, we tested different options for setting root probabilities. Our final models included: (a) Asymmetric rate reconstruction with fixed root probabilities using the 'yang' option (Yang, 2006), (b) Asymmetric rate reconstruction fixing root to oviparous, (c) Asymmetric rate reconstruction with empirically estimated transition probabilities following Dollo's law of (almost) no reversals to oviparity and (d-h) hidden rates model (HMM) using a single to five rate categories (Beaulieu et al., 2012). From this set of models, we chose the one with the lowest AICc to proceed and extracted the nodes where a transition occurred (Table S1). Transitions were included in downstream analyses when the likelihood differential between the oviparous (coded as 1) ancestor and the viviparous (coded as 0) descendant was larger than 0.5 (= probability that a transition from oviparity to viviparity occurred in that node >50%). In addition, we performed separate analyses using less conservative differentials (0.3, 0.4, 0.45) and more conservative differentials (0.75, 0.8, 0.9). Reversals to oviparity from viviparity involve re-evolving egg-shells and were thought to be unlikely or impossible (Pincheira-Donoso et al., 2013;Shine, 2005;Sites et al., 2011); there are few reported convincing cases (Lynch & Wagner, 2010;Recknagel et al., 2018). However, the potential for reversals is still disputed and recent studies have proposed that reversals to oviparity have occurred more frequently but still rarely (Blackburn, 2015;Fenwick et al., 2012;King & Lee, 2015;Pyron & Burbrink, 2014). Therefore, we did not include potential transitions from viviparity to oviparity in our analysis.
We then tested whether transition time estimates for viviparity reflect the background rate of lineage diversification in squamates.
First, we explored the diversification pattern using a lineagethrough-time plot visualizing the number of squamate lineages against time using the recent time-calibrated phylogeny by Zheng and Wiens (2016) with the 'ltt.plot' function in ape using log transformed species numbers. Second, a binary trait was simulated 50 times based on that phylogeny with a transition rate to viviparity of 0.0015 and a much less likely transition rate to oviparity of 0.0001 using the 'rTraitDisc' function in ape (Paradis et al., 2004). The root was fixed to oviparous. These transition rates were inferred from empirically assessing a range of values (for oviparity to viviparity: 0.01, 0.005, 0.002, 0.0015, 0.001) and comparing the average number of transitions and proportion of viviparous species known in squamates (Pyron & Burbrink, 2014;Zheng & Wiens, 2016) (see Table S1 to compare the simulations to empirical data in number of transitions and proportion of viviparous species). After the trait simulation, 'corHMM' was used to reconstruct ancestral traits with three rate categories and the root fixed to oviparous as above. The ggtree package in R (Yu et al., 2017) was then used to infer binary trait transitions, based on a differential in trait likelihood larger than 50% and 80% between any parent and child node. Background rate and simulated transitions were compared with the empirical transitions to viviparity. The difference between simulated and empirical transitions per million year bin was calculated, normalizing the simulated transitions to the total empirical transitions with a likelihood differential of 50% and 80%, respectively. Lastly, we estimated the transition rate to viviparity from the number of transitions per million year time bin divided by the number of branching events in oviparous species that did not result in a transition to viviparity.

| Statistical analysis
To examine whether transitions to viviparity are associated with temperature, we used fossil oxygen isotope (δ 18 O) measurements from benthic oceanic foraminifera, a proxy for estimating palaeoclimatic temperature (Zachos et al., 2008). Measurements were averaged in bins of one million years. From these temperature estimates, we developed two different summaries: (1) the average temperature per million years (T) and (2) the relative temperature change as variance over a window of 10 million years (ΔT L ). The two temperature summaries were not correlated (Pearson correlation coefficients between −0.035 and 0.117; p > .1). Autocorrelation in temperature data should not be a problem for our analysis in principle, as there is no intrinsic causal relationship. However, we performed a test for autocorrelation for our best statistical model and found that there was no autocorrelation detected in the residuals using the 'acf' function in R. Accordingly, the regression analysis between lagged residuals (using a lag of one, which had the highest score in acf) was not significant (t = 1.794, p = .09). Therefore, we conclude that autocorrelation is not an issue in our statistical analyses.
We statistically tested for association between the two temperature parameters and the number of evolutionary transitions using R (R Core Team, 2015). Because the data were 'zero inflated' (58% of observations = 0, i.e. when the million year bin included no transitions to viviparity), we applied a generalized linear model with a negative binomial distribution (Crawley, 2007). Transition time estimates included the following: (a) conservative estimates of transition time (likelihood differential = 80%) inferred from the Zheng and Wiens (Zheng & Wiens, 2016) timetree, (b) less conservative estimates of transition time (likelihood differential = 50%), (c) the 'Dollo' model assuming highly unlikely transitions from viviparity to oviparity (likelihood differential = 80%) and (d) transitions inferred from simulations of binary trait evolution with conserved (80%) and (e) less conserved (50%) likelihood differentials, models correcting empirical transitions for simulated transition assuming (f) conserved (80%) and (g) less conserved (50%) likelihood differentials, and a transition rate measured as empirical transitions to viviparity per speciation node per million years with a (h) conserved and (i) less conserved likelihood differential. All models included an interaction term between T and ΔT L and were gradually simplified by AICc value and comparison to the null model. The model with the lowest AICc was chosen as the best model.

| Early Miocene burst in the origins of viviparity
We found that the HMM model with three rate categories best reconstructed the evolution of viviparity across squamates (Table S1).
The worst performing model was the 'Dollo' model, which assumed reversals to oviparity to be extremely unlikely. Most transitions from oviparity to viviparity in squamate reptiles occurred within the last 66 million years (Table 1), with only up to maximally three transitions occurring earlier than that. Depending on the likelihood differential of a transition, we found between 59 and 126 transitions in total (Table 1). The majority of transitions had a high likelihood differential between the oviparous parent node and the viviparous child node (e.g. likelihood differential cut-off at 50%: mean = 0.868, median = 0.925). With the least conservative cut-off likelihood differential of 30%, the mean was still relatively high with 0.755. Transitions tended to occur most frequently within the last 25 million years (>73% of all transitions, irrespective of cut-off). In general, using more conservative likelihood differentials led to transitions shifted towards more recent times (Table 1). The peak of transitions was between 11 and 16 mya, with 25.0%-28.8% of all transitions occurring in this time period.
On geological scales, origins of viviparity increased markedly in frequency from the late-Oligocene (around 25 mya) and decreased in the late-Miocene (around 5 mya) (Figures 1 and 2). The timing of transitions is broadly congruent across lizards and snakes, coinciding in the early/mid Miocene epoch; a peak around 10 to 25 mya in lizards (50% likelihood differential: 31 of 70 transitions; average = 21.3 mya) and snakes (50% likelihood differential: 14 of 29 transitions; average = 18.1 mya). There was no difference in transition rate between lizards and snakes: 70.9% of squamates are lizards, and the percentage of transitions occurring in lizards relative to all transitions within squamates was similar (HMM 50%: 71.0%; HMM 80%: 71.4%).
We found that the background rate of lineage diversification in squamates, based on the comprehensive phylogeny of 4,162 taxa (Zheng & Wiens, 2016), closely matched that of exponential growth through time during the last 100 million years (Figures 1 and S1).
Simulations of binary trait transitions generally matched the empirical values found for the number of viviparous species with regard to proportions: in simulations, on average 24.7% of squamate species were viviparous, which is close to the real data estimates of 19.8% of species (Zheng & Wiens, 2016). The larger average proportion of viviparous species was due to a few simulations in which transitions to viviparity had occurred very early in squamate evolution, rendering large clades of the phylogeny to be viviparous and reversing the pattern that oviparity is most common to a predominance of viviparity. The exclusion of those simulations reduces the proportion of viviparous species to 18.1%, matching the proportions of real data more closely.
Therefore, the simulations overall closely matched the empirical data matching the number of viviparous and oviparous species, but differed in when these transitions occurred (Table 1), and in particular regarding the number of transitions in early squamate evolution and transitions during the past 10-25 mya (Figures 1, S1 and S2). In simulations, between four (80% likelihood differential) and six (50% likelihood differential) transitions to viviparity occurred earlier than 65 million years ago, whereas empirical transitions varied between one and three in that time, respectively.
In addition, the empirically inferred transitions to viviparity differ in slope ( Figure 1) and shape ( Figure S1) from the binary trait simulations, with a lower number of early transitions to viviparity, and a steeper slope representing an increase approximately 10 to 23 million years ago (Figure 1, Figure S1). This is reflected by a significant difference between empirical and simulated transitions in their distribution (Kolmogorov-Smirnov test: p < .001). Distributions of empirical transitions were independent of using alternative likelihood differential cut-offs (Kolmogorov-Smirnov test: p < .94). Note: Different likelihood thresholds are applied for counting transitions (e.g. 50% threshold: differential between a parent and child node in likelihood of being viviparous is more than 0.5). Using lower thresholds (= including less likely transitions) leads to lower mean likelihoods across all transitions. The time estimate at which 50% of all transitions occurred (t 50% ) and 75% (t 75% ) are recorded, counting from present to past.

TA B L E 1 Empirical and simulated transitions to viviparity across squamates
In summary, based on the differences between empirical and simulated transition times, we found that empirical transitions to viviparity occurred less frequently during early squamate evolution, and more frequently than expected during the Oligocene between 10 and 25 million years ago ( Figure S2).

| Origins of viviparity during cold and stable temperatures?
We found that the transitions from oviparity to viviparity identified in the time-calibrated phylogeny were statistically associated with relatively cold and stable global climatic conditions (Figure 2;  (Figure 2; Table 2). This is consistent with cold climate (low temperatures) and environmental stability (low variance in long-term temperatures) associated with the evolution of viviparity. In analyses using the average mid-branch point instead of the timepoint of diversification as the viviparity transition point, cold temperature remained a significant association (not shown). However, simulations of the binary trait based on the squamate phylogeny showed a similar pattern, with associations in the same direction (model 4 and 5,  Figure S2). Indeed, when controlling for the simulated trait, empirical transitions to viviparity were no longer associated with any climatic variable (model 6 and 7, Table 2; Figure S2). Moreover,

| Palaeoclimate and origins of viviparity
Our results find support for all three aspects of our nonexclusive hypotheses: cold climate, stability and background diversification rates. We found that transitions from oviparity to viviparity across the squamate tree occurred most frequently during the past 10 to 25 million years. The most pronounced increase in transition frequency relative to background diversification occurred in the late-Oligocene and decreased again in the mid-Miocene, peaking at around 11 to 16 million years ago, in the early-/mid-Miocene ( Figure S2).
Our findings are consistent with the cold-climate hypothesis because we found that transitions to viviparity were statistically associated with historically relatively cold climate ( Figure 2; Table 2). This is in line with recent studies linking the prevalence of contemporaneous viviparous squamates to colder climates (Esquerré et al., 2019;Feldman et al., 2015;Pincheira-Donoso et al., 2017;Pyron & Burbrink, 2014;Scharf et al., 2015).
Although here conducted at palaeoclimatic and global scales, our results broadly agree with previous, lineage-specific studies that found a connection between transitions to viviparity and historically cold climatic conditions. For example, in vipers transitions to viviparity have been associated with global cooling during the Cenozoic period, and in the Oligocene speciation rates for viviparous species increased substantially relative to oviparous species F I G U R E 2 Frequency of transitions from oviparity to viviparity and palaeoclimate. (a) Empirically estimated transitions from oviparity to viviparity in squamates (dark grey) and the ratio of oxygen isotope (δ18O) as proxy for global mean temperature (light grey) are displayed per million years from 65 million years ago to present. Temperature (y-axis) and number of transitions (z-axis) are shown as smoothed lines (span[λ] = 0.25). Oxygen isotopes from (Zachos et al., 2008) were used to determine palaeotemperature (negative relationship). Geological epochs are indicated below. In (b), transitions of a simulated binary trait are compared to global temperature estimates through time. Simulated transition time estimates include mean values for 50 independent replicates (Lynch, 2009). In iguanid lizards, transitions to viviparity occurred within the last 65 million years but mostly during and following the cooling during the Eocene (35.2-52.6 Mya), predating the Pleistocene-Pliocene glaciations of the last 4 million years (Schulte & Moreno-Roark, 2010). In phrynosomatid lizards, transitions to viviparity were found to be correlated with cold summer temperatures (Lambert & Wiens, 2013).
A challenge with inferring the discrete timing of past evolutionary transitions is that the transition could theoretically have happened at any point on the branch leading to the species with the derived acquired trait (Duchen et al., 2017). We note that if viviparity evolved on average at the midpoint of the branch leading to the viviparous species rather than at the divergence point, this shifts origins of viviparity to be more recent, leading to a stronger association with cold climate.
Our results are also consistent with origins of viviparity being associated with periods of stable climate. Compared with cold climate, less attention has been given to the importance of environmental stability, especially at larger evolutionary time scales (i.e. speciation/ diversification times over millions of years) (Svenning et al., 2015). Climate reconstructions show that there was limited temperature variation from the Oligocene to the mid-Miocene (Figure 2), which was generally characterized by a relatively stable climate (Flynn et al., 1995). This period of stability was followed by a dramatic drop in temperature from the mid-Miocene. Our models of transition rates corrected for a simulated trait showed larger transition frequencies during this time; however, there was no significant relationship across the whole phylogeny and time ( Figure S2, Table 2). Although compatible with extrinsic influences of cold and stable climate, we find that a simulated binary trait across the squamate phylogeny showed similar statistical associations. Indeed, when the empirically estimated transitions to viviparity were corrected for the number of diversification events or for a binary simulated trait, the association with climate was no longer present. This suggests that due to their natural correlation, these two processes cannot be disentangled.

| Lineage diversification and transitions to viviparity
Our study shows that origins of viviparity are consistent with a role of stable and colder climate, but binary trait simulations showed a similar association with colder climates. Transitions to viviparity coincided with diversification of squamate reptiles in general, casting doubt on whether climate alone is enough to either explain or manifest the evolution of viviparity. Our result is consistent with a previous study concluding that transitions to viviparity had not occurred more frequently than expected during cold epochs within the last 65 million years (King & Lee, 2015). One reason for the association between origins of viviparity and diversification could be that transitions to viviparity favour diversification. This has been previously suggested for squamate reptiles (Lambert & Wiens, 2013), although extinction rates were also higher for viviparous species (King & Lee, 2015;Pyron & Burbrink, 2014 For each of the nine tested models, only the model with optimal performance (by AICc) is shown. Explanatory variables are the two climatic parameters average temperature per million year (T) and long-term change in temperature (ΔT L ) and their interaction.
Therefore, when coupled with a lineage having a genetic propensity to evolve phenotypic diversity, lineage diversification may itself intrinsically lead to trait diversification given ecological opportunity (Blount et al., 2012;Elmer et al., 2014;Lovette et al., 2002;Reaney et al., 2018). For example, birds have repeatedly evolved adaptation to altitude that has resulted in speciation. Molecular convergence underpinning this adaptation has been found to be high but dependent on the genetic background of the lineages (Natarajan et al., 2016).
A similar association between molecular convergence and phylogenetic signal has been observed in other studies (Agrawal, 2017;Christin et al., 2015;Dobler et al., 2012;Ng & Smith, 2016;Ujvari et al., 2015). Therefore, mechanistically the association between lineage diversification and convergent trait diversification may underlie a genetic predisposition for the evolution of certain traits. We speculate that such a combination of factors play a role in origins of squamate viviparity, though its validation would require future research.

| CON CLUS ION
In summary, our study links palaeoclimatic conditions to transitional events from oviparity to viviparity in the evolutionary past on a broad phylogenetic scale. Importantly, we find that transitions to viviparity in squamates closely followed background diversification rate, meaning that diversification and trait evolution cannot be separated. Consequently, our findings showed that viviparity evolved during stable and long-lasting cold climatic conditions, this is correlated with intrinsic lineage diversification, and the effects cannot be separated. In light of other recent studies (Ma et al., 2018;Pincheira-Donoso et al., 2017), the ultimate causes for transitions from oviparity to viviparity in squamate reptiles appear to be complex and encompass a wide range of environmental and intrinsic factors.

More generally, this is in line with other famous examples of trait
diversification-such as seen in cichlids or Darwin's finches-with an implication of intrinsic and extrinsic factors shaping the diversification process (Stroud & Losos, 2016;Tokita et al., 2017;Wagner et al., 2012;Xiong et al., 2021).

ACK N OWLED G EM ENTS
We thank J. Zachos for providing oxygen isotope data. We thank R. Biek, W. Harvey and A. Crawford for project discussion, and anonymous referees for comments that improved the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The paper analyses publicly available data.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.13886.