DOI: 10.5061/dryad.7h44j0zxr
title: Data for: Cessation of grazing causes biodiversity loss and homogenization of soil food webs

description: 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.
description: <strong>Study
 sites. </strong>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 (&gt;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 <em>Calluna
 vulgaris, Vaccinium myrtillus, Erica tetralix </em>and<em>
 Erica cinerea</em>. All sites were visited and sampled once, between
 28th of April and 7th June 2015. Fenced grazing exclosures varied in size
 from 25 m<sup>2</sup>–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<sub>2</sub>SO<sub>4</sub> 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 <strong>Nematode
 communities.</strong> 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. <br> <strong>Microbial communities.
 </strong>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. <strong>Microarthropod communities.
 </strong>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). <strong>Diversity
 calculations.</strong> 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.
description: # Data for: Cessation of grazing causes biodiversity loss and
 homogenization of soil food webs
 [https://doi.org/10.5061/dryad.7h44j0zxr](https://doi.org/10.5061/dryad.7h44j0zxr) Dear Reader, The 7 Excel files in this DRYAD repository constitute an addition to the data supplied in the supplemental material of this paper. Please note that most files consist of more than 1 sheet. Below we have described the content of these sheets Corresponding names of the sheets are denoted between brackets below (e.g. ("summary of values"). The file entitled "All\_diversity\_and\_abiotic\_data.xlsx" contains the total numbers for each of the measured plots for all the variables of interest for which the methods are described in the method section of the paper. * All column names of the overview sheet ("ALL\_DATA") are explained in a second sheet, sheet 2, entitled "Metadata". Info on the column names can be found in column B of the sheet "Metadata". The file entitled "ANALYSIS\_veg\_species\_&\_funct\_groups.xlsx" contains an overview of the species found in each of the plots, the changes in species diversity for plants and the changes in rarity. * Sheet 1: contains all the total cover for each plant species in each of the plots (Column A-D: Plot Characteristics, Column E: Plant name, Column F Total cover, Column G-Z: cover for each of the plant species recorded in each of the plots * Sheet 2: plant species divided in vascular plants (Columns B-CA) and other recorded species (Columns CC-CO) * Sheet 3: vascular plant species divided in rare (Columns B-BR) and widespread (Columns BU-CC). Summation of Rare species in Column BS and common species in Column CD. * Sheet 4: For each of the sites, a response ratio was calculated for the rare (Column C) and common species (Column F). The file entitled "Bact\_rarity\_analysis.xlsx", has the data and the analyses that were carried out to calculate changes in relative rarity in response to grazing. * Sheet 1: ("Bact OTU") contains the abundance of each bacterial OTU (rows 2-10337 contain all OTUs) in each of the plots (column B-CS). Rows are sorted based on commonness of OTUs (column CT). * Sheet 2: ("rare") Similar structure to sheet 1 but columns & rows transposed; contains only the rare species, present in &lt;25% of the plots. * Sheet 3: ("common") Similar structure to sheet 1 but columns & rows transposed; contains only the common species, present in 25-75% of the plots. * Sheet 4: ("widespread") Similar structure to sheet 1 but columns & rows transposed; contains, only the widespread species, present in &gt;75% of the plots. * Sheet 5: ("summary rare common widespread") summary of abundance classes (rare, common, widespread) of bacterial OTUs on sheet 2-4 * Sheet 6 ("response ratio"): summary of abundance of bacterial OTUs, taken from on sheets 2-4. Calculations on the response ratios for each of the rarity classes are in column D,F and H according to the description in the methods of the paper. In the files "Fung\_rarity\_analysis.xlsx" and "Protozoa\_rarity\_analysis.xlsx" one can find all the analyses that were carried out to calculate changes in its relative rarity in response to grazing. The structure of these files is identical to the file "Bact\_rarity\_analysis", please use the explanation above. The file entitled "Overall\_Species Richness\_analyses.xlsx" presents the changes in total diversity between grazed and ungrazed and response ratios for the groups of soil organisms are calculated. * Sheet 1("sp. rich. all groups") contains the species richness data for each of the groups of study. Unique code for each site, composed of the name of the site, whether or not it was grazed and the plot nr. Specifically, columns names can be explained as follows: * Column A-C: plot characteristics -name of side, and whether it was grazed or not, replicate plot nr. * Column D-J: richness of the various groups investigated in this study: * Fung\_OTUs: Number of fungal OTUs observed in this plot * ProtozoaOTUs: * Number of protozoa OTUsobserved in this plot * Bacteria\_OTUs: Number of bacterial OTUs observed in this plot * Mite\_family\_richness: Number of mite families observed in this plot * Collembola\_family\_richness: Number of collembola families observed in this plot * Nematode\_richness: Number of nematode taxa observed in this plot * Plant\_species\_richness: Number of plant species observed in this plot. * Sheet 2 ("analysis"): For each of the species groups, a response ratio was calculated for the rare (Column C) and common species (Column F). * Response ratios for species richness are calculated in this sheet, for each of the species groups, according to the method description in the paper.