Stochastic process drives the dissimilarity in biodiversity patterns between Pinus kwangtungensis coniferous forest and evergreen deciduous broad-leaved mixed forest in karst area

Study area and plot establishment

This study was conducted in Maolan National Nature Reserve, Yunnan-Guizhou Plateau, Southwest China. The reserve has a typical subtropical climate, with an annual temperature of 15.1 °C and annual rainfall of 1,374 mm. The highest temperature occurs in July (mean temperature of 15.1 °C), while the lowest is in January (5.2 °C). More than 80% of precipitation concentrates from May to October. Unlike typical subtropical evergreen broad-leaved forests, the karst area forms a different vegetation type with mixed evergreen and deciduous broad-leaved species due to its more heterogeneous habitats and poorer soil conditions. However, the evergreen coniferous forest featured by P. kwangtungensis usually distributes on the upper slopes due to the relatively less soil.

However, due to the greater than 50 m distance between each plot, which aims to avoid the confusing impacts of spatial autocorrelation, and P. kwangtungensis is normally distributed in habitats with higher altitudes or deeper slopes, we only established nine FDPs for both CF and EDBM.

To explore the difference in species diversity pattern between CF and EDBM and find out the driving factors of the difference, we received the field survey approval from the Maolan National Nature Reserve Administration of Guizhou Province and a series of FDPs were established in each vegetation type from September to November 2021. Since the least distance between any two sampling plots should be more than 50 m to avoid the impacts from spatial autocorrelation and ecotone, we established nine 20 m × 20 m FDPs in centralized distribution area of CF. In addition, to ensure the uniformity of repetition, we also established nine 20 m × 20 m FDPs in EDBM. Thus, a total of 18 20 m × 20 m FDPs were established in the study area (Table 1). In each FDP, the information of locations, elevation and soil conditions was recorded. All individuals in the plots were mapped, tagged, and identified to species, Furthermore, the diameter at breast height was measured to calculate species importance value (IV) (Eq. (1)), and determine vegetation community types. (1)${\displaystyle \mathrm{IV}=\left(R_f+R_d+R_{do}\right)/3}$ where R_f is the relative frequency, R_d is the relative density, R_do is the relative significance.

To explore how the abiotic factors influence the species composition, we collected the top layer of soil (0–20 cm) without the litter at five points in each FDP, and the soil samples were mixed to determine the soil physicochemical properties. A total of five indices, including soil pH value, soil water content, soil total nitrogen content, soil total phosphorus content, and soil exchangeable calcium content were measured. The soil pH was obtained by measuring soil suspension through a pH meter (PH500T, INESA, Shanghai, China). The soil total nitrogen content and soil exchangeable calcium content were determined through Elemental Analyzers (UNICUBE trace, Elementar, Langenselbold, Germany). The soil total phosphorus content was determined spectrophotometrically at 700 nm by a continuous flow automated analyzer (AA3, Bran+Luebbe, Norderstedt, Germany).

Differences in biodiversity patterns between CF and EDBM

Using the Margalef index, Shannon-Wiener index (Eq. (2)), Pielou’s evenness index (Eq. (3)), and rarefied richness (sampled 10 individuals), a one-way analysis of variance (ANOVA) was performed to test the differences in α-diversity between CF and EDBM. Shannon-Wiener index (H) can indicate the diversity pattern of a community due to its sensitivity to variations in species richness and evenness of species abundance. Pielou’s evenness index (E) is usually used to reflect the evenness of species abundance. Rarefied richness is used to characterize the species richness by controlling the confusing effects of sample size (Heck Jr, Van Belle & Simberloff, 1975). Additionally, the skewness of log-transformed species abundance depicts SAD symmetry and better quantifies the shape of SAD (Magurran, 2005). Usually, negative skewness refers to a higher proportion of rare species, while positive skewness reflects strong species dominance compared with a lognormal SAD (Fig. 1).

Metric multidimensional scaling (MDS) analysis was conducted on obtained SADs to examine the differences in species composition between communities (Norden et al., 2009). The similarity in species composition between types of FDPs was evaluated through the Chao-Jaccard abundance-based estimator (Chao et al., 2005). This approach can find a stable solution using several random starts, and standardizes the scaling in the result by a principal components rotation (Oksanen et al., 2013). This analysis was performed by using the “metaMDS” of the package “vegan” (Oksanen et al., 2013) in R3.4.2 (R Development Core Team, 2022). (2)${\displaystyle H=-\sum _{i=1}^S{P}_i\ln P_i}$

where s is the total number of species in the plot; P_i is the relative abundance of species i. (3)${\displaystyle E=H/H_{\max }}$

where H is the actual Shannon-Wiener index while H_max is the maximum of it, calculated by ln s.

Examination of the first-order assembly process between CF and EDBM

Following the PER-SIMPER framework proposed by Gibert & Escarguel (2019), this study examines the first-order assembly process influencing the dissimilarity in SADs between CF and EDBM. Based on Clarke (1993) SIMPER analysis, Gibert & Escarguel (2019) developed the PER-SIMPER analysis, which introduces a random permutation-based step. SIMPER analysis is a distance-based procedure that computes the relative contribution of each taxon to the OAD among groups of taxonomic assemblages. The Bray–Curtis coefficient for abundance was employed in SIMPER analysis to quantify the contribution of each taxon to the dissimilarity among groups of taxonomic assemblages (Clarke, 1993; Legendre & Legendre, 1998). Utilizing the empirical SIMPER pattern of a taxon locality occurrence data set, the PER-SIMPER analyzes 1,000 independently randomized occurrence data tables by fixing column or row to form the null SIMPER profile distribution (95% confidence intervals) of dispersal- and niche-assembly process. In addition, both sample richness and taxon uniquity remained unchanged under a maximally constrained null model. A third null SIMPER profile distribution was generated, corresponding to the null hypothesis that both dispersal- and niche-assembly processes dominate the taxonomic assembly. The null distribution closest to (if not including) the empirical SIMPER profile indicates the primary driver contributing to the observed taxonomic occurrence structuring among groups (details about the null permutation can be found in the paper of Gibert & Escarguel (2019)). For each permutation model, the empirical SIMPER profile is compared with the null SIMPER profiles by computing the logarithm of the sum of squared deviations between two profiles (E is as follows): ${\displaystyle E={\mathrm{Log}}_{10}\left(\sum _{i=1}^{i=p}{\left({{\overline{\gamma }}_i}_{\left(\mathrm{null}\right)}-{{\overline{\gamma }}_i}_{\left(\mathrm{obs}\right)}\right)}^2\right)}$

where ${\overline{\gamma }}_i=\frac{{\delta }_i}{\overline{\delta }}$ is the contribution of the ratio of δ to the mean of δ. Lower E indicates a closer similarity between the two compared SIMPER profiles. Thus, a PER-SIMPER analysis was conducted on SADs between CF and EDBM, using the “PER-SIMPER” function in R3.4.2 software (R Development Core Team, 2022).

Determination of abiotic factors on α-diversity

Redundancy analysis (RDA) is a multivariate statistical approach based on principal component analysis that is used to examine the correlation between explanatory variables (i.e., environmental factors) and response variables (i.e., species composition). Abiotic factors could play a crucial various role in shaping community structure (Xu et al., 2018). Therefore, RDA was used to examine the influence of soil physicochemical properties on dissimilarities in species composition between CF and EDBM in this study. This analysis was performed with the “vegan” package in R3.4.2 (R Development Core Team, 2022).

Difference in α-diversity between CF and EDBM

CF had a total of 46 species and was dominated by P. kwangtungensis, with an IV of 0.43. EDBM had a total of 128 species and was dominated by Boniodendron minus, with an IV of 0.06 (File S2). By comparing species diversity between CF and EDBM, significant differences were found in the Shannon-Wiener index, Margalef richness, Pielou evenness, rarefied richness, and skewness (Fig. 2), indicating significant differences in both components (species richness and evenness) of α-diversity between CF and EDBM. Even after controlling the impacts of abundance, the species richness of EDBM was still significantly higher than that of CF (Fig. 2D). In addition, the EDBM showed higher species dominance compared to CF (Fig. 2E). MDS result also showed that the species composition of CF was significantly different from that of EDBM (Fig. 2F).

First-order assembly processes driving dissimilarity in species composition between CF and EDBM

Results of PER-SIMPER analysis showed that the lowest E-values for dispersal- and niche-assembly profile, indicating that both niche-assembly and dispersal-assembly processes drove the dissimilarity in species composition between CF and EDBM (Fig. 3). Additionally, compared with niche-assembly process, the dispersal-assembly process had a lower E-values, indicating a more contribution of dispersal-assembly process to the dissimilarity in species composition.

Abiotic factors influencing the α-diversity

Results of RDA showed that the soil physicochemical properties determined the difference in α-diversity between CF and EDBM (Fig. 4). Additionally, the soil exchangeable calcium content, soil total phosphorus content, soil water content, and soil pH value had influence on the difference in α-diversity. Specifically, soil exchangeable calcium content and soil total phosphorus content showed a positive relationship with Margalef index, and a negative relationship with the other α-diversity indices. Soil total nitrogen had a relatively minor impact on the community composition.

Difference in α-diversity between CF and regional climax community

By comparing the CF with a typical regional climax community (EDBM), it was found that both the species richness and evenness of the CF were significantly lower (Fig. 2). The average individual density of the CF (0.297 individuals/m²) was also lower than that of the EDBM (0.377 individuals/m²). According to the species–energy hypothesis (Akatov et al., 2023), more individuals would increase species accumulation within the same regional scope, which might be one possible reason for the higher species richness of EDBM compared to CF (Chu et al., 2019; McGill, 2011). However, the rarefied richness of the EDBM was notably higher than that of the CF, indicating that the species richness of the CF remained lower even after controlling the cumulative effects of individuals. The biotic and abiotic filtering frameworks may play an important role (Münkemüller et al., 2020). The dominance of the CF was markedly lower than that of the EDBM, showing a high dominance of the EDBM (Fig. 3). The higher abundance and larger size of individuals compressed the ecological niche space of other species, indicating a much higher competitive ability of P. kwangtungensis compared to others, which limited colonization and survival of other species and formed monodominant species pattern (Stanley Harpole & Tilman, 2006; Wang & Cui, 2023). Similar patterns have been found in many evergreen coniferous forests, such as naturally occurring Pinus roxburghii forests in tropical regions, as well as artificially planted forests such as Pinus massoniana, Cryptomeria japonica, and Cunninghamia lanceolata, where evergreen coniferous species dominate and exhibit low species richness patterns (Sloan, Zimmerman & Sabat, 2007; Wang & Cui, 2023; Zang et al., 2021b). Meanwhile, the lack of subdominant species in CF allowed the colonization of rare species, contributing more to the difference in species composition (Lamanna et al., 2017; Zhang et al., 2014). This phenomenon could also be a primary factor leading to the scarcity of species under such biotic conditions. Some researchers have pointed out that evergreen coniferous species commonly exhibit allelopathy, which could inhibit the establishment of other tree species, thereby reducing species richness (Ehlers, Charpentier & Grøndahl, 2013). Other studies also showed that poor light resources and soil conditions in evergreen coniferous forests could also be limiting factors for species diversity patterns. For example, low light resource quality in evergreen coniferous forests limits the establishment of shade-tolerant tree species (Loke & Chisholm, 2023; Wang & Cui, 2023). This study showed that the CF had almost no shade-tolerant pioneer tree species, which supported the view to some extent. In addition, the CF is shade-intolerant, similar to deciduous broad-leaved tree species in eco-strategy. Therefore, the presence of CF somewhat restricts the ecological niche space of deciduous broad-leaved tree species, reducing species richness (especially in broad-leaved tree species (Cui et al., 2012; Wang & Cui, 2023)). In contrast, the EDBM is mixed with both evergreen and deciduous broad-leaved tree species. More deciduous broad-leaved tree species can increase species richness and evenness, resulting in a higher level of species diversity. Such diversity patterns also emphasized the importance of considering reducing deciduous tree species when performing P. kwangtungensis conservation to reduce the severe competition for lights. High habitat heterogeneity can also change species diversity patterns (Bar-Massada, Kent & Carmel, 2014; Brown et al., 2013). The habitat heterogeneity hypothesis suggests that increasing heterogeneity leads to an increase in small habitat types that can support plant growth and survival. The differing quality of small habitats can support the survival of different life forms, thereby maintaining a higher level of diversity (Bar-Massada, Kent & Carmel, 2014; Franklin et al., 2013). This study found that although both the CF and the EDBM were formed on limestone (as the parent rock), the microhabitat of the CF was relatively homogeneous without the presence of gullies, exposed rocks, or surface patterns. The EDBM exhibited higher levels of horizontal habitat heterogeneity, which provided a non-biological environmental foundation for maintaining diversity patterns and further promoted greater species diversity (Bar-Massada, Kent & Carmel, 2014; Brown et al., 2013; Franklin et al., 2013; Wang & Cui, 2023).

First-order assembly process driving the dissimilarity in species composition

The results showed that both the CF and the EDBM were primarily driven by dispersal- and niche-assembly processes, with the stochastic process playing a slightly larger role in determining species composition differences (Fig. 3). This indicated that the community assembly of CF may be mainly controlled by stochastic processes relative to deterministic processes. The CF community exhibited monodominance with a sharp decline in the number of individuals of other species (Fig. 2E). Both results indicated that apart from P. kwangtungensis, the number of individuals of other tree species was extremely few, with a significant variation in presence between plots (Fig. 2F). This reflected the uncertainty of companion species in the CF, indicating strong stochasticity on the occurrence of other species (Matthews et al., 2019; McGill, 2011). From the perspective of species occurrence, neutral processes dominate the probability of occurrence of most species, especially the rare species. These results could also be driven by demographic, phylogenetic and biogeographic histories of individuals, populations and species. In contrast, the EDBM are mostly common species sharing the same regional species pool (Chen et al., 2019). Functional traits adapt to variations in abiotic factors, while convergence of a trait value suggests co-occurring species often appeared in similar abiotic conditions, leading to stronger habitat filtering (Pappas, Fatichi & Burlando, 2016; Xu et al., 2018). Therefore, changes in functional traits and abiotic factors may by one of the factors in community assembly process. Maolan National Nature Reserve is a typical karst area where strong habitat filtering effects and habitat heterogeneity restrict the dispersal and colonization of most species (Gu et al., 2019; Wang et al., 2023). When species disperse into this region, which species become companion species is largely randomized in CF. Thus, although the α-diversity of CF is totally different from that of the EDBM, no certain evidence was found to support the necessity of such patterns.

Abiotic factors driving the difference in α-diversity

The results of this study showed that the soil exchangeable calcium content and soil total phosphorus content were the main factors driving the difference in α-diversity between the CF and EDBM. Previous studies showed that the phosphorus content of soils was the limiting factor on species diversity in tropical or subtropical forests (Wang et al., 2023; Zang et al., 2021a). However, our results indicated that the soil total nitrogen content had no effect on the regional difference in species diversity among different vegetation types (Fig. 4). The varied speeds of biological cycles in different soil nutrients might be one possible reason (Alvarez-Clare, Mack & Brooks, 2013; Wang et al., 2023). The view of phosphorus limitation was mainly due to the relatively slow phosphorus cycle (Alvarez-Clare, Mack & Brooks, 2013). Previous studies showed that the soil available phosphorus came from rock decay, which was relatively slow and limited the speed of phosphorus cycling (Laliberté et al., 2015; Zotz & Asshoff, 2010). More phosphorus was fixed in the plant tissue with increasing biomass, decreasing the available phosphorus content (Alvarez-Clare, Mack & Brooks, 2013). Additional studies suggested that the ${\mathrm{H}}_2{\mathrm{PO}}_4^{2-}$ tended to form insoluble complexes with Al³⁺ or Fe³⁺ under acidic soils (Laliberté et al., 2015; Zang et al., 2021a). However, our previous studies on karst restoration showed no significant phosphorus limitation, as nitrogen cycling may be more limiting to nutrient content than phosphorus cycling (Laliberté et al., 2015; Turner, Brenes-Arguedas & Condit, 2018; Wang et al., 2023). Such nitrogen limitation also supported the idea that the EDBM still performed potential succession compared with subtropical evergreen broad-leaved forests. Thus, the soil nitrogen content might dominate the soil fertility, driving the species diversity patterns. In addition, higher calcium content has been one key component of soil physicochemical properties in karst areas. A growing body of studies showed that diversity had a higher dependence on the soil calcium content (Batalha et al., 2015; Guo et al., 2019; Guo et al., 2017; Guo et al., 2015), which is also clearly demonstrated in this study.