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.