Forests store substantial amounts of soil organic carbon (SOC), but SOC stocks differ strongly between forest ecosystems dominated by arbuscular mycorrhizal (AM) or ectomycorrhizal (EcM) fungi. In temperate forests, nearly all tree species associate with either AM or EcM fungi, but it is unclear if variation in SOC stocks is linked to the dominance of AM vs. EcM trees at local scales. However, SOC stocks are also influenced by many other factors, including plant diversity, plant traits, soil properties and microbial community composition, that can vary substantially at small spatial scales. Thus, elucidating how tree mycorrhizal associations interact with other drivers of SOC formation at local scales could improve our understanding of forest SOC storage. Using multivariate data from a 25-ha temperate forest plot, we hypothesized that SOC stocks would be greater in AM-dominated than in EcM-dominated subplots, due to the rapid decomposition of labile litter inputs from AM trees. We expected that variation in SOC stocks would be explained by interactions between AM-associations, plant diversity or traits, and soil microbial communities. We also accounted for differences in soil abiotic conditions and used piecewise structural equation models to identify potential causal pathways. Forest stand attributes played a primary role in predicting forest SOC stocks (45% relative contribution), compared to microbial community composition (26%) and abiotic factors (20%). Forest SOC stocks increased with plant and litterfall diversity and root diameter but declined with soil clay-silt content. Importantly, forest SOC stocks were lower in subplots with a high ratio of EcM to saprotrophic fungi and subplots with high predominance of bacterial functions related to carbon and nitrogen cycling. Our models suggest that smaller SOC stocks in stands with high EcM dominance are indirectly linked to lower tree species diversity, greater litter carbon input and distinct soil microbial community composition. Our findings highlight that tree mycorrhizal type also influences SOC storage at local scales but stand-level relationships between SOC stocks and mycorrhizal dominance differ from regional and global patterns. We show that mycorrhizal associations regulate complex relationships between producers, decomposers, and soil properties at small spatial scales. Considering these linkages could therefore improve estimates and management of temperate forest SOC stocks.

Study site and experimental design

The research area was situated within the Changbai Nature Reserve in Northeast China (41°42 ′ to 42°26 ′ N and 127°42 ′ to 128°17 ′ E; Figure 1a). We conducted our sampling within the 25-ha Changbaishan (CBS) Forest Dynamics Plot in the core zone of Changbai Nature Reserve (Figure 1b,c), which is one of the CTFS-ForestGEO global forest monitoring sites (http://www.forestgeo.si.edu). The climate in the study area is temperate, with average annual temperature and precipitation of 2.8°C and 700 mm, respectively (Wang et al. 2012). The elevation within the study site ranges from 791.8 to 809.5 m.a.s.l. All fieldwork was conducted with the authorization of the Jilin Changbaishan Forest Ecosystem National Observation and Research Station.

The 25-ha plot is mixed forest dominated by late-successional deciduous broadleaved Korean pine (Pinus koraiensis) on slightly acidic dark brown forest soil (mollisol according to States & Service 1975), which formed over granite and basalt. The most common tree species are P. koraiensis, Tilia amurensis, Quercus mongolica, Fraxinus mandshurica, Ulmus japonica, and Acer mono (Figure 1a; Hao et al. 2007). All individual trees with diameter at breast height (DBH) ≥1 cm were tagged, identified and measured in 2004, 2009, 2014 and 2019. According to the initial census data from 2004, there are 38,902 individuals belonging to 52 species, 32 genera and 18 families (Wang et al. 2010). To assess local-scale variability, the 25-ha plot was divided into 625 subplots (20 × 20 m) following a standard field protocol (Condit et al. 1995). We calculated the slope, aspect, and convexity of each subplot using the mean elevation at each corner. (Harms et al. 2001). As the influence of plants on soil properties generally corresponds to plant canopy cover (Zinke 1962; Waring et al. 2015), we randomly selected 120 subplots for sampling that were at least 20 m distant from the an adjacent sampling plot or the plot edge (Figure 1d). Each subplot contained 21-91 individual trees, allowing us to assess linkages between tree diversity, soil microbial communities, and soil properties at a relevant spatial scale (Kershaw Jr et al. 2016).

Soil organic carbon and physicochemical properties

We measured organic C stocks in the surface horizon of the mineral soil (0-10 cm), as this horizon contains the highest proportion of fine roots and fungal hyphae (Fisk et al. 2004; Lindahl et al. 2007; Makita et al. 2011; Voříšková et al. 2014), and interactions between biotic and abiotic factors are most likely to influence SOC formation and storage. To estimate soil C stocks, we collected soil samples at 30-m intervals throughout the plot on 1–4 October 2019 following methods in John et al. (2007). To capture the differences in SOC stocks at smaller scales, two more sample locations (2, 5, or 15 m) were chosen at a random cardinal or intercardinal direction from each base point (Figure S2; John et al. 2007; Yavitt et al. 2009). We first removed any large debris, and we used a 3.8-cm diameter soil drill to take three subsamples to a depth of 10 cm. The three sub samples were then mixed thoroughly to create a representative composite sample. Five sampling locations were inaccessible due to the presence of large roots, tree trunk, or stumps. As a result, we obtained 967 soil samples, which were air-dried and subsequently sieved using a 2 mm mesh screen, while the roots were carefully extracted.

To determine soil bulk density, an additional sample was collected at each sampling location, using a ring core (100 cm³ core volume). The soil was dried at 105 °C for 48 hours and soil bulk density was determined by dividing the dry weight of the soil by the volume of the core (100 cm³). Gravel content (>2 mm) was then determined by sieving the soil samples (2 mm mesh).

Soil pH in water (1:1 ratio of soil to solution) was determined using a pH meter (PHB-5, Beckman, California). SOC content was determined by titration with a strong oxidizing agent (K₂Cr₂O₇) in the presence of H₂SO₄ (Nelson & Sommers 1983). Soil texture (sand, silt, and clay content) was determined by laser diffraction (Mastersizer 2000, Malvern Instruments, Malvern, England).

As we only sampled one soil horizon (0–10 cm), SOC stocks (kg m^-2) were calculated using the formula in (Yang et al. 2010):

                           (1)

where  is gravel content at sample location i (%),  is soil bulk density in the surface layer (g cm³), and  is organic carbon content (g kg⁻¹). Finally, we interpolated the physicochemical properties for each subplot (20 × 20 m) using ordinary kriging (John et al. 2007).

Forest stand and litterfall characteristics

Forest stand characteristics in each 20 × 20 m subplot were derived from the 2019 plot census data. Tree diversity and biomass were represented by species richness (Tree SR) and basal area (Tree BA). Tree mycorrhizal associations (AM or EcM) were determined based on published records (Wang & Qiu 2006; Soudzilovskaia et al. 2022), and online resources (http://mycorrhizas.info/index.html), which resulted in 39 AM tree species, and 9 EcM tree species. EcM dominance was determined as the proportion of EcM-associated basal area compared to the overall basal area, expressed as a percentage (EcM%) in each subplot. Thus, low EcM tree dominance indicates subplots with high AM-dominance.

To quantify forest structure and estimate crown complementarity, we computed stand structural complexity attributes based on unmanned aerial vehicle (UAV) lidar data in August 2022 (GreenValley International LiHawk system; (Ma et al. 2022). The system has a RIEGL VUX-1 UAV laser scanner that can detect objects up to 1000 m away, and rapidly acquires data (550 kHz) using a narrow near-infrared laser beam. The collected UAV lidar data were preprocessed using the LiDAR360 software (GreenValley International Inc., Beijing, China), which denoises, filters, and normalizes the data to remove noise from the point clouds and classify ground points; a digital terrain model (DTM) is then generated from the ground points. Three metrics of structural complexity were extracted from the normalized lidar point clouds: stand cover (horizontal) was calculated as the percentage of vegetation points to the total number of lidar points (Ma et al. 2017); stand height (vertical) was calculated as the mean surface height of the forest canopy; and the internal structural coefficient of variability (ISCV) was calculated as a structural complexity attribute using the standard deviation and mean for height in each subplot.

To test our assumptions of litter quality and quantity, we collected freshly fallen leaf litter in a 0.5 m² litter trap (0.2 cm mesh, 1 m above the ground) at the central points of the 120 subplots from September to December 2019. We identified and separated the litter to species level, air-dried the samples at 80°C for 12 hours and weighed them (±0.01 g). Then, we determined litter carbon and nitrogen concentrations for each species by combustion (Elementar Vario El III analyzer, Elementar, Hanau, Germany). We calculated the annual amounts of deposited litter C and litter N in each subplot using species-specific litter mass multiplied by their litter C and N contents, respectively, and calculated the litter C:N ratio. Subplot litter diversity was expressed as the total number of leaf litter species found in each trap.

We measured root functional traits for 5 to 20 individuals per tree species within the whole 25-ha plot (Pérez-Harguindeguy et al. 2013). Fine roots (< 2 mm in diameter) were collected, carefully cleaned, and then scanned and analyzed using WinRHIZO (Regent Software, Canada). The roots were subsequently dried to constant weight at 80°C to determine dry weight. Average root density (RD in mm) and specific root length (SRL in cm g^-1) were calculated for each species. For each subplot, we then calculated community-weighted mean (CWM; Garnier et al. 2004) root traits, root diversity, and functional dispersion (FDis), defined as the weighted variance of the trait values within the subplots (Laliberté & Legendre 2010). Community weighting was based on tree species BA in each subplot.

Soil microbial community characteristics

To characterize soil microbial communities, we randomly selected two soil sampling sites in each of the 120 subplots following Yuan et al. (2020) (Figure 1b). At each sampling site, surface litter was first removed, then five soil cores (3.8 cm diameter, 10 cm deep) were collected, pooled to make one composite sample per sampling point, and stored at −80°C for microbial analyses. Laboratory analyses were performed on two composite samples per subplot, but data were averaged across the sampling points to give mean soil microbial diversity and community composition per subplot for statistical analysis.

We extracted soil bacterial and fungal DNA from 0.25 g and 0.5 g fresh soil, respectively, using PowerSoil DNA Isolation Kits (MoBio Laboratories, Carlsbad, CA, U.S.A.). We used agarose gels to check the quality and size of the extracted DNA before storing the samples at −20 °C until further use. For bacterial 16S rRNA gene amplification by PCR, the V4~V5 region was targeted with the primer pair 515F/907R (Tamaki et al. 2011). For fungal ITS region amplification, the primers ITS1F and ITS2F were utilized. (GeneAmp 9700, ABI, USA). PCR reactions were conducted in triplicate under standard conditions following Zhang et al. (2022b). The PCR products were sequenced using 300PE MiSeq (Illumina, San Diego, CA, USA), at Shanghai Majorbio Bio-pharm Biotechnology Co., Ltd (Shanghai, China).

Subsequently, raw reads were processed using QIIME software (www.qiime.org) by demultiplexing, quality filtering, and OTU clustering with a 97% similarity threshold using UPARSE. Chimeric sequences were detected and removed using UCHIME (Edgar 2013). A total of 28,872 fungal and 17,463 bacterial effective sequences were obtained from 240 soil samples across 120 subplots. Taxonomic assignments were performed using the UNITE database for fungal sequences and the Greengenes corset database for bacterial sequences, based on a confidence threshold of 70%.

We used functional annotation to categorize mycorrhizal types and predict microbial functional groups based on OTUs. For fungi, we classified fungal sequences into three trophic groups: symbiotroph, saprotroph, and pathotroph using FUNGuild (Nguyen et al. 2016). Symbiotroph OTUs were subdivided into EcM fungi, AM fungi, and other symbiotic fungi. OTUs assigned multiple trophic strategies were excluded from further analysis. For each soil sample, we calculated the relative abundance of each fungal trophic group (Table S1). The ratio of EcM fungi to saprotrophic fungi (ECM:SF ratio) was used to represent the functional composition of the fungal community (Gallart et al. 2018). For bacteria, we annotated OTUs with functional information using the FAPROTAX database (Louca et al. 2016). We were able to assign the majority of annotated species to groups associated with ecologically significant functions, such as the C, N, and sulphur (S) biogeochemical cycles and the degradation of organic matter. From these, we selected 11 bacterial functions involved in soil C and N cycling (Table S1), and used the first axis of a principal component analysis (PCA axis 1, explaining 46% of the variation; Table S2) to represent the overall C and N cycling function of the bacterial community. The Shannon-Wiener index was used to calculate the fungal and bacterial diversity of each sample based on OTUs.

Replication statement

Statistical analyses

We performed all data analyses in R 4.4.1 (R Development Core Team, 2024), using the “vegan” package (Oksanen et al. 2020) to calculate indices for tree species and soil microbial diversity, the “FD” package (Laliberté et al. 2014) for root functional trait analyses, the “MuMIn” package (Bartoń 2020) for model averaging, and the “piecewiseSEM” package (Lefcheck 2016) for structural equation model. Before analysis, all explanatory variables were standardized by centering and dividing them by two standard deviations using the rescale function in the “arm” package (Gelman 2011).

To examine the relationships between SOC stocks and stand or microbial community attributes (tree diversity, stand structure or functional composition), and account for the effects of environment factors (topography and soil physicochemical properties), we performed multiple linear regressions. As aspect is a circular variable, we used sin (aspect) and cos(aspect) in regression models, and slope (usually < 10°) was log-transformed. We limited collinearity issues among explanatory variables by removing variables with a variance inflation factor > 5, based on the strength of their relationship with SOC stocks. For example, we removed SRL and FDis_RD due to their high collinearity with RD and FDis_SRL but weaker relationship with SOC stocks. Using this approach, we incorporated a total of six parameters related to environmental conditions, ten forest stand attributes, and four attributes of the soil microbial community, as predictors in the multiple linear regression models to estimate the SOC stocks. Next, we conducted regression analysis on each predictor using all possible subsets and chose the best model based on the Akaike Information Criterion adjusted for small sample sizes (AICc). If the DAICc was less than 2 units, we averaged the models using the dredge function in the MuMIn package.

We identified direct and indirect links between SOC stocks and plant, microbial, or environmental variables by constructing piecewise structural equation model (pSEM) based on the initial conceptual framework. Predictors were selected based on the results of the multiple linear regression and conceptual models, and subplot was fitted as a random effect. We assessed spatial autocorrelation (SAC) of the residuals for each path in the SEM by calculating Moran's I (Moran.I function in the ape package; Paradis & Schliep 2018) but spatially corrected pSEMs failed to converge. We therefore present the results of the uncorrected pSEM, while acknowledging that the presence of SAC could inflate the strength of some of the relationships between variables. The final model fit was evaluated using a Fisher's C value, non-significant p value (>.05) and Akaike information criterion (AIC). Indirect linkages between each predictor and SOC stocks were determined by multiplying the standardized direct effects of the predictor on SOC stocks via mediators for each relevant path. The total indirect effect of a predictor was then determined by summing all its indirect effects. To quantify the relative contribution of different predictors to SOC stocks, we divided the beta coefficient of each predictor by the sum of the beta coefficients of all predictors to provide a measure of relative importance.