1. Introduction

Mountain regions include only about 25% of the world's land area (excluding Antarctica) but host about 87% of terrestrial global biodiversity (Rahbek et al., 2019). Many mountain regions are home to many small-ranged species and centers of biodiversity hotspots (Myers et al., 2000). Mountain regions located in tropical and subtropical latitudes have high species richness partly because they are topographically heterogeneous and environmentally diverse. Because many biodiversity hotspots are located in mountain regions (Myers et al., 2000), understanding the underlying mechanisms driving species assembly along elevational gradients is of great interest in ecology and biogeography.

Climate is a major determinant of species diversity (Ricklefs, 2004). Some major climatic variables (e.g., temperature) that determine species distributions change with increasing elevation in the same way as they change with increasing latitude. However, for a given range of climatic variation that is found across a latitudinal gradient of over 1000 km, it can be found across an elevational gradient of only 1 km in a mountain region. Thus, climate diversity is more tightly packed across elevational gradients than across latitudinal gradients. For a given climatic gradient, because geographic distance is much shorter along an elevational gradient than along a latitudinal gradient, dispersal limitation plays a lesser role in driving species distributions along elevational gradients. Thus, species distributions are expected to better reflect their climatic niches along an elevational gradient than along a latitudinal gradient (Qian et al., 2023).

The Himalaya is one of the world's hotspots of biodiversity (Myers et al., 2000). The elevational gradient of the central Himalaya in Nepal is one of the longest elevational gradients in the world. It spans over 8700 m and includes a series of life zones from tropical forests at low elevations to nival areas at high elevations. Thus, the elevational gradient in Nepal is one of the best elevational gradients for addressing climate-based questions of biodiversity.

Species richness of the elevational gradient in Nepal has been well studied for angiosperms, pteridophytes, and bryophytes (e.g., Bhattarai et al., 2004; Grau et al., 2007; Qian et al., 2022; Rana et al., 2019). Because species composition of a biological community is determined by the interplay of ecological and evolutionary processes (Ricklefs, 1987), investigating spatial patterns and climatic correlates of phylogenetic structure (e.g., phylogenetic diversity, phylogenetic dispersion or relatedness) of biological assemblages over an elevational gradient is important to better understanding the assembly of species along the elevational gradient. Previous studies have investigated phylogenetic structure of angiosperms and ferns across elevational gradients in the Himalaya (e.g., Qian et al., 2019, 2023; Rana et al., 2019; Li et al., 2022). However, no studies have investigated spatial patterns and climatic correlates of phylogenetic structure in mosses along any elevational gradient in the Himalaya. Because mosses are important constituents of vegetation in the Himalaya, knowledge of geographic patterns and ecological causes of phylogenetic structure of mosses along elevational gradients in the Himalaya is critical to understanding the assembly of plant diversity in general, and moss diversity in particular, in the Himalaya. Here, we fill this critical knowledge gap.

In this study, we investigate the relationships between phylogenetic structure and climatic variables in moss assemblages along the elevational gradient of the central Himalaya located in Nepal. Specifically, we investigate the relationships of metrics of phylogenetic structure with elevation and climatic variables representing mean temperature and precipitation conditions, climate extremes, and climate seasonality. We assess the relative importance of temperature-related versus precipitation-related variables and the relative importance of climate extreme versus climate seasonality variables on the phylogenetic structure of moss assemblages along the elevational gradient.

2. Materials and methods 2.1. Study area

The area of this study is the Himalayan range that is located within the boundary of Nepal (80°04′–88°12′ E and 26°22′–30°27′ N). It covers 147, 516 km² and an elevational gradient ranging from 60 to 8848 m above sea level. The elevation gradient of Nepal was divided into vertical belts (N = 89), each of which had 100 m elevation, with the lowest belt located between 0 and 100 m above sea level.

2.2. Species distributional data

Distributional data of moss species in Nepal were compiled by Grau et al. (2007) based on the data published in Kattel and Adhikari (1992). We used the package U.Taxonstand (Zhang and Qian, 2023) to standardize botanical nomenclature of mosses based on the database Bryophyte Nomenclator (Brinda and Atwood, 2023). We combined infraspecific taxa with their respective species. As a result, 406 moss species in 196 genera and 63 families were included in this study. Each species was assigned to all vertical belts which were entirely or partially located within the elevation range of the species, as in previous studies (e.g., Patterson et al., 1996; Grau et al., 2007; Qian et al., 2019, 2022). Mosses in Nepal were distributed in elevational belts ranging from 100 m to 5100 m.

2.3. Phylogenetic tree

We used the phylogeny of mosses reported by Sanbonmatsu and Spalink (2022) as a backbone to generate a phylogenetic tree for the moss species in the present study. We used the package U.Taxonstand (Zhang and Qian, 2023) to standardize botanical nomenclature of mosses in the phylogeny based on Bryophyte Nomenclator (Brinda and Atwood, 2023). Duplicate species were removed. All the 63 families of mosses in our study were present in the phylogeny backbone. Of the 196 genera in this study, 190 (97%) were present in the phylogeny. We used the functions build.nodes.1 and Scenario 3 in the package U.PhyloMaker (Jin and Qian, 2023; also see Jin and Qian, 2022) to graft the remaining six genera onto the phylogeny based on most recent evidence about their phylogenetic position (Table S1). Of the 406 species present in this study, 292 (72%) were present in the phylogeny backbone. For those species of our study that were absent from the phylogeny backbone, we used the functions build.nodes.1 and Scenario 3 in the package U.PhyloMaker (Jin and Qian, 2023) to add them to the phylogeny. With the functions build.nodes.1 and Scenario 3 (Jin and Qian, 2023), U.PhyloMaker placed the nodes evenly between dated nodes and terminals within the family or genus (i.e., the node of the branches of genus and species combined within a family was set at 2/3 of the length of the branch of the family, and the node of the branches of species within a genus was set at 1/2 of the length of the branch of the genus). This approach, or a similar one, has been commonly used to generate phylogenies for plant species (e.g., Sandel et al., 2020; Cai et al., 2023). Previous studies (e.g., Lehtonen et al., 2015; Qian and Jin, 2021) showed that phylogenetic structure metrics, including those used in this study (see below for details), estimated on a tree resolved at the genus level are highly similar to those estimated on a tree resolved at the species level.

2.4. Metrics of phylogenetic diversity and dispersion

Faith's (1992) phylogenetic diversity (PD) and Webb et al.'s (2008) mean pairwise distance (MPD) and mean nearest taxon distance (MNTD) are commonly used metrics of phylogenetic diversity, and their standardized effect sizes (PD_ses, MPD_ses and MNTD_ses, respectively) are commonly used metrics quantifying phylogenetic dispersion or relatedness of species assemblages (Cadotte and Davies, 2016; Qian et al., 2021). Higher phylogenetic dispersion reflects lower phylogenetic relatedness or clustering. We used these metrics to quantify phylogenetic diversity and dispersion in the present study. PD represents the sum of the phylogenetic branch lengths linking the species in an assemblage (Faith, 1992). MPD quantifies the mean phylogenetic distance between all pairwise species in an assemblage, and MNTD quantifies the mean distance between each species within an assemblage and its closest relative (Webb et al., 2008). PD is positively correlated with species richness whereas MPD and MNTD are mathematically independent of species richness. PD_ses, MPD_ses and MNTD_ses were calculated using the following formula: X_ses = (X_obs − mean(X_null))/sd(X_null), where X_obs is the observed value of PD, MPD or MNTD, X_null is the average expected value for randomized assemblages, and sd(X_null) is the standard deviation of expected values among randomized assemblages. MPD_ses quantifies the more basal structure of the phylogenetic tree, while MNTD_ses and PD_ses quantify the more terminal structure of the phylogenetic tree (Mazel et al., 2016). Thus, PD_ses, MPD_ses and MNTD_ses measure the phylogenetic dispersion of assemblages at different evolutionary depths. For each of these three metrics of standardized effect sizes, a positive value reflects relative phylogenetic overdispersion of species (i.e., relatively low phylogenetic relatedness) while a negative value reflects relative phylogenetic clustering of species (i.e., relatively high phylogenetic relatedness). We used the package PhyloMeasures (Tsirogiannis and Sandel, 2016) to calculate the phylogenetic metrics.

2.5. Climate data

The following six climatic variables have been considered to influence phylogenetic diversity and dispersion of plants (Weigelt et al; Lu et al., 2018; Qian et al., 2023): mean annual temperature, minimum temperature of the coldest month, temperature seasonality, annual precipitation, precipitation during the driest month, and precipitation seasonality. We related phylogenetic metrics to these climatic variables. We obtained data of these climatic variables from the CHELSA climate database at the 30-arc-second resolution (https://chelsa-climate.org/bioclim) (Karger et al., 2017). The aforementioned six climatic variables corresponded to bio1, bio6, bio4, bio12, bio14, and bio15, respectively, in the CHELSA database. We calculated the mean value of each of the climate variables for each 100-m elevational belt. We considered bio1, bio6, bio4 a set of temperature-related variables, and bio12, bio14, and bio15 a set of precipitation-related variables; we considered bio6 and bio14 a set of climate extreme variables, and bio4 and bio15 a set of climate seasonality variables.

2.6. Data analysis

We excluded elevational belts with fewer than five moss species, as in previous studies (e.g., Rana et al., 2019; Qian et al., 2023), because species-poor assemblages may have extreme values for some metrics of phylogenetic structure (Fritz and Rahbek, 2012) and may thus produce results that are unreliable (Kamilar and Guidi, 2010). The elevational belts included in the final analyses of this study and species richness in each of the elevational belts were shown in Fig. 2. We graphically explored the relationship of each of the six phylogenetic metrics with elevation. Considering that area per elevational belt varied among elevational belts, which would affect both species richness due to the species−area relationship (i.e., the number of species increase with area; Rosenzweig, 1995) and Faith's phylogenetic diversity due to the strong association between species richness and phylogenetic diversity, we followed previous studies (e.g., Qian et al., 2022) to calculate area-corrected species richness and phylogenetic diversity, by dividing the number of species or phylogenetic diversity in an elevation belt by the log₁₀-transformed area of the belt.

Considering that the three metrics of standardized effect size were richness-standardized, our climate-based analyses focused on these metrics. Furthermore, because PD_ses was strongly correlated with MNTD_ses (Pearson's r = 0.92 in our data), suggesting that their relationships with climatic variables would be very similar, and because both MNTD and MPD calculate averages, accordingly, our climate-based analyses further focused on MPD_ses and MNTD_ses.

We conducted two sets of variation partitioning analysis (Legendre and Legendre, 2012) for each phylogenetic metric: one set determined whether temperature-related variables (bio1, bio6, bio4) or precipitation-related variables (bio12, bio14, bio15) have a stronger effect on the phylogenetic metric, the other set determined whether climate extreme variables (bio6, bio14) or climate seasonality variables (bio4, bio15) have a stronger influence on the phylogenetic metric. We used the ordinary least squares regression to conduct variation partitioning analyses. To account for non-linear relationships, we used a third-order polynomial model in each regression (i.e., linear, quadratic and cubic terms of each climatic variable under consideration were included in each model). We used SYSTAT v.7 (Wilkinson et al., 1992) for correlation and regression analyses.

3. Results

PD peaked at the mid elevation (about 2500 m) of the moss elevational gradient of Nepal (Fig. 1, Fig. 2), a pattern similar to species richness of the elevational gradient (Fig. 1). This pattern held when variation in area among elevational belts was accounted for (Fig. 1). In contrast, MNTD tended to be lowest near the mid elevation (Fig. 2b). MPD generally increased with elevation (Fig. 2c).

Both PD_ses and MNTD_ses decreased from the lowest elevation to about 1500 m, and then increased with increasing elevation (Fig. 2d and e). Thus, the elevational gradient of PD_ses and MNTD_ses can be divided into two segments. MPD_ses tended to show a triphasic (zig-zag) pattern comprising three segments along the elevational gradient: it increased with elevation from the lowest elevational belt up to ca. 2000 m, decreased with elevation up to ca. 3600 m, and then increased with elevation again (Fig. 2f).

The amount of variation in each phylogenetic metric that was explained by each climatic variable varied greatly among phylogenetic metrics and among climatic variables (Table 1). When phylogenetic diversity metrics were considered, temperature seasonality (bio4) was the strongest effect on both PD and MNTD whereas mean annual temperature (bio1) was the strongest effect on MPD (Table 1). Each of the six climatic variables examined in this study explained, on average, 67%, 59% and 75% of the variation in PD, MNTD and MPD, respectively. When the standardized effect sizes of these three phylogenetic diversity metrics were considered, temperature seasonality also was the strongest effect on both PD_ses and MNTD_ses, and mean annual temperature also was the strongest effect on MPD_ses (Table 1). Each of the six climatic variables examined in this study explained, on average, 68%, 37% and 67% of the variation in PD_ses, MNTD_ses and MPD_ses, respectively (Table 1).

When MNTD_ses and MPD_ses were focused and the relationship between either phylogenetic metric and each of the six climatic variables was assessed for each of their elevational segments noted above, temperature seasonality was the strongest effect on MNTD_ses in both elevational segments, and mean annual temperature was the strongest effect on MPD_ses in all the three elevational segments (Table 2). Each of the six climatic variables explained, on average, 84% and 74% of the variation in MNTD_ses in the low and high elevational segments, respectively, and 56%, 81% and 75% of the variation in MPD_ses in the low, middle and high elevational segments, respectively (Table 2).

When the whole elevational gradient was considered, MNTD_ses and MPD_ses both tended to increase with increasing elevation, and the increasing rate for MNTD_ses was greater than that for MPD_ses (regression coefficient β = 0.496 and 0.316, respectively). For the two elevational segments of MNTD_ses, the slope of decreasing MNTD_ses with increasing elevation for the low segment was about twice as steep as the slope of increasing MNTD_ses with increasing elevation for the high segment (β = −2.005 and 1.100, respectively; Fig. 3). For the three elevational segments of MPD_ses, the slope of increasing MPD_ses with increasing elevation for the low segment was shallower than the slope of decreasing MPD_ses with increasing elevation for the middle segment (β = 0.436 and −1.208, respectively); the slope of increasing MPD_ses with increasing elevation for the high segment was not significant (P = 0.325; Fig. 3).

When MNTD_ses or MPD_ses was related to the three temperature-related variables and the three precipitation-related variables simultaneously, these six climatic variables explained 87.7% and 94.4% of the variation in MNTD_ses and MPD_ses, respectively; the variation independently explained by the temperature-related variables was greater than that independently explained by precipitation-related variables by a factor of over 9 in both MNTD_ses and MPD_ses (Fig. 4a). When MNTD_ses or MPD_ses was related to the two climate extreme variables and the two climate seasonality variables simultaneously, the four variables explained 85.4% and 51.4% of the variation in MNTD_ses and MPD_ses, respectively; the variation independently explained by the climate seasonality variables was greater than that independently explained by the climate extreme variables in both MNTD_ses and MPD_ses (Fig. 4b).

4. Discussion

We explored patterns in three metrics of phylogenetic diversity (PD, MPD, MNTD) and their standardized effect sizes (PD_ses, MPD_ses, MNTD_ses, respectively) in moss assemblages across the elevational gradient of Nepal, and found that geographic patterns and climatic correlates of different phylogenetic metrics differ substantially in some cases (Fig. 2; Table 1). Among the three metrics of phylogenetic diversity, PD is hump-shaped, MNTD is U-shaped, and MPD tend to be monotonic (Fig. 2). The elevational patterns of PD and MNTD for mosses are similar to those of their counterparts in ferns across the elevational gradient of Nepal (compare Fig. 2a and b of this study with figs. S1a and g of Qian et al., 2023). The hump-shaped pattern of PD is likely driven by the hump-shaped pattern of species richness, because PD is strongly and positively correlated with species richness (r = 0.97). Because MNTD represents the average of the tip branch lengths of the closest taxa within an assemblage, the combination of the higher species richness and lower MNTD at middle elevations, compared with low and high elevations, indicates that divergence times of closely related species are younger in the middle of the elevational gradient than those at the lower and higher elevations for both mosses and ferns in Nepal. However, the elevational patterns of MPD are opposite between mosses and ferns over the elevational gradient: with increasing elevation, MPD tended to increase for mosses (Fig. 2c, this study) but tended to decrease for ferns (fig. S1d of Qian et al., 2023). Because MPD emphasizes deep phylogenetic relationships, the opposite elevational patterns between mosses and ferns suggest that deep evolutionary histories of these two groups of plants have differently influenced phylogenetic structure of their respective current species assemblages. The increase of MPD towards high elevations for mosses and low species richness at high elevations suggest that at high elevations there exist only a few species from very scattered lineages across the phylogeny. As a result, assemblages at high elevations are built up of distantly related species whereas in assemblages at low elevations, each clade contains several species, which reduces MPD (Qian et al., 2023).

When phylogenetic metrics standardized for species richness (i.e., PD_ses, MPD_ses, MNTD_ses) were considered, their geographic patterns are, to a large degree, similar. For example, all the three metrics tended to increase with increasing elevation across the elevational gradient. In particular, the pattern of PD_ses is similar to that of MNTD_ses (Fig. 2d–f). These led us to focus on MPD_ses and MNTD_ses in more detailed climate-based analyses. Our subsequent discussion focuses on these two phylogenetic metrics.

MPD_ses and MNTD_ses are larger at low elevations than at high elevations (Fig. 3). Furthermore, MPD_ses and MNTD_ses decrease with mean annual temperature and annual precipitation (Table 1). These patterns for mosses are, in general, contrary to those observed for angiosperms and ferns across the elevational gradient of Nepal (Qian et al., 2019, 2023). These patterns are also inconsistent with the tropical niche conservatism hypothesis, which was developed to explain the variation of phylogenetic diversity and dispersion in species of the lineages originating and diversifying in tropical climates. According to the tropical niche conservatism hypothesis, due to phylogenetic niche conservatism (i.e., the tendency of species to retain ancestral ecological distributions; Wiens and Donoghue, 2004), phylogenetic dispersion is lower, and thus phylogenetic relatedness is higher, for species assemblages in areas with colder and drier climates (Wiens and Donoghue, 2004). This hypothesis has been supported by many studies on angiosperms and polypod ferns (Lu et al., 2018; Qian et al., 2021), which originated and diversified in tropical climates (Takhtajan, 1969; Qian et al., 2021). The inconsistency of the patterns of mosses across the temperature and precipitation gradient of Nepal with the tropical niche conservatism hypothesis might suggest that mosses originated in extra-tropical climates. Indeed, previous studies have suggested that mosses originated in temperate climates (Shaw et al., 2005; Wu et al., 2022). Nevertheless, the patterns of MPD_ses and MNTD_ses oberserved in mosses across the elevational gradient of Nepal are consistent with the 'out of tropical lowlands' hypothesis (Qian and Ricklefs, 2016), which was proposed to explain patterns of phylogenetic structure in angiosperms along elevational gradients in tropical regions.

When more closely looking at the patterns of MPD_ses and MNTD_ses in mosses across the elevational gradient of Nepal, the elevational gradient of MNTD_ses can be roughly divided into two segments forming a "V" shape (i.e., first decreasing and then increasing with elevation; Fig. 3b) while that of MPD_ses can be roughly divided into three segments forming a zig-zag shape (first increasing, then decreasing, and finally increasing again with elevation; Fig. 3d). The V-shaped pattern of MNTD_ses may be due to that divergence times of closely related species are younger at mid elevations than those at the lower and higher elevations, as noted above. Interestingly, the general trend of MPD_ses in mosses is increasing with elevation whereas the general trend of MPD_ses in angiosperms is decreasing with elevation (Qian et al., 2019; note that NRI, which is equivalent to −1 × MPD_ses, was used in their study), both groups of plants have shown a zig-zag pattern in MPD_ses. Furthermore, a similar zig-zag pattern was found in MNTD_ses of ferns across the elevational gradient of Nepal (Qian et al., 2023). More interestingly, such a zig-zag pattern has been also observed in angiosperm assemblages along an elevational gradient in the eastern Himalaya (Li et al., 2022). Qian et al. (2019) attribute the causes of the zig-zag pattern of phylogenetic relatedness in angiosperms to the interplay of geophysical and eco-evolutionary processes. Specifically, they hypothesized that the patterns observed at low, middle, and high elevations resulted from plate tectonics, niche conservatism, and trait convergence, respectively. It is not certain whether the mechanisms proposed by Qian et al. (2019) to explain the zig-zag pattern of MPD_ses in angiosperm are applicable to explain the zig-zag pattern of MPD_ses observed in this study for mosses.

Our study found that temperature-related variables explained more variation in MPD_ses and MNTD_ses in mosses than did precipitation-related variables. Our finding is contrary to the finding of Qian and Chen (2016) that species richness of mosses is more strongly correlated with annual precipitation than with mean annual temperature in China. However, our finding for the relative importance of temperature- and precipitation-related variables on MPD_ses in mosses is consistent with the finding of Qian et al. (2023) for the relative importance of temperature- and precipitation-related variables on MPD_ses in ferns in Nepal. For MNTD_ses in ferns in Nepal, Qian et al. (2023) found that precipitation-related variables explained more variation than did temperature-related variables, which is inconsistent with our finding for MNTD_ses in mosses.

Our study found that climate seasonality variables explained more variation in MPD_ses and MNTD_ses in mosses than did climate extreme variables. This finding for MPD_ses in mosses is consistent with the finding of Qian et al. (2023) for MPD_ses in ferns in Nepal; however, Qian et al. (2023) found that these two groups of climatic variables explained about the same amount of the variation in MNTD_ses for ferns in Nepal. Qian and Chen (2016) found that in China, temperature seasonality explained more variation of species richness in mosses than did minimum temperature of the coldest month, which is consistent with our finding that climate seasonality variables have a stronger effect on MPD_ses and MNTD_ses in mosses than do climate extreme variables.

In conclusion, our study reveals that the composition of moss assemblages along the elevational gradient in Nepal shows strong signatures of evolutionary processes that cannot be linked to patterns of species richness alone. For a given climatic variable, different metrics of phylogenetic structure show different spatial and climatic patterns. However, all metrics measuring phylogenetic dispersion show that phylogenetic dispersion of species in moss assemblages tend to increase with increasing elevation and decreasing temperature. These patterns are not only contrary to those for angiosperms and ferns across the elevational gradient of Nepal but also contrary to the tropical niche conservatism hypothesis. However, the finding of our study that the pattern of MPD_ses is in a zig-zag pattern is consistent to those of Qian et al. (2019, 2023) for angiosperms and ferns across the elevational gradient of Nepal. Our study shows that temperature-related variables and climate seasonality variables are more important drivers of phylogenetic dispersion in mosses in Nepal, compared with precipitation-related variables and climate extreme variables, respectively.

Acknowledgements

We thank two anonymous reviewers for their helpful comments.

Data accessibility

Moss distribution data used in this study have been published and cited in this article. Climate data are available at the CHELSA climate database (https://chelsa-climate.org/bioclim).

CRediT authorship contribution statement

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

Declaration of competing interest

The authors have no 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.07.005.