Plant diversity effects on herbivory are related to soil biodiversity and plant chemistry

Insect herbivory is a key process in ecosystem functioning. While theory predicts that plant diversity modulates herbivory, the mechanistic links remain unclear. We postulated that the plant metabolome mechanistically links plant diversity and herbivory. In late summer and in spring, we assessed individual plant above‐ground herbivory rates and metabolomes of seven plant species in experimental plant communities varying in plant species diversity and resource acquisition strategies. In the same communities, we also measured plant individual biomass as well as soil microbial and nematode community composition. Herbivory rates decreased with increasing plant species richness. Path modelling revealed that plant species richness and community resource acquisition strategy correlated with soil community composition. In particular, changes in nematode community composition were related to plant metabolome composition and thereby herbivory rates. Synthesis. These results suggest that soil community composition plays an important role in reducing herbivory rates with increasing plant diversity by changing plant metabolomes.


| INTRODUC TI ON
Insect herbivory is an essential ecosystem process that can remove substantial amounts of biomass in grasslands within a single season (Meyer et al., 2017;Seabloom et al., 2017). Plant diversity can influence the abundance and diversity of insect herbivores (Haddad et al., 2001;Hertzog et al., 2016) as well as herbivory rates (Ebeling, Meyer, et al., 2014;Wan et al., 2020). Plant traits, in particular those associated with resource acquisition and competition, are considered to provide mechanistic links between plant diversity and herbivory (Loranger et al., 2013). These traits can reflect functional aspects of light interception, deep soil nutrient and water use, or resource use along a seasonal gradient (Ebeling, Pompe, et al., 2014;Marr et al., 2021). Whereas these traits may predict the performance of plant species in their niches, they also display plasticity in response to plant soil feedbacks and soil legacies (Delory et al., 2021;Xi et al., 2021). However, a recent study shows that the commonly used morphological and physiological traits only explain up to 12.7% of the variance in herbivory (van der Plas et al., 2020).
Plant chemical composition may be a better predictor for individual herbivory, because many herbivores use plant metabolites to locate their host (Agrawal & Weber, 2015) while plants use metabolites to defend themselves (van Dam & van der Meijden, 2011). By using the plant's metabolome, that is, the composition of all metabolites produced by an individual plant (Oliver et al., 1998), as an additional functional plant trait, we may gain deeper insights in the molecular mechanisms underlying differences in herbivory.
Several factors may explain differences in herbivory rates across plant diversity gradients. Higher levels of plant diversity may increase niche diversity by increasing spatial heterogeneity and the variety of food sources, thus, supporting more insect herbivores and increasing community-level herbivory rates (Ebeling, Meyer, et al., 2014). At individual plant level, however, increased plant diversity may lead to dilution effects which decrease herbivory, as it will be more difficult for specialized insect herbivores to localize their host plant (Castagneyrol et al., 2014;Finch & Collier, 2000).
Indeed, in a previous study conducted in the Jena Experiment, individual herbivory decreased with increasing plant species richness (Scherber et al., 2006). Lastly, the abundance of predatory and parasitoid arthropods, which can reduce herbivore populations and thus plant community-level and individual herbivory via top-down control, is commonly higher in more diverse plant communities (Haddad et al., 2009;Hines et al., 2015;Schuldt et al., 2019;Wan et al., 2020).
In addition, differences in herbivory across plant species richness gradients may also be explained through changes in plant chemistry.
Plant metabolomes change in response to (a)biotic variation. This metabolomic response to environmental conditions co-determines the defensive status of a plant (van Dam & van der Meijden, 2011).
For instance, plants increase the synthesis of defensive metabolites following an attack by herbivores (Bezemer & van Dam, 2005;Karban & Baldwin, 1997). These induced responses can change defences both locally, that is, in the attacked tissue, and systemically, that is, throughout the plant (van Dam & Heil, 2011). In addition to herbivory-induced changes, plant diversity itself can affect plant metabolomes. Plant-plant interactions can alter the metabolome through competition, which may induce the production of volatile (Baldwin et al., 2006) and nonvolatile allelopathic compounds (Fernandez et al., 2016). Seen the broad biological activity spectrum of plant metabolites, these changes are likely to affect herbivory rates (Broz et al., 2010). Lastly, soil legacy effects, which may result from systemically induced changes triggered by soil biota, such as microbes and nematodes (van Dam & Heil, 2011;Wondafrash et al., 2013) can also affect plant metabolomes (Ristok et al., 2019).
Taken together, the plant metabolome both affects and reflects above-and below-ground interactions with insect herbivores, other plants and soil biota, in a species-specific and context-dependent way (Bezemer & van Dam, 2005;Ristok et al., 2019). Hence, we argue that measuring plant metabolomes will provide novel insights into the relationship between plant diversity and above-ground herbivory.
The aim of our study was to jointly analyse the relationships between plant diversity, soil biota communities, plant metabolomes, and above-ground herbivores to provide a mechanistic framework for above-below-ground multitrophic interactions in grasslands.
Here, we analysed the metabolomes and individual plant herbivory of three grass and four forb species in experimental plant communities manipulated to vary in spatial or temporal resource acquisition traits (Ebeling, Pompe, et al., 2014). Our species selection covered a range of functional traits related to resource acquisition (Ebeling, Pompe, et al., 2014) and included both grasses and forbs, because their metabolomes and response to the abiotic and biotic environment may differ (Dietz et al., 2019(Dietz et al., , 2020Huberty et al., 2020).
All plants were grown in 34 experimental plant communities that varied in plant diversity, that is, species richness and functional trait diversity (Ebeling, Pompe, et al., 2014). We tested if and how plant diversity alters the secondary metabolome and how this relates to herbivory. We hypothesized that (1) plant species richness and the resource acquisition strategy of the plant community affect individual plant herbivory. Moreover, we calculated partial-least-squares path models to explore if (2) plant species richness and plant community resource acquisition traits directly or indirectly, via the soil biota, relate to the plant's metabolome and thereby may explain variation in herbivory. Our hypotheses are based on observations that plant species richness affects soil community composition  and that differences in soil biota can affect the plant's metabolome and thereby herbivory (Huberty et al., 2020;Ristok et al., 2019). Additional paths in our models accounted for relationships between the soil microbial community and soil nematode community (Dong & Zhang, 2006) as well as for potential direct relationships between plant species richness and individual plant herbivory rates (e.g. due to dilution effects; Castagneyrol et al., 2014;Scherber et al., 2006). We inferred similar paths for functional trait composition, accounting for observations in which plant communities containing tall-statured species with large leaves and deep roots increased individual plant herbivory, as these species may have provided more niches for insect herbivores (Loranger et al., 2012).
We also modelled the relationship of growth and flowering time with herbivory as plant chemistry is known to change with ontogeny (Barton & Koricheva, 2010;Boege, 2005). Lastly, we included the relationship between plant biomass and herbivory in our path models, because soil community composition can affect plant biomass, which in turn may affect herbivory, whereby larger plants may incur more herbivory (Windig, 1993).
We show that increasing plant species richness reduces individual plant herbivory. Furthermore, our path models suggest that soil community composition, especially the composition of the nematode community, and plant individual metabolomes are key players in the relationship between plant diversity and above-ground herbivory.

| Experimental design
The trait-based experiment (Ebeling, Pompe, et al., 2014) was established in 2010 within the 'Jena Experiment' (www.the-jenaexper iment.de) field site, Thuringia, Germany; 50°55′N, 11°35′E, 130 m a.s.l. (Roscher et al., 2004; see Appendix S1 for details). We sampled 34 plots that differed in plant species richness (1, 2, 4, and 8 species) and plant functional trait dissimilarity (see Table S1 for information on the plant community of each plot). 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. All plots were arranged in three blocks, mown in June and September, and weeded three times per year.

| Secondary metabolome sampling and sample processing
We sampled twice under different environmental conditions to account for seasonal variation in the plants' metabolomes. Initially, we sampled above-ground biomass of eight common central European grassland species (grasses: Anthoxanthum odoratum L., Furthermore, we excluded seven samples (4 samples of P. pratense, 2 samples of R. acris, 1 sample of G. pratense) due to contamination during sample processing (i.e. the final number of analysed samples = 443). We harvested the shoot biomass by cutting the plants c. 1 cm above 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) with slight changes (see Appendix S1). In short, we extracted 20 mg dried ground plant tissue of each sample in 1 ml of extraction buffer (methanol/50 mM acetate buffer, pH 4.8; 50/50 [v/v]). The samples were homogenized for 5 min at 30 Hz using a ball mill (Retsch mixer mill MM 400), and subsequently centrifuged (20,000 g, 10 min, 4°C). The supernatant was collected in a 2 ml Eppendorf tube. We repeated the extraction procedure with the remaining pellet and combined the supernatant with the first one. We centrifuged (20,000 g, 5 min, 4°C) all extracts, transferred 200 μl to an HPLC vial and added 800 μl extraction buffer, resulting in a 1:5 dilution.
Metabolites were analysed on a liquid chromatography quadrupole time-of-flight mass spectrometer (LC-qToF-MS; Bruker maXis impact HD; Bruker Daltonik) with an electrospray ionization source operated in negative mode (Appendix S1).

| LC-MS data processing and metabolite prediction
The LC-MS data are presented as a list of features described by mass-to-charge ratios, retention times, and intensities. We processed LC-MS data as in Ristok et al. (2019) with minor changes (see Appendix S1). 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.fiehn lab.ucdav is.edu/) spectral database. We used a mass tolerance of 0.5 D for comparison. Furthermore, 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.cambr idges oft.com).

| 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 6 June 2016. We pooled and homogenized five subsamples per plot to account for spatial heterogeneity. We sieved soil samples to 2 mm. We stored one part at -20°C for phospholipid fatty acid analysis and the other part at 4°C for nematode extraction.

| Nematode extraction and identification
We extracted nematodes from 25 g fresh soil using a modified Baermann method (Ruess, 1995;Wagner et al., 2015). We counted all nematodes at 100× magnification and identified at least 100 randomly chosen nematodes (if available) at 400× magnification using a Leica DMI 4000B light microscope. Nematodes were identified to genus or family level following Bongers (1994). We classified all nematodes into plant feeders, fungal feeders, bacterial feeders, predators and omnivores. Moreover, we assigned all nematodes a c-p score (colonization-persistence gradient) that ranged from 1 to 5 (Bongers & Bongers, 1998). Finally, we combined the trophic group and c-p score to create functional nematode guilds as a proxy for nematode community structure Ferris et al., 2001).
Based on earlier studies in the same experiment (Beugnon et al., 2019;Steinauer et al., 2017), we calculated community mean scores (CMS) to represent resource acquisition strategy (spatial and temporal). We based our CMS calculations on the original PCA species scores calculated when the experiment was designed (Ebeling, Pompe, et al., 2014;Fischer et al., 2016), and 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 strategy (Ebeling, Pompe, et al., 2014). Plots with high community mean scores on the first PCA axis (CMS_PCA1) were mostly dominated by tall-statured species with deep 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 species (Fischer et al., 2016).
We tested our first hypothesis by calculating linear mixed effects models. We fitted herbivory rate (log-transformed) as the response variable. As predictor variables, we fitted sampling campaign (categorical; August 2015 or May 2016), plant functional group identity (categorical; grass or forb), and either plant species richness (metric; 1, 2, 4 or 8) or either CMS_PCA1 or CMS_PCA2 (metric), as well as the two-way and three-way interactions.
We tested for the overall and pairwise differences in shoot metabolome composition among the different sown plant species richness levels by calculating permutational multivariate analyses of variance using distance matrices. We log +1 transformed the metabolite intensity data to achieve multivariate normality, and used Bray-Curtis dissimilarity to calculate the distance matrices. All analyses were permuted 9999 times. Each analysis was species-specific and sampling campaign-specific. We were not able to calculate pairwise comparisons of the metabolome composition between plants grown in monoculture (lowest plant species richness level) and in the highest diversity plot (8 species mixture). This is due to the experimental design (see Table S1). For each species, only one monoculture plot was present. In addition, there was only one 8-species plot. This meant there were not enough replicates to run permutational multivariate analyses of variance and, as such, the pairwise comparisons between monoculture and the 8-species plot were excluded from the analyses.
In addition, we calculated two metrics of metabolite diversity: (a) the richness of secondary metabolites, that is, the number of metabolites within a plant individual; and (b) the Shannon diversity of secondary metabolites, that is the abundance-weighted diversity of metabolites expressed as the exponential of the Shannon-Weaver index (Hill, 1973) based on plant individual-level metabolite intensities. We also calculated community-weighted mean (CWM) trait values for each trait considered in the design of the Trait-Based-Experiment . Here, we based the calculations on the relative species-specific cover for each plant community.
We calculated linear mixed effects models to test for the effect of sown plant species richness or CWM trait values on the richness or Shannon diversity of secondary metabolites. We fitted either the richness or the Shannon diversity of secondary metabolites as response variables. As predictor variables, we fitted sampling campaign, plant functional group identity, and either sown plant species richness or each of the CWM traits separately (metric, scaled), as well as the two-way and three-way interactions.
Finally, we tested for the relationship between richness or Shannon diversity of secondary metabolites and plant-individual herbivory using linear mixed effects models. We fitted herbivory rate (log-transformed) as the response variable. As predictor variables, we fitted sampling campaign, plant functional group identity, and either richness or the Shannon diversity of secondary metabolites, as well as the two-way and three-way interactions.
For all linear mixed effects models calculated in our study, we fitted plot nested in block and species identity as independent random effects. We performed backwards model simplification, first removing nonsignificant interactions and then nonsignificant predictors, until the change in Akaike Information Criterion (AIC) was <2. Finally, the most parsimonious model with the lowest AIC was chosen. All linear mixed effects models were based on restrictedmaximum likelihood estimation and Type I analysis of variance with Satterthwaite approximation for degrees of freedom.
We hypothesized direct links from the experimental design variables plant species richness and resource acquisition traits to microbial and nematode community composition (De Deyn & Van der Putten, 2005;Strecker et al., 2016), as well as to plant individual biomass , plant metabolome (Scherling et al., 2010), and individual plant herbivory (Scherber et al., 2006; for details on latent variables see Table 1). In addition, we hypothesized links from the microbial to the nematode community composition (Dong & Zhang, 2006) as well as from either soil biota community to plant biomass (van der Putten et al., 2013) and to the composition of the plant metabolome (Huberty et al., 2020;Ristok et al., 2019). Furthermore, we hypothesized links from plant biomass to metabolome (de Jong, 1995;Fernandez et al., 2016) and herbivory (Barnes et al., 2020). Finally, we hypothesized a link from plant metabolome to herbivory (van Dam & van der TA B L E 1 Description of the latent and observed variables used in the PLS-PM models

Latent variable Description
Resource acquisition traits The latent variable represents the community-weighted mean (CWM) trait values of maximum plant height, leaf area, rooting depth, root length density, growth start, and flowering start. We based the calculations of the CWM trait values on the relative species-specific cover for each plant community. For instance, a positive relationship of resource acquisition traits with the latent variable biomass means higher values in the CWM traits correlated with a higher plant biomass

Microbial community
The latent variable represents the microbial biomass for gram-negative bacteria, gram-positive bacteria, undefined bacteria, saprophytic fungi, and arbuscular mycorrhizal fungi. For instance, a negative relationship of microbial community with the latent variable biomass indicates that a higher biomass of soil bacteria and fungi correlated with a lower individual plant biomass Nematode community The latent variable represents the relative abundance of functional nematode guilds. A functional nematode guild is the combined information of trophic guild and colonizer-persistence score. For instance, a positive relationship of nematode community with the latent variable metabolome means a greater relative abundance of certain functional guilds correlated with a higher concentration of metabolites

Metabolome
The latent variable represents the abundance of secondary plant metabolites. For instance, a negative relationship of metabolome with the latent variable herbivory means that a higher abundance of metabolites is correlated with a lower herbivore damage on individual plants

Species richness
The observed variable represents the plot level sown plant species richness from 1 to 8

Biomass
The observed variable represents the above-ground dry biomass of individual plants

Herbivory
The observed variable represents the herbivory rate on individual plants Meijden, 2011; Figure 4a). We calculated three separate PLS-PMs:

| The effects of plant species richness and resource acquisition strategy
Herbivory rates decreased significantly with plant species richness (Table S2; Figure 1a) and were lower in May 2016 than in August 2015 (Table S2; Figure 1). Together, plant species richness and sampling campaign explained 26% (marginal R 2 value, hereafter R 2 marg ) of the total variation in herbivory rates.
We found no significant effect of the community's spatial resource acquisition strategy on herbivory rates (Table S2). In other words, the relative abundance of tall-statured species with large leaves and deep roots did not significantly affect the plant-individual herbivory rates (Figure 1b). However, when we tested for the effect of community spatial resource acquisition strategy within plant functional groups, we found a significant effect for grasses (Table S3). Grasses growing in communities that predominantly contained tall-statured plants with large leaves suffered lower herbivory rates ( Figure S1a; R 2 marg = 0.35). In contrast, we observed no significant effect of spatial resource acquisition strategy on individual herbivory in forbs (Table S3; Figure S2a). When we tested for the effect of the temporal resource acquisition strategy of the plant communities on herbivory rates, we only found a marginally significant relationship (Table S2). More specifically, plant-individual herbivory tended to be greater in communities containing mostly later growing and flowering species (Figure 1c). Separate analyses for each plant functional group showed that this effect was significant for grasses (R 2 marg = 0.36; Table S3; Figure S1b), but not for forbs (Table S3; Figure S2b).
We also tested for the effects of plant species richness on plant metabolome composition and on metabolite diversity, that is, richness or Shannon diversity of secondary metabolites. We found a significant effect of plant species richness on the metabolome composition across most plant species in at least one sampling campaign, except for Holcus lanatus ( Table 2). Consecutive pairwise comparisons revealed that the metabolome of plants grown in monocultures most frequently differed from the metabolome of plants grown in more diverse communities (Table S4).
We did not find an effect of plant species richness on the richness of secondary metabolites (Figure 2a). Rather, we observed an effect of plant functional group identity on the richness of secondary metabolites (F 1,5 = 6.69; p = 0.049; Table S5a). Forb species had significantly more secondary metabolites (396 ± 8.1; mean ± SE) than grass species (311 ± 8.5; mean ± SE). Together with sampling campaign, plant functional group identity could explain 68% (R 2 marg ) of the total variation in secondary metabolite richness. We also discovered that the effect of plant species richness on the Shannon diversity of secondary metabolites depended on the functional group identity (F 1,138 = 5.35; p = 0.022; Table S5a). Increasing plant species richness increased the Shannon diversity of secondary metabolites in grasses, while it reduced the Shannon diversity of secondary metabolites in forbs (Figure 2b). Sampling campaign, plant functional group identity, and plant species richness together explained 49% (R 2 marg ) of the total variation in the Shannon diversity of secondary metabolites.
Moreover, we analysed the extent to which resource acquisition strategy affected metabolite diversity. We found that some traits associated with spatial resource acquisition strategy can increase or decrease the richness of metabolites dependent on sampling campaign and functional group identity (Table S5) (Table S6).
We observed a similar significant effect of plant height on the richness of metabolites in forbs (Table S6; (Table S6). In addition, we discovered that the community-weighted mean of leaf area, a trait associated with spatial resource acquisition, and the community-weighted mean of flowering onset, a trait associated with temporal resource acquisition, had similar effects on the Shannon diversity of metabolites (Table S5)

TA B L E 2
Overall differences in the plant species-specific shoot metabolome composition between plant species richness levels in different sampling campaigns. The plant species richness levels were 1, 2, 4, and 8 species. Statistical parameters resulting from permutational multivariate analyses of variance (PERMANOVA) using distance matrices. We used bray-Curtis dissimilarity matrices and 9999 permutations. Significant differences (p < 0.05) are given in bold community composition was positively related to plant metabolomes, which itself was negatively related to plant individual herbivory (for details of each latent variable see Table 1; for all direct, indirect, and total path coefficients see Table S7). Our most parsimonious model predicted 34% of the total variation in the secondary metabolome, and 22% of the total variation in individual herbivory. Plant species richness was negatively correlated with the relative abundance of predatory, omnivorous and plant feeding nematodes ( Figure S3).
The spatial resource acquisition trait plant height was positively correlated with the relative abundance of predators, fungivores, omnivores, and plant feeders. In contrast, leaf area was negatively correlated with the relative abundance of bacterial feeders, fungivores, omnivores and plant feeders. Rooting depth and root length density negatively correlated with bacterial feeders, predatory nematodes and omnivores, but positively correlated with plant feeders.
Conversely, the temporal resource acquisition traits growth starting date and flowering start were negatively correlated with the relative abundance of plant feeders ( Figure S3). In addition, we extracted the 100 most important metabolite mass spectra that characterized the metabolome, that is, the metabolites with the strongest positive correlation with the latent variable 'metabolome'. We could assign molecular formulas and structures to 13 mass spectra (Table S8; Figures S4-S16). These metabolites were mainly phenolics, their precursors or their derivatives, which are all products of the shikimic acid pathway. Moreover, these compounds are known to respond to phytopathogenic nematode infection (Ohri & Pannu, 2010) and play a role in plant-herbivore interactions (Whitehead et al., 2021).
As part of the nematode community composition, the relative abundance of bacterial feeders, predators, omnivores, and plant feeders showed the strongest positive correlations with the concentration of the assigned metabolites. Especially sinapic acid, a flavonol, the chlorogenic acid dimers, and quinic acid, were negatively correlated with plant herbivory ( Figure S17).
Our full-model PLS-PM also indicated that plant herbivory was positively related to community-weighted resource acquisition traits and plant individual biomass (Figure 4b). Plant height, growth starting date, and flowering start were most strongly positively correlated with the latent variable 'resource acquisition traits'. Neither trait was individually correlated with plant herbivory ( Figure S3).
However, our path model suggests a synergistic effect on plant herbivory, that is, plant communities of tall-growing species with late growth and flowering start may increase plant individual herbivory.
Plant biomass was negatively related to resource acquisition traits, plant species richness, nematode community composition, and plant metabolome. Lastly, the microbial community composition was positively related to plant species richness and negatively related to resource acquisition traits. We performed sensitivity analyses and calculated two alternative full-model PLS-PMs: (a) a path model F I G U R E 2 The richness and Shannon diversity of secondary metabolites in response to (a, b) sown species richness across all samples, as well as in response to selected community-weighted mean traits (c-e, i, j) in grasses only, and (f-h, k, l) in forbs only of the Trait-Based Experiment. For clarity, the placement of the symbols corresponding to the plant functional group identity in panels 'a, b' have been slightly shifted along the x-axis. In panels 'a, b' the relationship in grasses is displayed in circles and orange colour, and the relationship in forbs is displayed in triangles and green colour. In panels 'c-l' the relationship in August 2015 is displayed in circles and red colour, and the relationship in May 2016 is displayed in triangles and black colour. Significant relationships are displayed by solid lines. Nonsignificant and marginally significant relationships are displayed by dashed lines. Regression line estimates are based on linear mixed effect models with plot nested in block, and species identity as independent random terms (Tables S5 and S6). The shaded band displays the standard error. CWM, community-weighted mean. ( (Table S2). The shaded band displays the standard error.

F I G U R E 4
Hypothesis-based conceptual partial-least-squares path model (a) as well as path model including data across both sampling campaigns and all plant species (b), only across all grasses (c), and only across all forbs (d). Species richness represents the plot-level sown plant species richness. Resource acquisition traits represent the community-weighted mean traits maximum plant height, leaf area, rooting depth, root length density, growth starting date, and flowering start. Microbial community represents PLFA-based estimates on plot-level gram-negative, gram-positive, and undefined bacteria, as well as arbuscular mycorrhizal fungi and all other fungi abundance. Nematode community represents plot-level summed relative abundance of functional nematode guilds, that is, bacterial-feeding, carnivorous, fungalfeeding, omnivorous and plant-feeding. Biomass represents plant-individual above-ground dry biomass. Metabolome represents plantindividual secondary metabolite composition. Herbivory represents plant-individual herbivory rate expressed as the proportion of damaged leaves to the total number of leaves. All data is scaled. Variables taken at plot level are highlighted by a grey-shaded background. Variables taken at the plant-individual level are highlighted by a white-shaded background. Black arrows display significantly positive relationships. Red arrows display significantly negative relationships. Numbers on arrows are path coefficients. Numbers in the round boxes display the explained variation (R 2 ).
that links herbivory to plant metabolome, which would account for herbivore-induced responses (GoF = 0.14; Figure Table S7). In the forbs-only model, the plant metabolome was negatively related to nematode community composition, but positively related to biomass and herbivory. Moreover, the resource acquisition traits were not related to any other latent variable, plant individual biomass was not related to herbivory, but plant species richness was negatively related to herbivory (see Table S7). Both functional group models, however,  van Dam & Heil, 2011) our study yields novel insights by highlighting how below-ground communities may shape plant metabolomes, thereby becoming a significant driver of above-ground herbivory.
The abundance, diversity and community structure of soil biota are commonly determined by the species identity and traits of individual plants as well as the plant community diversity (Bezemer et al., 2010;Lange et al., 2015;Strecker et al., 2016). Accordingly, our path model showed that plant species richness and variation in resource acquisition-related functional traits can explain variation in soil microbial and nematode community composition. The relationship between plant species richness and microbial community composition is likely due to an increased and more diverse influx of organic matter in the form of rhizodeposits Lange et al., 2015;Steinauer et al., 2016). The observed negative relationship between resource acquisition traits and most functional nematode guilds was mainly driven by community-weighted plant height, growth starting time, and flowering onset. This suggests that the abundance and seasonality of resource influx from the plant community into the soil determines nematode community structure (Yeates, 1999). In contrast, rooting depth and root length density were positively correlated with phytophagous nematodes, suggesting that phytophagous nematodes can also be affected by root architecture (Yeates, 1999). Changes in soil community composition were related to significant changes in plant metabolomes.
Specifically, the abundance of bacterial feeders, predators and phytophagous nematodes positively correlated to the concentration of defence-related metabolites in individual plants. Bacterial-feeding nematodes contribute to the mineralization of nitrogen in the soil, which supports plant growth and potentially the synthesis of defence-related metabolites (Freckman & Caswell, 1985). Predatory nematodes control plant parasitic nematodes, thus also indirectly supporting plant growth (Freckman & Caswell, 1985 Similarly, also root herbivory may impact shoot metabolomes (Bezemer & van Dam, 2005).
While an aim of our study was to test if the plant's metabolome can explain variation in herbivory, we could not disentangle potential effects of herbivory on the plant's metabolome. Because we analysed field plants, it is likely that the metabolomes result from several (a)biotic interactions. If anything, this adds realism to our results, as in nature above-ground herbivores often encounter plants induced by other interactors (van Dam & Heil, 2011). In all models, the relationship between biomass and metabolome was maintained.
The negative relationships found in the full and grasses model, may support the hypothesis that larger plants produced less defence, because they can tolerate biomass loss to herbivory and prioritize growth over defence production (de Jong, 1995). In forbs, however, there was a positive relationship between biomass and metabolome.
This might point to the fact that flowering forbs are commonly larger and produce more and different metabolites to protect their reproductive organs (McKey, 1979). Additional experiments are necessary to test the hypothesis that the growth-defence trade-off varies with ontogeny and between plant functional groups.
Overall, sampling campaign had a strong effect on plantindividual herbivory, the richness and Shannon diversity of secondary metabolites, and the relationship among these variables as well as with resource acquisition traits. While this potential seasonal effect is certainly interesting, our experimental and sampling design does not allow for mechanistic or causal interpretation. Our data are limited because (1) we have no repeated seasonal measurements, (2) we sampled at the end of one growing season and at the beginning of the next growing season, which could have led to differences in herbivore community and herbivory between sampling campaigns (Meyer et al., 2017), and (3) leaf traits and the plant's metabolome can vary within and between years (Peters et al., 2018). Dedicated experiments that repeat sampling throughout the season (e.g. Marr et al., 2021) and in multiple consecutive years are necessary to analyse the seasonal variation in the relationship between the plant's metabolome, resource acquisition traits, and herbivory.
We also discovered that the metabolite diversity in grasses and forbs varied differently to changes in resource acquisitionassociated community-weighted traits. These contrasting responses are likely due to differences in defensive strategies. Grasses possess silica crystals providing mechanical protection from herbivory (Massey & Hartley, 2009), while forbs invest in carbon-based defences, such as phenolics (Cooke & Leishman, 2012;Larson, 1988).
Moreover, grasses and forbs differ in their associations with soil biota, such as the symbiosis with mycorrhizal fungi, which can contribute to diverging metabolomic responses (Chialva et al., 2018;Ristok et al., 2019). While the difference between grasses and forbs was not the focus of our study and we only analysed a small subset in each functional group, our results stress the importance of including functional group identity to improve predictive models analysing plant-herbivore interactions.
While the present experiment provides novel insights into above-below-ground relationships, additional experiments should manipulate and disentangle the individual and interactive roles of plant and soil biodiversity in driving changes in plant metabolomes and herbivory rates (Peters et al., 2018). Such studies should be conducted in the presence and absence of above-ground herbivores (Seabloom et al., 2017), to assess if above-ground herbivory modulates plant and soil biodiversity effects on the metabolome, for example via induced responses (Peters et al., 2018). Preferably, these studies should include specialist and generalist herbivores as well as different feeding types (Mithöfer & Boland, 2008).
Taken together, the present study provides support for the existence of tight relationships between plant diversity, soil biota communities, plant metabolomes, and above-ground herbivores. Our results especially suggest that the plant metabolome is an important functional trait (Walker et al., 2022) that can aid to explain more variation (22%) in herbivory than commonly used morphological and physiological traits (on average 12.7%; van der Plas et al., 2020).
By including metabolomic analyses, we advanced our knowledge on the potential mechanisms linking plant diversity and herbivory rates via changes in plant metabolomes (Peters et al., 2018). In addition, we highlight that the soil nematode and microbial communities shape above-ground interactions and that season and plant functional group identity should be considered when analysing such relationships. Our study creates a framework for future experimental research which can further illuminate the underlying mechanisms through targeted and independent manipulation of the plant, soil biota, and herbivore community. It thereby expands our capability to better characterize the complex nature of multitrophic interactions above and below the ground.