Global determinants and conservation of evolutionary and geographic rarity in land vertebrates.

Deciphering global trends in phylogenetic endemism is crucial for understanding broad-scale evolutionary patterns and the conservation of key elements of biodiversity. However, knowledge to date on global phylogenetic endemism and its determinants has been lacking. Here, we conduct the first global analysis of phylogenetic endemism patterns of land vertebrates (>30,000 species), their environmental correlates, and threats. We found that low temperature seasonality and high topographic heterogeneity were the main global determinants of phylogenetic endemism. While phylogenetic endemism hotspots cover 22% of Earth, these regions currently have a high human footprint, low natural land cover, minimal protection, and will be greatly affected by climate change. Evolutionarily unique, narrow-range species are crucial for sustaining biodiversity in the face of environmental change. Our global study advances the current understanding of this imperilled yet previously overlooked facet of biodiversity.

INTRODUCTION

Human-induced biodiversity declines are increasing and affect all corners of Earth (–). Current species extinction rates are estimated to be two to three orders of magnitude higher than “background” extinction rates (). Limited conservation resources make identifying key areas of biodiversity and their potential threats a paramount research and conservation agenda (–). Much scientific attention is directed toward rare, unique, and at-risk species (). Such species have irreplaceable ecological and evolutionary characteristics and a grim future outlook (, , ). Consequently, understanding the global geographical distribution of range-restricted species (endemics) and the mechanisms that create and maintain endemism centers is vital as these locations hold facets of biodiversity not represented elsewhere on the planet (, , ). Moreover, endemics’ spatial restriction makes them especially vulnerable to extinction due to different causes such as hunting, invasive species, climate change, and land-use change (, , ). Most historical tetrapod extinctions were of endemic species, e.g., (, ).

The concept of endemism has recently been placed within an evolutionary perspective in what is known as phylogenetic endemism (, ). Phylogenetic endemism is a spatial measure of relative phylogenetic diversity in an assemblage weighted by range size (). It is usually measured as the sum of branch length of descendant species in an assemblage weighted by the descendant range sizes (, ). Phylogenetic endemism thus extends the concept of phylogenetic diversity to highlight regions containing the most geographically rare and evolutionarily unique lineages (, ). Evolutionarily unique species are likely to be associated with unique ecological, morphological, behavioral, and functional attributes (, ). However, to-date global explorations of endemism (, , ) have yet to incorporate phylogenetic aspects (). Further, areas with high species richness or endemism are not always congruent with high phylogenetic endemism (, ). Therefore, phylogenetic endemism can be an important complementary metric of biodiversity in conservation prioritization (–).

Regions of high phylogenetic endemism are also potentially key to understand evolutionary mechanisms underpinning biodiversity patterns. For example, it can suggest what makes some species and clades radiate and spread geographically, while others remain confined in space over geological time (, ). Isolation, for example, on islands, promotes geographic rarity (, ) and consequently phylogenetic endemism. Furthermore, some environmental factors can facilitate the origination of phylogenetic endemism, while other factors could promote its persistence (see Table 1 for a list of main hypotheses of phylogenetic endemism determinants). For instance, spatial heterogeneity in topography, or habitat complexity (Table 1, hypothesis 1), offers ecological niche opportunities that promote diversification of geographically restricted lineages (–). Moreover, long-term climatic stability has been shown to aid the maintenance of geographic rarity (Table 1, hypothesis 2). Specifically, long-term climatically stable regions act as “museums” for endemic species as regions with unstable climates in the past lost many of their small-ranged species (, ) or select for species with high dispersal abilities (to enable tracking climates) (). Therefore, identifying regions where endemic species have been maintained for eons (i.e., paleoendemism) or have recently originated (i.e., neoendemism) could improve our understanding of the evolution of geographic rarity (, , ). Nevertheless, explicit tests of these hypotheses at the global level and how they differ between taxonomic groups have been missing.

Table 1.

Hypotheses and related predictors of (phylogenetic) endemism in tetrapods.

#	Factor	Variables	Hypothesis
1	Heterogeneity/barriers	Mean elevation range,number of soil types, andnumber of ecoregions	Heterogeneity, e.g., topography (28),generates dispersal barriers orpromote specialization, creatingsmaller species ranges (28–30)
2	Climatic stability	Temperature and precipitationvelocity since theLast Glacial Maximum	Long-term, less stable climate selectsfor species with high dispersalabilities (to enable them to trackchanges rapidly), reducingendemism (31) or species withsmall ranges in climaticallyunstable regions may face highextinction rate (31)
3	Climatic seasonality	Temperature andprecipitation seasonality	Seasonal variation in climaticconditions selects for species withbroader climatic tolerances andlarge ranges—Rapoport’s rule andthe climatic variability hypothesis (42).
4	Climatic uniqueness	Principal components of 19bioclimatic variables	Regions with distinct climaticconditions than those of thesurrounding area select for localadaption (i.e., specialization) andlimits range expansion (15).
5	Energy availability	Net primary productivity	High-energy availability may increasespecialization and opportunity forspeciation promoting endemism(20, 34) but simultaneously canreduce phylogenetic diversity byallowing the coexistence of closelyrelated species (49).

Here, we took advantage of the recently available, almost-complete distributions and phylogenies of terrestrial tetrapods (30,123 species) and addressed the following questions: (i) Where are the evolutionarily unique narrow-ranged species located, and how congruent are these regions for different classes of tetrapods? (ii) What are the main determinants of phylogenetic endemism? Lastly, (iii) how well are hotspots of phylogenetic endemism protected, and what are the anthropogenic threats they are facing now and are projected to face in the future?

RESULTS AND DISCUSSION

Global patterns of phylogenetic endemism

Using species distribution maps and molecular phylogenies, we mapped global phylogenetic endemism patterns for terrestrial tetrapods at a resolution of 96.5 × 96.5 km equal-area grid cells. We found a strong latitudinal pattern of phylogenetic endemism in terrestrial tetrapods, with regions of high phylogenetic endemism concentrated in the tropics and in temperate regions of the Southern Hemisphere (Fig. 1 and fig. S1). This latitudinal pattern of phylogenetic endemism arises as a result of a strong correlation between phylogenetic endemism and species richness (fig. S2). After correcting for species richness by taking residuals from regressions between per-grid-cell richness and phylogenetic endemism, we found a nonlinear relationship between phylogenetic endemism and latitude (Figs. 1 and 2; henceforth, we refer to phylogenetic endemism to those values corrected for species richness). Phylogenetic endemism is relatively low at high latitudes in the Northern Hemisphere and generally increases southward, peaking at high latitudes in the Southern Hemisphere (Figs. 1 and 2). Regions of high phylogenetic endemism with low species richness are predominantly located on islands and along coastal or peninsular regions of continents. They are often found in regions bordering the most species-rich areas (Figs. 1 and 2).

Fig. 1.

Global patterns of phylogenetic endemism and congruence of phylogenetic endemism hotspots for tetrapods.

(A) Spatial and (B) latitudinal patterns of phylogenetic endemism (uncorrected for richness; red). (C) Spatial and (D) latitudinal patterns of phylogenetic endemism (corrected for richness; blue). Note that the values of phylogenetic endemism uncorrected for richness are log-transformed in graph (B), and the fitted line is from loess regression. (E) Congruence map of phylogenetic endemism hotspots. Hotspots are selected according to randomization analysis [categorical analysis of neo- and paleoendemism (CANAPE)] conducted separately for each tetrapod class. Hotspots are colored according to whether cells are covered by a single class (blue, mammals; red, birds; purple, amphibians; and yellow, squamate reptiles) or more than one class (peach, any two; brown, any three; and black, all four class). The inset doughnut graph represents the percentage of unique significant hotspot cells for each tetrapod class and the percentage of significant hotspot cells shared by more than one tetrapod class. See fig. S4 for a detailed exploration of the percentage of overlap between classes.

Fig. 2.

Global patterns of phylogenetic endemism for tetrapod classes.

Global patterns of phylogenetic endemism uncorrected for richness (A, D, G, and J), richness corrected phylogenetic endemism (B, E, H, and K), and their latitudinal patterns for mammals (A to C), birds (D to F), amphibians (G to I), and squamate reptiles (J to L). Note that the values of phylogenetic endemism uncorrected for richness are log-transformed in graph (C,F,I and L), and the fitted line is from loess regression.

The montane speciation model proposes that diversification centers in the mountains frequently “feed” species to the surrounding lowlands (, ). Our findings that phylogenetic endemism is especially high in coastal areas adjacent to mountain systems (e.g., west of the Andes in Neotropics, in the Western Ghats in India, and east of the Great Dividing Range in Australia) are consistent with this model, suggesting that the mountains may restrict the ranges of these lowland coastal species. Furthermore, species range size is constrained by the available land area (). Therefore, the high concentration of phylogenetic endemism at southern latitudes (Figs. 1 and 2) can be explained by the general decline in the continental area toward the southern tips of the southern continents. We examined this pattern by developing a metric that combined grid cell proximity to coastlines and mountains and compared it to phylogenetic endemism (see the Supplementary Materials). We found a positive, albeit weak, association between phylogenetic endemism and grid cells’ proximity to both coastlines and mountains across all classes (fig. S3).

Phylogenetic endemism hotspot congruence across tetrapods

To identify hotspots of phylogenetic endemism for each taxon, we used a randomization procedure developed by Mishler et al. (). This method identifies hotspots of phylogenetic endemism under appropriate null expectations by accounting for the strong phylogenetic endemism—species richness relationship (, ). Hotspots are designated if the empirical phylogenetic endemism value was higher than 975 random values in the simulations (see Materials and Methods for more details). Of 17,080 terrestrial grid cells, 3797 (22%) were identified as hotspots of endemism for at least one of the four tetrapod classes (Fig. 1E and fig. S4). Almost half of these grid cell hotspots (49.1% of 3797 grid cells, i.e., 1867 grid cells) belong to a single class, and the rest are shared by 2 to 4. We found that phylogenetic endemism hotspots of mammals and birds overlap less (i.e., have more unique grid cells: 770 and 639 cells or 20.3 and 16.8% of the total 3797 grid cells for mammal and birds, respectively) with hotspots for all the classes combined, compared to amphibians (2.8% or 108 grid cells) and squamates (reptiles henceforth; 9.2% or 350 grid cells; Fig. 1E).

In total, 9.1% of the phylogenetic endemism hotspots (344 of 3797 cells) are shared among all four tetrapod classes (Fig. 1E). This value is 16 times higher than the 0.057% overlap expected if phylogenetic endemism hotspots of different classes were distributed randomly (averaged from 1 million simulations; see Materials and Methods). Regions of high phylogenetic endemism congruence (i.e., where phylogenetic endemism is high in several groups) are found in the Caribbean islands, southwestern Mexico, along the Andes, the Cape region in Africa, eastern Madagascar, Sri Lanka and the southern Western Ghats, Vietnam, New Guinea, Australia east of the Great Dividing Range, the Philippines, and Sarawak (Fig. 1E).

While phylogenetic endemism hotspots are clustered in space, they lack congruence across tetrapod classes. This suggests that processes underlying the generation and maintenance of phylogenetic endemism centers are potentially different across tetrapod classes (see the next section). Nevertheless, cross-taxon congruence and endemism are scale-dependent (, ). When repeating our analyses at the ecoregion scale (847 ecoregions globally), we found higher cross-taxon congruence (fig. S5). Conversely, we expect lower cross-taxa congruence when analyzed at finer resolutions (, , ). Consequently, when setting conservation actions, a scale explicit framework should be used (, , ).

We then explored which of the taxa contributes more to the overall phylogenetic endemism of tetrapods. We used a global tetrapod tree to compare the top 5% phylogenetic endemism (corrected for richness) grid cells of all tetrapods with the top 5% phylogenetic endemism (corrected for richness) grid cells of each taxon. We found that reptiles contributed relatively more to tetrapod phylogenetic endemism than other classes (at least 1.5 times higher than other classes; fig. S6). This confirms previous research highlighting the importance of reptile evolutionary history to global phylogenetic endemism ().

Spatial determinants of phylogenetic endemism

We tested our hypotheses relating to the effects of the different geoclimatic predictors on phylogenetic endemism (Table 1, hypotheses 1 to 5) or their combination simultaneously using regression models. We conducted tests on phylogenetic endemism values corrected for species richness separately for (i) mainland regions, (ii) island (landmasses smaller than Australia), and (iii) all regions with insular-endemism (yes/no) as a predictor. As we were predominantly interested in drivers of endemism in mainland regions, we present here the results of the first test and report the results of the other analyses in the Supplementary Materials (referring to them when they differ from the mainland only tests). As we found significant (α = 0.05) spatial autocorrelation in the residuals of the fitted ordinary least squares linear models (table S1), we conducted a spatial simultaneous autoregressive (SAR) analysis, which we present here. We consistently found that models that included more than one group of predictors—heterogeneity, climatic stability (since the Last Glacial Maximum), climate seasonality, climatic uniqueness, and available energy; see Materials and Methods and Table 1—perform better than models with a single predictor (tables S2 to S6). Consequently, we report the results of model-averaged multivariate models (at α = 0.05). Top-ranked models mostly included temperature seasonality, elevation range, and net primary productivity (tables S2 to S6). In models including islands as an additional predictor, insularity was always the second most important factor (fig. S7).

Current temperature seasonality (which is associated strongly with latitude) was the strongest predictor of phylogenetic endemism, which peaked in a-seasonal regions across all classes (especially birds; Fig. 3). Species inhabiting high latitudes (where temperature seasonality is high) are thought to have large ranges (Rapoport’s rule), potentially due to their adaptation to seasonal climates (the climate variability hypothesis 3) () and therefore lower phylogenetic endemism (). While we did find high phylogenetic endemism in high southern latitudes (which contradicts the overall trends), (i) these regions have relatively smaller areas than northern high latitudes, and (ii) the overall seasonality in high southern latitudes is lower than in the northern high latitudes. The prominence of seasonality as a predictor globally was also maintained when islands were included as a predictor and was the most important predictor also for islands separately (fig. S7). For amphibians, we found that climatic uniqueness (hypothesis 4) had a stronger positive effect on phylogenetic endemism than for other groups (Fig. 3). These results are consistent with the unique physiological requirements of amphibians—especially their sensitivity to low precipitation ().

Fig. 3.

Determinants of global phylogenetic endemism in terrestrial tetrapods.

Estimates (standardized slopes) and SEs were obtained from model-averaged SAR analysis. Nonsignificant estimates (P > 0.05) are colored red. Results from nonspatial linear regression analyses are presented in fig. S8. LGM, Last Glacial Maximum.

Elevation range (a proxy of habitat heterogeneity; Table 1, hypothesis 1) was the second strongest predictor and significantly positively correlated with phylogenetic endemism for amphibians, mammals, and reptiles (standardized coefficient: mammals, 0.04; birds, <0.001; amphibians, 0.23; reptiles, 0.08; Fig. 3). This is consistent with previous research regarding the role of elevation in promoting endemism (, ). The relative importance of topographic heterogeneity and temperature seasonality differed across classes and may stem from the lower dispersal ability of these taxa. Large homeotherms with high dispersal ability are expected to be less affected by local barriers (heterogeneity) (). In contrast, amphibians and reptiles have, on average, smaller body sizes and distributions () and lower dispersive abilities and consequently could potentially be affected by factors that vary at a finer spatial resolution (, ). Models predicting phylogenetic endemism across zoogeographic regions (fig. S10) consistently identified elevation range as the strongest predictor of phylogenetic endemism for reptiles and amphibians (compared to other predictors). No such trends were found in birds. For mammals, elevation range was only an important predictor in the Neotropic and Afrotropic realms (fig. S10).

Net primary productivity was positively associated with phylogenetic endemism uncorrected for species richness for all classes combined and each class separately (fig. S9). This mirrors trends regarding net primary productivity’s association with species richness (, ). After correcting for richness, net primary productivity had a negative effect on phylogenetic endemism for all classes combined and for birds (Fig. 3). However, phylogenetic endemism corrected for richness was positively associated with net primary productivity for reptiles and was nonsignificant for mammals and amphibians (Fig. 3). In general, the negative association between net primary productivity and phylogenetic endemism supports recent findings that highly productive sites potentially promote the coexistence of closely related species due to denser niche packing (). For reptiles, however, the occurrence of many large ranged species in less productive regions, e.g., deserts (), relative to the prevalence of range-restricted species in highly productive regions (e.g., tropical forests) can explain the positive trend.

When we included grid cell island status as a predictor, we found that the links between phylogenetic endemism and net primary productivity are unchanged (fig. S7A). In island-only grid cells, we found that phylogenetic endemism was negatively associated with net primary productivity for all classes combined and for mammals (nonsignificant predictor for the other groups; fig. S7B). We further found net primary productivity to be negatively correlated with phylogenetic endemism for all classes combined and for each class (including reptiles) in a model not accounting for spatial autocorrelation (fig. S8). Hence, there are small regions of high reptile phylogenetic endemism that have high net primary productivity and also others with low phylogenetic endemism and low net primary productivity, while large clusters of reptile phylogenetic endemism show the opposite trend. This may be attributed to the more local and patchy nature of reptile distributions compared with mammals and birds ().

We found that phylogenetic endemism was consistently higher in regions that experienced lower temperature velocity since the Last Glacial Maxima across all classes but only in a model that did not account for spatial autocorrelation (fig. S8). After correcting for spatial autocorrelation, temperature velocity was not significantly associated with global phylogenetic endemism; see also (). However, temperature velocity was significant in the Australasian, Palearctic, and Nearctic realms for all classes and in the analysis of island-only cells (figs. S7B and S10). The role of paleoclimate in shaping current biodiversity patterns has recently gained prominence (), but the relative role of historical versus current climate remains debated (, ). Our results highlight the dominance of current climates in shaping phylogenetic endemism patterns, with past climate probably only acting at the regional level (, , ). The influence of paleoclimate on phylogenetic endemism can be explained by two main mechanisms. Climatically unstable areas may select for species with high dispersal ability (and hence large ranges) that are better at tracking changing climates (). Therefore, climatically stable regions may have acted as refugia for narrow ranged species (). A related hypothesis postulates that new species seldom evolve in unstable climatic regions due to limited opportunities for specialization (). We found some support for the former hypothesis in mammals and reptiles where paleoendemism hotspots have lower temperature velocity since the Last Glacial Maximum than neoendemism hotspots (fig. S11). Beyond these, most other predictors had little effect on phylogenetic endemism across taxa (Fig. 3).

Hotspots of paleo- and neoendemism and the importance of topographic heterogeneity

Regions containing many endemic species can be divided according to species evolutionary age (neo- versus paleoendemism) (, ). Paleoendemics are species that diverged long ago (compared to the rest of their clade) and have small ranges either because their ranges contracted or because they never expanded from their source of origin (). Hence, paleoendemics have long branch lengths and high phylogenetic endemism (). Conversely, neoendemics are species that originated recently (i.e., have short branch length) and are (still) confined to narrow ranges. To understand spatial patterns and drivers of neo and paleo endemism, we used a categorical analysis [categorical analysis of neo- and paleoendemism (CANAPE)] (). We found that most hotspots of phylogenetic endemism (identified using randomizations; see Materials and Methods) contained lineages of both recent and ancient origins (mixed endemism; i.e., hotspots that cannot be defined as either paleo- or neoendemism). These accounted for 86% of hotspot cells for mammals, 76% for birds, 96% for amphibians, and 87% for reptiles and are predominantly concentrated in the tropics (Fig. 4, A to D, and fig. S12A). High species diversity in the tropics has been proposed to arise either because the tropics are a “cradle” (due to the high rates of speciation) or a museum of diversity (due to low rates of extinction) (). Our findings suggest that the tropics serve as both (), with most mixed endemism hotspots occurring in or around tropical mountain regions and which are associated with both high climatic stability and topographic heterogeneity (fig. S12B).

Fig. 4.

Global distribution of endemism types and importance of topographic heterogeneity.

Centers of neoendemism (Neo), paleoendemism (Paleo), mixed endemism (Mixed), or not significant regions (Not Sig) for (A) mammals, (B) birds, (C) amphibians, and (D) squamate reptiles from the CANAPE analysis. (E to H) Boxplot showing differences in elevation range calculate per grid cell for endemism types across the taxon class. Results from Fisher-Pitman permutation test (10,000 permutations; α = 0.05 after “holm” correction) are presented only for comparison between neo- versus paleoendemism hotspots. Other comparisons are presented in figs. S11 and S12. ***P < 0.001; **P < 0.01.

We did not find overarching patterns for either neo- or paleoendemic hotspots across all classes (Fig. 4, A to D). We found neoendemic hotspots in the southern Andes for reptiles; Central America for amphibians and mammals; Patagonia for mammals; and the Andes, the Hengduan, and Asir Mountains for birds. Paleoendemics were predominantly found in Morocco and Oman for reptiles; Namibia, Hispaniola, Somalia, and south-west Australia for mammals; and Madagascar, Papua, Borneo, south-east Australia, New Zealand, and northern Europe for birds. For amphibians, no such distinct paleoendemic hotspots could be identified (Fig. 4C). These results suggest that different mechanisms give rise to trends in different taxa. In accordance with this, an in-depth evaluation of factors that shape neo- versus paleoendemism hotspots did not reveal a consistent determinant of endemism types among tetrapods (except for elevation, see below) in either univariate or multivariate analyses (fig. S11 and S12B).

Nevertheless, we found that the mean of grid cell elevation range in neoendemism hotspots was significantly higher than in paleoendemic hotspots (all classes in univariate analysis and reptiles and birds in the multivariate analysis; Fig. 4, E to H, and fig. S11). Therefore, these results support the notion that mountains have greater importance as cradles of diversity (, –), as also shown recently for Neotropical orchids (), and in simulations of South-American biodiversity highlighting the importance of the Andean slopes as cradles (). However, the mechanisms underlying the generation of neoendemism hotspots across taxa can be idiosyncratic. For mammals and reptiles, hotspots in Patagonia and the southern Andes may have arisen from recent local radiations of few lineages: subterranean rodents from the genus Ctenomys in mammals () and lizards from the genera Liolaemus and Phymaturus in reptiles (). For birds living in the tropical Andes, local radiations in several lineages may have also been promoted by the mountainous terrain (). However, it is possible that neoendemic birds in the Asir Mountains may have evolved in this region when it became a climatic refugium during the Pleistocene. In general, these results highlight the importance of mountains in shaping biodiversity patterns (, ).

Assessment of protection and threats to phylogenetic endemism hotspots

Protected areas are a key tool to tackle the biodiversity crisis (). To understand the current protection status of tetrapod phylogenetic endemism, we calculated the percentage overlap of the global network of protected areas (protectedplanet.net) with each grid cells. We then tested whether this overlap percentage significantly differed between grid cells designated as phylogenetic endemism hotspots and nonhotspots (see above). This was analyzed either for all protected area categories or just those designated for biodiversity conservation—categories I to IV that we present in the Supplementary Materials (fig. S13). The median per-grid-cell overlap of protected areas for all classes combined is higher for hotspot grid cells than nonhotspot cells (2.74% versus 0.48%; Fig. 5A). However, for mammals, amphibians, and reptiles, the protection level of hotspots was similar to that of nonhotspot cells (Fig. 5B). The Neotropics stood out as having lower protection for their phylogenetically endemic hotspots (compared to nonhotspots; Fig. 5A). Moreover, in the Afrotropics and Australasia, hotspots have the same level of protection as nonhotspots (Fig. 5A). The median protected area coverage of phylogenetic endemism hotspots per se on a per-grid-cell basis is low, 2.7% at the global level and less than 10% for all taxon classes (Fig. 5, A and B), the minimum threshold target set to safeguard biodiversity (). This low per-grid coverage may stem from the fact that protected areas conventionally target flagship species or species richness and neglect regions where endemic species are concentrated, e.g., islands and mountains (, ). Hotspot regions that require protection mainly occur along the southern Andes, Horn of Africa, Southern Africa, and the Solomon Islands (Fig. 5A). We further compared the overall overlap of protected areas for either all hotspots or all nonhotspot grid cells combined together. We found that globally hotspots are better covered than nonhotspots (11.36% versus 10.67%, respectively). Overall, our study highlights the urgent need to incorporate species evolutionary history in the current conservation priority setting (, , ). Our analyses are carried out at a coarse resolution—the most appropriate at a global scale intended to highlight global patterns (). Nevertheless, most conservation actions are conducted at smaller scales, and actual phylogenetic endemism hotspots are probably smaller than one grid cell.

Fig. 5.

Assessment of current and future threats to phylogenetic endemism hotspots (i.e., either neo-, paleo-, or mixed endemism) across all tetrapod classes.

(A and B) The proportion of overlap with protected areas [all protected area (PA) categories] comparing hotspot (H) and nonhotspot (NH) areas, (C) human footprint index, (D) number of threatened species as classified by the International Union for Conservation of Nature (IUCN) per grid, (E) the difference in the proportion of natural landcover (both primary and secondary forest) between current (2015) and future (2100) time periods, climate change vulnerability assed as the ratio of future (2070) to past velocity (Last Glacial Maximum) of climate change for (F) temperature and (G) precipitation. For estimates involving future threats, the Representative Concentration Pathway (RCP) 4.5 scenario was used. In the boxes, white horizontal lines represent medians; lower and upper bounds of the box correspond to first and third quartiles, respectively. Significance codes indicate differences between the hotspot and nonhotspot cells from the two sample Fisher-Pitman permutation test (10,000 permutations): ***P < 0.001; **P < 0.01; *P < 0.05; NS, not significant.

As phylogenetic endemism hotspots are less covered by protected areas, we further explored whether they are also threatened by human modification and face high species extinction risk. For this, we explored the overlap of phylogenetic endemism hotspots with the human footprint index to quantify the exposure of hotspots to human disturbances. Hotspots are currently more affected by humans than nonhotspot regions (Fisher-Pitman permutation test; P < 0.001 for all taxa; Fig. 5C). We further found that phylogenetic endemism hotspots hosted more species classified by the International Union for Conservation of Nature (IUCN) as threatened [vulnerable (VU), endangered (EN), and critically endangered (CE)] than nonhotspots areas (Fisher-Pitman permutation test; P < 0.001 for all taxa; Fig. 5D). This is expected because Red List assessments rely strongly on small-range size as a determinant of threat (iucn.org). Together, these results suggest that the evolutionarily unique range-restricted taxa are under high extinction risk and face high impact by human modifications ().

We then used future projections of land use and climate to understand the potential threats to phylogenetic endemism hotspots in the future. We explored the overlap of hotspots with natural (primary and secondary forests, savannahs, deserts, etc.) versus modified land cover classes, both currently and in the future, and compared these proportions to nonhotspot regions. We found that, currently, hotspots have a lower proportion of natural land cover than nonhotspot areas (for mammals, birds, and reptiles; fig. S14). When exploring the trends of losses and gains to natural land cover between 2015 and 2100 [under a moderate greenhouse gas mitigation scenario, Representative Concentration Pathway (RCP) 4.5; results under other scenarios are presented in fig. S15], we found that hotspots are not projected to gain more natural land cover by 2100 (mammals, birds, and reptiles; Fisher-Pitman permutation test; P > 0.05; Fig. 5E). We compared the ratio of future to the past velocity of climate (2013 to 2070 and the Last Glacial Maximum to the beginning of the 21st century, respectively) across hotspot and nonhotspot areas. For all tetrapod classes, this ratio was significantly higher in hotspots than in nonhotspot regions for temperature velocity (Fisher-Pitman permutation test; P < 0.001; Fig. 5F). For precipitation, the ratio of future to past velocity was significantly lower in hotspots compared to nonhotspot regions for all classes, except for amphibians (Fisher-Pitman permutation test; P < 0.001; Fig. 5G).

High spatial variation in climate within hotspots (due to mountains) is thought to buffer them from the impact of climate change (, ). However, our results suggest that hotspots might be facing faster changes in temperature than surrounding areas. These results support recent findings that despite these hotspots served as a refugium due to long-term climatic stability in the past, they are most vulnerable to human-induced climate change (, ). Moreover, many of the phylogenetic endemism hotspots that we identified are in islands and coastal areas, and dispersal of species to climatically suitable regions may be limited, particularly in fragmented landscapes dominated by humans (). Overall, these results suggest that the implementation of policies to curtail greenhouse gas emissions () and establishing conservation measures that also account for global change is paramount (). The latter can be achieved by incorporating connectivity of future habitats in protected area establishment design () or translocation of species to regions with suitable future climates () when this is a viable option.

Summary

Understanding where and how biodiversity originated, and is being maintained, is one of the long-standing endeavors of biogeographers, evolutionary biologists, and conservationists. Species found nowhere else that are also evolutionarily distinct are of particular fascination and represent a crucial yet overlooked facet of biodiversity. Our study integrates patterns, mechanisms, and conservation aspects of phylogenetic endemism for more than 30,000 land vertebrates globally. Above and beyond general richness patterns, global phylogenetic endemism is predominantly correlated with current climatic seasonality and topographic heterogeneity across most tetrapod classes. Our findings highlight that there are many paths to become a phylogenetic endemism hotspot that is idiosyncratic to geographic location, taxonomic class, and the evolutionary age of taxa. We emphasize the importance of including phylogenetic endemism in conservation decision-making and provide several conservation implications of incorporating this new dimension. Our global assessment revealed that most of the phylogenetic endemism hotspots have minimum protection and face a greater risk from both current threats and projected climate change. In light of these findings, proper action to tackle the ongoing environmental change needs to be taken to ensure the persistence of thousands of narrow-ranged species that represent millions of years of unique evolutionary history.

MATERIALS AND METHODS

Species distribution and phylogenetic data

We obtained vector range maps of all mammals and amphibians available from the IUCN (iucn.org; version 6.2), from BirdLife for birds (version 3.0) (), and an updated version of squamates distribution (GARD internal version 1.5) (). Resident and native distribution ranges of species were gridded using letsR package () onto grid cells of dimension 96.5 km by 96.5 km (~1° at 30°N/S) with a Behrmann equal area projection (). This resolution has been suggested to be less prone to incur the false presence of species at the global level, which was our scale of exploration (). When assigning hotspots for regional or local conservation initiatives, analyses at finer spatial resolution should be attempted (). We used near-complete phylogeny of tetrapods from published sources, mammals (), birds (), amphibians (), and squamates (). We matched phylogeny and distribution data using online taxonomic databases accounting for synonyms (taxonomic list provided in the supporting dataset). We based all our analysis on 100 randomly drawn trees from the posterior distribution of phylogenetic trees and summarize the results with median values. We grafted the species-level trees of the four classes onto the backbone tree of Tetrapoda obtained from the TimeTree website (www.timetree.org/) and scaled them to match the time of divergence for the major classes using the backbone tree. Since we are interested in hotspots of phylogenetic endemism and their relevance to conservation, we omitted grid cells with few species overall (<5 species, e.g., small islands). The final dataset contained 5370 mammals, 9206 birds, 6335 amphibians, and 9212 reptile species, which accounted for 84% of currently recognized tetrapods.

Data analyses

All analyses were conducted either in Rstudio () or ArcMap (version 10.7.1) (). We accessed the BIODIVERSE software (version 3.0.0) () in R using the “Biodiverse pipeline” (https://github.com/NunzioKnerr/biodiverse_pipeline) to estimate phylogenetic endemism.

Phylogenetic endemism

Phylogenetic endemism was calculated using the formula below based on the method described by Rosauer et al. ()where C is the set of branches representing the minimum spanning path joining all taxa in a grid cell to the root of the phylogenetic tree; L is the length of a branch, c, from the spanning path C; and R is the global range size of all descendant species of branch c (as the union of combined ranges), measured as the total number of grid cells that they occupy throughout their range without repeated count of grid ().

Significance test for phylogenetic endemism

Since regions with high species richness are likely to hold more endemic species than expected by chance (), we performed randomization tests to identify grid cells that have significantly higher or lower observed phylogenetic endemism than expected given the richness and range size of species per grid cell (). For this, we retained the same number of species per grid cell and range size of each taxon and randomly assigned species to each grid cell based solely on the number of species it contains. To maintain the continuity of range for each species in randomization trials, we used the “rand_struct” option in BIODIVERSE (). We performed 999 such randomization iterations and compared empirical values to random values. We designated cells as containing significantly high phylogenetic endemism if the empirical value was higher than 975 random values (). Likewise, cells were classified as significantly very high (> 99%, i.e., in the top 10), significantly low (<2.5%), or significantly very low (<1%) phylogenetic endemism [see () for more details]. All other cells were considered as having phylogenetic endemism not significantly different than expected by chance. These cell categorizations were further used to identify significant hotspots of neo-, paleo-, and mixed endemism hotspots (below).

Categorical analysis of neo- and paleoendemism

We performed CANAPE analysis in the BIODIVERSE software to define the type of endemism of a given cell (), based on its evolutionary age. We defined “new endemism hotspots” (“neoendemism”) as cells with shorter branch lengths than predicted by chance (i.e., where endemism is due to recent radiations), “old-endemic hotspots” (“paleoendemism,” i.e., more species with long branch lengths than expected by chance), or a mixture of both—mixed endemism (). For this, we estimated relative phylogenetic endemism. Relative phylogenetic endemism is the ratio of phylogenetic endemism calculated using the original tree (PEoriginal) divided by the phylogenetic endemism measured on a tree whose branch length of all species in the entire phylogeny is made equal but retaining the topology and number of species (PEequal) (). Therefore, a significance test of measured relative phylogenetic endemism will tell us whether the branches of all species present in each cell are disproportionately longer or shorter than expected

The significance of relative phylogenetic endemism values was tested using randomizations as described above. We designated cells with significantly high relative phylogenetic endemism as characterized by paleoendemism if their phylogenetic endemism is higher than PEequal in the 97.5% of randomizations and as hotspots of neoendemism if PEoriginal is lower than PEequal in the 97.5% of the randomizations. Cells were classified as mixed endemism if both PEoriginal and PEequal are significantly high (97.5% of the randomizations; mixture of paleo-, neo-, and nonendemic), but relative phylogenetic endemism is not significantly high or low (). Last, if neither PEoriginal nor PEequal is significantly high, then cells were classified as not being endemism hotspots. We further simulated the overlap of phylogenetic endemism hotspots of all tetrapod classes based on null expectation. For this, we randomly sampled the same number of hotspot cells for each taxon and calculated the overlap of these cells by chance. We repeated the randomization process 1,000,000 times and calculated the mean number of overlapping cells, which we then compared to our empirical data.

Geoclimatic predictors of phylogenetic endemism

We compiled a set of nine geographic and climatic predictors that have been shown to influence endemism (Table 1). For all analyses, we obtained predictor variables in their finest available spatial resolution and calculated their mean values to match our grid’s resolution. We used BIO4 and BIO15 variables, which represent seasonality in temperature and precipitation, respectively, from the Climatologies at high resolution for the earth’s land surface areas (CHELSA) dataset at a resolution of 30 arc sec averaged for the period 1979–2013 (). Next, we measured climatic uniqueness as the distance in climatic space between each focal cell and its neighboring cells. This was done by computing the principal components (PC) analysis across all 19 bioclimatic variables obtained from the CHELSA dataset () and using the first two PCs (PC1 and PC2, which explained more than 75% of the variance). We then calculated the mean Euclidean distance between the focal cell and all the neighboring cells for PC1 and PC2. If the focal cell is climatically very different from the surrounding regions, then the estimated distance should be high for the focal cell. We estimated these values for cells within a radius of approximately 96.5 (3 × 3), 193 (5 × 5), and 289.5 (7 × 7) km from the focal cell. Results of this analysis were qualitatively similar irrespective of neighbor distance, and we present here results only for a neighbor distance of 96.5 km.

Briefly, we quantified climate change velocity as the spatial displacement in climate per unit time, i.e., the ratio of the temporal to spatial climate component following (). First, the temporal component was estimated as the absolute difference between the current and past climate during the Last Glacial Maximum (as degree Celsius per year). The spatial component of climate velocity, expressed as the displacement of climate over distance (as degree Celsius per kilometer), was calculated for the current climate using the slope function in package SDMTools (). We then divided the spatial component by the temporal component to obtain the velocity (kilometer per year). Before this, slope values equal to zero were adjusted to 0.001 to prevent dividing the temporal component by zero (). Both the current and the Last Glacial Maximum climate values were obtained from the CHELSA dataset (). Historic climate velocity measures were calculated for all available general circulation climate models separately (table S7). For each grid cell, we accounted for climate model uncertainty by using the median climate velocity estimated from all the seven models at a resolution of 30 arc sec. Velocity values were then averaged for each 1° grid cell separately for two climatic measures: mean annual temperature and mean annual precipitation.

We used the grid cell elevation range (maximum-minimum) as a proxy for spatial heterogeneity in altitude. We obtained elevation data from the STRM 4.1 at a resolution of 30 arc sec (). The number of unique categorical major soil types and ecoregions was estimated per grid cell using the harmonized soil dataset () and Terrestrial Ecosystem of the World database ecoregion map (), respectively. Net primary productivity values were obtained from the Socioeconomic Data and Applications Centre at a resolution of 0.25° × 0.25° ().

Assessment of threats to phylogenetic endemism hotspots

Protected area data were obtained from World Database on Protected Areas (WDPA) and was cleaned following WDPA guidelines (https://protectedplanet.net/, accessed February 2021). Protected area coverage was calculated as the percentage of overlap area with each grid cell (after correcting for available land area in each cell). Next, we obtained the human footprint index for the period 2009 at a resolution of 30 arc sec from Venter et al. (). The number of threatened species per grid cell was estimated as the sum of species classified as CR, EN, and VU in the grid cell for which IUCN assessment data were available ().

Future climate velocity values were estimated by using the same approach as described for historic climate velocity (i.e., since the Last Glacial Maximum). Mean annual temperature and precipitation data (i.e.., BIO1 and BIO12) for 2070 at a resolution of 30 arc sec were obtained from the CHELSA dataset from 32 general circulation climate models (table S7) (). Climate velocity for each grid cell was estimated as the median values for all the 32 climate models. To evaluate the vulnerability of phylogenetic endemism hotspots to future land-use change, we calculated the difference in the proportion of area occupied by natural vegetation (primary and secondary forests, savannahs, deserts, etc.) to modified land use between the years 2015 and 2100 for each grid cell. This value ranges between −1 and 1. A positive value indicates a loss of natural land cover, and a negative value indicates a gain of natural land cover in the future. The projected land use for the representative period 2015 and 2100 was obtained from Land-Use Harmonization database at a resolution of 0.25° × 0.25° (version LUH2 v2f) (). Both climate and land-use data were obtained for two emission levels, RCPs 4.5 and 8.5, to highlight “middle of the road” and “high-end” case scenarios of climate change.

Statistical analyses

For analysis involving determinants of phylogenetic endemism, most predictor variables had a moderate correlation among them (Spearman ρ < 0.7), and the variation inflation factors for all the predictors were lower than 5 in our linear models suggesting low multicollinearity. Therefore, we retained all nine predictor variables in the analyses. Phylogenetic endemism values were natural log-transformed to satisfy assumptions about the normality of residual error distribution. All predictor variables were scaled and centered to standardize the data (i.e., to estimate effect size). SAR error analysis (SARerror), which accounts for spatial autocorrelation, was performed using the spdep package in R (). Our SARerror models used a conservative row standardized weighing scheme “W.” Initial exploration of models that accounted for spatial autocorrelation within distances equivalent to 1° (96.5 km), 2°, 4°, 6°, 8°, and 10° at the equator revealed that a spatial error distance of 193 km (2°) best accounted for spatial structure in the analyzed data and, hence, was used throughout. We used five groups of models that included (i) each individual predictor (nine models), (ii) group of predictors [three models; climate seasonality (temperature and precipitation), climate velocity (temperature and precipitation), and heterogeneity (number of soils, number of ecoregions, and elevation range)], (iii) all nine predictors (one model), (iv) intermediate models included in backward stepwise elimination of nonsignificant terms (seven models), and (v) null model without any predictor (one model). We then used the model averaging procedure using the second-order Akaike information criterion (AICc) and Akaike weights to compute model-averaged coefficients and confidence intervals from models with ΔAICc < 7 relative to the best model (). Model averaging was performed using MuMIn package ().

We performed a two-sample Fisher-Pitman permutation test (10,000 permutations) using the coin package () to compare differences in current and future threats between hotspots and nonhotspot regions. We adjusted P values using the method holm to account for multiple comparisons.