This study was carried out in the entire forest area of South Korea; by the end of 2020, the forests occupied about 63% of the total area of the country (Korea Forest Service, 2022). About 81% of the forest area is comprised of stands between 31 and 50 years old in forest age. The survey also showed that the forest area was made up of 37% coniferous, 32% broadleaved, and 26% mixed forests (Korea Forest Service, 2022). For the coniferous forests, Pinus densiflora Siebold and Zucc., Pinus rigida Mill., and Pinus thunbergii Parl, which all belong to the Pinaceae family, were the most abundant tree species. Quercus mongolica Fisch. Ex Ledeh., Quercus variabilis Blume, and Quercus serrata Murray, which all belong to the Fagaceae family, were the most abundant tree species for the broadleaved forests, while Q. mongolica Fisch. Ex Ledeh. and P. densiflora Siebold and Zucc. were the most abundant for the mixed forests. Under temperate climate zones, the mean annual temperature (MAT) of South Korea ranges between 7 and 15°C and the mean annual precipitation (MAP) ranges between 1,011 and 1,921 mm (Korean Meteorological Administration, 2021).

In this study, the NFI data collected from 2016 to 2020 were used to analyze the drivers that control AGB in each stand type (i.e., broadleaved, coniferous, and mixed stands) and total stands (i.e., the combination of with coniferous, broadleaved, and mixed stands). In the NFI, plots are located at the intersection of grids dividing South Korea into 4 km × 4 km areas and one plot consists of four 400-m² circular subplots (Korea Forest Service, 2017). In other words, the area of a circular subplot is 0.04 ha, and the area of a plot is 0.16 ha (Supplementary Figure 1). The 2,070 plots belonging to natural forests were finally selected and comprised 394 plots (19.0%) of coniferous forests, 829 plots (40.0%) of broadleaved forests, and 847 plots (40.9%) of mixed forests (Figure 2). In the vegetation survey, the tree species name and DBH were measured for all individuals over 6 cm DBH within the four subplots. Overall, the mean age of the forests ranged from 4.2 to 154.25 years (Supplementary Table 1). As proposed in previous studies, species-specific allometric equations based on DBH were used to calculate the AGB of each plot by the summation of AGBs of all individuals (Supplementary Table 2).

Species, phylogenetic, and functional trait diversity; CWM trait values; and stand structural attributes (Fotis et al., 2018; Ali et al., 2019b; Chun et al., 2020) were used as biotic drivers. For species diversity, species richness (SR), Shannon H index (SD), and Pielou species evenness (SE) were used. Phylogenetic trees for each of the three stand types and for the total stands were formed to quantify phylogenetic diversity using the investigated tree species and a megatree approach in R’s “V.PhyloMaker” package, as suggested in earlier studies (Jin and Qian, 2019; Chun et al., 2020). In all, we quantified phylogenetic species diversity (PD), phylogenetic SR, phylogenetic species clustering, phylogenetic species variability, phylogenetic species evenness (PSE), standard effective size of phylogenetic diversity (SES PD), standardized effect size of mean pairwise phylogenetic distance (SES MPD), and standardized effect size of mean nearest taxon distances (SES MNTD; Helmus et al., 2007; Chun et al., 2020) for each stand type.

Seven functional traits were selected to calculate functional trait diversity and CWM trait values, including leaf dry matter content (LDMC, g), specific leaf area (SLA, mm² g^–1), wood density (WD, g cm^–3), leaf nitrogen content (N, mg g^–1), leaf phosphorus content (P, mg g^–1), seed mass (mg), and maximum tree height (H, m), which are recognized as being important traits for the survival and growth of the plant and directly correlated with the functioning of the ecosystem such as AGB and stand dynamics as demonstrated in previous studies (Poorter et al., 2017; Yuan et al., 2018; Ali et al., 2019a; Chun et al., 2020). We obtained the functional traits from the published literature and open data sources or obtained the data by directly analyzing the leaf and wood samples (i.e., leaves for LDMC, SLA, N, P, and wood pieces for WD) of tree species at Dr. Chang-Bae Lee’s Laboratory following standardized protocols (Perez-Harguindeguy et al., 2013). We adapted the mean value at the same genus level for the species when species level data were absent for a species.

Functional dispersion, functional richness (FRic), functional divergence (FDiv), and functional evenness as a proxy for functional trait diversity were quantified using the R “FD” package (Laliberté et al., 2014) in R version 4.2.3 (R Development Core Team, 2023) for seven functional traits. To quantify CWM values, each trait value was calculated by weighing plot mean trait values by the relative basal area for each species (Chun et al., 2020). To describe the structural characteristics of the stand, we used the mean DBH (DBH mean), standard deviation of DBHs (DBH STD), coefficient of variation of DBHs (DBH CV), Shannon H index of DBH classes with 2-cm-diameter interval (DBH SD), skewness of DBH (DBH Skew), and stem density (Stand density) in each plot.

For each stand type and total stands, three climatic and five topographic drivers and stand age were selected to quantify the effect of abiotic drivers on AGB. The climatic drivers including MAT, MAP, and climate moisture index (CMI) were quantified in each plot from WorldClim,¹ and Global Aridity Index and Potential Evapotranspiration Climate Database.² In our study plots, MAT ranges from 5.6 to 16°C, MAP ranges from 618 to 3,870 mm, and CMI ranges from 0.1 to 0.68. Topographic drivers such as elevation, slope, topographic position index, topographic wetness index, and topographic ruggedness index were constructed using a digital elevation model with ArcGIS Pro version 3.1.0 spatial analyst tools. Annual rings were collected using cores from the five dominant trees in each subplot and then averaged across each plot as the stand age for that plot. The comparisons of the biotic and abiotic drivers and AGB along the three stand types and total stands are shown in Supplementary Figures 2–6.

To improve normality and linearity, the abiotic and biotic variables and AGB were log or square root transformed before implementing statistical procedures. The variables were standardized to unify the different units (Chun et al., 2020). For all variables, Pearson correlation analysis was conducted to remove highly correlated (| r| ≥ 0.7) explanatory variables (Supplementary Figures 7–10). In addition, to check for multicollinearity, we calculated the variance inflation factor (VIF) of each variable in the multiple regression model (Graham, 2003). Typically, the explanatory power of the model is reduced by a VIF value above 10, but in this study, multicollinearity had no effect on the results of the model as the VIF values were <3 in all models.

To compare the spatial and non-spatial model fitness, we calculated generalized least squares models (Legendre and Legendre, 2012). Spatial models include the latitudinal and longitudinal coordinates, whereas non-spatial models do not contain the coordinates for each plot. By comparing the suitability of the spatial and non-spatial models within the Akaike Information Criterion (AIC), the effects of spatial autocorrelation were found in several variables (Supplementary Table 3). Hence, we applied latitudinal and longitudinal coordinates for further multimodel inference tests and piecewise structural equation modeling (pSEM; Lefcheck, 2016) to explain spatial autocorrelation.

To select the most significant variable for AGB in each stand type and total stands, we employed multimodel inference tests as proposed in previous studies (Chun et al., 2020). As a result, the variables within the highest standardized regression coefficients (β) for AGB in each category of abiotic and biotic drivers were selected (Supplementary Figure 11).

Based on the basic conceptual model (Figure 1), we conducted pSEM (Lefcheck, 2016) to assess the multiple pathways between abiotic and biotic drivers of AGB. In this study, we implemented four pSEMs including three major stands and one total stand: (1) coniferous stand model, (2) broadleaved stand model, (3) mixed stand model, and (4) total stand model. For each model, we first constructed a model that included all pathways and then we removed the insignificant ones (Yi et al., 2021). Fisher’s C statistic, P-value, and AIC were calculated for all models to evaluate fitness.

In three major and total stands, we found that a higher variation in tree size with increasing stand age led to a higher AGB. Based on niche complementarity effects, it appears that tree size variation acts as a primary regulatory mechanism for promoting positive effects on the AGB (Zhang and Chen, 2015; Ullah et al., 2021). For instance, variations in tree size occurring within and between different species facilitate higher packing densities of diverse canopy heights, promoting the capture of above ground light, and resulting in improved light use efficiency (Zhang and Chen, 2015). In forest ecosystems, tree stands pass through successional stages over time, with the early stages characterized by a higher density of smaller, faster growing trees, and the later stages consisting of larger, slower growing trees (Lienard et al., 2015). As the stand progresses through these stages, the variation in tree size increases, leading to a wider range of tree sizes within niche complementarity and consequently an increase in AGB (Zhang and Chen, 2015; Matsuo et al., 2021). In temperate forests, most successional series require a minimum of 100 years to progress to a climax stage (Choung et al., 2020). However, most forest stands in South Korea are less than 50 years old (Korea Forest Service, 2022) due to historical forest degradation events such as the Japanese colonization of the Korean Peninsula and the Korean War, and large areas have been reforested in a short time with the three main forest types: coniferous, broadleaved, and mixed stands (Korea Forest Research Institute, 2016). In our study, across all forest types, the increase in tree size variation in the early stages of forest succession is the main factor that drives the increase in AGB, but over time, coniferous stands will decline and be replaced by the next stand type (i.e., Quercus broadleaved forests; Choung et al., 2020), so it is crucial to consider not only stand structural diversity but also biotic and abiotic drivers for increasing AGB in forest types at each forest succession stage (Qianwen et al., 2022).

In coniferous stands, we found that as stand age increased, SR decreased and CWM LDMC increased. High leaf LDMC means few spaces between the cells and resistance of the mesophyll to the diffusing gases (Bussotti and Pollastrini, 2015), which conserves resources by increasing leaf dry matter at higher construction cost to reduce water loss in coniferous stands (Qin and Shangguan, 2019). This suggests that the dominant species within specific functional trait values play an influential role in species-poor coniferous stands, in accordance with the selection effect (Poorter et al., 2015; Ullah et al., 2021). CMI, as a climate driver, had a positive indirect effect on AGB through DBH STD. Since variability in potential evapotranspiration and precipitation indicates the prominent role of soil moisture as a driver of tree growth and site productivity, the availability of moisture from higher levels of precipitation could promote greater tree growth and subsequently result in greater tree size heterogeneity in coniferous stands (Pompa-García et al., 2021). We also found that elevation increases AGB in coniferous stands. Forest carbon uptake is believed to increase along an elevation–climate environmental gradient, where greater resource availability (water and soil nutrients) occurs in temperate forests (Chun et al., 2020; Yuan et al., 2021). This is consistent with our findings of a high correlation between CMI and elevation in coniferous stands (Supplementary Figure 7). When considering the negative relationship between elevation and CWM LDMC, earlier studies have revealed that with increasing drought stress, increasing values of LDMC were observed in different tree species along an ecological gradient (Bussotti and Pollastrini, 2015). Thus, elevation, which is highly correlated with precipitation, would have had a negative direct effect on CWM LDMC in coniferous stands. For these reasons, biotic drivers and AGB are regulated by abiotic conditions, among which tree size variation and certain functional traits increase AGB through complex mechanisms in temperate coniferous stands.

Our results showed that SES MPD had a negative direct effect on AGB in temperate broadleaved stands. Generally, low SES MPD indicates that the communities are composed of species that are more phylogenetically clustered and have relatively recent common ancestors (Edwards et al., 2017). Due to the evolutionary conservatism of many species traits within a lineage, the presence of a phylogenetically clustered community structure is considered to be the most important evidence for the verification of the phylogenetic effect of environmental filtering (Asefa et al., 2017). Our results showed that MAT had a negative effect on SES MPD. Based on previous studies that have shown a strong correlation between high temperatures and growth rates and leaf nitrogen productivity in temperate broadleaved forests (Tang et al., 2018; Chen et al., 2023), these climatic conditions have created an environment conducive to specific phylogenetic groups (Ali et al., 2020), resulting in increased AGB in broadleaved stands where phylogenetic clustering occurred.

We also found that CWM P had a direct negative effect on AGB. Generally, leaf phosphorus is essential for tree metabolism and influences the distribution of photosynthate by participating in the transport of plant carbohydrates, which controls biomass accumulation (Tang et al., 2018; Yu et al., 2022). In this study, a lower CWM P was associated with a higher AGB in broadleaved forest stands. This is probably because broadleaved stands with low leaf phosphorus content strategically allocate available phosphorus resources to maximize AGB production, resulting in lower leaf phosphorus concentrations but higher total biomass, thereby maximizing their competitive advantage (Yu et al., 2022).

Our study concludes that warm climate conditions in broadleaved stands promote the formation of phylogenetically clustered communities through environmental filtering effects and that these communities employ resource allocation strategies that contribute to their high AGB in temperate forests (Liu et al., 2023).

We found that elevation had positive direct effects on AGB and CWM P, whereas it had a negative effect on PSE in the plots of mixed stands. Previous studies have shown that broadleaved tree species are known to have higher phosphorus content in their leaves than coniferous tree species (Yuan and Chen, 2009; Tang et al., 2018). Coniferous tree species mostly have a smaller SLA compared to broadleaved tree species, which limits their ability to absorb nutrients such as phosphorus (Khan et al., 2022). In contrast, broadleaved tree species tend to absorb and accumulate more phosphorus in their tissues due to their faster growth rate and higher leaf turnover, which can lead to higher AGB (Rahman and Tsukamoto, 2013). However, our results showed that forest stands dominated by broadleaved tree species with higher phosphorous content at high elevations could reduce AGB in temperate mixed stands. In natural forests, the growth of trees can be constrained by various site-specific factors including temperature, water availability, and the accessibility of available nutrients (Elser et al., 2007). In nutrient-poor environments, coniferous tree species exhibit higher growth rates compared to broadleaved tree species despite their lower biomass growth and photosynthetic rates (Sawada et al., 2016). Moreover, under less favorable conditions, such as at higher elevations, coniferous tree species can outcompete broadleaved tree species due to their higher survival rate, higher growth rate, and more efficient use of nutrients (Augusto et al., 2015). Thus, forest stands dominated by deciduous tree species with high leaf phosphorous content at higher elevations would have a low AGB in mixed stands. Moreover, we found that PSE decreased as elevation increased in mixed stands. As mentioned earlier, a higher elevation could create a more favorable environment for certain tree species to grow, environmental filtering could occur under these conditions, and phylogenetically associated species with unequal abundance could contribute to high AGB in mixed stands (Kunwar et al., 2021; Qin et al., 2022).

We also found that CMI and stand age had positive direct and indirect effects on AGB through CWM P and DBH STD in temperate mixed stands. Since elevation and CMI are positively correlated but also negatively correlated with MAT, mixed stands with enough CMI and at high elevation would have been dominated by coniferous tree species with lower leaf phosphorus content (Zhou et al., 2018). Moreover, coniferous tree species that adapt to cold and harsh environments would have better resource utilization, resulting in greater tree size variation (Ali et al., 2020). However, our pSEM results showed that CMI had a negative direct effect on AGB, which is in contrast to the results of the bivariate analysis (Figure 4C). This is expected to be due to the negative indirect effect of CMI on AGB masking the direct positive effect.

Under these conditions, as the stand age increases, mixed stands dominated by taxonomically related coniferous species that exhibit environmental filtering effects were expected to have a higher AGB in temperate mixed stands.

In total stands, we found that slope had a positive direct effect on AGB (Figure 3D). In general, low slopes tend to increase AGB in forests because gentle slopes tend to reduce soil erosion and nutrient leaching. However, about 51% of the total area of forests in South Korea are steep with slopes of 20° or more (National Geography Information Institute, 2020), and the average slope across our study plots was 28.7° (Supplementary Table 1). The forest stands distributed on the higher slopes could avoid the effects of human activities and could be better preserved, contributing to a higher plant growth rate that promotes biomass accumulation (Xu et al., 2013). We also found that slope and CMI had positive direct effects on CWM P. Generally, steep slopes are often washed by rainfall, making it difficult for vegetation to grow due to nutrient limitations and faster nutrient leaching; however, if suitable moisture is present, this can be beneficial for tree growth (Wang et al., 2020). Steep slopes induce plants to adapt by conserving and efficiently utilizing available nutrients such as phosphorus and allocating a greater proportion of these nutrients to leaves for growth and metabolic processes (Yan et al., 2016). This adaptive strategy leads to higher leaf phosphorus content as plants optimize nutrient utilization to overcome nutrient constraints on steep slopes.

In total stands, stand age had positive direct and indirect effects on AGB through DBH STD and CWM P. This is identical to the results of the coniferous forest pSEM. In line with the selection effect, this suggests that forest stands dominated by productive species with low leaf phosphorus content play an influential role in temperate forests (Poorter et al., 2015; Ullah et al., 2021). A consistent mechanism across the stand types was suggested by the observation of an increase in tree size variation with increasing stand age and suitable environmental conditions for plant growth in the three major stands in temperate forests (Koch et al., 2022).