Divergent national-scale trends of microbial and animal biodiversity revealed across diverse temperate soil ecosystems.

Soil biota accounts for ~25% of global biodiversity and is vital to nutrient cycling and primary production. There is growing momentum to study total belowground biodiversity across large ecological scales to understand how habitat and soil properties shape belowground communities. Microbial and animal components of belowground communities follow divergent responses to soil properties and land use intensification; however, it is unclear whether this extends across heterogeneous ecosystems. Here, a national-scale metabarcoding analysis of 436 locations across 7 different temperate ecosystems shows that belowground animal and microbial (bacteria, archaea, fungi, and protists) richness follow divergent trends, whereas β-diversity does not. Animal richness is governed by intensive land use and unaffected by soil properties, while microbial richness was driven by environmental properties across land uses. Our findings demonstrate that established divergent patterns of belowground microbial and animal diversity are consistent across heterogeneous land uses and are detectable using a standardised metabarcoding approach.

Introduction

Soil n class="Species">biota, iclass="Chemical">ncludiclass="Chemical">ng bacteria, archaea, protists, fuclass="Chemical">ngi, aclass="Chemical">nd aclass="Chemical">nimals, uclass="Chemical">nderpiclass="Chemical">n globally importaclass="Chemical">nt ecosystem fuclass="Chemical">nctioclass="Chemical">ns. Fuclass="Chemical">ndameclass="Chemical">ntal fuclass="Chemical">nctioclass="Chemical">ns of soil commuclass="Chemical">nities iclass="Chemical">nclude class="Chemical">nutrieclass="Chemical">nt aclass="Chemical">nd hydrological cycliclass="Chemical">ng, decompositioclass="Chemical">n, pollutioclass="Chemical">n mitigatioclass="Chemical">n, aclass="Chemical">nd supporticlass="Chemical">ng terrestrial primary productioclass="Chemical">n, which are iclass="Chemical">nextricably liclass="Chemical">nked to global food security, climate regulatioclass="Chemical">n, aclass="Chemical">nd other ecosystem services[1,2]. Nevertheless, uclass="Chemical">ntil receclass="Chemical">ntly, characterisiclass="Chemical">ng soil biodiversity (popularly referred to as a ‘black box’) has beeclass="Chemical">n coclass="Chemical">nstraiclass="Chemical">ned by our iclass="Chemical">nability to ideclass="Chemical">ntify typically iclass="Chemical">ntractable levels of diversity usiclass="Chemical">ng either traditioclass="Chemical">nal or molecular approaches. High-throughput sequeclass="Chemical">nciclass="Chemical">ng has however resulted iclass="Chemical">n a step chaclass="Chemical">nge, facilitaticlass="Chemical">ng the characterisatioclass="Chemical">n of bacteria[3-7], archaea[6-8], fuclass="Chemical">ngi[9,10], protists[11-13], aclass="Chemical">nd aclass="Chemical">nimals[14] withiclass="Chemical">n the belowgrouclass="Chemical">nd biosphere. Iclass="Chemical">ncreasiclass="Chemical">ngly, efforts have beeclass="Chemical">n made to iclass="Chemical">nvestigate the total biodiversity of the soil biosphere across large ecological[15-17] aclass="Chemical">nd taxoclass="Chemical">nomic scales[15,16,18,19].

Understanding the response of the total soil biosphere to changes in land use and environmental drivers has become an important research focus in regional soil monitoring programmes[15,16,19] and in small-scale field[20,21] and mesocosm experiments[18,20]. Yet despite the move towards unified study of soil n class="Species">biota, fuclass="Chemical">ndameclass="Chemical">ntal challeclass="Chemical">nges of techclass="Chemical">nique aclass="Chemical">nd scale remaiclass="Chemical">n. Ofteclass="Chemical">n, such studies require the comparisoclass="Chemical">n of soil class="Chemical">n class="Species">biota metrics captured through both traditional and modern molecular techniques[15,19-21]. To our knowledge, relatively few studies have attempted to assess all components of belowground communities using a multi-marker metabarcoding approach[22].

There is mounting evidence that the microbial and animal fractions of soil communities may respond differentially to land use change. Microbial richness increases[15], whereas n class="Disease">richness of soil fauna declines iclass="Chemical">n respoclass="Chemical">nse to more iclass="Chemical">nteclass="Chemical">nse laclass="Chemical">nd use[15,23,24]. However, these ficlass="Chemical">ndiclass="Chemical">ngs come from relatively homogeclass="Chemical">nous laclass="Chemical">ndscapes, such as grasslaclass="Chemical">nds[15]. It is uclass="Chemical">nclear whether the differeclass="Chemical">ntial respoclass="Chemical">nses of soil microbes aclass="Chemical">nd fauclass="Chemical">na exteclass="Chemical">nd across heterogeclass="Chemical">neous laclass="Chemical">nd uses. For example, across heterogeclass="Chemical">neous laclass="Chemical">ndscapes of Wales, UK, α-diversity of mesofauclass="Chemical">na is both lowest iclass="Chemical">n agricultural aclass="Chemical">nd bog systems, which are the most- aclass="Chemical">nd least-iclass="Chemical">nteclass="Chemical">nsively maclass="Chemical">naged systems iclass="Chemical">n the couclass="Chemical">ntry, respectively[23]. Chaclass="Chemical">nges iclass="Chemical">n soil properties may further dictate decliclass="Chemical">nes of commoclass="Chemical">n soil fauclass="Chemical">na iclass="Chemical">n low-iclass="Chemical">nteclass="Chemical">nsity laclass="Chemical">nd uses. Therefore, it is critical to assess whether the positive effect of iclass="Chemical">ncreasiclass="Chemical">ng laclass="Chemical">nd use iclass="Chemical">nteclass="Chemical">nsity oclass="Chemical">n microbial richclass="Chemical">ness is coclass="Chemical">nsisteclass="Chemical">nt across regioclass="Chemical">ns made up of markedly diverse ecosystems aclass="Chemical">nd laclass="Chemical">nd uses. Similarly, the importaclass="Chemical">nce of iclass="Chemical">ndividual soil properties iclass="Chemical">n shapiclass="Chemical">ng belowgrouclass="Chemical">nd commuclass="Chemical">nities has also proveclass="Chemical">n difficult to diseclass="Chemical">ntaclass="Chemical">ngle. Maclass="Chemical">ny studies have democlass="Chemical">nstrated the coclass="Chemical">nsisteclass="Chemical">nt domiclass="Chemical">naclass="Chemical">nce of pH iclass="Chemical">n shapiclass="Chemical">ng belowgrouclass="Chemical">nd commuclass="Chemical">nity compositioclass="Chemical">n at class="Chemical">natioclass="Chemical">nal[23,25-28] aclass="Chemical">nd global scales[4,5,9,29]. However, climatic factors[9,30] aclass="Chemical">nd other soil properties, iclass="Chemical">ncludiclass="Chemical">ng class="Chemical">n class="Disease">organic matter, nitrogen (N) availability, and the carbon (C)-to-N ratio[9], are also recognised as important drivers of belowground community composition yet consistent trends remain elusive[30]. Therefore it is unclear whether the total soil biosphere responds to changes in land use and soil properties in the same manner across heterogeneous landscapes.

Here, we sought to assess whether divergent responses to land use and soil properties in the microbial and animal fractions of soil communities persist across heterogeneous systems at the national-scale using a standardised metabarcoding approach. We present a national-scale analysis of soil biodiversity across Wales, UK, from the micro-to-macro scale including all major groups of soil microbes in addition to animals, from 436 sites over 2 years across a diverse array of oceanic-temperate ecosystems, including grasslands, forests, bogs, and managed systems. Biotic metrics come from high-throughput sequencing of prokaryotic, fungal, microbial eukaryotic, and soil animal communities using 16S, ITS, and 18S rRNA marker genes; these are complemented by an extensive suite of co-located abiotic soil properties and vegetation cover data. Specifically, we investigate how richness and β-diversity of all major fractions of subterranean life respond to land use type and prevailing soil properties (e.g. n class="Disease">organic matter, pH, aclass="Chemical">nd N) to explore which liclass="Chemical">neages play a democlass="Chemical">nstrable role iclass="Chemical">n determiclass="Chemical">niclass="Chemical">ng belowgrouclass="Chemical">nd commuclass="Chemical">nity structures across large aclass="Chemical">nd complex ecological gradieclass="Chemical">nts. Our results democlass="Chemical">nstrate that across a gradieclass="Chemical">nt of heterogeclass="Chemical">neous laclass="Chemical">nd uses, richclass="Chemical">ness of soil aclass="Chemical">nimals is goverclass="Chemical">ned more by laclass="Chemical">nd use regime rather thaclass="Chemical">n iclass="Chemical">ntriclass="Chemical">nsic soil properties. Iclass="Chemical">n coclass="Chemical">ntrast, microbial richclass="Chemical">ness is driveclass="Chemical">n by soil properties aclass="Chemical">nd democlass="Chemical">nstrates a largely liclass="Chemical">near treclass="Chemical">nd of decreasiclass="Chemical">ng richclass="Chemical">ness aloclass="Chemical">ng a productivity gradieclass="Chemical">nt of laclass="Chemical">nd use based oclass="Chemical">n decreasiclass="Chemical">ng soil class="Chemical">nutrieclass="Chemical">nt availability.

Results

Sequencing results

Illumina sequencing and environmental data were collected from across Wales as part of the Glastir Monitoring and Evaluation Programme (n class="Chemical">GMEP)[31]. Sample sites were categorised iclass="Chemical">nto Aggregate Vegetatioclass="Chemical">n Classes (AVCs) based oclass="Chemical">n placlass="Chemical">nt species assessmeclass="Chemical">nts usiclass="Chemical">ng established criteria (see Supplemeclass="Chemical">ntary Note 1). Aclass="Chemical">n explaclass="Chemical">natioclass="Chemical">n of the compositioclass="Chemical">n of AVCs is described iclass="Chemical">n Supplemeclass="Chemical">ntary Table 1. Briefly, the 7 AVCs used iclass="Chemical">n the curreclass="Chemical">nt study were established by clustericlass="Chemical">ng samples based oclass="Chemical">n aclass="Chemical">n assessmeclass="Chemical">nt of vegetatioclass="Chemical">n data usiclass="Chemical">ng a detreclass="Chemical">nded correspoclass="Chemical">ndeclass="Chemical">nce aclass="Chemical">nalysis[32]. The ordiclass="Chemical">natioclass="Chemical">n of the detreclass="Chemical">nded correspoclass="Chemical">ndeclass="Chemical">nce aclass="Chemical">nalysis has showclass="Chemical">n that the laclass="Chemical">nd use categories follow a gradieclass="Chemical">nt of soil class="Chemical">nutrieclass="Chemical">nt coclass="Chemical">nteclass="Chemical">nt[32] from which soil productivity aclass="Chemical">nd maclass="Chemical">nagemeclass="Chemical">nt iclass="Chemical">nteclass="Chemical">nsity caclass="Chemical">n also be iclass="Chemical">nferred (see Supplemeclass="Chemical">ntary Note 1 aclass="Chemical">nd Supplemeclass="Chemical">ntary Table 1). The AVCs iclass="Chemical">n desceclass="Chemical">ndiclass="Chemical">ng order of productivity are crops/weeds, fertile grasslaclass="Chemical">nd, iclass="Chemical">nfertile grasslaclass="Chemical">nd, lowlaclass="Chemical">nd wood, uplaclass="Chemical">nd wood, moorlaclass="Chemical">nd grass-mosaic, aclass="Chemical">nd heath/bog.

In total, 29,690 bacterial and 156 archaeal operational taxonomic units (OTUs) were identified from 16S reads. Overall, the most abundant class was Alphaproteobacteria (Fig. 1a). Proportional abundances (OTU n/total × 100) of Acidobacteria increased in less-productive land use types from its lowest in crops/weeds to its highest in heath/bog AVCs. In contrast, abundances of Actinobacteria followed the exact opposite trend, as did Spartobacteria and Bacilli (Fig. 2a). For archaea, Nitrososphaeria was the most abundant class overall (Fig. 1d); however, the proportion of Thermoplasmata became dominant in less productive AVCs (Fig. 2d).

Fig. 1

Sankey diagrams of proportional abundances of OTUs from all samples for major soil biota groups. Arms denote proportions of OTUs at the class-level for a bacteria; b fungi; of major lineages of c protists; class-level for d archaea; and at the phylum-level for e animals. For information on how this figure was created, please see Supplementary Methods

Fig. 2

Proportionate abundances of OTUs for major soil biota groups within each Aggregate Vegetation Class. Land uses are ordered from most (crops/weeds) to least (heath/bog) using the same divisions as Fig. 1 for a bacteria; b fungi; c protists; d archaea; and e animals

Sankey diagrams of proportional abundances of OTUs from all samples for major soil n class="Species">biota groups. Arms declass="Chemical">note proportioclass="Chemical">ns of OTUs at the class-level for a bacteria; b fuclass="Chemical">ngi; of major liclass="Chemical">neages of c protists; class-level for d archaea; aclass="Chemical">nd at the phylum-level for e aclass="Chemical">nimals. For iclass="Chemical">nformatioclass="Chemical">n oclass="Chemical">n how this figure was created, please see Supplemeclass="Chemical">ntary Methods

Proportionate abundances of OTUs for major soil n class="Species">biota groups withiclass="Chemical">n each Aggregate Vegetatioclass="Chemical">n Class. Laclass="Chemical">nd uses are ordered from most (crops/weeds) to least (heath/bog) usiclass="Chemical">ng the same divisioclass="Chemical">ns as Fig. 1 for a bacteria; b fuclass="Chemical">ngi; c protists; d archaea; aclass="Chemical">nd e aclass="Chemical">nimals

There were 7582 OTUs recovered from ITS1 sequences. Agaricomycetes were the most abundant class of fungi overall. There was also a large proportion of Sordariomycetes (Fig. 1b). Proportionate abundances of Sordariomycetes and Agaricomycetes followed contrasting trends, with the dominance of the former replaced by the later in lower productivity AVCs (Fig. 2b).

In total, 8683 protist OTUs were recovered from the 18S reads. Chloroplastida (green n class="Species">algae) was by far the most abuclass="Chemical">ndaclass="Chemical">nt protist group, followed by Rhizaria, class="Chemical">n class="Species">Stramenopiles, and then Alveolates (Fig. 1c). Green algae, largely comprised of unidentified sequences (Supplementary Fig. 1a), were least abundant in crops/weed and heath/bog sites (Fig. 2c). Proportions of Rhizaria were relatively constant across AVCs (Fig. 2c) and entirely comprised of Cercozoa (Supplementary Fig. 1b). Among Stramenopiles, proportions of Ochrophyta were also largely consistent, while those of Oomycetes and Bicosoecida followed contrasting trends across the productivity gradient of AVCs, declining and increasing, respectively (Supplementary Fig. 1c). Ciliates were the most common Alveolates in most AVCs; however, the proportion of Apicomplexa was greater in the lowland wood and grassland AVCs (Supplementary Fig. 1d). The proportion of Amoebozoa was surprisingly low (Fig. 1c), potentially due to primer bias in our study when compared to other studies[12,15]. Across AVCs Tublulinea was consistently dominant among the Amoebozoa, though divergent trends in Gracilipodida and Discosea can be seen along the productivity/intensity gradient (Supplementary Fig. 1e).

In the animal dataset, 1138 OTUs were recovered. Nematode OTUs were the most abundant animal group across all samples (Fig. 1e). Annelids and arthropods followed opposing trends in proportionate abundance, increasing and decreasing respectively, across the productivity gradient. Proportions of Platyhelminthes and Tardigrades also increased in less-productive AVCs (Fig. 2e).

Effect of land use on belowground richness

We found significant differences in biodiversity trends across land use types. There was a marked shift along the productivity gradient of crops/weeds-to-heath/bog in all organismal groups, except animals (Fig. 3). Significant differences in the mean richness of bacterial OTUs were prominent (F6,264 = 78.47, p < 0.0001) following ANOVA. Bacterial richness decreased in AVCs across the productivity gradient with highest values in the most productive crops/weeds and grasslands and lowest in the low productivity land uses (i.e. moorland grass-mosaic, heath/bog) (Fig. 3a). The same trend was also observed in fungi (F6,248 = 48.98, p < 0.001; Fig. 3b), and protists (F6,249 = 59.86, p < 0.001; Fig. 3c). For individual pair-wise comparisons see Supplementary Note 4. Richness of archaeal OTUs had an opposing trend to that of other microbial groups. Archaeal OTU richness was significantly lower (F6,185 = 24.37, p < 0.001) in higher-productivity AVCs and highest in the least-productive land-use types (Fig. 3d). In the crops/weeds, AVC richness of archaeal OTUs was significantly lower than upland wood (p = 0.01), moorland grass-mosaic (p = 0.005), and heath/bog sites (p < 0.001) based on Tukey’s post hoc tests, with the remaining land uses displaying intermediate OTU richness values.

Fig. 3

Boxplots of OTU richness for each organismal group. Richness of a bacteria; b fungi; c protists; d archaea; e animals are plotted against Aggregate Vegetation Class ordered from most (crops/weeds) to least (heath/bog) productive. Boxes are bounded on the first and third quartiles; horizontal lines denote medians. Black dots are outliers beyond the whiskers, which denote 1.5× the interquartile range. Source data are provided as a Source Data file

Animal OTU richness did not follow the trends observed in microbial communities. Differences observed with ANOVA were significant (F6,244 = 6.25, p < 0.001) but plateaued after the grassland AVCs, as opposed to the sloped trend of microbial groups across the productivity gradient (Fig. 3e). Richness in the infertile grasslands was significantly greater than in crops/weeds (p = 0.008), heath/bog (p = 0.003), and upland wood (p = 0.02) based on Tukey’s post hoc tests. Richness was lowest in the most intensively management crops/weeds sites and was shown to be significantly lower than richness of lowland woods (p = 0.04) with Tukey’s test. Collectively, the results demonstrate a strong divergence between the richness of animal and microbial communities across all AVCs.

Relationships of richness between organismal groups

Bacterial richness from the total dataset was significantly correlated with all other organismal groups (Supplementary Table 2). Such relationships were positive between bacterial richness and richness of fungi, protists, and animals. Similarly, there was a positive relationship between protistan richness and both fungal and animal richness. However, n class="Disease">archaeal richness democlass="Chemical">nstrated sigclass="Chemical">nificaclass="Chemical">nt, but class="Chemical">negative correlatioclass="Chemical">ns with all orgaclass="Chemical">nisms except aclass="Chemical">nimals. Iclass="Chemical">ndeed aclass="Chemical">nimal richclass="Chemical">ness (measured by metabarcodiclass="Chemical">ng) was oclass="Chemical">nly sigclass="Chemical">nificaclass="Chemical">ntly correlated with aclass="Chemical">nimals (measured by taxoclass="Chemical">nomic assessmeclass="Chemical">nt; Table 1) aclass="Chemical">nd protists (Supplemeclass="Chemical">ntary Table 2).

Table 1

Results of partial least squares regressions for soil biota against soil properties for richness

Soil and environmental variables	Taxon
	Bacteria	Archaea	Fungi	Protists	Animals
Total Ca	1.14 (R2 = 0.44***)	1.21 (R2 = 0.13***)	0.44	1.3 (R2 = 0.35***)	0.9
Total Na	0.93	0.89	0.93	0.8	1.18
C:N ratiob	1.45 (R2 = 0.41***)	1.31 (R2 = 0.09***)	1.64 (R2 = 0.28***)	1.67 (R2 = 0.35***)	0.1
Total P (mg kg−1)b	0.35	0.59	0.7	0.85	0.67
Organic matter (% LOI)a	1.47 (R2 = 0.5***)	1.27 (R2 = 0.14***)	1.13 (R2 = 0.29***)	1.27 (R2 = 0.35***)	1.08
pH (CaCl2)	1.98 (R2 = 0.51***)	1.68 (R2 = 0.25***)	1.52 (R2 = 0.23***)	1.56 (R2 = 0.33***)	0.9
Soil water repellencya,c	1.31 (R2 = 0.2***)	0.9	1.23 (R2 = 0.13***)	0.93	0.98
Volumetric water content (m3 m3 −1)	0.36	1.33 (R2 = 0.13***)	0.6	0.41	0.4
Soil bound water (g water g dry soil−1)	1.25 (R2 = 0.41***)	0.83	1.08 (R2 = 0.26***)	1.23 (R2 = 0.31***)	0.63
Rock volume (mL)	0.25	0.61	0.64	0.27	1.3
Bulk density (g cm3 −1)	1.39 (R2 = 0.44***)	1.43 (R2 = 0.18***)	1.41 (R2 = 0.29***)	1.5 (R2 = 0.35***)	1.39
Clay content (%)d	0.85	1.19 (R2 = 0.1***)	0.84	1.14 (R2 = 0.09***)	0.05
Sand content (%)d	0.45	0.16	0.6	0.51	0.78
Elevation (m)	1.66 (R2 = 0.42***)	1.7 (R2 = 0.27***)	1.68 (R2 = 0.22***)	1.65 (R2 = 0.36***)	0.57
Mean annual precipitation (mL)	1.08 (R2 = 0.25***)	1.75 (R2 = 0.3***)	1.44 (R2 = 0.18***)	1.48 (R2 = 0.27***)	0.46
Temperature (°C)	0.51	0.5	0.56	0.58	0.35
Collembolae	0.34	0.06	0.41	0.17	1.14 (R2 = 0.03***)
Mitese	0.49	0.2	1.17 (R2 = 0.03***)	0.23	1.74 (R2 = 0.08***)
Total mesofaunae	0.44	0.1	1.03 (R2 = 0.01*)	0.15	1.71 (R2 = 0.08***)

Positive relationships are written in bold and negative relationships are written in italics

aLog10-transformation

bSquare-root-transformation

cSoil water repellency was derived from median water drop penetration times (s)

dAitchison’s log-ratio transformation

eLog10 plus 1 transformation

***p < 0.001; **0.001 > p < 0.01; *0.01 > p < 0.05, and blank indicates p > 0.05

Results of partial least squares regressions for soil n class="Species">biota agaiclass="Chemical">nst soil properties for richclass="Chemical">ness

cSoil n class="Chemical">water repelleclass="Chemical">ncy was derived from mediaclass="Chemical">n class="Chemical">n class="Chemical">water drop penetration times (s)

Relationships between richness and environmental variables

Partial least squares (n class="Disease">PLS) regressioclass="Chemical">ns democlass="Chemical">nstrated that the divergeclass="Chemical">nce observed betweeclass="Chemical">n aclass="Chemical">nimal aclass="Chemical">nd microbial commuclass="Chemical">nities may be due to the effects of soil properties. No soil properties were sigclass="Chemical">nificaclass="Chemical">ntly correlated with richclass="Chemical">ness of soil aclass="Chemical">nimal OTUs (Table 1). Coclass="Chemical">nversely, there were stroclass="Chemical">ng relatioclass="Chemical">nships betweeclass="Chemical">n microbial richclass="Chemical">ness aclass="Chemical">nd a raclass="Chemical">nge of soil properties. However, although microbes were iclass="Chemical">nflueclass="Chemical">nced by the same eclass="Chemical">nviroclass="Chemical">nmeclass="Chemical">ntal variables, there were disticlass="Chemical">nct patterclass="Chemical">ns withiclass="Chemical">n each group. For example, while pH was the best predictor of bacterial richclass="Chemical">ness, it was raclass="Chemical">nked as secoclass="Chemical">nd for fuclass="Chemical">ngi aclass="Chemical">nd protists aclass="Chemical">nd third for archaea. Bulk declass="Chemical">nsity aclass="Chemical">nd C:N ratio were also major drivers of richclass="Chemical">ness across all microbial groups. Elevatioclass="Chemical">n (here closely liclass="Chemical">nked with precipitatioclass="Chemical">n aclass="Chemical">nd class="Chemical">n class="Disease">organic matter content) was the most important environmental variable in relation to archaea and protist richness. Organic matter and bulk density were strong predictors of fungal OTU richness. All environmental properties that had positive relationships with OTU richness of bacteria, fungi, and protists had negative relationships with archaea.

Community structure (β-diversity) across land uses

Non-metric multidimensional scaling (NMDS) using Bray–Curtis distances showed consistent differences in β-diversity between AVCs across all organismal groups. Plots show tight clustering of the crops/weeds, fertile grassland, and infertile grassland AVCs, whereas the other AVCs form a more dispersed organismal assemblage (Fig. 4 for bacteria and Supplementary Figs. 2–5). Results of PERMANOVAs were significant across all groups and analyses of dispersion were also significant (Fig. 4 for bacteria and Supplementary Figs. 2–5) for all groups except for the dispersion of animals (F6,401 = 0.67, p = 0.68) owing to the wide range of sample numbers within each AVC (Supplementary Fig. 5). We also found that this clustering was present using constrained canonical analyses of principle components (CAP) ordinations for each organismal group (Supplementary Figs. 6–10).

Fig. 4

Plot of the non-metric dimensional scaling ordination (stress = 0.06) of bacterial community composition across GMEP sites. Samples are coloured by Aggregate Vegetation Class. Results of PERMANOVA (F6,427 = 30.76, p = 0.001) and dispersion of variances of groups (F6,427 = 10.97, p = 0.001) were significant

Plot of the non-metric dimensional scaling ordination (stress = 0.06) of bacterial community composition across n class="Chemical">GMEP sites. Samples are coloured by Aggregate Vegetatioclass="Chemical">n Class. Results of PERMANOVA (F6,427 = 30.76, p = 0.001) aclass="Chemical">nd dispersioclass="Chemical">n of variaclass="Chemical">nces of groups (F6,427 = 10.97, p = 0.001) were sigclass="Chemical">nificaclass="Chemical">nt

pH was the best predictor of β-diversity from linear fitting for all soil organisms (Table 2 and Supplementary Tables 3–6). The n class="Chemical">carbon-to-class="Chemical">n class="Chemical">nitrogen (C:N) ratio was the second most important variable in all major groups except animals. Mean C:N values were higher in the crops/weeds and grassland AVCs and lower in the remaining land use types (Supplementary Table 6). Mean pH values and C:N ratios (Supplementary Table 6) reflect the distribution of points in NMDS plots, with tight groupings observed in the crops/weeds and grasslands AVCs and increasingly more spread out groupings in all other AVCs as pH values decreased and became more varied (Fig. 4 for bacteria and Supplementary Figs. 2–5). Across all groups, all or nearly all variables were significant following linear fitting; however, most were only weakly correlated with β-diversity values. Other important variables varied in their ranked importance, including elevation, mean annual precipitation, organic matter content, total C, bulk density, volumetric water content, and clay content of soil (Table 2 and Supplementary Tables 3–6). The results of linear model fitting for CAP ordinations, though not identical (Supplementary Tables 7–11), were highly related to those of the NMDS ordinations (Supplementary Fig. 11).

Table 2

Summary of relationships amongst environmental factors and bacteria communities

Soil and environmental variables	R 2	Correlation
		Axis 1	Axis 2	Axis 3
pH (CaCl2)	0.71***	−	−	+
C:N ratioa	0.52***	+	−	+
Volumetric water content (m3 m3 −1)	0.49***	+	−	+
Bulk density (g cm3 −1)	0.47***	−	+	−
Organic matter (% LOI)b	0.46***	+	−	+
Elevation (m)	0.45***	+	−	−
Mean annual precipitation (mL)	0.43***	+	−	−
Total Cb	0.39***	+	−	+
Clay content (%)c	0.33***	−	+	−
Soil bound water (g water g dry soil−1)	0.31***	+	−	+
Soil water repellencyb,d	0.27***	+	−	−
Total N (%)b	0.26***	+	−	+
Sand content (%)c	0.21***	+	+	+
Collembolae	0.09***	−	+	−
Mitese	0.06***	+	+	−
Total P (mg kg−1)a	0.06***	−	−	−
Total mesofaunae	0.06***	+	+	−
Rock volume (mL)	0.05**	−	+	−
Temperature (°C)	0.03*	+	+	−

+/− signify the direction of association between each variable and respective NMDS axes

aSquare-root-transformation

bLog10-transformation

cAitchison’s log-ratio transformation

dSoil water repellency was derived from median water drop penetration times (s)

eLog10 plus 1 transformation

***p < 0.001; **0.001 > p < 0.01; *0.01 > p < 0.05, and blank indicates p > 0.05

n class="Chemical">dSoil water repelleclass="Chemical">ncy was derived from mediaclass="Chemical">n class="Chemical">n class="Chemical">water drop penetration times (s)

Discussion

High-throughput sequencing of the biosphere amongst heterogeneous soils revealed both expected and novel relationships between soil organisms and environmental drivers. The richness of microbes and animals had notable contrasting trends across land use types. The richness of microbial communities was strongly influenced by both land use and environmental variables, especially pH, C:N ratio, elevation, n class="Disease">organic matter, aclass="Chemical">nd aclass="Chemical">nclass="Chemical">nual precipitatioclass="Chemical">n. Coclass="Chemical">nversely, we fouclass="Chemical">nd class="Chemical">no sigclass="Chemical">nificaclass="Chemical">nt associatioclass="Chemical">ns betweeclass="Chemical">n measured eclass="Chemical">nviroclass="Chemical">nmeclass="Chemical">ntal variables aclass="Chemical">nd aclass="Chemical">nimal richclass="Chemical">ness, which was class="Chemical">negatively impacted by higher iclass="Chemical">nteclass="Chemical">nsity laclass="Chemical">nd use, suggesticlass="Chemical">ng that richclass="Chemical">ness patterclass="Chemical">ns of microbial aclass="Chemical">nd macrobial life fractioclass="Chemical">ns adhere to differeclass="Chemical">nt ecological determiclass="Chemical">naclass="Chemical">nts. For β-diversity, pH was by far the most importaclass="Chemical">nt eclass="Chemical">nviroclass="Chemical">nmeclass="Chemical">ntal variable iclass="Chemical">n shapiclass="Chemical">ng commuclass="Chemical">nity compositioclass="Chemical">n of all orgaclass="Chemical">nismal groups, yet other drivers were attributable for iclass="Chemical">nflueclass="Chemical">nciclass="Chemical">ng patterclass="Chemical">ns of α-diversity.

Our findings demonstrate that diverging trends between soil microbes and fauna extend across distinct, heterogeneous land uses. Furthermore, we build on the work of Gossner et al.[15] by demonstrating that microbial richness, with the exception of archaea, increases with greater land use intensity across heterogeneous ecosystems at the national-scale. The divergence between microbes and animals at this scale is supported by previous findings from French soils[17,25]. Across France, bacterial richness[17] and biomass[25] were strongly linked to belowground environmental properties but largely unaffected by aboveground climatic variables, which commonly influence animal and plant biogeography[25,30]. Our findings show that richness of fungi and protists also follow this trend—whereas archaea follow an opposing trend to all other groups.

There are several mechanisms that may explain the relationship between higher microbial richness and intensifying anthropogenic disturbance. One explanation is that consistent nutrient inputs from fertilisers and disturbance under tillage stimulate high α-diversity in these areas[16]. Indeed higher α-diversity has been observed in cropping systems than in forest or grassland sites for both bacteria[16,17] and fungi[16]. Interestingly, high microbial richness in more productive land use types (e.g. arable) may illustrate the intermediate disturbance hypothesis (IDH) within soil ecosystems. Under the IDH, as outlined by Connell[33], diversity reaches its highest levels where succession has been interrupted by intermittent disturbance events. In our sites, microbial richness was highest in AVCs concurrent to disturbances (augmented by nutrient inputs) from agricultural interventions such as fertilisation, tilling, clearing, and the cultivation of livestock. However, it is also possible that the high diversity observed in the grassland and especially in agricultural land uses stems from organisms that have entered a dormant state after disturbance-induced changes to their environment[13,34]. Disturbance pressures can also lead to high bacterial diversity through the reduction in dominant OTUs, which are replaced by a wide range of weaker competitors. It has been demonstrated that α-bacterial diversity is greater in the phyllosphere of ivy in urban habitats associated with more anthropogenic stressors than in less disturbed sites[35]. Our findings suggest that the phenomenon of greater species richness resulting from the addition of nutrients and non-equilibrium dynamics induced by disturbance may extend to across all microbial groups, with the possible exception of archaea.

Richness of all microbial groups, except archaea, followed the land use productivity/management intensity gradient[32] with higher richness in the highly productive and more disturbed grasslands and arable sites and lower richness in the least productive, relatively undisturbed upland heath/bog sites. Changes within bacterial and fungal communities reflected expected within-community changes following the shift in soil nutrient quality across land uses. Actinobacteria[36] and Sordariomycetes[37] are known to dominate bacterial and fungal communities in high productivity grasslands as witnessed here. In contrast, Acidobacteria increased in proportion in low productivity, highly acidic AVCs as expected based on previous studies from the UK[27] and across the globe[7]. Likewise, the greater proportion of Agaricomycetes OTUs in low productivity AVCs is intuitive as many Agaricomycete fungi are common in bogs and related low-productivity habitats across Wales[38].

Protists have been chronically overlooked in European soil monitoring programmes (but see ref. [28]), as extracting trends of protist diversity across land uses is difficult. For example, Gossner et al.[15] were not able to show changes in richness across all protists with land use intensification. We demonstrate that protistan richness follows the trends of bacteria and fungi across land uses, with the highest richness levels in arable land. As with other microbes, there is evidence of increased protist richness at the mesocosm[39] and field[40] level, in response to fertiliser addition. Furthermore, in German grassland soils, protist richness has been shown to increase with land use intensity[41]. Our results show that an association between intensification and protistan richness extends across the national-scale over multiple land uses.

Unlike other microbes, n class="Disease">archaeal richness was greatest iclass="Chemical">n low productivity AVCs aclass="Chemical">nd lowest iclass="Chemical">n highly productive sites (Fig. 3d). Furthermore, our uclass="Chemical">nderstaclass="Chemical">ndiclass="Chemical">ng of the exteclass="Chemical">nt of soil archaeal diversity aclass="Chemical">nd its fuclass="Chemical">nctioclass="Chemical">nal capabilities is coclass="Chemical">nticlass="Chemical">nually iclass="Chemical">ncreasiclass="Chemical">ng[6-8]. Receclass="Chemical">nt research has revealed maclass="Chemical">ny liclass="Chemical">neages of Thaumarchaeota are crucial liclass="Chemical">nks iclass="Chemical">n the N cycle aclass="Chemical">nd methaclass="Chemical">nogeclass="Chemical">nesis iclass="Chemical">n soils[7,8]. class="Chemical">n class="Disease">Archaeal richness was highest in the moorland grass-mosaic and heath/bog AVCs, likely due to the specialised nature of acidophilic lineages. In particular, the Thaumarchaeota[42] and Thermoplasmata[43] are known to proliferate (Fig. 2d) under reduced competition from bacteria.

Animal richness did not change linearly with land use and was not strongly influenced by environmental variables. Our molecular analysis of soil eDNA supports recent findings by George et al.[23] based on morphological assessments of coincident soil mesofauna. Both the present work and George et al.[23] demonstrated that animal richness and abundance were lowest in land uses associated with more intensive management. Animal richness peaked in infertile grasslands and was lowest in crops/weeds sites (Fig. 3e). Agricultural disturbance negatively affects soil n class="Disease">faunal richness aclass="Chemical">nd diversity across large geographic scales[14,23,24]. However, iclass="Chemical">n the low-productivity laclass="Chemical">nd uses, although proportioclass="Chemical">nal abuclass="Chemical">ndaclass="Chemical">nces of arthropod taxa decliclass="Chemical">ned similarly to the ficlass="Chemical">ndiclass="Chemical">ngs of George et al.[23], overall richclass="Chemical">ness was class="Chemical">not as stroclass="Chemical">ngly affected due to aclass="Chemical">n iclass="Chemical">ncrease iclass="Chemical">n fractioclass="Chemical">ns of Aclass="Chemical">nclass="Chemical">nelids, Platyhelmiclass="Chemical">nthes, aclass="Chemical">nd Tardigrades. Such aclass="Chemical">n iclass="Chemical">ncrease iclass="Chemical">n the peat-rich, low-disturbaclass="Chemical">nce, higher elevatioclass="Chemical">n sites is rather iclass="Chemical">ntuitive siclass="Chemical">nce Aclass="Chemical">nclass="Chemical">nelids, Platyhelmiclass="Chemical">nthes, aclass="Chemical">nd Tardigrades are susceptible to desiccatioclass="Chemical">n aclass="Chemical">nd require moist habitats to be active compoclass="Chemical">neclass="Chemical">nts of the soil commuclass="Chemical">nity[44,45]. As soil aclass="Chemical">nimals still exhibited expected lower diversity treclass="Chemical">nds iclass="Chemical">n more iclass="Chemical">nteclass="Chemical">nsively maclass="Chemical">naged laclass="Chemical">nd uses[15,23,24], there are further opportuclass="Chemical">nities for research iclass="Chemical">nto uclass="Chemical">nderstaclass="Chemical">ndiclass="Chemical">ng the mechaclass="Chemical">nisms uclass="Chemical">nderlyiclass="Chemical">ng the divergeclass="Chemical">nt richclass="Chemical">ness treclass="Chemical">nds betweeclass="Chemical">n microscopic aclass="Chemical">nimals aclass="Chemical">nd the rest of soil commuclass="Chemical">nities.

Soil pH, as evidenced by ordination results, was the most important environmental variable in our study for β-diversity and in most cases richness as has been previously observed across the UK[27,28] and at larger national[25,26] and continental scales[4-6]. pH has been implicated with driving n class="Disease">richness of soil Archaea[42,43] aclass="Chemical">nd is the most importaclass="Chemical">nt driver of protist commuclass="Chemical">nities iclass="Chemical">n the UK[28]. However, pH oclass="Chemical">nly plays a margiclass="Chemical">nal role iclass="Chemical">n shapiclass="Chemical">ng soil protist commuclass="Chemical">nities globally[11]. Likewise, pH is a poor predictor of global fuclass="Chemical">ngal biogeography, yet is a good predictor of ectomycorrhizal fuclass="Chemical">ngal richclass="Chemical">ness[9], which may coclass="Chemical">ntribute to the Agaricomycetes OTUs observed iclass="Chemical">n the preseclass="Chemical">nt study. Nevertheless, it is importaclass="Chemical">nt to ackclass="Chemical">nowledge the iclass="Chemical">ncoclass="Chemical">nsisteclass="Chemical">nt class="Chemical">nature of correlatioclass="Chemical">ns betweeclass="Chemical">n microbial biodiversity aclass="Chemical">nd pH, poteclass="Chemical">ntially due to variatioclass="Chemical">ns iclass="Chemical">n soil properties occurriclass="Chemical">ng at scales that do class="Chemical">not aligclass="Chemical">n with large-scale soil surveys[30].

We also observed a strong effect of C:N ratio in determining richness of microbes and β-diversity of all organismal groups, as has been observed in bacterial[27] and protistan[28] β-diversity across Britain and some fungi globally[9]. Yet C:N ratio is often co-correlated with other soil properties including bulk density, total C, n class="Disease">organic matter, elevatioclass="Chemical">n, aclass="Chemical">nd meaclass="Chemical">n aclass="Chemical">nclass="Chemical">nual precipitatioclass="Chemical">n. Diseclass="Chemical">ntaclass="Chemical">ngliclass="Chemical">ng such related variables is difficult; despite usiclass="Chemical">ng class="Chemical">n class="Disease">PLS analyses[46] we could not disentangle co-correlated soil properties. For example, AVCs such as moorland grass-mosaic and heath/bog generally had higher elevation, mean annual precipitation, C:N ratio, and both total C and N (Supplementary Table 12) owing to their less-disturbed, upland location, and often peat-rich soils. Higher C:N ratios are indicative of lower-quality soils[47] and have historically been associated with a shift in microbial biomass from bacterial to fungal dominance[48]. Our results suggest that, with the exception of archaea, microbial richness is equally susceptible to the effect of soil quality degradation. According to our results, archaea, on the contrary, appear to be well adapted to habitats with lower nutrient quality.

We observed strong relationships between soil properties and microbial, but not animal richness. We suspect this is due to the direct effects of soil properties on microbes. For example, shifts in pH towards either a more alkaline or acidic condition inhibit the ability of most non-specialised bacteria to uptake nutrients from their environment[26]. In addition the quality of soil nutrients, as discussed previously, was likely a strong determinant of available nutrient resources and therefore total richness of microbes. We also found strong relationships between soil properties and β-diversity and across all organismal groups. These relationships between Bray–Curtis dissimilarities and soil properties demonstrate that more dissimilar belowground communities correlate positively with indicators of better quality soils across the breadth of soil n class="Species">biota (Supplemeclass="Chemical">ntary Table 6). However, associatioclass="Chemical">ns betweeclass="Chemical">n class="Chemical">nutrieclass="Chemical">nt quality aclass="Chemical">nd aclass="Chemical">nimal commuclass="Chemical">nity compositioclass="Chemical">n are likely the result of class="Chemical">nutrieclass="Chemical">nts iclass="Chemical">nflueclass="Chemical">nciclass="Chemical">ng the compositioclass="Chemical">n of the abovegrouclass="Chemical">nd placlass="Chemical">nt commuclass="Chemical">nity[49] rather thaclass="Chemical">n direct iclass="Chemical">nteractioclass="Chemical">ns with aclass="Chemical">nimals. Furthermore, aclass="Chemical">nimals are more vagile thaclass="Chemical">n microbes aclass="Chemical">nd caclass="Chemical">n actively seek out microhabitats with better resources[50], limiticlass="Chemical">ng the direct impact of soil properties oclass="Chemical">n aclass="Chemical">nimal richclass="Chemical">ness.

Using an extensive soil sampling programme and metabarcoding, we present perhaps the most comprehensive assessment of the belowground diversity in Europe. Despite uncertainties on the ability of environmental DNA methods using small soil volumes to accurately characterise communities of larger organisms[51], we were still able to detect key differences in larger organisms (i.e. animals) across land uses. Our results highlight the complexity of belowground ecology by demonstrating a divergence of patterns of richness between soil fauna and microorganisms at a national-level. We show that microbial richness is strongly influenced by soil properties in a near-uniform manner, whereas animal richness is not. Rather, animal richness is likely driven by changes in aboveground communities that stem from intensive land use management, while microbial richness was affected by soil properties in addition to land use. A particularly interesting outcome of our analyses is the near-uniform trend of declining microbial richness along a gradient of decreasing land use productivity/management intensity. The data therefore suggest that soil properties strongly affect bacteria, fungi, and protists in a similar manner, whereby n class="Disease">richness decreases with soil quality; whereas archaea showed aclass="Chemical">n opposiclass="Chemical">ng treclass="Chemical">nd with iclass="Chemical">ncreasiclass="Chemical">ng richclass="Chemical">ness as productivity decliclass="Chemical">ned. The richclass="Chemical">ness of aclass="Chemical">nimal OTUs, oclass="Chemical">n the coclass="Chemical">ntrary, was class="Chemical">not affected by soil properties although β-diversity was. Although ofteclass="Chemical">n coclass="Chemical">nsidered as ecological ‘black boxes’, soils coclass="Chemical">nticlass="Chemical">nue to provide uclass="Chemical">nique aclass="Chemical">nd cohereclass="Chemical">nt iclass="Chemical">nsights iclass="Chemical">nto the differeclass="Chemical">nces betweeclass="Chemical">n iclass="Chemical">ntercoclass="Chemical">nclass="Chemical">nected microbial aclass="Chemical">nd macrobial assemblages. Our ficlass="Chemical">ndiclass="Chemical">ngs also highlight the importaclass="Chemical">nce of the dyclass="Chemical">namics betweeclass="Chemical">n biotic aclass="Chemical">nd abiotic processes that drive the orgaclass="Chemical">nisatioclass="Chemical">n of belowgrouclass="Chemical">nd biological diversity.

Methods

Sampling

Soil samples were collected between late spring and early autumn in 2013 and 2014 as part of n class="Chemical">GMEP (Supplemeclass="Chemical">ntary Note 2), established to moclass="Chemical">nitor the Welsh Goverclass="Chemical">nmeclass="Chemical">nt’s agri-eclass="Chemical">nviroclass="Chemical">nmeclass="Chemical">nt scheme, Glastir. The scheme covered aclass="Chemical">n area of 3263 km2 with 4911 laclass="Chemical">ndowclass="Chemical">ners[31]. Briefly, surveyors collected samples from raclass="Chemical">ndomly selected 1 km2 squares with up to 3 locatioclass="Chemical">ns withiclass="Chemical">n squares, followiclass="Chemical">ng protocols established by the UK Couclass="Chemical">ntryside Survey[52]. As described previously, habitat withiclass="Chemical">n plots was classified usiclass="Chemical">ng placlass="Chemical">nt species assessmeclass="Chemical">nts iclass="Chemical">nto oclass="Chemical">ne of seveclass="Chemical">n AVCs[32]: crops/weeds (class="Chemical">n = 9), fertile grasslaclass="Chemical">nd (class="Chemical">n = 98), iclass="Chemical">nfertile grasslaclass="Chemical">nd (class="Chemical">n = 162), lowlaclass="Chemical">nd wood (class="Chemical">n = 17), uplaclass="Chemical">nd wood (class="Chemical">n = 44), moorlaclass="Chemical">nd-grass mosaic (class="Chemical">n = 54), aclass="Chemical">nd heath/bog (class="Chemical">n = 52) (Supplemeclass="Chemical">ntary Note 1; Supplemeclass="Chemical">ntary Table 1). Soil type was derived from the Natioclass="Chemical">nal Soil Map[53] (Supplemeclass="Chemical">ntary Note 3; Supplemeclass="Chemical">ntary Table 13). class="Chemical">n class="Disease">Organic matter content was classified by loss-on-ignition (LOI) following the protocols of the 2007 Countryside Survey[51].

A total of 436 cores were collected from 1 km2 squares, with up to 3 samples coming from an individual square based on a randomised sampling design. Cores were transported to the Centre for Ecology and Hydrology, Bangor, UK, and stored at −80 °C until DNA extraction. Soil physical and chemical properties were taken from 4 cm diameter by 15 cm deep cores co-located with the high-throughput sequencing cores. These included total C (%), N (%), P (mg kg−1), n class="Disease">organic matter (% LOI), pH (measured iclass="Chemical">n 0.01 M class="Chemical">n class="Chemical">CaCl2), mean soil water repellency (median water drop penetration time in seconds), bulk density (g cm3 −1), volume of rocks (cm3), soil bound water (g water g dry soil−1), volumetric water content (m3 m3 −1), as well as clay and sand content (%) of soil. Abundances of mesofauna collected as part of GMEP were taken from George et al.[23] and geographic data including grid eastings, northings, and elevation were also included in our analyses. For complete details on chemical analyses, see Emmett et al.[51]. Temperature (°C) and mean annual precipitation (mL) were extracted from the Climate Hydrology and Ecology research Support System dataset[54]. Mean values for each variable were recorded for each AVC (Supplementary Table 12) and soil properties were normalised where appropriate.

Soil texture data were measured by laser granulometry with a LS320 13 analyser (Beckman-Coulter). We subsampled approximately 0.5 g of soil taken from 15 cm cores by manual quartering and removed organic C using n class="Chemical">H2O2 aclass="Chemical">nd theclass="Chemical">n traclass="Chemical">nsferred the sample iclass="Chemical">nto 250 mL bottles, added 5 mL of 5% Calgoclass="Chemical">n® aclass="Chemical">nd shook overclass="Chemical">night at 240 rpm. Bottles were emptied maclass="Chemical">nually iclass="Chemical">nto the laser diffractioclass="Chemical">n iclass="Chemical">nstrumeclass="Chemical">nt for measuriclass="Chemical">ng particle size distributioclass="Chemical">n. Full Mie theory was used to obtaiclass="Chemical">n a particle size distributioclass="Chemical">n from the raw measuremeclass="Chemical">nt data, with the real refractive iclass="Chemical">ndex set to 1.55 aclass="Chemical">nd the absorptioclass="Chemical">n coefficieclass="Chemical">nt at 0.1 as iclass="Chemical">n Özer et al.[55]. The cut-off poiclass="Chemical">nts for clay, silt, aclass="Chemical">nd saclass="Chemical">nd were 2.2, 63, aclass="Chemical">nd 2000 μm, respectively. Clay aclass="Chemical">nd saclass="Chemical">nd perceclass="Chemical">ntages were selected for subsequeclass="Chemical">nt aclass="Chemical">nalyses aclass="Chemical">nd class="Chemical">normalised usiclass="Chemical">ng Aitchisoclass="Chemical">n’s log-ratio traclass="Chemical">nsformatioclass="Chemical">n.

DNA extraction

Soils were homogenised by passing through a sterilised 2 mm n class="Chemical">stainless steel sieve. Sieves were sterilised betweeclass="Chemical">n samples by riclass="Chemical">nsiclass="Chemical">ng uclass="Chemical">nder the tap class="Chemical">n class="Chemical">water using high flow, applying Vircon laboratory disinfectant and UV-treating each side for 5 min DNA was extracted by mechanical lysis and the homogenisation step performed in triplicate from 0.25 g of soil per sample using a PowerLyzer PowerSoil DNA Isolation Kit (MO-BIO). Pre-treatment with 750 μL of 1 M CaCO3 following Sagova-Mareckova et al.[56] was performed as it was shown to improve PCR performances, especially for acidic soils. Extracted DNA was stored at −20 °C until amplicon library preparation began. To check for contamination in sieves 3 negative control DNA extractions were completed and an additional 2 negative control kit extractions were performed using the same technique but without the CaCO3 solution.

Primer selection and PCR protocols for library preparation

Amplicon libraries were created using primers for rRNA marker genes, specifically for the V4 region of the 16S rDNA gene targeting bacteria and archaea (515F/806R)[57], ITS1 targeting fungi (ITS5/5.8S_fungi)[58], and the V4 region of the 18S rDNA gene (TAReuk454FWD1/TAReukREV3)[59] targeting a wide range of, but not all, eukaryotic organisms. We used a two-step PCR following protocols devised in conjunction with the Liverpool Centre for Genome Research. Amplification of amplicon libraries was run in triplicate on DNA Engine Tetrad® 2 Peltier Thermal Cycler (n class="Disease">BIO-RAD Laboratories) aclass="Chemical">nd thermocycliclass="Chemical">ng parameters for each PCR started with 98 °C for 30 s aclass="Chemical">nd termiclass="Chemical">nated with 72 °C for 10 miclass="Chemical">n for ficlass="Chemical">nal exteclass="Chemical">nsioclass="Chemical">n aclass="Chemical">nd held at 4 °C for a ficlass="Chemical">nal 10 miclass="Chemical">n For the 16S locus, first-rouclass="Chemical">nd PCR amplificatioclass="Chemical">n followed 10 cycles of 98 °C for 10 s; 50 °C for 30 s; 72 °C for 30 s. For ITS1, there were 15 cycles of 98 °C for 10 s; 58 °C for 30 s; 72 °C for 30 s. For 18S there were 15 cycles at 98 °C for 10 s; 50 °C for 30 s; 72 °C for 30 s. Twelve μL of each first-rouclass="Chemical">nd PCR product was mixed with 0.1 μL of exoclass="Chemical">nuclease I, 0.2 μL thermoseclass="Chemical">nsitive alkaliclass="Chemical">ne phosphatase, aclass="Chemical">nd 0.7 μL of class="Chemical">n class="Chemical">water and cleaned in the thermocycler with a programme of 37 °C for 15 min and 74 °C for 15 min and held at 4 °C. Addition of Illumina Nextera XT 384-way indexing primers to the cleaned first-round PCR products were amplified following a single protocol which started with initial denaturation at 98 °C for 3 min; 15 cycles of 95 °C for 30 s; 55 °C for 30 s; 72 °C for 30 s; final extension at 72 °C for 5 min and held at 4 °C. Twenty-five μL of second-round PCR products were purified with an equal amount of AMPure XP beads (Beckman Coulter). Library preparation for 2013 samples was conducted at Bangor University. Illumina sequencing for both years and library preparation for 2014 samples were conducted at the Liverpool Centre for Genome Research.

Bioinformatics

Bioinformatics analyses were performed on the Supercomputing Wales cluster. A total of 130,219,260, 104,276,828, and 98,999,009 raw reads were recovered from the 16S, ITS1, and 18S sequences, respectively. Illumina adapters were trimmed from sequences using Cutadapt[60] with 10% level mismatch for removal. Sequences were then de-multiplexed, filtered, quality-checked, and clustered using a combination of USEARCH v. 7.0[61] and VSEARCH v. 2.3.2[62]. Open-reference clustering (97% sequence similarity) of OTUs was performed using VSEARCH; all other steps were conducted with USEARCH. Sequences with a maximum error greater than 1 and shorter than 200 bp were removed following the merging of forward and reverse reads for 16S and ITS1 sequences. A cut-off of 250 bp was used for 18S sequences, according to higher quality scores. There were 15,202,313 (16S), 7,242,508 (ITS1), and 9,163,754 (18S) cleaned reads left at the end of these steps. Sequences were sorted and those that only appeared once in the dataset were removed. Briefly, filtered sequences were matched first against a number of different reference databases: Greengenes 13.8[63], UNITE 7.2[64], and n class="Disease">SILVA 128[65] for 16S, ITS1, 18S, respectively. Teclass="Chemical">n perceclass="Chemical">nt of sequeclass="Chemical">nces that failed to match were clustered de class="Chemical">novo aclass="Chemical">nd used as a class="Chemical">new refereclass="Chemical">nce database for failed sequeclass="Chemical">nces. Sequeclass="Chemical">nces that failed to match with the de class="Chemical">novo database were subsequeclass="Chemical">ntly also clustered de class="Chemical">novo. All clusters were collated aclass="Chemical">nd chimeras were removed usiclass="Chemical">ng the uchime_ref commaclass="Chemical">nd iclass="Chemical">n VSEARCH.

Chimera-free clusters and taxonomy assignment were used to create an OTU table with QIIME v. 1.9.1[66] using RDP[67] methodology with the GreenGenes database for 16S and UNITE database for ITS1 data. Taxonomy was assigned to the 18S OTU table using BLAST[68] against the n class="Disease">SILVA database aclass="Chemical">nd OTUs appeariclass="Chemical">ng oclass="Chemical">nly oclass="Chemical">nce or iclass="Chemical">n oclass="Chemical">nly 1 sample were removed from each OTU table.

Newick trees were constructed for the 16S and 18S tables using 80% identity thresholds. The trees were combined with their respective OTU tables as part of analyses using the R package n class="Chemical">phyloseq[69], removiclass="Chemical">ng OTUs that did class="Chemical">not appear iclass="Chemical">n both the tree aclass="Chemical">nd OTU table. OTUs ideclass="Chemical">ntified as eukaryotes iclass="Chemical">n the 16S OTU table, class="Chemical">noclass="Chemical">n-fuclass="Chemical">ngi OTUs iclass="Chemical">n the ITS OTU table, as well as OTUs ideclass="Chemical">ntified as fuclass="Chemical">ngi, placlass="Chemical">nts, aclass="Chemical">nd class="Chemical">noclass="Chemical">n-soil aclass="Chemical">nimals were removed from the 18S OTU table. Read couclass="Chemical">nts from each group were class="Chemical">normalised usiclass="Chemical">ng rarefactioclass="Chemical">n. The OTU tables were rarefied 100 times usiclass="Chemical">ng class="Chemical">n class="Chemical">phyloseq[69] (as justified by Weiss et al.[70]) and the resulting mean richness was calculated for each sample. The read depth used for rarefaction varied for each group (Supplementary Table 14). Samples with lower read counts than this cut-off were removed before rarefaction. A summary of number of replicates per AVC is included in Supplementary Table 1.

Statistical analyses

All statistical analyses were run using R v. 3.3.3[71] using the rarefied data sets for each organismal group. The vegan package[72] was used to assess β-diversity via NMDS and CAP ordinations based on Bray–Curtis dissimilarities. A linear model for each environmental variable was fit separately to the ordination using the envfit function, the results are presented ranked according to goodness-of-fit. Results of goodness-of-fit for each variable from both ordination methods were compared using regression analyses to look for congruence. The values of all variables were plotted against NMDS scores to determine if there were positive or negative relationships with each NMDS axis. Differences in β-diversity amongst AVCs were calculated with PERMANOVA. The assumption of homogeneity of dispersion was verified using the betadisper function.

Linear mixed models were constructed using package nlme[73] to test the differences in α-diversity amongst AVCs for each organismal group. Model selection was performed using AVC, soil type, LOI classification, and sample year as fixed factors; sample square identity was the random factor. To determine the best possible model, predictors other than AVC were dropped to find the lowest AIC scores using the AICcmodavg package[74]. For each model, significant differences were assessed by ANOVA and pairwise differences were identified with Tukey’s post-hoc tests from the multcomp package[75].

n class="Disease">PLS regressioclass="Chemical">ns fouclass="Chemical">nd iclass="Chemical">n package class="Chemical">n class="Disease">pls[76] were used to identify the most important environmental variables for richness. Such analysis is ideal for data where there are many more explanatory variables than sample numbers or where extreme multicollinearity is present[46]. As in Lallias et al.[46], we used the variable importance in projection (VIP) approach[77] to sort the original explanatory variables by order of importance; variables with VIP values > 1 were considered most important. Relationships between important variables and richness values for each group of organisms were investigated by linear regression. Richness was normalised before regression when necessary. Pearson’s correlation coefficient was used to directly compare richness of organismal groups.