TRACING THE PATTERNS OF NON-MARINE TURTLE RICHNESS FROM THE TRIASSIC TO THE PALAEOGENE: FROM ORIGIN TO GLOBAL SPREAD

: Turtles are key components of modern vertebrate faunas and their diversity and distributions are likely to be affected by anthropogenic climate change. However, there is limited baseline data on turtle taxonomic richness through time or assessment of their past responses to global environmental change. We used the extensive Triassic – Palaeogene (252 – 223 Ma) fossil record of terrestrial and freshwater turtles to investigate diversity patterns, ﬁnding substantial variation in richness through time and between continents. Globally, turtle richness was low from their Triassic origin until the Late Jurassic. There is strong evidence for high richness in the earliest Cretaceous of Europe, becoming especially high following the Cretaceous Thermal Maximum and declining in all continents by the end-Cre-taceous. At the K – Pg boundary, South American richness levels changed little while North American richness increased, becoming very high during the earliest Palaeogene (Danian). Informative data are lacking elsewhere for this time period. However, the Selandian – Thanetian interval, approximately 5 myr after the K – Pg mass extinction, shows low turtle richness in Asia, Europe and South America, suggesting that the occurrence of exceptional turtle richness in the post-extinction Paleocene fauna of North America is not globally representative. Richness decreased over the Eocene – Oligocene boundary in North America but increased to its greatest known level for Europe, implying very different responses to dramatic climatic shifts. Time series regressions suggest number of formations sampled and palaeotemperature are the primary inﬂuencers of face-value richness counts, but additional factors not tested here may also be involved. analyses suggest that turtle richness dynamics are complex and do not correspond strongly or simply to common drivers of richness, such as palaeotemperature and non-marine area. The number of turtle-bearing formations was expected to rank highly in the model comparisons, as a large number of formations producing turtles should be more likely to produce specimens that can be identiﬁed to genus level. When this is removed, the best models are the tetrapod-bearing collections (TBCs) only and tetrapod-bearing formations (TBFs) only models, though the R 2 values are very low (see Cleary et al. 2020). Given that all of the best models selected here included a sampling proxy,

T E S T U D I N A T A , consisting of terrapins, tortoises and turtles (collectively referred to as 'turtles' herein), spans a time period from the Late Triassic to the modern day (Joyce 2017). Today, it is a geographically widespread clade with approximately 90 genera containing 350 species (Uetz et al. 2018) and is divided into two major lineages: Cryptodira (hidden-neck turtles) and Pleurodira (side-neck turtles). Extant pleurodires are currently found only in the southern hemisphere, but ranged into the northern hemisphere in the past, and are restricted to aquatic or semi-aquatic lifestyles. Cryptodires, on the other hand, are common in both hemispheres and occupy a wider range of aquatic and terrestrial niches, including the marine realm (Pough et al. 2013). Both extant lineages extend back into the Jurassic, alongside various stem-group taxa that diverged prior to the cryptodire/ pleurodire split (Joyce 2017) and have a rich and widespread fossil record. Turtle origins are heavily debated, with multiple theories on their origin and closest relatives (see Joyce 2015; Schoch & Sues 2019).
Turtles have been abundant and widespread components of non-marine vertebrate faunas since the Middle Jurassic and their ectothermic physiology and broad geographical distribution make them a key clade for studying the effects of both deep-time and anthropogenic climate change on vertebrates (Ara ujo et al. 2006;Foufopoulos et al. 2011;Urban et al. 2014;Waterson et al. 2016). For example, turtle sex ratios are determined by temperature and are expected to be heavily influenced by future climate warming (Bickford et al. 2010). The ability of individual taxa to cope with rapid climate change varies, but ectothermic vertebrates appear to have relatively poor plasticity with regard to thermal tolerances (Gunderson & Stillman 2015).
Several regional studies have assessed face-value turtle richness at key moments in time, often focusing on major extinction events or climatic changes (Hutchison & Archibald 1986;Hutchison 1992;Hirayama et al. 2000;Holroyd et al. 2001Holroyd et al. , 2014Corsini et al. 2011). However, besides inclusion in whole-tetrapod analyses of richness (Kalmar & Currie 2010;Benson et al. 2013Benson et al. , 2016Close et al. 2017Close et al. , 2019, only three studies have examined the turtle fossil record in detail using methods that account for or alleviate geological and anthropogenic sampling biases, all of which deal with Mesozoic non-marine turtles (Nicholson et al. 2015(Nicholson et al. , 2016Vlachos et al. 2018).
Here, we investigate non-marine turtle richness through time using an expanded, updated dataset for non-marine turtle distribution and richness that extends from the Late Triassic to the Palaeogene (237-223 Ma). This expanded interval allows more detailed comparisons with long-term datasets that have been generated for other non-marine tetrapod clades (Newham et al. 2014;Mannion et al. 2015;Cleary et al. 2018). Similar patterns might identify common causes, whereas notable differences could indicate clade-specific responses to long term environmental change. We generate sampling-corrected taxon-richness curves and discuss some of the driving factors that might have influenced these richness patterns, using these data to assess problems with both the turtle fossil record and the fossil record as a whole.

Data
Occurrences of Testudinata from the Triassic-Palaeogene were downloaded from the Paleobiology Database (PBDB; https://www.paleobiodb.org) on 13 January 2020. Our dataset builds on that presented by Nicholson et al. (2015Nicholson et al. ( , 2016, with further additions and corrections to taxonomic records by PAH and other contributors to the PBDB; further details of individual contributions are listed in the acknowledgements. After downloading, these data were filtered to remove any remaining marine and ichnofossil occurrences as well as those not identifiable to generic level. Additionally, we ensured dating was consistent between occurrences from the same time intervals (see Cleary et al. (2020) for the script covering both processes). We followed the definition of Testudinata used by Joyce (2017), as a clade originating from the first amniote to produce a fully developed turtle shell, and thus excluded Odontochelys and Pappochelys from this study. The resulting dataset consists of 2989 occurrences representing 308 genera.
We recognize that there are problems with using genera rather than species, particularly with 'wastebasket taxa' that have been used widely in the turtle fossil record (e.g. 'Trionyx'). Wastebasket genera are, at worst, a conservative underestimation of potential turtle richness in each time bin (indicating at least one genus of that group, but possibly representing more). Furthermore, species taxonomy is poorly resolved in many clades or species cannot be readily distinguished from the available fossils. However, genera can be assigned with more confidence in many assemblages (Holroyd & Hutchison 2002). For the latter reason, we included genus-level occurrences that are specifically indeterminate (>1000 occurrences). Occurrences that can be identified to higher taxonomic levels only (e.g. family or clade) were excluded because they do not reliably represent a single taxon or a comparable rank.

Subsampling
We used shareholder quorum subsampling (SQS; Alroy 2010a, b) to produce subsampled genus richness estimate curves for non-marine turtles at global and continental levels. SQS uses taxon occurrence frequencies, rather than absolute counts of occurrences, to determine when the subsampling quota (= quorum) has been met. These frequencies are the proportion of total occurrences that each taxon accounts for in a time bin. Occurrences are sampled until a desired summed frequency quorum is met, with each taxon that is drawn contributing its frequency only once per bin. This would result in fair and comparable sampling of the abundance-frequency distribution only if all intervals/regions were equally sampled to begin with. However, the observed abundance-frequency distribution represents only a subset of the true distribution. Therefore, Good's u, the proportion of single-occurrence taxa ('singletons'), is used to estimate the coverage of a time bin's dataset (i.e. how much of the 'true' underlying richness has been recovered). By dividing a chosen quorum by each time bin's value of u, one can estimate the 'true' coverage of each bin and subsample more intensively when coverage is low (for further information see Alroy 2010a). Consequently, if a time bin's estimated coverage is lower than the chosen quorum level, then that bin will be excluded from the analysis as the quorum cannot be met.
SQS was performed using the v.4.3 Perl script (see Cleary et al. 2020) on global and continental-level turtle occurrences at quorum levels of 0.4-0.7. Quorum level 0.4 is considered to be the lowest level at which a sufficient level of standing diversity can be recovered (Alroy 2010b;Benson et al. 2016). Some standard stratigraphic stages were merged to reduce variance in durations between bins (Table 1), resulting in interval durations of approximately 9 myr. Occurrences were only included in the analyses if their stratigraphic ages were known with sufficient certainty to be assigned to a single time bin. Exclusions for each continental region ranged from fewer than 12% of occurrences (Indo-Madagascar, North America, Europe) to very high due to high stratigraphic uncertainty (Asia, 34%; Oceania, 64%, but only 11 collections in total). To account for the broad or poorly known dating of some Asian localities, we also included individual points for the Jehol Group (plotted between bins K2 and K3) and Nemegt Formation (plotted between K7 and K8) biotas.
In order to avoid problems with small sample sizes when using Good's u, points generated from fewer than five collections were excluded from plots, except for one figure with these points included for reference (see Cleary et al. 2020). All SQS curves are annotated with the numbers of collections that contribute to each point in order to display underlying data quality.
We also used classical rarefaction (CR), which uses absolute counts of occurrences rather than coverage when subsampling. Using absolute counts with a set quota ('uniform sampling') has been criticized, as this can dampen genuine diversity signals (Alroy 2010b). Nevertheless, we use this method as a comparison to SQS to see if results vary due to the application of different subsampling techniques. Results for these analyses can be found in Cleary et al. (2020).
As a visual aid to our discussion of geographical biases, we plotted the palaeolatitudinal distributions of every occurrence used in our analyses by major clade (Cryptodira, Pleurodira and Paracryptodira). Palaeolatitudes were calculated by the PBDB using tectonic plate rotation data from the GPlates model (Wright et al. 2013).

Comparing richness to proxies for sampling and environmental factors
To examine relationships between face-value occurrence data and proxies for sampling and environmental variables, we used generalized least-squares regressions (GLS). Sampling proxies considered herein are the numbers of non-marine tetrapod-bearing collections (TBCs), tetrapod-bearing formations (TBFs) and turtle-bearing formations (TuBFs), which were collated from PBDB tetrapod data (59 823 occurrences downloaded on 20 January 2020; see Cleary et al. 2020). Environmental proxies consist of sea-level (Miller et al. 2005), non-marine area (Smith et al. 2004) and palaeotemperature (d 18 O record;Prokoph et al. 2008vide Mannion et al. 2015. We interpolated palaeotemperature data into 0.1 myr intervals to facilitate subsequent aggregation into time bin averages (see Cleary et al. 2020). These data were only available from the Jurassic onwards, so our GLS analyses include bins J4 to Pg5 (Table 1). We added time bin duration as a non-optional variable to determine if it was influencing our model outcomes, following Marx & Uhen (2010).
For GLS we fitted autoregressive (AR) models of order one to all possible combinations of the above proxies. Each model included only one of our sampling proxies, in order to limit model complexity. We compared the models using Akaike's information criterion for small sample sizes (AICc). We calculated Akaike weights to identify the best model from those tested; AIC weights reward goodness of fit but penalize a higher number of variables in a model (Burnham & Anderson 2002). We manually calculated R 2 values using the method of Nagelkerke (1991), which provides an estimate of how much variance in face-value counts of richness is explained by each model. Sampling proxies and face-value richness were ln-transformed prior to analyses to ensure homoskedasticity and normality of the models' residuals; T A B L E 1 . Time bins of roughly 9 myr each used in richness and GLS analyses, and the standard stratigraphic intervals they correspond to.

RESULTS
Global face-value and subsampled curves Figure 1 shows face-value, counted genera richness compared against subsampled global richness using SQS at a quorum of 0.4. In the face-value curve, apparent richness increased steadily from a very low level in the Late Triassic following the first appearance of the clade until the mid-Cretaceous, after which it sharply rose until the end of the Mesozoic. There was a sharp decline over the Cretaceous-Palaeogene (K-Pg) boundary followed by a sequential pattern of peaks and troughs in the Palaeogene. Palaeogene richness remained high compared to pre-Late Cretaceous levels. The subsampled curve has a number of marked differences, including a peak in richness in the late Early Cretaceous (K4; Albian). The latest Cretaceous (Campanian-Maastrichtian; bins K7-8) is flatter, lacking the strong peak present in the face-value data, and relative richness is roughly the same or lower than in the late Early Cretaceous. By contrast with the face-value curve there was a slight increase in richness across the K-Pg boundary, and while the peaks and troughs remain in the Palaeogene their scale is uneven compared to the facevalue curve. The quora of bins J1-J4 and K3 are lower than the chosen quorum (0.4) and do not appear in the figure. At higher quorum levels, the latest Cretaceous (K7-K8) and late Eocene (Pg4) had the highest peaks of apparent global richness (though some bins drop out due to poor sampling; see Cleary et al. 2020, fig. 1). Classical rarefaction results are very similar to the SQS results but the signal is dampened, and hence we focus our discussion on the results of SQS (for CR see Cleary et al. 2020, figs 2-3).

Continental-level subsampling
Subsampled richness for individual continents at quora of 0.4-0.7 (Figs 2-3; Cleary et al. 2020, figs 5-6) reveals substantial variation in richness among regions, which indicates that the global level curves mask marked geographical differences in pattern. Due to poor sampling, robust estimates of regional richness cannot be recovered prior to the Late Jurassic even at a quorum of 0.4. The only exception to this is the confident recovery of low richness in the Late Triassic of Europe. Former Gondwanan continents are also poorly sampled. Of these, subsampled richness estimates can only be recovered for South America for a few time intervals and Indo-Madagascar for one. To demonstrate the differences between genuinely poor sampling and sampling intensively in one place, a face-value continental curve is included for reference in Cleary et al. (2020, fig. 7).
In the Late Jurassic, we recovered an intermediate level of richness for Asia in K5 (Callovian-Oxfordian) and low North American and high European richness in K6 (Kimmeridgian-Tithonian; Fig. 2). When including points with fewer than five collections (Cleary et al. 2020, fig. 4) we found that J6 of Asia had a low level of richness comparable to North America.
Over the J-K boundary and into the Early Cretaceous, richness decreased in Europe though it remained at relatively intermediate levels. Sampling is poor in the middle-late Early Cretaceous for most continents. The occurrences from the Jehol Group that fit into a single time bin contribute a large portion of the richness in Asia during the late Early Cretaceous, and subsampling this biota alone suggests that it was characterized by intermediate richness. Early Cretaceous Asia has a low number of localities contributing data, however. Asia in bin K4 (Albian) had a slightly higher level of richness to bin K2. Bin K3 is too poorly sampled to assess any sequential patterns between the bins surrounding it. North America had an increase in richness from bins K3 to K4, but richness remained low. South America in bin K4 also had low richness levels (although the estimate is based on data from only eight localities). Overall, these findings suggest substantial regional variation in the genus richness of turtles during the Early Cretaceous.
Late Cretaceous richness patterns can be reconstructed for North America, Asia, Europe, South America and Indo-Madagascar. At a quorum of 0.4 we recover a nearcomplete richness curve for the entire Late Cretaceous of Asia, missing bin K7 (Campanian). Richness increased from bins K4 to K5 (Albian to Cenomanian) then decreased to lower levels in bin K6 (Turonian-Santonian). After a gap due to poor sampling, Asian richness in the Maastrichtian (K8) was at intermediate-to-high levels.
North America had intermediate to high richness throughout the latest Cretaceous (K6-K8), with a higher richness in bin K7 compared to K6 and K8. In the Maastrichtian, richness levels are roughly the same as Asia. In Europe, sampling is too poor to estimate richness during K6, but in K7-K8 richness was initially intermediate and decreased into K8. The Maastrichtian (K8) of South America and the well-sampled separate Nemegt Formation of Asia (dated as Campanian/Maastrichtian) exhibit low richness, similar to that in Europe. Indo-Madagascar had an intermediate level of K8 richness that was slightly higher than in South America, Asia or Europe.
Subsampled richness increased among North American turtles across the K-Pg boundary. Richness appears to have neither increased nor decreased in South America but cannot be estimated at a quorum greater than 0.5 (Cleary et al. 2020, figs 5-6). Subsampled North American turtle richness was at an all-time high in Pg0 (Danian). The richness of other regions cannot be estimated during this time. Richness declined from the early to late Paleocene (Pg0 to Pg1, Danian-Thanetian) in both North America and South America, though richness remained relatively high in the North America compared to South America. Richness in Pg1 of Asia and Europe was low, approximately similar to that in South America, and was comparable to their lowest levels in the Mesozoic. This provides a striking contrast between Pg0-1 of North America and the other continents.
A return to high richness levels occurred in Pg2 (Ypresian) in Europe and Asia, while North America decreased slightly compared to Pg1. These continents had roughly comparable richness at this time, with Europe's slightly higher. Richness cannot be estimated for other continents. F I G . 1 . Non-marine turtle global genera richness from Triassic-Palaeogene using: A, uncorrected (face-value) richness; B, subsampled (SQS) richness at a quorum of 0.4. Numbers indicate collections drawn for each time bin subsampled, as an indicator of underlying data quality. Unfilled point is contributed to by < 10 collections.
North America maintained roughly the same level of high richness until bin Pg4 (Bartonian-Priabonian), in which richness peaked to a very high level. Europe, in comparison, showed a pattern of alternating increase and decrease akin to the global data in bins Pg2-4. By Pg4, Europe had a high richness that was slightly lower than in Pg2, but lower than the very high North American richness in Pg4. Asia is recovered again in this bin, with a similar richness to Europe. At the Eocene-Oligocene boundary (E-O; Pg4-5) there was a fall in richness for North America, whereas in Europe there was a sharp rise to a very high level. No other continents are recovered for this time period, except a very low richness in South America (only three collections; Cleary et al. 2020, fig. 4), as they are too poorly sampled.
Patterns obtained at a quorum of 0.5 ( Fig. 3) are broadly similar to those at 0.4, but with several key differences. First, one bin is too poorly sampled to yield richness estimates at a quorum of 0.5, K3 (Aptian) of North America, but this does not affect the overall curves very much.
Fewer regions/intervals provide richness estimates at a quorum of 0.6 (Cleary et al. 2020, fig. 5), for which the following bins do not provide sufficiently informative data: the Jehol biota of Asia (K2/3), K4 and K8 of South America and K4, K6 and Pg3 of North America. This removes South America's pre-Cretaceous record and the K-Pg boundary pattern meaning that very few meaningful results were obtained from the southern hemisphere at this quorum level. At a quorum of 0.7 (Cleary et al. 2020, fig. 6) many more bins drop out, including most of Asia (except K4, K6 and Pg1), J6, Pg1 and Pg5 of Europe, Pg4 of North America, K8 of Indo-Madagascar and Pg0 of South America.
We also examined differences in the richness of the major testudinate clades (Cryptodira, Pleurodira and Paracryptodira) through time . These clades were assigned to taxa manually (by PAH and TJC; see script in Cleary et al. 2020). Pleurodira, which is less species-rich and confined to freshwater in the southern hemisphere today, is predictably more poorly sampled (Fig. 5A), while cryptodires contribute most of the richness record for non-marine turtles (Fig. 4). The majority of paracryptodires are known from North America, and thus are well sampled only on that continent (Fig. 5B). The Palaeogene part of the North American genus richness curve and the entire Asian curve appear to have been driven almost entirely by cryptodires with minimal influence from paracryptodires, while the South American and Indo-Madagascan data are dominated by pleurodires. The European record had a mixture of both cryptodires and pleurodires, with more input from the former.

Generalized least-squares analyses
Of the combinations of explanatory variables that we tested with generalized least-squares regression (GLS ;  Table 2), eight models were found to be a better fit than the null model. The best four of these were turtle-bearing formations (TuBFs) only (AICc = 9.27), followed by single combinations of turtle-bearing formations with d 18 O (AICc = 13.47), non-marine area (NMA; AICc = 19.47) and sea level (AICc = 21.62), respectively. Tetrapod-bearing collections and formations (TBCs and TBFs) trailed just behind these. From these results, the number of turtle-bearing formations sampled appears to have had the greatest influence on generic richness through time. Bin duration was not a significant variable in any of the best models, and the differences in the ordering of the best models after the top two is minimal, so the analyses using bin length as an extra compulsory variable are presented in Cleary et al. (2020).
While the TuBFs-only model is recovered as the best model, with a significant p value (p < 0.01), the model's R 2 value is not especially high (0.67). This is also true of the second and third-best models (TuBFs + d 18 O, R 2 = 0.65; TuBFs + NMA, R 2 = 0.50). Consequently, there are likely to be other variables that were not tested in these analyses that contribute to turtle generic diversity patterns.

DISCUSSION
Turtle richness and the effects of subsampling and regional analysis Global face-value genus counts show steadily increasing richness from the first appearance of turtles up to the 'middle' Cretaceous, then a more abrupt increase to the Maastrichtian (K8). This was followed by a decrease over the K-Pg boundary and a series of fluctuations throughout the Palaeogene, with peaks in richness during Pg2 (Ypresian), Pg4 (Bartonian-Priabonian) and Pg5 (Oligocene; Fig. 1). However, many features of the global facevalue curve changed drastically when subsampling was applied, particularly during the Mesozoic, and especially when examining regional datasets. These results illustrate the pitfalls of using 'global' diversity patterns extracted uncritically from the fossil record, and we caution strongly against using global diversity estimates unless the underlying regional structure has also been thoroughly explored. Many supposedly 'global' patterns, both for face-value and subsampled richness, are in fact caused by F I G . 5 . Subsampled genera richness at a quorum of 0.4: A, Pleurodira; B, Paracryptodira. AF, Africa; AS, Asia; EU, Europe; IM, Indo-Madagascar; NAm, North America; SA, South America. Numbers indicate number of collections drawn for each time bin subsampled, as an indicator of underlying data quality. Unfilled points are contributed to by <10 collections. Colour online. richness patterns present in only one or two samplingdominant regions plus a combination of other areas that are not individually well-sampled enough to provide subsampled richness estimates at a quorum of 0.4. Indeed, a strong relationship between subsampled richness and the geographical spread of fossil localities has been demonstrated for tetrapods (Barnosky et al. 2005;Benson et al. 2016;Close et al. 2017). Because of these issues, our interpretations will focus on regional subsampled diversity estimates.
We were not able to characterize richness patterns in many regions during the Early-Middle Jurassic and Early Cretaceous (Figs 1B, 2). The turtle fossil record is poorly sampled during these intervals (Nicholson et al. 2015) and further field exploration will be required to at least partially address this imbalance. We recovered high richness in K4 (Albian) and K7-8 (Campanian-Maastrichtian) from a combined Asia and North America, with some input from other continents, particularly in the latter bins. The transition across the K-Pg boundary is primarily observable in North America (with some data from South America), meaning that the 'global' curve for this time is driven almost entirely by the North American record. Whilst the global curve shows fluctuating increases and decreases in richness for the Palaeogene from bins Pg0-4, the underlying, northern hemisphere, continents have differing patterns of richness that are not captured by an amalgamated 'global' curve.

Mesozoic richness patterns and major clade diversifications
Triassic. Turtle diversity during the Triassic was low compared to many later intervals in all the analyses. The fossil record for this time is sparse and consists of stem turtles only, which had their first definite appearances in the Norian (Tr4; Joyce 2017). We recovered a European Rhaetian (Tr5) point at a quorum of 0.4 with slightly lower richness than Tr4 but this is excluded from Figure 2 due to the very low number of collections that contribute to it (two). Similarly low levels of subsampled richness from less than the threshold number of collections were also recovered for Asia in Tr4 (see Cleary et al. 2020, fig. 4, which includes points supported by < 5 collections). The low level of richness in Tr4 persists up to a quorum of 0.7 (see Cleary et al. 2020, fig. 6), and although the European record is the only one robust enough to be recovered in our analyses, by the Norian (Tr4) turtles had already achieved a cosmopolitan distribution (Joyce 2017). This supports the earlier observation that poor sampling is obscuring much about the origins and early history of turtles (Anquetin et al. 2009;Joyce 2017), but also provides evidence for a global spread of early turtles that appears separate from taxonomic diversification.
Jurassic. The Early-Middle Jurassic is very poorly sampled providing insufficient information to estimate richness, even at a quorum of 0.4. This limits what can be understood about the origin of the turtle crown group, which was thought to occur in the Middle-Late Jurassic, based on the first appearance of geographically widespread crown group fossils in the Late Jurassic (Sterli 2008;Anquetin et al. 2009). These observations suggest that turtles could have been widespread and abundant  (Newham et al. 2014). It is unclear whether this sampling gap is due to genuinely low diversity or to a lack of sampling from suitable environments (e.g. arid places unsuitable for freshwater aquatic turtles).
The late Middle and Late Jurassic witnessed various lineage-splitting events that gave rise to six distinct groups of crown and stem turtles. These lineages developed, perhaps through vicariance, in three different geographical areas, with one crown and stem lineage in each. In Euramerica, these crown and stem lineages were Paracryptodira (North America) and Helochelydridae (= Solemydidae; Europe); in Asia they were Cryptodira and Sichuanchelyidae; and Pleurodira and Meiolaniiformes were predominantly Gondwanan in distribution (Joyce et al. 2016; Joyce 2017). However, it should be noted that Paracryptodira might belong to the cryptodire stem (Joyce, 2017), which could modify interpretation of this biogeographical pattern. Due to data deficiency for the Early and Middle Jurassic, our analyses uncover no clear richness patterns until the Late Jurassic (J5-J6), when we find divergent patterns among geographical regions, indicating relatively low richness in North America, intermediate levels in Asia and high richness in Europe. Unfortunately, no meaningful diversity estimates are available from the Late Jurassic of the southern hemisphere, which means that we cannot comment on the origins or early diversification of Pleurodira or compare our results to those of Vlachos et al. (2018), who examined the diversity of South American turtles from the Triassic to present day based on face-value counts combined with ghost lineage diversity.
Cretaceous. The Jurassic-Cretaceous boundary is poorly sampled for all continents except Europe, but here richness decreased, between the Kimmeridgian-Tithonian (J6) and the Berriasian-Valanginian (K1), although the apparent decrease is small and richness remained at intermediate-to-high levels that were comparable to those seen in the Late Cretaceous and Palaeogene. High sampling in Europe at this time is due in part to the very productive Wealden and Purbeck group localities in the UK. During the Early Cretaceous, more derived members of Eucryptodira and Paracryptodira flourished (P erez-Garc ıa 2012; P erez-Garc ıa et al. 2014) and the Xinjiangchelyidae made their first appearance in Europe at this time (probably immigrating from Asia; P erez-Garc ıa et al. 2015). Unfortunately, the record is too poorly sampled to allow comprehensive richness estimates from all continents, preventing us from observing the potential effects of these The separately plotted Jehol Group (between K2 and K3) had comparable, but lower, richness to that of the cumulative Asian record at this time. This is consistent with the fact that this Konservat Lagerst€ atte makes a large contribution to what we know of the turtle record from this time, and highlights potential issues with sampling discrepancies between areas with vastly different preservation potentials. South American richness was at a similar but slightly lower level to that of the Jehol biota during the subsequent Albian bin (K4); here, richness was composed largely of pleurodires (Pelomedusoides) from the Santana Formation of Brazil, as South America is otherwise very poorly sampled in the Early Cretaceous. At this time in North America, new lineages, such as the endemic Baenidae, established themselves (Hirayama et al. 2000; Joyce & Lyson 2015) alongside taxa that persisted from the Late Jurassic (e.g. pleurosternids, solemydids). North American richness increased from bin K3-K4 (Aptian-Albian), but it is difficult to observe what occurred between here and the higher richness in the Late Cretaceous due to poor sampling in K5 (Cenomanian). Overall, the points recovered from the Early Cretaceous that are robust up to a quorum of 0.6 (K1 Europe, K2 and K4 Asia), are in areas with either spectacular preservation or that have been subject to intensive study. These points highlight the continued diversification and persistence of non-marine turtles despite the overall low level of sampling indicated by our results.
In the Late Cretaceous we can estimate richness for more intervals from the northern hemisphere, particularly in Asia and North America. The Cenomanian (K5) suffers from poor sampling apart from a very low richness recovered for Europe, and an increased richness in Asia compared to K4, for which most of our data is derived from the Khodzhakul Formation in Central Asia. This diversity decreased into K6 (Turonian-Santonian) in Asia, while North America had a much higher richness than in earlier bins. The increase in North American richness coincides with the Cretaceous Thermal Maximum (CTM;Forster et al. 2007). The occurrence of turtles at high latitudes in the Turonian (71°N; see Fig. 6) suggests that the higher temperatures allowed for geographical range expansion, as turtles tend to have a conservative latitudinal range compared to other terrestrial tetrapods, such as mammals and most dinosaurs, at other times (Nicholson et al. 2015(Nicholson et al. , 2016. Given this, it is surprising that Asian diversity decreased in bin K6 (the sampled area, mostly Central Asia, remains very similar), although in later Cretaceous intervals it returned to high levels. It is likely, however, that temperature is not the primary environmental factor driving turtle distribution and that other factors accompanying temperature (e.g. precipitation or habitat heterogeneity) may have had differing magnitudes of influence on different continents (Iverson 1992;Holroyd & Hutchison 2002;Waterson et al. 2016).
During K7 (Campanian) there is a gap in the Asian richness curve due to poor sampling, although in the cryptodire-only curve (Fig. 4) (Fig. 2). The Nemegt Formation (Campanian-Maastrichtian; for which richness was estimated separately) is intensively sampled, with 44 collections, but richness is lower than in Maastrichtian Asia, but comparable to K8 in Europe and South America. Richness remained high in the Maastrichtian of Asia compared to other continents, and matched North America, and a diverse fauna including land-adapted clades thrived (e.g. nanhsiungchelyids).
European subsampled richness decreased between the Campanian and Maastrichtian (K7-K8). Campanian richness was lower than that of North America, but still relatively high, as non-marine pleurodires (bothremydids) appeared in the northern hemisphere during the Santonian and quickly became the dominant turtle clade in Late Cretaceous Europe. It has been proposed that this was influenced by the CTM, as modern pleurodires have stricter environmental requirements than cryptodires, and a period of warming may have facilitated their dispersal to other continents via short-range oceanic rafting from Africa (Nicholson et al. 2016; P erez-Garc ıa 2017; Ferreira et al. 2018). Nevertheless, unlike other turtles, pleurodires do not reach particularly high latitudes at this time (Fig. 6). This is consistent with their hypothesized narrower thermal tolerances, although their expansion to other continents might also be explained by other environmental or sampling factors (the mid-Cretaceous European record is very poor).
Richness decreased in the Maastrichtian, which could be attributed to poor sampling, but the countries sampled in this period are broadly similar to those sampled in the Campanian (with the addition of Romania), so it might be that a genuine decrease occurred at this time. It is also possible that there are issues with undiagnostic material or wastebasket taxa. In Romania, for example, the Maastrichtian turtle fauna recovered so far is monogeneric, consisting solely of the meiolaniiform Kallokibotion, but there is evidence of other turtles unidentifiable to the generic level (e.g. occurrences of 'Dortokidae indet.' from the Maastrichtian of Teleac; Rabi et al. 2013).
South America and Indo-Madagascar are sufficiently sampled in the Maastrichtian (up to a quorum of 0.6 for Indo-Madagascar) to estimate subsampled genus richness. South American richness was slightly lower than that in the Maastrichtian of Europe, while Indo-Madagascar richness was at intermediate levels between those of South America, Europe and Asia, with the pleurodire clades Pelomedusoides and Podocnemidoidae dominating all assemblages. Other Gondwanan continents are too poorly sampled to estimate Maastrichtian turtle richness.
North America exhibits high richness during K6-K8, following the CTM. Trionychids first appeared in the North American record during the Santonian (Brinkman et al. 2017), as pleurodires did in Europe, perhaps again due to the warmer conditions facilitating immigration from Asia via high latitude dispersal routes. The terrestrial Nanhsiungchelyidae also appeared here in the Late Cretaceous (Basilemys; Holroyd & Hutchison 2002). Regional differences within North America are observed in the face-value record between the Campanian and Maastrichtian (Brinkman 2003; Holroyd et al. 2014) but the data from the PBDB makes it difficult to corroborate these observations with subsampled curves. The majority of North American Campanian occurrences are from the south-western United States (Utah, New Mexico) plus Alberta, Canada, whilst the Maastrichtian records primarily come from the more northerly United States (Montana, Wyoming, the Dakotas) with some overlap of states between the two stages. This means we are comparing two different geographical areas that may have had very different habitat suitability for turtles. The decrease in Maastrichtian richness may reflect the fact that we are comparing a region closer to the edge of the turtle range (K8) with a region much further south (with some input from the edge in Canada; K7). Regardless, the decrease in richness differs from North American lepidosaurs, which increased from K7-8 (Cleary et al.

Results compared to Nicholson et al. (2015)
The results from the Mesozoic portion of our study differ from those of Nicholson et al. (2015) largely due to changes in data availability; we excluded points recovered from regions/intervals with fewer than five turtle-bearing collections and added richness estimates representing the Early Cretaceous Jehol biota of China and the Late Cretaceous Nemegt Formation of Mongolia. Our richness estimates for regions/intervals based on fewer than five collections (see Cleary et al. 2020, fig. 4) differ, a finding that matches the expectation that richness estimates based on inadequate information are highly sensitive to the addition of new data Among these data-poor regions/intervals, we both recovered the same Triassic richness for Europe and Asia, and the same K7-K8 decrease for Europe, but outside of this the patterns recovered are different or have missing data making comparisons impossible. Some of these differences are due to changes in the stratigraphic ages of various localities since the earlier study. For example, the Khodzhakul Formation of Uzbekistan is now dated as Cenomanian (Averianov & Archibald 2005) whereas in the Nicholson et al. older PBDB data occurrences were dated Albian, Cenomanian or spanned both; F I G . 6 . Palaeolatitudinal distribution of non-marine turtle occurrences from Triassic-Palaeogene; light grey circles indicate other non-marine tetrapod occurrences. the change in allocated time bin resulted in a reduction of richness in bin K4 (Albian) and an increase in K5 (Cenomanian). Also, the ability to provide a subsampled diversity estimate in North American bin K6 (Santonian) reflects differences in dating the Canadian Milk River Formation in the PBDB between their dataset and ours (previously Campanian, now Santonian, following Payenberg et al. 2002).
Palaeogene richness patterns and major clade diversifications K-Pg boundary and Paleocene. The K-Pg mass extinction has been intensely studied for its effects on global richness, most famously for the demise of the non-avian dinosaurs (for a review, see Brusatte et al. 2015). Most previous studies suggested that turtles had very high rates of survival, at least in North America and possibly elsewhere (Hutchison 1982;Hutchison & Archibald 1986;Sheehan & Hansen 1986;Hutchison & Holroyd 2003;Lyson et al. 2011;Holroyd et al. 2014).
In North America, which has the best continuous rock record across the K-Pg boundary, we found a definitive increase in richness between the Maastrichtian (K8) and Danian (Pg0). The earliest richness in the Puercan (= Danian) Tullock Member of the Fort Union Formation of Montana was previously observed at face-value as being just as high as in the preceding Maastrichtian Hell Creek Formation (Holroyd et al. 2014). Minimum survivorship was 75% for the northern Great Plains region, with 18 of 24 genus-level lineages confidently known to have crossed the K-Pg boundary. Holroyd et al. (2014) highlighted that stratigraphic uncertainty means that we cannot tell if certain turtle genera made it to the boundary or were extinct long prior (perhaps even hundreds of thousands of years prior) to the K-Pg event so that true survivorship could be much higher. Nevertheless, a large proportion of North American Cretaceous genera are found in the earliest Paleocene, including many baenid and trionychid genera as well as members of smaller clades such as kinosternoids, pleurosternids and adocids. Our subsampled curves are consistent with high turtle survivorship, and with diversification during the earliest Paleocene in North America leading to increased postextinction richness. This post-Cretaceous increase concurs with results obtained for non-marine tetrapods by Benson et al. (2016) and Close et al. (2017), which were mainly driven by large increases in mammalian richness.
In South America there was essentially no difference between pre-and post-boundary richness and the same podocnemidoid genera occur in bins on either side of the boundary in our data. However, the number of collections contributing to the latest Cretaceous (K8) is low (disallowing estimates of richness above a quorum of 0.5) and the geographical spread of the data differs substantially between the two bins. This is exacerbated by formations whose dates span too long a time period, for example the Adamantina Formation (late Coniacianlate Maastrichtian; Castro et al. 2018) or La Colonia Formation (Campanian-Maastrichtian, possibly also Danian;Gasparini et al. 2015) and could not be included in this analysis. Once the ages of these formations are better constrained it may be possible to observe richness patterns more precisely but based on the evidence currently available, it appears that turtle richness changed little over the K-Pg boundary in South America. The pattern from subsampled data differs from that found in the face-value richness and phylogenetic diversity curves produced by Vlachos et al. (2018). In the face-value curve there is an increase in South American richness, but there is a decrease in phylogenetic diversity (represented by the number of branches on a time-calibrated phylogenetic tree). As neither of these curves were corrected for sampling biases and use different methodologies it is difficult to compare them to our results further.
Asia is poorly sampled across the K-Pg boundary. If we include estimates generated from fewer than five collections (Cleary et al. 2020, fig. 4) then Asia shows a drop. However, we caution that this inference is based upon an insufficiently informative data pool drawn from a small number of collections. Extinctions across the K-Pg boundary included two land turtle clades, the Nanhsiungchelyidae and Sichuanchelyidae (Danilov & Syromyatnikova 2008;Joyce et al. 2016), which were previously common in Asia. It is likely that ecological differences played a large part in turtle extinction selectivity at the boundary and that aquatic turtles had much higher survivorship, as noted by other authors studying turtles and (similarly aquatic) non-marine crocodylians (Hutchison 1982;Hutchison & Archibald 1986;Brochu 2003;Bronzati et al. 2015;Mannion et al. 2015). Increased survivorship of aquatic clades has been posited to be due to differences in trophic interactions compared to those for fully terrestrial clades. After the bolide impact the atmosphere would have filled with dust and debris, obscuring sunlight and leading to a reduction in primary productivity (Alvarez 1983;Kaiho et al. 2016). Modern aquatic habitats, however, rely more on detrital-based food webs rather than in situ primary productivity, and so aquatic environments may have been buffered against the short and long-term environmental changes following this event (Sheehan & Hansen 1986;Robertson et al. 2013;Kaiho et al. 2016). No alternative hypotheses have been proposed so far to explain the selective survival of freshwater species during the K-Pg mass extinction event.
Later, in the Paleocene (Pg1), richness was low in Asia, Europe and South America, representing all regions for which informative data are available except for North America. This low diversity was driven by the extinction of major clades that had survived the K-Pg boundary (e.g. compsemydids, pleurosternids, adocids, bothremydids) coupled with relatively few first appearances of new clades in this interval. The reasons for high richness in the Paleocene of North America compared to all other sampled regions are not clear, but might relate to differences in the environments available or diversification dynamics between regions. Some modern lineages began to appear during this interval, including testudinoids in both North America (platysternids; Vlachos 2018) and Asia (testudinids, the extant tortoise clade; Claude & Tong 2004).
Eocene. In the early Eocene (Pg2) turtles had richness levels close to, if not higher than, those of the Late Cretaceous. Only the northern hemisphere continents are wellsampled enough to track from this point forwards in our subsampled curves but Asia is still poorly sampled, obscuring the continuity of richness patterns in this region. Higher temperatures during and following the Paleocene-Eocene Thermal Maximum (PETM) opened high latitude dispersal routes (e.g. Beringia) for interchange between Eurasia and North America (Godinot & de Broin 2003), which corresponds with the high recovered richness of Pg2. In Asia, testudinoids continued to diversify and, based on current fossil evidence, dispersed to other continents from there, though whether tortoises truly originated in this region is debated (Hofmeyr et al. 2017).
In Europe, warmer temperatures appear to have facilitated more dispersals of pleurodires from Gondwana via Africa (P erez-Garc ıa 2017). Carettochelyids (a clade of soft-shelled turtles), which were endemic to Asia during the Cretaceous (but are restricted to Australasia today), appeared in Europe and North America in the Ypresian (Joyce 2014). North American assemblages continued to include high turtle richness. The warm climate enabled turtles to colonize high latitudes for the first time since the Late Cretaceous and in North America~10 clades of turtles reached at least 75°N alongside squamates, amphibians and crocodylians (Eberle & Greenwood 2012;Nicholson et al. 2016: see also Fig. 6). For comparison, the modern maximum latitude of extant turtles in North America is approximately 58°N (Chrysemys picta; Angielczyk et al. 2015). It is unfortunate that the southern hemisphere is so poorly sampled at this time, as other terrestrial vertebrates appear to have had similar, southerly range expansions due to warmer temperatures during the PETM in South America (e.g. lepidosaurs; Head et al. 2009;Albino 2011). Given modern pleurodire environmental requirements (P erez-Garc ıa 2017) it seems likely that they might have responded favourably to warmer conditions in this region, though further evidence is required to test this hypothesis.
Turtle richness appears to decrease slightly in North America during the early middle Eocene (Lutetian, Pg3), whilst in Europe it decreases sharply, in a pattern very similar to that for non-marine lepidosaurs (Cleary et al. 2018). For the European data, the countries being sampled differ a little from the earlier Eocene and so the reason for this decline could be a change in the geographical areas being sampled. As theorized for lepidosaurs (Cleary et al. 2018), a decrease could also be linked to the temperature decline that occurred at the end of the early Eocene Climatic Optimum (EECO; Woodburne et al. 2009). Mannion et al. (2015 found the inverse pattern for non-marine crocodylians, with Pg3 containing a higher richness of taxa than the two surrounding bins. In North America, turtle decline in the Lutetian (containing most of the Uintan North American Land Mammal Age (NALMA)) is very minimal (possibly negligible), and might have been due to a decrease in aquatic environment availability and increasing aridity, which was suggested by Hutchison (1992) based on the reduction in carapace size of aquatic turtles in the Uintan, whereas tortoises (testudinids) remained mostly unaffected. A decrease in collections in Pg3 relative to Pg2 is caused partially by the difficulties in correlating NALMAs with standard European stratigraphic stages. For example, the Uintan NALMA overlaps both the Lutetian and Bartonian stages. If an occurrence lacks sufficient biostratigraphical refinement for assignment to a standard stage those occurrences were not included in these analyses.
High richness is recovered in the late Eocene (Pg4) on both continents, with European richness bring roughly the same as during Pg2 and the Late Jurassic. Asia is well-sampled enough to be recovered again, with richness levels comparable to those in Europe. There are fewer monogeneric localities in Europe in Pg4 compared to Pg3, and numerous collections are known from the UK, but the geographical spread of collections does not change drastically. In North America, baenids made their last appearances in the fossil record in the late Uintan NALMA (early in Pg4) and are assumed to have become extinct by the end of the Eocene (Joyce & Lyson 2015), while testudinids became more common. This coincided with increasing aridity throughout the later Eocene (Hutchison 1992). Precipitation (and the associated availability of aquatic habitats) is the most important factor determining species richness in extant freshwater turtles (Iverson 1992;Waterson et al. 2016) and increasing aridity or changes in precipitation patterns could have had differential effects on the distribution and richness of aquatic versus terrestrial turtles. The interaction of various environmental variables is complicated, however, and aridity alone may not necessarily cause aquatic taxa to decline in richness if sufficient riparian or paludal habitats remain.
Oligocene and the Grande Coupure. The Eocene-Oligocene boundary (Pg4-Pg5) has been the subject of intense interest due to the Grande Coupure ('great break'), an extinction or turnover event recognized for many groups of terrestrial vertebrates, particularly mammals (Stehlin 1909;Hooker 1989;Meng & McKenna 1998;Costa et al. 2011;Eronen et al. 2015 It coincides with a sudden decrease in temperature, increased seasonality and the onset of Antarctic glaciation (Ivany et al. 2000;Zachos et al. 2001;Liu et al. 2009). In our analyses, the turtle richness of North America greatly decreased, but that of Europe increased during this event.
Relatively few Oligocene deposits have been sampled in Asia, so examining the effects of this climatic shift is impossible with our current data. It appears, however, that large parts of the continent experienced the same transition from warmer temperate conditions to an arid or semi-arid cooler climate with a short rainy season and fewer permanent bodies of water (Meng & McKenna 1998;B€ ohme 2007;Sun et al. 2014). This probably affected vertebrate assemblages and B€ ohme (2007) noted that late Eocene Mongolian faunas contained many aquatic forms that do not appear in the Oligocene. Nevertheless, it will be interesting to investigate the effects of the Grande Coupure further when more data are available from Asian Oligocene localities.
In Europe, where the Grande Coupure was first recognized (Stehlin 1909;Hooker 1989) pleurodires are poorly sampled, which is frustrating as the stricter environmental requirements of extant members would make an interesting comparison. Cryptodires increased in richness, contrary to what is expected during a dramatic shift in climate (Fig. 4). Little attention has been paid to European turtle responses to the supposed event, but their overall richness increased to its highest for the entire Triassic-Palaeogene interval. Geoemydids and trionychids were still relatively abundant even in the early Oligocene, which suggests that European aquatic turtles were not affected by these global environmental changes to the same extent as other major vertebrate groups.
North American turtle faunas underwent various changes at this time, with emydids and trionychids disappearing from the Western Interior after very rare occurrences in the early Oligocene, leaving testudiniddominated assemblages in the continental interior. They did not go extinct, as some specimens from this time exist in museum collections (PAH, pers. obs.) but are certainly less common. Testudinid body size decreased across the Eocene-Oligocene boundary (Hutchison 1992), corresponding with global cooling and drying (Hutchison 1982), but temperatures were still relatively warm and testudinid geographical ranges remained unchanged. The differences in responses between continents could indicate a number of possibilities, including that Oligocene sampling is much poorer. There is a lack of turtle-bearing localities in the late Oligocene of North America (Arikareean NALMA; Joyce & Bourque 2016) and those that do exist are often poorly constrained temporally, which could exacerbate the steep drop in richness found here compared to the situation elsewhere in the northern hemisphere.
Conversely, or perhaps additionally, the sudden global temperature drop could have caused migration (or extirpation) of the aquatic taxa out of the sampled areas to more southerly latitudes. In North America, for example, kinosternids, chelydrids, emydids and trionychids are found widely in Recent faunas, and it is hypothesized that these turtles persisted at lower latitudes (e.g. in the Gulf Coast region) until conditions were more favourable further north (Hutchison 1992). Evidence for a brief warmer period in the late Oligocene is suggested by the presence of pleurodires in South Carolina and Florida (Weems & Knight 2012; Bourque 2016), for example. This was also theorized to be the case for North American squamates (Smith 2011) andFerreira et al. (2018) found a reduction in the phyletic diversity of freshwater pleurodires versus coastal and marine taxa at this time, which lends more credence to this hypothesis. However, the lack of more southerly localities currently limits testing of these hypotheses.
Another possible factor which should be explored with tetrapods on a broad scale, is whether the North Atlantic current creates a more equable climate in Europe than in the North American interior during the Oligocene. Nevertheless, the climatic changes across the Eocene-Oligocene boundary probably limited the geographical ranges of turtle taxa with stricter temperature and environmental tolerances (e.g. pleurodires) and allowed those with wider temperature tolerances and semi-terrestrial to terrestrial life habits to spread (testudinoids).

Correlations of turtle richness with potential driving variables
We compared models with various combinations of explanatory variables (sampling, temperature, land area, sea-level) to face-value turtle richness using generalized least-squares regressions. Of these, eight models were found to be a better fit than the null model based on AICc values ( Table 2). The best four are turtle-bearing formations (TuBFs) only, followed by single combinations of turtle-bearing formations with d 18 O, non-marine area (NMA) and sea level. Even though they were chosen as the best models, many had relatively low R 2 values (between 0.3 and 0.5), and so do not provide a complete explanation of the variance in fossil taxon counts.
Our analyses suggest that turtle richness dynamics are complex and do not correspond strongly or simply to common drivers of richness, such as palaeotemperature and non-marine area. The number of turtle-bearing formations was expected to rank highly in the model comparisons, as a large number of formations producing turtles should be more likely to produce specimens that can be identified to genus level. When this is removed, the best models are the tetrapod-bearing collections (TBCs) only and tetrapod-bearing formations (TBFs) only models, though the R 2 values are very low (see Cleary et al. 2020). Given that all of the best models selected here included a sampling proxy, it is clear that sampling has a sizable effect on the patterns recovered from the turtle record, as expected from other analyses of the fossil record. The number of fossil-bearing collections, as a sampling proxy, is recovered as the best explanatory variable, or as a component of a model alongside non-marine area, temperature and/or Lagerst€ atten, for a number of vertebrate groups including sauropodomorph dinosaurs richness. Despite being ectotherms, turtles do not respond as strongly to temperature as expected for a clade reliant on temperature for many of their life habits (e.g. reproduction). North American crocodylian richness also did not correspond closely with changes in palaeotemperature through time (Mannion et al. 2015). It is likely, as discussed above, that precipitation and the availability of freshwater environments were more important to nonmarine turtle (and possibly crocodylian) richness than temperature influences (as in modern turtle faunas; Iverson 1992), since a high proportion of Triassic-Palaeogene turtles were aquatic, but these factors are hard to measure directly from the fossil record. Ecological niche modelling of the Maastrichtian versus the Recent found that precipitation drove the niche expansion of aquatic turtles, while thermal limits provided important constraints on niche ranges (Waterson et al. 2016), but additional modelling is required to determine the influence of precipitation on niche occupancy and breadth in other time slices.

CONCLUSION
Temperature, rainfall and habitat heterogeneity are thought to drive patterns of richness in turtles both now and in deep time (Iverson 1992;Holroyd & Hutchison 2002;Waterson et al. 2016). However, our analyses demonstrate that observed turtle richness through time is also significantly influenced by sampling biases that must be accounted for when interpreting the temporal and spatial patterns of turtle diversification and extinction events. Our analyses show where sampling is not yet sufficient to extract macroevolutionary signals, but also identifies previously unrecognized richness peaks and shared long-term patterns of diversity change within and between regions.
Turtle richness was low and poorly sampled until the middle-Late Cretaceous, with the exception of a few instances such as the Late Jurassic to Early Cretaceous of Europe. Richness was high following the Cretaceous Thermal Maximum, and it is tempting to associate this with increased opportunities for immigration and dispersal events at higher latitudes. However, in all recovered continents with a sequential K7-K8 pattern, richness declined to a lower level by the Maastrichtian (K8). Non-marine turtles appear to have been relatively unaffected by the K-Pg mass extinction, besides the extinction of most land turtles, with high survivorship in North America and stable richness across the boundary in South America. We find different results to those of Nicholson et al. (2015) for Mesozoic turtles, resulting from major recent taxonomic and stratigraphic revisions, as well as the addition of new taxa to the Paleobiology Database. In the Palaeogene, richness was only recoverable for northern hemisphere continents after the Paleocene (Pg1). North American richness remained high but decreased slightly until a sharp increase in the late Eocene. European richness fluctuated in a series of peaks and troughs, corresponding to those recovered for lepidosaurs and suggesting similar driving factors for both groups, which might be climate, sampling factors not corrected for by SQS, or a combination of these. Asia, where recovered, had a high richness. At the Eocene-Oligocene boundary North America decreased in richness while Europe increased, casting doubt on the existence of the Grande Coupure event in turtles. GLS analyses suggest that sampling (particularly the number of turtle-bearing formations) has some influence on the recovered face-value richness through time, but much variation remains unexplained by sampling alone and further investigation is required. Ecological niche models provide an alternative method for examining potential driving factors, like precipitation, which are difficult to extract directly from the fossil record. However, even in well-sampled intervals fossil localities represent only a portion of the total land area for that time period (e.g. the Western Interior representing all of North America), and evenness can affect richness estimates strongly. It is therefore crucial to examine richness at local scales, in order to avoid extrapolating what could be purely local richness patterns to continental or global scales. We must rigorously assess the data that we do have, make new field collections and increase assessments of untapped museum collections to improve our datasets for the future.
These caveats aside, our data indicate that responses to current and future anthropogenic warming will probably differ among major turtle clades if the past is an accurate guide to the present and future. Pleurodires, for example, will probably benefit from a warmer climate, as this seems to have been a factor that facilitated their dispersal on several occasions in Earth history and extant taxa are known to prefer warmer temperatures. However, as freshwater turtles they are also reliant on rainfall-driven habitat availability, and so clade responses to current climatic warming will depend on the interaction of numerous factors. Ultimately turtle survivability will depend on the rapidity of their adaptation to changing conditions and their ability (and opportunity) to migrate to locales within their temperature tolerances. The data we have gathered on past richness provides a useful benchmark with which to model and forecast future responses, and the incorporation of new data in the future will only increase its utility in elucidating the factors that most influence turtle richness. funding