1. Introduction

Species range sizes vary greatly (Rosenzweig, 1995). Range-restricted species are often endemic species, and vulnerable to extinction. Species that are confined to specific small areas are typically endemic and face a higher risk of extinction. Endemic species, which have limited geographical ranges within certain regions, often exhibit overlapping distributions, potentially indicating unique historical events in those areas (Mishler et al., 2014). Such species, when geographically concentrated, can establish centers of endemism. Endemism is a reflection of the processes of speciation, extinction, and dispersal shaping the distribution of species (Sandel et al., 2020). Consequently, elucidating the global patterns of endemic species and identifying the factors and mechanisms that lead to the formation of endemism centers are critical objectives in the field of biogeography (Qian et al., 2024a, 2024b).

Most prior research on endemism is taxon-centric, analyzing endemism based on species or genera. However, these taxon-based measures of endemism, known as taxonomic endemism, fail to encompass the evolutionary history and diversification patterns of clades over time (Rosauer et al., 2009; Mishler et al., 2014; Sandel et al., 2020). It is therefore crucial to employ metrics that reflect the evolutionary distinctiveness and phylogenetic relationships of species in endemism studies. Increasingly, studies are adopting a metric known as phylogenetic endemism, which considers the range sizes of clades rather than individual taxa (Mishler et al., 2014; Sandel et al., 2020; Cai et al., 2023; Qian et al., 2024a, 2024b). This approach is particularly valuable because regions with high phylogenetic endemism, often identified as centers of paleo-endemism, contain a greater number of evolutionarily unique lineages with limited geographic distributions compared to regions with low phylogenetic endemism (Cai et al., 2023). Consequently, incorporating evolutionary history and phylogenetic uniqueness into assessments of endemism allows for the identification of both paleo-endemism and neo-endemism centers (Cai et al., 2023; Qian et al., 2024a, 2024b).

Several factors are recognized as primary influences on endemism patterns, including present and historical climatic conditions, as well as the topographic and environmental variability (heterogeneity) within regions (Sandel et al., 2020; Qian et al., 2024a, 2024b). Areas with higher current temperatures and precipitation levels are observed to harbor a greater number of species with restricted ranges and elevated levels of endemism (Myers et al., 2000). Additionally, regions that have maintained a stable climate over long periods are likely to have acted as refuges for species with limited ranges during climatic shifts, thereby supporting higher endemism (Enquist et al., 2019). Topographic diversity also plays a dual role in fostering endemism. On one hand, it facilitates the emergence of neo-endemics by promoting diversification through habitat isolation and ecological specialization (Fjeldsa and Lovett, 1997). On the other hand, it aids in the preservation of paleo-endemics by enabling species to locally track suitable climatic conditions during periods of climate change (Carnaval and Moritz, 2008; Rosauer and Jetz, 2015). Consequently, regions with greater topographic variability tend to exhibit higher endemism (Sandel et al., 2020; Cai et al., 2023). Understanding the relative significance of these factors is essential for comprehending the geographic variations in endemism.

For plants at both global and regional scales, previous studies have investigated geographic variation of endemism, locations of the centers of paleo-endemism and neo-endemism, and the relationships between environmental factors and endemism for vascular plants. For example, at a global scale, Qian et al., 2024a, Qian et al., 2024b, and Cai et al. (2023) have explored geographic patterns and environmental correlates of endemism in angiosperms, ferns, and seed plants, respectively. However, for bryophytes, such studies are lacking. The present study aims to bridge this knowledge gap.

Liverworts, comprising approximately 7500 species (Laenen et al., 2018), represent one of the major groups of bryophytes. Collart et al. (2021) published a database with liverwort distributions worldwide. This study is structured around three main objectives: (1) to investigate the geographic patterns of taxonomic and phylogenetic endemism in liverworts in various regional floras globally; (2) to examine the correlations between liverwort endemism and environmental factors, including current climatic conditions, Quaternary climate changes, and variations in topography and environment; and (3) to pinpoint regions characterized by paleo-endemism and neo-endemism of liverworts.

2. Materials and methods 2.1. Operational geographic units and liverwort floras

The present study included 390 operational geographic units (hereafter regions) worldwide (Fig. S1). Small oceanic islands and Antarctica were excluded. Species lists of liverworts in these regions were obtained from Collart et al. (2021), which provided a species list for each of the regions used in our study. In total, 5971 species of liverworts were included.

2.2. Reconstruction of phylogenetic tree

In this study, a phylogenetic tree was constructed for the species analyzed. We utilized the time-calibrated phylogeny from Laenen et al. (2014), which establishes the earliest possible origin of the root node at 475 million years ago, as the backbone of our phylogenetic tree. The functions build.nodes.1 and Scenario 3 from the U.PhyloMaker package (Jin and Qian, 2023; see also Jin and Qian, 2022) were employed to integrate species from our dataset into this genus-level phylogenetic backbone. Prior research, such as that by Lehtonen et al. (2015) and Qian and Jin (2021), has demonstrated a strong correlation between phylogenetic metrics calculated from species-level and genus-level resolved phylogenetic trees. This methodology aligns with approaches used in other studies examining phylogenetic endemism for other plant groups (e.g., Qian et al., 2024a, b).

2.3. Endemism metrics

In this study, the range size of a species was defined as either the total number of regions in which the species was distributed or the total area of the regions (see below for details). We quantified taxonomic endemism (TE), also referred to as weighted endemism (Linder, 2001), for a specific region by summing 1/range size for all species within that region (i.e., we calculated 1/range size for each species and then summed resulting values of all species within the region), following the methodology outlined by Linder (2001) and Sandel et al. (2020). Additionally, we assessed the phylogenetic endemism (PE) of a geographic area by calculating the total phylogenetic branch length covered by taxa within the area, adjusting each branch length by dividing it by the global range size of its descendant clade, as in Rosauer et al. (2009) and Sandel et al. (2020). To enable comparisons between TE and PE, we normalized these values using the formula (x_i – minimum)/(maximum – minimum), where i = 1, 2, 3, … 390; this transformation allows TE and PE to vary from 0 to 1 across the 390 regions analyzed in this study. Recognizing that endemism metrics based solely on the number of geographic units might underestimate endemism in smaller geographic areas and that endemism metrics based solely on the total area of geographic units might underestimate endemism in larger ones, we computed TE and PE using these two approaches, as done in prior studies (e.g., Sandel et al., 2020; Cai et al., 2023b; Qian et al., 2024a, b). We then averaged the values obtained from both methods to determine the TE and PE for each region, mitigating potential biases associated with either method (Qian et al., 2024a, 2024b). The calculations of TE and PE were performed using the PDcalc package (available at https://github.com/davidnipperess/PDcalc).

In this investigation, we employed the categorical analysis of neo- and paleo-endemism (CANAPE; Mishler et al., 2014) to evaluate relative phylogenetic endemism (RPE) in each of the 390 regions, thereby identifying global centers of neo- and paleo-endemism in liverworts. Regions exhibiting significantly high RPE were classified as centers of paleo-endemism, which typically harbor ancient taxa that were once more widespread but are now confined to specific areas. Conversely, regions with notably low RPE were identified as centers of neo-endemism, often containing taxa that have recently diverged and remain endemic due to limited dispersal or migration to other locations. Additionally, regions displaying non-significant RPE but significantly high PE derived from both original and comparative phylogenetic trees were designated as centers of mixed endemism, representing centers of both neo- and paleo-endemism. The CANAPE methodology, detailed in Mishler et al. (2014), has also been applied to discern global endemism centers in angiosperms and ferns (Qian et al., 2024a, 2024b), as well as seed plants (Cai et al., 2023). For our calculations of PE and RPE significance, we utilized the canaper package (Nitta et al., 2023), and used the 'curveball' null model (Strona et al., 2014) using 1000 randomizations.

2.4. Effects of environment on endemism

Prior research has demonstrated that factors such as the current climatic conditions, historical climate changes since the Last Glacial Maximum, and environmental heterogeneity within a region significantly influence endemism in vascular plants (Cai et al., 2023; Qian et al., 2024a, 2024b). Therefore, our study analyzed the correlation between endemism and these variables across the 390 regions. We used the following six variables to represent current climatic variables: minimum temperature of the coldest month (T_min), mean annual temperature (T_mean), temperature seasonality (T_seas), precipitation during the driest month (P_min), annual precipitation (P_mean), and precipitation seasonality (P_seas). We used the following two variables to represent historical climatic variables: differences in annual precipitation and mean annual temperature between the Last Glacial Maximum and the present as a precipitation anomaly (P_anom) and temperature anomaly (T_anom), respectively. We assessed environmental variability through the standard deviations of elevation (E_sd), mean annual temperature (T_sd), and annual precipitation (P_sd). We used climate data at the 30-arc-second resolution available at the CHELSA climate database (https://chelsa-climate.org/bioclim) (Karger et al., 2017). For each region, the mean value of all values of a given variable within the region was calculated.

2.5. Statistical analysis

Given the strong correlation between taxonomic endemism (TE) and phylogenetic endemism (PE) (Pearson's r = 0.99), which suggests that statistical outcomes for TE and PE would be virtually identical, and considering the greater informativeness of PE, we opted not to perform statistical analyses for TE. Instead, our analyses concentrated on PE and relative phylogenetic endemism (RPE). Initially, we regressed each endemism metric against the eleven environmental variables and evaluated the relationships between pairwise variables using the standardized regression coefficients (Legendre and Legendre, 2012). Subsequently, we implemented variation partitioning analyses (Legendre and Legendre, 2012) to decompose the explained variance in each endemism metric into distinct components, utilizing adjusted coefficients of determination from simultaneous autoregressive error models (Kissling and Carl, 2008). For example, when current climatic variables and historical climate change variables were considered, variation partitioning analyses assessed the relative contributions of current climatic variables and historical climate change variables to the endemism metrics. Each variation partitioning analysis separated the explained variance into three portions: variance uniquely attributable to current climatic variables, variance uniquely attributable to historical climate change variables, and variance jointly explained by both variable sets. Additional variation partitioning analyses were conducted to compare the effects of current climatic variables against environmental variability variables, and historical climate change variables against environmental variability variables. We used simultaneous autoregressive (SAR) models of the spatial error type (Kissling and Carl, 2008) to determine standardized regression coefficients of the explanatory variable(s) and coefficients of determination in the regressions. Comprehensive variation partitioning analyses incorporating all three groups of environmental variables were also performed to segregate the explained variance into independent and shared effects. The aforementioned analyses were conducted using the software packages Spatial Analysis in Macroecology (Rangel et al., 2010) and SYSTAT (Wilkinson et al., 1992).

Principal component analysis (PCA) and t-test (based on Bonferroni-adjusted probability) were conducted to determine if the significant centers of neo-endemism, paleo-endemism, and mixed endemism can be distinguished by the environmental variables used in this study. PCA and t-test were conducted using PC-ORD (McCune and Mefford, 1999) and SYSTAT (Wilkinson et al., 1992), respectively.

3. Results 3.1. Geographic patterns of endemism

Phylogenetic endemism (PE) exhibited a strong correlation with taxonomic endemism (TE), as indicated by a Pearson correlation coefficient of 0.99. Consequently, the geographic distributions of PE and TE were remarkably similar, as illustrated in Fig. 1a and b. Regions that displayed significantly high values of TE and PE were predominantly found in tropical Asia (extending from India to Malesia), Madagascar, eastern Australia, and the Andes, as shown in Fig. 1. Areas characterized by notably low TE and PE included temperate regions of Eurasia and North America, parts of Africa, and eastern South America, also depicted in Fig. 1.

Geographic patterns of RPE were generally similar to those of TE and PE in the Southern Hemisphere, but differed substantially from those of TE and PE in the Northern Hemisphere (compare Fig. 2a with Fig. 1a and b). Regions with significantly high RPE were located in southern South America and tropical Asia, whereas regions with significantly low RPE were primarily located in northern South America and southern Africa (Fig. 2b).

Our analysis utilizing the CANAPE approach pinpointed 38 regions as significant centers of liverwort phylogenetic endemism. These included 7 regions of paleo-endemism, 7 of neo-endemism, and 24 exhibiting characteristics of mixed endemism, as depicted in Fig. 2c. The centers of paleo-endemism were primarily located in southern South America, tropical Asia, and New Zealand. In contrast, the main centers of neo-endemism were situated in southern Africa. The majority of regions demonstrating a combination of both paleo- and neo-endemism were found extending from Mexico to the tropical Andes, from India to Malesia, and across much of Australia, as illustrated in Fig. 2c.

3.2. Relationships of phylogenetic endemism metrics with environmental variables

Of the 11 environmental variables, variability in annual precipitation (P_sd) was most strongly correlated with PE, and it was followed by annual precipitation (P_mean) and precipitation during the driest month (P_min) (Fig. 3a). In contrast, precipitation seasonality (P_seas) and precipitation anomaly (P_anom) were not significantly correlated with PE (Fig. 3a). For RPE, of the 11 environmental variables, temperature seasonality (T_seas) was the strongest correlate; it was followed by P_seas (Fig. S2). While PE was correlated positively with minimum temperature of the coldest month (T_min), mean annual temperature (T_mean), annual precipitation (P_mean), and precipitation during the driest month (P_min), and negatively with temperature seasonality (T_seas) and temperature anomaly (T_anom), the opposite patterns were observed in RPE (Fig. S2). While PE was significantly and positively associated with all the three variables measuring environmental heterogeneity, no significant associations between RPE and the three variables were found (Fig. S2).

Current climatic variables accounted for 17.7% of the variation in PE and 25.5% in RPE. Historical climatic change variables explained 2.2% and 7.3% of the variation in PE and RPE, respectively. Environmental variability variables explained 24.8% of the variation in PE and less than 0.1% in RPE. Current climatic variables and historical climate change variables together explained 18.1% of the variation in PE and 25.1% in RPE, with the majority of the explained variation attributed independently to current climatic variables, as shown in Fig. 4a and Fig. S3a. The combination of current climatic variables and environmental variability variables accounted for 26.7% and 25.9% of the variation in PE and RPE, respectively. In the case of PE, this variation was largely explained jointly by these variables, as illustrated in Fig. 4b, whereas for RPE, it was explained independently by current climatic variables, as depicted in Fig. S3b. Furthermore, historical climatic change variables and environmental variability variables together explained 24.5% of the variation in PE and 6.3% in RPE; in the case of PE, the majority of the explained variation was attributed independently to environmental variability variables as shown in Fig. 4c; in the case of RPE, the majority of the explained variation was attributed independently to historical climatic change variables, as shown in Fig. S3c.

The three categories of environmental variables together accounted for 32.1% of the variation in PE and 28.4% in RPE. Independently, current climatic variables, historical climatic variables, and environmental variability variables explained 5.0%, 0.4%, and 8.9% of the variation in PE, respectively, as depicted in Fig. 5. For RPE, these variables independently explained 18.7%, 0.3%, and 1.2% of the variation, respectively, as shown in Fig. S4. These groups of environmental variables jointly accounted for 2.8% of the variation in PE and 0.6% in RPE, as illustrated in Figs. 5 and S4.

3.3. Discrimination of endemism centers across current climate gradient

The first two principal components derived from the six current climatic variables in the 38 endemism centers explained 85.7% of the variation in the current climatic variables (59.3% by PC1, 26.4% by PC2). Climatic conditions of the three categories of endemism centers (paleo-endemism, neo-endemism, and mixed endemism) did not differ on the PC1 axis (t-test, P > 0.50 in all the three pairs of endemism center categories) but differed significantly from one another on the PC2 axis in all the three pairwise comparisons (P < 0.05) (Fig. 6). Mean annual temperature and precipitation seasonality had stronger associations with PC2 than any other four variables (Table S1), compared with other current climatic variables. Paleo-endemism centers tended to be located in areas with lower mean annual temperature and lower precipitation seasonality whereas neo-endemism centers tended to be located in areas with higher mean annual temperature and higher precipitation seasonality, with mixed endemism centers tending to be located between paleo- and neo-endemism centers on the PC2 axis (Fig. 6).

4. Discussion

The present study demonstrated that geographic patterns of phylogenetic endemism (PE) and taxonomic endemism (TE) in liverworts across the world are greatly consistent (Fig. 1). This finding is in agreement with findings reported in previous studies for global angiosperms and ferns (Qian et al., 2024a, 2024b). Strong correlations between PE and TE are found in other groups of organisms across the world (e.g., trees, Sandel et al., 2020; angiosperm genera, Qian et al., 2024b; ferns, Qian et al., 2024a). This study found that regions with high liverwort endemism are located in tropical latitudes and the Southern Hemisphere. This finding is also highly consistent with those for angiosperms and ferns (Qian et al., 2024a, 2024b), suggesting that the same set of main mechanisms and factors might have shaped the patterns of endemism in liverworts observed in this study and in other groups of plants observed in previous studies. Temperate latitudes in the Northern Hemisphere have a lower endemism than those in the Southern Hemisphere. This difference might have partly resulted from differences in land areas between the two hemispheres. Temperate latitudes in the Northern Hemisphere have a much larger land area than those in the Southern Hemisphere. Furthermore, land area in the Northern Hemisphere tends to increase with latitude whereas land area in the Southern Hemisphere tends to decrease with latitude. Thus, the setting of the landmasses in the Southern Hemisphere favors to produce smaller ranges (Orme et al., 2006), which would lead to higher endemism (Sandel et al., 2020), compared with that in the Northern Hemisphere.

Our study identified 38 centers of endemism for liverworts, nearly all of which are located in either tropical latitudes or temperate latitudes in the Southern Hemisphere. The fact that most of the mixed endemism centers (i.e., having both neo-endemics and paleo-endemics) reflect that both diversification and persistence have played a crucial role in creating endemism centers. The endemism centers identified for liverworts largely overlap with those for angiosperms (Qian et al., 2024b), and also overlap with those for ferns to some degree (Qian et al., 2024a). Unlike Asia and South America which each have all the three categories of endemism centers, Africa only has neo-endemism centers for liverworts. For those endemism centers identified for all three groups of plants (liverworts in this study, angiosperms and ferns in Qian et al., 2024a, b), there is substantial inconsistency among the three plant groups in terms of whether a center is identified as a paleo-, neo-, or mixed endemism center. For example, Venezuela is a neo-endemism center for liverworts, a paleo-endemism center for ferns (Qian et al., 2024a), and a mixed endemism center for angiosperms (Qian et al., 2024b). This inconsistency reflects, at least to some degree, different evolutionary histories in different plant groups.

Our research demonstrates a positive correlation between phylogenetic endemism and topographic heterogeneity, aligning with prior findings on other plant groups such as ferns, angiosperms, and seed plants (Qian et al., 2024a, 2024b; Cai et al., 2023). Topographic heterogeneity often introduces dispersal barriers, potentially enhancing speciation rates and fostering the emergence of small-ranged lineages, thereby increasing endemism and contributing to neo-endemism (Quintero and Jetz, 2018). Additionally, regions with varied topography provide opportunities for species to adapt to climate changes over short altitudinal ranges, potentially lowering extinction risks and aiding in the preservation of ancient lineages, thus supporting paleo-endemism (Jump et al., 2009). Our findings suggest that areas with diverse topography tend to harbor both neo- and paleo-endemics, corroborating the observation that most identified centers of phylogenetic endemism in this study are mixed endemism centers located in topographically diverse areas.

Our study revealed that areas with more humid and warmer climates have higher phylogenetic endemism, but precipitation tended to have a greater effect on phylogenetic endemism than temperature. For example, among the six current climatic variables examined in this study, annual precipitation (P_mean) is most strongly correlated with phylogenetic endemism, and it is followed by precipitation during the driest month (P_min) (Fig. 3). This finding for phylogenetic endemism in liverworts is consistent with that observed by Qian et al. (2024a) for phylogenetic endemism in ferns, in which P_mean is a stronger correlate of phylogenetic endemism than mean annual temperature (T_mean) and minimum temperature of the coldest month (T_min), but is inconsistent with that in angiosperms (Qian et al., 2024b), for which P_mean and P_min have weaker effects, compared with T_mean and T_min. However, our finding is consistent with Qian and Chen (2016b) who found that liverwort species richness is less strongly correlated with mean annual temperature than with annual precipitation in China. Stronger effect of precipitation on species richness and endemism in liverworts may be related to the poikilohydric nature of liverworts, whose growth strongly depends on water availability (León-Vargas et al., 2006; Proctor et al., 2007).

Quaternary temperature change is negatively associated with phylogenetic endemism in liverworts, although such an association does not exist for Quaternary precipitation change. These results were also observed for phylogenetic endemism in angiosperms and ferns worldwide (Qian et al., 2024a, 2024b). Great temperature anomalies during Quaternary glacial cycles are likely linked to increased species extinctions (Sandel et al., 2011), which may result in reduced neo-endemism. Conversely, regions with lower temperature anomalies might have experienced greater climatic stability during glacial cycles, facilitating the evolution of lineages with narrow physiological tolerances and decreasing the extinction risk for small-ranged lineages. Consequently, climatically stable areas could serve as refuges for endemic taxa amid climate fluctuations (Dynesius and Jansson, 2000; Jansson, 2003), thereby supporting the survival of ancient lineages (paleo-endemics).

Current climatic variables, historical climate change variables, and environmental variability together explained about one-third of the variation in phylogenetic endemism in liverworts. This amount of explained variation is much lower than that reported in previous studies on phylogenetic endemism in angiosperms and ferns (Qian et al., 2024a, 2024b). The current climatic variables alone accounted for about 60% of the variation in angiosperm phylogenetic endemism (Qian et al., 2024b), and 55% of the variation in fern phylogenetic endemism (Qian et al., 2024a). Our study found that environmental variability within region explained independently more variation in phylogenetic endemism than independently either current climatic variables or historical climatic variables. However, about the half of the explained variation in liverwort phylogenetic endemism was jointly explained by current climatic variables and environmental variability; this suggests that these two sets of variables act jointly, to a marked degree, in determining patterns of liverwort endemism. The finding of our study that historical climate change variables explained less variation than current climatic variables for liverworts is consistent with the finding of Qian et al. (2024b) for angiosperms and the finding of Qian et al. (2024a) for ferns. Future studies should investigate why current climatic variables have such a lower power in explaining geographic variation in phylogenetic endemism in liverworts, compared with other groups of plants, and why historical climate change variables explain much less variation in phylogenetic endemism than current climatic variables in all the three groups of plants (i.e., liverworts, angiosperms, ferns).

Acknowledgements

We thank anonymous reviewers for their helpful comments.

Data accessibility statement

Liverwort species distribution data that were used in the present study are available at the Figshare Repository at https://doi.org/10.6084/m9.figshare.22587199. Climate data that were used in the present study are available at the CHELSA (https://chelsa-climate.org/bioclim).

CRediT authorship contribution statement

Hong Qian: Writing – review & editing, Writing – original draft, Investigation, Formal analysis, Data curation, Conceptualization. Shenhua Qian: Writing – review & editing, Formal analysis.

Declaration of competing interest

The authors do not have competing interest to declare.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2024.08.004.