Material and methods
Experimental design. This experiment was set up as a fully factorial design containing all combinations of three levels of soil biodiversity (low, moderate, high) and four levels of plant species richness (1, 4, 8, 12 species) (Fig. 1 and Supplementary Methods, Table S1). We created the four- and eight- plant species treatments, which were replicated 7 times, by randomly selecting plant species from a 12-species pool for each replicate. For each replicate of plant diversity treatments, there were identical species composition among low, moderate and high soil biodiversity treatments to avoid a confounding effect of plant community composition and soil biodiversity treatments. There were 12 treatments and a total of 123 grassland microcosms. Twelve typical species from a local grassland in Brandenburg, Germany, were used in this study: four grasses (Holcus lanatus , Anthoxanthum odoratum , Lolium perenne and Festuca rubra ), four herbs (Daucus carota ,Achillea millefolium , Hieracium pilosella andPlantago lanceolata ) and four legumes (Trifolium repens ,Vicia cracca , Medicago lupulina and Lotus corniculatus ).
Soil inoculum preparation. We used the dilution-to-extinction approach (Yan et al. 2015; Hol et al. 2015; Roger et al. 2016; Maron et al. 2018; Domeignoz-Horta et al. 2020) to create high, moderate and low soil biodiversity (Fig. 1 and Supplementary Methods). Fresh field soil was collected from the top 20 cm of the local grassland in Brandenburg (Rillig et al. 2010). The fresh field soil was treated as the undiluted treatment (100). We mixed 200 g dw (dry weight) of fresh soil with 1800 g dw sterilized soil (autoclaving for 90 min at 121°C) to obtain the 10-1 dilution treatment, and the 200 g of the 10-1 dilution was mixed with 1800 g dw sterilized soil to obtain the 10-2 dilution treatment. This procedure was repeated several times to reach 10-3 and 10-6 dilutions. A plastic bag was filled with 2000 g dw of the 100, 10-3 or 10-6 dilutions, and then received sterile water to reach the initial soil moisture. To recover soil microbial biomass, all soil dilutions were incubated in a dark room (20°C) for two months, until no differences among soil dilutions were observed using the substrate-induced respiration method (Fig. 2a). The 100, 10-3 and 10-6dilutions of soil inocula were treated as high, moderate and low soil biodiversity treatments, respectively.
Microcosms establishment and sampling. Each grassland microcosm (22.5 cm diameter and 16.5 cm height) was established by filling with 5.5 L (6.8 kg) of a sterilized 1:1 sand:field soil mixture (autoclaving for 90 min at 121°C) and 200 g of soil inoculum. Microcosm with 1, 4, 8 or 12 plant species richness treatment received 24, 6, 3 or 2 seedlings of each species, respectively. All microcosms were maintained in a climate-controlled greenhouse.
We simulated environmental variation in precipitation by inducing three wet-dry cycles (Fig. 1). Each microcosm was watered twice weekly to maintain gravimetric soil moisture of 12-18% by weight during the wet periods. During the dry periods, each microcosm received 300 mL of water only when most plant species started to wilt. Each of wet and dry periods lasted at least for 8 weeks. At the end of each period, all plant shoots were harvested by cutting at 5 cm above the soil surface, sorted by species, oven-dried for 48 h at 70°C and weighed. In total, there were six such harvests, including shoot biomass from 17,712 plant individuals. Based on the maximum nutrients removed by the first harvest, 400 mL of Hoagland nutrient solution was added to each microcosm after each harvest.
Soil fungal and bacterial diversity . After the final harvest, 100 g of fresh soil was collected from each of the 12-species microcosms in each soil biodiversity treatment. Fresh soil samples were stored at -80°C for DNA extraction to evaluate the establishment of the soil microbial community during the experimental period. DNA from each soil inocula and fresh soil samples of the final harvest was extracted from 250 mg soil, using DNeasy PowerMax Soil Kits (MoBio Laboratories Inc., Carlsbad, CA, USA), following manufacturer’s instructions. Soil fungal and bacterial diversity were determined following Illumina MiSeq high-throughput sequencing with fITS7 and ITS4 for fungi and 515f and 806r for bacteria (Fierer et al. 2005; Ihrmark et al.2012). Sequencing raw data were processed using DADA2 (Callahan et al. 2016) to obtain denoised, chimera-free, non-singleton amplicon sequence variants (ASVs). Soil microbial diversity was indicated by the sum of bacterial and fungal ASVs in each soil biodiversity treatment. Details of sequence processing can be found in the Supplementary Methods.
Temporal stability of biomass production . The temporal stability of plant community biomass production in each microcosm is defined as µ /σ , where µ is the temporal mean of biomass production during three rounds of wet–dry cycles, and σis the standard deviation of biomass production across this experimental period (Tilman et al. 2006). We partitioned the temporal stability into community-wide plant species asynchrony and weighted population variance (Thibaut & Connolly 2013). Plant species asynchrony (Loreau & de Mazancourt 2008) was calculated as 1 – φ = 1 –σ 2/(\(\sum_{i=1}^{s}\sigma\)i)2, whereφis species synchrony, σ2 the temporal variance of biomass production, σi is the standard deviation in the shoot biomass of species i in a community withS species. If 1 – φ equals 0, species fluctuate synchronously, indicating that there is no compensatory effect; if 1 – φ is higher than 0, species fluctuate asynchronously, indicating that there are compensatory effects (Loreau & de Mazancourt 2008; Song & Yu 2015). Asynchrony of plant functional groups was calculated by replacing species i by each plant functional group, which indicates the strength of asynchronous responses at the plant functional group level. Weighted population variance (CVpop ) is defined asCV pop = (\(\sum_{i=1}^{s}\sigma\)i)/µ . A decrease in population variance can promote temporal stability by reducing the variability of populations (Thibaut & Connolly 2013).
Functional composition and diversity . To better understand the role of plant diversity, we calculated four functional diversity indices for each plant community using ten leaf and root traits relevant to drought tolerance: functional richness, functional evenness, functional divergence, and functional dispersion (Villéger et al. 2008; Laliberté & Legendre 2010; Fischer et al. 2016; Liu et al. 2018; Lozano & Rillig 2020; Freschet et al. 2020; Mahautet al. 2020) (See Supplementary Methods for details). Principal component analysis (PCA) of community-weighted means of ten traits was employed, using the PCA function in FactoMineR (Lêet al. 2008). The first axis of PCA analysis was treated as an index of the community-weighted mean (CWM) of all traits (Cravenet al. 2018; Mahaut et al. 2020). Functional diversity indices of each harvest were calculated, and then were averaged across the six harvest time points.
Statistical analysis . We used linear models to test the effects of the dilution-to-extinction approach on soil microbial, fungal and bacterial diversity, and microbial biomass. Microbial data were log-transformed before analyses to meet the normality assumption. We tested whether soil biodiversity interacts with plant diversity to influence response variables using a linear model. Plant species richness was log-transformed to represent the biodiversity effect which typically saturates with increasing species richness due to increasing redundancy (Hautieret al. 2014). For each response variable, a weighted least squares regression was applied to take heteroskedasticity into account. Since variance differed with plant species richness, each sample was weighted with an inverse of the variance at each plant species richness level. We evaluated the importance of the interaction between soil biodiversity and plant diversity treatments based on model comparison using Akaike information criterion (AIC) that compares the additive model (soil biodiversity + plant species richness) and the interactive model (the additive model with their interaction term).
We calculated multitrophic biodiversity as an average of standardized plant species richness and the number of soil microbial ASVs (Allanet al. 2014), to investigate the effect of multitrophic biodiversity on stability-related indices. Compared with polynomial, quadratic and exponential models in regression analysis, linear models fitted the data better and were employed to test the relationships between multitrophic biodiversity and stability-related indices.
We tested whether soil biodiversity affects interactions among plant functional groups by comparing responses of plant functional groups in the monocultures and species-mixed communities. In the monocultures, a mixed-effect model with harvest time point and plant species as random factors was employed to investigate the effect of soil biodiversity loss on shoot biomass of individual plant functional group. By excluding monocultures, we further tested the effect of soil biodiversity and plant species richness on shoot biomass of plant functional groups using a mixed-effect model with harvest time point as a random factor. These analyses were conducted by sub-setting data during wet and dry periods, separately. All data analyses were performed in the software R4.0.0 (R -Core-Team 2020). For data handling, analysis and visualization, we used the packages dplyr , MASS ,nlme , ggplot2 , patchwork , ggridges ,reshape2 , AICcmodavg and cowplot . Details of data analysis can be found in the Supplementary Methods.