A global phylogeny of turtles reveals a burst of climate-associated diversification on continental margins

Phylogeny.

We developed a well-sampled phylogenetic hypothesis for the extant Testudines. Our taxon sampling includes all 14 living families, 90 of the 92 living genera (98%), and 279 of the 348 living species (80%) of turtles (Table 1). When possible, we included two (occasionally more) individuals per species to allow an assessment of species-level monophyly, resulting in a full dataset of 591 individual turtles.

The tree recapitulates several previously identified relationships in the deep history of the group (Figs. 1 and 2 and SI Appendix, Fig. S1) (5, 19–23). We recovered the early divergence between the Cryptodira (hidden-necked turtles) and Pleurodira (side-necked turtles), with maximal support values (Bayesian posterior probabilities = 1); within Cryptodira, we recovered softshell turtles (Trionychia) as the sister group to Durocryptodira (hard-shelled cryptodires); and, within Pleurodira, Chelidae and Pelomedusoides are reciprocally monophyletic (Figs. 1 and 2 and SI Appendix, Fig. S1). We recovered all 14 currently recognized families as monophyletic, with interfamilial relationships that are congruent with recent molecular analyses. At more recent levels of divergence, we recovered several notable relationships that we discuss in greater detail in the SI Appendix.

Divergence Times and Diversification.

Crown turtles likely originated in the Late Triassic ∼210 mya (4, 5). The paleontological record for the group is rich and one of the best known of any major vertebrate lineage (10, 51, 52). Using a series of well-characterized fossil calibrations from Joyce et al. (4) and a version of the molecular dataset that was reduced to one individual per species, we coestimated phylogeny and divergence times for the living turtles (Figs. 1 and 2). This analysis indicates that the two major lineages within the clade, Cryptodira and Pleurodira, diverged in the Late Triassic (median 208 mya; 95% highest posterior density [HPD] 183 mya to 236 mya), followed by the Early Jurassic divergence between the softshell Trionychia (Carettochelyidae and Trionychidae) and the remaining cryptodires 182 mya (95% HPD 161 mya to 200 mya). Most recognized families diverged in the Late Cretaceous or Early Paleogene (4, 5, 46, 53, 54).

Some earlier analyses of the fossil record suggest that climate has played a role in shaping diversity through geological time (10, 47, 49, 50), albeit in different ways and at different scales (also see ref. 48). The number of fossil genera increased in both number and total occupied latitudinal range in the lower Cretaceous (145 mya to 100 mya) and again following the Cretaceous Thermal Maximum in the Late Cretaceous 85 mya to 90 mya (10), raising the possibility that diversification was mediated by an increase in climatically suitable geographic area. Cranial disparity of turtles also accumulated gradually across this time period and accelerated toward the present, although this analysis did not detect a particularly strong correlation with climate (48).

These paleontological results provide an interesting comparison to the timing of diversification we found in the phylogeny of extant turtles and highlight the significant role that extinction has played in shaping their present-day diversity. Our analysis identified no detectable change in diversification rate during the Cretaceous (Fig. 2B). While the fossil record documents a buildup in diversity during the Middle Cretaceous, our divergence time analysis indicates that most living turtle diversity is much younger, with most speciation events occurring within the last 50 mya, well after the end of the Cretaceous (65 mya). The increased diversification observed in the fossil record by Nicholson et al. (10) is at least partly associated with lineages that have no living descendants and therefore do not register in the molecular phylogeny of extant diversity. Taken together, these results document a dramatic turnover in turtle diversity through geological time, with both recent speciation and more ancient extinction contributing to the tree we observe today.

The new chronogram provides evidence for an increasing diversification rate beginning roughly 50 mya during the Eocene (Fig. 2 B and C). The estimated net diversification rate (speciation minus extinction) roughly triples beginning in the Eocene and remains elevated toward the present day (Fig. 2B). Given recent concerns about the identifiability of diversification shift models (55), we also confirmed that a qualitatively similar, but perhaps somewhat earlier, shift was recovered using a state-dependent modeling framework (Fig. 2D). Bayes factors (BF) among alternative rate shift models indicate positive (2 ln BF greater than 2 and less than 6; ref. 56) support for a shift that occurs at the Eocene−Oligocene transition 34 mya (Fig. 2C). This period coincides with large-scale declines in global sea level, which began falling from a peak of more than 150 m above modern levels ∼50 mya, and experienced a dramatic decline by the end-Eocene (57–60). It also raises the possibility that the increasing availability of suitable habitat along continental margins contributed to increased diversification of modern turtles. As lineages expanded into newly exposed continental margins, they may have rapidly built up species numbers through allopatric isolation in low-lying areas or among new isolated river basins. These environmental changes have been shown to be an important force in structuring diversity of other organisms (61). Previous work has also implicated incised river basins and sea level changes as a potential factor that drives a region of high endemicity in turtles (38, 62).

To test this hypothesis, we asked whether rates of diversification were faster on newly exposed continental margins, as would be expected if declining sea levels created ecological opportunity that led to speciation. Using the most recently available estimates of species ranges (11) and reconstructed paleocoastlines for the Eocene−Oligocene transition (59, 63), we classified each species as currently occupying 1) low-lying coastal areas that had not yet emerged before the Eocene−Oligocene transition (“coastal”), 2) inland habitats exposed and available prior to the transition (“inland”), or 3) both (“widespread”) (Fig. 3). We then fit a series of 10 models of geography-dependent (or independent) speciation and extinction using the new chronogram and range data (Table 2). These models vary from scenarios where diversification rate changes may occur, but are independent of geographical range, to scenarios where diversification rates depend on geographical range either alone or in combination with other unobserved (“hidden”) traits.

We found that the best-fitting models were those where geographic range was linked to variation in diversification rates (Table 2). Parameter estimates from these models indicate that speciation rates were higher in coastal compared to inland or widespread range lineages (Fig. 4 and SI Appendix, Fig. S2). The two best-fitting models together make up 0.951 of the Akaike information criterion (AIC) weight (0.535 and 0.416, respectively) and vary only in whether extirpation is modeled independently and whether jumps in ranges are allowed (i.e., from “coastal” to “inland,” rather than evolving through a “widespread” intermediate ancestor). The best-fitting model is the simpler of the two and includes neither separate extirpation nor jumps. Conversely, all the models we tested that treat range evolution and diversification as independent processes were much poorer fits. The best of this group had an AIC 4.79 units lower than the pair of best-fitting models, and an AIC weight of 0.048. All remaining models in the set had AIC weights below 0.0001. These results are not sensitive to uncertainty in topology or divergence time estimates, as repeated analyses across a sample of 100 chronograms drawn from the posterior recovered qualitatively similar results (Table 2).

While models that link geographic range to diversification fit the data much better than those that do not, we also found substantial model uncertainty between the two best-fitting models. Rather than rely on results from only the single “best” model, we averaged rates across all 10 models using their AIC weights. The weighted average parameter estimates again indicate speciation rates that are threefold higher in coastal regions than inland regions (Fig. 4) and increase through time regardless of the modeling framework (compare Fig. 2 B and D). Ancestral range estimates (as weighted averages across 10 models) indicate that nearly all coastal species transitioned to coastal areas after the end of the Eocene (dashed line in Fig. 4A), and that this pattern is common across families and global regions. If newly exposed land areas served as new habitat for chelonian range expansion and speciation, it stands to reason that this signature may be detectable in ancestral range estimates. Here, the models are agnostic about timing of transitions among areas, but they still recover transitions into and speciation within newly exposed continental margins during this time period. The exceptions to this are the four lineages (within the Carretochelyidae, Geoemydidae, Podocnemidae, and Dermatemydidae) that have moderate weight for coastal ranges well before the Eocene−Oligocene transition (Fig. 4A). Three of these are single-taxon-long branches where the reconstructions are necessarily equivocal. The remaining lineage, the diverse Geoemydidae, had equivocal support for coastal ranges before the Eocene−Oligocene transition, which increases to strong support at the time of the transition (Fig. 4A, nodes 1 to 8). The model-averaged reconstructed state probabilities for these nodes are provided in Table 3. This uncertainty could arise from the geographic distribution of the Geoemydidae. Most of the family’s diversity occurs on the continental mainland and the islands of Southeast Asia, an area with an extremely complicated geological history. This may lead to a relatively greater level of uncertainty in paleogeographic reconstructions within the region, which would then propagate as error in our categorization of present-day species ranges. Because of this, reconstructions for these early geoemydid nodes could reasonably be interpreted as either modest support against our hypothesis or the result of error in classifying the ranges of species, given the complexity and uncertainty surrounding the area’s geology. Whichever the explanation, if we restrict our attention to range reconstructions where the data provide a strong statistical signal, all lineages including Geoemydidae conform to the hypothesis.

Reconciling Fossils and Phylogenies.

Some paleontological studies of turtle diversity have recovered an association between climate and diversity, with warm wet regions supporting higher diversity. The fossil record also documents a consistent decrease in diversity over the Eocene−Oligocene transition coincident with global cooling, at least at regional scales in western North America (50), a trend that sharply contrasts with the increasing rate of speciation and accumulation of diversity that we recover during the same period globally. The resolution of this contrast lies in the differences in geographic focus between approaches, and highlights the value of large-scale diversification analyses. In particular, some of the most influential paleontological studies have focused on the rich fossil beds that occur in western North America (in the vicinity of present-day Wyoming), an area that underwent aridification across the Eocene−Oligocene transition as sea levels declined. Given this, the paleontological findings actually accord well with our findings. Declining sea levels drove climate-associated habitat loss in some areas, including western North America where particularly significant fossil beds for turtles now occur. However, it expanded habitat in other regions, including present-day coastal margins that contain far fewer significant fossil-bearing sites. In accordance with these results, the paleontological studies also suggest that differences in patterns of diversification are regional. The same drying that resulted in aridification and extinction in western North America resulted in exposed continental margins and speciation in other areas of the world. As a case in point, southeastern North America, which is now a global hotspot of turtle diversity (64), is partly assembled from the same lineages (e.g., Apalone, Chelydra, Chrysemys, and other Emydidae) that disappeared from western North America during the Eocene−Oligocene transition. Relying on fossils or molecular information alone would miss these different regional signals, yielding an incomplete picture of diversification dynamics at odds with the more complex pattern that they jointly reveal.

Taxon and Data Sampling.

We follow the taxonomy from the most recent compilation from the Turtle Taxonomy Working Group, which recognizes 348 extant species organized into 92 genera and 14 families (11). Our taxon sampling comprises 279 of these 348 species or 80% of the extant turtles of the world (Table 1). We include representatives of all 14 families and 90 of the 92 genera (98%). When possible, we included two (occasionally more) individuals/species to allow a minimal assessment of species-level monophyly, which led to a full dataset of 591 individual turtles (Dataset S1). We also included GenBank sequences from chicken (Gallus gallus) and alligator (Alligator mississippiensis) as outgroups (66, 67). We generated new nuclear DNA (nuDNA) sequence data from 15 nuclear loci (AHR, AIING, BDNF, BMP2, HMGB2, HNFL, NB22519, PSMC1, PAX1P1, R35, RAG-1, TB01, TB29, TB73, ZFHX1B; refs. 22, 26, 37, 64, 66–72). To do so, we extracted DNA from blood or soft tissue samples using a salt extraction protocol (73) and amplified PCR products for all loci in 20-µL-volume reactions using the locus-specific annealing temperatures, extension times, and primers for these markers developed in our earlier studies (23, 30, 74). PCR products were checked on 1% agarose gels and sequenced in both directions using Sanger sequencing by Beckman Coulter Genomics (https://www.beckman.com).

Data Matrix.

We assembled a data matrix containing 6,029 new sequences generated for this study (68%), 2,004 additional sequences from GenBank (∼23%, most of which were submitted previously by our group), and 832 missing sequences (∼9.4%). The 15 nuDNA markers ranged in size from 593 bp to 1,488 bp (mean = 904 bp), and totaled 13,559 bp of aligned sequence data for 593 individuals, including outgroups. All sequences are available in GenBank (see Dataset S1 for accession numbers), and the data matrix is available in the Dryad repository associated with this study (https://doi.org/10.5061/dryad.jh9w0vt8w) (75).

Phylogenetic Analyses.

We used Geneious v5.1 (76) to translate exons and check for the presence of unexpected stop codons, frame shifts, and other indicators that pseudogenes or other unwanted products may have been unintentionally sequenced. Alignments were carried out using the MAFFT software (77) plugin for Geneious (76). We used jModelTest (78) to select models of nucleotide substitution for each locus prior to the phylogenetic analyses. For the 593-taxon dataset (591 turtles plus two outgroups), we performed Bayesian phylogenetic analyses on the entire 15-locus concatenated and partitioned dataset using MrBayes v3.2 (79). Bayesian analyses consisted of 10 independent runs, each comprising one cold chain and three incrementally heated chains (temperature parameter = 0.1), that ran for up to 10,000,000 generations each, recording the state of the cold chain every 1,000 generations. We examined the Markov chain Monte Carlo (MCMC) samples in Tracer (80) and Are We There Yet? (81) to ensure that all chains appeared to be mixing adequately, that all were sampling from the same target distribution, and that a large sample of all parameters (Effective Sample Size > 200 after burn-in) in the analysis had been obtained. We discarded the first 25% of samples as burn-in, provided that the chains had apparently reached stationarity prior to that point, and combined and summarized the posterior samples using LogCombiner and TreeAnnotator v1.8.2 (82).

Divergence Time Estimates.

We estimated divergence times with BEAST v.1.8.2 (82) using the uncorrelated lognormal relaxed clock model (83) with a diffuse gamma prior on both the mean clock rate and SD, the birth−death stochastic branching process prior, and default priors for remaining parameters of the model. We reduced the dataset to one representative per species or subspecies (288 total) to provide a better match to the birth−death prior (which does not account for extensive intraspecific sampling) and used 22 uniform fossil calibration priors that have been previously defined and justified (4) to calibrate the clock. We estimated divergence times under the settings and priors described above using independent Hasegawa-Kishino-Yano substitution models for each partition (which was necessary to obtain adequate MCMC performance). We performed nine replicated analyses each using a single MCMC chain that was allowed to run for at least 50,000,000 generations, recording the state of the chain every 5,000 generations. We examined the MCMC samples as above to check for adequate mixing, convergence, and sampling, removed the first 25% of samples as burn-in, and combined results from the nine runs. We then calculated the mean divergence time and corresponding 95% HPD interval for each node in the phylogeny across the replicates.

Diversification Analysis.

We investigated the role that changing rates of speciation and extinction may play in explaining the standing diversity of turtles. Some previous studies that focus on fossil species have found increased periods of diversification in association with changes in global temperatures (10, 47). We investigated whether this signal emerged in the phylogenetic tree of extant species by investigating models of tree-wide shifts in diversification rate. Using reversible jump MCMC (rjMCMC), we averaged over models of temporal shifts in speciation and extinction rate using the maximum clade credibility time-calibrated phylogeny. These analyses were implemented using the CoMET functions of the R package TESS (84, 85). We used empirical hyperpriors on the rates of speciation and extinction, set the prior on the expected number of diversification rate shifts to 1, and assumed uniform taxon sampling for the 80% of species that we sampled. We started the MCMC using the software’s autotuning and autostopping rules to terminate the chain after 500 effective samples had been achieved, repeated the analysis four times, and checked that each individual run had large effective sample sizes, a stable Geweke statistic near zero across all time periods in the analysis, and the replicated analyses all converged to similar results. We ran additional analyses varying the lognormal priors on speciation and extinction rates, the expected number of rate shifts, and the taxon sampling probabilities to ensure that our qualitative results were not sensitive to these settings. We repeated these analyses targeting the joint prior probability distribution to ensure that the posterior samples strongly diverged from the prior and that the data (the tree and divergence times) were informative with respect to the parameters of the model. We then allowed the rjMCMC to explore models that allowed for mass extinction events, identifying no strong evidence for these under any settings. We repeated these analyses on a sample of 100 trees drawn randomly from the posterior distribution of chronograms to check that the results were not sensitive to uncertainty in topology or divergence times. Finally, because recent work has raised concerns about the reliability of analyses under such models (55), we treated these results cautiously and rechecked them using the Geographic Hidden State Speciation and Extinction (GeoHiSSE) modeling framework.

We further investigated whether diversification dynamics are associated with geographic ranges by fitting a set of 10 models from the GeoSSE and GeoHiSSE modeling framework (86) to the tree and range data for turtles, assuming uniform taxon sampling for the 81% of turtles that we sampled in coastal areas, 90% in widespread areas, and 77% for inland areas. These models vary according to whether diversification is linked to range, whether range contractions are treated separately from extirpation, and whether the models allow diversification rate variation associated with unobserved (or hidden) states (Table 2). We fit the models to the tree and range data, reconstructed ancestral range under each model, and calculated AICs to identify the best supported models in the set using the R package HiSSE (87). We then averaged estimates of diversification parameters and ancestral ranges across the set of models using AIC weights following Caetano et al. (86). We repeated these analyses across a posterior sample of 100 chronograms to check that the results were robust to uncertainty in topology and timing of divergence. To check the results of the CoMET analysis, we also examined mean diversification rate variation across time by dividing the tree into 20 equally sized time bins and calculating the average net diversification rate across the nodes in each bin.