Short-term effects of non-grazing on plants, soil biota and aboveground-belowground links in Atlantic mountain grasslands

Mountain grasslands in the Iberian Peninsula are the result of extensive grazing. However, a progressive abandonment of grazing activity is currently observed in the study region. The objective of this work was to evaluate the short-term (2 years) effects of non-grazing on the diversity and composition of plants, soil microorganisms (prokaryotes, fungi, arbuscular mycorrhiza), mesofauna, macrofauna and aboveground-belowground links, through the study of 16 grazed vs. non-grazed areas in Atlantic grasslands located in the Basque Country (Spain). Sites were divided between 4 habitat types with different elevation, pasture productivity, vegetation type and parent material. Herbivores appeared to influence plant community composition, contributing to increase aboveground diversity, while having unequal effects on belowground communities depending on the organisms analysed. This may be explained by the different habitat and trophic level of each soil organism, which may be more or less affected by the predominating negative effects of grazing, such as soil compaction, and only partially compensated by other positive effects. Finally, habitat type appeared to be the strongest influence on both above- and belowground communities, also influencing the effect of the absence of grazing.

In mountain grasslands of the Iberian Peninsula, grazing has been the most important economic activity since the Neolithic 1 . These semi-natural grasslands deliver many ecosystem services, such as primary productivity and nutrient cycling, food production, preservation of genetic diversity, water storage, pest regulation, pollination, recreation, etc. However, noteworthy changes in management practices have taken place in these extensive livestock systems, including shorter grazing time in higher mountain pastures, lower stocking rate and less shepherd control 2 . These trends have also been detected in many parts of Europe; if they continue, a progressive mediumor long-term grazing abandonment is expected 3 . In particular, in the Basque Country, management practices are leading to the abandonment of grazing activity in many mountain areas.
Belowground responses to livestock grazing are less widely investigated than aboveground responses 2,4,5 . In this respect, soil microbial communities are known to play a crucial role in soil processes and the delivery of vital ecosystem services 6,7 ; besides, they can influence plant community composition and productivity 8 . Among soil fungi, arbuscular mycorrhiza are thought to help plants tolerate grazing 9 . On the other hand, the activity and composition of soil microbial communities can also be influenced by (i) predator-prey interactions with small soil invertebrates and (ii) habitat formation due to macrofaunal activity 10 . These interactions are strongly influenced by, and may react in unison to, the presence of large herbivores. Herbivores exert this influence by (i) regulating the quantity and quality of organic matter and resource inputs, (ii) altering the functional composition of vegetation and (iii) physical disturbance 11 . These effects are likely to change depending on the edaphoclimatic conditions present in the grazed area, which, in turn, affect the structure and function of the biological communities present in each site. Unfortunately, there is little research and data integration regarding above-belowground responses and links, and few studies have simultaneously investigated the effect of grazing/non-grazing on plant species composition, soil macrofauna, mesofauna, microorganisms, and soil abiotic properties 12 . In this respect, one of the most significant challenges currently faced by ecologists is probably the enormous complexity and heterogeneity of the soil foodweb and its micro-habitats 13 , although technological advancements in DNA sequencing have opened up for large-scale studies of microbial community structure.
The main objective of this study was to assess the short-term (2 years) effects of non-grazing on plant communities, soil biological communities and the aboveground-belowground links between these two communities in Atlantic mountain grasslands. We hypothesized that cessation of grazing would either decrease or increase belowground activity and diversity in different ways for different soil organisms. Increased diversity may indicate negative grazing effects from compaction, while a decrease may indicate positive effects from e.g. increased availability of labile carbon due to fecal matter from grazers and increased production of plant exudates due to compensatory growth 11 . Which of these effects will dominate would depend on the habitat and trophic position of the community studied. Climatic and edaphic factors will also influence the effects 12 , with less fertile soils being less likely to show positive effects from grazing 11 .
Specific questions investigated in this study were: (Q1) Does cessation of grazing affect plant communities and various soil biological communities differently? (Q2) Are the links between these communities different in the presence vs. the absence of grazing? (Q3) Do the abovementioned effects vary depending on bedrock and habitat type, with a special focus on the fertility of the site? (Q4) Do traditional methods, such as community-level physiological profiles (CLPPs) and arbuscular mycorrhizal (AM) spore morphotyping, provide similar responses as those offered by molecular techniques, such as DNA sequencing, regarding the effects studied? To this end, we compared grazed to non-grazed sites in an Atlantic grassland area, performing a comprehensive analysis of soil biological communities through the determination of (i) basal and substrate-induced respiration; (ii) types of macrofauna; (iii) the "Qualitá Biologica del Suelo" -Soil Biological Quality (QBS) from mesofauna communities; (iv) community-level physiological profiles (CLPPs) of cultivable bacteria through Biolog Ecoplates TM ; (v) diversity of morphotypes of AM spores; and (vi) prokaryotic and fungal diversity through 16 S rRNA and ITS amplicon sequencing, respectively. In total, grassland sites belonging to four distinct categories were studied (i.e., mountain and valley grasslands, each on siliceous as well as on calcareous bedrock) in order to cover the variability existing in the studied area.

Results
Relatively infertile mountainous sites with high soil acidity and a pH between 3.9 and 5.5 dominated the sites studied, while soil pH ranged between 5.6 and 6.9 in more fertile valley sites. None of the chemical parameters related to soil nutrient levels (i.e., soil organic matter-SOM, total nitrogen, Olsen phosphorus and extractable potassium) nor pH were significantly affected by grazing (Table 1). However, soil compaction decreased significantly (39% in average) in the absence of grazing and, concomitantly, infiltration time significantly decreased (48% in average). In turn, soil humidity significantly increased (12% in average) in the absence of grazing (Table 1).
Regarding biological parameters (Table 2), root depth and root abundance were not significantly affected by grazing. However, the number of macrofauna types and CO 2 emission values increased significantly in the absence of grazing (39% and 36% in average, respectively), while the mesofauna QBS index decreased significantly (15% in average).
The diversity indexes for all community composition datasets are listed in Table 3. In all plots, a total of 73 different plant species, 22,219 prokaryotic OTUS, 8,037 fungal OTUs and 18 AM spore morphotypes were identified. Plant richness and Shannon's diversity index decreased significantly in the absence of grazing, while prokaryotic richness instead increased significantly. Other diversity indexes, i.e. fungal and glomeromycota diversity (from ITS amplicon sequencing data), CLPPs of bacteria and AM spore morphotypes, did not change significantly. However, when performing this analysis considering only the more infertile mountainous sites, the richness of AM spore morphotypes increased significantly in the absence of grazing (t-value = 2.212, p-value = 0.038). Instead, when considering only the valley sites, Shannon's diversity from CLPPs of bacteria decreased significantly (t-value = −2.777, p-value = 0.039).
Supplementary Figure S1 shows the distribution of the 20 most abundant (A) plant, (B) prokaryotic and (C) fungal taxa detected. Plant communities appeared to vary considerably both between mountain vs. valley plots and between non-grazed vs. grazed plots (Supplementary Figure S1A) Table 1. Physicochemical parameters (average ± standard deviation; n = 4). In the right (italics), summary of the linear mixed-effect models. G = grazed; NG = non-grazed.
( Supplementary Table S1). Instead, prokaryotic communities varied less between habitats at the order level and none of the top 20 orders showed significant differences when comparing grazed vs. non-grazed plots (Supplementary Figure S1B). Fungal communities were less even at order level (Supplementary Figure S1C) and the most abundant taxa differed greatly depending on habitat. For example, Archaeorhizomycetales constituted 36 and 33% in 6230a and 6230c mountain habitats, respectively, but only 6 and 4% in 6170 mountain and valley plots, respectively. Further, Mortierelalles constituted 23, 21 and 27% in 6170, 6230c and valley plots, respectively, but only 4% in 6230a habitat. Finally, Agaricales constituted 11 and 8% in 6170 and valley plots, respectively, but only 1% in 6230a and 6230c habitats. From these top 20 orders, only Chaetothyriales and Eurotiales showed significantly higher relative abundance values in the absence of grazing. In order to study the effect of grazing at the community composition level, pCCAs were performed with datasets of plants, CLPPs of bacteria, AM spore morphotypes, and amplicon sequence data of 16S rRNA and ITS, including grazing as an explanatory variable and site as a covariable. Out of these analyses, only that using plant community composition as response variable showed statistically significant effects according to the Monte Carlo permutation test ( Supplementary Fig. S2). Therefore, we proceeded to study the effect of other variables besides grazing (such as habitat and bedrock) on community composition. Valley habitats (compared to others) accounted for 50% of the explained variation in plant composition, while the Type 6170 mountain habitat accounted for 34% of the explained variation (pseudo-F = 7.3, p-value = 0.002; Fig. 1a). These two habitats presented a fairly different plant community composition. Instead, the effect of grazing, even if statistically significant, was limited to 10% of the explained variation. Both for prokaryotes ( Fig. 1b; pseudo-F = 23.7, p-value = 0.002) and fungi ( Fig. 1c; pseudo-F = 9.4, p-value = 0.002), most of the explained variation in community composition data was also attributable to habitat type. As with plant communities, the valley and Type 6170 habitats accounted for the largest shares of the explained variation. Grazing was not significant in explaining variation for neither prokaryotic nor fungal community composition. Further, RDAs indicated higher prokaryotic and, to a lesser extent, fungal richness in valley plots.
Supplementary Figure S3 show similar RDA plots performed with Glomeromycota community composition at the species level, CLPPs of bacteria and AM spore morphotype abundances. These three community datasets coherently agree with the results obtained with the broader and deeper amplicon sequence analyses, pointing mainly at valley habitat as the variable explaining most of the data variation and as the habitat that showed higher richness values. Also, grazing contributed significantly to 15% to the variation of CLPPs of bacteria obtained from Biolog EcoPlates TM .
Other soil physicochemical and biological parameters analyzed here are shown in Fig. 2a (pseudo-F = 10, p-value = 0.002). Most of the variation in these parameters was due to type of habitat, although also grazing  Table 2. Biological parameters (average ± standard deviation; n = 4). In the right (italics), summary of the linear mixed-effect models. SIR = substrate-induced respiration; G = grazed; NG = non-grazed.  Table 3. Biodiversity indexes in different community datasets (average ± standard deviation; n = 4). In the right (italics), summary of the linear mixed-effect models. G = grazed; NG = non-grazed.
regime and type of bedrock were significant variables. Soil respiration (basal and substrate-induced), total nitrogen, SOM, root depth and abundance, and soil humidity positively correlated with the 6230a mountain habitat, and pH was negatively correlated. When including all community diversity indexes as response variables ( Fig. 2b; pseudo-F = 9.6, p-value = 0.002), the type of habitat (valley in particular) accounted for the greatest variability, but grazing was also significant. All belowground diversity indexes calculated here, except for H' of AM spore morphotypes, were positively associated towards the positive region of RDA-1 and the valley habitat. Instead, plant diversity indexes were positively associated towards the positive region of RDA-2, the mountain 6170 habitat and the grazed sites.
The compositions of almost all community datasets correlated strongly with each other, as indicated by Mantel tests (Supplementary Table S2). Further, the strength of correlations did not change to a great extent depending on whether only grazed or non-grazed plots were considered. However, if considering diversity indexes instead of community composition, only belowground communities correlated significantly to each other ( Table 4). The only exceptions to this was plant richness correlating with prokaryotic richness in grazed plots and plant richness with AM spore morphotype richness in non-grazed plots. In general, diversity indexes obtained from CLPPs prokaryotic community composition at the order level, and (c) fungal community composition at the order level, as response variables, with grazing, bedrock and habitat as explanatory variables. Only the 15 response variables with best fit are shown. Variation explained by each axis is indicated between brackets. The explanatory variables that appear are those that significantly explained the variation in response data following forward selection; the contribution percentage of each variable is shown between brackets. correlated significantly with those obtained from amplicon data. With few exceptions, the diversity of AM spore morphotypes did not correlate significantly with any other parameter.

Discussion
Effect of grazing on soil parameters, diversity indexes and composition (Q1). Two years after establishing the non-grazed areas, no consequences were observed at the soil nutrient level. However, the reduction of soil compaction in the absence of grazing was very clear, explaining its higher water infiltration capacity. In fact, soil compaction was severe in some grazed areas (i.e., a pressure of over 3 MPa was necessary to penetrate the soil; data not shown). This suggests that the effect of grazing was largely negative, in accordance to Bardgett & Wardle 11 for relatively infertile soils. Moreover, a higher rate of CO 2 emissions was observed in non-grazed areas. This is also consistent with a previous study were grazing abandonment was also simulated by fencing; after five years, soil compaction decreased, values of enzyme activities and microbial biomass were reduced, and CO 2 emission values increased 2 . Jensen et al. 14 also detected a decrease in CO 2 emissions due to compaction in a grazed soil and in a soil dedicated to corn production. In our case, this increase in non-grazed areas could be due to both increased circulation of water and gases through the soil and the accumulation of plant necromass. However, Klumpp et al. 15 showed that grassland systems heavily disturbed through grazing foster soil carbon loss.
In our study, macrofauna richness increased in non-grazed areas, but the mesofauna QBS index decreased. Related results of grazing-induced effects on soil biota abundance and activity can be found in the literature: in a study by Bardgett et al. 16 , the absence of grazing decreased the abundance of microarthropods. Bardgett and Leemans 17 also showed that two years without grazing reduced soil microbial biomass and activity, due to the cessation of liming on ungrazed grassland, but also due to the removal of readily utilizable substrate in the form of sheep excreta and urine and changes in the quantity and quality of root exudates.
The increase in percentage cover of Erica vagans in the absence of grazing may indicate an incipient succession towards heathland already after two years of exclosure. Apart from that, the increase of the unpalatable Brachypodium pinnatum in ungrazed grasslands has been known for a long time (e.g., Lousley 18 ). Besides, this species belongs to a group of grasses that can dramatically reduce biodiversity as it colonizes the territory 19 . This Figure 2. Biplot of the redundancy analysis performed with (a) physicochemical and biological parameters measured, and (b) all diversity indexes calculated, as response variables, with grazing, bedrock and habitat as explanatory variables. Variation explained by each axis is indicated between brackets. The explanatory variables that appear are those that significantly explained the variation in response data following forward selection; the contribution percentage of each variable is shown between brackets. might partly explain the lower plant richness and diversity in the absence of grazing, as well as the lower percentage cover of species like Agrostis curtisii, Trifolium repens, Danthonia decumbens and Carex caryophillea.
Despite the decrease in plant diversity, prokaryotic richness increased in the absence of grazing as did the abundance of some fungal taxa (e.g., Chaetothyriales and Eurotiales). By contrast, Aldezabal et al. 2 observed no significant differences in bacterial richness between excluded and grazed plots. However, these authors used DGGE community profiling, which is generally less sensitive to changes in community diversity than the method used here. The increase in microbial diversity observed here may be related to the reduction of soil compaction and the higher values of water infiltration capacity found in the absence of grazing. Another reason could be the cessation of vegetation removal and, hence, the increase of litter accumulation in non-grazed areas. Changes in organic matter dynamics and nutrient supply are likely to have a profound influence on microbial community structure, as well as on plant physiological responses to grazing such as changes in root exudates or plant litter quality 20,21 .
As opposed to alpha diversity, we could not detect any significant effect of grazing on soil prokaryotic and fungal community composition. Soil is very heterogeneous, and each microorganism is circumscribed by its own immediate environment 22 . It might be that indirect changes that occur at a macroscale, such as the modification of the grazing regime, do not affect to such a great extent to these organisms living in microhabitats. Also, the spatial resolution of soil biological responses is higher than that of aboveground plant responses, implying that a more intensive sampling strategy might be needed to reveal microbial responses. Besides, two years of absence of grazing might be a too short period of time to produce consistent changes on soil microbial community  Table 4. Pearson's correlations to compare richness (S) and Shannon's (H') diversity for every composition dataset. Bold: statistically significant results (p-value < 0.05).
composition. The strong influence of habitat types may also have obscured more subtle changes in microbial composition in our study. From a meta-analysis, Andriuzzi and Wall 12 found that the responses of soil biota abundance and activity to herbivores vary by (i) climate, (ii) ecosystem type and (iii) herbivore identity. Our contribution to this list of variables would be (iv) the type of soil organisms studied. According to the results obtained here, we can speculate that the mesofauna are less negatively affected by compaction and the absence of plant litter than the macrofauna. In respect to the microbial communities, it seems that herbivore activity affected more bacterial communities than fungal ones in general, as their richness increased in the absence of grazing. This may be related to the fact that highly grazed grasslands are dominated by bacteria 15 . Links between aboveground and belowground prokaryotic and fungal communities (Q2). In terms of diversity, there was a strong correlation between belowground prokaryotic and fungal communities (consistently with Lanzén et al. 23 ). However, plant and prokaryotic richness only correlated in grazed plots. This might indicate a stronger dependence of plants on belowground prokaryotes under a stress situation, such as the presence of grazing by large herbivores.
As opposed to alpha-diversity, a strong compositional correlation was found among all communities. Several authors have found a correlation between plant and soil microbial communities, such as Millard and Singh 24 who suggested that while plant community composition can drive fungal community composition, bacterial community structure is more influenced by SOM quality, hence driven indirectly by plant community composition via plant biomass.
Finally, the strength of these links in terms of alpha-diversity and composition between all studied communities was not found to be affected by grazing.

Effect of other variables on the measured parameters (Q3).
Rather than grazing, habitat type (and divergence between valley and mountain habitats) better explained differences in community structure. Valley habitat plots were consistently more diverse, especially for prokaryotes. Further, fungal communities (Archaeorhizomycetales, Mortierelalles and Agaricales orders in particular) differed greatly depending on habitat type. This is expected, since aboveground plant and belowground microbial communities are structured according to the edaphoclimatic conditions of a particular site 25 .
Interestingly, the different responses observed to the absence of grazing in the highly acidic mountain sites (increased richness of spore morphotypes) and in the more productive valley sites (reduced Shannon's diversity from CLPPs) agree with the general predictive framework proposed by Bardgett & Wardle 11 , in which grazing increase soil biological activity in fertile ecosystems, and decrease it in less fertile ones. The significantly increased richness of spore morphotypes found here, and the rather subtle effects detected for other fungal measurements, may also support the hypothesis that predicts stronger aboveground effect on specialist soil organisms 26 such as AM fungi, that are intimately associated with plant roots.

Complementarity of methods (Q4).
The determination of CLPPs through Biolog EcoPlates TM appeared more sensitive for detecting community changes caused by grazing abandonment. This may be because it is a method that can detect functional and, therefore, more ecologically relevant aspects of soil microbial communities. This is consistent with a study by Yang et al. 27 , using the functional gene array GeoChip 4.0, which indicated that the functional structure of soil microbial communities was very sensitive to the impact of livestock grazing. However, a disadvantage of Biolog TM is the exclusion of organisms that cannot be cultivated, and a bias to fast-growing opportunists that can. However, we recommend the use of Biolog EcoPlates TM to complement the molecular determinations of soil microbial diversity, as they have previously shown repeatability, discriminating power, and sensitivity for many environmental factors 28 .
Determination of diversity of AM spore morphotypes (H' in particular) showed divergent results compared to those obtained with amplicon sequencing. Few studies have concurrently analysed mycorrhizal diversity in different propagule types 29 . It was not surprising that amplicon data would provide a very different picture of Glomeromycota diversity, since it includes all DNA in a soil sample including which would form spores as well as mycelium, whereas morphological characterisation focuses only on sporulation during a certain season 30 . Morphological characterisation of arbuscular mycorrhiza is a time-consuming, expertise-requiring technique that, based on our results, we cannot recommend to complement molecular methods since it did not add any complementary insights into the communities studied.

Conclusions
In response to Q1, the cessation of grazing affected plant communities, but also various communities of soil organism in a different way, probably due to the different habitat and trophic level of each community. In the short term, the absence of large herbivores decreased plant diversity and significantly changed its species composition (increasing the Erica vagans shrub and the unpalatable Brachypodium pinnatum, and decreasing others). Meanwhile, although an increase of macrofauna richness, prokaryotic richness and the abundance of some fungal taxa (Chaetothyriales and Eurotiales), as well as a decrease of the mesofauna QBS index, was detected, the overall belowground microbial community composition did not appear significantly affected. This might be because soil microorganisms operate at much smaller spatial scales and are only indirectly affected by grazing.
In general, the links between plant communities and soil biological communities did not differ in the presence vs. the absence of grazing (Q2). Prokaryotic and plant diversity significantly correlated in grazed but not in non-grazed plots, although grazing did not appear to modify the strength of the links between aboveground and belowground communities. There was an influence of bedrock and habitat type on the observed effects, since more beneficial and negative effects of the absence of grazing were observed in the more fertile and infertile sites, respectively (Q3). Habitat type (especially mountain vs. valley) appeared to represent the strongest influence on both above-and belowground communities.
Finally, our results support the use of Biolog EcoPlates TM as a complementary technique to more advanced molecular approach when studying soil microbial communities (Q4).
As a general lesson, and with the final objective of elucidating the consequences of disturbances such as grazing or its absence for ecosystem functioning, the results obtained here bring out the need to consider the complex interactions between aboveground and belowground communities and the effects produced on different soil taxa.  (Table 5). Four permanent exclosures of 10 × 10 m were placed in each site (Table 5).

Study area and experimental design. This study was carried out in Gorbeia
In situ vegetation and soil analyses. Two years after the establishment of the exclosures, plant richness and percentage cover were determined by randomly locating a 50 × 50 cm square (5 times per sampling area; an average of these measurements was used hereinafter). Root abundance and root depth were measured in situ following Mijangos et al. 32 . Compaction was measured using a digital penetrometer (Rimik CP40II) to record the pressure needed to penetrate the soil over a 0-75 cm depth profile. Water infiltration capacity was determined by quantifying the time needed to absorb 230 ml of water poured into a 10 cm diameter tube fixed to previously water-saturated soil. Soil CO 2 emissions were determined with an infrared gas analyser (EGM-4, PP Systems, Amesbury, MA, USA) linked to a cylindrical soil respiration chamber (SRC-1, PP-Systems) as described in Biau et al. 33 . Finally, morphotypes of macrofauna (e.g., worms, cockroaches, woodlice, millipedes, centipedes, earwigs, ants, termites, grasshoppers, beetles, true bugs, spiders, snails, cicadas) were also counted within a block of soil 30 cm deep × 25 cm long × 25 cm wide 23,32 . All these measurements were performed inside and outside (in an area of similar size immediately next to the exclosures) the 16 exclosure areas (32 plots in total). Soil sampling and analyses. On the same day that in situ measurements were taken, a composite soil sample of topsoil (0-10 cm) was collected using a core soil sampler (25 mm diameter; 15-20 cores taken at random throughout the whole sampling area) and sieved to < 2 mm (except for the determination of mesofauna). Subsequently, soil humidity was determined. For the measurements of physicochemical parameters, soil samples were air-dried at ambient temperature until constant weight. Soil organic matter content, pH, total nitrogen, Olsen phosphorus and extractable potassium were measured according to standard methods 34 .
For biological parameters, soils were stored fresh at 4 °C for a maximum of one month until analysis. Basal and substrate-induced respiration (SIR) were measured following ISO 16072 Norm 35 and ISO 17155 Norm 36 , respectively. The mesofauna QBS index, where microarthropods are separated according to convergent morphological features to evaluate the level of adaptation to the soil environment, was calculated following Parisi et al. 37 . CLPPs of cultivable heterotrophic bacteria were determined with Biolog Ecoplates TM,38 , where measurements corresponding to an incubation time of 48 h were chosen for further calculations, following Epelde et al. 39 . Arbuscular mycorrhizal (AM) spores were isolated from 50 g dry weigh (DW) of soil from each sample by wet sieving and decanting. After manually selecting all the AM spores present under a microscope (40x), they were mounted on slides in either polyvinyl alcohol -lactic acid -glycerol (PVLG) or a mixture of PVLG -Melzer reagent. Spore characteristics were observed under a microscope at 200x and 400x for the establishment and quantification of morphotypes.
Soil sub-samples for molecular analysis were stored at −20 °C. DNA extraction and amplicon library preparation were carried out as described in Lanzén et al. 23 . Briefly, PowerSoil DNA Isolation kits were used for extraction (Mo-Bio Laboratories, Carlsbad, CA, USA); primers 515F (CTGNCAGCMGCCGCGGTAA) and 806R (GGACTACHVGGGTWTCTAAT) used for amplification of prokaryotic communities, and ITS1F (CTTGGTCATTTAGAGGAAGTAA 40 ) and ITS2R (GCTGCGTTCTTCATCGATGC 41 ) for fungi, using a dual-indexed adapter-linked protocol based on Caporaso et al. 42 . Library preparation and Illumina MiSeq V2 sequencing was carried out at StarSEQ GmbH, Mainz, Germany. Sequence data has been deposited to the European Nucleotide Archive with study accession PRJEB9654, under sample group "Biopasto".
Sequence data analysis. Sequence data processing and taxonomic classification was performed as described in Lanzén et al. 43 . Briefly, overlapped read-pairs from 16 S rRNA amplicons were quality-filtered and overlapped using usearch 44 and truncated to a length of 253 nt, and ITS read-pairs trimmed to remove the reverse primer using cutadapt 45 . All quality-filtered overlapped sequences from 16S rRNA and ITS amplicons were clustered into OTUs using a 3% divergence cut-off using vsearch 46 . Representative prokaryotic and fungal OTU sequences were aligned to the SilvaMod v106 and UNITE reference databases, respectively, using blastn (v.2.2.25 + task megablast), and taxonomically classified using CREST with default parameters 47 . The relative taxon abundances derived by CREST at order rank were used in further analysis. This approach was motivated by the fact that sequence-based relative abundance has been demonstrated to provide meaningful semi-quantitative information when comparing community structure between samples 48 , in spite of being affected by issues such as ribosomal copy number variability and preferential amplification. Taxa of lower ranks detected that lacked child nodes at order level were manually included into the respective dataset.
Statistical analysis. In every community composition dataset (plants, CLPPs of bacteria, AM spore morphotypes, amplicon sequences of 16S rRNA and ITS), two measures of diversity were calculated: species richness (S) and Shannon's diversity index (H') 49 . In the case of CLPPs of bacteria, the number of utilized substrates (i.e., the number of substrates with an absorbance value > 0.25) was used equivalent to S; H' was calculated considering absorbance values at each well as equivalent to species abundance. For amplicon sequence data, rarefied richness estimates, interpolating the expected richness at the lowest sample-specific sequencing depth, were used to compensate for variation in read depth across samples, using the R package vegan 50 . Pearson correlation analyses were carried out in order to compare these diversity indexes.
To compare grazed vs. non-grazed plots on all measured plant and soil parameters, linear mixed-effect models were used as implemented in the R package nlme (lme function 51 ), including grazing as fixed effect and site as random effect. In order to study the influence of grazing at the community composition level, partial canonical correspondence analyses (pCCA) were performed with compositions as response variables, and including grazing as an explanatory variable and site as a covariable. Likewise, in order to study the influence of other explanatory variables besides grazing, such as habitat and bedrock, redundancy analyses (RDA) were performed following forward selection of the subset of explanatory variables that significantly explained the variation in community composition and other response variables. The statistical significance of all canonical axes was also tested by means of non-parametric Monte Carlo permutation test. All multivariate analyses were done using Canoco 5 52 . Mantel-tests were performed to compare community composition between Bray-Curtis similarity matrices, using the program PAST 53 .