Small-scale agricultural grassland management can affect soil fungal community structure as much as continental scale geographic patterns.

A European transect was established, ranging from Sweden to the Azores, to determine the relative influence of geographic factors and agricultural small-scale management on the grassland soil microbiome. Within each of five countries (factor 'Country'), which maximized a range of geographic factors, two differing growth condition regions (factor 'GCR') were selected: a favorable region with conditions allowing for high plant biomass production and a contrasting less favorable region with a markedly lower potential. Within each region, grasslands of contrasting management intensities (factor 'MI') were defined: intensive and extensive, from which soil samples were collected. Across the transect, 'MI' was a strong differentiator of fungal community structure, having a comparable effect to continental scale geographic factors ('Country'). 'MI' was also a highly significant driver of bacterial community structure, but 'Country' was clearly the stronger driver. For both, 'GCR' was the weakest driver. Also at the regional level, strong effects of MI occurred on various measures of the soil microbiome (i.e. OTU richness, management-associated indicator OTUs), though the effects were largely regional-specific. Our results illustrate the decisive influence of grassland MI on soil microbial community structure, over both regional and continental scales, and, thus, highlight the importance of preserving rare extensive grasslands.

INTRODUCTION

Relatively, more studies in the existing literature have examined the effect of environmental factors (e.g. annual temperature and precipitation, soil pH) in driving the structure of the grassland soil microbiome across large spatial/ continental scales (e.g. Chen et al. 2015; Bahram et al. 2018; Xue et al. 2018; Plassart et al. 2019). In contrast, studies looking into the effects of grassland management intensity over comparable scales are currently missing. This is a pertinent research gap, particularly in light of the sharp decline in permanent, extensively managed grasslands over the past century. Before the green revolution, traditional grassland farming with low levels of farmyard manure application and low numbers of annual utilizations created such grasslands and represent some the most biodiverse habitats in temperate Europe (Wesche et al. 2012). They can contain a species-rich, compositionally complex plant community (Peter et al. 2009; van Dobben et al. 2017).

Land use intensification represents one of the principle drivers of global change (Erb et al. 2017), and extensively managed grasslands are threatened at both ends of the management intensity spectrum, through intensifcation and abandonment. Over the course of the 20th century, the intensification of European grasslands has resulted in their plant communities becoming comparatively simple, with a greatly reduced species richness (Wesche et al. 2012). This was achieved through high nutrient fertilizer inputs, greatly increased number of annual utilizations (i.e. cutting and/or grazing events) and even through the sowing of grass monocultures or two to four-species mixtures. This has resulted in a multitrophic homogenization of grasslands communities, including soil microbial communities (Gossner et al. 2016). On the other hand, abandonment of extensive grasslands threatens what remains of this grassland type (through reforestation) and its associated ecological value (Le Clec'h et al. 2019; Zehnder et al. 2020). This has now resulted in extensive, species rich grasslands being one of the most endangered habitats in Europe, with the need of their conservation representing an urgent societal challenge (Habel et al. 2013). Promoting and preserving this grassland type, and the traditional farming practices which preserve them, within the agricultural landscape is now a major aim of the common agricultural policy (CAP) and habitats directive (2014–2020) of the European Union (EU). The primary aim of extensive grasslands under such schemes is to reduce nutrient fertilizer application, while concomitantly enhancing their capacity for biodiversity promotion (Kampmann et al. 2012).

The decline of extensively managed grasslands has important implications for soil microbial biodiversity patterns, as grassland management intensity has been shown to be a strong determinant of soil microbial communities (de Vries et al. 2012b; Meyer et al. 2013; Sayer et al. 2013; Szukics et al. 2019). Due to very limited nutrient fertilizer inputs, extensively managed grasslands are typically considered nutrient-poor, with such environments being favorable to soil fungi (van der Heijden, Bardgett and van Straalen 2008). Indeed, management practices such as reduced fertilizer application, followed by reduced grazing and cutting events, have been shown to promote the abundance of fungal biomass (de Vries et al. 2007). By contrast, increased grassland management intensity has been shown to promote bacterially dominated soil microbial communities (Xun et al. 2018). This difference may be evoked by distinct mechanisms, such as differences in the quantity and quality of plant-derived substrates (i.e. plants residues and root exudates) entering the soil matrix due to contrasting plant community structures (Fox, Lüscher and Widmer 2020). Additionally, plant traits associated with intensively managed grasslands, e.g. high growth rate and competitive ability for nutrients, have been shown to promote bacterial abundance in grasslands (Orwin et al. 2010).

The majority of reported studies into the effects of grassland management intensity on the soil microbiome are at the experimental field or plot scale, while few have examined the influence of grassland management intensity at the regional or continental scale (de Vries et al. 2012b; Szukics et al. 2019; Felipe-Lucia et al. 2020; Andrade-Linares et al. 2021). Thus, whether management intensity or geographic pattern (and the associated differences in pedo-climatic conditions which determine growth) is the stronger driver of grassland soil microbial community structure remains unclear. Microbial diversity in soil has been shown to strongly exhibit geographic patterns (Martiny et al. 2006), with environmental factors being shown to be the strongest driver of bacterial community structure at the landscape scale (Mayerhofer et al. 2021). In contrast, soil bacterial diversity has been shown to be largely independent from geography, but highly influenced by plant diversity, a major differentiator between intensively and extensively managed grasslands (Fierer and Jackson 2006). While, in vineyard sites, it was shown that land use had a strong influence on fungal community structure, while geography was the stronger driver of the bacterial community (Coller et al. 2019).

Additionally, a more comprehensive understanding of the effect of grassland management intensity on the taxonomic and functional composition of the soil microbiome is also needed. Specific fungal taxa, which are known to be associated with ‘unimproved’ European grasslands include the Hygrocybe, Entoloma, Clavaria and Geoglossum (McHugh et al. 2001). More broadly, high soil nitrogen (N) content has been shown to suppress soil fungi (Fox, Lüscher and Widmer 2020). The abundance of bacterial and archeal genes involved in key steps in the soil N cycle, fixation, nitrification and denitrification have been shown to be altered by grassland management intensity (Meyer et al. 2013). Additionally, temperate European intensive grasslands typically receive fertilizer inputs in the form of slurry and manure, potentially introducing exogenous microorganisms into the soil matrix, such as bacteria of the phyla Firmicutes and Bacteriodetes (Abubaker et al. 2013).

Previous studies have examined various aspects of grassland soil biology at the European scale, especially comparing it with both arable and forest soils, as well as the associated environmental drivers (Bouffaud et al. 2016; Plassart et al. 2019), while studies comparing different grassland management types are missing. This study aimed to compare the effect size of grassland management intensity (a composite term comprising many factors) against the benchmark of environmental factors (also a composite term). To achieve this, our sampling strategy attempted to maximize the environmental conditions and management practices employed. The study did not aim to disentangle the variation explained by individual factors. The employed study design is not ideal for partitioning the effects of individual factors, given how regionally specific grassland management is in Europe (Gilhaus et al. 2017). Our study specifically focused on the influence of grassland management intensity on both fungal and bacterial community structure across both a continental and regional scale.

Specifically, the study had two research hypotheses:

Geographic factors (i.e. location and environmental variables), and not management intensity, are stronger determinants of both soil fungal and bacterial community structures in grasslands at the continental scale.

As differing management regimes will exist across the different transect regions, we hypothesize that there will be inconsistent effects of grassland management intensity across many regions of the continental transect.

MATERIALS AND METHODS

Survey design

A European transect, representing an agro-climatic gradient, was established comprising Sweden, Germany, Switzerland, Portugal (mainland) and the Azores (Fig. 1). For simplicity, these will, hereafter, be collectively referred to as Countries. A total of two separate growth condition regions (GCR, characterization see Table 1) within each country were selected, a favorable region (FR), where, due to environmental conditions, agricultural production was considered good for that country, and a contrasting less favorable region (LFR), where agricultural production was considered constrained by prevailing environmental conditions. There were 10 of these regions in total, which were varied in their biogeography, comprising the north (LFR, boreal) and south (FR, continental) of Sweden (regions contrast in temperature), the lowland (FR, continental) and mountain (LFR) areas of Germany and Switzerland (regions contrast in altitude), the north (FR) and south (LFR, both mediterranean) of Portugal (regions contrast in precipitation) and the north (LFR) and south (FR, both macronesia) of Azores (regions contrast in slope direction, Roekaerts 2002). Greater detail on the location of the study sampling sites is provided in Figure S1 (Supporting Information). The aim of this sampling strategy was 2-fold. Firstly, to extend the gradient by including contrasting environmental regions within each country. Secondly, splitting each country into two regions allowed for INT and EXT to be sampled within comparatively small regions, keeping the environmental conditions (and thus growth conditions) as comparable as possible and the confounding factors as low as possible.

Figure 1.

Overview of the European transect with the mean sampling area in a favorable region (FR, upright triangle) or less favorable region (LFR, downward triangle) in each country of the transect; Sweden (SE), Germany (DE), Switzerland (CH), Portugal (PT) and the Azores (AZ). Overlain (in color) are the biogeographic regions, to highlight the large variation in climatic/environmental conditions in the transect, while the dotted line was added to emphasize the transect scale. Greater detail on the sampling sites are provided in Figure S1 (Supporting Information).

Table 1.

Definition of growth condition region (GCR), FR or LFR, which was either mean annual temperature (°C), altitude (m), mean annual precipitation (mm) and slope direction. Also shown is the best practice recommendations for intensive (INT) and extensive (EXT) grasslands management intensities (MI) in each country of the transect (country and GCR codes in the legend of Fig. 1). This includes details on the expected utilizations (grazing was typically by cattle, sheep or some other ruminant) and nutrient application (kg/ha/yr, mineral, slurry or solid manure) in a year.

Country	GCR	Definition	MI	Utilizations/yr	Nutrient additions	Additional
SE	FR	∼7.5°C	INT	Grazed, 2–3 cuts	∼180 kg N/ha/yr (mineral + slurry)	Ley in a crop rotation (sown 2007– 2014)
	FR		EXT	Grazed, 2 cuts	None	Permanent grassland
	LFR	∼3°C	INT	Grazed, 2 cuts	∼55 kg N/ha/yr (mineral + slurry)	Ley in a crop rotation (sown 2012– 2014)
	LFR		EXT	Grazed, 2 cuts	None	Permanent grassland
DE	FR	<800 m	INT	Grazed, 5–6 cuts	∼230 kg N/ha/yr (mineral + slurry)	Permanent grassland
	FR		EXT	1–2 cuts	∼80 kg N/ha/yr (only solid manure)	Permanent grassland
	LFR	>900 m	INT	Grazed, 2 cuts	∼170 kg N/ha/yr(only slurry)	Permanent grassland
	LFR		EXT	1 cut	None	Permanent grassland
CH	FR	∼500 m	INT	5–7 cuts, may also be grazed	∼140 kg N/ha/yr (mineral + slurry)	Permanent grassland, first cut end of April
	FR		EXT	1–2 cuts, may also be grazed	None	Permanent grassland, first cut after 15th June
	LFR	∼1300 m	INT	3–4 cuts, may also be grazed	∼80 kg N/ha/yr (mineral + slurry)	Permanent grassland, first cut in May
	LFR		EXT	1 cut, may also be grazed	None	Permanent grassland, first cut after 15th July
PT	FR	>700 mm	INT	Grazed, no cuts	50–150 kg P/ha (mineral, at sowing)	Sown (2010–2015), irrigated
	FR		EXT	Grazed, no cuts	None	Semi-natural, minimal management
	LFR	<700 mm	INT	Grazed, no cuts	50–150 kg P/ha (mineral, at sowing)	Sown (2010–2015), irrigated
	LFR		EXT	Grazed, no cuts	None	Semi-natural, minimal management
AZ	FR	Southward facing slopes	INT	Grazed all year, 1 cut per month during summer	∼180 kg N/ha/yr mineral + slurry	3–5 years since last sowing
	FR		EXT	Summer grazing	None	Semi-natural pasture
	LFR	Northward facing slopes	INT	Grazed all year, 1 cut per month during summer	∼180 kg N/ha/yr mineral + slurry	3–5 years since last sowing
	LFR		EXT	Summer grazing	None	Semi-natural pasture

In each of the 10 GCR, 24 grassland sites across two contrasting grassland management intensities (MI) were chosen; 12 sites representing highly intensive (INT) managed grasslands and 12 sites representing extensively (EXT) managed grasslands. It should be noted that farmers adapt their grassland management practices to best suit regional conditions, therefore, the specific aspects of grassland management intensity will be regionally specific to each of the GCR in the transect. Though the primary difference between INT and EXT in all regions is the same, in that INT aims to increase the biomass yielding potential. Table 1 summarizes the agricultural practices characterizing the INT and EXT management type in each of the 10 GCR. We aimed at sampling the broadest range of, best practice, MI in each GCR, as the rational of this study was to compare the effect sizes of the biggest range occurring in environmental/climatic conditions (as ‘Country’ and ‘GCR’) with the effect size of the biggest range occurring in grassland management intensity (as ‘MI’). Generally, INT were characterized by high inputs of plant growth resources (N and/or P fertilizer, irrigation) combined with frequent and many utilizations throughout the growing season (i.e. cuttings and/or grazing events), while EXT received little/no resource input and had infrequent utilizations (EXT). In Portugal and Azores, the number of samples per region differed slightly. In Portugal, 8 INT and 14 EXT were sampled in FR and 16 INT and 10 EXT in LFR while in Azores, 10 INT and 8 EXT were sampled in FR and 12 INT and 6 EXT in LFR. Over the whole transect this resulted in 228 grassland sites sampled. A third grassland MI, low intensive (LI), characterized by a strongly reduced fertilization rate and infrequent utilizations, was also sampled, though for clarity purposes was removed from further analysis.

Soil samples were taken at the peak of the vegetation season. This represented the summer months of 2017 in Sweden, Germany, Switzerland and Azores. In Portugal, grassland phenology is distinctly different due to the Mediterranean climate, with peak biomass in spring. Therefore, the soil sampling in this country was delayed until the spring of 2018. In total, 228 sites were sampled as part of this study.

Soil sampling

A standardized field sampling design was employed for all sites of the transect. This structured study unit had four 4 m2 subplots, to get an optimum representation of the study site. A schematic of the sampling design employed is provided in Figure S2 (Supporting Information). From each corner of the subplots, a soil sample (ø 2.5 cm) was taken to a depth of 20 cm. The 16 collected cores were combined, homogenized and sieved (2 mm), and a subsample was frozen at −20°C, and later freeze-dried. All freeze-dried samples were then sent to Switzerland, which conducted the soil DNA extraction and later molecular work. There, ∼0.3 g of each sample was added to a 2 mL Eppendorf tube which contained 0.5 g of glass beads (ø 0.10–0.11 mm). Subsequently, 1.2 mL extraction buffer (0.2 M Na3PO4 at pH 8, 0.1 M NaCl, 50 mM EDTA, 0.2% CTAB) was added, and the tubes were vortexed and stored at −20°C.

Measurement and modeling of environmental variables

Total soil carbon (C) and nitrogen (N) was determined from air-dried, milled soil (<2 mm) using the FLASH 2000 Elemental analyzer (Thermo Scientific, Waltham, MA). For the determination of total soil P (Ptot), 300 mg of soil was burned (550°C, 3 h, Nabertherm B400) to remove the organic matter and the remaining ash was digested in 10 mL 13% HNO3 (1:4) and analysed using ICP-OES. Soil pH was measured in a 1:5 dH2O suspension. Soil texture (relative amounts of sand, silt and clay) was measured using the particle fractionation method, combining sieving and sedimentation (ISO 11277:2009). Climate variables were calculated using the CHELSA (Climatologies at high resolution for the earth's land surface areas) model (Karger et al. 2017). Variables relating to temperature (i.e. mean annual temperature, mean annual diurnal range, mean temperature of the warmest and coldest quarters) and precipitation (i.e. mean annual precipitation, precipitation of the wettest and driest quarters) were extracted.

DNA extraction and metabarcoding

DNA was extracted from each soil sample in three repeated extractions with a final phenol–chloroform extraction according to the procedure previously described (Bürgmann et al. 2001), alongside aforementioned modifications (Hartmann et al. 2005). DNA extracts were subsequently quantified using the PicoGreen dsDNA quantification assay (Invitrogen, Carlsbad, CA) and a Cary Eclipse fluorescence spectrophotometer (Varion Inc., Palo Alto, CA). Extracts were then diluted to a concentration of 5 ng/µL using sterile ddH2O. The primer pair ITS3 (5′ CAH CGA TGA AGA ACG YRG 3′) and ITS4 (5′ TCC TSC GCT TAT TGA TAT GC 3′) was used to PCR amplify the fungal internal transcribed spacer region (ITS2) of the rRNA operon (Tedersoo et al. 2014), while the primer pair 341F (5’ CCT AYG GGD BGC WSC AG 3’) and 806R (5’ GGA CTA CNV GGG THT CTA AT 3’) was used to PCR amplify the V3–V4 region of the bacterial 16S rRNA gene (Frey et al. 2016). The 5’ end of each forward primer was tagged with a CS1 adapter, with the reverse primer being tagged with CS2, which allowed multiplexing with the Fluidigm Access Array System (Fluidigm, San Francisco, CA). PCR amplification conditions were as described previously (Mayerhofer et al. 2017), with the modification that the reaction volume was 20 μL and 35 cycles of amplification were used. Negative PCR controls were also conducted to check for contamination. Each reaction was repeated four times, pooled and sent to the Génome Québec Innovation Center (Montréal, Canada) for PE-300 sequencing using the Illumina Miseq v3 platform (Illumina Inc., San Diego, CA). No PCR product could be obtained from five of the samples, due to a low quantity of DNA extracted, or to humic inhibitors being present in the sample (Matheson et al. 2010). These included four from Sweden (one each from EXT in FR and LFR, two from INT in FR) and one from Germany (INT in LFR), resulting in 223 sites being analysed.

Sequence processing and taxonomic classification

The processing of obtained sequences was done using a customized pipeline in UPARSE implemented in USEARCH version 9 (Edgar 2010) with the specifications described previously (Frey et al. 2016; Gschwend et al. 2020). The following comprised the main steps of the bioinformatic workflow; merging of paired ends, quality filtering (maximum expected error per sequence ≤1), removal of singletons, operational taxonomic units (OTU) centroids clustering at 97% identity and target verification using METAXA2 (Bengtsson-Palme et al. 2015) and ITSx (Bengtsson-Palme et al. 2013) for prokaryotic and eukaryotic sequences respectively. Taxonomic classification was performed using the Bayesian classifier implemented in MOTHUR version 1.36.1 (Schloss et al. 2009) with a custom-made database extracted from NCBI GenBank (Benson et al. 2014) for eukaryotic sequences. The SILVA ribosomal RNA gene database (version 132; Quast et al. 2012) was used as a reference for prokaryotic sequences. Non-fungal sequences were discarded from the eukaryotic dataset and retained fungal OTUs (fOTUs) were assigned using the UNITE reference database (version 7.2; Abarenkov, Henrik Nilsson and Larsson 2010). The minimum bootstrap value for a taxonomic assignment was set to 80%. Lastly, prokaryotic OTUs which were assigned to chloroplasts, mitochondria and archaea were removed to retain only bacterial OTUs for analyses.

Data analysis

All statistics were conducted in the R statistical software program, unless otherwise stated (R Development Core Team 2019). Both the fungal and bacterial OTU count tables were square-root transformed, converted to relative abundance and a pairwise Bray–Curtis dissimilarity matrix constructed (Hartmann et al. 2015). Canonical analysis of principal coordinates (CAP) was applied to assess the discriminatory power of MI on both fungal and bacterial community structures at the transect level using the function ‘CAPdiscrim’ in the R package ‘BiodiversityR’ (Anderson and Willis 2003). This CAP analysis revealed that LI straddled the two management extremes (Figure S3, Supporting Information), indicating the sampled sites were good representations of INT and EXT grasslands. Therefore, for clarity purposes, we focused all further comparative analysis between INT and EXT, and excluded LI. CAP analysis was also performed to discriminate the fungal and bacterial community structures, and the environmental variables, the latter scaled and centered and a Euclidean distance matrix constructed (Gschwend et al. 2020) between INT and EXT in both GCR within each country. The CAP analysis is a constrained ordination analysis which can be done on the basis of any distance or dissimilarity measure, and can uncover patterns of difference which are masked in an unconstrained ordination. To determine how the countries of the transect varied from each other purely in terms of their environmental properties, certain variables relating to soil chemistry were not considered (i.e. total % soil C and N, pH and Ptot), as these could be influenced by agricultural management. Pairwise Spearman correlations were conducted on the remaining variables, with the following being selected (Spearman rho ≥ 0.6); mean annual temperature, mean annual temperature diurnal range, mean annual precipitation, % sand, altitude and site gradient. The variables removed were; the mean temperature of the warmest and coldest quarters, mean precipitation of the wettest and driest quarters as well as the % soil silt and clay, as they were highly correlated to the variables already in the model. A principle component analysis (PCA) was then carried out using ‘prcomp’ in the stats package, on these scaled and centered variables. While the principle objective of measuring the environmental variables was to determine whether such differences within regions would compound the grassland management effect on soil microbial community structure, a distance-based linear modeling (DISTLM) analysis (AIC selection criterion) was still conducted to determine the influence of individual variables on both fungal and bacterial community structure. The influence of the selected environmental variables on both fungal and bacterial community structure were also visualized using the distance based redundancy analysis (dbRDA). Both the DISTLM and dbRDA analysis were conducted in the PRIMER-7 software.

Variation partitioning among ‘Country’, ‘GCR’ and ‘MI’ at the transect level was assessed via permutational multivariate analysis of variance (PERMANOVA, Monte Carlo, 9999 permutations, Anderson 2001) using a doubly-nested design (i.e. GCR nested within Country and MI nested within GCR) in the statistical software PRIMER-E (Version 7, Plymouth, UK). A nested (rather than factorial) statistical design was chosen, as the management types within each region of the transect are not equivalent. The centroid distance (in Euclidean ordination space) of the environmental variables, as well as the fungal and bacterial community structures, between ‘Country’, ‘GCR’ and ‘MI’ were calculated using the ‘betadisper’ function (Anderson, Ellingsen and McArdle 2006; Gschwend et al. 2020), and used to determine the effect size of each factor, along with the square-root of the component of variation (√CV) derived from the PERMANOVA model. The effect of GCR and MI within each country was tested via PERMANOVA, with a nested design (MI nested within GCR). A pairwise test in the PERMANOVA function was used to indicate statistical differences between INT and EXT in each region.

Fungal and bacterial community OTU richness was calculated using the ‘summary.single’ command in MOTHUR V1.36.1 (Schloss et al. 2009). Differing sequencing depth in the respective communities was corrected for, using an iterative subsampling approach (10 000 x) to the lowest sequencing number (15 769 for fungi and 10 894 for bacteria, Figure S4, Supporting Information). The effect of ‘MI’ on fungal and bacterial OTU richness, in each of the ten regions was statistically assessed using a t-test, this was also done to assess the effect of ‘GCR’ and ‘MI’ on each of the selected environmental variables, as well as those which had been removed from the PCA. The effect of ‘Country’ on these variables was assessed via analysis of variance (ANOVA). In instances where the assumption of normality was not met, the non-parametric Kruskal–Wallis test was used. The effect of MI on the relative abundance of the major fungal and bacterial phyla (> 0.1% rel. total sequence abundance) was measured as a log response ratio (Hedges, Gurevitch and Curtis 1999).

Finally, associations of fungal and bacterial OTUs (fOTU and bOTU, respectively) to MI both within each of the ten regions, and at the transect level, was determined via correlation based indicator species analysis (‘multiplatt’ function) with 9999 permutations in the package indicspecies (Cáceres and Legendre 2009). Prior to analysis, OTUs which were present on just one site were removed from the dataset. An indicator OTU was defined as those with an IndVal ≥ 0.7 and a P-value ≤ 0.05. Functional information on the MI indicator fOTU (at the transect level) was obtained through the FUNGuild software (Nguyen et al. 2016).

RESULTS

Soil physicochemical and climatic variables across the European transect

The factor ‘Country’ was by far the strongest driver of site differences in environmental variables across the transect (mean centroid distance = 3.086, P ≤ 0.001) followed by ‘GCR’ (1.320, P ≤ 0.001 Fig. 2A, Table 2a). This indicates that the sampling design did incorporate a large range of environmental and growth conditions, as was the intention. The environmental variables that most strongly differentiated the countries of the transect, were temperature (country mean ranged from 5.22°C in Sweden to 16.43°C in Portugal), precipitation (652 mm in Sweden to 1484 mm in Germany) and altitude (111 m in Sweden to 917 m in Germany, all Table 3). Azores could be clearly differentiated from the other countries, while Germany and Switzerland were the most similar to each other (Fig. 2A). The factor ‘GCR’ had a strong effect on site differences (P ≤ 0.001, mean centroid distances of 1.945, 1.250 and 2.461, respectively, Fig. 2A, Table 4a) in Sweden, Germany and Switzerland, while it had no significant influence in Portugal and Azores (P > 0.05, mean centroid distance = 0.761 and 0.182, respectively). The mean annual temperature and precipitation significantly differed between the GCR of each country (all at least P ≤ 0.01), with the exception of Azores, and Germany for the latter. Moreover, the GCR of Germany and Switzerland significantly differed in their altitude and site inclination, with both variables increasing in the LFR (Table 3, all at least P ≤ 0.05).

Figure 2.

(A) Principal component analysis (PCA) displaying how the selected environmental variables differentiate the different regions of the transect, i.e. the five countries; SE (green), DE (blue), CH (red), PT (black) and AZ (purple), and their constituent GCR regions; FR (dark color) LFR (light color). Selected environmental variables were: mean annual temperature, mean annual diurnal temperature range, mean annual precipitation, % soil sand, altitude and site inclination. (B and C) Differentiation of both fungal (B) and bacterial (C) community structure in the different regions. Circles indicate soil samples that were correctly classified, while triangles indicate samples that were incorrectly classified to its sampled region in a leave-one-out reclassification test. Reclassification and constrained ordination are based on canonical analyses of principal coordinates (CAP) when maximizing differences between regions. Axes show linear discriminants (LD). Squares represent the community centroid for each of the 10 regions. Also displayed in each panel is the community centroid distance for each variable (Euclidian distance) between each country (left) and between the GCR in each country (right). Explanations for abbreviations can be found in Fig. 1 and Table 1

Table 2.

The effect of Country, GCR and MI on the measured environmental variables (i.e. selected soil physicochemical and climatic variables) and on both the fungal and bacterial community structure in the transect, as tested by nested PERMANOVA analysis. At the survey level, the PERMANOVA model was doubly nested, with GCR nested within country and MI within GCR. Displayed is the degree of freedom (df), the mean sum of squares (MS), the Pseudo F value (Pseudo F), square-root of the component of variation (√CV), the mean distance (Mean distance) and distance range (Distance range). Significance codes: ‘***’ P< 0.001, ‘*’ P < 0.05 and ‘ns’P > 0.05. Explanations for abbreviations can be found in Fig. 1 and Table 1.

	Factor	df	MS	Pseudo F	√CV	Mean distance	Distance range
(a) Environment
	Country	4	206.640	7.502***	2.035	3.086	1.815–3.759
	GCR (Country)	5	27.637	4.407***	0.992	1.320	0.182–2.461
	MI (GCR)	10	6.285	4.760***	0.675	0.959	0.202–1.832
	Residual	203	1.320		1.149
(b) Fungi
	Country	4	4.020	4.310***	0.247	0.429	0.280–0.493
	GCR (Country)	5	0.935	0.992ns	0.000	0.273	0.142–0.383
	MI (GCR)	10	0.945	4.505***	0.233	0.415	0.209–0.555
	Residual	203	0.210		0.482
(c) Bacteria
	Country	4	3.377	4.921***	0.249	0.388	0.229–0.499
	GCR (Country)	5	0.688	1.479*	0.101	0.227	0.115–0.339
	MI (GCR)	10	0.466	3.367***	0.174	0.289	0.141–0.448
	Residual	203	0.138		0.372

Table 3.

Difference (mean ± se) in selected environmental variables between countries, GCR and MI across the transect. Variables include mean annual temperature (Temperature, °C), diurnal temperature range (Diurnal temp range, °C), mean annual precipitation (Precipitation, mm), % soil sand (Sand), altitude (Altitude, m) and slope inclination (Inclination, as a decimal). Also reported is the F-value from an ANOVA test, on the influence of country on each variable, along with the associated significance. Explanations for abbreviations can be found in Fig. 1 and Table 1. Significance code: ‘***’ P ≤ 0.001 ‘**’ P ≤ 0.01 ‘*’ P ≤ 0.05 ‘⁑’ P ≤ 0.10 and ‘ns’ P > 0.10.

	Temperature	Diurnal temp range	Precipitation	Sand	Altitude	Inclination
Country	172.790***	157.550***	163.910***	71.643***	163.680***	95.06***
SE	5.22°C (±0.335)	5.12°C (±0.116)	652 mm (±10.27)	43% (±3.46)	111 m (±6.32)	0.03 (±0.003)
FR∼LFR	***	***	***	***	ns	ns
INT_FR	7.39°C (±0.097)	4.79°C (±0.236)	728 mm (±12.02)	62% (±3.55)	124 m (±12.16)	0.02 (±0.004)
EXT_FR	7.57°C (±0.098) ns	4.36°C (±0.225)ns	700 mm (±13.10)ns	58% (±2.03)ns	103 m (±12.12)ns	0.02 (±0.005) ns
INT_LFR	3.12°C (±0.107)	5.65°C (±0.065)	602 mm (±6.74)	25% (±5.73)	110 m (±12.87)	0.03 (±0.008)
EXT_LFR	3.17°C (±0.116) ns	5.61°C (±0.070)ns	590 mm (±8.29)ns	29% (±6.53)ns	108 m (±14.08)ns	0.03 (±0.007) ns
DE	7.11°C (±0.170)	8.10°C (±0.009)	1484 mm (±42.04)	25%(±1.66)	917 m (±29.09)	0.11 (±0.019)
FR∼LFR	***	***	⁑	***	***	*
INT_FR	7.85°C (±0.065)	8.06°C (±0.016)	1395 mm (±62.38)	27% (±2.24)	777 m (±10.91)	0.04 (±0.010)
EXT_FR	8.06°C (±0.086)⁑	8.04°C (±0.015) ns	1332 mm (±67.35)ns	33% (±2.81) ns	747 m (±17.04)ns	0.13 (±0.051)ns
INT_LFR	6.82°C (±0.153)	8.14°C (±0.007)	1521 mm (±57.91)	24% (±3.79)	978 m (±24.28)	0.08 (±0.018)
EXT_LFR	5.66°C (±0.336)*	8.16°C (±0.005)*	1690 mm (±104.67)ns	15% (±2.08)*	1169 m (±49.33)**	0.19 (±0.040)*
CH	7.83°C (±0.320)	7.80°C (±0.015)	974 mm (±35.00)	35% (±1.79)	876 m (±60.91)	0.25 (±0.024)
FR∼LFR	***	***	**	***	***	*
INT_FR	9.86°C (±0.133)	7.74°C (±0.011)	1077 mm (±50.47)	30% (±3.20)	491 m (±17.35)	0.12 (±0.024)
EXT_FR	9.82°C (±0.137)ns	7.73°C (±0.013) ns	1090 mm (±38.70)ns	25% (±3.46) ns	482 m (±20.50)ns	0.24 (±0.040)*
INT_LFR	5.94°C (±0.407)	7.86°C (±0.031)	848 mm (±72.26)	42% (±2.62)	1226 m (±68.76)	0.19 (±0.033)
EXT_LFR	5.69°C (±0.271)ns	7.86°C (±0.031) ns	880 mm (±85.06)ns	42% (±2.41) ns	1306 m (±49.56)ns	0.43 (±0.043)***
PT	16.43°C (±0.134)	7.82°C (±0.281)	726 mm (±17.37)	58% (±3.17)	185 m (±21.45)	0.03 (±0.006)
FR∼LFR	**	***	***	ns	ns	ns
INT_FR	15.83°C (±0.455)	8.24°C (±0.538)	840 mm (±45.37)	68% (±3.40)	227 m (±77.96)	0.01 (±0.005)
EXT_FR	16.07°C (±0.281)ns	8.88°C (±0.323) ns	789 mm (±34.42)ns	53% (±6.55)*	229 m (±46.68)ns	0.03 (±0.009)*
INT_LFR	17.09°C (±0.028)	6.55°C (±0.562)	647 mm (±15.48)	66% (±5.57)	102 m (±15.51)	0.02 (±0.005)
EXT_LFR	16.37°C (±0.177)***	8.01°C (±0.520)ns	672 mm (±5.33)ns	45% (±5.58)*	222 m (±29.70)**	0.08 (±0.015)***
AZ	14.64°C (±0.199)	1.212°C (±0.0003)	1284 mm (±19.78)	46% (±1.15)	574 m (±35.43)	0.10 (±0.012)
FR∼LFR	ns	**	ns	ns	ns	ns
INT_FR	15.17°C (±0.276)	1.212°C (±0.0005)	1269 mm (±34.14)	46% (±2.75)	465 m (±41.02)	0.08 (±0.017)
EXT_FR	13.45°C (±0.207)***	1.210°C (±0.0002) ns	1362 mm (±18.67)*	46% (±1.76)ns	793 m (±46.06)***	0.14 (±0.023)⁑
INT_LFR	15.64°C (±0.258)	1.213°C (±0.0006)	1184 mm (±30.84)	45% (±2.00)	400 m (±42.07)	0.07 (±0.021)
EXT_LFR	13.49°C (±0.091)***	1.212°C (±0.0005)ns	1373 mm (±26.77)***	48% (±1.79)ns	789 m (±25.71)***	0.13 (±0.034) *

Table 4.

The effect of both GCR (FR and LFR) and MI within each of the five countries of the transect; SE, DE, CH, PT and AZ, as measured by nested PERMANOVA. At the country level, MI is nested within GCR. Displayed is df, the Pseudo F value, the P-value. Also displayed are separate pairwise PERMANONA tests displaying the differences (or lack thereof) in community structure between INT and EXT (INT ∼ EXT) in each of the ten regions of the transect. Explanations for abbreviations can be found in Fig. 1 and Table 1 and 2. Significance codes: ‘***’ P< 0.001 ‘**’ P < 0.01 ‘*’ P < 0.05 and ‘ns’ P > 0.05.N

nested PERMANOVA				Pairwise PERMANOVA
	Factor	df	MS	Pseudo F	√CV	GCR	Factor	t-value
(a) Environment
SE	GCR	1	121.790	53.030***	2.326	FR	MI	1.210ns
	MI (GCR)	2	2.295	0.698ns	0.000	LFR	MI	0.357ns
	Residual	40	3.286		1.800
DE	GCR	1	100.100	7.239*	1.917	FR	MI	1.370ns
	MI (GCR)	2	13.833	4.088**	0.943	LFR	MI	2.418**
	Residual	43	3.384		1.840
CH	GCR	1	132.740	15.599**	2.275	FR	MI	1.496⁑
	MI (GCR)	2	8.509	2.831**	0.677	LFR	MI	1.767*
	Residual	44	3.006		1.734
PT	GCR	1	39.064	2.063 ns	0.950	FR	MI	1.012ns
	MI (GCR)	2	19.063	4.139**	1.134	LFR	MI	3.106***
	Residual	44	4.606		2.146
AZ	GCR	1	6.996	0.200 ns	0.000	FR	MI	2.228**
	MI (GCR)	2	35.096	8.380***	1.898	LFR	MI	3.654***
	Residual	32	4.188		2.059
(b) Fungi
SE	GCR	1	1.403	1.927*	0.176	FR	MI	2.185***
	MI (GCR)	2	0.729	3.745***	0.221	LFR	MI	1.674**
	Residual	40	0.195		0.441
DE	GCR	1	0.833	1.367ns	0.098	FR	MI	1.159ns
	MI (GCR)	2	0.610	2.786***	0.183	LFR	MI	1.983***
	Residual	43	0.219		0.468
CH	GCR	1	1.757	1.884⁑	0.185	FR	MI	2.552***
	MI (GCR)	2	0.933	4.408***	0.245	LFR	MI	1.730**
	Residual	44	0.212		0.460
PT	GCR	1	0.531	0.457ns	0.000	FR	MI	2.117***
	MI (GCR)	2	1.169	5.091***	0.245	LFR	MI	2.394***
	Residual	44	0.230		0.495
AZ	GCR	1	0.152	0.118ns	0.000	FR	MI	2.701***
	MI (GCR)	2	1.284	6.882***	0.361	LFR	MI	2.542***
	Residual	32	0.187		0.430
(c) Bacteria
SE	GCR	1	1.243	3.322**	0.199	FR	MI	1.738**
	MI (GCR)	2	0.375	2.930***	0.150	LFR	MI	1.679***
	Residual	40	0.128		0.358
DE	GCR	1	0.350	1.461ns	0.069	FR	MI	0.980ns
	MI (GCR)	2	0.240	1.735⁑	0.093	LFR	MI	1.534*
	Residual	43	0.138		0.372
CH	GCR	1	1.383	4.574**	0.212	FR	MI	2.063***
	MI (GCR)	2	0.302	2.251*	0.118	LFR	MI	1.110ns
	Residual	44	0.134		0.366
PT	GCR	1	0.356	0.559ns	0.000	FR	MI	1.886**
	MI (GCR)	2	0.642	4.003***	0.176	LFR	MI	2.115***
	Residual	44	0.160		0.411
AZ	GCR	1	0.109	0.141ns	0.000	FR	MI	2.313***
	MI (GCR)	2	0.772	6.047***	0.276	LFR	MI	2.586***
	Residual	32	0.128		0.357

The factor ‘MI’ had the weakest influence on environmental variables across the transect (mean centroid distance = 0.959, P ≤ 0.001). Within countries, ‘MI’ had no significant impact in Sweden (Table 4a), while a significant effect was observed in Germany and Switzerland (both P ≤ 0.01), ‘MI’ was not as strong a driver as ‘GCR’ in those countries (Table 4a). In the southern countries of the transect (Portugal and Azores), ‘MI’ was the sole significant driver of the site differences in environmental variables (both P ≤ 0.01, Table 4a). In five of the 10 GCR, there was a significant difference between INT and EXT in terms of environmental variables (all at least P ≤ 0.05, Table 4a). Site inclination was the most common variable, which significantly differed between INT and EXT, especially in Switzerland where there was a sharp increase in the EXT sites (INT = 0.12 and 0.19, EXT = 0.24 and 0.43, in the FR and LFR, respectively, Table 3). The variable Ptot, which could be influenced directly by MI, was significantly higher in INT, as compared to EXT, in six of the ten GCR, including both GCR of Switzerland and Azores (all P ≤ 0.05, Table S1, Supporting Information).

Microbiome structures across the European transect

A total of 11 735 475 high-quality partial fungal ITS sequences were obtained, with a mean Good's coverage per sample of 0.99, and yielding 19 330 fungal OTUs (fOTU). The two factors ‘Country’ and ‘MI’ were strong drivers of fungal community structure (P ≤ 0.001) with comparable effect sizes (mean centroid distances of 0.429 for ‘Country’ and 0.415 for ‘MI’, Table 2b). The same pattern was evident for √CV (‘Country’ = 0.247, ‘MI’ = 0.233; Table 2b). In contrast, ‘GCR’ had no effect on fungal community structure when analysed over the whole gradient (P = 0.525, Table 2b). For fungal community structure, there was a clear separation of the southern countries, Portugal and Azores, from the other countries of the transect (i.e. Sweden, Germany and Switzerland, centroid distance = 0.430–0.493, Fig. 2B) as well as from each other (centroid distance = 0.468). Within countries, the factor ‘GCR’ was a significant driver of the fungal community structure only in Sweden (Table 4b), although it just failed to reach significance in Switzerland, (P = 0.056), while it was clearly insignificant for Germany, Portugal, and Azores. Within every single country, ‘MI’ was a highly significant driver of fungal community structure (P ≤ 0.001, Table 4b). Looking at the results from the DISTLM analysis, the temperature variables of Temperature (proportion = 0.083, P ≤ 0.001) and Diurnal temp range (proportion = 0.058, P ≤ 0.001) were the strongest environmental drivers of fungal community structure (Table S1a and Figure S5a, Supporting Information).

For bacteria, 9 414 029 high quality partial 16S rRNA gene sequences were obtained (mean Good's coverage per sample of 0.96), yielding 30 803 bacterial OTUs (bOTU). ‘Country’ was clearly the strongest driver of bacterial community structure at the transect level, (P ≤ 0.001, mean centroid distance = 0.388, Table 2c) followed by ‘MI’ (P ≤ 0.001, mean centroid distance = 0.289), while ‘GCR’ was clearly the weakest driver (P = 0.011, mean centroid distance = 0.227). The bacterial community structure of Portugal and Azores strongly separated from each other (centroid distance = 0.449, Fig. 2C), as well as from Germany (centroid distance = 0.414 and 0.439, respectively) and Switzerland (centroid distance = 0.398 and 0.499, respectively), though their separation from Sweden was not as strong (centroid distance = 0.380 and 0.324, respectively). As seen with fungi, within countries, ‘GCR’ was a significant driver of bacterial community structure only in Sweden and Switzerland (Table 4c), while it was not significant in Germany, Portugal, and Azores. While a strong driver at the transect level, within countries, however, ‘MI’ was the weaker driver of bacterial community structure in Sweden (√CV = 0.150) and Switzerland (√CV = 0.118) as compared to GCR (√CV = 0.199 and √CV = 0.212, respectively, Table 4c), though it was, by far, the stronger driver in the southern countries of Portugal and Azores. Like with the fungi, the temperature variables of Diurnal temp range (proportion = 0.097, P ≤ 0.001) and Temperature (proportion = 0.064, P ≤ 0.001) were the strongest of the selected environmental variables in driving bacterial community structure (Table S1b and Figure S5b, Supporting Information). Among the environmental variables which could be influenced by grassland management intensity, soil pH was, by far, the strongest driver of soil bacterial community structure (proportion = 0.194, P ≤ 0.001, Table S1b, Supporting Information).

Effect of management intensity on microbiome structures at the regional level

At the regional level within each country, there was a significant difference (mostly P ≤ 0.001) in the fungal community structure between INT and EXT in nine of the ten regions (Table 4b and Fig. 3F–K). The distance of community centroids in fungal community structure between INT and EXT was greatest in the GCRs of the two southern countries, Azores (FR = 0.547 and LFR = 0.555) and Portugal (FR = 0.471 and LFR = 0.443). In the remaining regions with a significant difference between INT and EXT, centroid distances ranged from 0.314 to 0.430 (Fig. 3F–K).

Figure 3.

Differentiation of selected environmental variables (A–E), as well as the fungal (F–K) and bacterial (L–P) community structure in INT (red) and EXT (blue) in each country of the transect; SE, DE, CH, PT and AZ, and in their constituent GCR; FR (dark color) and LFR (light color). Circles indicate soil samples that were correctly classified, while triangles indicate samples that were incorrectly classified to its sampled region ∼ management intensity in a leave-one-out reclassification test. Reclassification and constrained ordination are based on canonical analyses of principal coordinates (CAP) when maximizing differences between management intensities (MI). Squares represent the community centroid for each MI in each GCR. The number within each graph represents the community centroid distance (Euclidian) between INT and EXT. Black and grey numbers represent the difference between INT and EXT in the FR and LFR region within each country, respectively. These centroid distances (for each variable) are directly comparable across panels, and across graphs. Axes show linear discriminants (LD). Explanations for abbreviations can be found in Fig. 1 and Table 1.

There was a significant difference in bacterial community structure between INT and EXT (mostly P ≤ 0.01, Table 4c) in eight of the ten regions. The two exceptions were FR in Germany (P = 0.424), and LFR in Switzerland (P = 0.251). As with fungi, the centroid distance in the bacterial community structure between INT and EXT was greatest in the two southern countries, Azores (FR = 0.404 and LFR = 0.448) and Portugal (FR = 0.354 and LFR = 0.324). While in Sweden, Germany and Switzerland, the significant centroid distances in bacterial community structure between INT and EXT ranged from 0.141 to 0.277 (Fig. 3L–P).

In Sweden, Germany, and Switzerland, fungal OTU richness was higher in all six EXT locations, compared to INT, with increases being significant in four of the six locations (P ≤ 0.05, Table 5). In the southern countries of the transect (i.e. Portugal and Azores), fungal OTU richness responded in the opposite manner at three of the four remaining locations, i.e. FR in Portugal as well as both FR and LFR in Azores, with EXT having a significantly lower value (P ≤ 0.05), as compared to INT. The only significant (P ≤ 0.05) MI-associated differences in bacterial OTU richness were in the same three regions (FR in Portugal, both FR and LFR in Azores) where EXT harbored significantly lower bacterial OTU richness than INT (Table 5).

Table 5.

The average (±1 se) fungal OTU richness in the two grassland MI; INT and EXT in each of the 10 regions of the transect. Also displayed is the result of an ANOVA test to determine differences in both variables between INT and EXT. Significance codes: ‘***’ P< 0.001 ‘**’ P < 0.01 ‘*’ P < 0.05 and ‘ns’ P > 0.05. Additionally, the number of indicator (IndVal ≥ 0.7, P-value ≤ 0.05) fungal and bacterial OTUs (fOTU and bOTU, respectively) for the five countries of the transect in their constituent GCR. The percentage of f and bOTU which are shared between the different GCR of the transect are shown in Tables S4 and S5, Supporting Information, respectively. Explanations for abbreviations can be found in Fig. 1 and Table 1.

	SE		DE		CH		PT		AZ
	FR	LFR	FR	LFR	FR	LFR	FR	LFR	FR	LFR
(a) Fungi
INT richness	638	618	604	631	697	721	813	690	633	666
	(±18.46)	(±23.35)	(±25.62)	(±22.37)	(±19.10)	(±26.69)	(±10.71)	(±20.09)	(±26.67)	(±39.62)
EXT richness	780	637	687	653	833	862	653	789	518	485
	(±37.38)**	(±23.56)ns	(±25.39)*	(±50.40)ns	(±16.48)***	(±48.74)*	(±33.60)***	(±44.83)⁑	(±39.38)*	(±40.11)**
Indicator OTU INT	152	85	9	163	169	120	426	163	150	310
shared FR∼LFR (%)	10.97		1.74		11.07		11.04		24.57
Indicator OTU EXT	144	62	17	83	257	120	71	296	267	226
shared FR∼LFR (%)	6.31		1.00		4.24		8.99		27.59
(b) Bacteria
INT OTU richness	2168	2373	2154	2013	2047	2039	2574	2374	2510	2448
	(±119.22)	(±144.48)	(±52.63)	(±66.33)	(±59.26)	(±51.23)	(±55.89)	(±89.48)	(±122.39)	(±141.52)
EXT OTU richness	2201	2218	2025	1799	1960	1924	2040	2228	2039	2050
	(±117.20)ns	(±114.01)ns	(±63.79)ns	(±105.66)ns	(±49.33)ns	(±47.75)ns	(±117.32)**	(±57.76)ns	(±33.38)**	(±111.70)*
Indicator OTU INT	358	344	49	313	291	151	1161	490	600	1277
shared FR∼LFR (%)	8.26		2.21		7.69		11.51		22.22
Indicator OTU EXT	145	198	26	76	228	43	191	687	796	783
shared FR∼LFR (%)	3.50		0.98		0.00		7.90		27.61

Taxonomic differences with management intensity at the regional level

None of the seven major fungal phyla (> 0.1% rel. total sequence abundance) responded to MI across all 10 GCR, though five did respond in at least two of the 10 GCR (Table S2, Supporting Information). Of the 15 major bacterial phyla detected (> 0.1% rel. total sequence abundance), all displayed a significant difference to MI in at least one of the 10 GCR (Table S2, Supporting Information), with the exception of the Chloroflexi. A total of two phyla, the Firmicutes and Nitrospirae, displayed significant but inconsistent directions of differences in their relative abundance according to MI in seven and six out of the 10 GCR, respectively. The relative abundance of both the Firmicutes and Nitrospirae was significantly lower in EXT, compared to INT, in five and four of the GCR, respectively, with a mean reduction of −63% in the Firmicutes and −59% in the Nitrospirae. However, the relative abundance of the Firmicutes was significantly higher in EXT in the FR of Sweden and Germany, in the latter case by +154% and the Nitrospirae had a mean relative abundance +71% higher in EXT in both regions of Azores.

Fungal indicator OTUs associated with either INT or EXT were found in each of the ten regions (9 OTU–310 OTU, Table 5), with the smallest number being found in the FR of Germany. A very low percentage of indicator fOTUs for both INT and EXT were shared between the two GCR within the same country (Mostly < 10%, Table 5) and likewise, between the GCR of different countries (Tables S3 and S4, Supporting Information). Azores was somewhat of an exception. Here, about 25% of the indicator fOTUs for INT or EXT were shared between both regions. As was also seen with indicator fOTUs, a small percentage (mostly >10%) of indicator bOTU for both INT and EXT were shared among GCR, both between and within countries (Table 5, Tables S3 and S4, Supporting Information). Azores again had the highest percentage of bOTU shared between its constituent GCR (∼25%, Table 5).

When looking for MI-associated indicator fOTUs at the transect level, 18 fOTUs were found (Table 6). Of the five indicator fOTUs detected for EXT at the transect level, most were functionally assigned as saprotrophs. A total of four of these fOTU were assigned to two genera, Clavaria (fOTU_636 and fOTU_647) and Leohumicola (fOTU_395 and fOTU_9761). A total of two of the indicator fOTUs for INT were functionally assigned as dung saprotrophs (OTUf_13 and OTUf_106). At the transect level, 22 bOTUs were found for INT (Table 6), but none were an indicator for EXT. All detected indicator OTU at the transect level were found in the majority of sites of their associated management types (Table 6).

Table 6.

Fungal and bacterial OTUs which were detected as unique indicator OTUs for either MI; INT and EXT through indicator species analysis (IndVal > 0.7, P ≤ 0.05). The taxonomy of each detected OTU is given at the lowest known taxonomic level:—k_: Kingdom, c_: Class, o_: Order, f_: Family, g_: genus and s_: species. Given is the mean relative abundance (%, ± 1 se) for each OTU (rel ab.) as well as the relative abundance (%) of each OTU in both INT and EXT. Also given is the number of sites of a given MI the OTU was detected in (INT and EXT had 115 and 108 sites, respectively). Significance code: ‘***’ P ≤ 0.001. Explanations for abbreviations can be found in Table 1.

Indicator OTU	Taxonomy	IndVal	MI	Rel. ab. (%) INT	Rel. ab. (%) EXT	Sites
(a) Fungi
fOTU_46	f__Nectriaceae	0.869***	INT	0.746 (±0.113)	0.132 (±0.034)	100/115
fOTU_134	s__Microascaceae sp	0.846***	INT	0.831 (±0.231)	0.041 (±0.015)	85/115
fOTU_13	s__Thelebolus globosus	0.830***	INT	0.824 (±0.237)	0.049 (±0.018)	87/115
fOTU_106	s__Podospora intestinacea	0.808***	INT	0.828 (±0.710)	0.044 (±0.016)	78/115
fOTU_431	s__Sordariales sp	0.792***	INT	0.775 (±0.132)	0.101 (±0.037)	80/115
fOTU_2295	s__Glomeraceae sp	0.790***	INT	0.721 (±0.115)	0.159 (±0.042)	87/115
fOTU_2815	k__Fungi	0.744***	INT	0.797 (±0.131)	0.077 (±0.023)	69/115
fOTU_576	o__Pleosporales	0.740***	INT	0.728 (±0.123)	0.150 (±0.035)	76/115
fOTU_966	o__Pleosporales	0.740***	INT	0.704 (±0.136)	0.176 (±0.046)	78/115
fOTU_873	o__Helotiales	0.735***	INT	0.689 (±0.112)	0.192 (±0.050)	79/115
fOTU_318	c__Sordariomycetes	0.708***	INT	0.723 (±0.206)	0.157 (±0.046)	69/115
fOTU_1580	k__Fungi	0.706***	INT	0.687 (±0.116)	0.194 (±0.046)	73/115
fOTU_506	o__Pleosporales	0.705***	INT	0.807 (±0.277)	0.066 (±0.022)	60/115
fOTU_9761	s__Leohumicola sp	0.785***	EXT	0.123 (±0.058)	0.795 (±0.230)	76/108
fOTU_647	s__Clavaria sp	0.745***	EXT	0.163 (±0.065)	0.752 (±0.221)	74/108
fOTU_636	g__Clavaria	0.711***	EXT	0.156 (±0.046)	0.760 (±0.163)	66/108
fOTU_395	g__Leohumicola	0.705***	EXT	0.166 (±0.072)	0.749 (±0.224)	64/108
fOTU_662	s__Alternaria hordeicola	0.703***	EXT	0.141 (±0.030)	0.776 (±0.272)	63/108
(b) Bacteria
bOTU_1640	g_Nakamurella	0.906***	INT	0.775 (±0.078)	0.101 (±0.015)	107/115
bOTU_809	g_Clostridium sensu stricto 1	0.852***	INT	0.700 (±0.086)	0.181 (±0.062)	104/115
bOTU_3139	g_Dokdonella	0.790***	INT	0.741 (±0.089)	0.137 (±0.021)	87/115
bOTU_1833	f_Gemmatimonadaceae	0.776***	INT	0.686 (±0.088)	0.195 (±0.037)	89/115
bOTU_808	f_Gemmatimonadaceae	0.772***	INT	0.686 (±0.072)	0.196 (±0.042)	89/115
bOTU_26 514	g_Mycobacterium	0.770***	INT	0.709 (±0.083)	0.171 (±0.029)	85/115
bOTU_556	f_Pedosphaeraceae	0.761***	INT	0.696 (±0.096)	0.185 (±0.043)	85/115
bOTU_844	g_Pir4 lineage	0.761***	INT	0.678 (±0.063)	0.204 (±0.042)	87/115
bOTU_2462	g_Dechloromonas	0.758***	INT	0.771 (±0.137)	0.105 (±0.032)	76/115
bOTU_1680	f_A4b	0.756***	INT	0.686 (±0.081)	0.196 (±0.041)	85/115
bOTU_4498	f_Pirellulaceae	0.755***	INT	0.694 (±0.073)	0.187 (±0.031)	83/115
bOTU_18 978	g_Candidatus Udaeobacter	0.755***	INT	0.739 (±0.098)	0.139 (±0.038)	78/115
bOTU_752	c_Subgroup_6	0.732***	INT	0.668 (±0.087)	0.214 (±0.045)	80/115
bOTU_22 078	f_Pirellulaceae	0.730***	INT	0.714 (±0.080)	0.166 (±0.032)	76/115
bOTU_5565	f_A4b	0.728***	INT	0.756 (±0.133)	0.121 (±0.019)	71/115
bOTU_3115	o_R7C24	0.727***	INT	0.733 (±0.123)	0.146 (±0.025)	74/115
bOTU_1901	g_Pirellula	0.720***	INT	0.689 (±0.069)	0.192 (±0.041)	77/115
bOTU_1203	g_Pir4 lineage	0.720***	INT	0.710 (±0.081)	0.170 (±0.035)	75/115
bOTU_2079	o_SBR1031	0.718***	INT	0.699 (±0.081)	0.181 (±0.038)	74/115
bOTU_15 989	g_Arenimonas	0.716***	INT	0.655 (±0.081)	0.228 (±0.050)	79/115
bOTU_23 857	o_Corynebacteriales	0.715***	INT	0.831 (±0.154)	0.041 (±0.011)	62/115
bOTU_5957	f_BIrii41	0.714***	INT	0.686 (±0.089)	0.196 (±0.041)	78/115

DISCUSSION

Management had a strong influence on microbial community structure at the European and local scale

Contrary to our hypothesis 1, ‘MI’ had a similar effect size as ‘Country’ on fungal community structure across the European transect. Remarkably, this result indicates that small-scale, regionally localized management practices, which have previously been shown to effect soil fungal communities (de Vries et al. 2007, 2012b), have a comparable effect size on fungal community structure (mean centroid distance = 0.415) as biogeographically diverse, continental-scale geographic factors (mean centroid distance = 0.429). An important point to consider is that the INT and EXT grasslands in each region were sampled comparatively close to each other, while the sites in each of the regions and countries were sampled far, and extremely far, from each other, respectively. If such short distances (resulting in comparable environmental growth conditions) would influence soil fungal community structures, this would strongly reduce the magnitude of the management effect observed, compared to both region and country. Thus, the effect size of management on soil fungal community structure seen in this study is quite impressive. There was also a significant effect of ‘MI’ on soil bacterial community structure at the transect level, though this effect was not as strong a driver as ‘Country’.

Effects of individual management variables/ practices have been studied under controlled conditions. For example, N and P mineral fertilizer addition has been shown to influence fungal and bacterial community composition over large spatial scales under experimental field conditions (Leff et al. 2015). Our study shows that effects of grassland management on soil microbial community structure are also observed in ‘on-farm’ conditions at the European scale. These changes in the soil microbial community structure as a result of MI are most likely the result of dissimilarities in the soil environment, specifically nutrient availability, and differing plant community structure and diversity in contrastingly managed grasslands (Orwin et al. 2010; de Vries et al. 2012a). Under field conditions, it is not possible to separate these compounding influences as they always act concertedly, though both, and their interactions, may contribute to the observed effects. The strong effect size of MI on the soil microbiome was even more impressive given the fact we did sample the management intensity range of agricultural land and excluded nature conservation areas. This highlights the ecological value of the rare EXT grasslands in the agricultural landscape.

In each country, the factor ‘MI’ had a stronger impact than ‘GCR’ on fungal community structure, which indicates that as differences among geographic patterns decrease, MI becomes the dominant influencing factor, confirming the continental scale finding at the regional scale. This is especially impressive in Sweden and Switzerland, given the large differences in environmental variables between their FR (e.g. mean annual temperature = 7.48°C and 9.84°C, respectively) and LFR (3.15°C and 5.82°C, respectively). The strong effect of management was also confirmed by the fact that ‘MI’ had a significant effect on fungal and bacterial community structure in nine and eight of the 10 regions, respectively (Table 4 and Fig. 3). Here, the effect of ‘MI’ was tested separately in these regions, thus the environmental conditions of the other regions could not have had an effect. Indeed, the effect of ‘MI’ was so strong at the regional level, that significance was reached even with the greatly reduced number of sites.

The causes of the MI effect on soil microbial community structure, e.g. differences in soil nutrient availability and consequentially plant community structure may be especially relevant in Switzerland, as it likely reflects the implementation of an agroecological scheme called ‘ecological compensation areas’, of which all sampled EXT sites in Switzerland were part of. Introduced already in 1993, the scheme aims to promote and protect biodiversity in agricultural areas, of which EXT grasslands are the dominant type (BLW 1998). The scheme has been shown to promote and protect plant biodiversity within the agricultural landscape, particularly in mountain regions (Peter et al. 2008; Kampmann et al. 2012; Ravetto Enri et al. 2020), and to have positive effects on insects (Albrecht et al. 2010). The results presented here demonstrate that they additionally harbor distinct fungal and bacterial community structures. This is most likely on account of them never having been intensified, as even short-term additions of soil nutrients to grassland can influence both the plant and soil microbial community structure for decades afterwards (Spiegelberger et al. 2006).

Even though the sampling design was optimized to sample INT and EXT under comparable environmental conditions by sampling nearby plots within small regions, the effect of MI was accompanied by significant differences in environmental variables between INT and EXT in five regions of the transect, the LFR region of Germany, Switzerland and Portugal, as well as both GCR of Azores. This is a consequence of farmers choosing a management scheme dependent on the differing abiotic conditions within their holding (Peter et al. 2008; Gilhaus et al. 2017). For example, in mountain regions (i.e. LFR Germany and Switzerland), the sharp difference in the site inclination between management types is reflective of EXT typically being on valley slopes, while INT are located on more easily accessible plots where mechanization is easier (Kampmann et al. 2008). Even these small, indirect MI effects could influence soil microbial community structure in grasslands, though we suspect that the direct effects of MI (e.g. fertilizer application) are more important. All the above-mentioned regions had higher Ptot in INT, suggesting that MI did directly influence the observed shift in microbial community structure in these regions, as P fertilization has been shown to do (Ikoyi, Fowler and Schmalenberger 2018; Ikoyi et al. 2020). Moreover, in four of the five remaining regions (i.e. FR of Sweden, Switzerland and Portugal, LFR of Sweden) there was a significant effect of MI on both fungal and bacterial community structure, even without the aforementioned indirect effects.

Geographic factors comparably influence both fungal and bacterial community structure

The effect of the factor ‘Country’ was comparable for both fungal and bacterial community structure (i.e. mean centroid distance = 0.429 and 0.388; √CV = 0.247 and 0.249, respectively). The clear contrast was the effect size of ‘MI’ (as discussed above), which was much greater on fungal than bacterial community structure (i.e. mean centroid distance = 0.415 and 0.289; √CV = 0.233 and 0.174, respectively). This would indicate that soil fungal communities are more strongly influenced by the availability of soil nutrient resources (Bahram et al. 2018; Fox, Lüscher and Widmer 2020), and that the validity of hypothesis 1 depends on the component of the soil microbiome considered.

Large differences in precipitation, as present across our transect (590–1690 mm), has previously been shown to be a strong determinant of both fungal and bacterial community structure in a large scale study of Mongolian grasslands (where precipitation ranged 104–412 mm, Chen et al. 2015). Other climatic factors, particularly temperature, has been shown to have a strong influence on soil fungal and bacterial diversity of other land use types, such as forests (Zhou et al. 2016). Particularly notable in our transect was the strong separation of the microbiome structures of the southern countries, Portugal and Azores, both from each other and from that of the other countries. Both Portugal and Azores represent quite different biogeographic regions, Mediterranean and Macaronesian respectively, from the other countries of the transect (Alpine, Boreal and Continental). Additionally, the isolation of Azores from mainland Europe (∼1190 km from the coast of Portugal), could be contributing to the Azores–Portugal difference.

This continental scale study not only generated a range of climatic conditions, but also of soil physicochemical properties (i.e. % soil sand in Portugal) that may further differentiate the microbial community structure between countries. Soil physicochemical properties have been shown to be strong drivers of soil fungal and bacterial communities in large scale studies (Chen et al. 2015; Bahram et al. 2018). Along a European transect, which included most of the countries studied here, it was shown that spatial variations in bacterial community structure and diversity were principally driven by soil physicochemical properties (Plassart et al. 2019).

We included GCR within countries as part of our transect to further expand the range of sampled environmental conditions and geographic patterns. The GCR within Sweden (centroid distance 1.95) and Switzerland (centroid distance 2.46) were especially biogeographically divergent, with the factor ‘GCR’ significantly determining microbial community structure in these two countries. The geographic distance between the FR (Continental) and LFR (Boreal) in Sweden was the largest within country distance of the transect at ∼1200 km. This distance is comparable to, or greater than, the distance between different countries of the transect, e.g. Portugal–Azores. This latitudinal distinction in Sweden was previously shown to influence the fungal community structure in Scots pine needles (Millberg, Boberg and Stenlid 2015). In Switzerland, the FR (Continental, ∼487 m) and LFR (mountain, ∼1266 m) differ sharply in altitude, and consequently their climate, and have previously been shown to differ in their grassland plant communities (Güsewell, Peter and Birrer 2012). There were much smaller and statistically not significant differences in the climatic variables in the GCR regions of the southern countries of Portugal (centroid distance 0.76) and Azores (centroid distance 0.18), which may not have been large enough to shape the soil microbiome among GCR.

Strong regional specificity in differences of the soil microbial taxa to management intensity

Strong regional specificity was seen in the particular response of the microbiome (i.e. OTU richness, abundance of fungal and bacterial phyla, management-associated indicator OTUs) to grassland MI, which is contrary to the second hypothesis of this study. This is in line with the studies second hypothesis and reflects on-farm conditions and the fact that optimal grassland management is adapted to suit regional conditions (Gilhaus et al. 2017). Thus the INT management can consequently address completely different plant limiting conditions, as was the case between Switzerland (N fertilization) and Portugal (irrigation), for example. Furthermore, there could be different regional interactions of grassland MI with soil type and physicochemical characteristics (Xue et al. 2018), climate (Chen et al. 2015; Bahram et al. 2018) and plant community diversity and composition (Gilhaus et al. 2017), as all differ among countries.

In four regions of the transect fungal OTU richness was significantly higher in EXT. This finding supports previous studies, which indicated that EXT grasslands promote fungal dominated soil microbial communities (de Vries et al. 2007, 2012b). However, the opposite trend was observed in three other regions in the southern part of the transect (i.e. FR of Portugal and Azores, LFR of Azores) and may represent a different interaction effect between the southern regions and MI to that observed in Sweden, Germany and Switzerland. This, along with the effect of MI on soil fungal biomass, warrants further research. There was also a lack of consistency in the differences in the abundance of fungal and bacterial phyla, with no single phyla having a uniform response to MI in all regions. Though a small number did respond in multiple regions, such as the Firmicutes, whose tendency for higher relative abundance in INT grasslands may reflect slurry or manure application (Abubaker et al. 2013).

There were many (in most cases hundreds) of indicator fOTUs and bOTUs for MI detected in each region. The extent of the lack of sharing of these MI-associated fOTU and bOTU was surprising, even being low between regions within the same country. This would indicate that abiotic conditions are strong determinants for the specific taxonomic differences of the soil microbiome to MI (Ma et al. 2018). This finding has important implications for soil biodiversity monitoring. Due to the strong regional-specificity of microbial indicators to grassland MI, monitoring efforts should focus on identifying and protecting EXT grasslands harboring regionally unique/rare microbial taxa. In addition, future studies should consider the temporal development of soil microbial diversity in differently managed grasslands, both within a growing season and across years, as this is a rarely considered aspect in soil monitoring schemes (van Leeuwen et al. 2017).

Despite the strong regionally specific response, a number of microbial indicators for grassland MI were identified over the whole European transect. Indicator fungal genera for EXT grasslands included the Clavaria (fOTU_647 and fOTU_636, found in 74 and 66 of the 108 EXT sites, respectively) and Leohumicola (fOTU_9761 and fOTU_395, found in 76 and 64 of the 108 EXT sites, respectively). Clavaria is a saprotrophic genus associated with ‘unimproved’ European grasslands with low available sulfur and phosphate, and is highly sensitive to the application of nutrients (McHugh et al. 2001). The Leohumicola genus has been identified as an ericoid mycorrhiza. The symbiosis of ericoid mycorrhiza with the Ericaceae family of plants has been shown to help them adapt to acidic and/or nutrient-poor conditions (Cairney and Meharg 2003; Hambleton, Nickerson and Seifert 2005). Fungal indicator genera for INT included Thelebolus (fOTU_13, found in 87 of the 115 INT sites) and Podospora (fOTU_106, found in 78 of the 115 INT sites), both have been reported to develop on dung samples (Richardson 2001), and likely benefits from dung application. All of the identified bacterial taxa were indicators for INT. One indicator (bOTU_18978, found on 78 of the 115 INT sites) was assigned to the genus Candidatus Udaeobacter within the Verrucomicrobia phylum. A representative of this genus, Ca. Udaeobacter copiosus was identified as a highly ubiquitous taxa, particularly in grasslands (Brewer et al. 2016). It was reported as having a relatively small genome (∼40% smaller than the mean bacterial genome size) and may be sacrificing metabolic versatility as a competitive strategy. This may be a particular advantage in INT grasslands, as they have greater soil nutrient availability and lability, compared to EXT grasslands.

CONCLUSION

Our results demonstrate the surprisingly strong effect that small-scale grassland MI has on microbial community structure. For fungi, it even was a comparable driver to geographic factors at the continental scale. While geographic factors had a similar effect size on both fungal and bacterial community structure, the effect of ‘MI’ was far stronger on the former. Importantly, these results establish that EXT grasslands harbor distinct fungal and bacterial communities in all examined regions, across a divergent climatic gradient. Despite the MI effects on individual microbial taxa being strongly regionally specific, some microbial indicators for grassland MI were detected at the transect level (e.g. the fungal genera Clavaria and Leohumicola). This study demonstrates the ecological importance of the rare EXT grasslands as a distinct habitat for the soil microbiota, providing additional evidence of their capacity for biodiversity promotion in agroecological programs.

DATA ACCESSIBILITY

The values of each variable for each site considered in this study can be found in Tables S6 and S7 (Supporting Information) in the supplementary materials. DNA sequences were deposited in the NCBI SRA archive with the project number PRJNA641340. Supplementary data associated with this article can be found, in the online version, at https://doi.org/….

FUNDING

This research was funded through the 2015–2016 BiodivERsA COFUND call for research proposals, with the national funders being the Swedish Research Council for Environment, Agriculture and Spatial Planning (FORMAS, grant 2016–01927), the German Research Foundation (DFG, grant RA1717/4–1), the Swiss National Science Foundation (SNF, grant 31BD30-172463) and the Governo Regional dos Açores through Fundo Regional para a Ciência e Tecnologia (FRCT, grant M1.1.C/2/001/2017).

ETHIC APPROVAL

The work described here did not require ethics approval.

Click here for additional data file.