Higher tree diversity increases soil microbial resistance to drought
Forest sites and soil sampling
The sites used in this study are part of a permanent network of existing mature forest plots across Europe established in 2011–2012 (see Baeten et al.54 for detailed descriptions). We included four sites ranging over a large climatic gradient: North Karelia (Finland), Białowieża (Poland), Râşca (Romania), and Colline Metallifere (Italy), which correspond to typical boreal forests, hemiboreal mixed broadleaved-coniferous, montane mixed beech, and Mediterranean thermophilous, respectively (Supplementary Table 1, Supplementary Fig. 1). At each site, we selected 30 m × 30 m forest plots dominated by either one tree species (monospecific stands) or by three co-dominating tree species, hereafter referred to as mixed stands, resulting in a total of 34 species combinations (species were considered co-dominant if they composed >15% of the stand; see Supplementary Data file 1 for plot and tree species information). Each site differed in total species numbers, species identity, and species combinations (Supplementary Table 1). There were two replicates per tree species for the monospecific plots of each site, except for Picea abies and Quercus robur, which were only replicated once and Betula pendula which had no mono-specific plot in Białowieża. There was a minimum of three mixed species plot replicates per site that were composed of any of the target species present at the site (Supplementary Table 1), i.e., the replicate mixed plots at each site did not necessarily have the same tree species combinations. There were 64 plots in total. The sampling design with the total plot number, their distribution over four forest ecosystems, and including a wide range of tree species is well suited to address the generality of our hypothesis that microbial responses to DRW cycles are modified by tree species mixing but poorly suited to identify site-specific patterns with plot numbers too limiting within specific sites for robust testing.
Within each plot, we selected five tree triplets, a triplet being a triangle of three tree individuals within a maximum distance of 8 m from each other and no obstructing tree individuals within the triangle. Each triplet was composed of either the same species in the monospecific stands (monospecific triplet) or the three tree species present in the mixed stands (mixed triplet). At the estimated tree individual size weighted (based on individual diameter at breast height) center within the triangle, we collected five soil cores from the topsoil (10 cm deep, 5.3 cm diameter) after the litter layer had been removed. The five soil cores were spaced at roughly 35 cm from each other circling the center point (approximate sampled area 50 cm × 50 cm). A depth of 10 cm was selected because it is the standard topsoil sampling depth in soil ecology, has the highest soil microbial activity, and is under the most influence from the plant community19. All soil cores from each sampling location (i.e., tree triplet) within a plot were then sieved together through a 2 mm sieve and air-dried immediately after sampling for transportation and experiment preparation.
Experimental design
The soils collected from the 64 forest plots at the four sites were split into six replicate microcosms, yielding a total of 384 microcosms that were housed at the Montpellier European Ecotron CNRS in Montpellier, France. Each microcosm contained 95 g dry weight of soil in a glass vial (soil volume 51–72 ml; air volume 259–279 ml), initially incubated at 80% of water holding capacity (WHC) using deionized water, 25 °C, no light, and 40% relative air humidity (the vials were covered with Parafilm® to allow gas exchange but to prevent soil desiccation) for 3 weeks to reactivate the microbial community (Supplementary Fig. 3). After this acclimation period, half of the microcosms (192, i.e., n = 3 per plot) was assigned to a drying-rewetting (DRW) treatment and the other half (192, i.e., n = 3 per plot) to a control treatment. Maximum microbial mineralization activity appears to be reached between 60% and 80% WHC55. We chose 80% to ensure soils were entirely and homogeneously humid; very sandy soils with a low WHC, such as those from the Polish site, were not completely wetted at the typically chosen 60% WHC. Each treatment replicate was housed in a 2 m3 individual growth chamber (n = 6). Within each chamber, the microcosms were randomly distributed on a single shelf and re-randomized weekly. The DRW treatment was defined as two DRW cycles while the soils in the control treatment were maintained at 80% WHC throughout the experiment (Supplementary Fig. 3). Water content was adjusted gravimetrically 2–3 times a week.
Due to the large latitudinal distribution and varying soil and climate conditions of the sites (Supplementary Table 1), the soil microbial communities do not necessarily have the same degree of drought history and adaptation56. We therefore applied a site-specific drought treatment representative of each of the four study sites, i.e., site-specific drought intensity and duration. We used the permanent wilting point as a water stress threshold indicator since there is not a known microbial equivalent. The permanent wilting point was measured using a pressure plate extractor (1500F2, Soilmoisture Equipment Corp., Santa Barbara, USA) at pF 4.2 (15.5 bar) for the plots with the fastest and slowest drying soils of each site. The soil drying speed, i.e., the number of days it took for the soil to dry from 80% WHC down to constant weight, was measured gravimetrically for each plot using a subsample of soil that was subsequently excluded from the experiment. We then averaged the permanent wilting point values per site and designated this average as the drought intensity: Colline Metallifere 11% H2O g−1 dry soil, Râşca 30%, Białowieża 12%, and North Karelia 12%. The beginning of the drought was considered the moment the soil water content arrived at this threshold. The drought duration was calculated using the forest drought history data from Grossiord et al.56 as the average annual number of days the relative extractable water (REW) dropped below 0.2 (unitless) over the 1997–2010 period. REW is the ratio of available soil water to maximum extractable water (i.e., WHC), ranging between the field capacity (REW 1.0) and the permanent wilting point (REW 0.0)56. Plants are in non-water limited conditions when REW is >0.4 and water limited when REW is 2 mm diameter) and fine roots (≤2 mm in diameter). Fine roots were further separated into tree and understory roots. Tree fine roots were further divided into dead (which are hollow, brittle, and dark-colored) and live fine roots, which were then sorted by species (based on distinct color, architecture, morphology, and mycorrhizal types) and subsequently further divided based on their functions into absorptive and transport roots66. Ectomycorrhizal root tips were counted on absorptive tree roots using a binocular. All absorptive tree fine-root samples were scanned with a flat-bed scanner (resolution of 800 dpi) and scans were then analyzed using WinRhizo (Regent Instruments, Quebec, Canada, 2009) to quantify root length, surface area, volume, and diameter. Coarse root samples were also scanned to obtain coarse root volume, which was used together with the stone mass to calculate fine-earth volume (cm−3) of each soil sample. Root samples were dried (minimum 72 h, 40 °C) and weighed. Carbon and nitrogen concentrations of milled absorptive fine-root samples were measured for samples pooled at the plot level using dry combustion (Elementar Vario El Cube). Absorptive root analysis results are provided in Supplementary Table 2. We chose absorptive root traits rather than the commonly used leaf traits to characterize functional trait characteristics of tree communities, because the majority of soil microorganisms are intimately associated to the rhizosphere and thus root traits67. CWM is a measure of the relative species abundance weighted trait values. FDis is a measure of the abundance weighted mean distance between the “trait space’ of individual species. Both indices were calculated with the same standard root chemical and morphological traits (Supplementary Table 2) using the R function ‘dbFD’ in the FD package (version 1.0-1268). Due to difficulties in differentiating between Quercus species in root samples from some of the Italian plots, we were unable to determine mean absorptive root trait values at the plot level. We therefore used mean root trait values at the site level calculated from the mono-specific stands. Although the root trait values were not at the plot level, we were still able to determine the CWM and FDis indices at plot level since the root traits values were reported to tree species relative abundance in each plot. The relative abundance of each tree species was calculated using the basal areas of the tree individuals used in the five plot tree triplets (three tree individuals per tree triplet). Within each plot, the basal areas of a tree species (including five or fifteen tree individuals depending on whether the plot was mixed or mono-specific, respectively) were summed and then reported to the total basal area of the 15 tree individual, giving the relative basal area of each tree species within each plot. In order to synthesize this data, we incorporated them in a principal component analysis (PCA) and extracted the first axis scores (explaining 52.8% of the variance; Supplementary Fig. 2). Although the evidence supporting a universal root economics spectrum (RES) for woody species is inconsistent69,70,71, we consider our PCA1 axis as an acquisitive to conservative trait gradient with lower scores represented acquisitive root traits (high N content, specific root length, and ectomycorrhizal colonization intensity) and higher scores represented conservative traits (large diameter and high tissue density). The FDis was calculated following Laliberté & Legendre64 based on all traits at the plot level. The mono-specific stands had a FDis value of zero, which limits FDis variability for half of the plots. Accordingly, there was one single FDis value per plot that was used in our statistical analyses.
Since soil microbial resistance and recovery are tied to soil parameters and resource availability17,18, we also included major topsoil parameters (0–10 cm) known to affect microbial activity and/or community composition (Supplementary Table 2) measured previously during the FunDivEurope project54 at the plot level. Similar to the CWM absorptive root traits, we incorporated the topsoil variables into a PCA using the function ‘prcomp’ from the factoextra package (version 1.0.662) and extracted the first axis scores (explaining 52.5% of the variance; Supplementary Fig. 2) for a synthetic soil parameter measure for each individual plot. High PC1 scores are associated with higher pH, carbon content, and clay content and lower bulk density, the inverse is correlated with low PC1 scores.
We used generalized mixed-effects linear models (two-sided) using the lme4 package (version 1.1–2172) to assess the effects of the DRW treatment and the influence of the tree species number on microbial C and N-related parameters. The root FDis, root CWM PC1, and soil PC1 variables were included with the treatment × tree species interaction as explanatory variables. For the response variables (instantaneous CO2 and N2O fluxes measured five times over the experiment and cumulative fluxes, DOC, and TDN leaching, qCO2, and resistance and recovery indices), extreme values were removed (±3 times the IQR of all values for each variable). The soil collection site and plot as well as the growth chambers used for the incubation were included as random variables with plot nested within site. We did not include any climatic variables from the different sites, because they were highly correlated to site, which was already a random effect in the model. The model structure was as follows: response variable ~ Root FDis + Root CWM PC1 + Soil PCA axis + Treatment * Tree species number * Flux measurement time + (1|Chamber) + (1 | Site/Plot). The “Flux measurement time” variable, which identifies the times the five flux measurements were taken (i.e., beginning, drought 1, rewetting 1, drought 2, rewetting 2; Supplementary Fig. 3), was used only in the models that looked at the temporal dynamics of CO2 and N2O fluxes. For the analysis of resistance and recovery indices, we did not keep the “Treatment” variable in the model since these indices were calculated using both the DRW and control treatment results (see above). Additionally, for the resistance and recovery indices, instead of a “Flux measurement time” variable, a “Cycle” variable was included to distinguish the microbial activity resistance and recovery of the first and second cycles; the “Cycle” result indicates the change between the first and second cycle. Model residuals were plotted to test for normality, and data was transformed (log2 or BoxCox) when normality was not met. We also verified for data homogeneity and model probability (Q–Q plots). In order to identify the most parsimonious model we used the R software (version 3.5.3) and the “dredge” function in the MuMIn package (version 1.43.673) which uses the lowest Akaike information criteria (AIC) to rank all possible models with all possible combinations of the explanatory variables in the full model.
The data presented here is tied to specific spatial and temporal ecological conditions (e.g., forest drought history, tree species presence, microbial community composition, and soil property heterogeneity) which are susceptible to change. This makes exact study replication challenging and underlines the importance of including a wide range of conditions (e.g., multiple forest types, tree species, tree species combinations, climatic conditions, and soil types) as done here in order to explore general, potentially reproducible, trends oppose to site-specific trends.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More
