Driving mechanisms of productivity stability vary with selective harvesting intensities in a mixed broad-leaved Korean pine forest

Abstract

Key message

We found that the stabilizing mechanisms for forest productivity varied across harvesting intensities in a mixed broad-leaved Korean pine forest. Effects of overyielding at high species richness and species asynchrony occurred only in unharvested and lightly harvested plots, whereas asymmetries between individuals of different size contributed significantly to stabilizing productivity when harvestings became intensive.

Context

Understanding the driving factors of forest ecosystem stability has become increasingly crucial in forest management. However, it remains unclear whether and how the stabilizing mechanisms of forest productivity might be influenced by management practices.

Aims

We related the temporal stability of aboveground biomass productivity to harvesting history. We further tested how three key driving mechanisms of stability might be modulated by selective harvesting intensities.

Methods

Based on a 10-year monitoring (five repeated tree inventories) of a mixed broad-leaved Korean pine forest in Northeastern China recovering from selective harvesting, we examined the relative importance of two diversity-dependent mechanisms (overyielding and species asynchrony) and one size-dependent mechanism (asymmetric growth) for productivity stability across a wide range of intensities (0–73.4% basal area removed).

Results

We found that selective harvesting significantly lowered the productivity stability, species asynchrony, and growth dominance coefficient. Growth dominance coefficient had an overall stronger effect on stability than species richness and asynchrony. Moreover, the strengths of stabilizing mechanisms varied across harvesting intensities: effects of overyielding at high species richness and species asynchrony were detected only in unharvested and lightly harvested plots, whereas the explanatory power of growth dominance coefficient outweighed the diversity-related variables when harvesting became intensive.

Conclusions

We emphasized the importance to consider both diversity- and size-related explanatory variables as potential mechanisms for the temporal stability of forest productivity. In fact, how growth is partitioned among trees of different species as well as sizes may co-determine the response of forest stability to disturbances.

1 Introduction

Forests play an essential role in hosting global terrestrial biodiversity and providing a wide range of ecosystem functions such as timber production and carbon sequestration. This role, however, can be critically modified by ongoing climate change and other disturbances (e.g., anthropogenic activities) (Grantham et al. 2020; Seidl et al. 2017). As one of the basic attributes of ecosystems, stability demonstrates the capacity of an ecosystem to defy change. One of the most frequently studied aspects of stability, temporal stability, measures the constancy of ecosystem functions, such as biomass production, over time (Tilman 1999). The temporal stability of biomass production is usually defined as the inverse of its coefficient of variation, i.e., the ratio between mean productivity (μ) and its variation in time (standard deviation, σ). Factors that affect μ and/or σ may modify the temporal stability of biomass production (Isbell et al. 2009).

Over the past two decades, temporal stability of biomass production has been frequently linked to biodiversity (e.g., Craven et al. 2018; De Boeck et al. 2018). Theoretically, species diversity may promote stability by increasing performance (increasing μ) or buffering variation (decreasing σ), corresponding to two mechanisms: (i) overyielding due to complementarity among species and (ii) asynchronous responses of species to environmental conditions (Hector et al. 2010). Compared with the “overyielding effect,” the buffering effect is much less explored, particularly in forests, because of the long lifespan of trees and their slow adjustment to environmental changes. However, asynchronous interannual dynamics in species productivity have been found to drive the stabilizing effect of diversity in some forests (Aussenac et al. 2019; Jucker et al. 2014). Tree species with different adaptation strategies can compensate for each other, a function of compensatory growth in the face of fluctuating environments (Cordonnier et al. 2018; del Río et al. 2017). While the important role of species asynchrony in forest stability has been recently highlighted, most evidences came from forest plantations that were relevant only to young development stages without seedling recruitment (Schnabel et al. 2021). Studies testing species asynchrony in natural forests, on the other hand, had the limitation of low tree species diversity, such as in western European temperate forests (Jourdan et al. 2021).

Disturbance history can have a profound impact on forest structure and functions, either directly by increasing mortality of individuals (or species) or indirectly by changing resource availability and heterogeneity. Selective harvesting represents a major disturbance characterized by the removal of single-tree or small groups. The opening of canopy gaps improves the tree species’ access to resources (e.g., light). As light availability becomes less limiting, the magnitude of species complementarity effect for light may decline (Lebourgeois et al. 2013; Toïgo et al. 2015). Besides, lower stand densities may weaken species interactions that trigger complementarity, because individuals are further apart and thus interact less (Feng et al. 2022; Forrester 2014). For these reasons, a decreasing trend of species complementarity due to harvesting may be expected (Fig. 1a). Such changes in species complementarity following harvesting could further modify the overyielding effect (affecting μ).

Fig. 1

Illustration of how selective harvesting may influence two diversity-dependent (a, b) and a size-dependent (c) regulatory mechanisms of the temporal stability of aboveground biomass productivity in multi-layered mixed forests. The width of the top arrows in each graph reflects the magnitude of the effects. In part a, along the gradient of increasing harvesting intensity (decreasing stand density), the magnitude of species complementarity for biomass production may show a decreasing (dashed blue line). These changes of complementarity further alter the overyielding effect. Part (b) shows how species with different shade tolerances and growth strategies respond to harvesting treatment. The solid lines indicate increased abundance of a pioneer (yellow) and a subordinate light-demanding (blue) tree species, while the dashed line indicates a decrease in the abundance of subcanopy shade-tolerate species (orange). The canopy dominant trees (green) do not show strong reactions (dotted green line). Such changes in community composition may modify the overall asynchronous species dynamics at the community level. Red X indicates tree mortality. Part c shows that selective harvesting reduces the degree of size-asymmetric resource competitions, benefiting the growth of small-sized trees (solid line) over large trees (dashed line), as indicated by a decline in the growth dominance coefficient (GDC)

Secondly, harvesting may influence stability via changing the degree of species asynchrony (affecting σ). In forest ecosystems, light availability is usually the key limiting abiotic factor for tree growth. The ability to optimize light capture contributes to the asynchrony of species dynamics via a diversity of light adaptation strategies. However, an increase in light availability in the understory may reduce the environmental heterogeneity and hence the diversity of tree shade tolerances (Dolezal et al. 2020). This lower level of niche segregation (negative species covariations) due to enhanced resource availability has been reported in grassland experiments, where nitrogen enrichment often reduced species asynchrony and led to decreased ecosystem stability (Lepš et al. 2018; Liu et al. 2019; Zhang et al. 2016). Moreover, disturbances tend to benefit more pioneer and subordinate light-demanding tree species that are temporarily competitively superior over slow-growing, shade-tolerant trees (Yano et al. 2021). However, these species with high growth rate generally display high temporal variability and similar responses to extrinsic factors (e.g., weather fluctuations, pests, and disease) (Bazzaz 1979; Coomes et al. 2014), leading to overall lower stability at the community level due to less asynchronous population dynamics (Fig. 1b). This is also the major reason why young (second-growth) forest patches, dominated by fast-growing shade-intolerant species, usually exhibit lower ecosystem stability compared with the late forest successional stage where conservative strategies take over (Musavi et al. 2017; Yuan et al. 2019; Zhang et al. 2012).

Lastly, as argued by Loreau and Mazancourt (2013), most studies only considered the different environmental responses among species while ignored the strong effect of tree size on biomass production stability. Asymmetry, in addition to asynchrony, may also stabilize ecosystem functioning. For example, asymmetric competition between upper-canopy stratum (25–30 m) and lower-canopy stratum (5–15 m) is a major driver of productivity stability in an old-growth forest (Dolezal et al. 2020). Selective harvesting reduces the degree of size-asymmetric resource competitions, and hence, overyielding can be asymmetrical: small-sized trees should benefit more than large trees, as indicated by a decline in the growth dominance coefficient (GDC), which has been used as a quantitative method to evaluate the consequence of silvicultural treatments on stand structure (Binkley 2004; Bradford et al. 2010; Forrester 2019; Soares et al. 2017) (Fig. 1c). Nevertheless, young and short-statured trees have inherently more variable growth rates and higher mortality risk (Jucker et al. 2014). Moreover, an increase of small trees leads to greater competition for moisture, making all of the trees more vulnerable to drought and even causing death of old-growth trees (Giles et al. 2021). Therefore, an increase in the contributions made by small trees to overall stand production may destabilize stability at the whole-stand level.

The broad-leaved Korean pine (Pinus koraiensis) forest (BKPF) is the dominant vegetation type in eastern Eurasian Continent, characterized by high biodiversity and abundance of high-quality timbers. Close-to-nature, uneven-aged silvicultural systems such as selective harvesting have been widely applied in BKPF after decades of excessive logging and poorly regulated management (Dai et al. 2009; Zhang et al. 2014). Although selective harvesting generally aims at promoting the growth of remaining plants, trees do not respond in the same way to environmental changes. The first decade after harvest represents a reorganization phase: demographic processes such as growth, mortality, and recruitment may change dramatically following harvesting, with size- and species-specific responses depending on social status, life history, shade tolerance, morphological plasticity and colonizing ability (de Avila et al. 2017; Dionisio et al. 2018; Kuehne et al. 2015). This period is a critical window because changes in community attributes (e.g., species composition, size structure, the spread of invasive species, and death of tree seedlings) can come quickly and the accumulation of biomass may be modulated by harvest-induced disturbance as well. Such altered forest dynamics during this transitional phase determine the direction and magnitude of subsequent forest change and may shape forest structure and composition for decades, even centuries (Rammer et al. 2021; Seidl and Turner 2022). A better understanding of stability and its driving mechanisms during this relatively short period of time would help predict long-term consequences of selective harvesting and contribute to enhance the suitability of harvesting management, balanced with its impacts as a disturbance itself.

In order to examine how the temporal stability of forest productivity (hereafter productivity stability) relates to stand disturbance history, and to elucidate factors that confer greater stability, here we conducted a decade-long detailed monitoring of individual trees in an experimental BKPF following selectively harvesting. Harvesting intensities varied substantially among patches within the forest, from unmanaged to intensively harvested (over 70% of initial basal area removed). Here we explored three potential mechanisms that may drive forest productivity stability following selective harvesting: two diversity-dependent mechanisms (overyielding and species asynchrony) and one size-dependent mechanism (asymmetric growth of trees). Specifically, we hypothesize that (i) with enhanced resource (e.g., light) availability by selective harvesting, the strength of complementarity among species may decrease, making the overyielding effect less prominent (Fig. 1a); (ii) harvesting may also alter community composition by favoring pioneer and subordinate light-demanding tree species with more rapid and synchronous growth, leading to reduced overall species asynchrony and its effect on stability (Fig. 1b); and (iii) as intense asymmetric competition for light prevails in the BKPF, shifts in growth dominance between large and small trees following harvesting may also be an important driver of stability (Fig. 1c).

2 Material and methods

2.1 Site description

The study site is part of an experimental manipulation of forest management in the Jiaohe Forest Experimental Zone in Jilin Province, Northeastern China. Forest management practices including selective harvesting were applied to examine their influences on forest ecosystem components. The site is characterized by a temperate continental monsoon climate, with the hottest and coldest months being July (with an average temperature of 21.7 °C) and January (with an average temperature of − 18.6 °C), respectively. The mean annual temperature is 3.8 °C and mean annual precipitation is 695.9 mm. Severe drought or wet events did not occur during the monitoring period (2011–2021) (Appendix Fig. 8). The soils are classified as dark-brown forest soils (FAO 2006). The typical forest type in this study area is a late successional uneven-aged broad-leaved Korean pine forest (BKPF). Most canopy trees are now more than 150 years old (Wang et al. 2016). Dominant tree species included Pinus koraiensis, Fraxinus mandshurica, Tilia mandshurica, Tilia amurense, Acer mono, Juglans mandshurica, and some pioneer species from the Betulaceae and Salicaceae families. The last recorded massive harvesting (> 50% stock extracted) occurred 60 years ago (Wang et al. 2016). The forest has been restricted from commercial logging and other human disturbances since the launch of China’s Natural Forest Protection Project.

2.2 Experimental design

Four sites (experimental blocks) were chosen from the Jiaohe Forest Experimental Zone (21.12 ha, 660 × 320 m). Sites were selected to represent similar topographic positions and community composition (Appendix Table 2). At each site, a 1-ha (100 × 100 m) permanent plot (working unit) was established. Among the four permanent plots, one was set as control without any disturbance. During October–November 2011, selective harvesting was applied to the other three plots according to the guidelines of “close-to-nature silviculture” (Lu 2006). In brief, according to the “target tree silviculture method,” identification of target trees (residual trees) was made based on several aspects including tree health, vigor, rareness, commercial values, competition with neighboring trees, and historical or cultural importance (Lu 2006). Harvesting aimed to promote the growth and quality of those target trees while limiting the competition from nontarget trees. Nontarget trees selected for removing covered a broad range of size classes and were not biased for particular species. In order to efficiently monitor the dynamics of community composition and biomass, following a standard field protocol (Center for Tropical Forest Science, Condit 1998) we split each plot into 25 subplots (20 × 20 m each). There were thus a total 100 subplots across the four sites. For each subplot, we quantified separately the harvesting intensity as the total basal area of the harvested trees. Among subplots, the observed harvesting intensities ranged from 0 to 73% of basal area (BA) removed, representing various levels of harvest intensity. We then assigned all subplots to one of the four harvesting intensity levels: HI0 (no harvesting); HI1 (lightly harvested, 0 < BA removed ≤ 20%); HI2 (moderately harvested, 20% < BA removed ≤ 40%); HI3 (heavily harvested, BA removed > 40%).

2.3 Temporal stability of aboveground biomass productivity

All trees that had a diameter at breast height (DBH) exceeding 1 cm were tagged, mapped in the subplots and their DBHs were recorded. Measurements on the inventoried trees occurred before and right after the selective harvesting in 2011. Re-measurements were taken in 2013, 2015, 2018, and 2021. A total of 2618 individuals belonging to 14 families and 32 species were inventoried.

The aboveground biomass of individual trees were quantified using a set of region- and species-specific allometric models obtained in the same study site (Appendix Table 3) (He et al. 2018). Dynamics of aboveground biomass were monitored at the subplot level. Aboveground biomass stock of each subplot was the sum of biomass of all trees (DBH > 1 cm) within the subplot. Aboveground biomass productivity (Mg ha⁻¹ year⁻¹) was estimated as the change of aboveground biomass stock between adjacent inventories, divided by time (four inventory intervals: 2011–2013, 2013–2015, 2015–2018, 2018–2021). Temporal stability of aboveground biomass productivity in each subplot was defined as μ/σ, where μ and σ represented the temporal mean and standard deviation (SD) of productivity over the four intervals, respectively (Lehman and Tilman 2000). More functionally, we classified species into four ecological groups: pioneer (PI), light-demanding upper-canopy (LU), (non-pioneer) light-demanding lower-canopy (LL), and shade-tolerant (ST) species according to their shade tolerance and growth strategies (Zhao et al. 2009) (Appendix Table 4), and calculate the stability of productivity for each group.

2.4 Species diversity and stand structural attributes

In each subplot, we measured species richness (number of tree species). Two stand structural attributes, the remaining standing biomass after harvest and the quadratic mean diameter (cm), were quantified as indicators of neighborhood competition and average tree size in a stand, respectively. We did not include soil variables as explanatory factors because the four sites were located within a relatively small region with similar topography and thus not significantly differed in site condition (Appendix Table 2). Therefore, here we better focused on the biotic factors as potential drivers of ecosystem stability. To explore the growth partitioning in relation to size structure and its relationship with stand productivity, the growth dominance coefficient (GDC) was also calculated. This approach has been proposed as an insightful examination of asymmetric competition and is frequently used to predict future development of size distribution in managed forests (Forrester 2019). In each subplot, GDC were calculated according to West (2014):

$$GDC=1-\sum_{i=1}^{n}({s}_{i}-{s}_{i-1})({\Delta }_{i}-{\Delta }_{i-1})$$

where n is the number of stems in a subplot, s_i is cumulative percent biomass of tree i, and Δ_i is cumulative percent increment in biomass of tree i. Values of GDC range between − 1 and 1. A positive GDC value indicates that the growth of larger trees is proportionally greater than that of smaller trees, and a negative value indicates that the relative contribution of smaller trees to total stand growth exceeds that of larger trees (Binkley 2004). Species richness, stand structural attributes, and GDC values were calculated separately for each inventory period, and then averaged to get the mean values for each subplot.

2.5 Species asynchrony

Species asynchrony was calculated for each subplot in order to understand how the aboveground biomass productivity of multiple species fluctuates differentially in time. The degree of species asynchrony was quantified as:

$$1-\varphi =1 - {\sigma }^{2}/ {\left(\sum_{i=1}^{S}{\sigma }_{i}\right)}^{2}$$

where $\varphi$ is species asynchrony, ${\sigma }^{2}$ is the temporal variance of community-level productivity, and ${\sigma }_{i}$ is the standard deviation of productivity of species i in the subplot with S species (Loreau and de Mazancourt 2008). This index attains zero when species fluctuations are perfectly synchronized and attains one if species biomass production are perfectly asynchronized.

2.6 Statistical analysis

To deal with the spatial autocorrelation in contiguous subplots, we fitted generalized least-square models (GLS) with and without spherical autocorrelation structure for each variable predicting stability, and selected the model with lower Akaike Information Criterion (AIC) (Chisholm et al. 2013; Dormann et al. 2007). Consistent with our previous studies conducted in these forests (Geng et al. 2021; Hao et al. 2020), non-spatial models showed lower AIC values compared with pairwise spatial models (Appendix Table 5). Therefore, we did subsequent analysis based on the original variables.

One-way analysis of variance (ANOVA) was used to compare the difference in aboveground biomass productivity and its temporal stability, species richness, species asynchrony, growth dominance coefficient, and the standing biomass after harvest among harvesting treatments. Linear regression was used to fit temporal stability of productivity with species richness and asynchrony, as well as size-related growth dominance. In particular, in order to clarify whether stabilizing effects are the result of increased mean (μ), decreased standard deviation (σ), or both, we regressed μ and σ of productivity against each predictor separately. Furthermore, to determine whether predictor-stability relationships vary among harvesting treatments, slopes and intercepts were also separately fitted for each harvesting intensity level (HI0, HI1, HI2, and HI3). Multiple linear regressions were fitted at the subplot level by including a nested random effect of subplots within plots (Zuur et al. 2009). All possible combinations of seven predictors (species richness, asynchrony, GDC, the standing biomass after harvest, quadratic mean diameter, aboveground biomass productivity, and harvesting intensity level) were evaluated in an all-subsets model selection procedure. Models were ranked according to their corrected Akaike information criterion (AIC_c) values and were supposed to be equally important if the difference in AIC_c was less than two units (Burnham and Anderson 2002). To avoid the issue of model uncertainty, parameter estimates were averaged across all candidate models (ΔAIC_c < 2) (Burnham and Anderson 2002). Multiple linear regressions were performed using the “lme4” package; model selection and averaging were performed using the “MuMIn” package (Barton 2018). Standardized coefficients of multiple regression were used to compare the effect magnitude of each driver on the stability of productivity at different harvesting levels. Furthermore, to clarify the direct and indirect effects of multiple driving factors on productivity stability as well as the interplay between the drivers, structural equation model (SEM) was used. We formulated a hypothetical causal model based on prior knowledge and the results from linear mixed model. Selective harvesting may directly influence stability (direct path) or indirectly via species richness, asynchrony, and growth dominance (indirect paths). The goodness-of-fit of the SEM was tested using chi-square test (χ²), the root square mean errors of approximation (RMSEA), the Akaike Information Criterion (AIC), the comparative fit index (CFI), and P values. SEM was performed using the package “lavaan” in R (Rosseel 2012). All statistical analyses were performed with R 4.0.5 (R Core Team 2021).

3 Results

3.1 Species richness and asynchrony, stand growth dominance, and productivity stability after selective harvesting

Aboveground biomass productivity was promoted by selective harvesting of all intensities (Fig. 2a). In contrast, the temporal stability of productivity decreased with harvesting intensity from low (HI1) to high (HI3). Moderately and heavily harvested subplots (HI2 and HI3) showed considerably lower stability compared with the undisturbed subplots (Fig. 2b). Harvesting had a negligible effect on species richness over the experiment period (Fig. 2c), suggesting that harvesting treatments did not introduce bias towards removing particular species. By contrast, the treatments significantly reduced species asynchrony and growth dominance coefficient (GDC), both of which decreased with increasing harvesting intensity from HI0 to HI3 (Fig. 2d, e). GDC was close to zero in HI0, whereas strong negative growth dominance was observed in HI2 and HI3, where smaller trees account for an increasingly disproportionate amount of the stand biomass production (Fig. 2e).

Fig. 2

Boxplots of aboveground biomass productivity (a) and its temporal stability (b), species richness (c), species asynchrony (d), growth dominance coefficient (e), and the standing biomass after harvest (f) across four harvesting intensity levels. Different letters indicate significant differences at P < 0.05. HI0, harvesting intensity (basal area removed) = 0; HI1, 0 < harvesting intensity ≤ 20%; HI2, 20 < harvesting intensity ≤ 40%; HI3, harvesting intensity > 60%

3.2 Drivers of productivity stability following selective harvesting

Species richness had a positive effect on the stability of productivity for the entire community (all HI levels pooled together) (Fig. 3a). This stabilizing effect of species richness was the result of increased μ (P < 0.001), suggesting an overyielding effect of tree diversity on productivity (Fig. 3b; Table 1). Species asynchrony also had a strong positive influence on temporal stability of productivity (Fig. 3d). In contrast to species richness, asynchrony stabilized biomass production by reducing σ (P = 0.036) (Fig. 3f; Table 1). Greater stability was associated with greater values of growth dominance coefficients (GDC) demonstrating the size-asymmetric competition (Fig. 3g), causing σ to decrease strongly (P < 0.001) (Fig. 3i; Table 1). The multiple regression analysis indicated that temporal stability of productivity was mostly driven by the standing biomass after harvest (β = 0.593) and GDC (β = 0.498), followed by species asynchrony (β = 0.303), species richness (β = 0.176), and quadratic mean diameter (QMD) (β = 0.135) (Appendix Table 6). The harvesting intensity was also a good predictor of productivity stability (β = 0.297) (Table 1 and Appendix Table 6). However, aboveground biomass productivity itself was a weak predictor (β =  − 0.008) (Appendix Table 6).

Fig. 3

Relationships between temporal stability of productivity (μ/σ, gray), the mean (μ, red), and SD (σ, blue) of productivity with species richness (a–c), species asynchrony (d–f), and growth dominance coefficient (g–i). Solid lines represent significant relationships at P < 0.05, while the dashed lines represent nonsignificant relationships (P > 0.05). The colored bands represent 95% confidence intervals. SD, standard deviation

Table 1 Model outputs of linear regressions testing the hypothesized drivers of temporal stability of productivity, the mean (μ) and SD (σ) of productivity

3.3 Driving mechanisms of productivity stability under different harvest treatments

Results of ANOVA showed that the interactive terms of predictor (species richness, species asynchrony, and GDC) and harvesting intensity were all significant (Appendix Table 7), suggesting that the stabilizing effects of each predictor on stability changed across levels of harvesting intensities. When the bivariate relationships between stability of aboveground biomass productivity and each predictor were further fitted separately for each harvesting level, different patterns emerged between intensity levels. Species richness showed a positive correlation with productivity stability in HI0 and HI1 plots, whereas this richness-stability correlation became no longer significant in HI2 and HI3 plots (Fig. 4a). Similarly, significant relationships between stability and species asynchrony were found only in HI0 and HI1 treatments (Fig. 4b). On the other hand, productivity stability increased with GDC and the standing biomass after harvest for all intensity levels (Fig. 4c, d).

Fig. 4

Temporal stability of productivity as a function of species richness (a), species asynchrony (b), growth dominance coefficient (c), and the standing biomass after harvest (d) for each harvesting intensity level separately. Solid lines represent significant relationship at P < 0.05, while the dashed lines represent nonsignificant relationships (P > 0.05). The colored bands represent 95% confidence intervals. HI0, harvesting intensity (basal area removed) = 0; HI1, 0 < harvesting intensity ≤ 20%; HI2, 20 < harvesting intensity ≤ 40%; HI3, harvesting intensity > 60%

Although species richness had an overall positive influence on productivity stability, the strengths (standardized slope) of the richness-stability relationships varied with harvesting intensity (Fig. 5a). Temporal stability of aboveground biomass productivity was positively correlated with species richness in HI0 (β = 0.231, P = 0.004) and HI1 treatments (β = 0.257, P = 0.004), but had no significant relationship at HI2 (β = 0.104, P = 0.854) and HI3 levels (β =  − 0.051, P = 0.803) (Appendix Table 6). Likewise, effect magnitude of species asynchrony declined with increasing harvesting intensity, from HI0 (β = 0.577, P < 0.001), HI1 (β = 0.411, P = 0.001), to HI2 (β = 0.212, P = 0.438) and HI3 (β = 0.177, P = 0.505) (Fig. 5b, Appendix Table 6). By contrast, the stabilizing role of GDC strengthened as harvesting intensified (β = 0.368, 0.425, 0.475, 0.602 from HI0 to HI3, respectively, all P < 0.05) (Fig. 5c, Appendix Table 6).

Fig. 5

Effect size (standardized coefficients of multiple regression) of species richness (a), species asynchrony (b), and growth dominance coefficient (c) on temporal stability of productivity across the selective harvesting levels. Vertical bars indicate standard errors. GDC, growth dominance coefficient. Other abbreviations are same as Fig. 4. *P < 0.05, **P < 0.01 and ***P < 0.001

In order to determine the complex relationships between stability and explanatory variables, we further applied structural equation (SEM) modelling to find out pathways linking harvesting intensity to forest stability. SEM showed that selective harvesting intensity directly reduced stability (standardized path coefficient =  − 0.42), and indirectly affected stability via species asynchrony (standardized path coefficient = 0.25) and growth dominance coefficient (standardized path coefficient = 0.44). In contrast, productivity was promoted by harvesting (standardized path coefficient = 0.27) (Fig. 6). However, the causal pathway between productivity and stability was nonsignificant, suggesting that harvesting influenced stability mainly through changes in the extent of species asynchrony and the asymmetries between large vs. small individuals, rather than through productivity (despite a productivity surge shortly after harvest), consistent with the results of multiple linear regressions (Appendix Table 6). Notably, species richness had a positive effect on stability (standardized path coefficient = 0.22) but was not detected as an important predictor of productivity (Fig. 6).

Fig. 6

Structural equation modelling relating selective harvesting intensity, species richness, species asynchrony, growth dominance coefficient, and aboveground biomass productivity to temporal stability of productivity. Solid red and blue arrows represent significant positive and negative influences, respectively, and grey dashed arrows represent nonsignificant influences. Standardized path coefficients are given next to each path. Chi-square (χ.²) = 1.771, P = 0.413; RMSEA (root mean square error of approximation) = 0.000, P = 0.210, CFI = 0.991

3.4 Difference in the productivity stability among ecological groups

By classifying species into different ecological groups, we found that at the population level light-demanding upper-canopy (LU) and shade-tolerant (ST) groups showed higher population productivity stability, compared with pioneer (PI) and non-pioneer light-demanding lower-canopy (LL) groups (Fig. 7). In particular, both μ and σ of PI were the highest among the four groups, demonstrating a high potential growth rate but a great temporal variability. In contrast, μ and σ of ST species were remarkably lower than other groups, representing a more conservative strategy of slower growth and less fluctuation over time (Fig. 7). When further examining the stability of each ecological group separately for each harvesting intensity, results showed that none of the stability and its two components (μ and σ) of LU varied among HI levels (Appendix Fig. 9). However, LL species were more stable after harvesting treatments, due to enhanced μ (and slightly changed σ). In contrast to LL, stability of PI decreased with intensity levels, because of an increase in μ and a greater increase in σ. Harvesting destabilized ST species by reducing its μ (Appendix Fig. 9).

Fig. 7

Temporal stability of productivity and its components (mean μ and standard deviation σ) among ecological groups. LU, light-demanding upper-canopy; LL, light-demanding lower-canopy (non-pioneer); PI, pioneer and ST, shade-tolerant species. Different letters indicate significant differences at P < 0.05

4 Discussion

Disturbance legacies can have strong effects on forest ecosystem trajectories and the stability of ecosystem functioning (Stone et al. 1996). We found a decreasing temporal stability of aboveground biomass productivity with harvesting intensities. In particular, harvesting treatment with the highest intensity showed the weakest temporal stability. Furthermore, the strengths of the three stabilizing mechanisms varied across the harvesting intensity gradient: overyielding (species richness) and negative species covariance (species asynchrony) effects occurred in undisturbed and lightly harvested treatments, whereas the degree of size-asymmetric partitioning of growth contributed strongly to productivity stability when harvesting disturbance was more severe.

4.1 The stabilizing effect of overyielding

Consistent with previous studies (Aussenac et al. 2019; del Río et al. 2017; Jucker et al. 2014; Schnabel et al. 2019), we found that overall species richness promoted productivity stability in our studied forests (Figs. 3 and 5). More specifically, species richness had a significant effect on temporal mean productivity (higher μ) while keeping the standard deviation of productivity (σ) relatively constant, pointing to an overyielding effect. Nevertheless, as observed in early studies, the relationship between species richness and stability may shift along environmental gradients. Species mixing showed a weak stabilizing effect under favorable climates but became significantly positive under limiting conditions (e.g., dry, cold and infertile) (Gao et al. 2021; Garcia-Palacios et al. 2018; Lebourgeois et al. 2013; Schnabel et al. 2019; Toïgo et al. 2015), in line with the stress gradient hypothesis (Lortie and Callaway 2006; Maestre et al. 2009). Similarly, with increased light (and other resources) availability by canopy opening, the complementarity processes would be less prominent, increasing the probability of finding neutral effects of species richness on stability. Thus, the stabilizing role of richness may decrease with harvesting intensity, supporting our first prediction (Fig. 1a).

4.2 Asynchronous growth of tree species

Species asynchrony, which is the compensatory dynamics among individual species, is generally attributed to inter-specific competition (Hautier et al. 2009; Tilman 1996). However, when there is a strong environmental driver that reduces competition, compensatory dynamics might be less common (Gonzalez and Loreau 2009). Indeed, in our study selective harvesting was likely to synchronize species growth as communities are less hierarchically competitive. The magnitude of negative species covariance declined because of the concordant responses of most species to light (and other resources) released by the creation of gaps. Further, harvesting disturbance promoted the abundance of pioneer and other light-demanding subordinate species (Bazzaz 1979), especially in the heavily harvested plots (Appendix Fig. 10). Although these exploitative species responded fast to forest gaps, they exhibited lower stability compared with shade-tolerant species (Dolezal et al. 2020; Yuan et al. 2019) (Fig. 6). Indeed, the reduced stability but increased abundance of pioneer species group (PI) may contribute to an overall lower community-level stability for the high harvesting intensity levels (HI2 and HI3). For example, in BKPFs Betula platyphylla is a typical pioneer species increasing dominance after disturbance. However, they are not constant in abundance and biomass increments (especially radial growth) over time and eventually outcompeted by pines and oaks in later stages. Thus, with increasing harvesting intensity, the marked change of community composition (sometimes referred to as compositional instability, Hillebrand et al. 2018) led to reduced asynchronous species dynamics at the community level, consistent with our second hypothesis (Fig. 1b).

It should also be noted that overyielding and species asynchrony effects are not mutually exclusive stabilizing mechanisms. Several studies have found that tree species richness increased community stability indirectly via enhanced asynchrony (Morin et al. 2014; Ouyang et al. 2021; Schnabel et al. 2021). We also found a positive correlation between species richness and asynchrony (Table 1, Appendix Fig. 11). In our study site, more diverse patches are more likely to show greater species asynchrony in response to changes in local light regimes, leading to a greater temporal stability in productivity.

4.3 The growth partitioning among tree sizes

In addition to species diversity, another level of complexity in BKPFs is the existence of varying individual tree size within a stand. Tree size reflects individual tree’s functional statutes and determines the competition for resources such as light, nutrients, and water among trees (Merlin et al. 2015). It has been argued that changes in tree size structure or asymmetric competition between canopy dominants and understory subordinates were important drivers of forest productivity and its stability (Dolezal et al. 2020; Sasaki and Lauenroth 2011; Zhang et al. 2015).

As expected, selective harvesting lowered the degree of asymmetric competition, indicated by an increased contribution made by small individuals to overall stand biomass production (a decline of GDC values). However, the increased growth partitioning of smaller trees was unstable in time because they are less resistant to external stresses (e.g., fire, windthrow, pathogen, and herbivores) and threats from neighboring trees, understory shrubs, and herbaceous species compared with large trees (reviewed in Gora and Esquivel-Muelbert 2021). Relationship between stability of biomass growth of individual tree and stem DBH emphasized a common positive trend across all species and for each species separately (Appendix Fig. 12), indicating that small trees generally had lower stability than large ones. As a result, the stabilizing role of large trees may become more essential in intensively harvested plots where smaller trees account for an increasingly disproportionate amount of the growth. Appendix Table 6 also shows that quadratic mean diameter was a significant predictor for HI2 and HI3, again demonstrating that as harvesting intensity increased communities with higher mean diameter tended to be more stable. Sasaki and Lauenroth (2011) found that asymmetric competition between canopy and subcanopy, or the presence of two distinct crown strata, affected the temporal stability of a Japanese mixed forest much more than species richness.

4.4 Limitations of the present study

Our work is subject to some limitations. Species asynchrony and temporal stability were estimated based on five repeated inventories over a relatively short period of time (10 years), while long-term continuous monitoring may be required to understand the stability-management relationship. The disturbance caused by harvest is generally transitory, varying with the degree of harvest intensity. Many medium- and long-term studies have shown that the harvesting effects on biomass accumulation decrease over time as vegetation regrew, and typically became undetectable after a decade (de Avila et al. 2017; Dionisio et al. 2018; Hu et al. 2020; Schwartz et al. 2012). As forests progressively return to the pre-harvest condition, it may be difficult to detect the impacts of harvesting on these parameters. Besides, selective harvesting abruptly releases resources and triggers changes in ecosystem processes. Forest structure and composition may be shaped by the early stand development after harvest (especially tree regeneration and mortality), making the short transitional phase disproportionally important for future trajectories of forests (Seidl and Turner 2022). Existing studies also found a good agreement between independent inventory and long-term observation (Wales et al. 2020). A minimum of four repeated inventories is recommend to analyze the stability of forest productivity (Musavi et al. 2017; Yuan et al. 2019), and it may not be necessary to collect data on an annual basis due to the tradeoff between high temporal resolution and limited spatial coverage (Ouyang et al. 2021). More importantly, we showed a clear decreasing trend of asynchrony, growth dominance, and productivity stability with harvesting intensity, implying that our sample size may not reduce the power of the study. Nevertheless, longer-term studies in stability and associated factors will provide a more complete understanding of the functioning of forest ecosystems.

5 Conclusion

Our study presents one of the first in-depth analysis to identify and disentangle the relative effects of diversity-dependent (species richness and asynchrony) and size-dependent (asymmetric growth partitioning) mechanisms for stabilizing aboveground biomass productivity following selective harvesting in a diverse temperate forest. We found that compared with unharvested plots, harvesting treatments were associated with a lower temporal stability of productivity. Another important finding is that productivity stability was mostly influenced by the pattern of stand growth dominance and its changes, and to a lesser extent, by asynchronous response of different species to harvesting. Thus, beyond taxonomic diversity, stand structural factors (size, age and canopy vertical heterogeneity) that increase the growth efficiency of the whole stand needs to be considered in order to maintain forests stable. In particular, because large trees contributed disproportionately to forest biomass stock and became increasingly crucial for stability in harvested plots, assessments focusing on the growth and mortality of large trees could be invaluable for predicting forest dynamics to disturbances (Gora and Esquivel-Muelbert 2021). Furthermore, since pioneer and some less abundant light-demanding species tended to destabilize total productivity by inherently more variable and temporally more synchronous biomass production, such tree species may be identified and limited, leaving the residual forests more stable.

Appendix

Table 2 Environmental and forest characteristics of the four sites in the mixed broad-leaved Korean pine forests. Data are mean ± SE. DBH, dimeter at breast height

Table 3 Region-specific allometric equations for the calculation of aboveground biomass (kg) for each tree species. DBH, the diameters at breast height (cm)

Table 4 Relative aboveground biomass and basal area of each species, as well as their ecological group. LU, light-demanding upper-canopy; LL, (non-pioneer) light-demanding lower-canopy (non-pioneer); PI, pioneer and ST, shade-tolerant species. LL/ST indicates species that tolerate both sun and shade

Table 5 Output from the fits of the generalized least-squares (GLS) models of each stability-predictor relationship. GDC, growth dominance coefficient

Table 6 Results of multiple regression of productivity stability for all subplots pooled together (a) and for each intensity level separately (b–e). Estimates in models with ΔAICc less than 2 from the top model were averaged. Standardized regression coefficients for each predictor (β), standard error (std. error), z test (z value) and significance level (P) are given

Table 7 Summary of analysis of variance (ANOVA) of the effects of harvesting intensity, species richness (a), species asynchrony (b) growth dominance coefficient (c), and the standing biomass after harvest (d) and the interactions between harvesting intensity and each predictor on the productivity stability. df, degrees of freedom; SS, sum of squares

Fig. 8

Monthly mean air temperature (blue curves) and precipitation (grey bars) during the study period

Fig. 9

Temporal stability of productivity and its components (mean μ and standard deviation σ) among ecological groups for each harvesting intensity level (HI0, HI1, HI2 and HI3). LU, light-demanding upper-canopy; LL, light-demanding lower-canopy (non-pioneer); PI, pioneer and ST, shade-tolerant species. Different letters indicate significant differences between intensity levels at P < 0.05

Fig. 10

Proportion of the abundance each functional group under each harvesting intensity level at the end of the study period. LU, light-demanding upper-canopy; LL, light-demanding lower-canopy (non-pioneer); PI, pioneer and ST, shade-tolerant species

Fig. 11

Species asynchrony as a function of species richness

Fig. 12

Relationship between temporal stability of biomass growth and stem DBH across all species (a) and for each common species individually (b)