Turtle species extinction across the Cretaceous/Paleogene boundary

Abstract

The last mass extinction event some 66 million years ago at the Late Cretaceous/Paleogene boundary caused the extinction of many clades, including the non-avian dinosaurs. Turtles, as well as several other vertebrate clades, survived. However, the debate about whether the diversity of turtles was affected during this event is still ongoing. Here, I calculate a global turtle diversity curve at the species level that shows that the diversity of turtle species was already in decline since the Campanian, before the extinction event, and was further reduced during the Danian. The sample coverage of turtle occurrences at the stage level is also calculated and discussed.

Non-technical Summary

Around 66 million years ago, a massive extinction event wiped out many species, including the dinosaurs. However, some animals, like turtles, managed to survive. Scientists have been debating whether this extinction event affected the variety of turtle species. This study creates a detailed curve showing the number of turtle species over time. It was found that the variety of turtles was already decreasing before the extinction event and continued to drop afterward. This suggests that the extinction event had a significant impact on turtle diversity, which had already been in decline.

Introduction

Since the seminal work of Raup and Sepkoski (1982) and the definition of the “big five” mass extinction events, we are now within historical phase 3 (sensu Marshall 2022), working on the recognition and understanding of diversity events in the history of life on Earth. During this phase, a number of analytical tools and techniques are at our disposal, including a constantly growing database of fossil occurrences, the Paleobiology Database (PBDB). Among the five or more mass extinction events, the most recent at the end of the Cretaceous (Cretaceous/Paleogene [K/Pg] boundary) has achieved mainstream status because of the extinction of the non-avian dinosaurs and the asteroid impact. In recent years, scientists have focused on studying the effect of this mass extinction event on different groups of vertebrates and invertebrates and explored possible declines in diversity before the impact of the asteroid. For example, some dinosaur clades had declined since the Campanian, before their final extinction (Brusatte et al. 2015). As these clades with reduced diversity involved mainly large-bodied herbivores, their decline could have affected the stability of local ecosystems (Brusatte et al. 2015). Bayesian models suggest a similar overall result because of higher extinction rates than speciation rates during the Late Cretaceous (Sakamoto et al. 2016; Condamine et al. 2021). These results indicate that dinosaur diversity responded to large-scale environmental changes, mainly because of falling temperatures during the cooler greenhouse climate of the Campanian–Maastrichtian, although other factors like geotectonic plate configurations, geological changes, sea-level fluctuations, or floristic turnover should be taken into account (Condamine et al. 2021 and references therein).

The effect of the last mass extinction on turtle diversity is still a matter of debate. The prevailing view is that the diversity of turtles had been largely unaffected or just slightly affected. Hutchison and Archibald (1986) studied turtle diversity with raw field data across the K/Pg boundary in eastern Montana (USA) and found only 16% extinction. Additional evidence from the same region allowed Holroyd et al. (2014) to estimate that 18 of 24 nonmarine turtle lineages managed to cross the K/Pg boundary (only 25% went extinct). Using phylogenetic diversity and ghost lineages, both Lyson and Joyce (2009) and Lyson et al. (2011) found that clades of baenid turtles from present-day North Dakota survived the K/Pg boundary extinction. In other local studies, Augustin et al. (2021) found differential survival in dortokid species after the end of the Cretaceous in present-day Romania. Cleary et al. (2020), with a detailed and state-of-the-art occurrence-based approach, found a clear increase in shareholder quorum subsampling (SQS)-corrected richness (at the generic level) across the K/Pg boundary, definitely increasing for North America and relatively stable for South America. But when diversity is counted at the species level, including estimations based on the phylogenetic diversity and ghost lineages, a significant drop in the diversity of turtles is observed, at least in South America (Vlachos et al. 2018). Recently, Pereira et al. (2024) found a phase of significant extinction in turtle species around the K/Pg boundary.

During the last decade, joint efforts of the turtle research community led by W. G. Joyce have thoroughly revised most of the testudinatan clades, allowing the availability of an unprecedented amount of updated taxonomic, anatomical, and distribution information on turtles (Joyce 2014, 2016, 2017; Lawver and Jackson 2014; Cadena and Joyce 2015; Joyce and Lyson 2015; Sterli 2015; Vitek and Joyce 2015; Joyce and Bourque 2016; Maniel and de la Fuente 2016; Anquetin et al. 2017; Georgalis and Joyce 2017; de la Fuente et al. 2018; Vlachos 2018; Joyce and Anquetin 2019; Georgalis et al. 2021). These works offer updated bibliographic and taxonomic information relating to turtles, with quite conservative subjective taxonomic opinions on the validity of turtle taxa, therefore allowing proper curation and completeness of the available PBDB data on turtles. Once this set of expert chapters is complete, it will be possible to perform an in-depth analysis at the turtle clade/family taxonomic level.

The present study aims to provide a global curve of turtle species diversity from the Norian to the present with a special focus on the K/Pg boundary, corrected against sampling bias and after considering the sample coverage in each stage. I focus on the following questions: Did the K/Pg mass extinction event affect turtle species? If so, was their specific diversity already in decline before this event? When did their diversity recover?

The Quality of the PBDB Data

A major difficulty of any quantitative analysis is the quality and completeness of the data analyzed. Here, I use the PBDB, a major resource of paleobiodiversity data initiated in 1999 and the result of collective and independent work of hundreds of researchers for more than two decades.

Currently, there are nearly 1300 references containing different types of information on turtles registered in the PBDB, spanning from Linnaeus's (1758) Systema Naturae published in 1758 to 25 turtle papers published in 2023 (Fig. 1A; also see Supplementary Information). Although the PBDB deals primarily with fossil occurrences, taxonomic information from non-paleontological papers is also recorded and helps greatly with classification, especially in diversity counts. During the first decade of the existence of the PBDB, the registration of turtle-related papers was inconsistent (Fig. 1B,C). Whereas at least 190 turtle papers were published from 1999 to 2009, only 234 turtle papers were registered during the same period in the PBDB, with a delay between publication date and inclusion within the database. For example, the average time for a paper published from 1999 to 2009 to be registered in the PBDB was 8.5 years. During the second decade of the PBDB, and especially during the last 5 years, combined efforts by several researchers have greatly improved the quality of the turtle information in the PBDB. Whereas 650 papers were added to the PBDB up to 2016 (18 years of PBDB data), another 621 papers have been registered since 2017, thus doubling the available information. The average delay of the time of inclusion of a paper in the PBDB has been dropping significantly for recent papers, which are generally included in the PBDB in the same year of their official publication (Fig. 1C).

Figure 1. The amount of turtle-related data from the publication of Systema Naturae until the end of 2023 and their inclusion in the Paleobiology Database (PBDB). A, The cumulative amount of different types of turtle information published across time (named genera, named species, references, collections, occurrences, taxonomic opinions). In the background and with a separate axis at the right, the number of recorded references per year is shown with dark gray vertical columns. B, Violin plot showing the difference between the year of publication and the year of the inclusion of a paper in the PBDB (considering only papers published after 1999, the starting year of the PBDB). More than half of the papers have been included in the database during the same year or the year after of publication. C, Linear regression of the time difference between year of publication and the inclusion of a paper in the PBDB across time. Each year, papers are included closer to their publication date at a faster rate.

In conclusion, I demonstrate that the published information on fossil turtles in the PBDB has increased dramatically in recent years, including hundreds of historical papers, while at the same it is kept up to date with recently published information.

Material and Methods

Data

All data were downloaded from the PBDB on October 23, 2023, for all regular taxa of the search term “Testudinata” at the species level (see “Data Availability Statement”), including both marine and nonmarine taxa as well as taxa with open nomenclature. In particular, the dataset used for this analysis consists of 4142 occurrences identified at the species level, coming from 2863 different collection points, with 2121 accepted valid species, based on information from 1356 different references.

Analyses

Two analyses were performed. The estimation of the sample coverage per time bin selected (in this case, per stage) was performed by entering the data into a spreadsheet and calculating the modified Good's coverage estimator from Chao and Jost (2012: equation 4a and references therein) that takes into account both singleton (i.e., species that are represented by exactly one occurrence in the time bin) and doubleton (i.e., species that are represented by exactly two occurrences in the time bin) occurrences. By “coverage” and sensu Chao and Jost (2012), I mean the sample coverage rarefaction, which estimates the extent to which sampling has covered the expected number of species. It makes sense to compare diversity counts in samples of similar levels of coverage. Diversity per time bin and at the species level was calculated with the SQS correction by using the recently published divDyn package (Kocsis et al. 2019, 2022) in R (R Core Team 2024) that calculates the SQS as in Alroy (2014) and default way in his Perl script also in R; please see Close et al. (2018) for further details. Default SQS uses the so-called three-timer correction (Alroy 2008; CSIB in the divDyn package [Koscis et al. 2022]) that counts taxa sampled in three consecutive time bins. However, it is not reasonable to use this correction at the species level and for stage-level time bins, as it is unrealistic to expect many species of turtles to be present in three consecutive time bins. Besides some species that could exceptionally be expected to extend across three stages and several millions of years, in paleontology, species with such extensive temporal distribution could also be those wastebasket taxa that represent species complexes that cannot be sufficiently distinguished. For example, from the 1156 species counted in this dataset, only 70 (6%) are present in three time bins, and 214 (18.5%) are present in more than one time bin. This analysis aims to capture extinctions especially in the time bins across the K/Pg boundary, and it would therefore be important to count all species present in a single time bin. Those taxa counted with the three-timer correction are few and would be exactly those taxa that have managed to survive across two or more stages or several millions of years. For these reasons, and because this analysis is made at the species level, the sampled-in-bin richness is used (divSIB in the divDyn package; Koscis et al. 2022), of course, SQS-corrected for the 0.8, 0.6, and 0.4 quorum levels. This diversity estimator simply counts all species sampled in each bin. The calculated proportional rates of origination and extinction for the raw diversity (extProp and oriProp) have been used to calculate the turnover between origination and extinction rates for each time bin.

Results and Discussion

Databases, Sampling, and Taxonomy

The latest global diversity analysis for turtles is that of Cleary et al. (2020), who found that SQS-corrected global turtle diversity based on counts of genera was relatively unaffected by the K/Pg mass extinction event. Comparing their results with the previous work of Nicholson et al. (2015; note that they do not address the K/Pg boundary, but they make an extensive analysis of the PBDB at the genus level), Cleary et al. (2020) identify two main reasons for differences: (1) the major taxonomic and stratigraphic revisions since 2015 and (2) a more complete and up-to-date dataset in the PBDB. Both reasons remain true and are taken into account herein. The global fossil record of turtles is still under revision. Still, most clades now have modern revisions with conservative estimations of alpha diversity as well as up-to-date higher-level taxonomy and relationships supported by phylogenetic analyses in most cases. At the same time, the turtle information in the PBDB is now two times more complete and up-to-date than in 2015, and another 120 new studies have been published since 2020.

However, another issue that affects counting comes from the choice of counting genera instead of species. Marshall (2022 and references therein) summarizes this point, showing how different extinction estimations are calculated by counting genera or species. Genera are artificial taxonomic human constructs (for further information, see Hendricks et al. 2014; Wiese et al. 2016 and references therein). The naming of genera and the inclusion of species vary per taxonomic group, the taxa being investigated (e.g., extant or extinct), the geological period in which the fossil record is found, and the subjective criteria of individual researchers and temporal trends in literature. An easy way to illustrate this issue is by plotting families according to the number of species and genera they contain (Fig. 2). In turtles, the ratio between genera and species is about 2:5 (0.42); only half of the named turtle genera are monospecific (349 out of 642 named genera, 54%). The ratio is even more skewed in diverse families that span the Cenozoic and contain extant species (e.g., Testudinidae, Geoemydidae, Emydidae). Therefore, counting genera in turtles would deflate counts. In that sense, the results presented herein are not directly comparable with previous results such as those in Nicholson et al. (2015) and Cleary et al. (2020), which both use generic counts. It is important to point out that subjective taxonomic opinions could affect also specific counts, based on different practices followed by researchers (i.e., lumping or splitting).

Figure 2. Comparison of total number of testudinatan genera and species by family in the current dataset. Families with only monospecific genera are plotted along the diagonal, whereas families with many genera that contain multiple species are plotted away from the diagonal. Turtle families, especially those with extant representatives, contain mostly genera with numerous species. This difference would affect generic and specific diversity counts.

Sample Coverage

There are 3827 occurrences of turtle species that can be safely included in the stage time bins of the analysis (Table 1). The time span represented by the Late Cretaceous and the Paleocene (= Cenomanian through Thanetian) contains 1247 turtle species occurrences. In particular, the time span represented by the Campanian through the Danian contains 1037 occurrences of valid turtle species. Therefore, one-third of the global fossil turtle occurrences are found across the K/Pg boundary. The Campanian and the Maastrichtian represent the stages with the highest sampling coverage, with 88% and 87%, respectively, followed by similarly elevated coverage in the Danian (82%). As such, the available material across the K/Pg boundary represents a good sample for estimating diversity. Given their similar levels of coverage, even the raw counts in these three stages would represent trustworthy species diversity values for ratio comparisons.

Table 1. Sample coverage of turtle fossil data per stage.

Other stages match or exceed the 80% coverage benchmark, albeit with varied numbers of occurrences: Late Pleistocene, Piacenzian, Lutetian, Ypresian, Selandian, Berriasian, and Kimmeridgian. The Toarcian and the Rhaetian show 100% coverage, but this is calculated based on minimal occurrences and should be ignored. Most of the remaining stages show very good sample coverage within the 60–80% range; only a few are poorly sampled, with coverage below 50% and 40%. The sample in most of the time bins shows high values of evenness, with the exception of a few poorly sampled stages.

Overall, the sample coverage analysis suggests that the available information at the species level would allow the calculation of an extensive diversity curve with bias correction at 60% of the sample. In contrast, around the K/Pg boundary, the correction could be elevated at 80% of the sample coverage (or quorum levels of 0.6 and 0.8, respectively).

Species Diversity

Based on the curve of the species diversity (Fig. 3), the highest diversity in the evolutionary history of turtles is observed in the Campanian (extant diversity excluded). During the entire Late Cretaceous, there has been positive turnover with comparatively higher origination than extinction rates. In both the Maastrichtian and the Danian, diversity at the 0.8 quorum is reduced by nearly one-quarter compared with the previous time bin, meaning that turtle species diversity in the Danian was half what it was during the Campanian. During the Maastrichtian–Danian, the rate turnover is always negative, with comparatively higher proportional extinction rates than origination rates. Based on these results, it seems clear that turtle diversity was affected by the K/Pg extinction event, but more importantly, global turtle species diversity was already in decline toward the end of the Cretaceous. By the Ypresian, a recovery of the species diversity to the Maastrichtian levels is observed, followed by a series of fluctuations during the rest of the Paleogene. In particular, comparatively higher proportional extinction rates are observed in the Lutetian, Priabonian, and Rupelian, whereas comparatively higher origination rates are observed in the Bartonian and the Chattian. As a result, in some stages (e.g., the Bartonian), the species diversity of turtles perhaps reached even lower levels than in the post-K/Pg times. The recovery toward the species diversity of the Late Cretaceous levels was achieved in the Neogene thanks to comparatively higher proportional origination rates. In the Pliocene–Pleistocene, comparatively higher proportional extinction rates are observed, but the overall counts are not reduced, probably because of the extant taxa (i.e., Pull of the Recent).

Figure 3. Global turtle species diversity across time, both raw counts (comparable only for stages with similar coverage) and shareholder quorum subsampling (SQS)-corrected counts for 0.8, 0.6, and 0.4 quora (or coverage). Generic raw and SQS-corrected counts for 0.6 quorum are added in gray, for comparison. The turnover between the calculated proportional origination and extinction rates is shown for each stage. In stages with negative turnover, extinction rates are higher than origination rates.

Before the Campanian, the rest of the stages of the Cretaceous are generally subsampled below 60%. Still, the 0.4 quorum curve indicates that the Albian global diversity of turtles could have been comparable to that of the Campanian. Then, another drop in diversity is also observed at the 0.4 and 0.6 quorum levels at the Jurassic/Cretaceous boundary, accompanied by higher extinction rates in the Tithonian. The calculated diversity is too low to draw general conclusions during the Late Triassic and the Early–Middle Jurassic.

As previous studies worked with generic counts, a comparison between specific and generic counts (shown in gray in Fig. 3) is made, based on the updated dataset. The curves of the raw diversity of species and genera are in good accordance during most of the time, showing differences mostly during the Neogene. The SQS-corrected curves are quite similar as well, and they also reconstruct the Maastrichtian–Danian drop, but not the Campanian–Maastrichtian drop seen in the specific curve. The generic curve is better sampled in the Early Cretaceous compared with the specific curve. A major diversity that is only recorded by the SQS-corrected specific diversity curve is the sharp increase during the Miocene, something that raw curves or SQS-corrected generic curves fail to capture. Therefore, the updated generic curve can also represent the majority of the overall changes in diversity across time but tends to flatten out peaks and drops in diversity that are otherwise seen as more marked in the specific diversity curves, especially in the Cenozoic. The currently available information and the progress made during the last decades in fossil turtle studies allows the construction of relatively complete specific diversity curves.

Conclusions

During the last two decades, a huge amount of new data on turtle fossil records has been published, accompanied by joint efforts to curate them with previously published information in open databases like the PBDB. This dataset now allows the construction of global diversity curves at the species level, which indicate that global turtle diversity was already in decline during the last 20 Myr of the Late Cretaceous, before the extinction event, and continued in decline after the extinction event. Compared with the historical peak of turtle diversity in the Campanian, the global turtle diversity was reduced by approximately half during the Danian, only to recover during the Cenozoic. These results differ from previous analyses not only because they are based on an updated dataset, but also because previous estimations are based on genus-level analyses, shown herein to affect generic counts in vertebrate groups like turtles that contain multispecific genera. Future steps should focus on more detailed, local studies that will test this historical and long-decline hypothesis to find whether differential trends can be observed at a continental and/or clade level.