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

Andrea S. Meseguer, Fabien L. Condamine

► To cite this version:

Andrea S. Meseguer, Fabien L. Condamine. Ancient tropical extinctions at high latitudes contributed to the latitudinal diversity gradient*. Evolution - International Journal of Organic Evolution, 2020, 74 (9), pp.1966-1987. 10.1111/evo.13967. hal-03049599v2

HAL Id: hal-03049599

https://hal.science/hal-03049599v2

Submitted on 1 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

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

Andrea S. Meseguer<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" /> ID and Fabien L. Condamine<fnref n="2" /> ID

Received April 15, 2019

Accepted March 21, 2020

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.

**KEY WORDS:** Biodiversity, climate change, extinction, fossils, Holarctic, tropics.

The current increase in species richness from the poles toward the equator, known as the latitudinal diversity gradient (LDG), is one of the most conspicuous patterns in ecology and evolution. This pattern has been described for microbes, insects, vertebrates, and plants, and for marine, freshwater, and terrestrial ecosystems (Willig et al. 2003; Hillebrand 2004; Novotny et al. 2006; Kreft and Jetz 2007; Fuhrman et al. 2008; Jenkins et al. 2013).

For decades, it has been thought that the modern steep LDG, with higher diversity concentrated at the equator, persisted

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),

<footnote n="1">¹INRA, UMR 1062 Centre de Biologie pour la Gestion des Populations (INRA | IRD | CIRAD | Montpellier SupAgro), Montferrier-sur-Lez, France</footnote>

<footnote n="2">²CNRS, UMR 5554 Institut des Sciences de l'Evolution de Montpellier (Université de Montpellier | CNRS | IRD | EPHE), Montpellier, France</footnote>

<footnote n="3">³Real Jardín Botánico de Madrid (RJB-CSIC), Madrid, Spain</footnote>

<footnote n="4">⁴E-mail: asanchezmeseguer@gmail.com</footnote>

<footnote>*This article corresponds to Lacroix R. (2020). Digest: The contribution of historical climate events in shaping the modern latitudinal diversity gradient of ancient reptiles. Evolution. https://doi.org/10.1111/evo.14069.</footnote>
<figure><img src="image_4.png" /></figure>

**Figure 1.** Changes in global temperatures and extension of the tropical belt during the Cenozoic, in relation with the shape of the LDG. Early Cenozoic global temperatures were higher than today and paratropical conditions extended over high latitudes. From the early Eocene climatic optimum (EECO; about 53–51 Ma), a global cooling trend intensified and culminated with the Pleistocene glaciations. Warm-equable regimes got then restricted to the equator. The LDG evolved following these global changes; during greenhouse periods diversity was similar across latitudes, such that the LDG flattened, whereas in cold periods diversity peaked at the equator (a steep LDG) (Mannion et al. 2014). The question mark denotes the focus of this study, which is to unveil the processes that mediated the transition between a flat and steep LDG. The relative temperature curve of the Cenozoic is adapted from (Zachos et al. 2008). Maps represent the extension of the tropical belt and Earth tectonic changes as derived from (Ziegler et al. 2003; Morley 2011). P = Pleistocene; Pli = Pliocene.

insects (Archibald et al. 2010, Archibald et al. 2013; Labandeira and Currano 2013), brachiopods (Krug and Patzkowsky 2007; Powell 2007; Powell et al. 2012), bivalves (Crame 2000, 2020), coral reefs (Leprieur et al. 2016), foraminifers (Fenton et al. 2016), crocodiles (Mannion et al. 2015), turtles (Nicholson et al. 2015, Nicholson et al. 2016), 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 Polarctic 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 ($\lambda$) and extinction ($\mu$) rates ($\lambda_h \approx \mu_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 Schlüter 2007), butterflies (Condamine et al. 2012), plants (Leslie et al. 2012; Kerkhoff et al. 2014), fishes (Siqueira et al. 2016), mammals (Weir and Schlüter 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 Schlüter 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 ($\lambda \approx \mu$) (Cardillo et al. 2005; Weir and Schlüter 2007; Condamine et al. 2012; Leslie et al. 2012; Pyron and Wiens 2013; Rolland et al.

2014, 2015; Rabosky 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 ($\mu > \lambda$), resulting in negative net diversification rates ($r = \lambda - \mu$; $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. (2010, 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}$; $\mu_h \le \lambda_h$) vs. “diversity loss” ($d_{eh} < d_{he}$; $\mu_h > \lambda_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 warm- and 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, Nicholson et al. 2016). All these lineages
<table><caption>**Table 1.** Predictions of the most common LDG hypotheses, including the *Asymmetric gradient of extinction and dispersal* (AGED) proposed here.</caption><thead><tr><th>1 Mechanism</th><th>Time (t)</th><th>Differences in r</th><th colspan="2">Differences in d</th><th>Both</th></tr></thead><tbody><tr><td>2. Rate</td><td>Time-constant models</td><td></td><td colspan="2"></td><td>Time-variable models</td></tr><tr><td>3. Source</td><td>Lower Holarctic species accumulation</td><td></td><td colspan="2"></td><td>Holarctic diversity loss</td></tr><tr><th>Hypothesis</th><th>Time for speciation</th><th>Cradle of diversity</th><th>Museum of diversity</th><th>Holarctic turnover</th><th>Into the tropics</th><th>Tropical niche conservatism</th><th>Out of the tropics</th><th>Asymmetric gradient of extinction and dispersal</th></tr><tr><th>Predictions</th><td>t<sub>top</sub> > t<sub>temp</sub><br>(r<sub>top</sub> = r<sub>temp</sub>)<br>(d<sub>temp</sub> = d<sub>temp</sub>)</td><td>r<sub>e</sub> > r<sub>h</sub><br>(λ<sub>e</sub> > λ)<br>d<sub>eh</sub> = d<sub>he</sub></td><td>r<sub>e</sub> > r<sub>h</sub><br>(μ<sub>e</sub> < μ<sub>h</sub>)<br>d<sub>eh</sub> = d<sub>he</sub></td><td>r<sub>h</sub> = 0<br>(μ<sub>h</sub> = λ<sub>h</sub>)<br>d<sub>eh</sub> = d<sub>he</sub></td><td>r<sub>e</sub> = r<sub>h</sub><br>d<sub>eh</sub> < d<sub>he</sub></td><td>r<sub>e</sub> = r<sub>h</sub><br>d<sub>eh</sub> > d<sub>he</sub></td><td>r<sub>e</sub> > r<sub>h</sub><br>(λ<sub>e</sub> > λ<sub>h</sub>)<br>(μ<sub>e</sub> < μ<sub>h</sub>)<br>d<sub>eh</sub> > d<sub>he</sub></td><td>Ice-greenhouse transition<br>Green-icehouse transition<br>r<sub>e</sub> = r<sub>h</sub><br>(μ<sub>h</sub> < λ<sub>h</sub>; r<sub>h</sub> > 0)<br>d<sub>eh</sub> > d<sub>he</sub><br>r<sub>e</sub> > r<sub>h</sub><br>(μ<sub>h</sub> > λ<sub>h</sub>; r<sub>h</sub> < 0)<br>d<sub>eh</sub> < d<sub>he</sub></td></tr></tbody></table>

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 (λ<sub>t</sub> > λ<sub>e</sub>) under the "*cradle of diversity*," and/or result from lower rates of extinction (μ<sub>t</sub> < μ<sub>e</sub>) 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 Schlüter 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.

Abbreviations: λ = speciation; μ = extinction; d = dispersal; t = time; e = equatorial; h = Holarctic; trop = Tropical biome; temp = Temperate biome.
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.

## Methods

### 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 assignment 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.

### ESTIMATION OF ORIGINATION AND EXTINCTION 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," $\lambda_H$) and the other with other equatorial and subtropical regions ("equator" or "E," $\lambda_E$), two extinction rates associated with the Holarctic ($\mu_H$) and the equator ($\mu_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%.

### ESTIMATION OF ORIGINATION AND EXTINCTION 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
4798 occurrences for 638 genera (120 extant and 518 extinct; 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 (temperature 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 ($\lambda^*$ and $\mu^*$) are the same to those estimated with phylogenetic data ($\lambda$ and $\mu$), that is, $\lambda^* = \lambda$, $\mu^* = \mu$; (ii) a Compatible model where parameters differ, but the differences could be explained by differences in speciation mode. In this model, $\lambda^*$, $\lambda$, $\mu$, and $\mu^*$ are constrained such that $\lambda^* - \lambda = \mu^* - \mu$ (i.e., equal net diversification rates) and $\lambda^* \ge \lambda$; (iii) an Incompatible model, where parameters $\lambda$, $\lambda^*$, $\mu$, and $\mu^*$ are allowed to take any value and thus differences in $\lambda$ and $\lambda^*$, as well as $\mu$ and $\mu^*$, cannot be explained by differences in speciation mode.

We estimated $\lambda$, $\mu$, $\lambda^*$, and $\mu^*$ simultaneously for the Testudines, Crocodilia, and Lepidosauria phylogenies pruned to the genus level, and the corresponding fossil data using maximum-likelihood 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 $\lambda^*$, $\mu^*$, $\lambda$, and $\mu$ obtained under the incompatible rates model to verify the conditions predicted by the compatible BDC model (i.e., $\lambda^* - \lambda = \mu^* - \mu$ and $\lambda^* \ge \lambda$) and assess the support for each model as suggested in previous studies (Silvestro et al. 2018).

## INFERRING ANCESTRAL GEOGRAPHIC DISTRIBUTION WITH PHYLOGENIES AND FOSSILS

We performed biogeographic analyses with the parametric likelihood method DEC (Ree and Smith 2008) 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.reptile-database.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).

## Results

### PHYLOGENY-BASED DIVERSIFICATION ANALYSES: ARE DIVERSIFICATION RATES HIGHER AT THE 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 (51/66 Ma and 23/34 Ma; Figs. S3 and S4).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Global pattern of turtle, lepidosaur, and crocodile diversification through time based on the fossil record. Origination (blue) and extinction (red) rates were estimated using time bins as defined by epochs of the geological timescale (on the top, main climatic periods are shown as follows: Greenhouse, Tran. = climatic transition, and I. = icehouse). Solid lines indicate mean posterior rates, whereas the shaded areas show 95% credibility intervals. Net diversification rates (green) are the difference between origination and extinction. The vertical lines indicate the boundaries between geological periods. Tr = Triassic; J = Jurassic; K = Cretaceous; Pg = Paleogene; Ng = Neogene.</figcaption></figure>

## 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
Holarctic vs. equatorial

<figure><img src="image_6.png" /></figure>

**Figure 3.** Global pattern of turtle, lepidosaur, and crocodile diversification between Holarctic and equatorial regions, based on the fossil record. Diversification dynamics are compared between fossils distributed in Holarctic and equatorial regions. Origination (blue) and extinction (red) rates were estimated using time bins as defined by the main climatic intervals since the Mesozoic (on the top, climatic periods are shown as follows: Greenhouse, Tran. = climatic transition, and I. = icehouse). Solid lines indicate mean posterior rates, whereas the shaded areas show 95% credibility intervals. Net diversification rates (green) are the difference between origination and extinction. The vertical lines indicate the boundaries between climatic intervals. Tr = Triassic; J = Jurassic; K = Cretaceous; Pg = Paleogene; Ng = Neogene.

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).

## ESTIMATIONS OF ANCESTRAL ORIGINS: DID GROUPS 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
Temperate vs. tropical

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Global pattern of turtle, lepidosaur, and crocodile diversification across temperate and tropical climates, based on the fossil record. Diversification dynamics is compared between fossils inhabiting under temperate and tropical macroclimates (see text). Origination (blue) and extinction (red) rates were estimated using time bins as defined by the main climatic intervals since the Mesozoic (on the top, climatic periods are shown as follows: Greenhouse, Tran. = climatic transition, and I. = icehouse). Solid lines indicate mean posterior rates, whereas shaded areas show 95% credibility intervals. Net diversification rates (green) are the difference between origination and extinction. Vertical lines show boundaries between climatic intervals. Tr = Triassic; J = Jurassic; K = Cretaceous; Pg = Paleogene; Ng = Neogene.</figcaption></figure>

(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
<table><caption>**Table 2.** Posterior (median) estimates of the preservation rate and heterogeneity parameter.</caption><thead><tr><th>Clade</th><th>Preservation rate</th><th>Heterogeneity parameter</th><th>Clade age (Ma)</th></tr></thead><tbody><tr><td>All crocodiles</td><td>1.2123 (1.0392–1.3431)</td><td>0.7566 (0.5207–0.89)</td><td>144.61 (140.99–152.14)</td></tr><tr><td>Holarctic crocodiles</td><td>1.1763 (1.0745–1.2881)</td><td>0.4439 (0.3653–0.5979)</td><td>145.81 (140.17–160.34)</td></tr><tr><td>Equatorial crocodiles</td><td>0.9774 (0.8264–1.1529)</td><td>0.7474 (0.5425–0.9496)</td><td>128.36 (125.03–129.43)</td></tr><tr><td>All turtles</td><td>0.7953 (0.754–0.8441)</td><td>0.6243 (0.5834–0.6692)</td><td>256.54 (248.8–275.06)</td></tr><tr><td>Holarctic turtles</td><td>0.8601 (0.7139–0.929)</td><td>0.6407 (0.5828–0.6996)</td><td>256.78 (250.68–273.94)</td></tr><tr><td>Equatorial turtles</td><td>0.897 (0.6289–1.1096)</td><td>0.3941 (0.176–0.4903)</td><td>145.2 (140.22–155.92)</td></tr><tr><td>All lepidosaurs</td><td>1.3108 (1.0566–1.5153)</td><td>0.6561 (0.6022–0.7169)</td><td>235.08 (228.69–247.45)</td></tr><tr><td>Holarctic lepidosaurs</td><td>1.5333 (1.116–1.6753)</td><td>0.6498 (0.6065–0.6973)</td><td>235.79 (228.05–249.92)</td></tr><tr><td>Equatorial lepidosaurs</td><td>0.3132 (0.2378–0.4074)</td><td>0.6934 (0.4226–1.037)</td><td>178.46 (173.73–182.08)</td></tr></tbody></table>

Preservation rate is the expected number of occurrences per lineage per million year and the heterogeneity parameter is the shape parameter of the Gamma distribution, as estimated under the PyRate MCMC model. Lower and upper bounds of the 95% credibility intervals are provided in parentheses.

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.

## Discussion

### TIME-VARIABLE EVOLUTIONARY PROCESSES 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 Schlüter 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
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Results from a joint Bayesian analysis of fossil and phylogenetic data for crocodiles, lepidosaurs, and turtles under the skyline BDC model. Posterior samples of speciation (blue) and extinction (red) rates are plotted against one another and calculated over three time intervals to account for rate heterogeneity; posterior samples of the two terms $\lambda^* - \lambda$ and $\mu^* - \mu$ are shown in black. Analyses were run assuming independent rates ($\lambda$ and $\mu$ for phylogeny; $\lambda^*$ and $\mu^*$ for fossils), and their joint posterior samples were used to assess which model (equal rates, compatible rates, or incompatible rates) best fit the data. The best model is indicated by the labels in the plots, and based on whether the posterior samples conform to the properties of the equal rate model ($\lambda = \lambda^*; \mu = \mu^*$), and the compatible rate BDC model: $\lambda < \lambda^*$, ($\lambda^* - \lambda$) = ($\mu^* - \mu$). P-values indicate whether phylogenetic and fossil parameter estimates differ significantly. For crocodiles (also for turtles), there is little phylogenetic information in speciation and extinction rates before to 23 Ma (as shown by the large spread of posterior values). For Lepidosauria, there is little fossil information, especially after 23 Ma.</figcaption></figure>

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 “Holartic diversity loss” scenarios to explain the LDG of turtles and crocodiles. Diversification rates estimated in the Holartic 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, 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 Holartic than at the equator, with extinction exceeding speciation in this region (i.e., Holartic 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, 2016;
<figure><img src="image_9.png" /><figcaption>B Fossil-informed biogeographical model</figcaption></figure>

<figure><img src="image_10.png" /></figure>

**Figure 6.** Biogeographic estimations of Testudines showing the effects of the incorporation of fossil information into biogeographic inference. (A) Results with DEC based on the distribution of extant taxa. (B) Results under the fossil-informed HFC (hard fossil constraint) model. Colored circles at tips and nodes represent current and ancestral ranges, respectively, whereas squares represent range inheritance scenarios. Colors correspond with the discrete areas in the legend. Black circles indicate fossil range constraints included in the analysis, with numbers corresponding with taxa in Table S8. The reconstruction under the soft fossil constraint (SFC; see text) model is presented in Figure S16. Tr = Triassic; J = Jurassic; K = Cretaceous; Pg = Paleogene; Ng = Neogene.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Estimated number of range extinction and dispersal events through time. Analyses were performed for turtles, lepidosaurs, and crocodiles under the unconstrained model (Unc.), based on present evidence only, and the fossil-based hard (HFC) and soft fossil constraint (SFC) biogeographic models. (A) Inferred number of range extinction events through time and across regions relative to the number of lineages currently distributed in each region. The black line represents the global mean temperature curve as modified from (Zachos et al. 2008). (B) Inferred number of dispersal events from the Holarctic into the equator (IntoEq) and out of the equatorial zone (OutEq), relative to the current number of lineages distributed in the Holarctic and equatorial zones, respectively. Tr = Triassic; J = Jurassic; K = Cretaceous; Pg = Paleogene; Ng = Neogene; Tran. = climatic transition; Ice. = icehouse.</figcaption></figure>

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, 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
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Prevalent evolutionary processes behind the latitudinal diversity gradient under the AGED model. It shows the change in evolutionary dynamics between Holarctic and equatorial regions through different climatic intervals: the greenhouse, icehouse, and transitions. For each climatic interval, inset plots represent the hypothetical distribution of species richness across latitudes (LDG shape).</figcaption></figure>

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 lepidosauurs 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).

## TOWARD AN INTEGRATIVE PHYLOGENETIC NICHE 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 lepidosauurs 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 Schlüter 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. temperate-adapted 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 lepidosauurs, 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).

## RECONCILING FOSSIL AND PHYLOGENETIC 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 lepidosauurs, 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 lepidosauurs, 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 lepidosauurs 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 lepidosauurs 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 diversity-loss 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 lepidosauurs.

## AUTHOR CONTRIBUTIONS

ASM and FLC designed the study and analyzed the data. ASM wrote the paper with contributions of FLC.
## ACKNOWLEDGMENTS

This preprint has been reviewed and recommended by Peer Community In Evolutionary Biology (https://doi.org/10.24072/pci.evolbiol.100068). The authors are very grateful to Drs. T. Ezard, J. Hortal, J. Arroyo, A. Mooers, and J. Calatayud, and to the various anonymous reviewers for comments and suggestions that greatly improved the study. Previous versions of the manuscript benefited from the comments of Drs. G. Mittlebach, E. Jousselin, and J. Rolland. Financial support was provided by a Marie-Curie FP7-COFUND (AgreenSkills fellowship-26719) grant to ASM and a Marie Curie FP7-IOF (project 627684 BIOMME) grant to FLC. This work benefited from an "Investissements d'Avenir" grant managed by the "Agence Nationale de la Recherche" (CEBA, ref. ANR-10-LABX-25-01).

## DATA ARCHIVING

The datasets supporting the results, the commands used in the analyses, and the Appendix 1 are stored in Dryad (https://doi.org/10.5061/dryad.zs7h44j5m).

## CONFLICT OF INTEREST

The authors declare no conflict of interest.

## LITERATURE CITED

Albino, A. M., and S. Brizuela. 2014. An overview of the South American fossil squamates. Anat. Rec. 297:349–368.

Alroy, J., M. Aberhan, D. J. Bottjer, M. Foote, F. T. Fürsich, P. J. Harries, A. J. W. Hendy, S. M. Holland, L. C. Ivany, W. Kiessling, et al. 2008. Phanerozoic trends in the global diversity of marine invertebrates. Science 321:97–100.

Angielczyk, K. D., R. W. Burroughs, and C. R. Feldman. 2015. Do turtles follow the rules? Latitudinal gradients in species richness, body size, and geographic range area of the world's turtles. J. Exp. Zool. Part B Mol. Dev. Evol. 324:270–294.

Archibald, S. B., W. H. Bossert, D. R. Greenwood, and B. D. Farrell. 2010. Seasonality, the latitudinal gradient of diversity, and Eocene insects. Paleobiology 36:374–398.

Archibald, S. B., D. R. Greenwood, and R. W. Mathewes. 2013. Seasonality, montane beta diversity, and Eocene insects: testing Janzen's dispersal hypothesis in an equable world. Palaeogeogr. Palaeoclimatol. Palaeoecol. 371:1–8.

Blois, J. L., and E. A. Hadly. 2009. Mammalian response to Cenozoic climatic change. Annu. Rev. Earth Planet. Sci. 37:181–208.

Brocklehurst, N., M. O. Day, B. S. Rubidge, and J. Fröbisch. 2017. Olson's Extinction and the latitudinal biodiversity gradient of tetrapods in the Permian. Proc. R. Soc. B Biol. Sci. 284:20170231.

Burin, G., L. R. V. Alencar, J. Chang, M. E. Alfaro, and T. B. Quental. 2019. How well can we estimate diversity dynamics for clades in diversity decline? Syst. Biol. 68:47–62.

Cardillo, M., C. D. L. Orme, and I. P. F. Owens. 2005. Testing for latitudinal bias in diversification rates: an example using new world birds. Ecology 86:2278–2287.

Coiffard, C., and B. Gomez. 2012. Influence of latitude and climate on spread, radiation and rise to dominance of early angiosperms during the Cretaceous in the Northern Hemisphere. Geol. Acta 10:181–188.

Condamine, F. L., F. A. H. Sperling, N. Wahlberg, J.-Y. Rasplus, and G. J. Kergoat. 2012. What causes latitudinal gradients in species diversity? Evolutionary processes and ecological constraints on swallowtail biodiversity. Ecol. Lett. 15:267–277.

Condamine, F. L., J. Rolland, S. Höhna, F. A. H. Sperling, and I. Sanmartín. 2018. Testing the role of the Red Queen and Court Jester as drivers of the macroevolution of Apollo butterflies. Syst. Biol. 67:940–964.

Crame, J. A. 2000. Evolution of taxonomic diversity gradients in the marine realm: evidence from the composition of recent bivalve faunas. Paleobiology 26:188–214.

—. 2001. Taxonomic diversity gradients through geological time. Divers. Distrib. 7:175–189.

—. 2020. Early Cenozoic evolution of the latitudinal diversity gradient. https://doi.org/10.1016/j.earscirev.2020.103090.

Davis, M. P., P. E. Midford, and W. Maddison. 2013. Exploring power and parameter estimation of the BiSSE method for analyzing species diversification. BMC Evol. Biol. 13:38.

De Celis, A., I. Narváez, and F. Ortega. 2019. Spatiotemporal palaeodiversity patterns of modern crocodiles (Crocodyliformes: Eusuchia). Zool. J. Linn. Soc. https://doi.org/10.1093/zoolinnean/zlz038.

Donoghue, M. J. 2008. A phylogenetic perspective on the distribution of plant diversity. Proc. Natl. Acad. Sci. 105:11549–11555.

Dunn, R. R., D. Agosti, A. N. Andersen, X. Arnan, C. A. Bruhl, X. Cerdá, A. M. Ellison, B. L. Fisher, M. C. Fitzpatrick, H. Gibb, et al. 2009. Climatic drivers of hemispheric asymmetry in global patterns of ant species richness. Ecol. Lett. 12:324–333.

Eiserhardt, W. L., F. Borchsenius, C. M. Plum, A. Ordonez, and J.-C. Svenning. 2015. Climate-driven extinctions shape the phylogenetic structure of temperate tree floras. Ecol. Lett. 18:263–272.

Farrell, B., C. Mitter, and D. J. Futuyma. 1992. Diversification at the insect-plant interface: insights from phylogenetics. Bioscience 42:34–42.

Fenton, I. S., P. N. Pearson, T. Dunkley Jones, A. Farnsworth, D. J. Lunt, P. Markwick, and A. Purvis. 2016. The impact of Cenozoic cooling on assemblage diversity in planktonic foraminifera. Philos. Trans. R. Soc. B Biol. Sci. 371:20150224.

FitzJohn, R. G. 2012. Diversitree: comparative phylogenetic analyses of diversification in R. Methods Ecol. Evol. 3:1084–1092.

FitzJohn, R. G., W. P. Maddison, and S. P. Otto. 2009. Estimating traitdependent speciation and extinction rates from incompletely resolved phylogenies. Syst. Biol. 58:595–611.

Frederiksen, N. O. 1988. Sporomorph biostratigraphy, floral changes, and paleoclimatology, Eocene and earliest Oligocene of the eastern Gulf Coast. U.S. Government Printing Office, Washington, DC.

Freyman, W., and S. Hoehna. 2018. The tempo of evolutionary decline in self-compatible plant lineages. Syst. Biol. 68:505–519.

Fritz, S. A., J. Schnitzler, J. T. Eronen, C. Hof, K. Böhning-Gaese, and C. H. Graham. 2013. Diversity in time and space: wanted dead and alive. Trends Ecol. Evol. 28:509–516.

Fuhrman, J. A., J. A. Steele, I. Hewson, M. S. Schwalbach, M. V. Brown, J. L. Green, and J. H. Brown. 2008. A latitudinal diversity gradient in planktonic marine bacteria. Proc. Natl. Acad. Sci. 105:7774–7778.

Ge, D., Z. Wen, L. Xia, Z. Zhang, M. Erbajeva, C. Huang, and Q. Yang. 2013. Evolutionary history of lagomorphs in response to global environmental change. PLoS ONE 8:e59668.

Gelman, A. 2004. Parameterization and Bayesian modeling. J. Am. Stat. Assoc. 99:537–545.

Goldberg, E. E., L. T. Lancaster, and R. H. Ree. 2011. Phylogenetic inference of reciprocal effects between geographic range evolution and diversification. Syst. Biol. 60:451–465.

Hawkins, B. A., J. A. F. Diniz-Filho, C. A. Jaramillo, and S. A. Soeller. 2006. Post-Eocene climate change, niche conservatism, and the latitudinal diversity gradient of New World birds. J. Biogeogr. 33:770–780.

Hillebrand, H. 2004. On the generality of the latitudinal diversity gradient. Am. Nat. 163:192–211.
Huang, D., E. E. Goldberg, and K. Roy. 2015. Fossils, phylogenies, and the challenge of preserving evolutionary history in the face of anthropogenic extinctions. Proc. Natl. Acad. Sci. 112:4909–4914.

Jablonski, D., K. Roy, and J. W. Valentine. 2006. Out of the tropics: evolutionary dynamics of the latitudinal diversity gradient. Science 314:102–106.

Jablonski, D., C. L. Belanger, S. K. Berke, S. Huang, A. Z. Krug, K. Roy, A. Tomasovych, and J. W. Valentine. 2013. Out of the tropics, but how? Fossils, bridge species, and thermal ranges in the dynamics of the marine latitudinal diversity gradient. Proc. Natl. Acad. Sci. 110:10487–10494.

Jaffe, A. L., G. J. Slater, and M. E. Alfaro. 2011. The evolution of island gigantism and body size variation in tortoises and turtles. Biol. Lett. 7:558–561.

Jansson, R., G. Rodríguez-Castañeda, and L. E. Harding. 2013. What can multiple phylogenies say about the latitudinal diversity gradient? A new look at the tropical conservatism, out of the tropics, and diversification rate hypotheses. Evolution 67:1741–1755.

Jenkins, C. N., S. L. Pimm, and L. N. Joppa. 2013. Global patterns of terrestrial vertebrate diversity and conservation. Proc. Natl. Acad. Sci. 110:E2602–E2610.

Joyce, W. G., M. Rabi, J. M. Clark, and X. Xu. 2016. A toothed turtle from the Late Jurassic of China and the global biogeographic history of turtles. BMC Evol. Biol. 16:236.

Kerkhoff, A. J., P. E. Moriarty, and M. D. Weiser. 2014. The latitudinal species richness gradient in New World woody angiosperms is consistent with the tropical conservatism hypothesis. Proc. Natl. Acad. Sci. 111:8125–8130.

Kindlmann, P., I. Schödelbauerová, and A. F. G. Dixon. 2007. Inverse latitudinal gradients in species diversity. Pp. 246–257 *in* D. Storch, P. Marquet, and J. Brown, eds. Scaling biodiversity. Cambridge Univ. Press, Cambridge, U.K.

Kreft, H., and W. Jetz. 2007. Global patterns and determinants of vascular plant diversity. Proc. Natl. Acad. Sci. 104:5925–5930.

Krug, A. Z., and M. E. Patzkowsky. 2007. Geographic variation in turnover and recovery from the Late Ordovician mass extinction. Paleobiology 33:435–454.

Labandeira, C. C., and E. D. Currano. 2013. The fossil record of plant-insect dynamics. Annu. Rev. Earth Planet. Sci. 41:287–311.

Latham, R. E., and R. E. Ricklefs. 1993. Continental comparisons of temperate-zone tree species diversity. Pp. 294–314 *in* R. E. Ricklefs and D. Schlüter, eds. Species diversity in ecological communities: historical and geographical perspectives. Univ. of Chicago Press, Chicago.

Leprieur, F., P. Descombes, T. Gaboriau, P. F. Cowman, V. Parravicini, M. Kulbicki, C. J. Melian, C. N. De Santana, C. Heine, D. Mouillot, et al. 2016. Plate tectonics drive tropical reef biodiversity dynamics. Nat. Commun. 7:11461.

Leslie, A. B., J. M. Beaulieu, H. S. Rai, P. R. Crane, M. J. Donoghue, and S. Mathews. 2012. Hemisphere-scale differences in conifer evolutionary dynamics. Proc. Natl. Acad. Sci. 109:16217–16221.

Maddison, W. P., P. E. Midford, and S. P. Otto. 2007. Estimating a binary character's effect on speciation and extinction. Syst. Biol. 56:701–710.

Mannion, P. D., R. B. J. Benson, P. Upchurch, R. J. Butler, M. T. Carrano, and P. M. Barrett. 2012. A temperate palaeodiversity peak in Mesozoic dinosaurs and evidence for Late Cretaceous geographical partitioning. Glob. Ecol. Biogeogr. 21:898–908.

Mannion, P. D., P. Upchurch, R. B. Benson, and A. Goswami. 2014. The latitudinal biodiversity gradient through deep time. Trends Ecol. Evol. 29:42–50.

Mannion, P. D., R. B. J. Benson, M. T. Carrano, J. P. Tennant, J. Judd, and R. J. Butler. 2015. Climate constrains the evolutionary history and biodiversity of crocodylians. Nat. Commun. 6:8438.

Marcot, J. D., D. L. Fox, and S. R. Niebuhr. 2016. Late Cenozoic onset of the latitudinal diversity gradient of North American mammals. Proc. Natl. Acad. Sci. 113:7189–7194.

Markwick, P. J. 1998. Fossil crocodilians as indicators of Late Cretaceous and Cenozoic climates: implications for using palaeontological data in reconstructing palaeoclimate. Palaeogeogr. Palaeoclimatol. Palaeoecol. 137:205–271.

Martin, P. S., and R. G. Klein. 1989. Quaternary extinctions: a prehistoric revolution. Univ. of Arizona Press, Tucson, AZ.

Mayr, G. 2016. The early Eocene birds of the Messel fossil site: a 48 million-year-old bird community adds a temporal perspective to the evolution of tropical avifaunas. Biol. Rev. 92:1174–1188.

Meseguer, A. S., J. M. Lobo, J. Cornuault, D. Beerling, B. R. Ruhfel, C. C. Davis, E. Jousselin, and I. Sanmartín. 2018. Reconstructing deep-time paleoclimate legacies in the clusioid Malpighiales unveils their role in the evolution and extinction of the boreotropical flora. Glob. Ecol. Biogeogr. 27:616–628.

Mittelbach, G. G., D. W. Schemske, H. V. Cornell, A. P. Allen, J. M. Brown, M. B. Bush, S. P. Harrison, A. H. Hurlbert, N. Knowlton, H. A. Lessios, et al. 2007. Evolution and the latitudinal diversity gradient: speciation, extinction and biogeography. Ecol. Lett. 10:315–331.

Moore, B. R., S.a. Smith, R. H. Ree, and M. J. Donoghue. 2008. Incorporating fossil data in biogeographic inference: a likelihood approach. Proceedings of the 52nd Annu. Syst. Symposium. Available at http://www.phylo.org/pdf_docs/2006_MoBot_biogeo.pdf.

Morley, R. J. 2011. Cretaceous and Tertiary climate change and the past distribution of megathermal rainforests. Pp. 1–31 *in* M. B. Bush and J.-R. Flenley, eds. tropical rainforest responses to climatic change. Springer, Berlin, Germany.

Morlon, H., T. L. Parsons, and J. B. Plotkin. 2011. Reconciling molecular phylogenies with the fossil record. Proc. Natl. Acad. Sci. 108:16327–16332.

Nee, S., R. M. May, and P. H. Harvey. 1994. The reconstructed evolutionary process. Philos. Trans. R. Soc. B Biol. Sci. 344:305–311.

Nicholson, D. B., P. A. Holroyd, R. B. J. Benson, and P. M. Barrett. 2015. Climate-mediated diversification of turtles in the Cretaceous. Nat. Commun. 6:7848.

Nicholson, D. B., P. A. Holroyd, P. Valdes, and P. M. Barrett. 2016. Latitudinal diversity gradients in mesozoic non-marine turtles. R. Soc. Open Sci. 3:160581.

Novotny, V., P. Drozd, S. E. Miller, M. Kulfan, M. Janda, Y. Basset, and G. D. Weiblen. 2006. Why are there so many species of herbivorous insects in tropical rainforests? Science 313:1115–1118.

Oaks, J. R. 2011. A time-calibrated species tree of crocodylia reveals a recent radiation of the true crocodiles. Evolution 65:3285–3297.

Ogg, J. G., F. M. Gradstein, and A. G. Smith. 2004. A geologic time scale 2004. Cambridge Univ. Press, Cambridge, U.K.

Peralta-Medina, E., and H. J. Falcon-Lang. 2012. Cretaceous forest composition and productivity inferred from a global fossil wood database. Geology 40:219–222.

Pereira, A. G., J. Sterli, F. R. R. Moreira, and C. G. Schrago. 2017. Multilocus phylogeny and statistical biogeography clarify the evolutionary history of major lineages of turtles. Mol. Phylogenet. Evol. 113:59–66.

Peterson, A. T., J. Soberón, and V. Sánchez-Cordero. 1999. Conservatism of ecological niches in evolutionary time. Science 285:1265–1267.

Pie, M. R., L. L. F. Campos, A. L. S. Meyer, and A. Duran. 2017. The evolution of climatic niches in squamate reptiles. Proc. R. Soc. B Biol. Sci. 284:1858.

Powell, M. G. 2007. Latitudinal diversity gradients for brachiopod genera during late Palaeozoic time: links between climate, biogeography and evolutionary rates. Glob. Ecol. Biogeogr. 16:519–528.
Powell, M. G., V. P. Beresford, and B. A. Colaianne. 2012. The latitudinal position of peak marine diversity in living and fossil biotas. J. Biogeogr. 39:1687–1694.

Pulido-Santacruz, P., and J. T. Weir. 2016. Extinction as a driver of avian latitudinal diversity gradients. Evolution 70:860–872.

Pybus, O. G., and P. H. Harvey. 2000. Testing macro-evolutionary models using incomplete molecular phylogenies. Proc. R. Soc. Lond. Ser. B Biol. Sci. 267:2267–2272.

Pyron, R. A. 2014. Temperate extinction in squamate reptiles and the roots of latitudinal diversity gradients. Glob. Ecol. Biogeogr. 23:1126–1134.

Pyron, R. A., and J. J. Wiens. 2013. Large-scale phylogenetic analyses reveal the causes of high tropical amphibian diversity. Proc. R. Soc. Lond. B Biol. Sci. 280:20131622.

Rabosky, D. L. 2010. Extinction rates should not be estimated from molecular phylogenies. Evolution 6:1816–1824.

Rabosky, D. L., J. Chang, P. O. Title, P. F. Cowman, L. Sallan, M. Friedman, K. Kaschner, C. Garilao, T. J. Near, and M. Coll. 2018. An inverse latitudinal gradient in speciation rate for marine fishes. Nature 559:392–395.

Rambaut, A., A. J. Drummond, D. Xie, G. Baele, and M. A. Suchard. 2018. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67:901–904.

Reddin, C. J., Á. T. Kocsis, and W. Kiessling. 2019. Climate change and the latitudinal selectivity of ancient marine extinctions. Paleobiology 45:70–84.

Ree, R. H., and I. Sanmartin. 2009. Prospects and challenges for parametric models in historical biogeographical inference. J. Biogeogr. 36:1211–1220.

Ree, R. H., and S. A. Smith. 2008. Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. Syst. Biol. 57:4–14.

Ricklefs, R. E. 2006. Global variation in the diversification rate of passerine birds. Ecology 87:2468–2478.

Rodrigues, J. F. M., and J. A. F. Diniz-Filho. 2016. Ecological opportunities, habitat, and past climatic fluctuations influenced the diversification of modern turtles. Mol. Phylogenet. Evol. 101:352–358.

Rodrigues, J. F. M., M. Á. Olalla-Tárraga, J. B. Iverson, T. S. Akre, and J. A. F. Diniz-Filho. 2017. Time and environment explain the current richness distribution of non-marine turtles worldwide. Ecography 40:1402–1411.

Rolland, J., F. L. Condamine, F. Jiguet, and H. Morlon. 2014. Faster speciation and reduced extinction in the tropics contribute to the mammalian latitudinal diversity gradient. PLoS Biol. 12:e1001775.

Rolland, J., F. L. Condamine, C. R. Beeravolu, F. Jiguet, and H. Morlon. 2015. Dispersal is a major driver of the latitudinal diversity gradient of Carcinova. Glob. Ecol. Biogeogr. 24:1059–1071.

Romdal, T. S., M. B. Araújo, and C. Rahbek. 2013. Life on a tropical planet: niche conservatism and the global diversity gradient. Glob. Ecol. Biogeogr. 22:344–350.

Rose, P. J., D. L. Fox, J. Marcot, and C. Badgley. 2011. Flat latitudinal gradient in Paleocene mammal richness suggests decoupling of climate and biodiversity. Geology 39:163–166.

Roy, K., and J. M. Pandolfi. 2005. Responses of marine species and ecosystems to past climate change. Pp. 160–175 in E. Thomas, ed. Climate change and biodiversity. Yale Univ. Press, New Haven, CT.

Sanmartín, I., and A. S. Meseguer. 2016. Extinction in phylogenetics and biogeography: from timetrees to patterns of biotic assemblage. Front. Genet. 7:35.

Saupe, E. E., A. Farnsworth, D. J. Lunt, N. Sagoo, K. V. Pham, and D. J. Field. 2019a. Climatic shifts drove major contractions in avian latitudinal dis-

tributions throughout the Cenozoic. Proc. Natl. Acad. Sci. U.S.A. 116:12895–12900.

Saupe, E. E., C. E. Myers, A. T. Peterson, J. Soberón, J. Singarayer, P. Valdes, and H. Qiao. 2019b. Spatio-temporal climate change contributes to latitudinal diversity gradients. Nat. Ecol. Evol. 3:1419–1429.

Shaffer, H. B., E. McCartney-Melstad, T. J. Near, G. G. Mount, and P. Q. Spinks. 2017. Phylogenomic analyses of 539 highly informative loci dates a fully resolved time tree for the major clades of living turtles (Testudines). Mol. Phylogenet. Evol. 115:7–15.

Shiono, T., B. Kusumoto, M. Yasuhara, and Y. Kubota. 2018. Roles of climate niche conservatism and range dynamics in woody plant diversity patterns through the Cenozoic. Glob. Ecol. Biogeogr. 27:865–874.

Silvestro, D., J. Schnitzler, L. H. Liow, A. Antonelli, and N. Salamin. 2014. Bayesian estimation of speciation and extinction from incomplete fossil occurrence data. Syst. Biol. 63:349–367.

Silvestro, D., B. Cascales-Miñana, C. D. Bacon, and A. Antonelli. 2015. Revisiting the origin and diversification of vascular plants through a comprehensive Bayesian analysis of the fossil record. New Phytol. 207:1469–8137.

Silvestro, D., R. C. M. Warnock, A. Gavryushkina, and T. Stadler. 2018. Closing the gap between palaeontological and neontological speciation and extinction rate estimates. Nat. Commun. 9:5237.

Silvestro, D., N. Salamin, A. Antonelli, and X. Meyer. 2019. Improved estimation of macroevolutionary rates from fossil data using a Bayesian framework. Paleobiology 45:546–570.

Siqueira, A. C., L. G. R. Oliveira-Santos, P. F. Cowman, and S. R. Floeter. 2016. Evolutionary processes underlying latitudinal differences in reef fish biodiversity. Glob. Ecol. Biogeogr. 25:1466–1476.

Sluijs, A., S. Schouten, M. Pagani, M. Woltering, H. Brinkhuis, J. S. S. Damsté, G. R. Dickens, M. Huber, G.-J. Reichart, R. Stein, et al. 2006. Subtropical Arctic Ocean temperatures during the Palaeocene/Eocene thermal maximum. Nature 441:610–613.

Smith, R. Y. 2012. and et al. Early Eocene plant diversity and dynamics in the Falkland flora, Okanagan Highlands, British Columbia, Canada. Palaeobiodivers. Palaeoenviron. 92:309–328.

Smith, S. A. 2009. Taking into account phylogenetic and divergence-time uncertainty in a parametric biogeographical analysis of the Northern Hemisphere plant clade Caprifoliaceae. J. Biogeogr. 36:2324–2337.

Spriggs, E. L., W. L. Clement, P. W. Sweeney, S. Madrinan, E. J. Edwards, M. J. Donoghue, and P. New. 2015. Temperate radiations and dying embers of a tropical past: the diversification of Viburnum. New Phytol. 207:1469–8137.

Stebbins, G. L. 1974. Flowering plants: evolution above the species level. Harvard Univ. Press, Cambridge, MA.

Stephens, P. R., and J. J. Wiens. 2003. Explaining species richness from continents to communities: the time-for-speciation effect in emydid turtles. Am. Nat. 161:112–128.

Tütken, T., and J. Absolon. 2015. Late Oligocene ambient temperatures reconstructed by stable isotope analysis of terrestrial and aquatic vertebrate fossils of Enspel, Germany. Palaeobiodivers. Palaeoenviron. 95:17–31.

Waterson, A. M., D. N. Schmidt, P. J. Valdes, P. A. Holroyd, D. B. Nicholson, A. Farnsworth, and P. M. Barrett. 2016. Modelling the climatic niche of turtles: a deep-time perspective. Proc. R. Soc. B Biol. Sci. 283:20161408.

Weir, J. T., and D. Schlüter. 2007. The latitudinal gradient in recent speciation and extinction rates of birds and mammals. Science 315:1574–1576.

Wiens, J., and M. Donoghue. 2004. Historical biogeography, ecology and species richness. Trends Ecol. Evol. 19:639–644.
Wiens, J. J. 2007. Global patterns of diversification and species richness in amphibians. Am. Nat. 70:S86–S106.

Wiens, J. J., D. D. Ackerly, A. P. Allen, B. L. Anacker, L. B. Buckley, H. V. Cornell, et al. 2010. Niche conservatism as an emerging principle in ecology and conservation biology. Ecology letters 13:1310–1324.

Wilf, P., C. C. Labandeira, K. R. Johnson, and N. R. Cuneo. 2005. Richness of plant-insect associations in Eocene Patagonia: a legacy for South American biodiversity. Proc. Natl. Acad. Sci. 10:8944–8948.

Willig, M. R., D. M. Kaufman, and R. D. Stevens. 2003. Latitudinal gradients of biodiversity: pattern, process, scale, and synthesis. Annu. Rev. Ecol. Evol. Syst. 34:273–309.

Xing, Y., R. E. Onstein, R. J. Carter, T. Stadler, and H. Peter Linder. 2014. Fossils and a large molecular phylogeny show that the evolution of

species richness, generic diversity, and turnover rates are disconnected. Evolution 68:2821–2832.

Zachos, J. C., G. R. Dickens, and R. E. Zeebe. 2008. An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics. Nature 451:279–283.

Ziegler, A., G. Eshel, P. M. Rees, T. Rothfus, D. Rowley, and D. Sunderlin. 2003. Tracing the tropics across land and sea: Permian to present. Lethaia 36:227–254.

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 S1.** Geographic distribution of extant turtle species based on the *Reptile Database* and the *IUCN Red List*.

**Table S2.** Geographic distribution of extant squamate species based on the *Reptile Database* and the *IUCN Red List*.

**Table S3.** Geographic distribution of extant crocodile species.

**Table S4–6.** Could be found on FigShare https://doi.org/10.6084/m9.fishshare.10078955.

**Table S7.** Speciation and extinction rates estimated from extant phylogenies and stratigraphic ranges of Crocodyles, Testudines and Lepidosauria under the equal, 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 S2a.** Latitudinal time-dependent diversification patterns of turtles inferred with BiSSE.td showing the diversification dynamics of turtles during the greenhouse (left panel), the climatic transition (central panel), and the coldhouse (right panel) intervals for Holarctic and Equatorial species.

**Figure S2b.** Latitudinal time-dependent diversification patterns of squamates inferred with BiSSE.td showing the diversification dynamics of turtles during the greenhouse (left panel), the climatic transition (central panel), and the coldhouse (right panel) intervals for Holarctic and Equatorial species.

**Figure S3a.** Latitudinal time-dependent diversification patterns of turtles inferred with BiSSE.td showing the diversification dynamics during the greenhouse (left panel), the climatic transition (central panel), and the coldhouse (right panel) intervals for Holarctic and Equatorial species. In this analysis, the climatic transition begins 66 Ma and ends 23 Ma.

**Figure S3b.** Latitudinal time-dependent diversification patterns of squamates inferred with BiSSE.td showing the diversification dynamics during the greenhouse (left panel), the climatic transition (central panel), and the coldhouse (right panel) intervals for Holarctic and Equatorial species.

**Figure S4a.** Latitudinal time-dependent diversification patterns of turtles inferred with BiSSE.td showing the diversification dynamics of turtles during the greenhouse (left panel), the climatic transition (central panel), and the icehouse (right panel) intervals for Holarctic and Equatorial species.

**Figure S4b.** Latitudinal time-dependent diversification pattern of squamates inferred with BiSSE.td showing the diversification dynamics of squamates during the greenhouse (left panel), the climatic transition (central panel), and icehouse (right panel) intervals 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.