Data analysis
Bioinformatic processing . Bioinformatic processing of the
sequence data was conducted for microarthropods (springtails and mites;
COI), bacteria (16S), and fungi and protists (18S) according to standard
procedures (see supplementary methods S6).
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. 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 (Oksanen et al. 2017). 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.
Statistical procedures . First, the effect of cessation of grazing
on biotic and abiotic properties were explored using a series of mixed
effect models where biotic and abiotic factors were included as response
variables (including aboveground biomass, litter depth, soil
temperature, electrical conductivity, available inorganic nitrogen,
mineralisation rates, pH and bulk density) and grazing treatment as a
fixed predictor. Second, to explore whether different grazing treatments
affected plant species composition, a nonmetric multidimensional scaling
analysis was conducted, using the R-package vegan (Oksanenet al. 2017). Linear mixed models were constructed to test the
effect of grazing treatment on α- and β-diversity for all groups of soil
organisms and plants, and included grazing/grazing removal treatment as
a fixed effect. Linear mixed models were constructed using R-packages
lme4 (Bates et al. 2015) and lmerTest to investigate the drivers
of α- and β-diversity for all groups of soil organisms. We included the
following predictors as fixed predictors: grazing treatment, pH, soil
carbon (% dw), depth of organic layer (cm), litter biomass
(g.m-2), aboveground biomass
(g.m-2), cover of grasses & herbs, average soil bulk
density (g.dm-3), soil moisture
(g.dm-3), plant richness (number of species per 2x2m
quadrat), mineral N (mg.kg-1) and two-way interactions
between each of these variables and the grazing/grazing removal
treatment. Site was included as a random predictor throughout. Sites
with very deep organic layers (n=8) were excluded from this analysis, as
we were unable to estimate the carbon storage of these sites. All
predictors were checked for multicollinearity and correlating predictors
were excluded from the analysis. Assumptions on homogeneity of variances
and normality of residuals were met for each of the models. For each
species group, a full model, consisting of all fixed effects and their
interaction with grazing removal was constructed, which was subsequently
simplified using the step() function in the lmerTest package (Kuznetsova
2017); all variables with p < 0.1 were retained in the final
model.