Data collection
We compiled the species list of global chelonians through a critical synthesis of four authoritative taxonomic frameworks: (1) the paper of Global Conservation Status of Turtles and Tortoises (Order Testudines) published in 2018 (360 species included)10, (2) Turtles of the World: Checklist and Atlas (9th Ed., 2021; 357 sp.)45, (3) the Reptile Database (version September 2024; 371 sp.)46, and 4) the IUCN Red List (version 2024-02; 272 sp.). For the Chelonoidis nigra complex, the Turtle Taxonomy Working Group (TTWG)45 and the Reptile Database46 recognize them as subspecies, whereas TFTSG10 elevates each to species rank. We adopted the TFTSG classification based on diagnostic morphological, behavioral, and genomic divergence identified by a recent study47. Seven recently proposed species (e.g., Chelodina ipudinapi, Elseya auramemoria) were excluded pending formal validation by TTWG. The final checklist comprises 378 validated species (Supplementary Data 1), with taxonomic decisions cross-verified through correspondence with Dr. Anders G.J. Rhodin.
The IUCN Red List evaluated the threat status of 272 chelonian species, and among which 138 species (>50%) were assessed before 201410. Based on the 2018 IUCN Red List, the TFTSG Red List 2018 assessed unevaluated tortoises and freshwater turtles and re-assessed these species with outdated assignments10. Overall, the threat status of 360 species was evaluated, including 7 Extinct, 187 threatened, and 136 non-threatened species. When compared, the threat status of 66 Data Deficient (DD) and Not Evaluated (NE) species in the IUCN Red List was evaluated, which included 45 Least Concern (LC) species, 5 Near Threatened (NT) species, 12 Vulnerable (VU) species, 2 Endangered (EN) species, and 2 Critically Endangered (CR) species10. Besides, TFTSG reevaluated 32 data-sufficient species assessed by the IUCN Red List, lifting the extinction risk of 16 species10. Considering the higher degree of data completeness of TFTSG over the IUCN Red List, we chose the TFTSG Red List as our primary dataset. We assigned the Red List status of species with a score to represent their extinction risk. Following previous studies11, we categorized the extinction risk in descending order as Extinct or Extinct in the Wild (EX = 5), CR = 4, EN = 3, VU = 2, NT = 1, and LC = 0.
In addition, since the TFTSG report has not updated the threat status since 2018, we took three additional steps to enhance the dataset. First, we updated the threat status of chelonians based on the TTWG 2021 report. Second, we updated the threat status of species that was revised by the IUCN Red List only after 2018, as these updates were also conducted by the IUCN-SSC TFTSG. Lastly, while the global population of green turtles (Chelonia mydas) is reported to be downlisted to LC in 20257, the status of other marine turtles recorded in the IUCN Red List is significantly outdated and does not accurately reflect their current risk status. Therefore, we averaged the risk status of their subpopulations as assessed by the IUCN-SSC Marine Turtle Specialist Group. Specifically, the loggerhead turtle (Caretta caretta) was assessed as VU (coded as 2) in 2015, with all ten subpopulations ranging in status from LC (0) to CR (4). The average risk status of these subpopulations was 1.7, and hence we rounded this number and assigned the threat status of loggerhead turtle as VU. However, for leatherback sea turtle (Dermochelys coriacea), its global threat status was assessed as VU in 2013. We assigned its threat status as CR (4) based on the average values of its seven regional assessments (assessed between 2013 and 2019), which was 3.8. Our assessment aligns with recent studies that leatherback turtles are experiencing declines in several regions24,48. Overall, the updated dataset included 335 data-sufficient species and was used in our analyses, comprising 7 Extinct, 73 Critically Endangered, 56 Endangered, 60 Vulnerable, 38 Near Threatened, and 101 Least Concern species (Fig. 1, Supplementary Data 1). (In sensitivity analyses, we used the IUCN or TFTSG Red List assessment, respectively as the response variable and found largely similar results, which are detailed in Supplementary Table 3).
We collected several traits that were shown by previous studies to be predictors of extinction risk in the studies of chelonians and other vertebrate groups (a summary in Supplementary Table 4), including activity time, diet, insular endemism, maximum lifespan (years), maximum body length, mean clutch size, and geographical range size. Specifically, we categorized activity time into nocturnal (0), crepuscular (1), diurnal (2), and cathemeral (3, day and night). Diurnal and cathemeral species have activity patterns that closely coincide with human activities. Secondly, previous studies have shown that carnivores are often at a higher risk of extinction than omnivores and herbivores23. Based on this, we quantified the chelonians as omnivores (0), herbivores (1), feeding on invertebrates (2), and feeding on vertebrates (3). Insular endemism was classified as either island endemic species (1) or not (0). Since species with a larger body size and a smaller geographical range size appear to have higher extinction risk49, we compiled the data on the maximum straight-line carapace length (mm) of adults to represent body size, and geographical range size as the extent of occurrence (km2). Species with smaller clutch sizes are likely to be prone to extinction due to a weak recovery potential from small population sizes50. Accordingly, we quantified clutch size as the mean number of eggs per clutch.
We recorded geographical range size from a recent study51 and supplemented it with information from the IUCN Red List2 and the TTWG 202145. Body size was collected from the TTWG 202145 and supplemented with the published species traits of chelonians52,53 as well as related publications54. The information on activity time, diet, insular endemism, maximum lifespan, and mean clutch size was derived from the two global datasets of the species traits52,53, the book written by Lovich and Gibbons55, the Reptile Database56, the IUCN Red List2, and related papers46. We collected trait data for all 378 species, except for a small proportion of missing values in the following categories: activity time (26 species), diet (66), maximum lifespan (51), and mean clutch size (23). We then imputed missing values using the impute function in R package funspace57, which fills in missing values (of continuous or categorical ones) using trait–trait correlations through a machine-learning process that accounts for phylogenetic relatedness. This methodology has been recommended in recent papers58, and this imputation procedure has been shown to produce reliable results with minimal error rates59. In this study, we run this gap-filling approach 100 times and performed the analyses on each imputed dataset60. We then pooled the results at the end using Rubin’s rule61, which adjusts the standard errors and provides more accurate parameter estimates by incorporating both within-dataset variance and between-dataset variance.
To represent the phylogenetic relationships among chelonians, we constructed a phylogenetic tree for chelonians based on Colston et al21. This tree recognizes 357 species of extant or recently extinct chelonians and contains all 335 data-sufficient species in the combined Red List. However, since we need to impute the missing values in species traits also for data-deficient species (20 sp.), we employed two different methods to address this issue. For species that were previously classified as subspecies or subpopulations, we positioned them alongside their sister species (the vast majority are in this situation. e.g., Galapagos giant tortoises). For some newly discovered species, we used the add.species.to.genus function (phytools)62 to incorporate the missing species midway along the branches of their respective genera. To assess whether the correlates of species extinction risk differed among habitats, we grouped chelonians into aquatic (freshwater and marine) and terrestrial species based on the classification of the IUCN Red List, the published species traits dataset, books, and related papers, detailed in Supplementary Data 1. Finally, we used the Reptile Database to align species names (i.e., matched the species names across all datasets) and ensure consistency across all datasets.
To determine the influence of external environmental factors on the extinction risk of chelonians, we collected four climatic variables (mean annual temperature, mean annual precipitation, temperature seasonality, and precipitation seasonality) by calculating the mean values of these bioclimatic data (WorldClim 1, 1970-2000)63 across species distribution ranges. In fact, many sea turtles have a global distribution across tropical, subtropical, and temperate regions, making it challenging to derive a climate metric64. Nevertheless, we measured these variables for sea turtles across their ranges, as temperature can influence the abundance and distribution of food resources in the oceans65, and to ensure consistency with terrestrial species. We did not consider other climatic variables from WorldClim website (e.g., minimum or maximum temperature) because they are not only highly correlated with the mean values but also highly variable between years. In addition, extreme weather events, such as heatwaves and droughts, are occurring with increasing frequency and intensity, which may lead to even greater biological consequences than those caused by changes in climate means29,30. Hence, we measured the Standardized Precipitation Evapotranspiration Index (SPEI) to represent extreme droughts66. This index is one of the most widely used tools for monitoring and quantifying meteorological droughts because it accounts for the impact of evaporative demand on drought severity. We averaged the historical period (1971-2000) of the SPEI raters (0.5° × 0.5°) and transformed it into a 5-km Mollweide projection raster. A lower SPEI value indicates a higher risk of drought exposure to species. Furthermore, we selected the Warm Spell Duration Index (WSDI) to represent heatwaves, which is defined as the annual count of days with at least six consecutive days when the daily maximum temperature exceeds the 90th percentile of maximum temperature for the corresponding calendar day, centered on a 5-day sliding window during the baseline period67. We used the same procedure to produce a 5-km WSDI raster. Lastly, we mapped these two extreme climates across species distribution ranges.
Additionally, we obtained the year of description as a variable45 since species described earlier may have a broader distribution, potentially leading to a reduced risk of extinction68. We grouped the chelonians into fourteen families based on the classification of the TTWG 202145. For biogeographic distributions, we made our results comparable to global biogeographical studies on chelonians69, by including Afrotropics, Australasia, Indomalaya, Nearctic, Neotropics, Palearctic, and Marine. We classified species into those realms by intersecting their ranges with the shapefile of the World Wildlife Fund for Nature (WWF) biogeographical realms. We calculated the proportion of species distribution range belonging to each biogeographical realm using ArcGIS 10.7. A species was assigned to a realm when at least 80% of its range overlapped with that realm, which left 14 unclassifiable species without any realm designation due to their widespread distribution (Supplementary Data 1).
Lastly, we calculated the functional distinctness of species using the functional traits, including activity time, body length, clutch size, diet, and lifespan. We computed the function distance matrix of these variables and then employed the distinctness function in the R package funrar70 in R 4.3.271 to quantify the average functional distance from a species to all other chelonians. It is computed as:
$${D}_{i}=\frac{{\sum }_{j=0,\,i\ne j\,}^{N}{d}_{{ij}}}{N-1}$$
(1)
Where Di is the functional distinctiveness of species i, N is the total number of species in the chelonian group, and dij is the functional distance between species i and j. The values range from 0 to 1, indicating the functional distinctness of species from low to high. Following Colston21, we calculated the evolutionary distinctness of species within the tree using the R package picante (“fair.proportion” metric)72. The metric quantifies the weighted sum of the branch lengths along the path from the root of an ultrametric tree to the tip, with weights assigned as the reciprocal of the number of tips that ultimately subtend that branch72. This analysis measures the evolutionary isolation of each species, with high values of evolutionary distinctness representing relict lineages on long branches72.
The populations of chelonians have been heavily threatened by overexploitation and trade, increasing habitat loss and fragmentation, environmental pollution, invasive species, and global climate change6,15. We gathered threat information for each species from several sources (Supplementary Data 1): the IUCN threat classification scheme, which provides a detailed description of the threats faced by each species but only covers ~ 270 species that they assessed; the Biology of Freshwater Turtles and Tortoises compiled by the IUCN/SSC TSTSG, which details the threats for 171 chelonian taxa; trade records from the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES); the reports by the Turtle Conservation Coalition (TCC)73; and books55. We reclassified threat types into five basic and recognized categories, including human exploitation (pet trade, human consumption, medicine), habitat destruction, pollution, invasive species, and climate change (Supplementary Table 5)6,8. Some other threats, such as roadkill or sacrificial purposes, were excluded due to minimal records. We constructed a species-by-threat matrix using these seven threat types, with 1 representing affected and 0 representing not affected. We derived the number of human threats for each species by summing the total records documented for a focal species (0 to 5).
Non-random extinction risk
Based on the extinction risk assessment of the combined Red List, we conducted a binomial distribution test to examine which families have a higher or lower proportion of threatened species than expected by chance alone. Overall, 335 data-sufficient species (from LC to EX) were grouped into two categories: at-risk/threatened species (VU, EN, CR, including EX) and non-threatened species (LC, NT). The null hypothesis posits that threatened species within each family occur randomly, implying that the proportion of at-risk species is equally distributed among families. In this case, the probability of K threatened species in a family with N species should follow a binomial distribution74, where the total proportion of threatened species in all families is 0.585 (196 divided by 335). Since the binomial distribution test was conducted independently for each family, we used the Dunn-Sidak method to calculate the adjusted critical value75. In addition, since the binomial test has relatively low statistical power when the taxonomic family is small76, we categorized small families into high-risk or low-risk groups when they contained 100% or 0% threatened or extinct species (those with > 0% and <100% were omitted). To identify the high-risk realms, we also performed the binomial distribution test for species in all biogeographical realms (including species in the unclassifiable group). Lastly, we linked the species families to the realms where they belong to identify the conservation priorities.
Correlates of extinction risk
To test for correlations of biological traits and external factors with the extinction risk of chelonians while accounting for the influence of shared ancestors among species, we conducted Phylogenetic Generalized Least Squares (PGLS) tests with Brownian motion (the Ornstein–Uhlenbeck model gave very similar results, Supplementary Table 6)77. The data analysis followed four steps. First, we employed univariate PGLS analysis to test the effect of each of the twelve ecological factors on the extinction risk category of species. Secondly, we selected the important variables with p < 0.1 in the first step and removed the variables with a variation inflation factor (VIF) > 5 to reduce the collinearity among predictor variables (all VIFs <5; Supplementary Table 7). Thirdly, we built models using various combinations of these variables and ranked all models using the information-theoretic approach based on Akaike’s information criterion corrected for small sample sizes78. We used the Akaike weight (wi) to compare model performance. If there were several competing models (ΔAICc ≤ 2), we conducted a model-averaging approach to calculate the relative importance of each variable within a 95% confidence interval76. The PGLS analysis was carried out using the package caper77.
We excluded DD and NE species from the univariate and multivariate PGLS analyses and log-transformed the continuous variables to achieve a normal distribution. We conducted the analyses using the combined Red List with categorical threat status ranging from 0 (LC) to 5 (EX). We conducted analyses for aquatic and terrestrial species since they generally show extensive variation in habitat preferences and life-history strategies. We also repeated all the above analyses for each realm using the combined Red List assessment (0–5) due to the availability of more adequate data for the analyses. To minimize the impact of any changes in the threat status of species on our findings, we further categorized the threat status into two groups: threatened (VU, EN, CR, EX, coded as 1) and non-threatened (NT, LC, coded as 0). We employed phylogenetic logistic regression models (PGLM, using the phylolm package)79 to identify the significant predictors of the probability of a species being threatened. We refrained from performing multivariable PGLS and model selection for species in the Palearctic region due to the limited sample size (15 sp.). We reported the standard errors, within-dataset variance, and between-dataset variance (<1.0 × 10−4 for most cases) across 100 datasets for both univariate and multivariate PGLS analyses in Supplementary Data 3 and 5.
The risk probability of DD and NE species
To predict the extinction risk probability of DD and NE species in the combined Red List, we combined the results from the binomial distribution test and the model-average analyses. At the family level, we grouped all species into high- (significantly higher than average), medium- (no difference from the average), or low-risk (significantly lower than average) family groups based on the results from the binomial distribution test. We then extracted the lists of DD and NE species belonging to each of the three groups. At the species level, we fitted a logistic regression model on the extinction risk of all data-sufficient chelonians using the significant variables identified through the model-average analysis following Senior et al.80. We then used this model fit to predict the probability of DD and NE species becoming threatened (1) or unthreatened (0) status using the predict.glm() function. The model predicted the probability of the DD and NE species having a threatened status (1), which ranged from 0 to 1. We then grouped these value estimates into three categories (0.75 ~ 1, 0.5 ~ 0.75, <0.5). The classification is highly artificial, while as suggested by Senior et al.80, a value > 0.75 suggests that species have a higher probability to be threatened in the future because of characteristics shared with congeneric species that have become rare or gone extinct. Lastly, we integrated the family-level and species-level results and formatted them into a matrix. We considered that the unassessed species belonging to moderate or high-risk families and having an estimated value ≥ 0.50 had a higher extinction risk.
The rates of trait evolution compared to future challenges
To predict how well chelonians can adapt to climate change, human expansion, and forest degradation in the next century, we calculated the rates at which traits were evolving and compared them to the rates of change of these threats, separately. We calculated the rate of trait evolution, which represents the variance among ln-transformed trait values per million years32. For two species, the evolutionary rate in felsens is:
$${ER}=\frac{{({x}_{1}-{x}_{2})}^{2}}{2t}$$
(2)
Where x1 and x2 represent the ln-transformed trait value of two species, t denotes the divergence time of t million years. A higher ER value indicates a faster rate of evolution32. The model considers that each trait is evolving at a constant rate with Brownian motion, so we used the BM model for all traits32.
To measure the rates of climate change, we first extracted historical climate data on mean annual temperature (1981–2000, with a spatial resolution of 2.5 minutes) for all non-marine species across their ranges. We then calculated the projected mean temperatures using three climate change models (ACCESS-CM2, MIROC6, MPI-ESM1-2-HR) for the period 2081–2100, under both the best-case (SSP1-2.6, Sustainability) and worst-case (SSP5-8.5, Fossil-fueled Development Scenario) climate change scenarios for each species. We calculated the rates of climate change per century for both scenarios using the aforementioned equation, and subsequently obtained the average and standard deviation values for each scenario. Similarly, we calculated the projected extreme droughts (SPEI) and heatwaves (WSDI) using the historical extreme climate data (1981–2000) and projected climate data (2081–2100) under the two climate change scenarios. We used the same climate change models and the equation above to determine the changes (Mean and SE) in the intensity of extreme droughts and heatwaves. Furthermore, to represent the change in intensity of human activities, we collected data on the world population and land-use projections of SSP1 (2100) and SSP5 (2100), as well as their base-year raster data of 2000 from the NASA Socioeconomic Data and Applications Center81 and Land-Use Harmonization82, respectively (only one model available). The Land-Use Harmonization provides seven land-use states (e.g., proportion of forested primary land, managed pasture, rangeland and urban land) and we only considered the changes in the total fraction of primary and secondary forested land types. We calculated the rate of population change and forest cover change per century in both the best-case and worst-case scenarios, respectively. To standardize all values for comparison, we divided the estimated evolutionary rate of each trait by 10,000 to convert the unit to centuries.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
