Experimental design
The Trait-Based Experiment23 was established in 2010 within the ‘Jena Experiment’ (www.the-jena-experiment.de) field site, Thuringia, Germany; 50°55’ N, 11°35’ E, 130 m a.s.l. 50. We sampled 34 plots that differed in plant species richness (1, 2, 4 and 8 species) and plant functional trait dissimilarity. The functional trait dissimilarity was based on traits that reflect spatial and temporal resource acquisition strategies. We chose plant height, leaf area, rooting depth, and root length density to reflect spatial resource acquisition. To reflect temporal resource acquisition, we chose growth starting date and flowering onset. The spatial resource acquisition strategy across plant communities ranged from small-statured species with small leaves and shallow roots to tall-statured species with large leaves and deep roots. The temporal resource acquisition strategy across plant communities ranged from early-growing and -flowering species to late-growing and -flowering species. The plots that were sampled, were arranged in three blocks, mown twice per year, and weeded three times per year. For more details on the design, see23.
Secondary metabolome sampling and sample processing
On 24 – 25 August 2015 and 31 May – 01 June 2016, we sampled aboveground biomass of seven common central European grassland species (grasses: Dactylis glomerata L., Holcus lanatus L., Phleum pratense L., forbs: Geranium pratense L., Leucanthemum vulgare (Vaill.) Lam., Plantago lanceolata L., and Ranunculus acris L.). We sampled 444 plant individuals across both sampling seasons and all species, i.e. three individuals per species and plot. One sample was lost during sample processing. We harvested the shoot biomass by cutting the plants ca. 1 cm above the ground and removed all inflorescences. All samples were taken between 15.00 and 19.00 h each sampling day to minimize diurnal variation. All samples were processed, extracted, and analysed according to Ristok et al. (2019)22 with slight changes (see Supplementary Material & Methods for a detailed description).
LC-MS data processing and metabolite prediction
We followed the LC-MS data processing protocol described in Ristok et al. (2019)22 with minor changes (see Supplementary Material & Methods for a detailed description). We predicted metabolite structures through the comparison of LC-MS data with literature references. We submitted high-resolution mass-to-charge values to the MassBank of North America (MoNA, http://mona.fiehnlab.ucdavis.edu/) spectral database. We used a mass tolerance of 0.5 D for comparison. In addition, we calculated high-resolution molecular weights, molecular formulae for putative molecular ions in neutral form, and particle weights for mass spectrometry generated fragments using ChemDraw Ultra 8.0 (www.cambridgesoft.com).
Leaf herbivory rate assessment
For each plant individual, we counted the total number of leaves and the number of leaves with herbivore damage. We noted and categorized herbivore damage for each damaged leaf per plant sample. Herbivore damage included signs of sucking, chewing, and mining on leaves. The damage categories were 1-10%, 10-20%, 20-30%, 30-40%, 40-50%, 50-60%, 60-70%, 70-80%, 80-90%, and 90-100%. We multiplied the number of leaves in each category with the damage level of that category (0.1 for 1-10%, 0.2 for 10-20%, etc.). Next, we took the sum across all categories. Finally, we calculated the herbivory rate for each sample by dividing the summed herbivory by the total number of leaves.
Soil sampling
In each plot, we took soil samples to a depth of 10 cm using a metal corer (diameter 2 cm) on 27 August 2015 and 06 June 2016. We pooled and homogenized five subsamples per plot to account for spatial heterogeneity. We sieved soil samples to 2mm. We stored one part at – 20 °C for phospholipid fatty acid analysis and the other part at 4 °C for nematode extraction.
Phospholipid fatty acid analysis
We determined the soil microbial community via phospholipid fatty acid (PLFA) analysis. We extracted PLFAs following the Frostegård et al. (1991)51 protocol with detailed description in Wagner et al. (2015)52. We analysed all samples on a gas chromatograph (see Supplementary Material & Methods for a detailed description). We used the following PLFA-markers, expressed as ng g-1 dry weight soil, as bacterial markers: (a) gram-negative bacteria: cy17:0 and cy19:0; (b) gram-positive bacteria: i15:0, a15:0, i16:0, and i17:0; and (c) widespread in bacteria: 16:1ω5 and 16:1ω7. We used the following PLFA markers as fungal markers: (a) saprophytic fungi: 18:1ω9t, 18:2ω6t, and 18:2ω6c; and (b) arbuscular mycorrhizal fungi: 20:1 52,53. We summed up all markers within a respective group to receive a representative value.
Nematode extraction and identification
We extracted nematodes from 25 g soil fresh weight using a modified Baermann method52,54. Then, we counted all nematodes at 100x magnification and identified at least 100 randomly chosen nematodes (if available) at 400x magnification using a Leica DMI 4000B light microscope. Nematodes were identified to genus or family level following Bongers (1994)55. We classified all nematodes according to their respective trophic group, i.e. plant feeders, fungal feeders, bacterial feeders, predators, and omnivores. In addition, we assigned all nematodes a c-p score (colonization-persistence gradient) that ranged from 1 to 5 56. Finally, we combined the information of the trophic group and c-p score to create functional nematode guilds57,58.
Statistical analysis
We analysed our data in the statistical software R v3.5 59 (http://www.r-project.org) using the packages ‘lme4’ 60, ‘lmerTest’ 61, ‘effects’ 62, and ‘plspm’ 63
To test for the effects of sown plant species richness and resource acquisition strategy (spatial and temporal), we calculated community mean scores (CMS) based on the original PCA species scores calculated when the experiment was designed23,64. We based our CMS calculations on the relative species-specific cover for each plant community recorded in August 2015 and May 2016, respectively. In short, the six functional traits plant height, leaf area, rooting depth, root length density, growth starting date, and flowering onset were analysed in a standardized PCA. The first PCA axis arranged species according to their spatial resource acquisition strategy. The second PCA axis arranged species according to their temporal resource acquisition strategy23. Plots with high community mean scores on the first PCA axis (CMS_PCA1) were mostly dominated by tall-statured species with deep sparse roots and large leaves. In contrast, plots with negative community mean scores on the first PCA axis contained a high proportion of small-statured species with dense shallow roots and small leaves. Plots with high community mean scores on the second PCA axis (CMS_PCA2) contained mostly late growing and late flowering species64.
We tested our first hypothesis that either the sown plant species richness or any of the resource acquisition strategies affect the plant-individual herbivory rate by calculating linear mixed effects models. We fitted the herbivory rate (log-transformed) as the response variable. As predictor variables, we fitted sampling season (categorical; August 2015 or May 2016), plant functional group identity (categorical; grass or forb), and either sown plant species richness (metric; 1, 2, 4 or 8) or either CMS_PCA1 or CMS_PCA2 (metric), as well as the three-way interaction. We fitted plot nested in block and species identity as independent random effects. Model simplification was achieved by model comparison using the Akaike Information Criterion (AIC). All linear mixed effects models were based on restricted-maximum likelihood estimation and Type I analysis of variance with the Satterthwaite approximation for degrees of freedom.
To test our second hypothesis that plant species richness and plant community resource acquisition traits affect directly or indirectly, via the soil biota, the plant metabolome and thereby herbivory, we calculated partial-least squares path models (PLS-PM)65. We hypothesized direct links from the experimental design variables sown plant species richness and resource acquisition traits to microbial and nematode community, as well as to plant individual biomass, plant metabolome, and herbivory (for a detailed description of each latent variable see Supplementary Table S4). In addition, we hypothesized links from the microbial and nematode community to plant biomass and the composition of the plant metabolome, as well as from plant biomass to secondary metabolome and herbivory. Finally, we hypothesized a link from plant metabolome to herbivory (Fig. 2a). To explore the effect of plant functional group identity, we calculated three separate PLS-PMs: (a) the full model using all data, i.e. across all seven plant species, all plots, and both sampling seasons; (b) the grasses-only model using only data of the three grass species; and (c) the forbs-only model using only data of the four forb species. All data were scaled and we used bootstrapping (n = 200) to calculate confidence intervals for the effect sizes within the path models. Finally, we simplified all path models, i.e. we reduced the number of moderators, until the most parsimonious solution was achieved65. Across all three models, the latent variable resource acquisition traits was characterized by, i.e. moderators with the highest weight relation as part of the linear combination65, the community-weighted means of leaf area, rooting depth, root length density, growth starting date, and flowering start. The latent variable soil microbial community was characterized by gram-negative, gram-positive, and general bacteria, as well as saprophytic The effects of the latent variable nematode community were mostly driven by bacterial feeders with a colonization-persistence score (c-p) of 1 and 3, predators (c-p 4 and 5), fungivores (c-p 4), omnivores (c-p 4), and plant feeders (c-p 2 and 3).