Forests are critical ecosystems to understand the global carbon budget, due to their carbon sequestration potential in both above- and belowground compartments, especially in species-rich forests. Soil carbon sequestration is strongly linked to soil microbial communities, and this link is mediated by the tree community, likely due to modifications of micro-environmental conditions (i.e., biotic conditions, soil properties, and microclimate). We studied soil carbon concentration and the soil microbial biomass of 180 local neighborhoods along a gradient of tree species richness ranging from 1 to 16 tree species per plot in a Chinese subtropical forest experiment (BEF-China). Tree productivity and different tree functional traits were measured at the neighborhood level. We tested the effects of tree productivity, functional trait identity and dissimilarity on soil carbon concentrations, and their mediation by the soil microbial biomass and micro-environmental conditions. Our analyses showed a strong positive correlation between soil microbial biomass and soil carbon concentrations. Besides, soil carbon concentration increased with tree productivity and tree root diameter while it decreased with litterfall C:N content. Moreover, tree productivity and tree functional traits (e.g. root fungal association and litterfall C:N ratio) modulated micro-environmental conditions with substantial consequences for soil microbial biomass. We also showed that soil history and topography should be considered in future experiments and tree plantations, as soil carbon concentrations were higher where historical (i.e., at the beginning of the experiment) carbon concentrations were high, themselves being strongly affected by the topography. Altogether, these results imply that the quantification of the different soil carbon pools is critical for understanding microbial community–soil carbon stock relationships and their dependence on tree diversity and micro-environmental conditions.

Study site

The study site is located in south-east China near the town of Xingangshan (Jiangxi Province, 29.08–29.11° N, 117.90–117.93° E). Our experimental site is part of the BEF-China experiment (site A, Bruelheide et al. 2014), and it was planted in 2009 after a clear-cut of the previous commercial plantation. The region is characterized by a subtropical climate with warm, rainy summers and cool, dry winters with a mean temperature of 16.7°C and a mean annual rainfall of 1.821 mm (Yang et al. 2013). Soils in the region are Cambisols and Cambisol derivatives, with Regosol on ridges and crests (Geißler et al. 2012; Scholten et al. 2017). The natural vegetation consists of species-rich broad-leaved forests dominated by Cyclobalanopsis glauca, Castanopsis eyrei, Daphniphyllum oldhamii, and Lithocarpus glaber (Bruelheide et al. 2011; 2014).

Study design

We selected 180 small-scale sample locations across five species richness levels (1, 2, 4, 8, and 16 species) per plot, according to the BEF China planting design (Appendix S1). These small-scale locations are local neighborhoods of trees defined as the twelve trees directly adjacent in the planting grid with two central trees (Appendix S1: Figs. S1-S2, see Tree Species Pairs in Trogisch et al. 2021). Each local neighborhood was replicated three times in each richness level when available (see “broken stick” design, Bruelheide et al. 2014).

Plot topography

To control for the topography effect on soil carbon concentration, topography measures were calculated from a digital elevation model (DEM). The DEM was interpolated in 2015 from elevation measurements with a differential global positioning system (DGPS) using the ordinary kriging algorithm and a cell size of 5 m x 5 m. Slope, altitude, plan curvature (Curv. PL), and profile curvature (Curv. PR) were calculated from the DEM at the plot level due to the low intra-plot variability (Scholten et al. 2017).

Microclimate modeling

The daily air temperature was recorded using 35 data loggers (HOBO® Pro v2, U23-001) installed at 1 m height in the center of 35 plots across the experiment, while a meteorological station was set up in the central part of the experimental site (see Appendix S2: Fig. S1 for more details; Bruelheide et al. 2014). To cover our full experimental area, the air temperature was modeled for all of our experimental plots using the available logger data. We modeled the temperature measurements of the 35 data loggers (i.e., daily minimum, mean, and maximum temperature) as a function of the meteorological station measurements (i.e., daily temperature, rainfall, and solar radiation), plot topography (i.e., latitude, longitude, altitude, orientation, slope, plot curvature, and mean annual solar radiation), forest vertical stratification (i.e. effective number of layers index, “ENL”, see below) and plot species richness (see Appendix S2 for more details). Spatio-temporal trends for the whole experiment were estimated using Gaussian radial basis functions (functions “auto-basis”, “eval_basis” from the R package FRK, see Appendix S2: Section S1 and Wikle, Zammit-Mangion, and Cressie 2019). Our model fits explained more than 90% of the loggers' temperature measurement variability. The fitted models were used to predict daily minimum, mean, and maximum temperature for all experimental plots with a standard error from 0°C to 2°C during our sampling period (see Appendix S2 for the complete procedure).

Field sampling

Our field measurements were performed from mid-August to the end of September 2018, before the main litterfall season (from September to December; Huang et al. 2017). To avoid spatio-temporal autocorrelation, each day another sampling area was randomly chosen. To test the effect of biotic conditions on soil microbial biomass and carbon concentrations, understory plant cover in each location was estimated on a five-level factorial scale from 'no understory plant' to 'mainly covered by understory plants'. Although this is a relatively coarse measurement, while being comparable with the Londo scale (Londo 1976), it allows considering the influence of understory vegetation which was shown to be of importance (Vockenhuber et al., 2011). We encourage subsequent studies to assess the understory vegetation in a more detailed way.

Starting from the center of the two central trees of the local neighborhood, we extracted two soil cores with 5 cm diameter and 10 cm depth, 5 cm away from the center (Appendix S1: Fig. S2). Two additional cores of the same dimensions were taken 20 cm away from the center in the direction of each tree. A composite soil sample was built for soil analyses from these four soil cores and sieved with a 2 mm mesh. As a first measure of the biotic environment, root fragments contained in the sieving residues were air-dried at 40°C for two days and weighed (± 0.01 g), while the composite soil samples were stored at -20°C.

To estimate the effect of biotic conditions and especially nutrient availability effect on soil microbial biomass and soil carbon concentration, the litter cover on the ground between the two central trees of each location was estimated on a five-level factorial scale from 'no-litter' to 'litter layer thicker than five centimeters'. Leaf litter was collected between the central trees from the ground excluding green understory plant residuals, air-dried at 40°C for two days, and milled to powder. Carbon and nitrogen concentrations were measured by micro-combustion from a subsample of 4 mg (Elementar Vario El III analyzer, Elementar, Hanau, Germany).

Soil analyses

Soil moisture was measured from a subset of 25 g soil by drying the soil at 40 °C for two days. A subsample was used to quantify soil pH in a 1:2.5 soil-water solution. Soil total nitrogen (TN) was determined on an auto-analyzer (SEAL Analytical GmbH, Norderstedt, Germany) using the Kjeldahl method (Bradstreet 1954). Soil total phosphorus (TP) was measured after wet digestion with H2SO4 and HClO4 (Parkinson and Allen 1975) using a UV-VIS spectrophotometer (UV2700, SHIMADZU, Japan). Soil total organic carbon (TOC) was measured by a TOC Analyzer (Liqui TOC II; Elementar Analysensysteme GmbH, Hanau, Germany). TOC in 2010 was quantified in a previous study (Scholten et al. 2017) at the plot level using the micro-combustion method (Elementar Vario El III analyzer, Elementar, Hanau, Germany).

Soil microbial biomass

Soil microbial biomass was measured using phospholipid fatty acid (PLFA) analysis. PLFAs were extracted from 5 g of frozen soil following Frostegård, Tunlid, and Bååth (1991). Biomarkers were assigned to microbial functional groups according to Ruess and Chamberlain (2010) using markers to assign bacteria (gram-positive bacteria: i15:0, a15:0, i16:0, i17:0; gram-negative bacteria: cy17:0, cy19:0; general bacteria markers: 16:1ω5; 16:1ω7), arbuscular mycorrhizal fungi (20:1ω9), and saprophytic and ectomycorrhizal fungi (18:1ω9 and 18:2ω6,9, see Appendix S3).

Tree biomass

Tree biomass was predicted for all neighbor trees using tree basal area (BA) and species-specific allometric relationships estimated on the two central trees. (1) Circumference at breast height (CBH) was measured in September 2018 for all trees in order to calculate the basal area of these trees. Tree height was measured for the two central trees using a laser meter (4.9 ± 2.1 m, PLR 50C, Robert Bosch GmbH, Gerlingen, Germany), and tree biomass was calculated following Huang et al. (2018). BA and biomass of the central trees were used to estimate species-specific allometric BA-biomass relationships (see Appendix S4). (3) These species-specific allometric relationships were used to calculate the neighborhood biomass (i.e., sum of the twelve surrounding trees’ biomass).

Leaf traits

For each tree species of the experiment, ten samples consisting of 10 to 25 pooled fresh leaves were collected across all diversity levels from mid-August to October 2018 (Davrinche and Haider 2021). Each sample was dried at 80°C for two days and milled for 5 min at 26 shakes per second. Carbon and nitrogen concentrations were measured by micro-combustion from a subsample of 5 mg (Elementar Vario El III analyzer, Elementar, Hanau, Germany).

Root traits

Root functional traits were measured from BEF-China Site A from September to October 2013 using two to three tree individuals per species per diversity level. First-order roots were collected, cleaned, scanned, and analyzed by WinRHIZO (Regent Software, Canada). After measurements, roots were air-dried at 60°C for two days and weighed. Average RD (in mm) and SRL (in m.g-1) were calculated from the measurements of each species at all species richness levels (Bu et al. 2017). The mycorrhizal status of the tree species was determined from the literature (Haug et al. 1994; Hawley and Dames 2004; Wang and Qiu 2006). The mycorrhizal status was assumed from the literature and confirmed by recent measurements in the same experiment (Singavarapu et al. 2021). However, intraspecific root functional trait variability can be high and may change over the course of an experiment depending on the biotic context (e.g., Zuppinger-Dingley et al. 2014), which could not be considered in the present study.

Root functional trait variables

We considered three functional root traits that are related to soil processes (Bardgett, Mommer, and de Vries 2014): root diameter (RD), specific root length (SRL), and mycorrhizal tree association (i.e., AM or EM). For each location, trait identity and diversity were calculated at the neighborhood level. At the neighborhood level, we calculated community-weighted means (CWM, Garnier et al. 2004) and functional dispersion (FDis) – defined as the weighted variance of the trait values within the neighborhood (Laliberté and Legendre 2010). All measures were weighted using tree BA. Calculations were made using the 'dbFD' function from the 'FD' package in R (Laliberté et al. 2014).

Forest vertical stratification

To quantify the forest vertical stratification and estimate crown complementarity, we computed the Effective Number of Layers (ENL, Ehbrecht et al. 2016) using terrestrial laser scanning measurements. A high ENL value indicates more evenly distributed layers, which can be an indication of higher crown complementarity and, thus, increase of canopy packing (Ehbrecht et al. 2016). A terrestrial laser scanning campaign took place in February-March of 2019 using a FARO Focus S120 and a FARO Focus X130 laser scanner (FARO Europe, Korntal-Münchingen, Germany; for more information see Perles–Garcia et al. 2021). The scanner was set up on a tripod at 1.3 m height in the center of each plot and a fully three-dimensional point cloud (360° x 305° field of view) with a spatial resolution of 6 mm at a distance of 10 m was acquired.

For each plot, the Effective Number of Layers (ENL) was computed. First, the scans were filtered using a statistical outlier removal filter (SOR, N=10, SD=3) in CloudCompare 2.9.1 software. Taking into account the dimensions of each plot (~667 m²), each point cloud was clipped in a 20 m square around the scan center (~400 m²). The point clouds were voxelized into a voxel grid of 5 cm voxels using the R package VoxR (Lecigne, Delagrange, and Messier 2018). Then, they were grouped in vertical slices of 50 cm and, for each slice, we quantified the proportion of filled voxels. The ENL was the result of calculating the inverse Simpson-Index, where n refers to the number of slices, calculated as (heightmax – heightmin ) / 50cm; and pi is the proportion of filled voxels of the ith slice.

Leaf litterfall measurement

From September to December 2018, the freshly fallen leaf litter between the two central trees of each location was collected in a 1 m2 litter trap (1 cm mesh). The collected litter was identified to species level, air-dried at 40°C for two days, and weighed (± 0.01 g). Annual amounts of litter carbon (i.e. "Clitterfall") and nitrogen (i.e. "Nlitterfall") deposited on the ground were calculated using species-specific leaf carbon and nitrogen contents and species-specific litter mass collected in the traps. We calculated the litterfall carbon to nitrogen ratio (CNlitterfall) from these measurements.