There is widespread concern that cessation of grazing in historically grazed ecosystems is causing biotic homogenization and biodiversity loss. We used 12 montane grassland sites along an 800-km north-south gradient across the United Kingdom, to test whether cessation of grazing affects local ɑ- and β-diversity of belowground food webs. We show cessation of grazing leads to strongly decreased ɑ-diversity of most groups of soil microbes and fauna, particularly of relatively rare taxa. In contrast, the β-diversity varied between groups of soil organisms. While most soil microbial communities exhibited increased homogenization after cessation of grazing, we observed decreased homogenization for soil fauna after cessation of grazing. Overall, our results indicate that exclusion of domesticated herbivores from historically grazed montane grasslands has far-ranging negative consequences for diversity of belowground food webs. This underscores the importance of grazers for maintaining the diversity of belowground communities, which play a central role in ecosystem functioning.

Study sites. We selected 12 montane grassland sites across an 800 km north-south gradient of the United Kingdom (Figure 1). Sites were selected based on the following criteria: 1) grasslands had never received inorganic fertilisers or herbicides; 2) grazer exclusion plots had to be present for at least 10 years; 3) sites needed to be sufficiently far apart (>10 km) to be considered independent from each other; 4) the main form of historical management at the site is extensive grazing by sheep. Sites were typically grazed by purebred sheep at stocking densities of 1–2 ewes per hectare per year, although historical variation in grazing pressure across sites has resulted in mosaics of vegetation with patches of short and tall grass, interspersed with patches of dwarf-shrubs dominated by Calluna vulgaris, Vaccinium myrtillus, Erica tetralix and Erica cinerea. All sites were visited and sampled once, between 28th of April and 7th June 2015. Fenced grazing exclosures varied in size from 25 m²–10.66 ha and in age since cessation of grazing from 10 to 65 years (Table 1, Table S1 and Figure 2). At the sites, we sampled 4 paired 5 x 5 m plots where extensive grazing had been excluded by fencing and 4 adjacent grazed plots (Figure 2). At one location, Exmoor (site 11, Figure 2), we sampled three paired sites, and in the Peak District (site 9, Figure 2) we sampled six sites; therefore, the total number of plots was 98, half grazed and half ungrazed, from 12 distinct sites. The elevation from the sites varied between 300 m and 700 m asl (Table 1), and differed in underlying geology, climatic conditions, soil characteristics and dominant plant species (Table S1–4). In each plot, we assessed the vegetation composition and biomass (see supplementary information for details). Soil samples were collected to determine soil abiotic properties (bulk density, water content, carbon (C), nitrogen (N), phosphorus (P) content, pH and the potential nitrogen-mineralization. To calculate soil bulk density and water content (% of initial fresh weight), the bulk density rings were weighed, dried at 70 °C for 72 h, and weighed again. Soil carbon (C), nitrogen (N) and phosphorus (P) content of the surface organic horizon and subsurface mineral layer were determined using a subsample of the composite soil sample. Soil C and N content was determined using an automated elemental analysis (vario EL Cube, Elementar Analysesysteme GmbH). Total P (ortho-P) was determined using the ashing method (Kuo 1996). Briefly, 2 grams of air-dried soil from the composite soil sample was placed in a muffle furnace at 550 °C for 4 hrs. After cooling down, samples were transferred to a centrifuge tube in which 50ml 0.5M H₂SO₄ was added, which was placed in a horizontal shaker for 16hrs at 120 rpm. After centrifuging (10 min at 1500 rpm) ortho-P was determined colourimetrically using an auto-analyser (Seal AA3, SEAL Analytical, UK).  Soil extractable nitrogen and soil pH was determined using a KCl-extraction on a 10g subsample from the composite soil sample using the method described by Berendse et al., (1989). In short, pH was determined 15 minutes after adding 25 ml MilliQ to the 10g sample (1:2.5 w/w sludge). After that, 25 ml of 2M KCl was added and samples were filtered and analysed colourimetrically using an auto-analyser (Seal AA3, SEAL Analytical, UK). Potential nitrogen-mineralization was determined by storing a second 200 g sample for 2 weeks at 20ºC, from which the amount of available N was extracted from a 10g subsample using the same method.

Assessment of the composition of the belowground community

Nematode communities. Nematodes were extracted from 200 g of composite soil sample using the elutriator - cotton wool filter method. Nematode suspensions were concentrated and DNA was extracted by a lysis buffer including mammalian DNA as an external standard to monitor losses due to sample handling and DNA purification. DNA extracts were purified using a glass-fibre column-based procedure.  Purified DNA extracts were stored at -20ºC. To assess overall nematode biodiversity across all sites, 5 ml aliquots of all purified extracts were combined. The resulting mixture was analysed by qPCR using 72 nematode taxon-specific primer sets (Table S5). Based on the outcome of the overall biodiversity assessment, 30 nematode taxa were selected for qPCR-based quantification in each of the 98 samples. Two additional qPCR primer sets were used: one primer set was used to assess total nematode densities per sample, and a second primer set was used to quantify DNA levels of the external standard. Quantitative PCR reactions were executed and Ct-values were converted to nematode counts per 200 g soil. 

Microbial communities. Genomic DNA of bacteria, fungi, and protists was extracted from 1.5 grams of each composite soil sample using the Qiagen DNeasy PowerSoil 96-well extraction method. DNA was amplified in triplicate using primers specific to targeted regions within either the 16S or 18S rRNA gene (for prokaryotic and eukaryotic analyses, respectively). A portion of the 16S rRNA gene was amplified using the archaeal- and bacterial-specific primer set 515f/806r [37]. This 16S rRNA gene primer set is designed to amplify the V4–V5 region of both Archaea and Bacteria, has few biases against specific taxa and accurately represents phylogenetic and taxonomic assignment of sequences [38]. The 18S rRNA gene was amplified using the eukaryotic-specific primer set F1391 (50-GTACACCGCCCGTC-30) and REukBr (50-TGATCCTTCTGCAGGTTCACCTAC-30). The 18S rRNA gene primer set is designed to amplify the hypervariable V9-region of eukaryotes, with a focus on microbial eukaryotic lineages, including both protists and fungi. Amplicons were sequenced on two lanes of a 2 x 151 bp sequencing run on the Illumina HiSeq 2500 operating in Rapid Run Mode.

Microarthropod communities. To determine the community composition of microarthropods (mites: Acari, and springtails: Collembola), the large intact soil core was extracted using the Tullgren extractors at Lancaster University. Per plot, the batch of extracted microarthropods was collected in 96% ethanol and further processed to enable DNA-based identification. DNA extraction was performed using the DNeasy Blood & Tissue kit. DNA was amplified using the MiteMinBarF7 and MiteMinBarR4 primers, which target a ~200 bp fragment located within the cytochrome oxidase subunit 1 (COI) region, and were specifically designed to cover a wide diversity of microarthropods in NW-European grasslands. At Aarhus University (Roskilde, Denmark), amplicons were prepared for in-house paired-end sequencing on an Illumina MiSeq platform, using the Nextera XT indexing kit (Illumina, San Diego, CA, USA). The resulting amplicon libraries were purified using HighPrep™ PCR (Magbio Genomics Inc., Gaithersburg, Maryland, US) beads, quantified and equimolarly pooled, upon sequencing using the 250 bp paired end MiSeq version 2 reagent kit (Illumina, San Diego, CA, USA).

Diversity calculations. Changes in species richness resulting from cessation of grazing were calculated per site (e.g. Lake District) for the different groups using the formula: Δ ɑ = (ɑ[ungrazerd]-ɑ[grazed])/ɑ[ungrazed], which gives a response ratio for the site where grazing was excluded. We refrained from using the Shannon- and Simpson’s indices as primer choice affects the relative abundance of phylotypes. Positive values represent an increase in species richness at a given site, negative values represent a decrease. To calculate β-diversity for each taxonomic group, we calculated the dissimilarity among all grazed and among all ungrazed plots at a given site, following Ferrier et al. 2007 to obtain this metric, we calculated a Bray-Curtis dissimilarity matrix within site, using the vegan package. We then averaged the dissimilarities per site to obtain an estimate for β-diversity for both treatments at each of the 12 locations. For example, the β-diversity of the grazed treatment in the Yorkshire Dales (which had 4 independent plots) was based on the average dissimilarity in each soil organismal group for all pairwise combinations of grazed plots 1, 2, 3, and 4, which were subsequently averaged. We used the same procedure for the ungrazed plots. Changes in community dissimilarity between grazing treatments were calculated as follows: Δdissimilarity = (dissimilarity[ungrazed] - dissimilarity[grazed]) / dissimilarity[ungrazed]. 

To investigate the effect of cessation of grazing on relatively rare, common and widespread taxa, we performed a separate analysis. For this, we divided taxa into three groups: taxa that were present in less than 25% of all plots (n=98) were considered relatively rare; taxa present in 25 to 75% of all plots were considered relatively common; and taxa present in more than 75% of all plots were considered relatively widespread (hereafter referred to as rare, common and widespread). We then investigated how species richness within each of the organismal groups responded to cessation of grazing using a similar method as explained above for ɑ-diversity.