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. Supplementary Information The online version contains supplementary material available at 10.1186/s40793-021-00396-9.


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][7][8]. Besides, intensive agricultural practices such as tillage, fertilization, and the use of pesticides can affect the soil microbiome [9][10][11]. Not only current but also past land-use can have effects that persist for decades, Romdhane et al. Environmental Microbiome (2022) 17:1 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][15][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 landuse intensification on microbial diversity have focused on either bacteria and/or fungi [17][18][19][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][22][23][24][25]. For example, the often neglected protists are important in agricultural ecosystems as they control microbial populations through predation [26][27][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][32][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 livedseasonal 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.

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 Agroecosystems, 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 m 2 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.

Sequencing and bioinformatic analysis
Sequence data from the 60 soil samples was analyzed using an in-house developed Python pipeline (https:// forge mia. inra. fr/ vasa/ illum iname tabar coding).
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 fulllength 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 PR 2 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.

Quantification of N-cycling bacterial and archaeal 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
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, lowabundance 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 Net-workAnalyser tool from Cytoscape was used to calculate network topological parameters (i.e. nodes, links, clustering coefficient and degree).

Land-use intensification impacts soil microbial communities
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).
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).

Identifying differentially abundant OTUs
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 Axis.  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).

Co-occurrence networks within microbial groups are influenced by land use
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.

Associations across microbial groups differ among land uses
To better understand how land use influences associations between bacteria, protist and fungi, we inferred cooccurrence 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).
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 landuse 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 (See figure on next page.) 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).  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 cooccurrence 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 rootassociated 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.