1. Introduction

Revealing the geographical distribution patterns and underlying mechanisms of biodiversity has long been a fundamental goal in ecological research (Rahbek et al., 2019; White et al., 2023). Elevational gradients on mountains have garnered considerable attention when studying biodiversity patterns because mountains offer high environmental heterogeneity over short distances (McCain, 2009; Peters et al., 2016; Rahbek et al., 2019; Wang et al., 2025). Karst mountains, which are distinct for combining the characteristics of arid zones and mountain ecosystems, exhibit significant environmental gradients, with a diverse set of microhabitats within a relatively small geographic area (Wang et al., 2019; Zhang et al., 2021; Zhao et al., 2024; Frei et al., 2025). However, the steep terrain of karst mountains poses a challenge for conducting systematic sampling across continuous elevational gradients (Bednar Jr, 2005; Wang et al., 2019; Zhang et al., 2021). This has limited our understanding of both biodiversity patterns and the mechanisms driving these patterns in karst mountain systems thus far.

Previous studies have shown that elevational biodiversity patterns in mountains are highly variable and are not universal (Rahbek et al., 2019). These differences arise from variations among ecosystem type, geographic region, and biological group (Peters et al., 2016; Ahmad et al., 2020; Ratier et al., 2023; Dai et al., 2024). Though common elevational patterns have been identified in certain mountain systems (Rahbek, 1995; Rahbek et al., 2019), species richness generally tends either to decline as elevation increases or to peak at mid-elevations (McCain, 2009; Ahmad et al., 2025). It remains uncertain whether these patterns apply to karst mountains, which are shaped by erosion and feature steeper, more fragmented terrain (Geekiyanage et al., 2019; Wang et al., 2019; Fu et al., 2023). Furthermore, most biodiversity studies in mountain ecosystems focus on a single diversity metric, such as species richness, and this may lead to incomplete or biased conclusions (Naeem et al., 2016; Soto-Navarro et al., 2021; Wang et al., 2022). To address this gap, it is essential to incorporate multidimensional diversity information (e.g., taxonomy, phylogeny, functional traits, and community structure) when aiming to comprehensively understand community assembly processes and their underlying mechanisms in shaping forest elevational diversity patterns within karst mountain systems (Naeem et al., 2016; Ahmad et al., 2025).

Plant community elevational diversity patterns are widely recognized as the result of multiple ecological processes related to community assembly, including climatic stress, nutrient availability, topography, and interspecific interactions (Franklin et al., 2013; Pérez-Toledo et al., 2022; Chen et al., 2024; Wan et al., 2024). For instance, resource availability influences species distribution by shaping plant life history strategies (e.g., encompassing phenology, seed dispersal, and individual growth) (Felton et al., 2021; Prochazka et al., 2024). Extreme low temperatures constrain the elevational range of certain species, while soil nitrogen and phosphorus levels restrict plant growth and development (Du et al., 2020; Körner, 2021). On the other hand, topographic barriers create biological refugia that mitigate stress from regional climate factors, strong winds, and other environmental pressures (Finocchiaro et al., 2023; Frei et al., 2023). Beyond these abiotic drivers, biotic interactions such as competition and mutualism among neighboring species can foster niche differentiation and influence community composition (Luo and Chen, 2011; Godoy et al., 2020; Blüthgen and Staab, 2024). The unique topography and habitat conditions of karst mountains may amplify these ecological processes, making them especially critical in shaping plant community elevational diversity patterns (Jucker et al., 2018; Zhang et al., 2021; Fu et al., 2023). Thus, a better understanding of the maintenance mechanisms underlying these elevational diversity patterns in karst mountains requires a comprehensive consideration of multiple ecological processes rather than focusing on a single perspective.

In this study, we integrated multidimensional diversity metrics, including taxonomy, phylogeny, functional traits, and community structure, and quantified ecological processes related to community assembly to explore the elevational patterns of diversity and the underlying mechanisms maintaining them in karst forest communities. This study aims to address the following questions: (1) What are the elevational patterns of multidimensional diversity in karst mountain forests? We hypothesized that the multidimensional diversity of karst montane forests may exhibit a hump-shaped pattern along the elevational gradient, peaking at mid-elevations. (2) Which ecological processes predominantly drive community assembly in karst mountain forests? We hypothesized that the assembly process may shift from environmental filtering to competitive exclusion with increasing elevation. (3) What is the relative importance of different ecological processes in shaping the elevational patterns of multidimensional diversity in karst mountain forests? We hypothesized that topographic filtering may be a key ecological process in shaping the elevational patterns of multidimensional diversity in karst mountain forests.

2. Materials and methods 2.1. Study site

This study was conducted at the National Nature Reserve of Jinfo Mountain (elevation ranging from 800 m to 2251 m), located at the transition zone between the Yunnan-Guizhou Plateau and the Sichuan Basin in southwest China (Fig. 1a). Jinfo Mountain is one of the biodiversity conservation priority areas of China (Ministry of Ecology and Environment of the People's Republic of China, 2015) and is also a key component of the Southern China Karst (UNESCO World Heritage, 2014). It is renowned for its unique topography and rich biodiversity.

The study area is characterized by a humid monsoon climate, with an annual average temperature of 8.3 ℃ (Fig. 1c) and annual precipitation of 1395.5 mm (Zhou et al., 2019). The area preserves relatively intact primary evergreen broad-leaved forests, thereby harboring many rare and endangered plants such as Cathaya argyrophylla, Ginkgo biloba, Davidia involucrata, Tetracentron sinense, and so on (Qian et al., 2016, 2018).

Jinfo Mountain exemplifies a typical Karst Table Mountain with its step-like terrain formations (Waltham, 2009). For this study, vegetation communities were classified along elevational gradients into three terraces (Ⅰ–Ⅲ). This classification was based on an integrative assessment of factors such as topographic features, anthropogenic influences, microclimatic conditions, and vegetation history (Table 1).

2.2. Field sampling

In 2021 and 2022, fourteen permanent plots (20 m × 20 m) were established on Jinfo mountain, 100 m apart and spanning 800–2100 m in elevation. Each plot was set up in a northerly direction (Fig. 1b). As for the selection of plots, we synthesized various variables such as community succession, anthropogenic disturbance, and natural disturbance to ensure the representativeness of the selected sample plots. Each plot was divided into four quadrats (10 m × 10 m), and species inventories were conducted in each sub quadrat in a clockwise direction (Fig. S1). The detailed method is as follows: the diameter at breast height (DBH) of all woody plants in each sub quadrat was measured using vernier calipers. All individuals with DBH ≥1 cm were tagged with an aluminum label displaying a unique number. The coordinates, DBH, height and crown size were recorded for each tagged individual within each quadrat. Each individual was also identified to species, and specimens were collected. The species Latin names were standardized according to The Flora of China (http://www.cn-flora.ac.cn). Overall, a total of 225 species representing 143 genera and 66 families were documented in our plots, including 3 bamboo species, and 8 vine species (Table S1 and Fig. S2).

2.3. Phylogeny tree construction

We used the R package "U.PhyloMaker" (Jin and Qian, 2023) to generate a species-level phylogenetic tree (Fig. S3) for the woody plants in this study. Specifically, we used the functions build.nodes.1 and Scenario 3 of U.PhyloMaker to generate the phylogenetic tree. For missing species in the U.PhyloMaker built-in mega phylogeny, we used the BLADJ method to randomly select a species from the identified genus to allow for tree construction (Webb et al., 2008; Zhang et al., 2021). Phylogenetic trees generated by U.PhyloMaker are very effective and have been widely used in studies on community phylogenetics (Zhang et al., 2021; Qian et al., 2023, 2024; Qian, 2025).

2.4. Collection of functional traits

Following the protocols of Zhou et al. (2022), we compiled the traits from the Flora of China (FOC: https://www.iplant.cn/frps) using plant Latin names. For species with missing traits in FOC, we conducted a literature search on databases of China National Knowledge Infrastructure and Web of Science using keywords "Latin names + trait", and extracted relevant traits from these literatures. Finally, a dataset of 12 functional traits was obtained that related to environmental adaptation, resource acquisition, and population regeneration. These included four categorical traits (lifeform, leaf habit, leaf texture, leaf margin) and eight continuous traits (leaf length, leaf width, petiole length, maximum tree height, fruit length, number of veins, flowering period, fruiting period) (Table S2).

To evaluate the trait-phylogeny relationship along with the changing elevation, we quantified the phylogenetic signal using Blomberg's K for continuous traits (Blomberg et al., 2003) and statistic D for binary traits (Fritz and Purvis, 2010). We utilized the phylosig functions in "phytools" package for Blomberg's K values, and phylo.d functions in "caper" package for Statistic D values. Finally, four continuous trait (maximum tree height, number of veins, petiole length, fruit length) and four binary traits (leaf habit, leaf texture, leaf margin) showed significant phylogenetic signals, thus indicating strong phylogenetic niche conservatism (Table S3).

2.5. Metrics of community diversity and structure

Four diversity metrics (taxonomic diversity, TD; phylogenetic diversity, PD; functional diversity, FD; structural diversity, SD), two phylogenetic structure metrics (the nearest taxon index, NTI; the net relatedness index, NRI), and two physical structure metrics (community maximum height, H_max; community density) were calculated based on abundance estimates of plant species across the elevation range. A summary of these indices, including their significance, formulae, software used, and reference(s), is given in Table 2.

Furthermore, to compute functional distances between species, a principal component analysis (PCA) on the standardized trait data was conducted (Villéger et al., 2008), and the first five PCA axes (67.39% of total inertia; Table S4) were used to define FD. We divided tree height into distinct classes using equal intervals and calculated the Shannon diversity of height based on these height classes (Cheng et al., 2024; Dong et al., 2024). The Shannon diversity of DBH (diameter at breast height) and canopy were calculated in the same way. Finally, the SD was estimated as the weighted average of these three metrics.

2.6. Environmental variables

Based on environmental monitoring and sampling from plots along the elevational gradients during 2021 and 2022, fifteen environmental variables were derived to quantify four ecological processes (topographic filtering, climatic stress, competition and soil nutrient limitation). These variables comprise six climate variables (annual mean air temperature, MAT; annual mean relative air humidity, RH; accumulated temperature above 5 ℃, AT; coldness index, CI; annual mean soil temperature, MST; annual mean soil moisture, MSM), five soil nutrient variables (soil pH; soil total nitrogen, TN; soil total phosphorus, TP; soil organic carbon, SOC; soil available phosphorus, AP), three topographic variables (elevation; aspect, slope), and one biological variable (stem basal area, BA). Detailed descriptions of each variable, including equipment for and method of data acquisition, sampling design, ecological significance, and relevant references are provided in Table S5.

We standardized all explanatory variables (SD = 1 and mean = 0) to compare the regression coefficients, and any variables with variance inflation factor (VIF) > 10 were excluded from the analysis to avoid multicollinearity (Table S6). Ultimately, we retained nine variables in our analysis: (Ⅰ) three micro-climate variables (AT, RH, and MSM); (Ⅱ) three soil variables (pH, AP, SOC); (Ⅲ) one biotic variable (BA); (Ⅳ) two topography variables (slope and aspect). Detailed information for all variables was provided in Table S7.

2.7. Statistical analyses

Elevational patterns of multidimensional diversity and community structure were assessed using first-, second-, and third-order polynomial regression. The polynomial regression with the lowest Akaike information criterion (AICc) value (Anderson, 2008), according to the "MuMIn" package (Bartoń, 2024), was considered the best-fitting function (Table S8). We also conducted a Principal Coordinates Analysis (PCoA) based on a Bray–Curtis dissimilarity matrix to compare community species composition across terraces and tested for significance using Permutational Multivariate Analysis of Variance (PERMANOVA; 999 permutations, α = 0.05) in R package "vegan". In addition, we checked for a correlation between each diversity and community structure parameter using the corr.test function in the R package "psych".

To disentangle the drivers of multidimensional diversity and community structure, we first conducted a simple ordinary least squares (OLS) regression using the "vegan" R package. This regression allowed us to examine potential associations between each predictor and diversity and community structure. Additionally, we also performed a multiple regression analysis to explore the multivariate explanatory power of predictors in shaping the elevational patterns of diversity and structure. In multiple regression models, subset procedures were performed using the dredge function in the "MuMIn" package. We used ΔAICc = 2 as the cut-off point and selected the model with the highest weight as the optimal model among the subset of models with ΔAICc < 2 (Table S9) (Chen et al., 2024). Furthermore, we used the "glmm.hp" package to calculate individual R² and the relative importance of each factor for the optimal model (Lai et al., 2022, 2023).

All statistical analyses were performed in R v.4.3.2 (R Core Team, 2023). Graphs were generated using the "ggplot2" package (Wickham, 2016).

3. Results 3.1. Elevational patterns of multidimensional diversity

TD, PD, FD, and SD all correlated significantly with elevation (Fig. 2). TD (_AdjustedR² = 0.59, P < 0.01), FD (_AdjustedR² = 0.60, P < 0.01), and SD (_AdjustedR² = 0.59, P < 0.01) followed similar left parabolic patterns along the elevational gradient (Fig. 2). Meanwhile, PD decreased linearly as elevation increased (_AdjustedR² = 0.47, P < 0.01; Fig. 2a–c, d).

Community diversity and species composition differed significantly across elevations (Figs. S4–S6). All multidimensional diversity was higher at Terrace Ⅰ (800–1100 m) than at Terrace Ⅲ (1800–2100 m) (Fig. S5). Furthermore, Terraces Ⅰ and Ⅱ were both dominated by evergreen trees, with some deciduous species present. Terrace Ⅲ was mainly composed of evergreen trees and shrubs (Table S1).

3.2. Elevational patterns of community structure

The relationship between NRI and elevation was not significant (_AdjustedR² = 0.095, P = 0.29, Fig. 3a), while NTI showed a significant left-skewed decline as elevation increased (_AdjustedR² = 0.58, P < 0.01, Fig. 3b). Both the NRI and NTI values were greater than 0 at elevations above 1900 m and less than zero at elevations above 1900 m.

Furthermore, H_max showed an S-shaped decline as elevation increased (_AdjustedR² = 0.72, P < 0.01), and the stem density showed an inverted S-shaped increase (_AdjustedR² = 0.87, P < 0.01; Fig. 3c and d). Community density and H_max did not differ significantly between Terrace Ⅰ and Terrace Ⅱ, though in Terrace Ⅲ, H_max was significantly lower and community density was significantly higher than in both other terraces.

3.3. Drivers of multidimensional diversity and community structure

The elevational patterns of multidimensional diversity and community structure were jointly shaped by multiple ecological processes with varying explanatory weights (Fig. 6). TD and FD shared similar drivers along the elevational gradients, including climatic stress, topographic filtering, interspecific competition, and nutrient limitation (Table 3 and Fig. 6). PD was mainly affected by climatic stress. SD, community density and community maximum height were mainly shaped by climate stress, followed closely by soil nutrient limitation (Fig. 6). NTI and NRI were mainly controlled by climatic stress, interspecific competition and topographic screening (Table 3 and Fig. 6). Furthermore, interspecific competition was also found to importantly influence community density (Fig. 5).

Significant positive correlations were found among multidimensional diversity (Fig. 4a and b). NTI was also significantly positively correlated with NRI and FD. In addition, community stem density showed a significantly negative correlation with structural diversity, maximum height of community, and NRI. The relationships of the remaining community metrics were not significant (Fig. 4a).

4. Discussion 4.1. The refugia effect of karst mountains on biodiversity

Our results indicate that the elevational patterns of multidimensional diversity are similar across metrics, with higher diversity observed at low-to-mid elevations (Fig. 2). These findings aligned well with those of previous studies in karst mountains (Rahbek, 2005; Zu et al., 2024), such as the studies at Yulong Snow Mountain (Luo et al., 2019), Daming Mountain (Wan et al., 2024), Fanjing Mountain (Zu et al., 2024) and so on. These patterns are generally explained by the mid-domain effect hypothesis (Moreno et al., 2008), the water-energy balance hypothesis (Huang et al., 2017) or the ambient energy hypothesis (Turner, 2004).

Interestingly, we found that the refugia effect, mediated by topography, may also explain these patterns, thereby providing a new perspective on forest elevational patterns in karst mountains (Selwood and Zimmer, 2020; Frei et al., 2023). The microhabitats of karst mountains may serve as "climate refugia" at a local scale (Selwood and Zimmer, 2020; Frei et al., 2023). These refugia can buffer against harsh environmental stresses (Frei et al., 2023, 2025), thereby providing stable niches for ancient relict plant species with narrow ecological amplitudes to persist, including Cathaya argyrophylla, Davidia involucrata, and Tetracentron sinense (Qian et al., 2016, 2018). In our study area, there were two geographically separated steep cliffs situated within mid-elevation zones (Fig. 1b). In fact, these cliffs acted as effective "topographic refugia" for native species, impeding the upward dispersal of nonnative species while also isolating the area from anthropogenic disturbances (e.g., understory grazing, bamboo shoot harvesting, tourist activities, and so on) (Tables S1 and S2). This dual effect ensures long-term habitat stability within medium elevations in karst mountains (Rull, 2009; Frei et al., 2023).

4.2. Mechanisms of community assembly along elevational gradients in karst mountains

Our results revealed that the phylogenetic structure of forest communities is clustered at low-to-mid elevations (Terraces Ⅰ and Ⅱ), whereas it is over-dispersed at high elevations (Terrace Ⅲ) (Figs. 3 and S6). These results contradict most previous studies where environmental filtering dominated community assembly at harsher, high elevations (Luo et al., 2019; Yao et al., 2019; Zhang et al., 2021). Some studies have also supported the competition exclusion hypothesis at high elevations, such as in the Rocky Mountains (Bryant et al., 2008), Baisha Mountain (Chun and Lee, 2018), and Lasha Mountain (Wan et al., 2024). Therefore, we propose that, despite the widespread influence of environmental filtering driven primarily by climatic stress, competitive exclusion may dominate community assembly at high elevations (Qian, 2018; Luo et al., 2019; Wan et al., 2024).

The underlying factors driving this mechanism could be nutrient limitation in karst mountains (Reich and Oleksyn, 2004; Sundqvist et al., 2013). Particularly in severely soil-eroded karst mountains, nutrient limitations may exacerbate competitive exclusion at high elevations (Fig. 6a), thus potentially resulting in an overdispersed phylogenetic structure (Dantas de Paula et al., 2021; Peng et al., 2025). Additionally, the clustered phylogenetic structure may be linked to anthropogenic disturbance (Luo et al., 2019; Ahmad et al., 2020). Because intensity understory grazing and tourist activities strengthen environmental filtering, likely leading to clustered phylogenetic structures (Zhang et al., 2021; Wan et al., 2024; Zu et al., 2024). On the other hand, topographic barriers may promote the colonization of closely-related species at low elevations, creating an "accumulation effect" that may result in clustered phylogenetic structures in karst mountains (Vasseur and Fox, 2007; Blonder et al., 2018).

Consequently, we inferred that community assembly would shift from environmental filtering to competitive exclusion as elevation increased. This is a process modulated by topographic filtering. Notably, community assembly processes are more complex than previously expected, particularly with regard to the coupling relationships among community ecological processes, and thus these require further investigation.

4.3. Relative importance of different ecological processes on elevational biodiversity in karst mountains

The elevational pattern of biodiversity arises from the combined effects of multiple ecological processes. Furthermore, these processes exhibit complex coupling relationships (Fig. 6). First, we found that climatic variables correlated strongly with both diversity and structural indices (Table 3 and Fig. 5), suggesting that climatic stress is a major driver of local biodiversity pattern in karst mountains (Zhang et al., 2021; Wang et al., 2022; Dai et al., 2024; Ahmad et al., 2025). This correlation can be attributed to climatic stress whereby high-heat-demand species typically persist at warmer, low elevations but are excluded from colder, high elevations (Peters et al., 2016; Gheyret et al., 2020; Carvalho-Rocha et al., 2023). Previous studies have shown that variations in microclimate could affect forest diversity along elevational gradients by altering local plant survival, growth, and mortality (Ding et al., 2019; De Pauw et al., 2021; Zhang et al., 2021). The highly heterogeneous terrain of karst mountains may allow for more diverse local microclimates (Geekiyanage et al., 2019; Finocchiaro et al., 2023). Thus, future researchers must place greater emphasis on the role of microclimatic variables in community assembly in karst mountains.

Secondly, we observed significant differences in all community indicators between Terrace Ⅲ and Terraces Ⅰ and Ⅱ (Figs. S5 and S6), along with a notable explanatory power of topographic variables for PD, TD, PD (Table 3 and Fig. 5). Because topographic changes can buffer vertical climatic variation and prompt the redistribution of soil nutrients, they can thereby also drive changes in community composition along elevational gradients (Hollunder et al., 2021; Wen et al., 2024). These findings support the hypothesis that long-term topographic filtering can promote differences in species composition between localities and convergence in functional traits within localities in karst mountains (Jucker et al., 2018; Hollunder et al., 2021). This further highlights the distinctive role of topography on shaping the elevational patterns of biodiversity in karst mountain ecosystems (Rull, 2009; Schmitt et al., 2021; Frei et al., 2023).

Third, we found that phosphorus limitation is one of the key nutrient limitations affecting community composition and plant growth in karst mountains (Table 3 and Fig. 5). This is consistent with results from other subtropical forests (Cui et al., 2022; Ma et al., 2024; Prochazka et al., 2024). Because soil available phosphorus can be directly absorbed and utilized by plants, it primarily determines community diversity and plant growth (Wassen et al., 2005; Cui et al., 2022). Furthermore, soil organic carbon also shares an equal role in community composition because it is closely linked to microbe-mediated biochemical processes, by which it indirectly regulates plant nutrient acquisition strategy (Felton et al., 2021; Prochazka et al., 2024; Song et al., 2025). Nevertheless, other available nutrients (such as available nitrogen, available carbon, and trace elements) also play a critical role in influencing community composition and plant growth (Harpole et al., 2011; Ding et al., 2019; Ma et al., 2024). Future studies should pay more attention to these available nutrients to elucidate community composition more comprehensively.

Fourth, we found that interspecific competition correlated with phylogenetic structure. At high elevations, the phylogenetic structure was over-dispersed (Figs. 3 and Fig. 5). This may reflect the competitive exclusion driven by resource limitation (Kunwar et al., 2021; He et al., 2023). Specifically, high-elevation nutrient constraints could intensify the competitive exclusion effect (Chen and Lewis, 2024; Gridzak et al., 2024; Schurman et al., 2024), which may result in increased phylogenetic distances between coexisting species and ultimately reduce species richness (Cadotte et al., 2013). Notably, the interspecific relationships beyond competition, such as mutualism and parasitism, play a crucial role in shaping biodiversity patterns along elevational gradients (Boucher et al., 1982; Hale et al., 2020).

5. Conclusions

This study demonstrates that simultaneously considering both multidimensional diversity and multiple ecological processes provides great value for comprehensively understanding species composition and distribution along elevational gradients in karst mountains. The present findings confirm that topography, a critical factor in shaping biodiversity elevational patterns in karst mountains, can regulate community assembly through complex couplings with other ecological processes. In addition, we discovered that community assembly may shift from environmental filtering to competitive exclusion, as elevation increased. Overall, these insights enhance our understanding of elevational biodiversity patterns and the underlying mechanisms driving their formation in karst mountains.

Acknowledgments

This research was supported by the National Natural Science Foundation of China [32071652] to Y.Y. and [32571936, 32301327] to L.H. We thank Wenqi Wang, Xue Ouyang, Liang Zeng, Zhuohao Cui, Jinxian Gu, Yin Zhu and Xun Hu for help with the field work. We thank the Mount Jinfo Scenic Area Management Committee for field survey permission and support. We thank the editors and reviewers for their comments and constructive suggestions on this work. This work is part of the BEST (Biodiversity along Elevational gradients: Shifts and Transitions) research network (https://best-mountains.org).

CRediT authorship contribution statement

Lihua Zhou: Writing – original draft & review & editing, Conceptualization, Investigation, Formal analysis, Visualization. Yuxiao Long: Investigation, Data curation, Review. Li Huang: Investigation, draft & review & editing. Siwei Hu: Investigation, Data curation. Min Luo, Jingwen Deng, Lisha Jing: Investigation, Data curation. Wenbo Mou: Data curation. Mingyue Pang: Conceptualization. Yongchuan Yang: Conceptualization, Supervision.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

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