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.