Understanding mechanisms stabilizing ecosystem functions, such as primary production, is crucial for forecasting global environmental responses. While biological diversity is expected to enhance stability through compensatory reactions to environmental changes, empirical evidence is lacking, especially in old-growth forests vital for biodiversity conservation and climate change mitigation. Moreover, whether increased niche complementarity and stronger intraspecific than interspecific competition are key mechanisms promoting compensatory dynamics and stabilizing ecosystem functions in diverse forests remains unexplored. This study investigates productivity and stability in temperate old-growth forests over 20 years at community and individual levels. Analyzing 4,380 trees in a 4-hectare plot in northern Japan with over 35 tree species, structural equation models evaluated the effects of biodiversity and average asynchrony in species fluctuations (compensatory dynamics) on productivity stability across 100 m² grid quadrats. Functional traits and taxonomic diversity represented species complementarity and the insurance effect. Temporal growth correlations between conspecific and heterospecific neighbors and neighborhood effects on growth performance indicated intra- and interspecific interactions at the individual level. Communities with greater stability exhibited higher diversity and asynchronous species fluctuations, suggesting that compensatory dynamics buffer community productivity against environmental variability. The inverse relationship between tree size variation and stability indicates that communities with less pronounced size and abundance hierarchies have more efficient compensatory mechanisms, ensuring stable forest functioning. The absence of negative temporal correlations in biomass production among heterospecific neighbors suggests the limited significance of interspecific competition in compensatory dynamics. Conversely, positive correlations among conspecific neighbors and their suppressed growth in dense conspecific patches highlight the importance of conspecific negative density-dependent mechanisms in sustaining tree species diversity and ensuring stable productivity. The study underscores the critical role of tree species richness in stabilizing ecosystem functioning via asynchronous growth in one of the world’s most diverse temperate forests. Stronger intraspecific than interspecific competition helps prevent single-species dominance, maintaining diversity and productivity stability. Despite occasional destabilization from size-asymmetric interspecific competition, species trait complementarity enhances stability by promoting overall biomass production. This study highlights the importance of overall diversity for the stability of forest productivity, with implications for nature conservation and ecosystem functionality.

The study was conducted within a 4-hectare plot (200 × 200 meters, divided into 1,600 5 x 5 m quadrats) established in 1996 in an old-growth stand within the 2,715-hectare Tomakomai Experimental Forest. All trees exceeding a height of 1.3 meters were tagged, and their girth at breast height of 1.3 m (GBH) was measured in 1996, 1998, 2000, 2002, 2004, 2006, 2009, 2012, and 2015. To comprehend their impact on average productivity and temporal stability, we investigated various factors, including taxonomic, phylogenetic, and functional (structural) diversity, compensatory dynamics (species asynchrony), stand density, and aboveground biomass. These associations were examined at a resolution of 100 m² by dividing the 4 ha plot into 400 quadrats, each measuring 10 m by 10 m (local communities).

The average productivity of community biomass in each quadrat was determined by analyzing changes in stand basal area increments from 1996 to 2015 (the sum of individual tree BAI calculated from remeasured GBH), while total stem density and aboveground biomass in each quadrat were estimated from all trees in 1996 and their stem basal area sums. We evaluated the temporal stability of productivity by computing the inverse of the squared temporal coefficient of variation (1/CV²). This metric indicates the ratio of the mean to the standard deviation of stand basal area growth over time (Wang and Loreau 2014). We quantified taxonomic, phylogenetic, and functional characteristics to examine the influence of multiple diversity facets on biomass productivity and its temporal stability directly and indirectly via elevated species asynchrony. Taxonomic diversity was assessed through species richness (tree species per quadrat). Phylogenetic diversity was computed from the studied tree species' phylogeny reconstructed from GenBank's nucleotide sequences. Functional diversity measures were computed for each quadrat from individually recorded stem girths, leaf traits, and branching patterns. Each tree's average relative growth rate (RGR) and its variation between 1996-2015 were calculated for each species in each quadrat from stem basal area changes (RGR = ln BA₂₀₁₅ – ln BA₁₉₉₆).

Furthermore, we used the fast–slow leaf economics spectrum (LES) to proxy for differences in tree leaf adaptive strategies. Finally, productivity stability may be influenced by how species respond to canopy disturbances and gap dynamics. To explore this, species were divided into two groups: those with monopodial branching, where the central leader shoot continues to grow while lateral branches stay subordinate, and those with sympodial branching, featuring a bifurcating pattern where one branch dominates. Sympodially branching trees are expected to be more adaptable to increased light levels following gap formation than monopodial species, which generally adopt more conservative crown formation strategies.

To assess how tree traits (mean and variance) impact community productivity and stability, we computed community-weighted mean (CWM) and functional diversity (FD) indices for each trait in every grid quadrant. Connections among community productivity and its temporal stability, species asynchrony, and diversity measures were tested by fitting a piecewise structural equation model to aggregate ecosystem properties in the 400 10 m x 10 m quadrats. We estimated a variance inflation factor (VIF) to assess whether multi-collinearity affected parameter estimates and only variables having VIF < 5 were included in the model. Fisher's C evaluated the goodness-of-fit for the whole model, and Nagelkerke pseudo-R² values were reported to show the explained variation (Lefcheck, 2016).

SEMs were fitted by ‘piecewiseSEM’ package version 2.3 (Lefcheck, 2016), VIF was computed by ‘car’ package version 3.0-12, species diversity was computed by ‘vegan’ package version 2.6-4, phylogenetic diversity by package ‘picante’ version 1.8.2, functional diversity by ‘FD’ package version 1.0-12.1 all in R version 4.3.1.