The study area encompasses the entirety of Northeastern China (Figure 2A) located between 40–53°N latitude and 120–135°E longitude. It includes the Heilongjiang, Jilin, and Liaoning provinces, as well as the eastern part of the Inner Mongolia Autonomous Region (Zhang and Liang, 2014). The forests in this area are predominantly found in three regions: Changbai Mountain (CB), Xiao Hingan Ling (XH), and Da Hingan Ling (DH) (Zhang and Liang, 2014) (Figure 2B). The Northeastern region of China (NE China) encompasses a forest area of more than 82.42 million ha² (Figure 2B). It encompasses a range of forest types, from temperate forests in the south to boreal forests in the far north, representing all of the major forest types in temperate East Asia (Zhang and Liang, 2014). We obtained 1 km-resolution spatial data of the area representing mean annual temperature (MAT), mean annual precipitation (MAP), and mean annual potential evapotranspiration (PET) for the period of 2000–2020 from the National Earth System Science Data Center, National Science & Technology Infrastructure of China dataset.¹ Data of soil physicochemical properties including pH value (pH), soil organic carbon (Soc), total nitrogen (Tn), total phosphorus (Tp), total potassium (Tk), bulk density (Bd), and coarse fragments (Cf) at various depths (0–5, 5–15, 15–30, 30–60, 60–100, and 100–200 cm) were obtained from the study by Liu et al. (2021) (90 m spatial resolution; available at National Earth System Science Data Center, http://www.geodata.cn/). We applied weightage to synthesize the soil data from each layer (0–5, 5–15, 15–30, 30–60, 60–100, 100–200 cm), considering that layer weight is equal to layer depth/200 cm. Finally, we summed the weighted values of each layer to generate the soil physicochemical property values of the sample plot.

This study utilized publicly available data from the National Ecology Science Data Center, National Science & Technology Infrastructure of China (NESDC, http://www.geodata.cn). This database comprises 122 sample plots of different forest types and disturbance classes that were established in Northeast China between 2015 and 2018 (Figure 2A). The forest layer with height (H) ≥ 1.3 m is designated as the tree layer, within which the species, status, diameter at breast height (DBH), and height of each tree are recorded in 36 quadrats of 5 × 5 m each. Ten quadrats were selected as survey plots to investigate understory wood plant diversity. The species, DBH, and height of each plant in the understory wood plant layer are investigated and recorded. We randomly set a 1 × 1 m understory herbaceous vegetation survey plot within the understory wood plant survey plot. The types, average height, coverage, and other relevant details of the understory herbaceous vegetation were recorded following the guidelines outlined in LY/T 3128–2019 (“The Regulations for Classification, Survey, and Mapping of Forest Vegetation, Forestry Industry Standard of the People’s Republic of China”).

Soil fauna were identified from three sieved soil core samples of 5 cm diameter each. Soil cores were taken from the upper 15 cm of soil in each plot, and the three cores were pooled as a composite sample for the plot (NESDC, http://www.geodata.cn). The samples collected in the field were placed on a sieve equipped with a funnel. The lower part of the funnel was connected to a device containing 75% alcohol by mass to collect and preserve the separated soil fauna. The separation time was 7 to 12 days. The device containing the soil sample was then poured into another glass bottle and heated in a water bath (Wu et al., 2021). After boiling the alcohol in the bottle, the suspended springtail sample sank and was removed for soil fauna sorting. Soil fauna were sorted by morphology, and the number of individuals was recorded before storing in centrifuge tubes containing 75% alcohol (Chen, 1983). The separated soil fauna were made into slide specimens and identified following slide specimen preparation processes such as fading, washing, and making (Xu et al., 2017). Species morphological characteristics were documented with photographs that were captured using a field microscope with high depth. Finally, two types and four trophic levels of soil fauna were collected from the soil under different forest types. The two types included large soil fauna and small and medium-sized soil fauna (such as Nematodes, Oribatida, and Collembola).

Soil fungi and bacteria were identified from sieved and freeze-dried soil core samples. Ten cores of 3 cm diameter each were taken from the upper 10 cm of soil per plot and pooled as the composite sample of the plot. The samples were then passed through a 2-mm soil sieve, placed into a ziplock bag, and cooled in a − 20°C refrigerator (NESDC, http://www.geodata.cn, Li P. et al., 2019). The PowerSoil DNA Isolation Kit (MoBio, United States) was used to extract total soil fungi DNA. The primers ITSF1-GCATCGATGAAGAACGCAG and ITSR1-TCCTCCGCTTATTGATATGC were employed to amplify the hypervariable region of fungal ITS2. The sequencing of paired-end 250 bp was conducted using the Illumina HiSeq platform (Li et al., 2020). The barcoded Illumina paired-end sequencing (BIPES) process was then employed for initial data processing to reverse the complementary assembly of the ITS2 sequence. The sequence was split into corresponding samples based on barcodes, and low-quality data such as primer mismatch, N-containing sequences, and excessively short sequences were eliminated. The UCHIME program was employed to screen out chimera sequences, and diversity analyses, such as clustering, were performed using the QIIME process. Operational taxonomic unit (OTU) clustering was conducted using the UCLUST method with a sequence similarity threshold of 97% (Li P. et al., 2019). The representative sequence for each OTU was selected using the highest frequency, and systematic classification was performed using the Amazon Relational Database Service (RDP) and UNITE databases.

The primers 515F/806R were similarly used to amplify the V4 hypervariable region of the bacterial 16S rRNA gene (Tamaki et al., 2011), and the above steps were repeated, except the SILVA database which was used to perform the systematic classification (Li et al., 2020). All measurements for trees, understory vegetation, and soil biodiversity reported in this study were taken from the same sample plot.

Tree and understory plant species diversity for each plot was determined by counting the unique species present in the tree and understory wood and herbaceous layers within each sample plot. Varieties within each species were aggregated (Zhou et al., 2017). In this study, we defined eight stand structure features, including β diversity of trees (β-Tree), stem abundance (SA), mean DBH (m-DBH), tree size variation (CV DBH), variation of height on crown base (CV HCB), height of crown (CH), stand basal area (SBA), basal area of deadwood (BAD), and the stand structural complexity index (SSCI). β-Tree measures the variation in species composition among sites, reflecting species turnover and nestedness (Baselga and Orme, 2012) and is calculated using the Jaccard dissimilarity index following the method stated by Zhang et al. (2022); SA represents stand density and is measured in individuals ha⁻¹; m-DBH is a proxy for stand age (Penone et al., 2019); CV DBH is calculated as the coefficient of variation for tree DBH within a given sample plot (Chu et al., 2019); CV HCB is used to assess the physical space under the tree canopy; CH is measured using the method stated by Liu X. Q. et al., 2022; SBA is calculated as the sum of individual stem basal area (Dolezal et al., 2020); BAD is the sum of deadwood basal area; and the aim of the SSCI to quantify the distribution of trees and their canopies in three-dimensional space as a function of leaf area and heterogeneity of biomass distribution (Ehbrecht et al., 2021). Each feature was calculated as the mean or sum of all individuals in the sample plot (Table 1) (Krebs et al., 2019).

Temporal stability of forest productivity is generally quantified as the ratio of the temporal mean of productivity to its standard deviation (SD; Tilman et al., 2006; Hautier et al., 2015; Craven et al., 2018; Liu S. et al., 2022). In our study, we used the peak Normalized Difference Vegetation Index (NDVI) from the period 2000–2020 as a proxy for aboveground plant biomass (following Yang et al. (2019) at the National Ecology Science Data Center of China²). The NDVI data had a spatial resolution of 30 × 30 m, which aligned with the area of our field survey sites. We calculated ecosystem temporal stability as the ratio of the mean NDVI from 2000–2020 to the SD of NDVI (Wu et al., 2023) (Figure 2B). The NDVI time scale is consistent with the climate data time scale.

We employed a piecewise structural equation model (SEM) as a linear mixed-effects model (Li Z. L. et al., 2019; Jin et al., 2022) to evaluate the direct and indirect effects of environmental factors and stand structure features on forest productivity stability, considering the effects of climate gradients as random. Response variables were summarized to the plot level for analysis. We then considered various alternative reduced models that shared the same causal structure as the initial model and were constructed by eliminating non-significant variables one at a time. The decision to remove a path was based on the performance of the overall model fit and the p-value for the path (p > 0.05) (Wang et al., 2023). The models were evaluated according to the following two criteria: (1) pathway significance (p < 0.05) and satisfactory fit (p > 0.05) and (2) the goodness of fit of the model (chi-square test (χ²) (p > 0.05) and 0 ≤ Fisher’s C/df ≤ 2). Finally, we removed non-significant paths with p > 0.05 in SEMs with satisfactory model fit and then reassessed the model fit (Chu et al., 2019).

Due to the high correlations between climate variables (MAT, MAP, and PET) (r = 0.782–0.945) and soil physicochemical properties (Soc, Tn, and Tk) (r = 0.689–0.856), we conducted principal component analyses (PCA) (Chu et al., 2019; Chen et al., 2021) and used the score of the first axis to represent climate (explaining 90.0% of the variation) and the scores of the first three axes to represent soil factors in each plot (explaining 77.4% of the variation) (Supplementary Figures S1, S2). We used separate analyses to account for collinearity among environmental variables and stand structure factors (Supplementary Figures S3–S6). We removed highly correlated variables (Pearson r > 0.60) and then performed principal component analysis (PCA) (Chu et al., 2019). To make model coefficients comparable in multi-predictor regressions, we normalized the response variables and all the predictors using the z-score (overall mean of 0 and SD of 1) before analyses (Penone et al., 2019; Chen et al., 2021; Zhang et al., 2022). To account for diversity–stability variation with multiple-taxon of relative importance, we calculated the proportions of variance components in a mixed-effects model as follows: lme (Stability – Climate_PCA1 + Soil_PCA1 + Soil_PCA2 + Soil_PCA3 + α-Tree + α-Shrub + α-Herbs + α-Fauna + α-Fungi + α-Bacteria + Stand structure, random = ~ 1|site regions (see Part 2.1)). Using R package glmm.hp (Lai et al., 2022), we calculated the contributions of these predictors (Liu C. et al., 2023), and the relative importance of the predictors was grouped into four identifiable variance fractions: climatic, soil, biodiversity, and stand structure (Table 1). Following the approach in a previous study (Carol Adair et al., 2018), we computed the gain or loss in biomass carbon and soil organic carbon stocks (via direct and indirect effects) by multiplying the unstandardized effect size by the range of that variable (Table 1) (Chen et al., 2023).

The analyses were carried out in RStudio 4.2.3 (Rstudio Inc) (R Core Team, 2023). The SEM analysis was implemented using the ‘piecewiseSEM 2.1.0’ package (Lefcheck, 2016), and LMM analysis was implemented using the ‘nlme 3.1–162’ package (Pinheiro et al., 2021) and ‘glmm. hp. 0.0–9’ package (Lai et al., 2022). All results were visualized with the ‘ggplot2’ package (Wickham, 2016).

Determining the contributions of abiotic and biotic factors to productivity stability can greatly improve current forest management and conservation strategies (Nabuurs et al., 2013). Our study provides the first quantitative assessment of the influence of understory vegetation and soil biodiversity on the stability of productivity in Chinese temperate forests. When considering both understory plant and belowground soil biota biodiversity, the best-performing model that includes soil and climatic variables accounts for 92% of the variation in forest productivity stability (Figure 3B). This value is larger than those in several previous large-scale studies that evaluated forest productivity stability using remote sensing to analyze grassland ecosystems (73% of variation in García-Palacios et al., 2018, and 43% in Chen et al., 2021) and studies that applied allometric modeling approaches to analyze China’s temperate and subtropical forest ecosystems (14% of variation in Qiao et al., 2022, and 50% of variation in Ouyang et al., 2021). Several previous studies have observed strong associations between tree species diversity and productivity stability in forest ecosystems (Chen et al., 2021; Qiao et al., 2022, 2023). Many studies have also shown that tree species diversity has a positive effect on productivity stability through a combination of processes such as species asynchrony and interactions (Jucker et al., 2014; Ouyang et al., 2021). Our study reveals that the influence of tree species diversity on productivity stability is mediated by stand structure. These findings are consistent with the study by Dolezal et al. (2020) who found that stand structural attributes significantly enhanced community productivity stability at small scales using stem increment data.

However, the results of this study do not indicate significant direct effects of understory woody and herbaceous plant, soil fungi, and bacteria species diversity on forest productivity stability.

A counterview of the biodiversity-stability relationship holds that dominant species play a major role in driving forest productivity stability (Grman et al., 2010). The stability of dominant species populations can strongly influence community biomass stability, especially in communities dominated by a small number of species (Ma et al., 2017). Although understory vegetation and soil biota exhibit high species diversity, it is worth noting that most of the understory plant species that appeared in response to our experimental treatments were classified as “rare,” contributing only a small fraction to community biomass stability compared to dominant species (Ma et al., 2017). In Northeast China, the biomass of understory vegetation ranges from 2.76 to 9.70 t∙ha⁻¹ (Jin et al., 2022), accounting for only 5–11% of the biomass of the entire forest ecosystem, the mean value of which is 85.24 t∙ha⁻¹ (Su et al., 2016). Large studies have shown that soil biota diversity and plant diversity both contribute to productivity stability in grasslands (Wagg et al., 2014; Chen et al., 2021). Soil biota affected productivity stability via multiple direct and indirect pathways through soil nutrient availability or plant species richness. In some cases, the association between productivity stability and soil biota species diversity is equivalent to that between productivity stability and plant species diversity in grassland ecosystems (Wu et al., 2023). However, the strength and direction of these effects vary between soil biota groups (Chen et al., 2021) and ecosystem types (Gross et al., 2014). One aspect, in particular, is likely to be especially relevant for determining which processes drive productivity stability: the fact that grassland community species change in relative abundance from year to year, whereas shifts in community composition occur much more slowly in forests (Jucker et al., 2014). Furthermore, there are strongly positive linkages between plant species diversity and soil biota species at local scales (Yang et al., 2014, 2016). Our research reveals that soil fauna influence the productivity stability of forest ecosystems through their interactions with tree species diversity at a large scale.

Our results indicate that CH (34%, p < 0.001), CV HCB (10%, p < 0.05), and SBA (7%, p < 0.01) are the most important drivers of productivity stability. These findings align with García-Palacios et al. (2018) who found that mean height was the most important variable for predicting productivity stability. Across the global network of drylands analyzed in the study, the dominance of particular species was associated with productivity stability, as indicated by the significant quadratic effect of mean plant height (García-Palacios et al., 2018). There are several reasons for this. First, high tree richness (α-diversity) promotes stand density and size difference due to species characteristics. These stand structure factors are critical component of ecosystem resilience, or the capacity of an ecosystem to persist through and re-organize after disturbance, adapt to shifting environmental conditions, and maintain basic structure and function over time (Gonzalez and Loreau, 2009; Churchill et al., 2013). Furthermore, these variables further promote CH, SBA, and CV HCB. Second, higher turnover in tree species diversity (β-diversity) enhances resilience by increasing the variation in canopy height and spatial asynchrony among communities, thereby stabilizing productivity according to the spatial insurance hypothesis (Hautier et al., 2020). Our results also indicate that stand basal area (SBA) and soil property (Tn) negatively affect productivity stability. Higher SBA may promote intraspecific competition and dieback, while higher Tn increases vegetation growth rate. In considering multiple environmental factors together, we found that climate is the ultimate driver of productivity stability in our models because it promotes tree species diversity. Variations in precipitation and temperature (Knapp and Smith, 2001; Dobrowski et al., 2013) create fine-scale mosaics of environmental conditions (Kane et al., 2019) and affect tree species diversity. Various tree species also contribute different heights, crown architecture, and structural complexity, which directly influence the ecological niche and micro-climate within forest communities. Climate factors not only control diversity and productivity but also influence the capacity for diversity to stabilize ecosystem function by altering the mechanisms that link diversity-structure-stability. These findings highlight the crucial roles of climate and diversity in determining forest ecosystem stability. However, it must be acknowledged that we did not calculate the gains of forest productivity stability brought about by changes in climate factors.