DOI: 10.5061/dryad.79cnp5j6n
title: Tree mycorrhizal associations regulate relationships between plant and microbial communities and soil organic carbon stocks at local scales in a temperate forest

description: 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.
description: <strong>Study site and experimental
 design</strong> 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<em> et al.</em> 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 (<em>Pinus koraiensis</em>) 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 <em>P. koraiensis</em>, <em>Tilia amurensis</em>, <em>Quercus mongolica</em>, <em>Fraxinus mandshurica</em>, <em>Ulmus japonica</em>, and <em>Acer mono </em>(Figure 1a; Hao<em> et al.</em> 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<em> et al.</em> 2010). To assess local-scale variability, the 25-ha plot was divided into 625 subplots (20 × 20 m) following a standard field protocol (Condit<em> et al.</em> 1995). We calculated the slope, aspect, and convexity of each subplot using the mean elevation at each corner. (Harms<em> et al.</em> 2001). As the influence of plants on soil properties generally corresponds to plant canopy cover (Zinke 1962; Waring<em> et al.</em> 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<em> et al.</em> 2016). <strong>Soil organic carbon and physicochemical properties</strong> 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<em> et al.</em> 2004; Lindahl<em> et al.</em> 2007; Makita<em> et al.</em> 2011; Voříšková<em> et al.</em> 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<em> et al.</em> (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<em> et al.</em> 2007; Yavitt<em> et al.</em> 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<sup>3</sup> 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<sup>3</sup>). Gravel content (&gt;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<sub>2</sub>Cr<sub>2</sub>O<sub>7</sub>) in the presence of H<sub>2</sub>SO<sub>4</sub> (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<sup>-2</sup>) were calculated using the formula in (Yang<em> et al.</em> 2010):                            (1) where  is gravel content at sample location <em>i</em> (%),  is soil bulk density in the surface layer (g cm<sup>3</sup>), and  is organic carbon content (g kg<sup>−1</sup>). Finally, we interpolated the physicochemical properties for each subplot (20 × 20 m) using ordinary kriging (John<em> et al.</em> 2007). <strong>Forest stand and litterfall characteristics</strong> 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<em> et al.</em> 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<em> et al.</em> 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<em> et al.</em> 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<sup>2</sup> 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<em> et al.</em> 2013). Fine roots (&lt; 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<sup>-1</sup>) were calculated for each species. For each subplot, we then calculated community-weighted mean (CWM; Garnier<em> et al.</em> 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. <strong>Soil microbial community characteristics</strong> To characterize soil microbial communities, we randomly selected two soil sampling sites in each of the 120 subplots following Yuan<em> et al.</em> (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<em> et al.</em> 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<em> et al.</em> (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<em> et al.</em> 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<em> et al.</em> 2018). For bacteria, we annotated OTUs with functional information using the FAPROTAX database (Louca<em> et al.</em> 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. <strong>Replication statement</strong> <strong>Scale of inference</strong> <strong>Scale at which the factor of interest is applied</strong> <strong>Number of replicates at the appropriate scale</strong> Forest Stand/subplot 120 subplots (20 × 20 m) along a natural environmental gradient <strong>Statistical analyses</strong> We performed all data analyses in R 4.4.1 (R Development Core Team, 2024), using the “vegan” package (Oksanen<em> et al.</em> 2020) to calculate indices for tree species and soil microbial diversity, the “FD” package (Laliberté<em> et al.</em> 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 <em>rescale</em> 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 &lt; 10°) was log-transformed. We limited collinearity issues among explanatory variables by removing variables with a variance inflation factor &gt; 5, based on the strength of their relationship with SOC stocks. For example, we removed SRL and FDis<sub>RD</sub> due to their high collinearity with RD and FDis<sub>SRL</sub> 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 <em>dredge</em> 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 (&gt;.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.
description: # Tree mycorrhizal associations regulate relationships between plant and
 microbial communities and soil organic carbon stocks at local scales in a
 temperate forest
 [https://doi.org/10.5061/dryad.79cnp5j6n](https://doi.org/10.5061/dryad.79cnp5j6n) ## Description of the data and file structure This study was conducted in a 25-hectare forest dynamics plot within the Changbai Mountain Nature Reserve in northeastern China. Following standard field survey methods, the plot was divided into 625 subplots (20 × 20 m), and 120 of them were randomly selected for sampling. Soil samples were collected from October 1–4, 2019, at 30 m intervals. At each sampling point, 2–3 additional locations were randomly selected for supplementary sampling. Soil was drilled to a depth of 10 cm, mixed to form composite samples, resulting in a total of 967 soil samples. Soil organic carbon, pH, texture, and other physicochemical properties were measured. Simultaneously, litterfall was collected within the 120 subplots to determine carbon and nitrogen contents, and UAV LiDAR was used to obtain forest structural data. Soil microbial communities were analyzed using 16S rRNA and ITS amplicon sequencing, with gene function predictions based on the FUNGuild and FAPROTAX databases. All data were processed following standardized methods and used for subsequent statistical analyses. ### Files and variables #### File: FE-2024-00530\_for\_Dryad.zip **Description:**  The full name and unit of the variables in the data are as follows : SOC stocks (kg m-2), soil organic carbon stocks; Slope, sin(Aspect), cos(Aspect),  Convexity (m) are topographical factors; Clay-silt (%) is the combined soil content of clay and silt; EcM% is the proportion of ectomycorrhizal (EcM) tree species based on basal area; Tree SR is the species richness of the tree community; Tree BA (m-2 per 400 m-2) is the total basal area of trees; Stand ISCV is internal structural complexity of the stand; Stand Cover and Stand Height (m) are the Canopy coverage and height; Litterfall SR is the species richness of leaf litter; Litter C  (g C m-2 yr-1) is the annual amount of carbon (C) from litter deposited on the soil surface; Litter C:N is the C to nitrogen (N) ratio in annual litterfall; FDis SRL and FDis RD are the functional dissimilarity of specific root length and root diameter; SRL and RD are the community-weighted mean of specific root length and root diameter; Fungal SWI and Bacterial SWI are the Shannon–Wiener indices for soil fungi and bacteria, respectively; EcM:SF is the ratio of ectomycorrhizal to saprophytic fungal abundance in the soil; Bac. Fun. (C:N) is the first axis of PCA analysis of 11 bacterial functional attributes, representing the abundance of bacteria associated with soil C and N cycling. ## Code/software 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).