Data for: Cessation of grazing causes biodiversity loss and homogenization of soil food webs
Cite this dataset
Schrama, Maarten et al. (2023). Data for: Cessation of grazing causes biodiversity loss and homogenization of soil food webs [Dataset]. Dryad. https://doi.org/10.5061/dryad.7h44j0zxr
Abstract
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.
README: Data for: Cessation of grazing causes biodiversity loss and homogenization of soil food webs
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 <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 >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.
Methods
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 m2–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 H2SO4 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.
Funding
Dutch Research Council, Award: Rubicon grant 2014.08
Biotechnology and Biological Sciences Research Council, Award: BB/L026406/1