Ancient tropical extinctions at high latitudes contributed to the latitudinal diversity gradient *

Global biodiversity currently peaks at the equator and decreases toward the poles. Growing fossil evidence suggest this hump‐shaped latitudinal diversity gradient (LDG) has not been persistent through time, with similar diversity across latitudes flattening out the LDG during past greenhouse periods. However, when and how diversity declined at high latitudes to generate the modern LDG remains an open question. Although diversity‐loss scenarios have been proposed, they remain mostly undemonstrated. We outline the “asymmetric gradient of extinction and dispersal” framework that contextualizes previous ideas behind the LDG under a time‐variable scenario. Using phylogenies and fossils of Testudines, Crocodilia, and Lepidosauria, we find that the hump‐shaped LDG could be explained by (1) disproportionate extinctions of high‐latitude tropical‐adapted clades when climate transitioned from greenhouse to icehouse, and (2) equator‐ward biotic dispersals tracking their climatic preferences when tropical biomes became restricted to the equator. Conversely, equivalent diversification rates across latitudes can account for the formation of an ancient flat LDG. The inclusion of fossils in macroevolutionary studies allows revealing time‐dependent extinction rates hardly detectable from phylogenies only. This study underscores that the prevailing evolutionary processes generating the LDG during greenhouses differed from those operating during icehouses.

throughout the Phanerozoic (the last 540 million years), even if the gradient was sometimes shallower (Mittelbach et al. 2007), based on published fossil record studies (Crame 2001;Alroy et al. 2008). However, the methodological limitations of fossil sampling have called this conclusion into question. Analyses controlling for sampling bias have suggested that, for many groups, the LDG was less marked in the past than it is today (i.e., with similar species diversity across latitudes) or even developed a paleotemperate peak during some periods (see Mannion et al. 2014 for a review). This sampling-corrected flattened LDG in deep time has been demonstrated for nonavian dinosaurs (Mannion et al. 2012), mammals (Rose et al. 2011;Marcot et al. 2016), birds (Saupe et al. 2019a), tetrapods (Brocklehurst et al. 2017),

Ice-free temperature (ºC)
Earth and biome history LDG pattern insects (Archibald et al. 2010, Archibald et al. 2013Labandeira and Currano 2013), brachiopods (Krug and Patzkowsky 2007;Powell 2007;Powell et al. 2012), bivalves (Crame 2000(Crame , 2020, coral reefs (Leprieur et al. 2016), foraminifers (Fenton et al. 2016), crocodiles (Mannion et al. 2015), turtles (Nicholson et al. 2015, or plants (Coiffard and Gomez 2012;Peralta-Medina and Falcon-Lang 2012;Shiono et al. 2018). The pattern emerging from fossil studies also suggests that steep LDGs, such as that currently observed, have been restricted to the relatively small number of short icehouse periods during the Earth's history: the Ordovician/Silurian, the Carbonifer-ous/Permian, and the Neogene. Most of the Phanerozoic has instead been characterized by warm greenhouse climates associated with a flatter LDG (Mannion et al. 2014;Marcot et al. 2016). Hence, the last change in the shape of the LDG hypothetically occurred during the greenhouse to icehouse transition of the Cenozoic, in the last 66 million years (Fig. 1).
Our perception on how diversity was latitudinally distributed in the past underlies LDG interpretations. Up to now, most studies have assumed the equator is the source of world diversity (Donoghue 2008;Jablonski et al. 2013) and diversity was always lower in the Holarctic than in the Equator (i.e., a steep LDG persisted in the deep time). The current shape of the LDG thus results from lower paces of diversity accumulation in the Holarctic than at the equator through time. Slower diversity accumulation in the Holarctic (hereafter "H") have been explained either by greater tropical diversification and limited dispersal out of the equatorial region (hereafter "equator" or "E") (r e > r h ; d eh > d he ) (Latham and Ricklefs 1993;Wiens and Donoghue 2004;Mittelbach et al. 2007;Rolland et al. 2014), or by high rates of turnover in the Holarctic, that is, similar high speciation (λ) and extinction (μ) rates (λ h ≈ μ h ; Table 1). This has been widely shown throughout the tree of life: amphibians (Wiens 2007;Pyron and Wiens 2013), birds (Cardillo et al. 2005;Ricklefs 2006;Weir and Schluter 2007), butterflies (Condamine et al. 2012), plants (Leslie et al. 2012;Kerkhoff et al. 2014), fishes (Siqueira et al. 2016), mammals (Weir and Schluter 2007;Rolland et al. 2014), or lepidosaurs (Pyron 2014). Contrary to the "slow Holarctic diversity accumulation" hypotheses, a scenario assuming the Holarctic was also a source of tropical diversity that flattened the LDG in the deep past could be considered. The current steep shape of the LDG would thus result from diversity lost in the Holarctic through evolutionary time, referred to as the "Holarctic diversity loss" hypothesis. The recent fossil investigations showing, for many lineages, similar diversity levels across latitudes in the past lend support to this scenario, suggesting we do not necessarily need to explain why diversity accumulated at slower rates in the Holarctic through time, but the question being how and when diversity was lost at high latitudes, giving rise to the current shape of the LDG (Mannion et al. 2014)?
Diversity losses in the Holarctic have been traditionally considered to underlie the LDG. They were initially attributed to Pleistocene glaciations (Martin and Klein 1989), but this hypothesis can be called into question because the LDG substantially predates the Pleistocene (Mittelbach et al. 2007). More ancient extinctions have also been considered (Latham and Ricklefs 1993;Markwick 1998;Roy and Pandolfi 2005;Hawkins et al. 2006;Weir and Schluter 2007;Dunn et al. 2009;Eiserhardt et al. 2015;Pulido-Santacruz and Weir 2016). For example, recent studies suggested the avian LDG resulted from the differential extirpation of older warm-adapted clades from the temperate regions newly formed in the Neogene (Hawkins et al. 2006;Pulido-Santacruz and Weir 2016;Saupe et al. 2019a). Pyron (2014) suggested that higher temperate extinction represents a dominant force for the origin of LDG in lepidosaurs. Wiens and Donoghue (2004) also proposed that range contractions (d eh < d he ) followed climate cooling.
Unfortunately, using phylogenies of extant taxa alone, studies on the LDG have not clearly demonstrated diversity losses in the Holarctic but instead high regional turnover (λ ≈ μ) (Cardillo et al. 2005;Weir and Schluter 2007;Condamine et al. 2012;Leslie et al. 2012;Pyron and Wiens 2013;Rolland et al. 2014Rolland et al. , 2015Rabosky et al. 2018). Nonetheless, high turnover in the Holarctic can only explain a slow accumulation of lineages, with one fauna being replaced by another, but does not explain diversity decline (i.e., a reduction in the net number of species) as elevated speciation rates counterbalance the effect of extinction. Diversity declines occur when extinction exceeds speciation (μ > λ), resulting in negative net diversification rates (r = λ − μ; r < 0). Accordingly, "diversity loss" hypotheses differ from "high turnover" scenarios.
The perceived difficulty for inferring negative diversification rates from phylogenetic data of extant species (Rabosky 2010;Burin et al. 2019) and the assumption that diversity levels always remained lower in the Holarctic than at the equator have resulted in "diversity loss" hypotheses being repeatedly proposed but seldom demonstrated using phylogenetic data (Pulido-Santacruz and Weir 2016). Meanwhile, numerous fossil investigations have detected signatures of extinction and diversity loss in the Northern Hemisphere. For instance, Archibald et al. (2010Archibald et al. ( , 2013 sampled insect diversity at an Eocene site in Canada, and in present-day temperate Massachusetts (USA) and tropical sites of Costa Rica. Insect diversity was higher at the Eocene paleotropical site than the modern temperate locality, and comparable to the modern-day tropical locality, suggesting that post-Eocene Nearctic insects have suffered great levels of extinction. This pattern is consistent with other studies on various taxonomic groups, including birds (Mayr 2016), invertebrates (Wilf et al. 2005), mammals (Blois and Hadly 2009;Rose et al. 2011;Marcot et al. 2016), and plants (Frederiksen 1988;Smith et al. 2012;Xing et al. 2014). Fossil studies are, however, generally restricted to a reduced geographic and temporal scale, which makes difficult to extrapolate local inferences of extinction in the global context of the LDG.
Here, we use comparative methods for both phylogenies and fossils to estimate the evolutionary processes behind the LDG of Testudines, Crocodilia, and Lepidosauria, and test the alternative predictions of two hypotheses: "slow diversity accumulation" (d eh > d he ; μ h ≤ λ h ) vs. "diversity loss" (d eh < d he ; μ h > λ h ) in the Holarctic. To this end, we propose to include a temporal component to study the LDG in which prevailing speciation, extinction, and dispersal dynamics may change between warmand cold-time intervals. The extant Crocodilia and Lepidosauria comprise mostly tropical-adapted species with a classic LDG pattern as shown by diversity peaks at equatorial latitudes (Markwick 1998;Pyron 2014). We investigate the origin of the LDG in subtropical taxa as well, by extending the study to Testudines, which display a hump-shaped gradient of diversity centered on subtropical latitudes (10°S-30°N) (Angielczyk et al. 2015). By contrast, the paleolatitudinal distribution of turtles was concentrated in the Holarctic (30-60°N) during the Cretaceous (Nicholson et al. 2015). All these lineages Main evolutionary hypotheses can be classified according to three criteria: (1) the mechanisms behind regional differences in species richness, including explanations based on evolutionary time (t), on dispersal (d), and on diversification (the composite value r = λ − μ). Explanations based on evolutionary time assume most groups originated in the tropics and had less time to diversify in the temperate regions (Stephens and Wiens 2003). Hypotheses focusing on the role of geographic movements (d), include the Tropical niche conservatism model assuming that most groups originated in the tropics and the LDG results from limited dispersal to the temperate regions (Farrell et al. 1992;Wiens and Donoghue 2004). The "Into the tropics" model assumes instead the LDG results from dispersals toward the equator (Condamine et al. 2012). Hypotheses that emphasize the LDG is generated by regional differences in net diversification rates, being higher in the tropics (Ricklefs 2006;Wiens 2007) assume the outstanding tropical diversity is the outcome of higher rates of speciation in the tropics than in the extra-tropical regions ( λ t > λ e ) under the "cradle of diversity," and/or result from lower rates of extinction (μ t < μ e ) under the "museum of diversity" (Stebbins 1974). The LDG could also result from higher turnover rates (i.e., higher λ and μ) in the Holarctic (Weir and Schluter 2007). Diversification and dispersal hypotheses are not mutually exclusive. In the "out of the tropics" model, the tropics are regarded as both cradle and museum, with lineages preferentially originating in the tropics and expanding into high latitudes (Jablonski et al. 2006). Hypotheses could also be classified according to (2) the rate at which processes acted through time; most studies assumed evolutionary processes acted constantly. In the AGED model, which extends the Tropical niche conservatism hypothesis (Wiens and Donoghue 2004), diversification/dispersal parameters vary for each temporal interval (the transition from icehouse to greenhouse, and vice versa).
(3) Hypotheses can be classified according to the origin of tropical diversity: "Lower Holarctic species accumulation" hypotheses assume the equator is the source of world diversity and species accumulated at slower rates on higher latitudes. Conversely, "Holarctic diversity loss" hypotheses assume the Holarctic was also a source of diversity that was lost during evolutionary history. are ancient and likely experienced climatic transitions during the early Cenozoic (Markwick 1998;Pyron 2014;Angielczyk et al. 2015;Mannion et al. 2015;Nicholson et al. 2015). They show contrasting patterns of species richness: turtles and crocodiles are species poor (350 and 25 species, respectively), whereas lepidosaurs include a large number of species (10,000 + species), and all have a rich fossil record extending back to the Triassic (Early Cretaceous for crocodiles), providing information about the variation of latitudinal species richness accumulation during their evolutionary history.

MOLECULAR PHYLOGENIES AND THE FOSSIL RECORD
A time-calibrated phylogeny for turtles (order Testudines) was obtained from Jaffe et al. (2011), including 233 species. We preferred this phylogeny over other more recent and slightly better sampled trees (Rodrigues and Diniz-Filho 2016) because the divergence time estimates are more consistent with recent estimates based on genomic datasets (Pereira et al. 2017;Shaffer et al. 2017). For scaled lizards (Lepidosauria), we retrieved the most comprehensive dated tree available, including 4161 species (Pyron 2014), and a complete phylogeny was obtained for crocodiles (Crocodilia) (Oaks 2011). Fossil occurrences were downloaded from the Paleobiology Database (https://paleobiodb.org/#/, last accessed 25 October 2017). We reduced potential biases in the taxonomic assignation of turtle, crocodile, and lepidosaur fossils, by compiling occurrence data at the genus level. We further cleaned the fossil datasets by checking for synonymies between taxa and for assignment to a particular genus or family on the basis of published results.

RATES WITH PHYLOGENIES
We ensured comparability with previous LDG studies and investigated possible differences between Holarctic and equatorial regions by combining the turtle and lepidosaur phylogenies with distributional data to fit trait-dependent diversification models in the Binary-state Speciation and Extinction (BiSSE) model (Maddison et al. 2007). We accounted for incomplete taxon sampling in the form of trait-specific global sampling fraction (FitzJohn et al. 2009). We did not use the geographic-state speciation and extinction model (Goldberg et al. 2011), which is appropriate for dealing with widespread species, because most of the species in our datasets were endemic to the Holarctic or equatorial regions (see below), and, for a character state to be considered in SSE models, it must account for at least 10% of the total diversity (Davis et al. 2013). We did not apply the BiSSE model to crocodiles, because simulations have shown that trees containing fewer than 300 species may have too weak a phylogenetic signal to generate sufficient statistical power (Davis et al. 2013).
We initially implemented a constant-rate BiSSE model in the R-package diversitree 0.9-7 (FitzJohn 2012) with six parameters: two speciation rates, one associated with the Holarctic ("H," λ H ) and the other with other equatorial and subtropical regions ("equator" or "E," λ E ), two extinction rates associated with the Holarctic (μ H ) and the equator (μ E ), and two transition rates (dispersal), one for the Holarctic to equator direction (q HE ), and the other for the equator to Holarctic direction (q EH ). We categorized each species as living in the equator or the Holarctic according to the threshold latitudes defining the tropics (23.4°N and 23.4°S). According to our distribution data, 84% of extant species of turtles, lepidosaurs, and crocodiles lives in the tropics, 15% in temperate regions, and 1% span both biomes. For turtles, there were 239 tropical species, 84 temperate, and six spanning both biomes (seven were marine). For lepidosaurs, 7955 tropical, 1337 temperate, and 124 spanning both biomes. Crocodiles had 23 tropical and two temperate species (Tables S1-S3).
We then used the time-dependent BiSSE (BiSSE.td) model, in which speciation, extinction, and dispersal rates are allowed to vary between regions and to change after shift times. We introduced two shift times to model different diversification dynamics among greenhouse, transitional, and icehouse periods. We assumed that a global warm tropical-like climate dominated the world from the origin of the clades until 51 Million years ago (Ma) (corresponding to the temperature peak of the Cenozoic). Thereafter, the climate progressively cooled until 23 Ma (the transitional period), when the climate definitively shifted to a temperate-like biome in the Holarctic (Ziegler et al. 2003;Zachos et al. 2008;Morley 2011). The climatic transition in the Cenozoic may have different temporal boundaries, with potential effects on the results. We thus applied the same model but with different combinations of shift times (we tested 51/66 Ma and 23/34 Ma for the upper and lower bounds of the climatic transition). We used a Markov Chain Monte Carlo (MCMC) approach to investigate the credibility intervals of the parameter estimates, with an exponential prior 1/(2r) and initiated the chain with the parameters obtained by maximum likelihood. We ran 10,000 MCMC steps, with a burn-in of 10%.

RATES WITH FOSSILS
We analyzed the three fossil records using a Bayesian model implemented in PyRate (Silvestro et al. 2019) for simultaneous inference of the temporal dynamics of origination and extinction, and of preservation rates (Silvestro et al. 2014). The turtle fossil dataset contains 4084 occurrences for 420 genera (65 extant and 355 extinct; Table S4). The lepidosaur fossil dataset comprises Table S5). The crocodile fossil dataset includes 1596 occurrences for 121 genera (nine extant and 112 extinct; Table S6). In this analysis, the preservation process is used to infer the individual origination and extinction times of each taxon from all fossil occurrences and an estimated preservation rate expressed as expected occurrences per taxon per million years. We followed a birth-death shift approach (Silvestro et al. 2015), also known as the Bayesian skyline model, and used a homogeneous Poisson process of preservation (-mHPP option). We used default settings and accounted for the variation of preservation rates across taxa, using a Gamma model with gamma-distributed rate heterogeneity (-mG option) and four rate categories to discretize the gamma distribution.
Given the large number of occurrences analyzed and the vast timescale considered, we dissected the birth-death process into time intervals, and estimated origination and extinction rates within these intervals. In one set of analyses, we defined time intervals using the geological epochs of the stratigraphic timescale (Ogg et al. 2004). In another set of analyses, we defined time intervals according to the major climatic periods characterizing the Cenozoic: the greenhouse world, the climatic transition, and the icehouse world, testing two alternatively boundaries for the climatic transition in the Cenozoic (34/23 Ma). We adopted this solution as an alternative to the algorithms implemented in the original PyRate software for joint estimation of the number of rate shifts and the times at which origination and extinction shift (Silvestro et al. 2014). The estimation of origination and extinction rates within fixed time intervals improved the mixing of the MCMC and made it possible to obtain an overview of the general trends in rate variation over a long timescale. Both the preservation and birth-death processes were modeled in continuous time but without being based on boundary crossings. One potential problem when fixing the number of rate shifts a priori is over-parameterization. The Bayesian skyline model overcame this problem by assuming that the rates of origination and extinction belonged to two families of parameters following a common prior distribution, with parameters estimated from the data with hyper-priors (Gelman 2004).
We ran PyRate for 10 million MCMC generations on each of the 10 randomly replicated datasets. We monitored chain mixing and effective sample sizes by examining the log files in Tracer 1.7 (Rambaut et al. 2018). After excluding the first 20% of the samples as a burn-in, we combined the posterior estimates of the origination and extinction rates across all replicates to generate plots of the change in rate over time. The rates of two adjacent intervals were considered significantly different if the mean of one lay outside the 95% credibility interval of the other, and vice versa.
In the context of the LDG, we performed additional analyses with different subsets of fossils, to separate the speciation, extinction, and preservation signals of different geographic regions (equator or Holarctic) and ecological conditions (temperate or tropical). For example, for turtles, we split the global fossil dataset into four subsets: one for the fossil genera occurring at the equator (429 occurrences), one for the fossils occurring in the Holarctic (3568 occurrences), one for the fossil genera considered to be adapted to temperate conditions (993 occurrences), and one for the fossils considered to be adapted to tropical conditions (2996 occurrences). We excluded the few fossil occurrences for the southern regions of the South Hemisphere (about 180) only in subset analyses, as they were poorly represented in our dataset. Note that a given fossil can be present in both the 'Holarctic" and "tropical" datasets. We encoded tropical/temperate preferences by considering macro-conditions in the Holarctic to be paratropical until the end of the Eocene, as previously reported (Ziegler et al. 2003;Morley 2011). We also assumed that taxa inhabiting the warm Holarctic were adapted to tropical-like conditions (i.e., a high global temperature, indicating probable adaptation to tropical climates). After the late Eocene, we categorized each species as living in the temperate biome as defined above. With these datasets, we reproduced the same PyRate analyses as for the whole dataset.
Finally, we examined the apparent incongruence between paleontological and phylogenetic estimates of diversification (see Results) using the birth-death chronospecies (BDC) model in PyRate (Silvestro et al. 2018). This model allows examining whether alternative speciation modes described for the fossil record, that is, budding, bifurcation, or anagenesis (Silvestro et al. 2018), are responsible for driving incongruences between fossil and phylogenetic estimates. We compared an (i) Equal rates model where diversification parameters estimated with stratigraphic data (λ * and μ * ) are the same to those estimated with phylogenetic data (λ and μ), that is, λ * = λ, μ * = μ; (ii) a Compatible model where parameters differ, but the differences could be explained by differences in speciation mode. In this model, λ * , λ, μ, and μ * are constrained such that λ * − λ = μ * − μ (i.e., equal net diversification rates) and λ * ≥ λ; (iii) an Incompatible model, where parameters λ, λ * , μ, and μ * are allowed to take any value and thus differences in λ and λ * , as well as μ and μ * , cannot be explained by differences in speciation mode.
We estimated λ, μ, λ * , and μ * simultaneously for the Testudines, Crocodilia, and Lepidosauria phylogenies pruned to the genus level, and the corresponding fossil data using maximumlikelihood optimization and assuming constant diversification rates through time. To assess support for the BDC model, we applied a likelihood ratio test (Silvestro et al. 2018), comparing which of the equal, compatible, or incompatible rates models are supported by our data.
Because the evolutionary histories of Testudines, Lepidosauria, and Crocodilia are long and exhibit great amount of temporal heterogeneity in both speciation and extinction rates (see Results), we also implemented a Bayesian skyline model with rate shifts at the climatic epoch boundaries defined above. We ran 10 million MCMC iterations to obtain joint posterior distributions of the stratigraphic and phylogenetic rates. We used the joint posterior samples of λ * , μ * , λ, and μ obtained under the incompatible rates model to verify the conditions predicted by the compatible BDC model (i.e., λ * − λ = μ * − μ and λ * ≥ λ) and assess the support for each model as suggested in previous studies (Silvestro et al. 2018).

WITH PHYLOGENIES AND FOSSILS
We performed biogeographic analyses with the parametric likelihood method DEC  using the fast C++ version (Smith 2009) (https://github.com/rhr/lagrange-cpp). Turtle, lepidosaur, and crocodile species distributions were obtained from online databases (www.iucnredlist.org and www.reptiledatabase.org). We also chose 23.4°N and 23.4°S as the threshold latitudes defining the tropics, and categorized each species as living in the Holarctic, in the southern temperate regions, or in the equatorial tropics and subtropical regions. We considered that all ranges comprising three areas could be considered an ancestral state (maxareas = 3).
We set up three different DEC analyses. We first ran DEC with no particular constraints, using only the distribution of extant species. We then performed DEC analyses including fossil information in the form of "fossil constraints" at certain nodes, according to the range of distribution of fossil occurrences assigned to a particular taxon during the relevant time frame. For example, the crown age of Carettochelyidae (Testudines) dates back to the Late Jurassic (150 Ma), and we set a geographic constraint on this node reflecting the distribution of all the Late Jurassic fossils attributed to Carettochelyidae. Similarly, for the origin of turtles (210 Ma), distribution constraints represent the range of Late Triassic fossils assigned to turtles.
We followed two different approaches to include fossil distributions. First, we used a soft fossil constraints (SFC) approach to incorporate fossil data into the anagenetic component of the likelihood framework. The direct impact of a given fossil is limited to the particular branch to which it has been assigned, although it may indirectly influence other branches. The inclusion of a fossil conditions the estimated geographic-transition probability matrix for that branch by imposing a spatiotemporal constraint on the simulation process. Only the simulations resulting in a geographic range including the area of fossil occurrence contribute to the geographic-range transition probability matrix for the branch concerned; simulations not meeting this constraint are discarded (Moore et al. 2008). This was achieved with existing functions in the C++ version of Lagrange, using the command "fossil." We consider this to be a "soft" constraint, because other areas different from that in which the fossil was found could be included in the ancestral states. In some cases, the SFC model may still overlook known fossil information. Second, we then implemented a hard fossil constraints (HFC) approach, where the estimation of ancestral areas was fixed to the location of fossils. For HFC, we used the command "fixnode." By fixing nodes to the distribution area of fossils, we assume fossil occurrences reflect the distribution of the ancestors. This is a strong assumption, but it makes it possible to recover all fossil ranges in the ancestral estimations. The real scenario probably lies somewhere between SFC and HFC inferences.
We then compared the timing and number of range extinction and dispersal events inferred with the three different biogeographic approaches through time. In DEC, extinction (range contraction) and dispersal (range expansion) events are modeled as stochastic processes occurring along the branches of the tree (Ree and Sanmartin 2009), with the probability of any extinction/dispersal event to occur being constant along the entire length of the branch. We therefore estimated the periods at which range extinction and dispersal occurred by dividing the phylogeny into intervals of 25 million years and calculating the number of branches for which extinction/dispersal was inferred crossing a particular time interval (the same branch could cross two continuous intervals).

EQUATOR?
Under the time-constant BiSSE model, net diversification rates for turtles were higher in the Holarctic than at the equator (Fig. S1A), but this difference was not significant, and rates of dispersal "into the equator" were 10 times higher than those "out of the equator." For lepidosaurs, a similar dispersal pattern was recovered, but net diversification rates were significantly higher in the equator (Fig. S1B). The time-variable BiSSE models, with shift times at 51 and 23 Ma, indicated that speciation and extinction rates of turtles were similar in the Holarctic and at the equator until the icehouse period, when Holarctic speciation increased. For lepidosaurs, speciation was lower in the Holarctic and extinction higher until the icehouse period, when Holarctic speciation increased and extinction decreased (Fig. S2). Dispersal was symmetric between regions (into the equator = out of the equator) during greenhouse periods, and asymmetric (into the equator > out of the equator) during the climatic transition and icehouse period. The same patterns were obtained assuming different combinations of shift times ( Crocodilia loss of diversity

FOSSIL-BASED DIVERSIFICATION ANALYSES: EVIDENCE FOR ANCIENT TROPICAL EXTINCTIONS?
We first inferred global diversification dynamics by analyzing the fossil datasets as a whole. For turtles, origination rates peaked during the Jurassic, subsequently decreasing until the present. Extinction rates were generally low and constant during the Mesozoic, but increased during the Jurassic and the Paleogene, resulting in negative net diversification during the Paleogene (Fig. 2). For lepidosaurs, origination rates peaked in the Jurassic and Late Cretaceous, whereas extinction increased steadily until the Late Cretaceous. In the Paleogene, net diversification approached zero, suggesting a high turnover. Crocodile origination peaked in the Early Cretaceous, subsequently decreasing toward the present, and extinction rates were generally low and constant. We also identified diversity losses in the Paleogene extending to the present, suggesting that crocodiles are still in a phase of declining diversity (Fig. 2).
Additional analyses with different subsets of the three fossil datasets separating origination and extinction signals between geographic regions (equator or Holarctic) and ecological conditions (temperate or tropical) showed that the diversity losses experienced by turtles and crocodiles during the Paleogene were mostly attributable to species living in the Holarctic and under tropical conditions (Figs. 3 and 4). The global diversity loss inferred for crocodiles during the Neogene was attributed to taxa living in both the Holarctic and equatorial regions (adapted to temperate and tropical conditions, respectively), providing further support for the hypothesis that this whole group is in decline.
For all groups, temperate taxa have been estimated to have high rates of diversification during the Oligocene, but lower rates during the Neogene. For the equatorial datasets, extinction and origination rates decreased over time, resulting in constant net diversification rates (except for lepidosaurs, which displayed a decrease in diversification during the Paleogene, followed by an increase during the Neogene). The same patterns were obtained for analyses with temporal boundaries defined according to the main climatic periods of the Cenozoic (Figs. S11-S13), and assuming a different combination of shift times to represent uncertainty on the temporal boundaries of the largest climatic oscillations (23/34 Ma; Figs. S5-S10). PyRate analyses also estimated similar preservation rates in the Holarctic and the equator for turtle and crocodile fossils, whereas preservation rates for lepidosaurs are much higher in the Holarctic than in the equator (Table 2). Finally, we found that crocodiles conform to the expectations of the compatible rates BDC model under a constant rate assumption and likelihood threshold values of 0.95 (they conform to the equal rates model at 0.99 threshold). For Testudines and Lepidosauria, the time-constant BDC model was rejected in favor for the incompatible rates model in all cases (Table S7). Relaxing the assumption of constant rates resulted in strong support for the compatible rates BDC model in turtles (P > 0.01 or 0.05) and the incompatible rates model in lepidosaurs (P < 0.01). Phylogenetic and fossil estimates of diversification for crocodiles are not significantly different during the first two time intervals (from the origin to 23 Ma) supporting the equal rates model (Fig. 5).

PREFERENTIALLY ORIGINATE AT THE EQUATOR?
Based on the unconstrained DEC analysis, we inferred an equatorial distribution for the deepest nodes for the turtles and lepidosaurs, whence these lineages colonized the other regions Ice. Tran.

Testudines
Lepidosaurs (Figs. 6A and S14). Crocodile ancestors were found to have been widespread during the Cretaceous, with an early vicariant speciation event separating Alligator in the Holarctic from the other genera of Alligatoridae in equatorial regions (Fig. S15). In contrast, our biogeographic estimates based on extant and fossil data yielded very different histories for the three groups (turtles: Figs. 6B and S16; lepidosaurs: Figs. S17-S18; and crocodiles: Figs. S19-S20), with Cretaceous and early Cenozoic ancestors of Testudines, Crocodilia, and Lepidosauria distributed in the Holarctic. We implemented 23 fossil constraints for turtles (Table  S8), 30 fossil constraints for lepidosaurs (Table S9), and eight for crocodiles (Table S10). Under the SFC model, turtles were found to have originated in the Northern Hemisphere (under the HFC model they were spread over both regions), whence lineages migrated toward the equator and southern regions (Fig. S16). Most dispersal therefore occurred "into the equator" (Fig. S21; Table S11). We also detected a larger number of geographic extinctions when fossil ranges were considered, predominantly for turtle lineages in the Holarctic (53 and 11 lineages disappeared from this region under the HFC and SFC models, respectively) and in southern temperate regions (nine in the HFC model; Fig.  S21; Table S12). The same pattern of Holarctic extinction was observed when the number of extinction/dispersal events was divided by the number of lineages currently distributed in each region (Fig. 7). The uncertainty associated to this estimation does not affect the overall result (Table S13 and Appendix 1). The most supported biogeographic scenarios in both SFC and HFC analyses also suggest that lepidosaur ancestors were widespread (Figs. S17-S18; uncertainty presented in Tables S14 and Appendix 1). During the greenhouse period, dispersal "into the equator" occurred at the same rate (or at a higher rate in the HFC model) than dispersal "out of the equator," and dispersal "out of the equator" prevailed thereafter ( Fig. S21; Table S11). Estimates of range extinction rates were high in this group under the unconstrained model, with 30 lineages extirpated from the Holarctic, two from southern temperate regions and 152 from the equator (Fig. S21; Table S12). Under fossil-informed models, the number of Holarctic extinctions increased (109 and 66 lineages in the HFC and SFC models, respectively), whereas the number of lineages extirpated from the equator was similar (144 and 109 in the HFC and SFC models, respectively; Fig. S21). When the number of events was controlled for the extant number of lineages distributed in each region, the number of Holarctic extinctions and dispersals "into the equator" increased dramatically, exceeding equatorial dispersal/extinctions (Fig. 7).
For crocodiles, analyses including fossil ranges showed that all deep nodes were distributed in the Holarctic (Figs. S19-20), and range extinctions were detected: four lineages disappeared from the Holarctic, three from southern temperate regions, and two from the equator (HFC model; Fig. S21; Tables S12). Only two lineages disappeared from the Holarctic in the SFC model. The same trends were observed after controlling the number of events for the current number of lineages in each region (Fig. 7). The uncertainty associated to this estimation does not affect the overall result (Table S15 and Appendix 1). A summary of all the results is presented in Table S17.

SHAPING THE LATITUDINAL DIVERSITY GRADIENT
Fossil investigations have shown that, at certain times during the Phanerozoic, the LDG has flattened, weakened, or developed a paleotemperate peak, with diversity at high latitudes being greater than currently for many groups (Mannion et al. 2014;Marcot et al. 2016). This observation has multiple consequences for the study of the LDG. The evolutionary mechanisms required to explain the formation of the current steep LDG are radically different whether we consider or not that high diversity levels existed previously in the Northern Hemisphere: if the pattern transitioned from flatten to steep, as suggested by fossils, then one hypothesis can argue for diversity loss in the Northern Hemisphere to explain the current LDG, via extinction or range contractions (Hawkins et al. 2006). If diversity was never elevated at high latitudes, then an alternative hypothesis can argue for slow accumulation of species in the Northern Hemisphere, due to limited dispersal to the Holarctic (Latham and Ricklefs 1993;Wiens and Donoghue 2004), high Holarctic turnover (Weir and Schluter 2007;Pyron 2014; Pulido-Santacruz and Weir 2016), or high rates of equatorial diversification (Ricklefs 2006;Wiens 2007;Jansson et al. 2013;Pyron and Wiens 2013;Rolland et al. 2014). Hypotheses related to "slow Holarctic diversity accumulation," however, cannot alone account for the formation of a flattened LDG, or for the transition from higher to lower diversity in the Holarctic.
Furthermore, although the processes shaping biodiversity vary over time and space, this has been largely overlooked in the context of the LDG, which has been generally explained in terms of time-constant uniform processes. That is, the parametrization of previous evolutionary models requires only one value per parameter (speciation/extinction/dispersal) and region (Holarctic and Equator) to explain the persistence of a steep LDG in the deep  time. Conversely, in our time-variable framework, these parameters adopt different values per region and through time: one value in each region during the greenhouse period, and different values during icehouses (Table 1). Doing this, our models identify gains and losses of tropical diversity at high latitudes, with prevailing speciation, extinction, and dispersal dynamics changing between warm and cold time intervals. Our time-variable fossil-based analyses support "Holarctic diversity loss" scenarios to explain the LDG of turtles and crocodiles. Diversification rates estimated in the Holarctic and equatorial regions were similar during the equable greenhouse period of the Cretaceous-early Cenozoic for all groups studied here (overlapping credibility intervals; Figs. 3 and S5-S10), consistent with the idea of the existence of a flattened LDG during this phase (Mannion et al. 2014;Marcot et al. 2016). We hypothesize that the expansion of tropical conditions to higher latitudes during greenhouse periods might have induced species diversification in the new paratropical areas (De Celis et al. 2019) and facilitated movements within the broad "tropical belt," such that tropical equatorial clades were able to disperse "out of the equator" (Jablonski et al. 2006(Jablonski et al. , 2013 (Fig. 8). By contrast, the contraction of the tropical biome following climate cooling provoked periods of declining diversity at high latitudes (Fig. 3), where climate change was more intensively felt, and mediated dispersal "into the equator" (Figs. 7 and 8). We found that diversification rates of turtles and crocodiles decreased in all regions during the transition to colder climates (Fig. 3)-they decreased since the Cretaceous in the analyses with time intervals defined by the main geological periods (Figs. S11-S13). The slowing of diversification was much stronger in the Holarctic than at the equator, with extinction exceeding speciation in this region (i.e., Holarctic diversity loss). In addition, using phylogenetic-based biogeographic models informed by fossils, we inferred that all groups had a widespread ancestral distribution that subsequently contracted toward the equator. This result is in agreement with previous fossil investigations for turtles (Nicholson et al. 2015;   Joyce et al. 2016) and crocodiles (Markwick 1998;Mannion et al. 2015;De Celis et al. 2019). Range contraction in our study started in the Cretaceous, intensifying during the late Paleogene cooling. They resulted from range extirpations at higher latitudes combined with "into the equator" dispersals (Condamine et al. 2012) (Figs. 6 and 7). Hence, our results suggest that climate change has likely driven the development of an "asymmetric gradient of extinction and dispersal" (AGED) within the tropical biome, and could have mediated the formation of a steep LDG (Fig. 8).
The AGED hypothesis reconciles previous contending ideas on the origin of the LDG by placing them in a temporal scenario (Table 1; Fig. 8). For instance, there is controversial support around the tropics being "cradle" or "museum of diversity" (Stebbins 1974), and dispersal prevailing "out of" (Jablonski et al. 2006(Jablonski et al. , 2013 or "into the tropics" (Condamine et al. 2012;Pyron 2014;Rolland et al. 2015). The interpretation of our results alternatively invokes the "museum of diversity" regarding the equatorial tropics as refuge during icehouse transitions, but also the "cradle of diversity" during greenhouse periods. Similarly, our hypothesis invokes "out of the equator" dispersals during greenhouse transitions and "into the equator" dispersals during icehouse transitions.
Support for the AGED hypothesis and for a "Holarctic diversity loss" scenario is mixed for lepidosaurs. On the one hand, we found similar diversification rates in the Holarctic and equator during the greenhouse period, and widespread ancestral distributions (Fig. 3). We also found a higher proportion of lepidosaur species that actually lost their ancestral Holarctic distribution and dispersed "into the equator" (Pyron 2014) than the other way around (Fig. 7), in agreement with the idea of an ancient flattened pattern shaped by extinction. On the other hand, we detected that during climate cooling, diversity losses of lepidosaurs occurred only in the equator (Fig. 3). However, equatorial estimates remain uncertain: diversity dynamics for species distributed at the equator have broad credibility intervals probably due to the poverty of the equatorial dataset in terms of the number of fossil lineages and the small number of records per lineage (Table S16). In the Holarctic, turnover rates were high during the transitional period to cold, indicating that species did disappear from high latitudes, but that a lepidosaur community got replaced by another. This result suggests the number of lepidosaur species may always have been unbalanced between regions with higher diversity in the equator. The high Holarctic turnover likely contributed to the maintenance of this pattern, together with the inferred temporal increases in diversification at the equator (Fig. 3), as previously hypothesized (Pyron 2014).

CONSERVATISM FRAMEWORK TO EXPLAIN THE LDG
Accumulating fossil, ecological and molecular evidence demonstrates that global climate changes over geological timescales could generate large-scale patterns of biodiversity (Mannion et al. 2014;Fenton et al. 2016;Saupe et al. 2019a). In the last decade, phylogenetic niche conservatism (PNC), that is, the tendency of species to retain their ancestral niches over time, emerged as a general principle to explain the effects of climate over diversity (Peterson et al. 1999). Wiens and Donoghue (2004) proposed PNC as a major explanation behind the LDG. The tropical niche conservatism (TNC) hypothesis posits that the difficulty of many tropical lineages to invade or persist in temperate environments determined the distribution of global diversity (Wiens and Donoghue 2004;Donoghue 2008). However, they were less specific about the mechanisms by which PNC shaped diversity. They considered time, limited dispersal, and also the contraction of the tropical belt (Wiens and Donoghue 2004;Donoghue 2008). They argued that tropical regions had a greater geographical extent in the past to explain why most taxa have tropical adaptations, but did not consider that diversity could have once reached equivalent levels across latitudes. This probably explains why extinction was not part of their original predictions. For example, Wiens et al. (2010) published a review on the topic where only limited dispersal explains the LDG. Moreover, diversity loss (and range contractions), if ever seen as an original prediction of the TNC, has not transcended in the literature as part of the TNC formulations. This is manifest in subsequent studies, which interpreted the LDG in terms of limited dispersal and PNC (Cardillo et al. 2005;Condamine et al. 2012;Kerkhoff et al. 2014;Rolland et al. 2014;Siqueira et al. 2016). For example, Kerkhoff et al. (2014) finds that high current tropical diversity results from a combination of (i) differential net diversification rates of tropical lineages due to larger cumulative area of tropical environments, (ii) greater time for diversification in tropical environments, and (iii) limited dispersal of tropical lineages into the temperate environments. We therefore argue that the AGED model represents an extension of the TNC hypothesis, by formalizing the mechanisms by which PNC shaped diversity (extinction and range contraction during icehouses, speciation, and range expansion during greenhouses) in a time-variable scenario and a varying LDG shape (Fig. 8).
Behind the AGED model prevails the PNC idea: when the tropical biome retreated toward the equator, most of the tropical-adapted taxa at high latitudes were unable to adapt and either went extinct or got their distribution restricted. The ancestors of turtles, lepidosaurs, and crocodiles were adapted to tropical conditions during the Late Cretaceous (Markwick 1998;Waterson et al. 2016;Pie et al. 2017). Our fossil-based diversification results when analyzing differences between taxa adapted to different climates indicate that extinction events were not random (Eiserhardt et al. 2015;Reddin et al. 2019), but instead preferentially affected taxa living in tropical-like climates at high latitudes (Fig. 4). Similar climatic-driven extirpation scenarios have been proposed for other vertebrates (Hawkins et al. 2006;Saupe et al. 2019b). This suggests that climate change could have imposed an asymmetric gradient of extinction and dispersal within the tropical biome due to PNC.
Lineages possessing or having evolved the appropriate adaptations to cope with cold climates diversified in the new temperate areas (Kindlmann et al. 2007;Meseguer et al. 2018). After the transition to temperate climates in the late Eocene, we found that diversification rates of turtles, crocodiles, and lepidosaurs living in temperate climatic conditions were significantly higher than those of tropical-adapted taxa living in Holarctic and equatorial regions (Fig. 4). New Neogene temperate habitats likely constituted an opportunity for diversification due to increased geographic ranges and ecological niches (Wiens 2007), eventually driving an inverse LDG for some groups (Kindlmann et al. 2007;Leslie et al. 2012). Several radiations following the appearance of temperate biomes have been identified in other groups, such as plants (Meseguer et al. 2018), mammals (Ge et al. 2013), or insects (Condamine et al. 2018). After this period, we estimated similar diversification rates between tropical and temperate lineages because speciation decreased dramatically in the temperate lineages of our focal groups (Fig. 4), possibly due to the effect of the Pleistocene glaciations. Thus, our study does not support lower rates of diversification under cold environments, in agreement with previous studies (Weir and Schluter 2007;Pyron 2014).
These results overall suggest that differences in species richness between geographic regions (i.e., the Holarctic vs. the equator) may be explained by asymmetric diversification and dispersal across regions and time. Differences in species richness between ecological types (i.e., tropical-vs. temperateadapted taxa) may be explained by the longer time available for tropical-adapted clades to diversify in tropical areas (Stephens and Wiens 2003). Turning up that the factors driving diversity differ between geographic regions and ecological types.
Here, we encoded tropical/temperate preferences by considering macro-conditions in the Holarctic to be paratropical until the end of the Eocene, as generally recognized (Ziegler et al. 2003;Sluijs et al. 2006;Morley 2011;and references therein), and temperate afterward. However, we acknowledge that the post-Eocene Holarctic could be seen as tropical like during warming events. We also assumed that taxa inhabiting the warm Holarctic were adapted to tropical-like conditions (i.e., a high global temperature, indicating probable adaptation to tropical cli-mates). For turtles, crocodiles, and lepidosaurs, this assumption is supported by Cenozoic climatic niche modeling (Waterson et al. 2016;Pie et al. 2017), stable isotope analyses, and other climate proxies (Markwick 1998;Tütken and Absolon 2015). These assumptions are, of course, oversimplifications that may introduce biases in the analyses, but we consider that general patterns may nevertheless emerge from such analyses (Romdal et al. 2013).

EVIDENCE
Our results demonstrate that the inclusion of fossils in macroevolutionary studies improves detecting ancient high-latitude extinctions and range extirpations, otherwise hardly detectable with molecular phylogenies only. The results exclusively based on extant species (time-constant and time-variable BiSSE models, and biogeographic analyses; Figs. 6, 7, S1-S4, and S14-S15) differed from the analyses including fossils in our study; they suggest higher levels of Holarctic diversification for turtles during Cenozoic cooling, together with an equatorial origin and recent invasion of high-latitude regions, resulting in less time for lineages to diversify in the Holarctic (Figs. 6, 7, and S1-S4). This is in agreement with the TNC hypothesis (Wiens and Donoghue 2004) and recent investigations (Pereira et al. 2017;Rodrigues et al. 2017). For crocodiles, they support the diversification hypothesis, with higher origination rates close to the equator and no effect on dispersal (Figs. 7 and S15). For lepidosaurs, they support the high temperate turnover and into the tropics hypotheses (Figs. 7, S1-S4, and S14), as suggested before (Pyron 2014). This conflict between extant and fossil evidence may extend beyond our study, pervading the LDG literature.
The observed incongruences between paleontological and neontological analyses could be attributed to the geographic and/or preservation biases of the fossil record (Albino and Brizuela 2014). We argue, however, that sampling artifacts on fossils have not biased diversification and biogeographic results for Testudines and Crocodilia. Preservation rates estimated for these lineages are similar for Holarctic and equatorial regions (overlapping credibility intervals; Table 2). PyRate has also been shown to correctly estimate the dynamics of speciation and extinction rates under low levels of preservation or severely incomplete taxon sampling (Silvestro et al. 2014). Nevertheless, because diversification dynamics was inferred independently for Holarctic and equatorial fossils, inferences of the Holarctic are robust and based on a large number of fossils (Table S16), indicating diversity loss during the climate transition.
In the case of lepidosaurs, preservation rates are higher in the Holarctic than at the equator, and thus we have a relative oversampling of Holarctic fossils (Table 2). Although this might have a lower impact on diversification estimates (as explained above), we cannot discard that higher preservation rates at higher latitudes potentially biased the biogeographic results for lepidosaurs toward nonequatorial origins, and artificially increased the estimates of "into the equator" dispersals. Nonetheless, the biogeographic reconstruction based on extant data only was not able to recover Holarctic distributions for lepidosaur ancestors, despite extensive fossil evidence suggesting a widespread (Holarctic and equatorial) distribution for these lineages (Table  S16; Fig. 7). Meanwhile, the reconstruction including fossils recovers widespread ancestors. The number of Holarctic and equatorial fossil constraints for lepidosaurs was relatively low given the size of the tree, but these constraints significantly increased the absolute number of Holarctic range extinctions (from 30 to 109) and "into the equator" dispersals (from 40 to 124) relative to estimates without such constraints. Meanwhile, the inclusion of fossil data did not alter the number of events estimated for equatorial taxa (Tables S11-S12). This finding suggests that there is an oversampling of equatorial species in the present (possibly due to extinction of Holarctic species in the past), which may lead to spurious increases of the probability of estimating dispersions "out of the equator," compared to "into the equator" if fossils are not considered. A better understanding of lepidosaur fossil taxonomy might facilitate the assignment of fossils on the tree, and the detection of additional high-latitude range extinctions not detected here.
Incongruences between paleontological and neontological analyses have also been attributed to the difficulty to estimate extinction rates from phylogenetic data. Although demonstrated mathematically (Nee et al. 1994;Sanmartín and Meseguer 2016), estimations of extinction in phylogenetic analyses are often inaccurate and highly underestimated (Rabosky 2010). This has been attributed to sampling biases (Pybus and Harvey 2000), lack of statistical power (Davis et al. 2013), violations of model assumptions (Rabosky 2010;Morlon et al. 2011), or substantial rate heterogeneity (Morlon et al. 2011;Sanmartín and Meseguer 2016). We cannot exclude that these artifacts affected our phylogenetic-based results. In the last years, however, studies proposed alternative explanations to biases in the data/methods. Conceptual differences on how ages, speciation, and extinction rates are estimated from paleontological data and molecular phylogenies could explain some of the incongruences (Huang et al. 2015;Silvestro et al. 2018). In our study, fossil and phylogenetic data for Crocodilia and Testudines fit the expectations of the BDC model, suggesting that discrepancies between phylogenetic and paleontological estimates are probably attributed to intrinsic differences in the fossil and phylogenetic properties due to bifurcating and anagenetic speciation events (Table S7; Fig. 5). These speciation modes are not detectable from phylogenetic data, which could alter evolutionary rates and contribute to an apparent incongruence between fossil and phylogenetic estimates (Silvestro et al. 2018). The BDC model was rejected for Lepidosauria, suggesting that differences between fossils and phylogenies may not be entirely explained by different speciation modes and are rather attributable to potential biases in the data/methods.
In previous LDG studies, difficulties to recover diversityloss scenarios could be attributed as well to the use of SSE models assuming that diversification parameters remain constant over time (Pyron and Wiens 2013;Pyron 2014;Rolland et al. 2015;Pulido-Santacruz and Weir 2016;Siqueira et al. 2016). Time-constant models are inadequate for testing complex scenarios where the processes underlying the LDG varied across latitudes and time. Moreover, the power of time-constant models for detecting negative diversification rates is questionable, because inferring negative diversification for the entire history of lineages conflicts with the fact that these groups are still extant. Testing the hypothesis of diversity declines thus requires the implementation of time-variable models (Morlon et al. 2011;Freyman and Hoehna 2018). When applied to the study of diversity patterns, time-variable SSE models have revealed extinction signatures in tropical clades (Spriggs et al. 2015). The incorporation of time-shifts into BiSSE improves but not completely reconciles paleontological and phylogenetic evidence in our study. These artifacts highlight the importance of combining fossils and phylogenies in macroevolutionary inferences (Fritz et al. 2013). The fossil record remains incomplete, but it provides the only direct evidence of the past diversity, and thus is fundamental in the study of extinction scenarios.

Conclusion
Speciation, extinction, and dispersal shape the LDG, but the contribution of these processes remains a debated topic in evolutionary ecology. Our results indicate these processes operated at different rates over time and space as climate changed. The current LDG of turtles and crocodilians can be explained by ancient high-latitude tropical diversity loss and range contractions as a consequence of the retraction of the tropical biome due to climate cooling. During greenhouse periods, equivalent diversification rates across latitudes prevailed, explaining the formation of a flattened LDG. Changes in global diversification and dispersal dynamics imposed by large-scale climatic transitions could represent a mechanism that shapes the LDG. This "asymmetric gradient of extinction and dispersal" hypothesis might account for the LDG of tropical-adapted groups that were once diverse at high latitudes, but might not be fully applicable to all organisms currently displaying a LDG, as shown here for lepidosaurs.

AUTHOR CONTRIBUTIONS
ASM and FLC designed the study and analyzed the data. ASM wrote the paper with contributions of FLC.
Associate Editor: T. Ezard Handling Editor: D. W. Hall

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.  Table S7. Speciation and extinction rates estimated from extant phylogenies and stratigraphic ranges of Crocodyles, Testudines and Lepidosauria under the equal, compatible and incompatible rates models assuming constant diversification rates over time Table S8. Fossil constraints used for the biogeographic analyses of turtles. Table S9. Fossil constraints used for the biogeographic analyses of Squamata. Table S10. Fossil constraints used for the biogeographic analyses of Crocodiles. Table S11. Estimated number of dispersals events through time from the Holarctic into the equator and out of the equatorial zone for Testudines, Squamata and Crocodiles under the unconstrained (Unc.), based on present evidence only, and the fossil-based hard (HFC) and soft fossil constraint (SFC) biogeographic models. Table S12. Estimated number of range-extinction events through time for Testudines, Squamata and Crocodilia under the unconstrained (Unc.), based on present evidence only, and the fossil-based hard (HFC) and soft fossil constraint (SFC) biogeographic models. Table S13. Uncertainty of the biogeographic reconstruction of Testudines inferred from Lagrange (DEC) under three different models (unconstrained, Unc.; soft fossil constraint, SFC; and hard fossil constraint model, HFC; see methods) on the 10 most basal nodes of the tree. Table S14. Uncertainty of the biogeographic reconstruction of Squamata inferred from Lagrange (DEC) under three different models (unconstrained, Unc.; soft fossil constraint, SFC; and hard fossil constraint model, HFC; see methods) on the 10 most basal nodes of the tree. Table S15. Uncertainty of the biogeographic reconstruction of Crocodilia inferred from Lagrange (DEC) under three different models (unconstrained, Unc.; soft fossil constraint, SFC; and hard fossil constraint model, HFC; see methods) on the 10 most basal nodes of the tree. Table S16. Number of fossil occurrences (occ.) per genera in Equatorial, Tropical, Temperate and Holarctic datasets for each group (occurrences for the Southern Hemisphere temperate regions not shown). Table S17. Results of the diversification and biogeographic analyses performed in this study for the transition from greenhouse to coldhouse climates. Figure S1a. Latitudinal diversification pattern of turtles inferred with BiSSE showing the estimates for speciation (a), extinction (b), transition (c) and net diversification (d) rates for Holarctic and Equatorial species. Figure S1b. Latitudinal diversification pattern of squamates inferred with BiSSE showing the estimates for speciation (a), extinction (b), transition (c) and net diversification (d) rates for Holarctic and Equatorial species.       Figure S5. Global pattern of turtle diversification based on the fossil record and estimated using time bins as defined by climatic epochs. Figure S6. Global pattern of Testudines diversification based on the fossil record and estimated using time bins as defined by climatic periods (51-33.7 Ma). Figure S7. Global pattern of squamates diversification based on the fossil record and estimated using time bins as defined by climatic epochs. Figure S8. Global pattern of Lepidosauria diversification based on the fossil record and estimated using time bins as defined by climatic periods (51-33.7 Ma). Figure S9. Global pattern of crocodiles diversification based on the fossil record and estimated using time bins as defined by climatic epochs. Figure S10. Global pattern of Crocodilia diversification based on the fossil record and estimated using time bins as defined by climatic periods (51-33.7 Ma). Figure S11. Global pattern of turtle diversification based on the fossil record and estimated using time bins as defined by geological periods. Figure S12. Global pattern of lepidosaurs diversification based on the fossil record and estimated using time bins as defined by geological epochs. Figure S13. Global pattern of crocodiles diversification based on the fossil record and estimated using time bins as defined by geological epochs. Figure S15. Biogeographic reconstruction of Crocodiles on DEC under the Unconstrained model. Figure S16. Biogeographical reconstruction of Testudines inferred with DEC including soft fossil constraints (SFC, see main text). Figure S19. Biogeographic reconstruction of Crocodiles showing the effects of the incorporation of fossil information into biogeographic inference under a soft fossil constraint (SFC, see main text) model. Figure S20. Biogeographic reconstruction of Crocodiles showing the effects of the incorporation of fossil information into biogeographic inference under a hard fossil constraint (HFC, see main text) model. Figure S21. Estimated number of range-extinction and dispersal events across regions and through time for Testudines, Lepidosauria and Crocodiles, relative to the current number of lineages in each group.