Plant-microbe interactions in response to grassland herbivory and nitrogen eutrophication

Plant-soil feedback is increasingly recognized as a vital framework to analyze multi-trophic interactions involving herbivores, plants and microbes


Introduction
Global change entails significant impacts on terrestrial ecosystems in the form of increased warming and increased loads of reactive nitrogen (N) from atmospheric deposition.A potential consequence of this is significantly larger insect outbreaks and increases of insect herbivory (Netherer and Schopf, 2010;Sommerfeld et al., 2018), due to changes in the secondary metabolites involved in plant defense systems (Throop and Lerdau, 2004;Jamieson et al., 2017), and to enhanced plant palatability.Invertebrate herbivores tend to target certain plant functional groups over others by preferentially feeding on fast-growing species such as grasses before perennial forbs, thereby altering the plant community composition of grassland ecosystems (Huntly, 1991).However, although preferentially fed upon, fast-growing grasses also benefit from herbivory through increased N-cycling and production of litter available for decomposition (Belovsky and Slade, 2000).Thus, under stable conditions plant herbivory induces several mutually antagonistic processes that might synchronously alter plant communities in opposite directions, although increased N loadings may also affect plant community composition and litter production.Plant-herbivore interactions are connected to the soil microbiome as plants and rhizosphere-associated microbiota live in intimate feedback with each other (Tedersoo et al., 2020).These feedbacks help to regulate plant nutrient uptake and plant defense mechanisms (Bever et al., 2012;Babikova et al., 2013).Indeed, accumulating evidence indicates that differing plant functional groups shape soil and root-associated microbial communities directly or indirectly by altering soil conditions over multiple plant generations (Bezemer et al., 2005;Bais et al., 2006;Bahram et al., 2020), and that this constitutes a significant driver of aboveground community responses to herbivory (Bezemer et al., 2005;Kardol et al., 2007;Eisenhauer et al., 2011).
Plants under recent herbivory attack tend to allocate more C into the rhizosphere in the forms of readily accessible C and root exudates (Schwachtje et al., 2006;Züst and Agrawal, 2017), thereby inducing substantial changes in the rhizosphere microbial community composition (Bardgett et al., 1998;Hamilton et al., 2008;Philippot et al., 2013).Over shorter time periods, this often leads to an increase in root-colonizing arbuscular mycorrhizal fungi (AMF) and beneficial soil-borne microbes (Van Wees et al., 2008;Nishida et al., 2009;Pineda et al., 2010;Tawaraya et al., 2012), and to significantly increased N-cycling (Bakker et al., 2004;Shahzad et al., 2015).However, as herbivory continues over time, C allocation into the rhizosphere is gradually reduced (Bardgett et al., 1998), which is expected to result in an overall depletion of rhizobacterial abundance and AMF colonization in line with the 'carbon-limitation hypothesis' (Wallace, 1987).While plant-AMF associations are typically mutually symbiotic, responses of AM plants to herbivory vary markedly according to herbivory feeding mode and diet breadth (Hartley and Gange, 2008;Koricheva et al., 2009).In general, AMF are thought to provide an overall increase in plant resistance to generalist leaf-chewing insects (Gehring and Whitham, 2003;Hartley and Gange, 2008), and to facilitate herbivory from e.g.phloem-sucking aphids (Koricheva et al., 2009).However, the ultimate outcome of herbivory on the root-and soil-associated microbes may depend on edaphic factors such as soil nutrient status (Barber et al., 2013;Bernaola and Stout, 2019) as well as on plant host type and functional group status (Bennett and Bever, 2007;Hartley and Gange, 2008;Koricheva et al., 2009;Middleton et al., 2015).
Prevailing plant-soil feedbacks under herbivory can be altered with increasing N levels in grasslands since deposition of reactive N significantly alters terrestrial ecosystems and plant-herbivore interactions (Matson et al., 2002;Bardgett and Wardle, 2010).In grasslands, N-eutrophication typically shifts microbial community composition by favoring copiotrophic taxa over slower-growing oligotrophic taxa (Fierer et al., 2012;Leff et al., 2015), and by decreased fungal:bacterial biomass ratios (Bardgett and McAlister, 1999;de Vries et al., 2006).Increased N levels may also change herbivore abundances and severity through direct changes in the plants' nutrient stoichiometry (Throop and Lerdau, 2004), or indirectly through induced changes in the plant community composition (Chen et al., 2019).So far, considerable effort has been dedicated to examining the effects of aboveground herbivory and N eutrophication in the plant-soil feedback context separately (e.g.Manning et al., 2008;van der Putten et al., 2016;Heinze and Joshi 2018;Heinze et al., 2020), but we are not aware of any study on the effects of both factors combined (but see Kos et al. (2015)).Results from Borgström et al. (2017) suggested that these factors combined lead to non-additive and unpredictable effects on the plant community composition, but there is scant evidence as to whether the same non-additive effects would prevail over additive, or indeed synergistic effects on the soil microbial community.
Here we used samples from a previous grassland mesocosm experiment (Borgström et al., 2017) to examine plant-soil interactions under herbivory and N eutrophication by analyzing soil bacterial and fungal assemblages as well as mycorrhizal association of two plant functional types.We hypothesized that (i) root AMF colonization would respond negatively to aboveground herbivory in line with the carbon-limitation hypothesis and that this would be further exacerbated by N-eutrophication, (ii) herbivory and N-eutrophication would entail direct but non-additive effects on the composition of soil bacteria and fungi, and on species associated with nutrient cycling and (iii) shifts in belowground microbial community composition would correlate with the aboveground changes in plant community composition observed by Borgström et al. (2017).

Experimental design and sampling
A grassland mesocosm experiment was established and run for three years (2013)(2014)(2015)(2016) by Borgström et al. (2017) in Krusenberg (59 • 44 ′ 27.9 ′′ N, 17 • 41 ′ 02.9 ′′ E), approximately 10 km south of Uppsala, Sweden.Soil at the site is characterized by a very acidic loamy sand of moraine deposits (pH ~ 4) with a clay content of 10-15% and a total C and N content of 1.45% and 0.08% respectively.The original experimental design comprised 64 mesocosms in a full factorial design with a control and three treatment factors with eight replicates: aboveground herbivory, belowground herbivory, and N addition.At the start of the experiment, nine grassland plant species were established in each 1 × 1 × 2 m mesocosm plot (forbs: Achillea millefolium L., Leucanthemum vulgare Lam., Plantago lanceolata L.; legumes: Lotus corniculatus L., Trifolium pratense L.; grasses: Agrostis capillaris L., Dactylis glomerata L., Lolium perenne L., Festuca rubra L.).Results from these experiments showed that the combined impact of herbivory and N-eutrophication resulted in a competitive advantage of forbs over grasses, and shifts to more fungal-driven decomposition as measured by fungal-feeding nematode abundance (Borgström et al. 2017(Borgström et al. , 2018)).Here, we only evaluated the effect of aboveground herbivory by the grasshopper Chorthippus albomarginatus with and without N fertilization on the soil microbial community in relation to the controls, sampling a total of 32 mesocosms in June 2017, four years after the initiation of the experiment.
A detailed account of the experimental design can be found in Borgström et al. (2017).Briefly, each enclosed mesocosm consisted of a wooden frame (dimensions 1 × 1 × 2 m) covered with anti-aphid fine-mesh net (mesh size 0.2 × 0.4 mm).Belowground, the enclosure was situated on a sheet metal base frame (dimensions 1 × 1 × 0.55 m) and covered with the same fine-mesh net at the bottom.The soil at each enclosure was dug out to a depth of 0.5 m prior to the establishment of the samples and were subsequently re-filled with soil from the field prior to planting.The N and combined herbivory and N treatment plots were fertilized yearly with a single application of 15 g of NH 4 NO 3 m − 2 in the form of a solution made from dissolvable pellets (2013) and distilled water or dissolvable pellets only (2014,2015), corresponding to a eutrophication rate of 40 kg N ha − 1 in addition to the estimated background rate for the region of around 8 kg N ha − 1 year − 1 .Thus, additional N was supplied a total of three times at the set rate over the course of the experiment, during the summer months (June or July).In the aboveground herbivory plots, 10 grasshoppers of the species Chorthippus albomarginatus were added in equal male-female proportion at the start of the experiment.Herbivory pressure of the treatment plots was evened out after one year by transferal of nymphs from high-density enclosures to low-density enclosures (for more details on the experimental design, see Borgström et al. (2017).

Soil sampling and chemical analysis
Three soil cores (3-cm-diameter, 10 cm depth) were sampled at random points within each mesocosm and pooled.The soil was stored at − 20 • C immediately after sampling.Prior to chemical analyses, subsamples (5 g) were thawed over 24 h in a fridge.These were subsequently analyzed for extractable P, K, Ca, and Mg through the ammonium lactate method (Egnér 1960, Swedish Standard [Svensk Standard] 1993, 1995) using ICP-AES.Soil pH was measured in a 1:5 soil:water suspension.

Root sampling and preparation
Out of the nine plant species used in the original experiment, root samples from a non-leguminous forb (Leucanthemum vulgare), and a grass (Agrostis capillaris), were collected from the 32 mesocosms.These represented dominant species within their respective functional group (Borgström et al., 2017).Three individual plants of each species were sampled by tracking to the rooting point and digging out intact root systems including surrounding soil with a small spade.The root samples were stored at − 20 • C immediately after sampling.Plant roots were subsequently sampled according to the method outlined by Bulgarelli et al. (2012), in which 3 cm root segments are cut, pooled from all three specimens, and subsequently washed with a sterile PBS buffer (Medicago AB; 140 mM NaCl, 2.7 mM KCl, 10 mM Phosphate buffer pH 7.4) under mechanical (20 min, 180 rpm), and manual shaking.The washing steps were repeated three times to ensure no soil particles were attached to the roots.In contrast to Bulgarelli et al. (2012) no root sonication was used, as this has been associated with root tissue destruction (Reinhold-Hurek et al., 2015;Richter-Heitmann et al., 2016).Roots were then stored in a 50% ethanol solution at 4 • C prior to staining.

DNA extraction, PCR amplification and library preparation
DNA was extracted from 2 g of freeze-dried and milled soil using the PowerMax Soil DNA Isolation Mini kit (Qiagen GmbH, Hilden, Germany) following the manufacturer's instructions.The extracted DNA was quality-checked on the basis of the 260/280 and 260/230 nm wavelength ratios using a NanoDrop™ (Thermo Scientific, Massachusetts, USA) and stored at − 80 • C until sequencing.For production of amplicons for sequencing, the universal prokaryote primers 515F (5 ′ -3 ′ , GTGYCAGCMGCCGCGGTAA) and 926R (5 ′ -3 ′ , CCGYCAATTYMTT-TRAGTTT) were used to amplify the 16S V4 subregion of rDNA (Walters et al., 2016).For eukaryotes, the ITS2 subregion was amplified with the ITS7ngs and ITS4ngsUni primers (Tedersoo and Lindahl, 2016).All primers were further tagged with a 10-12-base identifier code in order to recognize amplicons from each of the samples bioinformatically.DNA samples were amplified using the following conditions in three replicate runs: 95 • C for 15 min, followed by 26 cycles of 95 • C for 30 s, 50 • C for 30 s and 72 • C for 1 min with a final extension step at 72 • C for 10 min.
The 25 μl PCR mix consisted of 18 μl sterilized H2O, 5 μl 5 × HOT FIREPol Blend MasterMix 0.5 μl of each primer (20 μl) and 1 μl template DNA (final concentration of 400 nM).The amplicons from the replicates were pooled, purified using a purification kit containing agarose gel (FavourPrep Gel/PCR Purification mini Kit-300 Preps; Favourgen) and shipped for library preparation in the sequencing service facility of University of Tartu (the Estonian Biocenter).DNA libraries were sequenced on two runs using an Illumina MiSeq platform (2 × 250 bp paired-end chemistry).The raw sequences derived from all soil samples have been deposited at NCBI under accession PRJNA682012.
Known standards of DNA from the bacterial species Pseudomonas aeruginosa and the fungal species Pilidiom concavum with serial dilutions were used for generating standard curves and assays were run using a 96-well format on a CFX Connect Real-Time System (Bio-Rad).Duplicate reactions were prepared for each sample and blanks containing H 2 O instead of DNA template were used as negative controls.The DNA extracts were diluted for the reactions accordingly after running inhibition tests for each sample using M13 primers with a TOPO vector, to ensure absence of PCR inhibition.The qPCR cycling conditions were established as follows: initial DNA denaturation at 95 • C for 5 min, then 40 cycles each with of 95 • C for 15 s, annealing at 55 • C for 30 s, elongation at 72 • C for 45 s, and then signal reading at 80 • C for 5 s, followed by melt curve analyses.

Root AMF colonization
Root colonization of AMF was examined using the grid line intersection method (McGonigle et al., 1990).For this purpose, root samples were cleared in a 2.5% KOH solution, and subsequently stained with trypan blue in an acidic glycerol solution (Koske and Gemma, 1989).The stained roots were then mounted on microscope slides, and the colonization rates of mycorrhizal arbuscules, vesicles (if present) and hyphae were analyzed in a minimum of 100 random fields of intersection under 40X magnification.

Bioinformatics
We used the LotuS pipeline (Hildebrand et al., 2014) to quality-filter, demultiplex, and process the filtered reads into operational taxonomic units (OTUs).Chimera detection and removal was done using Uchime (Edgar et al., 2011), and all singletons and sequences shorter than 100 bp were discarded.Clustering of sequences was done using a de-novo clustering algorithm in UPARSE (Edgar, 2013) with a 97% similarity threshold, and taxonomy was assigned against the SILVA and UNITE databases for prokaryotic and fungal sequences respectively.All OTUs representing archaea (which comprised only 0.59% of OTUs), chloroplasts, eukaryotes and mitochondria were omitted from the prokaryote dataset, and OTUs not assigned to fungi were omitted from the ITS eukaryote dataset after OTU clustering.Contaminants were removed in both datasets.This resulted in a total of 358,101 and 598,188 reads, distributed over respectively 3528 and 1502 OTUs for the prokaryote and eukaryote datasets.The OTUs were subsequently assigned to guild and lineage (eukaryotes) using FUNGuild (Version 1.0; Nguyen et al., 2016) at default settings, and to functional groups (prokaryotes) using FAPROTAX (Version 1.2; Louca et al., 2016).We note that such functional assignments have limitations, as not all functions may be phylogenetically conserved.Relative abundances of OTUs, functional groups and guilds were calculated for each treatment (see below).

Statistical analysis
For analyzing diversity, samples were rarefied to minimum sequencing depths (6398 reads for the bacterial dataset, and 4716 reads for the fungal dataset).On the rarefied data, within-sample diversity (α-diversity) and species richness was measured by the Shannon and Chao1 indices respectively, and differences between treatments were evaluated with a one-way ANOVA and Tukey's HSD for pairwise comparisons.Due to the compositional nature of microbial community data derived from high-throughput sequencing, differences in relative abundances of taxonomic and functional groups between treatments were tested on rarefied data using a non-parametric Kruskal-Wallis test, with post hoc testing using pairwise Wilcoxon rank sum test for all possible group combinations followed by Benjamini-Hochberg multiple testing correction (Weiss et al., 2017).Similar to Carey et al. (2015), differences in microbial relative abundance for all treatments were also expressed in terms of their effect size relative to the control treatment, and calculated as %Effect = (Treatment-Control)/Control.The soil bacteria:fungi ratio was calculated based on results of the soil qPCR measurements, and differences between treatments were tested using one-way ANOVA.
Differences in bacterial and fungal community composition between treatments were tested using PERMANOVA with the adonis function with 999 permutations in the vegan package (Oksanen et al., 2019).As a first step of relating soil microbial community composition to plant functional groups (obtained from Borgström et al., 2017) and environmental variables, we calculated the biomass ratio relationships between all plant functional groups, along with plant evenness (E var ) for all samples.To visualize the relationships, we used global nonmetric multidimensional scaling (GNMDS) in the vegan package based on the T.R. Sveen et al.
following options: two dimensions, initial configurations = 100, maximum iterations = 200, and minimum stress improvement in each iteration = 10 − 7 .Furthermore, the importance of variables in explaining the variation in bacterial or fungal community composition was assessed using PERMANOVA on Hellinger-transformed data.To determine bacterial and fungal classes that may serve to discriminate between treatments, a generalized canonical discriminant analysis (GCDA) was performed on the class-level relative abundance matrix (Hellinger-transformed), using the candisc package (Friendly et al., 2020).
Mycorrhizal colonization was measured as the proportion of root length colonized by arbuscules, vesicles or hyphae, and hence a logistic regression model was used to evaluate the main effects of treatment and plant type and their interaction.Significance of main effects and interactions was examined through likelihood ratio tests (LRT) where models with and without the main effect and interaction were tested.Model validation was done by plotting residuals against fitted values for both models.Significant effects between treatment levels were evaluated with post-hoc tests using the multcomp package (Hothorn et al., 2008), whereas the phia package (De Rosario-Martinez, 2015) was used to assess significant interaction effects.All statistical analyses were carried out using R 4.01 (R Core Team 2019).

Soil properties across treatments
The treatments did not lead to any changes in the soil chemical properties investigated here, apart from a significant difference in K content between the herbivory and combined herbivory and N eutrophication treatment (ANOVA: F 3,28 = 3.4,P < 0.05; Table S1).

Root AMF abundance and colonization
AMF root colonization as assessed through staining and microscopy revealed an overall colonization level ranging between 0 and 8% with some evident outliers (Fig. 1), and significantly lower overall colonization of L. vulgare roots compared to A. capillaris roots (GLM: df = 1,8; χ 2 = 128.3,P < 0.001).Aboveground herbivory entailed no evident impacts on root AMF roots colonization, whereas N eutrophication had a negative impact compared to all other treatments (Tukey Contrasts: df = 1; P < 0.001 for all pairwise comparisons).However, this appeared to stem from the interaction between plant species and treatment (highly significant in all contrasts; Tukey Contrasts: df = 1; P < 0.001), where colonization of A. capillaris roots was reduced under the single impact of N eutrophication but not when combined with herbivory (Fig. 1).

Soil bacterial and fungal diversity, community structure and abundance
The combined effects of herbivory and N eutrophication led to a lower alpha diversity (Shannon's H) compared to each respective treatment for bacteria, but not for fungi (Fig. 2).Similarly, no differences were found in species richness (Chao1) among treatments for either bacteria or fungi.Further, there were no differences in community composition among treatments for either soil bacteria or fungi (Fig. S1; Table S2).Contrary to our expectations, the aboveground changes in plant community composition was not reflected by belowground changes (Fig. S1).Interestingly, PERMANOVA tests showed soil pH and  Ca to be the main determinants of both bacterial and fungal community composition (Table 1).In addition, there was a significant relationship between extractable P and bacterial community composition, whereas fungal communities appeared to correlate with the forb:legume biomass ratio and extractable K (Fig. S1b).There was no difference in bacterial and fungal abundances across the treatments (Fig. S2).Similarly, neither absolute abundances of bacteria and fungi nor the bacteria:fungi ratio showed any significant differences between treatments, but displayed around 100-fold higher bacterial than fungal abundance across the site (Fig. S2).

Differences in relative abundance between phyla
Due to the vast diversity and ecological differentiation in Proteobacteria (Spain et al., 2009;Nemergut et al., 2011), class was used instead of phylum level in the analysis.Dominant taxa in the bacterial dataset were Acidobacteria, Alphaproteobacteria, Actinobacteria and Gammaproteobacteria, together accounting for ~60% of total relative abundance (Fig. 3, S3a).Overall, the combined herbivory and N treatment entailed larger effect sizes on bacterial phyla than each respective single-impact treatment (Table S3), but none of these effects were statistically significant.For the most prevalent phyla, effect sizes were roughly similar in direction and magnitude among treatments.
However, the combined N-eutrophication and herbivory treatment appeared to impact some phyla non-additive negatively (Planctomycetes, Bacteroidetes), whereas others were impacted positively (Chloroflexi, Gemmatimonadetes).Soil fungi were dominated by Ascomycota and Basidiomycota (Fig. S3b), accounting for around 50% and 24% of total fungal reads, respectively.Glomeromycota was the fourth largest phylum and showed increased abundance under N-eutrophication treatments, but a negative effect size under the individual herbivory treatment (Table S3).Glomeromycota followed by Chytridiomycota were the main group discriminating herbivory treatment samples from those from control or N-eutrophication (Fig. 3).

Differences between guilds and functional groups
The classification of bacteria into functional groups highlighted some further impacts of herbivory and nitrogen on soil bacteria.Out of the 3528 OTUs, 929 (26%) were assigned to 32 functional groups, the dominant of which are shown in Fig. 4 and Table S4.Chemoheterotrophy was the dominant metabolic strategy in all bacterial communities, with aerobic chemoheterotrophy-associated bacteria accounting for 56-62% of functional group relative abundance (Fig. 4, Table S4).Compared to control and N-eutrophication treatment, herbivory and more strongly its combined effect with N-eutrophication reduced the relative abundance of bacteria associated with fermentation (P < 0.05) and nitrate reduction (P < 0.1), compared to all other treatments (Fig. 4).In addition, there were some indications that herbivory mediated a negative response of ureolytic bacteria to N addition, as well as increasing relative abundances of nitrogen-fixing bacteria while decreasing abundances of ammonia-oxidizing bacteria (Table S4).
Of the 1502 OTUs in the soil fungal dataset, a total of 702 (47% of OTUs and 41% of total reads) were assigned to functional guilds, and the remaining 800 (53%) were of unknown guild.The assigned OTUs were then further filtered into four main guilds stemming from unambigious guild-annotations in FUNGuild, yielding a total of 426 assigned OTUs at this level (Table S4).Saprotrophs was the most dominant group in all treatments, accounting for 34-38% of guild relative abundance, consisting of phyla that responded differently to the treatments (Fig. 4, Table S4).Thus, soil bulk AMF was the guild with the lowest relative abundances, accounting for 12-14% of total guild relative abundances, and showed a positive response to N-eutrophication but appeared negatively impacted by herbivory (Table S4).

Root AMF colonization
In this study we demonstrate that the AMF response, as measured by root colonization, to herbivory and N eutrophication is contingent on plant species and on the interactions of these.Host identity is known as an important determinant of AMF abundance and diversity (Scheublin et al., 2007;Becklin et al., 2012;Lugo et al., 2015), but recent evidence highlights large differences also between individual plants of a given species (van de Voorde et al., 2010;Lekberg and Waller, 2016), which was reflected also in the large within-treatment variance of our results.It is thus plausible that a higher number of replications, both within plant species but also with additional species for each plant functional group, would have increased the differences observed, as variation in colonization between individual plants is typically high.We nonetheless observed how N eutrophication drastically reduced colonization of A. capillaris roots, but not those of L. vulgare, whereas the inclusion of herbivores lessened this negative impact.Therefore, any boost in AMF colonization in line with the "cry-for-help" hypothesis -indicating that plants promote associations with mutualistic microbes under stressed conditions (Rudrappa et al., 2008) -or a decrease in line with the carbon-limit hypothesis (Wallace, 1987) needs to be considered in the light of plant species and perhaps functional types.In contrast with our initial hypothesis of a synergistic negative impact on root-colonizing AMF stemming from an aboveground generalist invertebrate herbivore in combination with N eutrophication, we instead found that herbivory entailed a positive mediation of the otherwise negative impacts caused by N eutrophication, but only for A. capillaris.Interestingly, the only larger meta-analysis conducted on the topic of AMF colonization and herbivory by Barto and Rillig (2010) found herbivory to entail reduced colonization of perennial grasses but not of forbs.They found, however, an overall increase in AMF colonization for both plant types when grown together in mixtures albeit with considerably fewer studies examined.Similarly, N-eutrophication at lower loads than those used here has been linked to reductions in AMF colonization across a range of ecosystems (Egerton-Warburton and Allen, 2000;Treseder, 2004;Mohan et al., 2014), and the indication seen here that aboveground herbivory might attenuate this negative response for some plant species or functional types thus warrants further research.For example, it remains to be investigated what underlying mechanisms drive the attenuating effect of herbivory.Although herbivory may reduce plant nutrient uptake, it may also result in compensatory growth (Fornoni, 2011), which requires more nutrients and perhaps an increased C-allocation belowground.
One potential explanation for the apparent positive mediation that herbivory entailed on AMF colonization of A. capillaris under N eutrophication observed here could lie in the shifts in biomass composition of plant functional groups observed at this study site (Borgström et al., 2017), and in the symbiotic role AMF can play in providing the host plant with phosphorus (Smith et al., 2011).As expected, Borgström et al. (2017) observed an increase in the biomass of both plant functional types under N eutrophication.Similarly, the herbivores showed a preference for feeding on grasses, which resulted in a concomitant decrease of grass biomass but increase in forb biomass, indicating a shift of competitive advantage.However, under the combined effects of herbivory and N eutrophication this shift was less pronounced and the root to shoot ratio of the whole plant community increased significantly, suggesting that much of the aboveground biomass lost across the entire plant community was compensated by stable or slight increases in the denser root systems of grasses (Kutschera, 1960).Greater root biomass and root:shoot ratio of grasses have been shown to correlate well with plant tolerance to herbivory (Strauss andAgrawal 1999, Stowe et al., 2000;Tiffin, 2000), and the induced re-location of biomass from exposed to inaccessible areas is a well-known plant defense mechanism in the face of herbivory (Orians et al., 2011;Thomas et al., 2017).Furthermore, despite the acidic soil of our experimental site, soil extractable P remained at the lower range of the 50 mg P kg − 1 considered optimum for plant growth in European grasslands (Janssens et al., 1998;Critchley et al., 2002).Thus, it is possible that the relatively stronger grass-AMF association observed here could stem from an increased plant uptake of P necessary to maintain the belowground biomass while suffering from losses of shoot biomass to herbivory.Ceulemans et al. (2019) found a negative relationship between both N eutrophication and P availability and AMF richness in a study of 40 European grassland sites.Their finding does not in principle contradict our results (Fig. 1) since the proliferation of a few abundant AMF species at the expense of many less abundant is an often-reported consequence of N eutrophication (e.g.Egerton-Warburton and Allen, 2000).The close correlation between P and the grass:forb biomass ratio observed here (Fig. S1a), along with the seemingly antagonistic impact of the combined elevated N and herbivory treatment on forb biomass observed by Borgström et al. (2017) at the study site would further corroborate this hypothesis.Nevertheless, as only one species per functional group was included in this study, further investigation is certainly needed to confirm whether this observation is due to functional groups rather than species-specific differences.
In summary, we found that herbivory moderated the negative effects of N eutrophication on AMF colonization, in line with an increased response to insect herbivory (Pineda et al., 2010;Tawaraya et al., 2012;Chen et al., 2017), but only for the studied grass species.This is at odds with our first hypothesis of additive negative effects of both impacts, and could have potential implications when predicting diversity and functional responses of grasslands to accelerated global change, although we do acknowledge that including more species of each respective plant functional group would have made the results more generalizable.

Shifts in microbial composition
Differences in microbial community composition in response to the individual or combined impacts of herbivory and N-eutrophication were smaller than expected.Interestingly, the inclusion of N-eutrophication induced significant decrease in the intrinsic diversity of the herbivoryaffected samples (Fig. 2a), without any concomitant significant differences in any dominant phyla (Fig. 3, S3).In a recent study by Eldridge et al. (2017), herbivory was shown to indirectly increase bacterial alpha diversity by suppressing dominant taxa on the phylum level, whereas globally, N-eutrophication appears to either decrease diversity by affecting certain phyla or to entail no discernible impact (Koyama et al., 2014;Leff et al., 2015).Thus, while the higher alpha diversity related to herbivory observed here would be expected, it is nevertheless surprising that herbivory in combination with N-eutrophication would lead to smaller within-sample diversity.
A fundamental tenet of the plant-soil feedback framework are the intertwined processes of biotic and abiotic factors as drivers of feedbacks.We expected to see differences in the soil microbial community composition reflecting shifts in metabolic strategies previously observed under N-eutrophication (Fierer et al., 2012;Ramirez et al., 2012), to interact with the often non-additive effects observed under herbivory (Bakker et al., 2004;Bardgett and Wardle, 2010).Therefore, the absence of differences in microbial community composition in the face of both individual and combined treatment effects was similarly surprising.In this context, a question arises: how well can yearly single applications of N purport to reproduce actual N deposition patterns, and what responses would these trigger?N deposition rates varies significantly across both months and years (Waldner et al., 2014).Yet, individual applications timed with rainfall or confined to the growing season appears to be a rule rather than exception when examining responses to simulated N deposition (e.g.Magill et al., 1997;Sinsabaugh et al., 2005;Eisenlord and Zak 2010;Frey et al., 2014).This could be of importance if the sampling is conducted in the immediate aftermath of application, as microbial communities tend to exhibit short-term responses to N addition that are not sustained over time (Fauci and Dick 1994;Stark et al., 2007).In the case of our study, however, sampling was conducted one year after the last of three N applications, after which any lasting responses would have presumably been incorporated into changes of the soil microbial community.
Overall, compared to fungal communities, bacterial communities responded more strongly to treatment effects, and some significant changes in the composition of functional groups could be discerned (Table S4).Interestingly, bacteria associated to fermentation and NO 3 reduction exhibited seemingly antagonistic responses to herbivory when combined with N eutrophication (Fig. 4).Increased N deposition, of which NO 3 is the dominating form, has been shown to increase nitrate reductase in grasslands (Morecroft et al., 1994) analogous to fertilization events (Herzog et al., 2015), whereas the opposite response is common in the face of environmental or herbivory-induced stress (Mattson, 1980;Louda and Collinge 1992).Here, the evident decreases in the bacterial functional group associated with NO 3 reduction indicates that a top-down control on the soil microbial community stemming from the significant increases in herbivory under N eutrophication observed by Borgström et al. (2017) might be dominating over any bottom-up effects originating from increased microbial N cycling.This could have significant ramifications for nutrient cycling and plant diversity in N-driven grasslands (Ceulemans et al., 2013) in the context of global change and increased N deposition along with larger insect outbreaks and herbivory, and the top-down control exerted might be additionally strengthened by a hotter climate and reductions in soil moisture (Ritchie, 2000).As for the similar treatment effects on bacteria associated with fermentation, we have not been able to link this to any known process relevant to the design used here, and further studies are Despite little variation across treatments, soil pH was significantly related to bacterial and fungal community composition across all samples (Table 1), which follows in line with a number of studies over large scales and natural settings (Lauber et al., 2009;Shen et al., 2013;Bahram et al., 2018).Soil pH, including in acid environment such as here, is a key driver and confining factor of microbial diversity and functioning in a variety of nutrient-poor ecosystems (Stark et al., 2014;Maestre et al., 2015;Zhalnina et al., 2015), and has been shown to exert a controlling influence over the diversity between soil microbial communities (Siciliano et al., 2014;Zeng et al., 2016).Further, the impacts of N-eutrophication on soil microbial community structure is sometimes only an indirect effect transmitted through alterations of the soil pH (Zhang et al., 2017), in line with the comparable levels of soil pH and microbial communities observed between treatments here (Table S1).
We similarly expected that the changes in plant functional composition observed in Borgström et al. (2017) in response to the individual aboveground herbivory and N eutrophication treatments at this study site would translate into changes in the soil microbiome, as has been previously shown in a number of greenhouse studies (Bezemer et al., 2005;Ma et al., 2017;Pineda et al., 2020).Yet, although aboveground herbivory significantly reduced the biomass of grass to the benefit of forbs and legumes and despite this pattern being somewhat modified by the inclusion of N-eutrophication (Borgström et al., 2017), these changes appeared to have a weak impact on the soil microbiome composition (Table 1).Zeglin et al. (2007) showed similar results for three grassland sites in the USA under N-enrichment when looking at microbial functioning, and pointed to the overall constraining factor that pH has for the soil microbiome despite significant changes in the aboveground plant community.Yet, as was shown by Marshall et al. (2011), and by Strecker et al. (2016), changes in the soil microbiome as a consequence of disturbances are slow to build, often taking several years or even decades to properly reflect.It is thus plausible that the observed changes in plant functional community could have more profound effect on soil microbiota over a longer period of time.

Conclusions
We found that aboveground herbivory and N-eutrophication have surprisingly comparable impacts on soil microbial community structure and potential functions.Joint impacts of these treatments were often stronger than each individually, suggestive of additive effects of biotic and abiotic disturbances on microbial communities.We also found that shifts in plant-mycorrhizal associations in response to disturbances depend on plant species, likely reflecting differential reliance of plant functional groups on their associated fungi for nutrient acquisition under stressed conditions.Our data suggest that the diversity and relative abundances of most microbial taxa persist in a disturbed local community over a short-term study even in the light of changes in the aboveground, and corresponding belowground, composition of plants.Yet, certain taxa and functional groups showed significant responses to the combined effect of herbivory and N-eutrophication, suggestive of possible consequences of such disturbances for some key taxa over shortterm periods.In addition, the variation in soil and vegetation across treatments explained a significant proportion of variation in the structure of soil microbiome.Thus, it is likely that in long-term, changes in such parameters exerted by herbivory and N-eutrophication would indirectly impact the soil microbiome.

Fig. 1 .
Fig. 1.Root AMF colonization of Agrostis capillaris and Leucanthemum vulgare associated with the plant roots across treatments (n = 8).Treatments: C = control, H = aboveground herbivory only, H.N = combined aboveground herbivory and N eutrophication and N = N eutrophication only.

Fig. 2 .
Fig. 2. Diversity and richness across treatments for a) bacteria, and b) fungi.Diversity and richness were calculated based on Shannon's H diversity index and Chao1 species richness index on rarefied data, respectively.** denotes significant differences (P < 0.05) based on pairwise Wilcoxon comparisons with Benjamin-Hochberg corrections for multiple testing.