In order to evaluate (Vrba, 1987, 1992) resource-use hypothesis in modern Testudines, our study was conducted on a global scale. Based on the time-calibrated phylogeny of Thomson et al. (2021) and geographic distribution data from Bonin et al. (2006) and Rhodin et al. (2021), supplemented by information from scientific literature for specific species (Fritz et al., 2008; Van Dijk and Rhodin, 2010; Thomson et al., 2011; Singkily et al., 2018; Caramaschi, 2021). We determined the biomes inhabited by each species, based on Walter (1970) biome classification, modified by Hernández Fernández (2001) ( Figure 1 ).

Following Hernández Fernández (2001) a biome was considered occupied by a species if it represented 15% or more of the species’ geographic range. In addition, due to the existence of species with extensive ranges, a biome was also considered occupied if the species inhabited 50% or more of a particular climatic domain, defined as a contiguous terrestrial area consisting of a single biome. Marine species were assessed based on the biomes associated with their primary breeding areas. These criteria allow us to represent the species’ adaptive capacities while maintaining their climatic specificity and facilitate comparison with previous works that used the same methodology (Hernández Fernández and Vrba, 2005a; Moreno Bofarull et al., 2008; Cantalapiedra et al., 2011; Gómez Cano et al., 2013; Menéndez et al., 2021; Gamboa et al., 2022; Hernández Fernández et al., 2022; Pelegrin et al., 2023).

We calculated the Biomic Specialization Index (BSI) for each species (Hernández Fernández and Vrba, 2005b), which represents the number of biomes occupied by the species. Once calculated, it allows the classification of species into three categories (Hernández Fernández and Vrba, 2005b): biome specialists, occupying a single biome (BSI = 1), moderate generalists (BSI between 2 and 4), and extreme generalists (BSI ≥ 5).

The biome occupancy data for all species, based on the taxonomy of Rhodin et al. (2021), were compiled into a matrix representing the presence (1) or absence (0) of each species in each of the 10 terrestrial biomes. Subsequently, 10,000 Monte Carlo simulations (Gotelli, 2000) were conducted to test for significant differences in the observed patterns of biome occupation compared to a random distribution of the species (Hernández Fernández and Vrba, 2005b; Hernández Fernández et al., 2022). Because the specific ecological characteristics of each biome affect species richness, there is no reason to consider the same number of species for all biomes in the null model (Jetz and Fine, 2012). The randomization was performed by randomly placing species in the biomes but limiting the species richness of the biomes according to the observed values. The significance (p-value) of the observed trends was assessed by comparing observed and simulated values (Hernández Fernández and Vrba, 2005b).

Given that families are significant evolutionary units (lineages) (Humphreys and Barraclough, 2014), Monte Carlo analyzes were conducted not only for all Testudines but also independently for their families. In order to ensure the robustness of the statistical analyses, only families with a minimum of 10 species were included. Specifically, seven families were analyzed (Chelidae, Geoemydidae, Emydidae, Kinosternidae, Pelomedusidae, Testudinidae, and Trionychidae), which represent half of the extant turtle families and 92% of all species.

Finally, to determine whether the observed pattern was reflected in species ecology, we conducted independent Monte Carlo analyses for terrestrial and freshwater species, excluding marine species due to their representation being fewer than 10. All analyses were conducted using the R software (R Core Team, 2024).

The relationship between diversification rate (DR) and the biome specialization, represented by the number of biomes inhabited by a species and the specific biome occupied (in the case of biomic specialist species), was assessed. Diversification rates for each species were calculated using the evol.distinct function of the R picante package (Kembel et al., 2010) based on the phylogeny from Thomson et al. (2021), which includes all 14 living families, 90 of the 92 living genera (98%), and 279 of the 357 living species (78%) of species in the group. This function estimates DR values by considering the number of splits and internodal distances of branches from each species or tip to the root of the tree, giving greater weight to more recent branches and splits (Jetz et al., 2012). Thus, to evaluate the statistical significance of the relationship between diversification rates and the number of biomes occupied by species, we employed a phylogenetic generalized least squares (PGLS) approach using the R nlme package (Pinheiro et al., 2023) under a Brownian evolution model (Freckleton and Harvey, 2006). This method accounts for the expected covariance of data, incorporating phylogenetic relationships among species (Mundry, 2014).

Given that many of the traits influencing species’ capability to inhabit different biomes are related to its anatomy, physiology or behavior, all of which are heritable traits (Cava et al., 2019), it is crucial to incorporate a phylogenetic perspective in statistical analyses. This approach allows to differentiate signals of evolutionary processes from potential phylogenetic conservatism (Gamboa et al., 2022; Pelegrin et al., 2023). In order to ascertain whether biome specialist species exhibit notable differences in diversification rates contingent on their occupation of disparate biomes, we conducted a phylogenetic analysis of variance (phyANOVA) to test for significant differences in DR among groups of specialist species adapted to different biomes using a Brownian motion model of evolution (Freckleton and Harvey, 2006). This analysis was performed using the R packages phytools (Revell, 2012) and geiger (function aov.phylo; Pennell et al., 2014).

Our findings indicate that turtles are present in eight of the ten terrestrial biomes, with the highest levels of occupancy observed in five of these (BSI= 5). The frequency distribution of BSI exhibited a pronounced rightward skew ( Figure 2 ), with a relatively low mean BSI (BSI = 1.56). In total, 205 species (57.7%) are confined to a single biome, while 113 species (31.9%) are distributed across two biomes. Furthermore, only 1.4% of Testudines (five species) can be classified as extreme biome generalists ( Supplementary Table S1 ).

Monte Carlo simulations for all species revealed a significantly higher proportion of biome specialist species than expected by chance. In contrast, the number of moderate biome generalist species was lower or not significantly different than expected from the simulations. Among extreme biome generalists, species with a BSI=5 were found in a higher number than anticipated from the simulations, while species with BSI=6 and above were not observed or found in the simulations ( Supplementary Table S1 ; Figure 2 ). In the same way, the Monte Carlo simulations by ecological groups showed the same trends for terrestrial and freshwater turtles as previously observed ( Supplementary Tables S4, S5 ), with a greater accumulation of specialists than expected by chance for both groups.

On the other hand, the distribution of species across biomes reveals the Tropical deciduous woodland as the richest biome, with 221 species, followed by the evergreen equatorial rainforest with 112 species. The remaining biomes exhibited less than 59 species, with the steppe recording the lowest number at just 10 species. It is noteworthy that neither the taiga nor the tundra biomes were home to any species ( Supplementary Table S3 ). We found that this distribution closely resembles that of biome specialist species, with a notable concentration of species in the Tropical deciduous woodland, followed by the Evergreen equatorial rainforest. The only distinction is the absence of biome specialists in the Steppe ( Figure 3 , Supplementary Table S1 ). On the other hand, terrestrial species exhibited a notably high occupation in the subtropical desert and sclerophyllous woodland and shrubland, comparable to their presence in evergreen equatorial rainforest and savanna ( Supplementary Table S6 ). In contrast, freshwater species exhibited the same trend as observed in previous analyses, with a higher accumulation of specialist in tropical biomes ( Supplementary Table S7 ).

Regarding the Monte Carlo simulations per family, six of the evaluated families (Chelidae, Geoemydidae, Emydidae, Kinosternidae, Pelomedusidae, Testudinidae) showed the expected pattern according to the resource-use hypothesis. However, the family Trionychidae showed an unexpected distribution, characterized by a significantly higher percentage of species than expected with BSI=2 and fewer biome specialist species than expected based on the Monte Carlo simulations ( Figure 4 ). Extreme generalist species were found among most of these families, although their occurrence was higher than expected by chance only in two families (Geoemydidae, Pelomedusidae). Additionally, we found that in all biomes, except for the subtropical desert, the broadleaf deciduous forest and steppe, the number of specialist species was higher than expected by chance according to Monte Carlo simulations ( Figure 3 ; Supplementary Table S3 ).

Among families, Emydidae exhibited the highest diversification rate. Conversely, Platysternidae and Carettochelyidae showed the lowest rates ( Figure 5 ). Upon calculating the diversification rates for all species in the phylogeny, the information was grouped based on their Biome Specialization Index (BSI). Biome specialists showed significantly higher diversification rates compared to moderate and extreme generalists ( Figure 6 , Table 1 ). The phylogenetic generalized least squares (PGLS) test showed that there is a strong correlation between DR and BSI (p = 0.014, r² adjusted = 0.886) ( Supplementary Table S9 ). Moreover, an increase in the standard error is observed as the BSI increases, due to the greater variability and sensitivity to fluctuations when working with smaller data sets (Camacho-Sandoval, 2007).

Despite the higher diversification rate for biome specialist species, their DR values were not significantly affected by the particular biome inhabited by these species; no significant differences were found between the DR of the biome specialist species with respect to the biomes they inhabit (phylANOVA p = 0.760; Supplementary Table S10 ), supporting the observations in Figure 7 , where the boxplots indicate a clear overlap between them.

Reptiles are ectotherms, known for a clear tendency towards endemism due to their low metabolic rate, which reduces their food requirements and decreases their activity levels, thereby diminishing their need for movement compared to endotherms (Hayden Bofill and Blom, 2024). Thus, due to metabolic temperature dependence, the diversification of biome specialist species should be more pronounced than in endotherms (Machac et al., 2012). The Monte Carlo simulations confirmed a tendency towards the overrepresentation of biome-specialist species, in accordance with the initial premise of the resource use hypothesis ( Supplementary Table S1 ; Figure 2 ). This trend was also consistent across most families, except for Trionychidae, which, instead of biome specialists, showed a significant accumulation of species inhabiting two biomes ( Figure 4 ). Essentially, while these families are phylogenetically distinct lineages, they share closely related morphologies derived from a conservative body plan, which could constrain the range of traits or adaptations to inhabit biomes with varying characteristics (Hu et al., 2019), being cold and aridity the main characteristics that seem to have limited them. These findings are in line with studies conducted in taxa such as mammals (Hernández Fernández and Vrba, 2005b; Moreno Bofarull et al., 2008; Cantalapiedra et al., 2011; Menéndez et al., 2021; Hernández Fernández et al., 2022) or butterflies (Gamboa et al., 2022). This indicates that historical climate changes are still a more relevant factor in the diversification processes observed in animal taxa than the presence of diverse thermoregulation mechanisms or patterns of geographical dispersal (Rolland et al., 2018; Pie et al., 2021). Ultimately this could mean that biomic specialization within a lineage is an intrinsic phenomenon associated to the geographical dynamics of the evolutionary process itself (Vrba, 1980; Gómez Cano et al., 2013), expressed across lineages independently of their morphological or phylogenetic characteristics, as is the case with niche conservatism or extinction events (Raup, 1994; Ferrière et al., 2004).

Although turtles are highly influenced by environmental temperature for their activity and survival, and depend on a certain degree of humidity for their integument (Labra et al., 2008), some species can inhabit up to five biomes ( Figure 2 ). This underscores the capacity of turtles to develop strategies or adaptations that have enabled them to thrive in a relatively broad range of environmental conditions. Nevertheless, this discovery represents a level of biome generalism much lower than that observed in butterflies (Gamboa et al., 2022), mammals (Hernández Fernández et al., 2022) or birds (Pelegrin et al., 2023), groups that include species inhabiting up to eight or even more biomes (e.g., Papilio polyxenes, Lasiurus cinereus or Falco peregrinus). Such difference possibly is related to diverse morphophysiological constraints associated with their limited morphological diversity, which is attributed mainly to the presence of the carapace (Claude et al., 2003). This result is not surprising given that these ectothermic organisms struggle to maintain an adequate body temperature in cold climates solely through behavioral thermoregulation. Turtles particularly lack a thick layer of fat or other insulation like plumage or fur found in birds and mammals respectively (Terrien et al., 2011). This absence of insulation significantly hinders their adaptation to cold biomes, such as the Taiga and Tundra, where they are notably absent.

Thomson et al. (2021) found a low net diversification rate from the emergence of the lineage up to the Eocene, suggesting a higher prevalence of generalist species before this period due to their lower tendency to diversify. This trend changed during the Eocene Climatic Optimum when global temperatures rose and altered the distribution of biomes, with tropical biomes expanding worldwide (Jaramillo, 2020). This new biome configuration benefited turtles, increasing the diversification of biome-specialist species and resulting in an adaptive radiation of the group during this period.

On the contrary, it is intriguing that despite the occurrence of moderate biome generalist species being lower than expected by chance, as posited by the resource-use hypothesis, the presence of species with a Biome Specialization Index (BSI) of 5 exceeded simulations for five out of seven families ( Figure 4 ). This finding could mean that the occupation of five different biomes reflects possibly an upper limit considering some colonization events from ancestral tropical biomes towards subtropical ones. Due to niche conservatism, the physiological and adaptive constraints made possible biome colonization in some preferential directions according to environmental gradients of rainfall or temperature (Pelegrin et al., 2023); for instance, biome colonization processes from subtropical desert to steppe (aridity gradient) or tropical deciduous woodland to evergreen rainforest (humidity). For this reason, biome colonization events in cold environments (such as taiga and tundra) were not possible.

Among all turtles, the family Trionychidae is a notable exception, since it exhibits a higher accumulation of species with BSI=2, along with a lower incidence of biome specialist species than expected by chance. This particularity may be attributed to the specific ecophysiology of this family. Almost all its species dive to the bottoms of rivers and lakes, where they camouflage using their flat, smooth shells (which vary in color according to the substrate in their habitat) while waiting for hunting their prey (Ferri and Soccini, 2019). The intimate connection with aquatic ecosystems entails the requirement for adaptations that enable the species to withstand more extreme temperature fluctuations (Chessman, 2020). These thermal tolerance adaptations should be reflected in other freshwater families. For example, this could explain the high adaptability of Emydidae to temperate biomes. However, the unique BSI distribution of Trionychidae compared to other freshwater families might be attributed to their great adaptation to aquatic life, which makes them highly skilled swimmers, as evidenced by their broad, webbed feet, or their bimodal respiration and even living in areas near to the coastline (Wang et al., 1989; Ferri and Soccini, 2019). However, these adaptations significantly limit their ability to thrive outside water (in comparison to other freshwater turtle families), largely confining them to the aquatic environments they inhabit. Considering that their primary habitats are tropical dry forests and tropical rainforests, biomes that are ecologically and geographically intertwined, often sharing borders, it is both common and expected that the water bodies these turtles occupy extend across both biomes, thus explaining their BSI distribution.

We observed a high level of specialization across biomes, with five of them having more specialists than would be expected by chance ( Figure 3 ). This has been documented in small mammals, which exhibit rapid specialization due to high reproductive rates, limited dispersal capabilities, and shorter generation times (Moreno Bofarull et al., 2008; Hernández Fernández et al., 2022). Turtles have long generation times (McGaugh, 2012) requiring a long time for the offspring to reach sexual maturity (between 5 and 50 years depending on the species), as well as a significant longevity of adult individuals (Shine and Iverson, 1995). But exhibit very variable reproductive rates across species, with some Testudinidae females producing less than 10 eggs (Epperson and Heise, 2003), and sea species with more than 40 offspring, typically around 90, and sometimes up to 200 each year (Hirth, 1980). Additionally, despite the average home range for turtles is larger than in small mammals, ranging from approximately 23 ha in terrestrial species to 32 ha in freshwater species (Slavenko et al., 2016), compared to 400-600 m² for rodents (Ostfeld, 1985; Ssuuna et al., 2023) and an average of 3 ha for lagomorphs (Hulbert et al., 1996), it is much smaller than the home ranges of large mammals, which can easily exceed 400 hectares and even reach over 30,000 hectares in some felids (Lindstedt et al., 1986). These small home ranges can explain the higher generation of biome specialists among turtles. Since species with small territories can maintain genetically viable populations in limited areas, they are able to survive in small biome fragments, promoting speciation processes in vicariant populations of biome specialist species (Moreno Bofarull et al., 2008; Sonne et al., 2016). Additionally, lower activity and energy requirements of ectotherms lead to high rates of endemism (Hayden Bofill and Blom, 2024), which is similarly linked to greater diversification of biome specialists.

The low proportion of specialist species in the subtropical desert and the steppe suggests that the conditions in these biomes are very challenging for the survival and diversification of turtles. In the subtropical desert, this is associated with high temperatures and low humidity, coupled with significant temperature variation between day and night. These extreme conditions are particularly difficult for ectotherms, whose survival and biological activity heavily depend on environmental temperatures, especially affecting juvenile organisms and thereby population stability (Griffis-Kyle, 2016). On the other hand, very low winter temperatures characterize the steppe, a feature it shares with taiga and tundra, posing a significant barrier to the survival and diversification of turtles in these biomes.

The absence of species in the taiga and tundra biomes is noteworthy. The potential for turtle colonization is constrained by their limited ability to withstand desiccation and low temperatures for long time periods. Turtle species from temperate zones spend the winter hibernating in underground shelters (Ultsch, 2006), while freshwater species often prefer to hibernate at the bottom of lakes, ponds, and even streams (Ultsch, 2006), as seen in species from the family Emydidae (Cadi et al., 2004). Some, like Chelydra serpentina, hibernate in shallow bodies of water, enabling them to reach the water surface by extending their necks (Meeks and Ultsch, 1990). Additionally, these species have developed two main adaptations to survive low temperatures and extend their range into colder areas. Freeze tolerance protects them at temperatures as low as approximately 4°C, while supercooling allows some species to survive temperatures down to -12°C (Costanzo et al., 1995). However, the latter strategy is susceptible to ice nucleation, leading to rapid and uncontrolled tissue freezing upon contact with environmental ice, which is often lethal (Costanzo et al., 2008). This increases mortality and impacts their ability to maintain stable populations in environments with long and very cold winters, such as tundra, taiga and the most continentalized areas of the steppe biome. In addition to temperature, many habitats for turtles in these cold biomes (especially for freshwater species) present hypoxic or anoxic conditions due to surface soil and water freezing, restricting the places where they can hibernate and thereby limiting their distribution (Ultsch, 2006). Juveniles are particularly vulnerable to these conditions, which can significantly affect population recruitment and limit the range expansion of species into colder biomes (Costanzo et al., 2008).

Finally, when examining the biome distribution of specialists by ecological group ( Supplementary Tables S6, S7 ), we observe a greater concentration of terrestrial species in biomes such as the subtropical desert and sclerophyllous woodland and shrubland. These are relatively young biomes, both more recent than the group’s adaptive radiation during the Eocene (Buerki et al., 2012; Guerrero et al., 2013). In contrast, freshwater species were more closely associated with tropical biomes, which are ancestral ecosystems that experienced significant expansion during the Eocene. This could suggest that terrestrial turtles may possess greater adaptive capacity than their freshwater counterparts, or that colonization processes into new biomes occur more rapidly in species adapted to terrestrial environments.

We found no significant relationship among the DR exhibited by the specialist inhabiting different biomes ( Figure 7 ; Supplementary Table S10 ). This is partly due to the high variability in data from temperate biomes, especially in temperate evergreen forests and broad-leaf deciduous forests. The lower number of specialist species in these biomes ( Supplementary Table S8 ) increases variability and sensitivity to fluctuations (Camacho-Sandoval, 2007). In combination with the fact that low morphological variability within the group could limits the range of adaptations that can be developed to thrive and diversify in a biome (Bergmann and Irschick, 2012; Belmaker and Jetz, 2015). This constraint hinders strategies involving significant movements, such as migrations (with the exception of marine species) (Southwood and Avens, 2010), or adaptations associated with environmental temperature insulation, among others. This suggests that the strategies or adaptations developed by turtles to cope with different environmental conditions across biomes may be similar or closely related, producing the observed patterns.

On the other hand, the evergreen equatorial rainforest and the tropical deciduous woodland exhibited a considerably low diversification rate given their high accumulation of specialist species ( Figure 7 ; Supplementary Table S3 ). These biomes provide ideal conditions for the diversification of reptiles in general due to their high temperatures, which accelerate metabolic rates and all associated processes (Head et al., 2009; Machac et al., 2012; Cadena, 2014; Jaramillo, 2020). Additionally, these are ancient biomes, with the evergreen equatorial rainforest appearing during the Late Cretaceous and the tropical deciduous woodland during the Early Eocene (Jaramillo et al., 2010; Jaramillo and Cárdenas, 2013; Carvalho et al., 2021). The high number of species present in these biomes today likely reflects significant historical diversification of specialists and possibly their climate niche conservatism associated with tropical forested biomes (Pelegrin et al., 2023). However, as species accumulate and niches become saturated, diversification has decreased. This explains why diversification rates are similar to, or in some cases lower than, those in areas with less favorable conditions but also lower species saturation. On the other hand, temperate biomes are relatively more recent. The youngest colonized biome is the sclerophyllous woodland and shrubland, which dates back to the Plio-Pleistocene (Buerki et al., 2012). This is followed by the taiga, originating in the late Miocene (DeVore and Pigg, 2013), and the steppe, which dates to the early Miocene (Wang, 2004). In other words, these biomes emerged after the major adaptive radiation event experienced by turtles during the Eocene (Jaramillo, 2020; Thomson et al., 2021). This, along with the saturation of previously mentioned tropical biomes and the reduction in diversification rates within the group after this period, may provide an additional explanation as to why no significant differences were found in the diversification rate (DR) of turtles across the various biomes they inhabit.

The PGLS test revealed a significant negative correlation between DR and BSI, indicating that species with lower BSIs tend to have higher DRs ( Supplementary Table S9 ; Figure 6 ), which aligns with the resource use hypothesis and is consistent with previous studies conducted on other vertebrate and insect groups (Pelegrin, 2016; Gamboa et al., 2022; Hernández Fernández et al., 2022). Finally, it is important to consider that due to the focus of this study, our results are restricted to phylogenetic data and do not incorporate fossil evidence. As a result, we do not delve into extinction events or their impact on the evolutionary history of turtles.

In conclusion, this study provides compelling evidence that the macroevolutionary trends of turtles align with the primary tenets of the resource use hypothesis. Particularly, the greater diversification and accumulation of biome specialist species within the lineage. However, they exhibit peculiarities that have diverged from previous studies in other groups. They did not show distinguishable diversification rates among biomes, nor did they demonstrate greater species accumulation in climatically extreme biomes. This might be associated with higher levels of biome conservatism associated to their small home ranges, ectothermy and low morphological variability, with the ecophysiological limitations of this group for the survival to freezing temperatures, or a combination of both factors. Ultimately, some lineages have demonstrated a remarkable ability to adapt and diversify in biomes with highly variable climatic conditions. These insights shed light on the complex interplay between turtles and their environment, highlighting the importance of considering both historical and ecological factors in understanding their evolutionary history.