Synchronous and Asynchronous Variation of Taxonomic and Phylogenetic Diversity During the Succession of Pinus kesiya var. langbianensis Forest in Yunnan, China

TL;DR

This study examines how species and evolutionary diversity change during the growth of a pine forest in Yunnan, China, showing they follow a cosine pattern with peak diversity in mid-succession.

Contribution

The study reveals asynchronous variation between taxonomic and phylogenetic diversity and highlights the importance of mid-successional stages for biodiversity conservation.

Findings

Taxonomic and phylogenetic diversity peak at mid- to late-successional stages, but phylogenetic diversity peaks earlier.

Community stability increases with stronger niche differentiation, while environmental filtering dominates mid-succession.

Mid-successional stages should be prioritized in forest management to maintain biodiversity.

Abstract

The changing patterns of taxonomic diversity (TD) and phylogenetic diversity (PD) during forest succession can provide a reference for optimizing forest ecosystem management. The widely distributed Pinus kesiya var. langbianensis forest (PKF) in subtropical Yunnan, China, has important ecological and economic values. However, little is known about species diversity patterns and driving factors during the pine forest succession. Adopting the “space‐for‐time‐substitution” (SFTS) approach and on the basis of community data from different successional stages, we investigated the dynamics of TD and PD across PKF succession by integrating environmental and spatial variables. The results show that both TD and PD follow a cosine pattern during succession, peaking at mid‐ to late‐successional stages, but with TD lagging behind PD. TD responds more rapidly to changes in dominant environmental factors than PD. Moreover, there is an asynchronous oscillation between taxonomic β diversity (TβD) and phylogenetic β diversity (PβD). TβD is consistently greater than PβD and increases monotonically throughout succession, whereas the PβD still fluctuates in a cosine pattern. TD and PD are equally important in maintaining community stability, and the community becomes increasingly homogeneous and stable. Notably, early and late successional stages are dominated by competitive exclusion, whereas environmental filtering prevailed at mid‐succession. The mean temperature of the driest quarter (bio9) plays an environmental filtering role in the community composition at the ecological scale, and the precipitation of the coldest quarter (bio19) shapes the phylogenetic structure by influencing the regional species pool at the evolutionary scale. Neutral and deterministic processes jointly govern β diversity, but niche differentiation has an increasing domination, which supports the “successional continuum hypothesis”. Spatial effects must be explicitly considered in SFTS‐based successional studies. The management should prioritize conserving mid‐successional stages (peak diversity) and balancing environmental heterogeneity with dispersal limitation. Ecological‐evolutionary assembly linkages should be considered in the pinewood sustainable utilization and management.

Based on community data and the space‐for‐time‐substitution (SFTS) approach, this study shows that both taxonomic and phylogenetic diversity follow a cosine pattern during succession, with PD peaking earlier than TD. Community stability increases as niche differentiation strengthens, with competitive exclusion dominating early and late stages, and environmental filtering prevailing mid‐succession. Mid‐successional stages are critical for biodiversity conservation and should be prioritized in forest management.

Figures

Structure and morphological characteristics of P. kesiya var. langbianensis forest (PKF) across the study area (Photographed by C. Y. Wang). (a) Community appearance; (b) trunk and understory vegetation; (c) crown structure; (d) branches and needles.

Research area and sampling plots of PKF.

α diversity, community stability index, and phylogenetic structure of PKF at different successional stages. AVD, average variation degree; NTI, nearest taxon index; PD, phylogenetic diversity.

Taxonomic β diversity (TβD) and phylogenetic β diversity (PβD) of PKF at different successional stages.

Correlation between diversity index and explanatory variables. Significance levels: p < 0.01, *p < 0.05.

CCA ordination results and permutation test of PKF.

Mantel test of taxonomic β diversity (TβD) and phylogenetic β diversity (PβD) at different successional stages.

NMDS ordination diagram and permutation test of taxonomic β diversity (TβD) and phylogenetic β diversity (PβD) in PKF's succession processes.

Variation partitioning of PKF community TβD and PβD on the basis of environmental factors and spatial factors. Env.—the partition in which environmental effects are explained separately, light blue; Spa.—the partition in which spatial effects are explained separately, light yellow; Interactive part refers to the partition in which environmental and spatial effects are explained together. PβD, phylogenetic β diversity; TβD, taxonomic β diversity.

Standardized α and β diversity in taxonomy and phylogeny of PKF during succession. PD, phylogenetic diversity; PβD, phylogenetic β diversity; TD, taxonomic diversity; TβD, taxonomic β diversity.

Response of taxonomic diversity (TD) and phylogenetic diversity (PD) to dominant environmental factors at different successional stages of PKF.

Relative explanatory ability of taxonomic diversity (TD) and phylogenetic diversity (PD) to community stability in PKF's succession.

Full text

Introduction

1

In the 18th century, European foresters first discussed the concept of forest succession (Spurr 1952). Succession has remained a central theme in ecology for more than a century (Walker and Del Moral 2003; Prach and Walker 2011; Letten et al. 2014; Meiners et al. 2015; Martínez‐Ramos et al. 2021; Poorter et al. 2024). Vegetation succession is a dynamic process of community assembly (Norden et al. 2009; Lebrija‐Trejos et al. 2010; Omidipour et al. 2021). As succession progresses, intra‐ or interspecific interactions intensify, the light‐use efficiency within the community increases, and resource allocation shifts. Consequently, community assembly patterns become increasingly heterogeneous among communities (Lohbeck et al. 2013, 2014; Purschke et al. 2013).

Two traditional theories have been proposed to explain community assembly during succession: the niche theory (Diamond 1975; Webb 2000; Webb et al. 2002) and the neutral theory (Hubbell 2011). The niche theory holds that community assembly is dominated by deterministic processes, such as environmental filtering or competitive exclusion (Chave 2004), whereas the neutral theory holds that community assembly is controlled by dispersal limitation or stochastic processes (Hubbell 2011). However, many researchers have increasingly recognized that these two processes coexist and complement each other in community assembly dynamics (Niu et al. 2009; Poorter et al. 2010), as their relative importance shifts across successional stages (Purschke et al. 2013; Lohbeck et al. 2014). There is still uncertainty in community assembly processes during succession. For example, the secondary and primary forests exhibit distinct assembly mechanisms at each stage (Bhaskar et al. 2014). The community assembly pattern of secondary forests is deterministic to a great extent, driven by the species' life history strategies and interspecific competition (Norden et al. 2012; Zhang et al. 2014; Lanta et al. 2023). In tropical forests, environmental filtering predominates (Whitfeld et al. 2012), whereas in subtropical forests, deterministic forces exert a weaker influence during succession (Mi et al. 2016).

Vegetation dynamics are accurately revealed by long‐term plot studies (Warming 1895; Lasky et al. 2014) and can also provide insights into continuous processes that cannot be captured by the space‐for‐time substitution (SFTS) approach (Smale et al. 2016; Fickert and Grüninger 2018). However, there are many challenges for conventional experiments and observations, requiring time frames well beyond instant results (Fukami and Wardle 2005; Johnson and Miyanishi 2008). The SFTS approach not only requires fewer human, material, and financial resources than the long‐term plot method but also produces reasonably reliable results (Li et al. 2016). It is appropriate for studying convergent secondary succession with low biodiversity, rapid species turnover, and infrequent, low‐intensity disturbances (Walker et al. 2010). As a result, the SFTS approach is widely used to study succession by many ecologists (Cowles 1899; Whittaker 1975; Pickett 1989; Blois et al. 2013; Song et al. 2019; Romanowski et al. 2021; Kharouba and Williams 2024; van der Sande et al. 2024). Then, will the two methods lead to significant differences in species diversity patterns in the process of succession? Subsequently, numerous studies have evaluated the reliability and applicability of the SFTS method by comparing it with long‐term monitoring. For example, Foster and Tilman (2000) demonstrated that SFTS‐derived data closely predicted temporal changes in species abundance, with discrepancies largely attributable to extreme drought events. Similarly, SFTS predictions of plant compositional changes in eastern North America were consistent with 72% of those derived from long‐term plots (Miao et al. 2018). These findings collectively highlight both the extensive use and practical validation of SFTS. More recently, the method has been widely applied to infer temporal dynamics in both taxonomic and phylogenetic diversity within communities (Li et al. 2015; Liu et al. 2016; Geng et al. 2022; Fu et al. 2023; Yuan et al. 2023).

Forest community structure and biodiversity exhibit uncertainty in succession (Howard and Lee 2003; Yuan et al. 2011; Marcilio‐Silva et al. 2016; Chai et al. 2023). Taxonomic diversity (TD) and phylogenetic diversity (PD) complement each other to provide an enhanced perspective of community assembly, which is driven by deterministic or stochastic processes in succession (Norden et al. 2009; Pausas and Verdu 2010; Letcher et al. 2012; Purschke et al. 2013; Jin et al. 2015; Jarzyna and Jetz 2016; Liu et al. 2016), as well as evolutionary or ecological variation (Webb et al. 2002; Cavender‐Bares et al. 2009; Norden et al. 2012). However, TD and PD show different patterns of change during succession (Purschke et al. 2013, 2017; Vellend et al. 2013; Chai et al. 2016; Yu et al. 2019), and have different effects on ecosystem functioning (Ouyang et al. 2016; Yuan et al. 2016). This suggests that they may capture different ecological processes. Meanwhile, phylogenetic community structure changes from clustering to overdispersion or randomization from early to late succession, or shows an opposite pattern (Letcher 2010; Norden et al. 2012; Stadler et al. 2017; Yu et al. 2019; Mastrogianni et al. 2024). TD exhibits uncertainty during succession, and four main patterns can be distinguished: TD can monotonically decrease or increase with ongoing succession, peak at middle successional stages, or may not show any temporal trend (Odum 1983; Howard and Lee 2003; Tilman et al. 2006; Baniya et al. 2009; Rolo et al. 2016). PD may increase significantly with the trajectory of community succession (Marcilio‐Silva et al. 2016; Muscarella et al. 2016; McKone et al. 2021; Karimi et al. 2022); decrease (Verdú et al. 2009; Li et al. 2015; Chai et al. 2016; Lazzaro et al. 2020; Kong et al. 2023); or show no significant trend over time (Bruelheide et al. 2011; Norden et al. 2012; Mi et al. 2016; Xu et al. 2020; Chai et al. 2023). Therefore, it is particularly noteworthy that studies on TD and PD patterns during succession have yielded a series of puzzling and contradictory conclusions. Such discrepancies may stem from differences in regional climate and biogeography (Chai et al. 2016), forest types (Bruelheide et al. 2011; Norden et al. 2012), types of succession (Muscarella et al. 2016), the completeness of coverage along the successional continuum (Karimi et al. 2022), and the scale effects associated with different research approaches (Miao et al. 2018). At present, there is still a lack of systematic studies comparing the patterns of TD and PD within the same successional framework. A study that jointly examines TD and PD across a full successional gradient while accounting for spatial effects in the SFTS framework is essential for reconciling these contradictory patterns and clarifying the mechanisms of community assembly.

The significance of spatial ecological structure has been widely acknowledged (Legendre 1993). Because ecological data are inherently collected across geographic space, spatial relationships among sampling sites often lead to spatial dependence in community patterns (Borcard et al. 2004). Although spatial eigenvector‐based approaches have been widely applied to detect spatial structures across ecological systems, including ferns in Amazonia (Borcard et al. 2004), benthic algae (Wu et al. 2021), and fish diversity in floodplain lakes (Virgilio et al. 2022), their role in succession studies remains underappreciated. Classical SFTS‐based succession studies (Parrish and Bazzaz 1982; Yuan et al. 2023) implicitly assume that spatial processes are negligible, leading to an underestimation of dispersal limitation and spatial structuring. Therefore, explicitly incorporating spatial factors is essential for accurately interpreting successional dynamics and avoiding misattribution of spatial patterns to temporal processes.

Forest succession is a fundamental process shaping terrestrial biodiversity worldwide (Crouzeilles et al. 2016). As one of the largest genera of conifers, pine (Pinus) forests are distributed across the entire Northern Hemisphere (Eckert and Hall 2006) and play a key role as pioneer communities in post‐disturbance landscapes (Sousa 1984). In particular, subtropical and tropical pine forests exhibit relatively rapid successional turnover, making them ideal model systems for investigating successional dynamics within observable time scales (Stephenson and van Mantgem 2005). Pinus kesiya var. langbianensis, Simao pine, is a geographical variant of P. kesiya Royle ex Gordon. Paleontological and molecular biological evidence support that this species may have differentiated southward from Pinus prekesiya with the continuous uplift of the Qinghai‐Tibet Plateau in the Neogene and the gradual strengthening of the South Asian monsoon climate (Xing et al. 2010). P. kesiya var. langbianensis is valuable for timber, pulp, and turpentine; more than 80% of turpentine production comes from it in Yunnan (Wu and Dang 1992). It is one of the typical pioneer tree species, spreading by winged seeds. Normally, P. kesiya var. langbianensis forest (PKF, Figure 1) will be replaced by monsoon evergreen broad‐leaved forests (Song et al. 2011). PKF is mainly distributed in the south subtropical and tropical lowland areas to the west of Ailao Mountain in Yunnan, China (Wu and Zhu 1987; Wen et al. 2010). As a typical secondary but long‐lasting forest, PKF stands are 1.48 × 10^6^ hm^2^, accounting for 6.5% of the forested area in Yunnan, and have a standing stock of about 1.50 × 10^9^ m^3^, giving it significant value for ecosystem service and carbon sequence (Wen et al. 2010).

Structure and morphological characteristics of P. kesiya var. langbianensis forest (PKF) across the study area (Photographed by C. Y. Wang). (a) Community appearance; (b) trunk and understory vegetation; (c) crown structure; (d) branches and needles.

Furthermore, PKF represents a transitional forest type, shifting from a relatively homogeneous coniferous community to a more heterogeneous broadleaved forest community. This transition reflects a fundamental ecological question concerning temporal changes in species composition, making PKF suitable for applying the SFTS approach to investigate long‐term successional processes (Tang et al. 2013; Ashton and Zhu 2020). Meanwhile, secondary pine forests such as PKF are increasingly becoming major components of global forest biomes and are facing challenges of adaptation and recovery under climate change (Zhou et al. 2024; Tudor et al. 2025). In this context, the PKF, which is distributed widely across the subtropical regions of southwestern China, provides a valuable model system for studying pine forest succession, with both regional specificity and global significance. Studies have been conducted on PKF's phylogeny (Zhao et al. 2018), aboveground biomass (Li et al. 2018; Fan et al. 2024), and carbon storage (Gao et al. 2021); however, few studies have explored the community assembly patterns throughout its successional process. To achieve a better understanding of the mechanisms underlying pine forest community assembly under climate change, we analyzed the changes in TD and PD during PKF succession on the basis of the SFTS approach. We aim to understand: (1) How similar are TD and PD at each successional phase, and how do potential driving factors shape them? (2) Do successional shifts in TD and PD show dissimilar community assembly processes? (3) What role does the spatial effect inherent in the SFTS approach play in explaining temporal variation? Although this study is geographically restricted to Yunnan, PKF shares similar ecological and evolutionary backgrounds with the subtropical pine forests widely distributed from Southeast Asia to South Asia. Therefore, the findings of this study not only enhance our understanding of the community assembly mechanisms of pine forests during succession, but also provide general insights into the successional dynamics of subtropical pine forests worldwide under the influence of climate change and human activities.

Materials and Methods

2

Study Area

2.1

This study was conducted in southwestern Yunnan Province, China (Figure 2), covering the geographic range from 21°05′–25°47′ N and 98°08′–102°25′ E. The region experiences a subtropical monsoon climate with distinct dry and wet seasons. The terrain slopes from north to south, resulting in notable topographic heterogeneity. Climatically, the area exhibits a large diurnal temperature range (mean annual temperature: 17.6°C) coupled with a narrow annual temperature variation. Annual precipitation averages 1270 mm, concentrated primarily during the wet season (Yang 1990). Precipitation exhibits pronounced seasonality, with approximately 85% of the annual amount occurring between May and October. The region has an accumulated annual temperature of 6400°C–7300°C and receives 2281–2453 h of sunshine per year. Soils are classified as acidic red soils, typical of subtropical monsoon climates.

Research area and sampling plots of PKF.

Sampling Design

2.2

From 2021 to 2023, we conducted a comprehensive survey of PKF in Yunnan using a stratified–spatially constrained sampling design on the basis of the typical plot method (Wang, Fang, et al. 2020). We first identified the continuous distribution area of PKF using Global Biodiversity Information Facility (GBIF, https://www.gbif.org/) and Chinese Virtual Herbarium (CVH, https://www.cvh.ac.cn/) databases, combined with forest resource inventory data and remote sensing imageries, and used the forest distribution as the primary stratum. On this basis, elevation, slope, aspect, and landform type were further considered to ensure that the plots covered the major environmental gradients. To minimize the influence of spatial autocorrelation, all sampling plots were separated by a minimum distance of at least 5 km. The size of each plot was defined as 20 m × 30 m. Each woody plant with a diameter at breast height (DBH) ≥ 2.5 cm and a height greater than 1.3 m was measured. The species name, individual count, base diameter, DBH, height, crown size, and growth status of the trees were recorded. Furthermore, we calculated tree density (ind./ha), tree cover (%), and stand height (m), which provide essential information for understanding community changes during forest succession (Table S1). As for shrub and herb species, five 5 m × 5 m and 1 m × 1 m subplots were established in the four corners and center of the main plot. Within these subplots, all plant species were identified, and for each species, base diameter, crown size, height, coverage, and individual number were measured. Additionally, any shrub or herb species present within the main plot but not encountered in the subplots were also recorded. A total of 144 plots were investigated in this study (Figure 2).

Classification of Successional Stages and Inference of Stand Age

2.3

Simao pine, as a typical pioneer tree species, often exhibits clustered and synchronous regeneration (Liu et al. 2024). Therefore, the age of the largest individual within a stand is generally considered a reliable proxy for the overall stand age (Alberdi et al. 2013). In this study, we employed the SFTS approach by selecting stands representing different successional stages to construct a successional sequence (Fu et al. 2023). The age range of each PKF stand was inferred on the basis of putative DBH–age relationships, supplemented by local forestry records, consulted with forestry or nature reserve managers, and interviewed with nearby villagers (Li 2010; Tang et al. 2020). Without doubt, this method involves a certain degree of uncertainty, mainly arising from microhabitat‐driven variation in individual growth and recall bias in interview‐based information. To assess and minimize these uncertainties, we implemented the following measures before conducting data analyses: (1) a plot was included only when the interview information was cross‐validated with the DBH structure of trees in the stand; and (2) considering the potential methodological errors, we set an age estimation uncertainty of ±3 years. Then five stage groups were classified according to the Age‐class and Age‐group Division Standards for Dominant Tree Species (LY/T 2908–2017, https://std.samr.gov.cn/), an official forestry industry standard promulgated by National Forestry and Grassland Administration of the People's Republic of China, that is, young forest (≤ 20 years old, N = 12), middle‐aged forest (21–30 years old, N = 27), near‐mature forest (31–40 years old, N = 37), mature forest (41–60 years old, N = 41), and overmature forest (≥ 61 years old, N = 27). N is the number of plots surveyed in each group. It is important to note that although the young forest stage was defined according to the standard criteria of LY/T 2908–2017 stands, no sapling‐dominated plots or early post‐establishment cohorts were encountered during the field survey. All individuals in this category had DBH values typically ≥ 10 cm, indicating that these stands had already passed the highly dynamic initial regeneration phase. Therefore, the earliest stage in this study represents an early developmental phase after the sapling stage rather than newly regenerated stands. The plot information is shown in Table S1. This stage classification ensures overall inference regarding community assembly patterns during succession. It has been widely demonstrated to be reliable and effective in similar studies of secondary forest succession (Abbas et al. 2019). Additionally, the pine forest plots of different age classes can be used as spatial substitutes for a “time series” under the premise that PKF is predominantly distributed across the study area and the climatic background remains consistent. The plots are applicable to quantify the law of niche distribution at different stages, which is suitable as a substitute index for the SFTS method in the processes of succession, thereby revealing the laws of community succession (Oliver and Larson 1996; Zhang et al. 2021; da Costa et al. 2025; Lanta et al. 2025).

Explanatory Variables

2.4

The bioclimate data were derived from the WorldClim database (http://www.worldclim.com/, accessed on 9 March 2024), including 19 bioclimatic factors. Topographic data (elevation, aspect, slope) and geographic coordinates were collected using handheld GPS receivers for all plots. Plot‐specific fire occurrences were recorded during field surveys and verified using local forest management archives. Monthly burned area data were sourced from the MODIS database (MCD64A1.061, https://developers.google.cn/earth‐engine/datasets/catalog/MODIS_061_MCD64A1, accessed on 9 March 2024). Burned area data were aggregated to annual averages using the GEE cloud computing platform (https://earthengine.google.com/, accessed on 10 March 2024) for subsequent analysis. During field surveys, we recorded plot‐level signatures of human disturbance. Annual Mean Human Activity Index data for 2000–2020 were retrieved from the Data Center of the Chinese Academy of Sciences (https://figshare.com/, accessed on 9 March 2024). Because the field‐based wildfire and human footprint data were collected as discrete ordinal variables, their applicability in quantitative analyses requiring continuous predictors was limited. To reduce the subjectivity and observer's bias inherent in discrete, plot‐level disturbance scoring and to retain more quantitative information, we additionally extracted continuous wildfire and human footprint layers from spatial databases and subsequently evaluated the correlation between the field‐measured disturbance levels and the corresponding database values. The results showed a significant and consistent correlation between the two datasets (Figure S1), indicating that the database variables reliably capture the disturbance gradients observed in the field. Therefore, we used the remote‐sensing‐based data of wildfire and human footprint. The spatial resolution of the above environmental data is 1 × 1 km. Using the GPS‐recorded latitude and longitude of each field plot, environmental data were extracted in ArcGIS (version 10.7). The explanatory variables were chosen as reliable proxies for key regional‐scale environmental gradients, also applicable at the plot scale. A critical feature of our sampling strategy is that the minimum distance between any two plots exceeds 5 km, ensuring that each plot is spatially independent from the 1 km grid environmental data. The plot network was effectively transformed into a discretized sampling of broad‐scale environmental gradients. Macroscale climatic variables represent long‐term environmental constraints that determine species' physiology, colonization, and successional trajectories, even at fine spatial scales. Thus, these variables do not represent microhabitat conditions; rather, they capture the distinct regional context of each plot, including climatic limitations, disturbance history, and the impacts of anthropogenic activities. These contextual factors filter the potential species pool and systematically influence the successional pathways observable at each plot. Such use of environmental variables to investigate regional climatic drivers of successional processes has been widely adopted (Janečková et al. 2024).

Because ecological and environmental data are derived from the Earth's surface, spatial autocorrelation among sampling points is often observed (Borcard et al. 2011). Spatial variables are represented by eigenvectors derived from the principal coordinates of the neighborhood matrix (PCNM) on the basis of community geographic coordinates (Borcard and Legendre 2002; Borcard et al. 2004). PCNM allows the spatial impact on local community development to be evaluated on multiple scales. In order to study the explanatory account of spatial effects, the spatial information data in our study were first transformed into the Cartesian coordinate system on the basis of ArcGIS software. We constructed a full model using redundancy analysis (RDA) and performed forward selection to obtain 16 principal component axes with significant positive eigenvalues for further variable screening (Table S2). A geographic distance matrix was also calculated using Euclidean transformation to serve as the explanatory variables of spatial characteristics. Spatial variables were generated using the “dbmem” function in the “adespatial” package (Guénard and Legendre 2022).

To avoid multicollinearity among explanatory variables, we first assessed pairwise relationships using Spearman's correlation (Figure S2). Subsequently, on the basis of canonical correspondence analysis (CCA) contribution rates and interpretability, we removed variables with low contributions and high intercorrelations. We applied a Hellinger transformation to the species‐abundance data. Then, we computed variance inflation factors (VIFs) and iteratively removed the variable with the highest VIF if it exceeded 5, repeating this process until all remaining VIFs were ≤ 5. The significance of both the overall model and individual environmental factors was assessed (Luo et al. 2023). Finally, 22 explanatory variables meeting the criteria (absolute Spearman's correlation coefficient < 0.8 and VIF ≤ 5) were retained (Bolded in Table S3).

Selection and Calculation of Diversity Metrics

2.5

To comprehensively characterize changes in community features and assembly mechanisms during PKF succession from multiple dimensions, we employed a multidimensional set of diversity metrics, including taxonomic, phylogenetic, and structural indices. These metrics capture different aspects of community structure and complement one another in revealing underlying patterns. The indices are introduced in detail in the following subsections.

Taxonomic Diversity Metrics

2.5.1

Species richness, Shannon–Wiener index, Simpson index, and Pielou evenness index were used to evaluate the species diversity of PKF at different successional stages. The above metrics describe the current species composition and abundance distribution of the community, reflecting species richness, patterns of individual abundance, and dominance structure. They are highly sensitive to species turnover and changes in abundance patterns during succession, and therefore provide complementary information on community structure that cannot be captured by phylogenetic metrics alone (Lanta et al. 2023). The calculation formulas for these diversity indices are as follows: (Simpson 1949; Spellerberg and Fedor 2003; Zhang et al. 2012):

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where S is the total number of species in the community; P _ i _ is the proportion of the total abundance of the ith species; N _ i _ is the sum of the number of a certain species i in the plots; N is the total number of individuals of this species in the community plots.

The Sørensen dissimilarity index was used to quantify taxonomic β diversity (TβD) of plant communities in the study area, because it is widely used in vegetation ecology and provides a robust measure of species turnover, particularly suitable for examining community replacement processes driven by succession (Baselga 2010; Duan et al. 2024). Because TβD does not incorporate phylogenetic relationships, it provides a critical reference against which phylogenetic β diversity (PβD) can be compared. This allows us to distinguish whether community differences arise primarily from species turnover itself or from phylogenetic divergence associated with successional differentiation. The equation is calculated as follows:

$$[eqn]$$

where A is the number of individuals of common species in community B and community C; B and C are the number of species individuals unique to community B and community C, respectively.

Community Stability Index

2.5.2

The average variation degree (AVD) was used to evaluate the community stability of PKF, and it is a stability metric that allows comparisons across different successional stages (Xun et al. 2021). It has been demonstrated in multiple studies to effectively capture fluctuations in community abundance, and it reflects community organization from both functional and structural perspectives (Xun et al. 2021). The AVD values were calculated by the following equation:

$$[eqn]$$

Here, xi denotes the rarefied abundance of the species in a single sample, x¯i is the mean value of xi, and δi signifies the standard deviation of rarefied abundances within each sample group. k denotes the number of samples within a sample group, and n represents the number of species in each sample group. A smaller AVD indicates higher community stability. In this study, we used 1 − AVD to represent community stability, such that higher values correspond to more stable communities.

Construction of Phylogenetic Tree

2.5.3

The “Plantlist” package was used to batch‐query the family and species information for all species recorded in the plots (Zhang et al. 2019). A phylogenetic tree was constructed for subsequent analysis using the “V.PhyloMaker” package (Jin and Qian 2019). Figure S3 shows the phylogenetic tree of the communities across the five successional stages.

Phylogenetic Diversity and Phylogenetic Structure

2.5.4

PD represents the sum of branch lengths in a phylogenetic tree connecting all species within a community, reflecting their shared evolutionary history. It serves as a fundamental metric for quantifying phylogenetic diversity (Faith 1992). The nearest taxon index (NTI) was used to infer community assembly processes (Webb et al. 2002). NTI quantifies phylogenetic community structure by measuring the clustering of closely related taxa at the tips of the phylogeny (Swenson et al. 2006). In this study, assuming phylogenetic conservatism of functional traits (Wiens et al. 2010; Li et al. 2022), we interpreted NTI values as follows: NTI > 0, phylogenetic clustering (closely related species co‐occur), suggesting environmental filtering dominance. NTI < 0, phylogenetic overdispersion (distantly related species coexist), suggesting competitive exclusion dominance. If NTI = 0, no significant deviation from random assembly, indicating stochastic processes (Webb et al. 2002). PD and NTI are primarily used to describe the evolutionary dimension of communities. The calculation was implemented in the “pd” function and the “ses.mntd” function of the “picante” package (Kembel et al. 2010). The calculation formula is as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where C is the minimal phylogenetic subtree connecting all species in the community, c represents a branch in C, and L _ c _ is the branch length of c. SES.PD calculates the standardized effect size of PD. PD_observed_ represents the observed PD in a community, PD_random_ is the mean PD from null communities, and sdPD_random_ is the standard deviation of PD in null communities. We utilized the “taxa.labels” null model and performed 999 iterations to satisfy statistical criteria (Jarzyna et al. 2021). MNTD_s_ represent the mean phylogenetic distance between each species and its closest relative in the observed community. When 999 null model simulations were used to generate random communities, MNTD_mds_ is the mean MNTD from null communities. SD (MNTD_mds_) stands for the standard deviation of MNTD in null communities (Webb et al. 2002).

The phylogenetic Sørensen index was used to quantify phylogenetic dissimilarity between communities (Stegen et al. 2012; Duan et al. 2024). It focuses on the sharing degree of phylogenetic branches among communities, representing PβD. The calculation formula is as follows:

$$[eqn]$$

where L _ A _ and L _ B _ are the total branch lengths of the two communities, and L _ C _ is the shared branch length.

Therefore, the combined use of taxonomic, phylogenetic, and community structural indices provides a multidimensional and complementary framework for understanding the community assembly mechanisms during PKF succession.

Data Analysis

2.6

The α and β diversity indices and AVD across five successional stages were compared using one‐way ANOVA, followed by Duncan's multiple‐range test to identify significant differences, with the significance level controlled at p = 0.05 for all multiple comparisons. To explore the relationship between α diversity and environmental factors, a correlation heatmap was generated on the basis of the Spearman correlation test using the “psych” package. We employed Spearman's correlation because of its robustness to non‐normal data distributions. Further, the detrended correspondence analysis (DCA) was used to judge the adoption mode of ordination analysis. When the gradient length of the first DCA axis exceeds 4.0, CCA is recommended; when it is below 3.0, RDA is more appropriate (Turktas et al. 2012). For intermediate gradient lengths (3.0–4.0), either CCA or RDA may be used (Bastow Wilson 2012). Finally, permutation tests (999 permutations) were performed on the first two ordination axes to assess their relationships with explanatory variables, on the basis of the results of DCA and CCA using the “decorana” and “anova.cca” functions in the “vegan” package. This data‐driven approach (DCA followed by CCA/RDA) ensures that the most appropriate unimodal or linear model is selected to reveal the relationships between community composition and environmental gradients.

To explore β diversity patterns during succession and their relationships with environmental factors, we first transformed β diversity indices into a distance matrix and constructed non‐metric multidimensional scaling (NMDS) ordination plots for both TβD and PβD. NMDS was selected for its ability to visualize complex community dissimilarities on the basis of any distance measure without assuming linear relationships. Statistical significance was assessed through 999 permutations of the Monte Carlo simulation. Furthermore, we used permutation tests to identify the key environmental factors that significantly influence β diversity, on the basis of the “envfit” function from the “vegan” package (Oksanen et al. 2019). Finally, Mantel tests were performed using the “mantel_test” function from the “linkET” package to explore the influence of environmental influences on β diversity at each successional stage (Tan et al. 2013). The Mantel test is appropriate here as it tests the correlation between two distance matrices, aligning with the β diversity data structure.

To investigate the influences of niche and neutral processes on β diversity patterns, we integrated bioclimatic factors, topographic factors, wildfire, and human footprint as environmental variables, and performed variation partitioning analysis (VPA) with spatial variables to quantify their relative contributions. VPA was specifically chosen to disentangle the pure and shared effects of environmental filtering (niche processes) and spatial dispersal (neutral processes) on community assembly. VPA was implemented using the “rdacca.hp” function from the “rdacca.hp” package (Lai et al. 2022). Additionally, we used linear regression to model how the explanatory power of environmental and spatial variables on β diversity changed throughout succession. Finally, standardized taxonomic and phylogenetic diversity indices were employed to examine variations in α and β diversity across successional stages, whereas polynomial regression was applied to evaluate relationships between diversity metrics and community stability. Polynomial regression was used to allow for the detection of potential non‐linear (e.g., hump‐shaped or U‐shaped) relationships between diversity and stability, which are theoretically expected during succession. In this study, box plots and polynomial fitting curves were generated using Origin 2024 software, and all other analyses were conducted in R version 4.4.1 (https://www.r‐project.org; R Core Team 2024).

Results

3

Changes in Community Structure and Taxonomic and Phylogenetic α Diversity of PKF During Succession

3.1

Community structure properties showed clear and significant differences among the five successional stages (Figure S4). Tree density exhibited a general decreasing trend from young to overmature forest (Figure S4a), with young forest having significantly higher tree densities compared with mature and overmature forest. Tree coverage also declined with succession (Figure S4b), with young and middle‐aged forest showing significantly higher canopy cover than later successional stages. In contrast, stand height increased markedly along the successional gradient (Figure S4c).

Taxonomic α diversity (TαD) of PKF significantly changed during succession. The species richness (SR) in near‐mature forest and mature forest was significantly higher than in young forest and middle‐aged forest (p < 0.05, Figure 3a), whereas the SR in overmature forest declined non‐significantly (p > 0.05, Figure 3a). The Shannon–Wiener index of young forest and middle‐aged forest was significantly lower than that of other stages (p < 0.05, Figure 3b). Throughout succession, the Simpson index decreased significantly, with young and middle‐aged forest showing higher values (p < 0.05, Figure 3c). The Pielou evenness index increased significantly and was lower in young and middle‐aged forest (p < 0.05, Figure 3d). Although community stability showed no significant differences among stages, it increased progressively during succession (Figure 3e,f).

α diversity, community stability index, and phylogenetic structure of PKF at different successional stages. AVD, average variation degree; NTI, nearest taxon index; PD, phylogenetic diversity.

Similar to SR patterns, both PD and SES.PD increased during the early succession and then declined. PD and SES.PD in the young forest and the middle‐aged forest were significantly lower than those in the mature forest (p < 0.05, Figure 3g,h). During succession, the NTI first increased and then decreased (Figure 3i). NTI values in young forest and middle‐aged forest were less than 0, in near‐mature forest were greater than 0, in mature forest were near 0, and in overmature forest were less than 0. This trajectory reflected a shift from phylogenetic overdispersion to clustering, then to randomness, and back to overdispersion. NTI dynamics demonstrated that different community assembly mechanisms dominated the local communities at different successional stages (p < 0.05, Figure 3i).

Taxonomic and Phylogenetic β Diversity of PKF During Succession

3.2

According to the Sørensen dissimilarity index of PKF at different successional stages (Figure 4a), TβD exhibited a monotonic increasing trend. In Figure 4b, PβD showed a fluctuating pattern, first increased and then decreased. PβD reached the maximum in the mature forest, which was significantly higher than at all other stages (p < 0.05). Moreover, PβD remained lower than TβD throughout all successional stages (Figure 4).

Taxonomic β diversity (TβD) and phylogenetic β diversity (PβD) of PKF at different successional stages.

Relationship Between α Diversity and Explanatory Variables

3.3

At different successional stages, the correlations between α diversity and explanatory variables were distinct (Figure 5). In the young forest, no α diversity or stability indices correlated significantly with environmental variables. In the middle‐aged forest, SR was negatively correlated with bio7 but positively correlated with bio14. The Shannon–Wiener index was positively correlated with bio14 and bio19. PD and SES.PD were positively correlated with slope. The stability index was positively correlated with bio3. In the near‐mature forest, PD and SES.PD were positively correlated with human footprint. NTI was positively correlated with altitude but negatively correlated with bio9. In the mature forest, SR and the Shannon–Wiener index were positively correlated with bio15, whereas the Simpson index was negatively correlated with bio15. PD and SES.PD were also positively correlated with bio15, and the stability index was positively correlated with bio3. In the overmature forest, SR was positively correlated with slope. The Shannon–Wiener index was positively correlated with human footprint, whereas the Simpson index was negatively correlated with both the human footprint and altitude, and the Pielou evenness index was negatively correlated with altitude. NTI was negatively correlated with bio9. The stability index was positively correlated with bio9, but negatively correlated with the human footprint. In general, the effects of environmental factors on forest diversity depended on the successional stage. Human impacts mainly occurred at later stages, and wildfire had no significant effect on TD and PD change during succession.

*Correlation between diversity index and explanatory variables. Significance levels: **p < 0.01, p < 0.05.

We first conducted DCA ordination and found that the first‐axis gradient length exceeded 4.0 (Table S4). Therefore, CCA ordination was used to explore the dominant drivers of local‐scale successional processes. The successional process exhibited a heterogeneous pattern that changed progressively under specific environmental drivers (Figure 6a). Permutation test results indicated that bio9 was the key environmental factor driving the successional dynamics of PKF (Figure 6b).

CCA ordination results and permutation test of PKF.

Relationship Between β Diversity and Explanatory Variables

3.4

On the basis of the Mantel test, we analyzed the dominant environmental factors influencing TβD and PβD at each successional stage (Figure 7; Table S5). The dominant environmental factors influencing β diversity varied across successional stages from young to overmature forest. In the young forest, β diversity was not significantly correlated with any environmental factor. In the middle‐aged forest, TβD was significantly correlated with bio7 and bio15, whereas PβD showed no significant correlations. In the near‐mature forest, TβD was significantly correlated with bio7, bio9, bio15, and wildfire. In the mature forest, TβD was significantly correlated with bio3 and bio15, whereas PβD was significantly correlated with bio9 and bio15. In the overmature forest, TβD was significantly correlated with bio7, bio9, and human footprint, and PβD was significantly correlated with bio9 and human footprint. The influence of environmental factors on TβD and PβD was variable among successional stages and shifted with forest maturity. At early stages, β diversity patterns were primarily shaped by stochastic processes, whereas at later stages, human disturbance became the dominant driver.

Mantel test of taxonomic β diversity (TβD) and phylogenetic β diversity (PβD) at different successional stages.

NMDS and permutation tests were used to analyze the dominant driving factors of β diversity during succession at the regional scale. For TβD, the different successional stages were divided into five distinct clusters and gradually shifted from the third quadrant to the first quadrant (Figure 8a). The permutation test results indicated that bio9 was the main environmental driver of TβD during succession (Figure 8b). As the succession progressed, the degree of phylogenetic differentiation among PKF slightly increased (Figure 8c), and bio19 emerged as the primary driver of PβD (Figure 8d). Other environmental factors showed no significant contribution to PβD.

NMDS ordination diagram and permutation test of taxonomic β diversity (TβD) and phylogenetic β diversity (PβD) in PKF's succession processes.

Effects of Spatial Factors and Environmental Factors on the Succession of PKF

3.5

VPA results showed that, except for the middle‐aged forest, the contribution of spatial factors to both TβD and PβD was greater than that of environmental factors at all other successional stages (Figure 9). Furthermore, by fitting the explanatory power of spatial factors and environmental factors at each successional stage, we found that spatial factors consistently explained more variation than environmental factors (Figure S5). However, the explanatory power of spatial factors gradually declined as succession progressed, whereas that of environmental factors increased significantly (Figure S5).

Variation partitioning of PKF community TβD and PβD on the basis of environmental factors and spatial factors. Env.—the partition in which environmental effects are explained separately, light blue; Spa.—the partition in which spatial effects are explained separately, light yellow; Interactive part refers to the partition in which environmental and spatial effects are explained together. PβD, phylogenetic β diversity; TβD, taxonomic β diversity.

Discussion

4

Changes of α and β Diversity in Succession Processes

4.1

The TD and PD of PKF both followed a cosine‐shaped pattern (Figures 3a,g and 10), a trajectory that transcends the traditional framework of linear or hump‐shaped curves and reveals diversity fluctuations driven by multiple stages and mechanisms, which can be implied as an extension of the classic “hump‐shaped diversity curve” (Geng et al. 2022). These results align with species diversity trends observed during the succession of abandoned lands in southern Illinois, USA (Bazzaz 1968; Kassen et al. 2000). At the early stage of succession, although niches were not saturated, the extremely high tree density likely intensified intraspecific competition for limited resources (Figure S4a), resulting in unstable community habitats, low niche differentiation, and less survival of pioneer species. This may explain the lower levels of both TD and PD (Chai et al. 2016; Fu et al. 2023). The priority effect led to a strong stochastic process (Figure 9; Viana et al. 2016), which explains the lack of significant relationships between species diversity and environmental factors (Figure 5). The succession progresses with decreasing tree density and increasing stand height (Figure S4a,c). The microclimatic and micro‐environmental conditions within communities gradually change, allowing more species to colonize. Intense niche differentiation promotes increases in both TD and PD, including the establishment of more shade‐tolerant species in the understory (Lebrija‐Trejos et al. 2016). In the later phases of succession, light becomes a more limiting factor. With the community structure shifting toward a taller and more open canopy (Figure S4b,c), light limitation is intensifying through deepening vertical shading, in consistent with poor understory light conditions during late succession (Baldissera et al. 2014). Meanwhile, niche saturation and limitations imposed by precipitation‐related factors intensify interspecific competition (Figure 7). There is the exclusion or elimination of some species by dominant ones (Selaya and Anten 2010), which in turn causes a decline in both TD and PD (Figure 3). During succession, gradual increases in species adaptation (Isbell et al. 2009) and niche differentiation (MacArthur and Levins 1967) lead to gradual increases in community evenness and stability (Figure 3e,f; Naeem and Li 1997).

Standardized α and β diversity in taxonomy and phylogeny of PKF during succession. PD, phylogenetic diversity; PβD, phylogenetic β diversity; TD, taxonomic diversity; TβD, taxonomic β diversity.

Further analysis revealed a synchronous time‐lag effect between TD and PD during succession (Figure 10a). Regarding β diversity, an asynchronous pattern was observed between TβD and PβD (Figure 10b), revealing a decoupled response of ecological and evolutionary processes to driving factors throughout succession. Successional changes directly affect TD (Figures 5 and 11). However, the existence of redundant species and the emergence of new lineages may strongly impact the phylogenetic structure of the region. These new lineages are the first to adapt and expand into new environments without immediately altering TD (Cavender‐Bares et al. 2006). With niche differentiation or the earlier filling of niche gaps (Cadotte et al. 2012), the available niches are rapidly occupied, leading to earlier changes in PD (Wang, Lv, et al. 2020, 2022). In contrast, TD reflects the number of coexisting species, which typically changes on an ecological timescale as it depends on the outcomes of ecological interactions and competition among species. This dependence usually results in TD changes lagging behind PD changes (Tilman 1994; HilleRisLambers et al. 2012). In terms of β diversity, TβD fluctuates more frequently, but many species replacements occur within closely related lineages, resulting in smaller changes in PβD (Graham and Fine 2008). Species turnover intensifies during the middle to late successional stages, driving the observed monotonic increase in TβD (Chase 2007). Generally, the regional species pool remains relatively stable, especially within the same plant community type (such as coniferous forests or grasslands). Even with species replacement, many of the newcomers still belong to closely related lineages, which constrains fluctuations in PβD (Cavender‐Bares et al. 2004). Therefore, the cosine fluctuation pattern of PβD reflects both the size and the evolutionary stability of the regional species pool in PKF (Figure 10b).

Response of taxonomic diversity (TD) and phylogenetic diversity (PD) to dominant environmental factors at different successional stages of PKF.

Inferences of Community Assembly in Succession Processes and Its Drivers

4.2

Given that phylogenetic clustering and overdispersion are respectively explained by environmental filtering and competitive exclusion (Emerson and Gillespie 2008), we infer that PKF's community assembly is primarily driven by competitive exclusion in the early stage, by environmental filtering in the mid‐successional stage, and again by competitive exclusion in the later stage (Figure 3i). This pattern extends the traditional framework of linear and unidirectional community assembly mechanisms in succession theory, revealing that environmental filtering and competitive exclusion are not simply alternative processes, but rather interact dynamically with changes in stand structure, resource availability, and evolutionary constraints. In the early successional stage, light transmittance under the forest is high, and the low species diversity may be due to the competition for light (Li, Peng, and Dang 2013). In the near‐mature forest, environmental conditions tend to be stable, and this stable environment strengthens the sorting for species with specific traits or lineages. On the other hand, both TD and PD reach their peak, suggesting the niche may be saturated. Therefore, environmental filtering plays a dominant role in shaping the phylogenetic structure (Figure 3i). In addition, the “intermediate disturbance” hypothesis may also play a role (Figure 5; Jiang et al. 2017). As succession progresses to the later stages, the community structure tends to stabilize (Figure 3f). High niche overlap leads to intense utilization of resources and increased competition pressure among species within communities (Webb et al. 2002; Kraft et al. 2008; Mayfield and Levine 2010; Fu et al. 2023). This pattern of community assembly during succession has been observed in the coniferous and broadleaved mixed forest of Changbai Mountain (Hou et al. 2017) and in tropical rain forests (Mo et al. 2013). However, on the Loess Plateau, environmental filtering is dominant in the later successional stages (Chai et al. 2016). This difference may be related to the aridity of the environment, specifically, a strong dependence on water availability.

Temperature‐related factors (bio9) directly affect variation in species composition during succession in PKF (Figures 6b and 8b). Temperature directly determines the colonization, growth, maturation, and survival of species (Sánchez‐González and López‐Mata 2005), greatly resulting in different species composition across regions or at different successional stages (Chen and Lewis 2024). P. kesiya var. langbianensis is a warm‐adapted coniferous tree species with a short evolutionary history (Gu and Li 1982); the mean annual temperature in its distribution area is 17.5°C–20.1°C, consistent with its biological characteristics (Li, Su, et al. 2013). In addition, the human footprint also plays a significant role in the successional processes (Figures 6b and 8b), which may be related to increased habitat fragmentation and the colonization of invasive alien species (Hoffmeister et al. 2005). PD reflects deeper evolutionary relationships shaped by paleoclimate. Sporopollenin records from the Ice Age indicate that global cooling was generally accompanied by decreased precipitation (Bartlein et al. 2011). However, because of the unique geographical location of Yunnan Province and its complex topography, the Quaternary Ice Age had less impact on the vegetation in this area (Bartlein et al. 2011; Tang et al. 2018). With the continuous uplift of the Qinghai‐Tibet Plateau in the Neogene and the gradual strengthening of the southern Asia monsoon climate, P. prekesiya began to differentiate into Pinus yunnanensis, which is gradually adapted to a more arid environment, and P. kesiya var. langbianensis, which persisted in more humid habitats (Xing et al. 2010). Precipitation‐dependent factors (bio19) continue to play a role in shaping the phylogenetic structure of PKF during contemporary succession (Figure 8d). This decoupling mechanistically explains how contemporary succession is simultaneously constrained by both immediate abiotic conditions and deep evolutionary history.

Explanations of Taxonomic Diversity and Phylogenetic Diversity on Community Stability

4.3

Pine forests are inherently pioneer communities that will eventually be replaced by broadleaved forests (Anyomi et al. 2022). Yet, they have persisted for millions of years, owing largely to periodic disturbances that reset successional trajectories (Maleki et al. 2019). Traditionally, community stability has been explained by TD alone (MacArthur 1955; Bezemer and Van Der Putten 2007), but this view may underestimate ecosystem complexity (Santos et al. 2014; Liu, Ma, et al. 2021; Liu, Zhang, et al. 2021; Wei et al. 2024). During the succession process of the PKF communities, the explanatory power of TD and PD for community stability exhibits a dynamic pattern of alternating dominance (Figure 12). Redundancy in TD and PD within a community contributes to stability across different levels (Sadeghinia et al. 2023). TD enhances communities' short‐term functional stability through niche differentiation (Liu et al. 2023), whereas PD confers long‐term adaptability and stress resistance via evolutionary diversity (Sadeghinia et al. 2023). Their complementarity enables the community to maintain strong stability facing environmental changes at different time scales (Symstad et al. 2003). This multidimensional stability mechanism during succession may provide a novel perspective on the long‐term persistence of pioneer pine forests.

Relative explanatory ability of taxonomic diversity (TD) and phylogenetic diversity (PD) to community stability in PKF's succession.

Relative Contribution by Spatial Variables and Environmental Variables

4.4

The significance of the spatial domain in ecological theory has been widely acknowledged (Legendre 1993; Borcard et al. 2004; Wu et al. 2021; Virgilio et al. 2022). Many ecologists mostly use the classic theory of SFTS in succession study (Parrish and Bazzaz 1982; Yuan et al. 2023). However, the spatial effects themselves are often overlooked. In this study, we divided the explanatory variables into two categories: environmental variables and spatial variables, and explored their effects on community structure from both species and phylogenetic perspectives during succession (He et al. 2022). Compared with spatial variables, environmental variables exert relatively less influence on species composition (Figure 9). However, as succession progresses, the influence of environmental variables gradually increases (Figure S5). In the early community, niches have not been fully occupied, random processes strongly influence community composition, and the environmental filtering effect is weak (Myers and Harms 2009). Species colonization is limited by the dispersal of pioneer species with close phylogenetic relationships (Weiher et al. 2011). The community stabilizes as succession progresses, and niches become saturated; then resource competitions are intensified. Species composition is increasingly determined by environmental adaptation, acting more strongly on specific evolutionary lineages (Cadotte and Tucker 2017). This pattern may reflect the “successional continuum hypothesis” (Gravel et al. 2006), which proposes that at the regional scale, both neutral processes and niche differentiation jointly shape community assembly, forming a continuum rather than a binary opposition. In addition, it is important to note that the unexplained variance at some succession stages in our study is high, underscoring the complexity of ecological processes (Legendre et al. 2009).

Limitation of the SFTS Method in Succession Studies

4.5

Although the SFTS method is widely used in many ecological studies, especially in studies of community succession (Blois et al. 2013; Kharouba and Williams 2024), this approach still has important limitations and should be applied with caution when interpreting results. Firstly, community spatial patterns may be shaped by both temporal and spatial processes (such as dispersal limitation and local interspecific interactions). Therefore, the SFTS approach may conflate temporal and spatial processes, making it difficult to determine whether observed community variation stems primarily from spatial structure or from temporal succession. To address this, our study introduces the “dbmem” function to disentangle temporal and spatial influences in the successional process (Figure 9). Secondly, SFTS often ignores the possible historical effects of each plot or community, such as the legacy effects of disturbance events (e.g., fire, human footprint) on the community structure. These historical events may not be correlated with temporal succession patterns along spatial gradients. Therefore, we explicitly incorporated two disturbance‐related variables into the analysis to minimize these potential sources of error. Our findings highlight that spatial variation should not be overlooked when applying the SFTS approach. However, the five successional stages formed a clear spatial continuum, and the patterns of community and environmental variation were highly consistent with theoretical expectations. Thereby, the applicability of the SFTS approach is supported in this context. Moreover, we carefully probed the inferred successional dynamics, acknowledging that our data represent a space‐for‐time comparison rather than a directly observed temporal sequence. While accounting for spatial connectivity and incorporating more detailed environmental variables and historical processes can help minimize biases, these measures may remain insufficient. The successional configuration of pine forests by SFTS will be helpful to conduct further validation and in‐depth analysis using long‐term time series data.

Suggestions on Management and Conservation of PKF

4.6

The resin secreted by P. kesiya var. langbianensis contains a variety of bioactive compounds such as monoterpenes, diterpenoids, flavonoids, and phenolic acids, which have antibacterial, anti‐inflammatory, antiviral, and antitumor activities (Liu et al. 2022). Moreover, as a fast‐growing tree species, it is highly valued for papermaking and timber production (Li, Peng, and Dang 2013). However, sustainable logging and resin production depend on the healthy development of PKF. Our results indicate that species richness is low and community stability is poor during early succession, and community assembly is dominated by competitive exclusion. So, forest management should focus on reducing competitive stress at early stages to promote community development. For example, community stability may be enhanced by selective thinning or removal of dominant individuals. At the middle stage, species diversity and stability reach higher levels, and environmental filtering plays a dominant role. Because this stage also coincides with the peak period for resin extraction, a regulated harvesting practice participated in by local stakeholders is necessary to maintain both high diversity and resource sustainability. The increasing diversity and recovery of PD at the late stage suggest that management should promote natural regeneration and, if necessary, artificially plant late‐successional species to maintain ecosystem resilience. Across all stages, establishing a long‐term monitoring system is crucial to link management practices directly to successional dynamics and to strengthen local community awareness of rational resin extraction and forestry resource protection. As a regionally important ecosystem, the protection of PKF is critical not only for biodiversity and ecosystem functioning but also for the local economy and community livelihoods. Overall, management practices integrated with successional patterns will ensure forest protection and sustainable usage.

Conclusion

5

This study integrates ecological and evolutionary perspectives and confirms the classic hump‐shaped pattern with diversity peaking in the mid‐successional stage. More importantly, we propose two novel patterns: a “synchronous time‐lag effect” between TD and PD and an “asynchronous oscillation” between TβD and PβD. The former reveals that ecological and evolutionary processes do not respond synchronously through time, whereas the latter shows that TβD increases monotonically during succession, whereas PβD exhibits a low‐amplitude cosine‐like fluctuation. These dual‐track patterns in subtropical pine forests provide new insights into the interplay between ecological filtering and evolutionary constraints. Next, we identify a periodic shift in the dominant community assembly mechanisms, that is, from early competitive exclusion to mid‐successional environmental filtering, and then back to competitive exclusion in the late stage. This oscillatory pattern expands the traditional view of succession as a linear transition and demonstrates that biotic and abiotic processes alternate repeatedly across temporal scales and jointly shape community structure. Furthermore, we reveal a multidimensional stability mechanism, suggesting that the long‐term persistence of pioneer communities may rely on the functional complementarity between taxonomic redundancy and phylogenetic potential pool across time. At the ecological scale, bio9 acts as an environmental filter shaping taxonomic composition, whereas at the evolutionary scale, bio19 influences the regional species pool and thus determines phylogenetic structure. Finally, our analytical framework explicitly quantifies the inherent spatial effects within the SFTS approach, demonstrating that spatial structure plays an essential role in explaining temporal variation. Long‐term monitoring data are still necessary to uncover the ecological–evolutionary differentiation and mechanisms driving succession.

Overall, by linking ecological processes with evolutionary constraints, the general understanding of pine forest succession is enriched theoretically, as well as the sustainable management of pine forest ecosystems is advanced by stage‐specific strategies.

Author Contributions

Xiaofan Wang: conceptualization (equal), data curation (equal), formal analysis (equal), investigation (equal), visualization (lead), writing – original draft (lead). Yunfei Ma: data curation (equal), investigation (equal). Biao Zhao: data curation (equal), investigation (equal). Dengpeng Chen: data curation (equal), investigation (equal). Yehong Luo: data curation (equal). Mingchun Peng: investigation (equal), writing – review and editing (equal). Yongping Li: investigation (equal), writing – review and editing (equal). Xinmao Zhou: investigation (equal), writing – review and editing (equal). Wen Chen: writing – review and editing (equal). Cindy Q. Tang: funding acquisition (equal), writing – review and editing (equal). Chongyun Wang: data curation (equal), funding acquisition (equal), investigation (equal), methodology (equal), writing – review and editing (equal).

Funding

This study received financial support from the Major Program for Basic Research Project of Yunnan Province, China (202101bc070002), and the Special Foundation for National Science and Technology Basic Resources Investigation of China (2019FY202302).

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Figure S1: Relationships between field‐recorded disturbance levels and remote‐sensing‐based disturbance indicators. Figure S2: Correlation between environmental factors and spatial factors. Figure S3: Phylogenetic tree of all species in Pinus kesiya var. langbianensis forest (PKF). Figure S4: Stand structural characteristics of PKF across five successional stages. Figure S5: Variation trend of explanatory rate of spatial factors and environmental factors in PKF succession processes. Table S1: Field investigation plots of Pinus kesiya var. langbianensis forest (PKF). Table S2: 16 MEM principal component axes with significant positive eigenvalues. Table S3: List of screened explanatory variables. Table S4: Results of DCA ordination parameters. Table S5: Mantel test of taxonomic β diversity (TβD) and phylogenetic β diversity (PβD) at different successional stages of PKF.

Data S1: ece372911‐sup‐0002‐DataS1.xlsx.