Land-use intensification differentially affects bacterial, fungal and protist communities and decreases microbiome network complexity.

BACKGROUND: Soil microbial communities are major drivers of cycling of soil nutrients that sustain plant growth and productivity. Yet, a holistic understanding of the impact of land-use intensification on the soil microbiome is still poorly understood. Here, we used a field experiment to investigate the long-term consequences of changes in land-use intensity based on cropping frequency (continuous cropping, alternating cropping with a temporary grassland, perennial grassland) on bacterial, protist and fungal communities as well as on their co-occurrence networks.
RESULTS: We showed that land use has a major impact on the structure and composition of bacterial, protist and fungal communities. Grassland and arable cropping differed markedly with many taxa differentiating between both land use types. The smallest differences in the microbiome were observed between temporary grassland and continuous cropping, which suggests lasting effects of the cropping system preceding the temporary grasslands. Land-use intensity also affected the bacterial co-occurrence networks with increased complexity in the perennial grassland comparing to the other land-use systems. Similarly, co-occurrence networks within microbial groups showed a higher connectivity in the perennial grasslands. Protists, particularly Rhizaria, dominated in soil microbial associations, as they showed a higher number of connections than bacteria and fungi in all land uses.
CONCLUSIONS: Our findings provide evidence of legacy effects of prior land use on the composition of the soil microbiome. Whatever the land use, network analyses highlighted the importance of protists as a key element of the soil microbiome that should be considered in future work. Altogether, this work provides a holistic perspective of the differential responses of various microbial groups and of their associations to agricultural intensification.

Background

Agricultural intensification represents one of the major causes of biodiversity loss in the twenty-first century [1, 2]. Intensive land use also threatens ecosystem multifunctionality [3] and increases negative environmental impact, for example by increased greenhouse gas emissions [4] and nitrogen losses through leaching [5]. Land intensification is most often associated with shifts in vegetation, which in return can affect the soil microbiome. It is well known that above-ground vegetation influences below-ground microorganisms through, for example, the exudation of organic carbon or the modification of soil properties in close vicinity to plant roots [6-8]. Besides, intensive agricultural practices such as tillage, fertilization, and the use of pesticides can affect the soil microbiome [9-11]. Not only current but also past land-use can have effects that persist for decades, affecting the microbiome in contemporary land use [12, 13]. For example, a history of tillage or no tillage 60 years before conversion to longleaf pine savanna was found to differentially affect the diversity and composition of bacterial and fungal communities [13]. Since soil microbes are major players of the biogeochemical cycles and therefore tremendously important for ecosystem functions [14-16], it is important to understand how microbial communities are affected by increased land-use intensity, and to what extent the responses are modulated by legacy effects of prior land use.

Although the soil microbiome is highly diverse and complex, most studies investigating the impact of land-use intensification on microbial diversity have focused on either bacteria and/or fungi [17-20]. However, microbes are not a mere collection of independent populations, but form complex communities and the direct and indirect interactions between taxa appear to play a crucial role for the assembly of microbial communities, with consequences for ecosystem functioning [21-25]. For example, the often neglected protists are important in agricultural ecosystems as they control microbial populations through predation [26-28] and influence ecosystem functions e.g. nutrient cycling [29]. In particular, by grazing on bacteria, phagotrophic protists are increasing the mobilization of bacterial nitrogen and are therefore contributing to nitrogen mineralization [29, 30]. This exemplifies the close interaction across soil microbes from different domains and that such interactions can affect soil functioning. An integrated assessment of the soil microbiome is thus needed for a comprehensive understanding of the impact of land-use intensification on soil microorganisms. Co-occurrence networks are often used to decipher relationships between microorganisms [31]. Although co-occurrence patterns might not reflect the true complexity of microbial interactions, this integrative approach can help to better understand the consequences of changes in land use on microbial communities [31-33]. For example, fungal networks were more connected in organic farming than in conventional or no-till farming systems, which was concomitant with a decrease in root colonization by arbuscular mycorrhizal fungi [17]. However, studies investigating how land-use intensification modifies co-occurrence patterns within and especially across domains are still scarce.

Here, we investigated how land-use intensification and the legacy effects of prior land-use affect the diversity, structure and co-occurrence patterns within and between bacterial, protist and fungal communities. For this purpose, we used a long-term (> 12 years) replicated field experiment, where a gradient of land-use intensity was established, with high (continuous cropping with annual crops), medium (alternating annual crops with a temporary grassland) and low (perennial grassland) intensity. We hypothesized that (i) changes in land-use will lead to distinct responses across microbial groups with smallest differences between temporary and perennial grasslands if there are no legacy effects (i.e. no impact of the continuous cropping preceding the temporary grassland), and (ii) co-occurrence patterns of the soil microbiome will shift along the gradient of land-use intensity with a higher network complexity in the perennial systems since the annual rotation of short lived-seasonal crops is less likely to support the co-evolution of multitrophic interactions within the soil microbiota. Because soil microbes are key players in the cycling of nitrogen, which is the major nutrient limiting primary production in terrestrial ecosystems [34], we also investigated whether land-use intensification affected the abundances of functional microbial guilds involved in N-cycling.

Methods

Experimental design and soil sampling

Soil samples were collected from a long-term field experiment SOERE-ACBB (Systems of Observation and Experimentation in Environmental Research in Agro-ecosystems, Biochemical cycles and Biodiversity) at the INRAE experimental station located at Lusignan, France (46°25′12.91″ N; 0°07′29.35″ E). The soil is a Cambisol with loamy texture (105 g kg−1 sand, 727 g kg−1 silt and 168 g kg−1 clay). Organic carbon content and nitrogen levels are shown in Additional file 1 (table S1). The experiment was established in 2005 in a randomized complete block design divided in four blocks containing plots with a surface of 4000 m2 each (Additional file 1: Figure S1). Each block comprises the following treatments: a three-year crop rotation of maize-wheat-barley (CC), a three-year temporary grassland alternated with the three-year crop rotation (TG) and a permanent grassland (PG). In the crop rotation, conventional tillage (mouldboard ploughing) and fertilizer management was applied, whereas in the grasslands, the timing and rate of nitrogen application was guided by a nitrogen nutrition index between 0.9 and 1.0 to provide the lowest nitrogen amount for potential plant production [35, 36]. The dominant vegetation of the temporary and permanent grassland was composed of Lolium perenne, Festuca arundinacea and Dactylis glomerata. The land-use intensity was defined by cropping frequency (continuous cropping > alternating cropping with a temporary grassland > perennial grassland).

Sampling was carried out in November 2017, which was before sowing for CC and TG, and therefore after the last year of the temporary grassland for TG. Five soil cores were collected randomly from each replicated plot using a soil corer of 20 cm depth and 5 cm diameter to obtain 20 samples from each treatment (N total = 60 samples). Plant debris and roots were removed from soil samples, which were then homogenized and sieved at 4 mm. Soil relative moisture was determined by drying 10 g of fresh soil at 105 °C for 24 h.

DNA was extracted from the 60 soil samples using the DNeasy PowerSoil-htp 96 well DNA isolation kit (Qiagen, France). Analysis of the diversity and composition of the bacterial, protist and fungal communities were performed for all DNA extracts using Illumina sequencing of the 16S rRNA and 18S rRNA genes and the fungal ITS1 region. Amplicons were generated for all DNA extracts in two steps. In the first step, the V3-V4 hypervariable region of the bacterial 16S rRNA gene was amplified by polymerase chain reaction (PCR) using the fusion primers U341F (5’-CCTACGGGRSGCAGCAG-3’) and 805R (5’-GACTACCAGGGTATCTAAT-3’)[37]. The V4 hypervariable region of 18S rRNA gene was amplified by PCR using the primers EK-565F (5′-GCAGTTAAAAAGCTCGTAGT-3′) and 18S-EUK-1134-R–UnonMet (5′-TTTAAGTTTCAGCCTTGCG-5′) [38, 39]. Fungal ITS1 region was amplified using the primers ITS1F (CTTGGTCATTTAGAGGAAGTAA) and ITS2 (GCTGCGTTCTTCATCGATGC) primers [40, 41]. Sequencing was performed on MiSeq (Illumina, 2 × 250 bp and 2 × 300 bp) using the MiSeq reagent kit v2 and v3 (500 and 600 cycles, respectively).

Sequence data from the 60 soil samples was analyzed using an in-house developed Python pipeline (https://forgemia.inra.fr/vasa/illuminametabarcoding). Briefly, 16S rRNA gene, 18S rRNA gene and ITS sequences were assembled using PEAR [42] with default settings. Further quality checks were conducted using the QIIME 1 pipeline [43] and short sequences were removed (< 400 bp for 16S, < 450 bp for 18S and < 300 bp for ITS). Reference based and de novo chimera detection, as well as OTU clustering were performed using VSEARCH [44] and the adequate reference databases (SILVA representative set of sequences for 16S and 18S rRNA, and UNITE’s reference dynamic dataset for ITS). The identity thresholds were set at 94% for 16S rRNA gene-based on replicate sequencing of a bacterial mock community that contains 40 different bacterial strains for which we have the full-length 16S rDNA sequences [45] and 97% for 18S rRNA gene and ITS. Representative sequences for each OTU were aligned using PyNAST [46] and MUSCLE for 16S rRNA and 18S rRNA, respectively. Phylogenetic trees were constructed using FastTree [47]. For 16S, taxonomy was assigned using UCLUST [48] and the SILVA reference database 132 [49]. For 18S, taxonomy was assigned using PR2 database 4.11.1 [50]. For ITS, the taxonomy assignment was performed using BLAST [51] and the UNITE reference database (v.7–08/2016 [52]). Raw sequences were deposited at the NCBI under the BioProjects PRJNA741976, PRJNA741982 and PRJNA742156.

Based on taxonomic assignments, we filtered out OTUs from the 18S rRNA gene sequences that were non-protist (i.e. OTUs belonging to Fungi, Streptophyta, Mollusca and Porifera) and OTUs from ITS sequences that were non-fungal (i.e. OTUs belonging to Cercozoa). In total, 2 319 835 bacteria sequences, 1 632 206 protist sequences and 3 244 474 fungal sequences were obtained and assigned to 6 999, 3 867 and 1 866 OTUs, respectively. Bacterial, protist and fungal α-diversity metrics (i.e. observed species, Simpson’s reciprocal, Shannon as well as Faith’s Phylogenetic Diversity PD for bacteria and protist [53]) were calculated based on rarefied OTU tables (23,000 sequences per sample for 16S rRNA, 9000 sequences per sample for 18S rRNA and 35,000 sequences per sample for ITS). Bray–Curtis dissimilarity matrix were also computed to detect variations in the structure of microbial communities.

The abundances of total bacterial and fungal microbial communities as well as that of nitrogen cycle microbial guilds were estimated by real-time quantitative PCR (qPCR) assays. For each land use, the five DNA extracts from each replicated plot were pooled in equimolar amounts which corresponded to 4 composite samples per land use used as templates for the qPCR assays (n = 4). Total bacterial and fungal communities were quantified using 16S rRNA and ITS primers as described by Muyzer et al. (1993) [54] and White et al. (1990) [41], respectively. Marker genes for nitrification (amoA in archaea and bacteria) and denitrification (nirK and nirS) were quantified as described previously [55, 56]. Finally, nifH and nrfA genes were used to quantify communities involved in nitrogen-fixation [57] and dissimilatory nitrate reduction to ammonium (DNRA) [58, 59], respectively. The qPCR reactions were carried out using a ViiA7 (Life Technologies, USA) in 15 µl reactions containing 7.5 µL of Takyon Master Mix (Eurogentec, France), 1 µM of each primer, 250 ng of T4 gene 32 (QBiogene, France) and 1 ng of DNA. Two independent runs were performed for each real time PCR assay. Standard curves were obtained using serial dilutions of linearized plasmids containing appropriated cloned targeted genes from bacterial strains or environmental clones. PCR efficiency for the different assays ranged from 77 to 101%. No template controls gave null or negligible values. Inhibition in qPCR assays was tested by mixing soil DNA extracts with either control plasmid DNA (pGEM-T Easy Vector, Promega, France) or water. No inhibition was detected with the amount of DNA used.

Statistical analyses were conducted using R statistical software version 3.4.1 [60]. Differences in gene copy number (16S rRNA, ITS, bacterial and archaeal amoA, nirK, nirS, nifH and nrfA) were tested using the Kruskal–Wallis test followed by Dunn’s multiple comparison test (adjusted p-value < 0.05). Differences in the microbial α-diversity indices were tested using ANOVA followed by Tukey's honestly significant difference (HSD) test (p-value < 0.05) using the agricolae package [61]. Normality and homogeneity of the distribution of residuals were verified and log-transformations were performed when necessary. Permutational multivariate analysis of variance (PermANOVA) was carried out on the Bray–Curtis dissimilarity distance matrices using adonis function implemented in the vegan package [62]. Pairwise post hoc tests were conducted using the function pairwise.adonis from the pairwiseAdonis package [63] with Holm corrections. As sequencing data are usually sparse, low-abundance OTUs were filtered out by keeping OTUs that are present at a threshold of 0.04% in all samples. We also discarded the OTUs that were not found in at least 10 out of the 60 samples. This filtering step allows reducing the zero counts in sequencing datasets, which can inflate the number of false positive for the differential abundances analysis and spurious correlation between OTUs in network analysis. This resulted in 472 OTUs for 16S rRNA, 341 OTUs 18S rRNA and 218 for ITS. Differential abundance analysis of microbial community composition was conducted by pairwise comparisons between land uses of filtered count matrices (n = 20) using the DESeq2 package (FDR-corrected p-value < 0.00001, Additional file 2) [64]. Significantly discriminant OTUs were illustrated by ternary plots using ggtern package [65] and Upset plot using UpsetR package [66].

Networks were inferred using a sparse multivariate Poisson log-normal (PLN) model with a latent Gaussian layer and an observed Poisson layer using the PLNmodels package (Additional file 2) [67]. The best network was selected using a Stability Approach to Regularization Selection (StARS) [68], which performs a random subsampling of the input data to select a network with low variability in the selected edges. A specific normalization with the TSS (Total Sum Scaling) method was performed in order to take into account the heterogeneity of sequencing depth within and between microbial groups. All networks were constructed using filtered count matrices (i.e. 472 OTUs for 16S rRNA, 341 OTUs 18S rRNA and 218 for ITS). Bacterial, protist and fungal networks were inferred separately for each land use (n = 20) and for visualization purpose only partial correlations with |ρ|> 0.1 were considered. Inter-domain networks were inferred using all microbial groups for each land use (n = 20), and for visualization purpose partial correlations with |ρ|> 0.08 were visualized. Networks were then visualized using the Cytoscape software [69]. The NetworkAnalyser tool from Cytoscape was used to calculate network topological parameters (i.e. nodes, links, clustering coefficient and degree).

Results

To characterize the soil bacterial, fungal and protist communities under continuous cropping (CC), alternating cropping with a temporary grassland (TG) and perennial grassland (PG), we used DNA metabarcoding targeting the bacterial, protist and fungal communities. Significant differences in α-diversity indices were observed for the bacterial community, with OTU richness, Shannon index, and phylogenetic diversity being higher in continuous cropping and temporary grassland than in the perennial grassland (Tukey’s test, p-value < 0.05, Additional file 1: Figure S2). By contrast, the α-diversity indices of the protist and fungal communities were similar across the three land-uses.

Comparison of β-diversity using Principal Coordinates Analysis (PCoA) of Bray–Curtis distances showed a clustering of samples according to management with 37%, 23% and 34% of the variance explained by the first two axes of the PCoA for bacterial, protist and fungal communities, respectively (Fig. 1). PermANOVA confirmed significant differences in the structure of the microbial communities between land uses (PermANOVA, p-value < 0.05). However, the structure of the bacterial and protist communities were more similar between continuous cropping and the temporary grassland than the perennial grassland (Additional file 1: Figure S3).

Fig. 1

Principal Coordinates Analysis (PCoA) of the Bray–Curtis distance matrices of a bacterial 16S rRNA gene fragments, b protist 18S rRNA gene fragments and (c) fungal ITS1 amplicons showing shifts in community structure between continuous cropping (CC), temporary grassland (TG) and perennial grassland (PG) land-use

The abundance of the total bacterial community, determined by quantification of the 16S rRNA gene, was significantly lower under continuous cropping compared to temporary and perennial grasslands (Dunn’s test, adjusted p-value < 0.05, Fig. 2). For the N-cycling guilds, the percentage of DNRA bacteria (nrfA) within the total bacterial community (nrfA/16S rRNA gene abundance) increased along the land-use intensity gradient (Dunn’s test, adjusted p-value < 0.05, Fig. 2). The percentage of AOB and of nirS-denitrifiers were also higher in the continuous cropping, but significant differences were not observed for other N-cycling guilds (Dunn’s test, adjusted p-value < 0.05).

Fig. 2

Abundances of a total bacteria (16S rRNA gene) and b fungal (ITS) in continuous cropping (CC), temporary grassland (TG) and perennial grassland (PG) soil samples (mean ± s.e.). The percentage of c ammonia-oxidizing archaea (AOA) to ammonia-oxidizing bacteria (AOB). Percentage of d AOB, e and f bacterial denitrifiers (nirK and nirS), g nitrogen-fixing communities (nifH) and h DNRA communities (nrfA) in the total bacterial community. Different letters above the bars indicate significant differences between land uses according to Dunn’s test (adjusted p-value < 0.05)

Across all samples, the dominant bacterial phyla were Proteobacteria (33%), Acidobacteria (21%) and Actinobacteria (17%). Protists were dominated by Cercozoa (45%), Chlorophyta (15%) and Stramenopiles (14%), whereas Mortierellomycotina (42%) and Sordariomycetes (22%) represented the most abundant fungal taxa (Additional file 1: Figure S4). The ternary plots showed the distribution of the most abundant OTUs, with those present at similar abundances in all samples being near the center of the ternary plots (Additional file 1: Figure S5). Differential abundance analysis based on the most abundant OTUs (472 OTUs for 16S rRNA, 341 OTUs 18S rRNA and 218 for ITS), identified 173 bacterial OTUs with significant changes in abundances between land-uses (FDR-corrected P < 0.00001) (Fig. 3a). The strongest differences were observed between continuous cropping and permanent grasslands, with 166 discriminant bacterial OTUs and 101 of them being significantly more abundant under perennial grasslands than continuous cropping (Fig. 3b). By contrast, only 28 OTUs differed significantly between continuous cropping and temporary grasslands, with 27 OTUs increasing in the latter. OTUs belonging to Actinobacteria were the most impacted by land-use intensity, showing a decrease in relative abundances from 7 to 3% from the perennial grassland to the continuous cropping (Fig. 3a and Additional file 1: Figure S4 a). For the protist, 79 out of the 341 dominant OTUs differed significantly between the three land uses (Fig. 3c). The majority were differentially abundant between perennial grassland and continuous cropping (Fig. 3d), with a decrease in the relative abundances of 55 OTUs from 32 to 10%, respectively. The most impacted OTUs belonged to Apicomplexa and Cercozoa (Fig. 3 c and Additional file 1: Figure S4b). For the fungal community, 95 out of 218 fungal OTUs were significantly impacted by the land management (Fig. 3e). Similar to the patterns observed for bacteria and protist, the largest differences were between the perennial grassland and the continuous cropping with 42 OTUs being more abundant under perennial grassland (e.g. Agaricomycetes and Mortierellomycotina), and 48 OTUs mainly assigned to Dothideomycetes having higher abundances under continuous cropping (Fig. 3f and Additional file 1: Figure S4c). We found that only 4 OTUs belonging to fungi and protist were significantly increasing along the land-use intensity gradient, (CC > TG > PG), while 12 OTUs that belonged to all three domains were decreasing (CC < TG < PG).

Fig. 3

Ternary plots representing the composition of the microbial community under different land uses. Differential OTU abundances between land uses identified by DESeq2 are represented for a bacterial, c protist and e fungal communities. The position of each circle on the axis represents the contribution of the indicated land use to the relative abundance of each OTU. The size of the circle indicates the mean frequencies of each OTU in all samples. The colors indicate the affiliation of OTUs at the phylum or class levels. Upset plots showing the differentially abundant OTUs between land uses for b bacterial, d protist and f fungal communities. Each vertical bar shows the number of differentially abundant OTUs that are shared or unique between land use comparisons

We inferred microbial association networks from continuous cropping, temporary and perennial grassland soil samples for the bacterial, protist and fungal communities using a recently developed sparse multivariate Poisson log-normal model [67]. The complexity of networks was estimated using network size (i.e. the number of nodes), links, the positive to negative links ratio in the networks, the clustering coefficient (i.e. the degree to which nodes are clustered) and the average degree (i.e. average links per node) (Table 1). In the bacterial networks, both the number of nodes and links gradually decreased along the land-use intensity gradient (Fig. 4a, Table 1). However, the highest ratio of positive to negative links in the temporary grassland was more than three times higher than that in the other land uses (CC = 2.33, TG = 7.45 and PG = 2.22), due to a decrease in negative associations in the temporary grassland with 13% of negative links comparing to 42% and 45% for continuous cropping and perennial grassland, respectively. For protists, the ratio of positive to negative links was also higher in the temporary grassland, whereas the other network properties were similar between land uses (i.e. number of nodes, links and average degree) (Fig. 4b and Table 1). The fungal networks showed a different trend with a higher complexity in the temporary grassland (Fig. 4c and Table 1). The fungal network also revealed an increase in the ratio of positive to negative links in the temporary grassland (CC = 3.20, TG = 5.21 and PG = 4.50), similar to bacterial and protist networks. Networks from the temporary grassland of the bacterial, protist and fungal communities exhibited the highest clustering coefficient.

Table 1

Properties of microbial co-occurrence networks

Microbial group	Bacteria			Protist			Fungi			All
Land use	CC	TG	PG	CC	TG	PG	CC	TG	PG	CC	TG	PG
Total number of OTUs	472	472	472	341	341	341	218	218	218	1031	1031	1031
Nodes	65	85	137	210	208	207	128	132	117	445	517	586
Links	60	93	161	326	330	315	168	211	143	778	971	985
Positive links	42	82	111	182	217	189	128	177	117	519	700	676
Negative links	18	11	50	144	113	126	40	34	26	259	271	309
Ratio ±	2.33	7.45	2.22	1.26	1.92	1.50	3.20	5.21	4.50	2.00	2.58	2.19
Avg. clustering coefficient	0.032	0.105	0.067	0.020	0.037	0.024	0.047	0.094	0.057	0.039	0.049	0.031
Avg. degree	1.85	2.19	2.35	3.10	3.17	3.04	2.63	3.20	2.44	3.50	3.76	3.36

Fig. 4

Co-occurrence networks of a bacterial, b protist and c fungal communities in continuous cropping, temporary grassland and perennial grassland. Nodes are colored according to their taxonomic affiliation at phylum and class levels. The size of the nodes is proportional to the number of links per node (i.e. degree). Link thickness is proportional to partial correlations between nodes and represents associative (black, ρ > 0.1) or exclusionary relationships (red, ρ < − 0.1)

To better understand how land use influences associations between bacteria, protist and fungi, we inferred co-occurrence networks including all three groups for each land use (Fig. 5a). The microbial network with all groups from the perennial grassland was the most complex with 586 nodes and 985 links, followed by the temporary grassland with 517 nodes and 971 links and continuous cropping network with 445 nodes and 778 links (Fig. 5b). By contrast, we found a higher clustering coefficient and average degree in the temporary grassland network indicating that the nodes were more connected than those of continuous cropping and perennial grassland (Table 1). The comparison between inter-group networks showed that 46.3% of nodes are shared between land uses. In contrast, only two microbial links (fungi-fungi and protist-protist) were shared between to the networks from the three land uses (Fig. 5d).

Fig. 5

a Co-occurrence networks across three groups in continuous cropping (CC), temporary grassland (TG) and perennial grassland (PG). Nodes are colored according to their taxonomic affiliation at phylum levels. The size of the nodes is proportional to the number of links per node (i.e. degree). Link thickness is proportional to partial correlations between nodes and represents associative (black, ρ > 0.08) or exclusionary relationships (red, ρ < -0.08). b Number of nodes per land-use for each microbial group. c Number of positive (black) and negative (red) links in inter-domain networks for each land-use type and within microbial domains. d The Venn Diagrams show the number of shared/unique nodes and links across inter-domains networks

To distinguish differences in taxa co-occurring among land uses, we compared the number of positive and negative links within and between microbial groups (Fig. 5c). Regardless of land use, microbial networks were dominated by fungi-protist (28%), protist-protist (27%) and fungi-fungi (24%) associations. Bacteria-protist and bacteria-bacteria associations represented only 11 and 3% of the total number of links, respectively. However, the positive to negative link ratio for the bacteria-bacteria associations was the highest in temporary grassland (CC = 16, TG = 29 and PG = 3.3). Associations across all three groups were mainly characterized by a higher number of links in temporary and perennial grasslands. This was the case for bacteria-Rhizaria, fungi-Rhizaria, fungi-Alveolata and fungi-Stramenopiles associations (Additional file 1: Figure S6).

Discussion

Investigations of belowground effects of land-use changes have often focused either on the bacterial or on the fungal community. However, as the main consumers of bacteria, protists influence the composition of the soil microbiome [26, 70]. This highlights the need to adopt more comprehensive approaches for holistic insights into the soil microbiome [24]. Here, we performed comparative and integrative analyses of the soil bacterial, fungal and protistan communities across a gradient of land-use intensity based on cropping frequency. We showed that β-diversity of all three domains was affected by land-use type, with significantly different communities in the continuous cropping, temporary grassland and perennial grassland. Little is known about the importance of land-use intensity in driving soil protists and there are large discrepancies between studies [27]. Although climatic factors have been described as the top predictor of composition of protist communities [71], our result indicate that protists can be affected by changes in land use to the same extent as bacterial and fungal communities. However, the α-diversity of the protist and fungal communities showed no significant differences across land uses, although the bacterial α-diversity did, with the highest diversity in continuous cropping and temporary grassland compared to perennial grassland. In line with our results, previous studies found that shifts to agricultural land-use increased bacterial diversity [72, 73]. This might be explained by the rotation of plant communities in these systems that likely led to a greater diversity in soil nutrients through root exudates and crop residues, ultimately increasing bacterial diversity [74]. This is congruent with the observed differences in the biochemical nature of the soil organic matter (SOM) between the long-term grasslands and continuous cropping systems studied here [75]. The relative abundances of N-cycling communities (i.e. nitrifiers, denitrifiers and DNRA bacteria) were the highest in the continuous cropping, suggesting that land-use can be a driver for these functional groups. Accordingly, previous studies reported significant effects of land use on the relative abundances of AOB and DNRA bacteria [55, 76]. However, since microbial communities involved in N-cycling processes resulting in nitrogen losses (denitrification and nitrification) as well as in nitrogen retention (DNRA) were enriched in the crop rotation, further work is needed to determine the fate of inorganic N in these systems.

To determine potential legacy effects of prior land use (i.e. continuous cropping followed by three years of temporary grassland), the soil samples were collected in all land uses when the cycle of the temporary grassland ended (TG). Without any effect of land-use history, it was expected that microbiomes would be more similar under temporary and perennial grasslands than under continuous cropping. Instead, the smallest differences in community structure for all domains were observed between temporary grassland and continuous cropping, which suggests lasting effects of the cropping system preceding the temporary grasslands. The impact of land-use history on microbial communities has been investigated in previous studies with large variation in the persistence of the legacy effects [77, 78]. For example, past arable farming land resulted in long-lasting legacies on forest microbial communities that were persisting over half a century after agricultural abandonment [13, 79]. An explanation for these long-lasting effects is that historical land use caused shifts in soil properties, such as C and N content or pH, which are important drivers of soil microbial communities and may require decades to recover. Consistent with this and our results, Crème et al. (2018) [75] observed that despite the insertion of three years of temporary grasslands within the crop rotation, the biogeochemical signature of the soil organic matter was more similar to the continuous cropping system than to the permanent grassland. Moreover, plant species can also affect differently soil microbial communities through a myriad of processes [8, 80] and create legacies that are detectable under subsequent plant communities [81], therefore also explaining the similarities between temporary grassland and continuous cropping [82]. Overall, we found that prior cropping systems that were at least three-year-old continued to play a role in shaping microbial community composition across groups, which highlight the importance to consider legacies of past land-use when assessing the impact of management practices. This is all the more important considering that such changes in soil communities can in turn affect subsequent plant communities [83].

Due to intrinsic differences in eco-evolutionary dynamics between perennial and annual systems, it can be hypothesized that plant-soil feedbacks can lead to more complex and more connected networks in permanent grassland. However, when inferring microbial co-occurrence networks for each group (i.e. bacteria, protist and fungi), only the bacterial network was more complex in the perennial grassland with a gradual increase in the number of nodes and links between bacterial OTUs along with decreasing land-use intensity (i.e. from CC to PG). This is in agreement with a previous study showing higher bacterial network complexity in grassland than cropping system [84]. Interestingly, perennial grasslands also had the lowest bacterial diversity, which suggests that changes in soil microbial diversity do not necessarily reflect changes in microbial networks [85]. In contrast to a previous study showing that agricultural intensification reduces fungal network complexity [17], we found a higher complexity in the temporary grassland compared to the perennial and continuous cropping systems. This discrepancy may be attributed to the fact that root-associated fungi were monitored by Banerjee et al. [17] whereas the present study focuses on soil communities. Differences in the structure and composition of the protist community across land-use types were not mirrored in the network complexity. However, when considering the link type, we found an increase in the positive to negative link ratio in the temporary grasslands for protists as well as the other two groups, indicating that positive associations between soil microorganisms were promoted. This pattern might be due to increased niche differentiation when alternating temporary grassland and crop rotation, which results in a decrease of competition between microbial species [86, 87]. On the opposite, the decrease in the positive to negative link ratio in the cropping systems suggests less effective establishment of cooperation rather than increased competition between microbial taxa. Alternatively, previous work showed that differences in community evenness are likely to also affect the positive edge percentage since less prevalent taxa tend to contribute more to negative edges [32]. In any case, further studies are required to decipher whether these changes in co-occurrence networks were induced by biotic interactions or differences in environmental filters among the three land-use types.

Finally, we inferred networks encompassing all three groups to obtain an integrated and holistic view of the soil microbiome along the land-use intensity gradient. Similar to the bacterial network, we found a higher complexity of the inter-group network under permanent grasslands, which supports our hypothesis that perennial systems lead to more connected microbial networks. This increase in complexity was mostly due to an increase in the number of bacterial and protist nodes as well as bacteria-bacteria, protist-protist and bacteria-protist links, despite that the individual protist networks were similar between land uses in contrast to those of the bacterial. While almost half of the nodes of the inter-domain networks were similar, links were distinct between land-use types, indicating that microbial associations are more sensitive than community structure to land-use intensification. That land-use management affected associations within the soil microbiome is important to consider given that recent studies have shown that ecosystem functioning could be related to microbial network complexity [33, 88]. Whatever the land-use, protists were dominant in the inter-group network (Fig. 5b), especially the Rhizaria (Cercozoa), which were highly connected to both bacteria and fungi. Although co-occurrence networks suffer of spurious correlations when the effects of habitat filtering are strong, they can also recapitulate possible interactions between microorganisms under certain conditions [89, 90]. While not all protists feed on other organisms, they are known as the main consumers of bacteria [29]. However, we found more negative associations between protists (Rhizaria but also Amoebozoa and Stramenopiles) and fungi than between protists and bacteria. Accordingly, recent work suggest that protist feeding on fungi might equally be important [91]. On the other hand, some fungi have developed trapping structures, such as adhesive spores, hyphae to capture soil-inhabiting microorganisms such as protists [92]. Fungi can also be protist parasites, when for example amoeba ingests the spores or conidia, resulting in its death [93]. These predatory/parasitic interactions could explain some of negative associations between fungi and protists observed in our study.

Conclusions

Using a holistic microbiome investigation of bacterial, fungal and protist communities in a long-term field experiment managed under different levels of land use intensity, we showed that land management only affected α-diversity of the bacterial community, with increased diversity in the continuous cropping system. However, we identified a clear shift in the structure and the composition of all communities in response to land use, in particular between the continuous cropping and the perennial grassland. Moreover, our results showed legacy effects of cropping on the structure of the soil microbiome that lasted after three years of temporary grasslands. This highlights that prior land use can shape the present-day community for multiple microbial groups across domains. The perennial grassland system led to more complex bacterial as well as inter-domain networks, which can have implication for the contribution of microbes to ecosystem multifunctionality [16]. Inter-domain networks also revealed the predominant role of the protist as key taxa in soil microbiome networks across all land-use types. Future work need to validate the importance of protists in shaping soil microbial communities, directly through biotic interactions and/or indirectly through changes in abiotic factors.

Additional file 1: Table S1. Total carbon and nitrogen in continuous cropping (CC), temporary grassland (TG) and perennial grassland (PG). Figure S1. Experimental design of the long-term observatory near Lusignan, France. (a) The experiment is set up as a randomized complete block design divided in four blocks comprising (b) a three year crop rotation of maize-wheat-barley (CC), a three-year temporary grassland alternated with the three-year crop rotation (TG) and a permanent grassland (PG). Figure S2. α-diversity of the microbial communities in continuous cropping (CC), temporary grassland (TG) and perennial grassland (PG). (a) Simpson’s reciprocal, (b) Shannon, (c) observed species and (d) Faith’s phylogenetic diversity (PD) indices are shown. Different letters above the boxplots indicate significant differences according to Tukey’s test (p-value < 0.05). Figure S3. Bray–Curtis distances between land uses (mean ± s.d.). Different letters above the bars indicate significant differences according to Tukey’s test (p-value < 0.05). Figure S4. Microbial composition in continuous cropping (CC), temporary grassland (TG) and perennial grassland (PG) of bacterial (a), protist (b) and fungal (c) communities. Relative abundances are shown at the phylum and class levels and expressed as a percentage of the total number of OTUs. Figure S5. Ternary plots representing the composition of the microbial community under different land uses. The distributions of the most abundant OTUs in the bacterial (a), protist (b) and fungal (c) communities are shown. The position of each circle on the axis represents the contribution of the indicated land use to the relative abundance of each OTU. The size of the circle indicates the mean frequencies of each OTU in all samples. The colors indicate the affiliation of OTUs at the phylum or class levels. Figure S6. Number of positive (black) and negative (red) links between bacteria (B) or fungi (F) and groups of protists.

Additional file 2: Script 1. R code for differential abundance analysis using DESeq2 package. Script 2. R code of network analysis using a sparse multivariate Poisson log-normal (PLN) model.