# Biogeography of Soil Bacterial Networks along a Gradient of Cropping Intensity

Battle Karimi, Samuel S. Dequiedt, Sébastien Terrat, Claudy C. Jolivet, Dominique D. Arrouays, Patrick Wincker, Corinne Cruaud, Antonio Bispo, Nicolas Chemidlin Prévost-Bouré, Lionel Ranjard

<figure><img src="image_1.png" /></figure>

To cite this version:

Battle Karimi, Samuel S. Dequiedt, Sébastien Terrat, Claudy C. Jolivet, Dominique D. Arrouays, et al.. Biogeography of Soil Bacterial Networks along a Gradient of Cropping Intensity. Scientific Reports, 2019, 9 (1), 10.1038/s41598-019-40422-y. hal-02093081

HAL Id: hal-02093081

https://institut-agro-dijon.hal.science/hal-02093081

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

OPEN

# Biogeography of Soil Bacterial Networks along a Gradient of Cropping Intensity

Battle Karimi<fnref n="1" />, Samuel Dequiedt<fnref n="1" />, Sébastien Terrat<fnref n="1" />, Claudy Jolivet², Dominique Arrouays², Patrick Wincker³, Corinne Cruaud³, Antonio Bispo²,⁴, Nicolas Chemidlin Prévost-Bouré<fnref n="1" /> & Lionel Ranjard<fnref n="1" />

Although land use drives soil bacterial diversity and community structure, little information about the bacterial interaction networks is available. Here, we investigated bacterial co-occurrence networks in soils under different types of land use (forests, grasslands, crops and vineyards) by sampling 1798 sites in the French Soil Quality Monitoring Network covering all of France. An increase in bacterial richness was observed from forests to vineyards, whereas network complexity respectively decreased from 16,430 links to 2,046. However, the ratio of positive to negative links within the bacterial networks ranged from 2.9 in forests to 5.5 in vineyards. Networks structure was centered on the most connected genera (called hub), which belonged to *Bacteroidetes* in forest and grassland soils, but to *Actinobacteria* in vineyard soils. Overall, our study revealed that soil perturbation due to intensive cropping reduces strongly the complexity of bacterial network although the richness is increased. Moreover, the hub genera within the bacterial community shifted from copiotrophic taxa in forest soils to more oligotrophic taxa in agricultural soils.

The response of soil bacterial communities to different environmental perturbations has been intensively investigated during the last decade. Studies have highlighted the significant role of soil characteristics<fnref n="1" />⁻³, plant communities, climate² and land use⁴ as drivers of the abundance, diversity and structure of soil microbial communities. However, the influence of these drivers on the number and intensity of interactions occurring between community members remains little known, partly because microbes, and most of their interactions, cannot be directly observed or measured⁵. Biological relationships or interactions are essential to ecosystem functions as previously demonstrated and discussed with regard to macroorganisms, as exemplified by plant-pollinators systems⁶,⁷, mycorhizal interactions⁸ and the regulation of biotic invasions⁹. Recently, a new metric has been developed to determine and assess the relationships occurring in microbial communities. This approach, based on evaluation of the co-abundance between taxa, highlights the positive and negative biological relationships, also known respectively as co-occurrence and co-exclusion¹⁰. Positive relationships in microbial communities can result from cooperative or co-dependence interactions, e.g. the facilitation of organic matter degradation, or the similarity of ecological niches between microorganisms. Negative relationships, in contrast, can be attributed to inhibitive or antagonistic interactions such as antibiosis or competition for resources¹¹. Thus, the co-occurrence approach contributes to a comprehensive evaluation of putative microbial relationships, independently of our technical capacity to observe and identify them¹⁰.

Analyzing the microbial co-occurrence network in soil could thus provide a promising way to improve our understanding of soil microbial community regulation, functioning and stability¹⁰,¹². However, the relationships between the microbial interaction network and biological function have not been clearly demonstrated in complex ecosystems. The only examples to date have been obtained on simplified microbial ecosystems such as biofilms, where the co-occurrence and interaction of different bacterial species enhance the biofilm biomass, tolerance of the community to physical (e.g. desiccation), chemical (e.g. sodium dodecyl sulfate), and biological stress (e.g. antimicrobial agents), microbial virulence and the degradation of organic compounds and pollutants¹³,¹⁴. Based on these demonstrations, it was assumed that, the diversity of microbial interactions in complex

Received: 20 November 2018

Accepted: 13 February 2019

Published online: 07 March 2019

<footnote n="1">¹Agroécologie, AgroSup Dijon, INRA, Université Bourgogne Franche-Comté, F- 21000, Dijon, France. ²INRA, US 1106, Unité INFOSOL, 45075, Orléans, France. ³CEA/Institut de Génomique/Genoscope, 91057, Evry cedex, France. ⁴ADEME, Service Agriculture et Forêt, 49000, Angers, France. Correspondence and requests for materials should be addressed to L.R. (email: lionel.ranjard@inra.fr)</footnote>

<footnote>Comment citer ce document :</footnote>

<footnote>Karimi, B., Dequiedt, S., Terrat, S., Jolivet, C., Arrouays, D., Wincker, P., Cruaud, C.,</footnote>

<footnote>SCIENTIFIC REPORTS | Preprint (2019) 9:3812 | https://doi.org/10.1038/s41598-019-40422-y
of Soil Bacterial Networks
along a Gradient of Cropping Intensity. Scientific Reports, 9 (1), DOI: 10.1038/s41598-019-40422-y</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Characterization of soil samples: (a) Map of the 1717 sampling sites used to compute the network replicates and classification for the four land uses encountered in France. (b–d) Distribution of the RMQS soils in the USDA soil texture triangle for each land use respectively. Color legend from red to blue indicates soil pH and the circle size represents the relative amount of organic carbon in the soil.</figcaption></figure>

ecosystems, such as soils, should affect biological functions at the ecosystem scale and that investigations of microbial network complexity might provide new insights into soil functioning and regulation<sup>5</sup>.

Recently, soil bacterial co-occurrence networks have been shown to be impacted by soil physico-chemical characteristics<sup>16,17</sup> and climate<sup>18</sup>, to vary with age of the plant cover<sup>15</sup> and to change with the soil land use<sup>19,20</sup>. Both studies by Morriën and Lupatini demonstrated that the least complex bacterial networks were found in cropped soils, i.e. those most disturbed by agricultural practices, as compared to pastured or forest soils. These studies provided the first proof of the sensitivity of microbial interactions networks to land management but were only carried out on a local spatial scale. Now large-scale demonstrations are needed to determine the genericity of these conclusions<sup>4</sup>. In this study, we investigated soil bacterial co-occurrence networks on the scale of France by clustering soils according to the main land uses encountered across this country (Fig. 1a). Data were obtained from the French Soil Quality Monitoring Network (RMQS), which represents the most intensive soil sampling
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Visualization of the most complex network among the 100 replicates for the 4 types of land use. The red edges represent the negative links and the green edges represent the positive links. The most complex network was the one with the most links, the highest connectance and the highest average degree.</figcaption></figure>

system on a wide spatial scale, due to its extensive area covered (5.5 10⁵ km²) and the high sampling resolution (about 2200 sites distributed along a systematic grid)⁴. The four land uses ranged from natural or semi-natural sites (forest and grassland) to agricultural sites (crops and vineyards), representing a gradient of cropping intensity, which integrated management of the plant cover and the degree of soil disturbance resulting from agricultural practices such as tillage, fertilization and pesticides inputs²¹ (see details in Material and Methods section). Previous studies based on RMQS data had revealed a lower microbial biomass but a higher bacterial taxonomic richness in vineyards and crop system than in forest and grassland soils³,²². Since these two microbial parameters were differently impacted by the land uses, investigating other microbial community parameters, such as the co-occurrence network, is a major concern to understand more comprehensively the response of microbial communities to the soil perturbations. Here, we evaluated and compared the structure of bacterial networks at the genus level by quantifying the complexity, the cohesion and the proportion of positive to negative relationships between land uses. We also identified the hub taxa for the bacterial network in each land use and related the ecological attributes of these taxa to the environmental and management context.

## Results

**Comparison of network structures between land uses.** Graphically, the networks were composed of nodes and connections (the links). The connections were the significant correlations between the nodes, which corresponded to the genera occurring in the soils under the respective land use. The connection could be positive (green color) or negative (red color). The more the node was connected, the nearer it was to the heart of the network. Visual comparison of the networks for each land use revealed a significant shift in structure ranging from a highly connected, tightly closed structure for forests to a sparse, open structure for vineyards (Fig. 2). The bacterial networks in forest soils formed a dense cluster. In grassland soils, the cluster seemed to split into two parts with several long chains extending from the clusters. In soils under crop systems, one of clusters split into several large satellites which remained connected. In vineyards soils, many of the links seemed to be lost and the satellites were smaller and less inter-connected. Statistical comparisons of the network structures between the land uses confirmed a highly significant decreasing gradient in network complexity with forest > grassland > crop

Comment to cite this document :

Karimi, B., Dequiedt, S., Terrat, S., Jolivet, C., Arrouays, D., Wincker, P., Cruaud, C.,
<table><caption>**Table 1.** Community metrics for the 4 land uses: 6 indices of bacterial diversity and 7 indices of bacterial network. Different letters indicate a significant statistical difference between land uses, based on a parametric variance analysis followed by a Tukey HSD post-hoc comparison (n = 100 for each land use). Letters indicate differences of values between land uses. They are provided by Tukey-HSD pairwise comparison.</caption><thead><tr><th rowspan="2">Community metrics</th><th rowspan="2"></th><th colspan="3">Forest</th><th colspan="3">Grassland</th><th colspan="3">Crop system</th><th colspan="3">Vineyards</th></tr><tr><th>Mean</th><th>sd</th><th>Diff.</th><th>Mean</th><th>sd</th><th>Diff.</th><th>Mean</th><th>sd</th><th>Diff.</th><th>Mean</th><th>sd</th><th>Differences</th></tr></thead><tbody><tr><td rowspan="6">Diversity</td><td>Richness in OTUs</td><td>1,083.3</td><td>224.0</td><td>c</td><td>1,238.1</td><td>157.1</td><td>b</td><td>1,305.9</td><td>167.9</td><td>a</td><td>1,342.8</td><td>152.3</td><td>a</td></tr><tr><td>Shannon's index for OTUs</td><td>5.27</td><td>0.39</td><td>c</td><td>5.41</td><td>0.32</td><td>b</td><td>5.50</td><td>0.31</td><td>a</td><td>5.57</td><td>0.23</td><td>a</td></tr><tr><td>Eveness for OTUs</td><td>0.18</td><td>0.04</td><td>ns</td><td>0.18</td><td>0.04</td><td>ns</td><td>0.19</td><td>0.04</td><td>ns</td><td>0.19</td><td>0.03</td><td>ns</td></tr><tr><td>Richness in genera</td><td>284.6</td><td>54.1</td><td>c</td><td>333.4</td><td>32.0</td><td>b</td><td>360.2</td><td>28.4</td><td>a</td><td>359.5</td><td>24.6</td><td>a</td></tr><tr><td>Shannon's index for genera</td><td>4.04</td><td>0.31</td><td>c</td><td>4.19</td><td>0.21</td><td>b</td><td>4.25</td><td>0.2</td><td>a</td><td>4.26</td><td>0.12</td><td>ab</td></tr><tr><td>Eveness for genera</td><td>0.21</td><td>0.03</td><td>ns</td><td>0.2</td><td>0.03</td><td>ns</td><td>0.2</td><td>0.03</td><td>ns</td><td>0.2</td><td>0.02</td><td>ns</td></tr><tr><td rowspan="6">Co-occurrence network</td><td>Number of links</td><td>16,430.8</td><td>3,109.3</td><td>a</td><td>5,846.4</td><td>2,208.0</td><td>b</td><td>4,122.6</td><td>1,138.8</td><td>c</td><td>2,046.7</td><td>281.3</td><td>d</td></tr><tr><td>Number of positive links</td><td>12,042.4</td><td>2174.9</td><td>a</td><td>4,539.5</td><td>1,917.4</td><td>b</td><td>3,101.2</td><td>697.3</td><td>c</td><td>1,794.0</td><td>172.2</td><td>d</td></tr><tr><td>Number of negative links</td><td>4,388.4</td><td>1,123.6</td><td>a</td><td>1,306.8</td><td>486.5</td><td>b</td><td>1,021.4</td><td>475.8</td><td>c</td><td>342.7</td><td>115.3</td><td>d</td></tr><tr><td>Positive:Negative links</td><td>2.9</td><td>0.7</td><td>c</td><td>3.7</td><td>1.8</td><td>b</td><td>3.5</td><td>1.2</td><td>b</td><td>5.5</td><td>1.7</td><td>a</td></tr><tr><td>Connectance</td><td>0.008</td><td>0.002</td><td>a</td><td>0.003</td><td>0.001</td><td>b</td><td>0.002</td><td>0.0005</td><td>c</td><td>0.001</td><td>0.0001</td><td>d</td></tr><tr><td>Average path length</td><td>0.139</td><td>0.007</td><td>a</td><td>0.129</td><td>0.019</td><td>b</td><td>0.116</td><td>0.016</td><td>c</td><td>0.080</td><td>0.009</td><td>d</td></tr><tr><td></td><td>Average degree</td><td>16.2</td><td>3.1</td><td>a</td><td>5.8</td><td>2.2</td><td>b</td><td>4.1</td><td>1.1</td><td>c</td><td>1.0</td><td>0.3</td><td>d</td></tr></tbody></table>

system > vineyards soils (Table 1). The number of links and the connectance were progressively reduced by 87% from forests to vineyards. Similarly, the average degree of the networks and the average path length decreased from forest to vineyards by 94% and 42%, respectively. In addition, the average number of links by genus was 16.2 in forest soils, 5.8 in grassland soils, 4.1 in crop system soils and only 2.0 in vineyards soils. A significant decrease in the number of positive and negative links was also observed from forest to vineyards (a loss of 85% and 92%, respectively). However, the greater decrease of negative links in vineyards, as compared to forests, led to a higher positive to negative ratio for vineyards. The positive to negative ratio increased from 2.9 in forest soils to 5.5 in vineyards soils.

To investigate the assumption that network structure and complexity was directly related to community diversity, the richness, the Shannon index and the Shannon evenness in genera and in OTUs were also compared between land uses (Table 1). For both genus and OTU level, a significant increasing gradient (of about 25%, p-value < 0.001) was recorded along the cropping intensity gradient, the lowest richness occurring in forest soils (285 genera and 1083 OTUs on average) and the greatest richness in crop systems and especially in vineyards (360 genera and 1324 OTUs on average). The Shannon indices followed the same gradient as the richness, when the Shannon evenness' did not differ significantly between land use.

**Identification of hub genera.** The hub genera are those genera which exhibit the most links within a network, in other words the highest degree. Here, the hub genera were identified by statistically comparing the degrees obtained from the 100 replicates of network per land use, using a non-parametric pairwise test. The main hubs, their average degrees based on the 100 replicates and their rank within each land use are provided in Table 2. The quantitative analysis of these hubs highlighted a decrease of 80% in the degree (the number of links of the genera) from forests to vineyards. In forest, the primary hub was *Pirellula* with 337.7 links, whereas in grassland and crop systems, the primary hub was *Acidicaldus* with 183.9 links and 161.7 links, respectively, and *Bradyrhizobium* in vineyards with 73.9 links. *Pirellula* was also found among the 20 main hubs in grassland and vineyards soils whereas *Acidicaldus* was the third hub in vineyards and *Bradyrhizobium* was also among the 20 main hubs under crop systems. Interestingly, some of the major hubs were highly specific to a single type of land use, e.g. *Flavobacterium* found exclusively in forest soils and *Bradyrhizobium* found exclusively in cropped soils.

## Discussion

**Global network structure changes between land uses.** Visual analysis of the bacterial networks showed changes in the bacterial networks along the gradient of cropping intensity (Fig. 2). The soil bacterial networks seemed to progressively split and to shift from a single block with highly connected structure in forest soils to an open structure, full of small slightly inter-connected satellites in vineyards soils. This network structure in vineyards soils suggested an increasing isolation of several bacterial genera, which interacted little with the rest of the community. This particular evolution in network structure has not been described in previous studies of the impact of land uses or soil parameters. This might be explained by the lack of network visualization¹⁹,²⁰ or the partial representation of networks²³,²⁴.

Regarding network metrics, the number of links and the connectance quantify the direct relationships and the network complexity. The average path length includes the indirect relationships, sometimes producing long chains with several intermediaries and translates the network cohesion. Finally, the positive to negative ratio

Comment cite this document:

Karimi, B., Dequiedt, S., Terrat, S., Jolivet, C., Arrouays, D., Wincker, P., Cruaud, C.,
<table><caption>**Table 2.** Hub genera for the 4 land uses: the table summarizes the specific hubs for each land use and indicates the average degree of the genera for each of the land uses, the statistical difference with the other genera in the respective land use (given by a Kruskal-Wallis pairwise test) and, in parenthesis, the ranking of each genus among the most connected genera for the land use. Letters indicate differences of degree between nodes within each land use. They are provided by Pairwise Kruskal-Wallis Rank Test corrected by Bonferroni method.</caption><thead><tr><th rowspan="2">Top Hub Genera</th><th rowspan="2">Phylum-level classification of genus</th><th colspan="2">Forest</th><th colspan="2">Grassland</th><th colspan="2">Crop system</th><th colspan="2">Vineyards</th></tr><tr><th>Mean</th><th>Differences</th><th>Mean</th><th>Differences</th><th>Mean</th><th>Differences</th><th>Mean</th><th>Differences</th></tr></thead><tbody><tr><td>Flavobacterium</td><td>Bacteroidetes</td><td>367.4</td><td><b>ab (2)</b></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Blastochloris</td><td>Alphaproteobacteria</td><td>361.0</td><td>abcd (4)</td><td>161.1</td><td>abcd (6)</td><td></td><td></td><td></td><td></td></tr><tr><td>Terrimonas</td><td>Bacteroidetes</td><td>366.8</td><td><b>abc (3)</b></td><td>170.6</td><td><b>ab (2)</b></td><td>133.2</td><td>bc (4)</td><td></td><td></td></tr><tr><td>Pirellula</td><td>Planctomycetes</td><td>377.7</td><td><b>a (1)</b></td><td>136.6</td><td>efg (14)</td><td></td><td></td><td>43.8</td><td>gh (11)</td></tr><tr><td>Geothrix</td><td>Acidobacteria</td><td>358.5</td><td>abcd (5)</td><td>155.0</td><td>bcde (7)</td><td>128.2</td><td>bcd (6)</td><td>42.6</td><td>gh (12)</td></tr><tr><td>Stella</td><td>Alphaproteobacteria</td><td>342.4</td><td>bcd (16)</td><td>170.2</td><td><b>ab (3)</b></td><td>151.0</td><td><b>ab (2)</b></td><td>63.9</td><td>abc (4)</td></tr><tr><td>Acidicaldus</td><td>Alphaproteobacteria</td><td></td><td></td><td>183.9</td><td><b>a (1)</b></td><td>161.7</td><td><b>a (1)</b></td><td>73.2</td><td><b>ab (3)</b></td></tr><tr><td>Isosphaera</td><td>Planctomycetes</td><td></td><td></td><td>162.7</td><td>abc (5)</td><td>149.9</td><td><b>ab (3)</b></td><td>53.7</td><td>cdef (7)</td></tr><tr><td>Frigoribacterium</td><td>Actinobacteria</td><td></td><td></td><td>144.7</td><td>cdefg (11)</td><td>128.9</td><td>bcd (5)</td><td>58.5</td><td>cde (6)</td></tr><tr><td>Catenulispora</td><td>Actinobacteria</td><td></td><td></td><td>126.2</td><td>g (18)</td><td>106.7</td><td>defg (16)</td><td>75.7</td><td><b>a (2)</b></td></tr><tr><td>Bradyrhizobium</td><td>Alphaproteobacteria</td><td></td><td></td><td></td><td></td><td>112.8</td><td>cdefg (12)</td><td>73.9</td><td><b>a (1)</b></td></tr></tbody></table>

indicates the balance between facilitative and inhibitive relationships within the network¹¹. In our study, all these metrics revealed a significant shift in bacterial network structure along the gradient of cropping intensity. The progressive decrease in the number of links, the connectance and the average path length observed from forest to vineyards soils, demonstrated a loss of network complexity and cohesion as agricultural practices were intensified, and soils disturbed. Surprisingly, this result contrasted with the significant increase in bacterial richness observed across RMQS soils, where the sequence was: forests < grasslands < crop systems = vineyards³ (Table 1). Theoretically, the presence of more taxa within a community would imply more potential interactions. Despite the increased number of potential links, the most taxa-rich communities may not necessarily exhibit the most complex interaction networks¹⁹. This has already been observed in the case of soil disturbance. In Karimi *et al.*, transplanting of samples from rural conditions to urban and industrial sites led to the loss of 40% of the links after 8 months, without any decrease of the microbial diversity²⁵. Similarly, in another study¹⁶, the contamination of a forest soil used as a tailings dump decreased the number of links in the bacterial co-occurrence network by 25% compared to the pristine soil while the bacterial richness was increased by 20%. Interestingly, our results confirmed that the pool of taxa and the biotic relationships within the community responded independently to environmental perturbation¹⁵.

The observed decrease in network complexity and cohesion supports the hypothesis that cropping may enhance the isolation of bacterial taxa, as previously suggested by visual observation. Isolation can be enhanced in three ways. Firstly, it can result from the loss of microbial biomass and therefore of cells. This hypothesis has been supported by the decreasing microbial molecular biomass in crop systems and vineyards compared to forest and grassland²². If the total number of cells in soil is reduced, the probability of each cell encountering another and interacting, whether directly or indirectly, may be also reduced. Secondly, isolation can be caused by the stimulated metabolic independence of microorganisms in agricultural soils. Soil disturbance, produced by agricultural practices (such as tillage), is known to increase soil bacterial diversity and to preferentially stimulate some opportunistic and/or pathogenic populations²⁶. These populations are self-sufficient and do not need to interact with others to be metabolically efficient and increase their fitness, which can therefore lead to even greater isolation of other bacterial populations. Finally, the isolation can be related to the weak spatial connectivity occurring between soil ecological niches in disturbed soils²⁷. At the microscale, the structure of tilled soils is more homogeneous, and the pores are less connected than in soils under minimum or without tillage²⁸. This can further induce the physical isolation of bacterial taxa in disturbed soils²⁹. As microbial co-occurrence networks translate both the biotic interactions and niche-sharing³⁰, the lower level of microbial networks cohesion (found in crop systems and vineyards) might be associated with the less spatially-structured distribution of ecological niches in such soils due to their frequent disturbance by tillage or fertilization. Altogether, these results suggest that forest and grassland soils represent a mosaic of connected ecological niches that are fully complete and shared by non-opportunistic bacterial genera when the soil structure and trophic resources are non-limiting³¹,³² (Fig. 2a,b). On the contrary, agricultural soils, especially vineyards, would consist of a mosaic of habitats with weakly-connected, partially-filled niches where the bacterial genera tend to be spatially and metabolically isolated by an altered soil structure and limited trophic resources³³,³⁴ (Fig. 2c,d). Unlike the other metrics, the ratio between positive and negative links increased from forest soils to vineyards soils. Although both positive and negative relationships decreased, the negative ones decreased more rapidly. This suggests that inhibition, rather than cooperation between bacteria taxa, was more affected by the cropping intensity. The weak availability of the soil organic matter in agricultural soils (Fig. 1) represented harsh conditions which could lead to selection of the most adapted taxa. Taxa with ecological attributes such as oligotrophy, ability to degrade recalcitrant matter or anaerobic metabolism, would be able to share this soil habitat³⁰,³⁵ and thus to maintain some of the positive relationships at the scale of bacterial network³⁰. This supports the concept of resource-driven co-occurrence patterns³⁶.
**Shift in hubs networks along the gradient of cropping intensity.** The hubs represent the most connected taxa within a network<sup>36,37</sup>. Most hub genera identified in our study (Table 2) were consistent with a recent review of microbial keystone taxa in different ecosystems based on 37 studies published between 2011 and 2017<sup>36</sup>. *Rhizobiales* (like *Blastochloris* genus in our study) and *Terrimonas* genus have previously been recorded as hub/keystone taxa in forest and woodlands. In grasslands, bacteria belonging the phyla Bacteroidetes, Planctomycetes, Acidobacteria and Actinobacteria were also identified as hubs. And in agricultural soils, *Rhodospirillales* (*Stella* and *Acidicaldus* here), *Rhizobiales* (*Bradyrhizobium* here) and *Actinomycetales* (*Frigoribacterium* and *Catenulispora* here) were keystone taxa. Due to their high level of connectivity with other community members, these taxa could have a key role in the soil bacterial community of these land uses<sup>36</sup>.

Given the known ecological attributes of the bacterial hubs identified here, the progressive shift of these hubs along the gradient of cropping intensity might be partly due to the soil trophic level. In cropped soils (vineyards and crop systems), demonstrated to be the poorest of the RMQS soils in terms of organic carbon and nutrient contents<sup>32</sup>, the major contribution of *Bradyrhizobium*, *Acidicaldus* and *Isosphaera* to the network may be explained by their oligotrophic nature<sup>38-40</sup>. Rhizobia, such as *Bradyrhizobium*, are also known to develop highly specific symbioses with their host plant<sup>41</sup>, which might enhance their co-occurrence with other bacterial populations associated with the rhizosphere. Conversely, Bacteroidetes genera (*Flavobacterium* and *Terrimonas*) were the hubs of bacterial networks in forest and grassland soils, the richest of the RMQS soils in terms of organic matter content<sup>32</sup>. The bacterial phylum Bacteroidetes is known to be copiotrophic and more abundant in forest soils and therefore to interact strongly with other taxa<sup>42,43</sup>. However, few ecological data concerning the biotic interactions and ecological attributes of most of these genera are currently available, which therefore limits our ability to draw conclusions.

Interestingly, our results also revealed that *Flavobacterium* and *Terrimonas* were, at the same time, network hubs and two of the most abundant genera in forest soils<sup>43</sup>. Even if the most connected taxa might also be expected to be the most abundant, this hypothesis was not corroborated for the other land uses. The respective hub genera in grassland, crop system and vineyards soils had previously been shown to be only slightly abundant in these soils<sup>43</sup>. The identification of both abundant and hub taxa provides complementary information about the bacterial community ecology. It also reinforces the need to investigate different parameters describing soil microbial communities<sup>11</sup> so as to better predict their responses to the changing environment and propose the most appropriate soil management practices.

**Robustness and limitation of the sampling and analytical strategies.** This work is based on the most intensive, without *a priori*, soil sampling survey (about two thousand soil samples) focusing on a nation-wide scale. The sampling design covered the major environmental variability across a 550,000 km² area<sup>4</sup>, thus precluding the bias associated with *a priori* samplings. This enable us to construct robust networks with a reasonable number of repetitions to verify our assumption concerning the effect of the cropping intensity. Compared to other studies based on 3 to 20 samples, our sampling design provides a cornerstone for robust analysis and conclusions about soil bacterial biogeography. The observed variability in network structure and metrics, (computed from 100 replicates) within each land use, was probably due to the variability in pedoclimatic conditions, i.e. the interacting soil characteristics, and/or agricultural practices recorded for each land use<sup>4</sup>. Nevertheless, according to the multi-dimensional plots (Fig. 1 and Supplementary Fig. 1), the environmental heterogeneity was similar for all four land uses. Moreover, as the local environment determined the microbial biomass<sup>1</sup>, the bacterial diversity<sup>3</sup> and the distribution of taxa<sup>43</sup>, the network structure could also be influenced by the local spatial effect. In our sub sampling, the 30 sites covered an area ranging from 200,000 km² to 550,000 km² across a territory of 551,500 km². However, the area covered had no impact on the network structure, whatever the land use (Supplementary Fig. 2 and Supplementary Table 1). The networks based on sites located closer together were not different from the others, which confirmed that the local environmental filters did not affect the bacterial network. We therefore assumed that the results obtained for the bacterial co-occurrence networks were only due to the different land uses and to the cropping intensities that they represent.

Our molecular analytical strategy is known to be highly robust<sup>44</sup>, although numerous biases inherent in amplification library preparation such as DNA extraction<sup>45</sup>, amplification<sup>46</sup>, sequencing and inference of the patterns of organism abundance from library data pertaining to relative abundance, are also well-known<sup>47</sup>. Analyses were conducted in a consistent manner to remove errors due to sequencing and chimeras, and the datasets were refied to the same sampling depth (i.e. 10,000 reads per sample), so that relative changes in microbial taxonomic composition levels could be compared across samples, even if the biodiversity sampling was not exhaustive<sup>47</sup>. Nevertheless, the comparison of our results with those from other studies is mainly limited by our choice of 16S rRNA primers, which were designed to specifically target both bacteria and archaea diversity. The results obtained by using these primers were able to reveal taxonomic groups (for example *Holophaga*)<sup>43</sup>, rarely detected in previous soil studies. In addition to using these different primers, we also chose a finer but more time-consuming method of taxonomic assignment than the approach currently used in QIIME: Instead of assigning the seed sequence of each OTU, all reads in the dataset were individually assigned. This led to changes in relative abundance of the taxonomic groups within the community, probably in the actual structure of the co-occurrence network and finally the revelation of hub *genera*, probably underestimated in other studies.

**Conclusions**

Altogether our study demonstrated that soil bacterial co-occurrence networks are different between land use types and are strongly shaped by the cropping intensity. We hypothesized that changes in the bacterial network would occur mainly in response to shifts in the heterogeneity and connectivity of the mosaic of microbial habitats as well as to the availability of C-substrates. Beyond the classical information obtained from bacterial richness or whole taxonomic composition, co-occurrence network analysis provides complementary insights in biotic

Comment cite this document :

Karimi, B., Dequied, S., Terrat, S., Jolivet, C., Arrouays, D., Wincker, P., Cruaud, C.,
interactions and niche connectivity, which could have repercussions on community stability⁴⁸ and soil functioning⁴⁹. Further investigations must now be focused on the dynamics of these networks to evaluate their response to environmental perturbations and also on the fungal community to obtain an overview of the soil microbial network.

## Material and Methods

**Sampling Design.** Soil samples were obtained from the French Soil Quality Monitoring Network ("Réseau de Mesures de la Qualité des Sols", RMQS) which is a soil monitoring network based on a 16 km regular grid across the 550 000 km² French territory¹,³. The soil profile, site environment, climatic factors, location, vegetation and land management for the 2173 RMQS monitoring sites have already been described. All samples were collected between 2002 and 2009. All sites have been geo-positioned with a precision of <0.5 m and the soil profile, site environment, climatic factors, vegetation and land use described¹. In the middle of each 16 × 16 km square, 25 individual core samples were taken from the topsoil (0–30 cm) using an unaligned sampling design within an area of 20 × 20 m. These core samples were bulked to obtain a composite sample for each RMQS site. The soil samples were gently air-dried, sieved to 2 mm and then stored at -40 °C before analysis. Physico-chemical parameters were measured for each composite soil, e.g. particle-size distribution, pH water, organic C, N, C/N ratio, soluble P contents, calcareous, cation exchange capacity (CEC) or exchangeable cations (Ca, Mg). Physical and chemical analyses are available for 2,131 soils and were performed by the Soil Analysis Laboratory of INRA (Arras, France, http://www.lille.inra.fr/las). Available climatic data for the RMQS were monthly rainfall, evapotranspiration and temperature at each node of a 12 × 12 km² grid, averaged for the 1992–2004 period. These climatic data were obtained by interpolating observational data using the SAFRAN model. The RMQS site-specific data were linked to the climatic data by finding the closest node within the 12 × 12 km² climatic grid for each RMQS site. Land cover was recorded according to the coarse level of the CORINE Land Cover classification (IFEN, http://www.statistiques.developpement-durable.gouv.fr/donnees-ligne/li/2539/0/base-donnees-geographique-corine-land-cover-clc.html), which consists of a rough descriptive classification into five classes: forest, croplands, grasslands, perennial crops (corresponding to vineyards and orchards) and others. All these data were available in the DONESOL database¹,³².

**Characterization of the gradient of cropping intensity.** The gradient of cropping intensity was composed of four land uses that were characterized by the management and turnover of the plant cover, the duration of the management plan, and the frequency and intensity of interventions which disturbed the soil: (i) The "forest" land use grouped together broad-leaved forests, coniferous forests, poplar grove and mixed forests, which could be natural or planted. The management plans were often established over 10 to 150 years with few perturbations until the wood cutting. These forest soils are not tillied, and few engines pass over them. There is no fertilization and no pesticide input. Turnover of the plant cover is on the scale of decades. (ii) The "grassland" land use grouped together 6 to 10-years old seeded grassland, >10-years seeded or natural grassland, unproductive permanent grasslands and meadows. Some of these grasslands may be seasonally pastured. Tillage is rare. In France, less than 30% of grasslands receive organic manure and 40% on average receive low doses of mineral fertilizer (50 kgN/ha on average). Only 6% of grasslands receive pesticides (data from the Ministry of Agriculture, 2011). Turnover of the plant cover exceeds 6 years. (iii) The "crop system" land use grouped together annual crops in monoculture, annual crops associated with permanent crops, complex cultivation patterns, market gardening and ornamental horticulture, or fallow land. These systems are strongly managed to ensure satisfactory yields. During the year, the different interventions concern soil tillage, seedling, fertilization and pesticides applications. Even if the frequency and intensity of interventions can vary considerably across the crop systems, they take place every year. Turnover of the plant cover is annual. (iv) The "vineyards" land use grouped together vineyards, orchards and others perennial tree crops. Like the crop systems, the annual interventions concern soil tillage, fertilization and pesticides applications. However, the frequency and intensity of pesticides input are greater than for annual crop systems due to the perennial implantation of the vines and trees.

Overall, the forest soils are the least disturbed and the vineyards soils are the most disturbed by the agricultural practice and interventions over time. These agricultural practices have also affected the soil nutrient status by determining the amount and the quality of organic matter. Thus, these land uses ranged from natural or semi-natural sites (forest and grassland) to agricultural sites (crop and vineyards). All soil samples were therefore classified according to the land use of the site: 492 forest soil samples, 464 grassland soil samples, 740 cropped soil samples and 36 vineyard or orchard soil samples.

**Molecular characterization of bacterial communities.** The analysis of these 2173 samples required the rigorous standardization of the range of different molecular tools involved in soil DNA extraction and sequencing technology. The following protocols were applied.

**Soil DNA extraction and purification.** Microbial DNA was extracted and purified from 1 g of each of the 2173 composite soils (composed of a bulk of 25 individual core soils) sampled at each RMQS site, using the previously-described GnS-GII procedure⁴⁵. Crude DNA extracts were quantified by agarose gel electrophoresis stained with ethidium bromide, using calf thymus DNA as standard curve¹. Crude DNA was then purified using a MinElute gel extraction kit (Qiagen, France) and quantified using a QuantiFluor staining kit (Promega, USA), prior to further investigations.

**PCR amplification and pyrosequencing of 16S rRNA gene sequences.** A 16S rRNA gene fragment targeting the V3-V4 regions to characterize bacterial diversity was amplified using the primers F479 (5'-CAGCMGYGCNGTAANAC-3') and R888 (5'-CCGYCAATTCMTTTRAGT-3') with the method described previously⁴⁵. From the 2173 DNA soil

Comment citere ce document :

Karimi, B., Dequiedt, S., Terrat, S., Jolivet, C., Arrouays, D., Wincker, P., Cruaud, C.,
samples, 2132 soil samples were successfully amplified. The PCR products were then purified using a MinElute PCR purification kit (Qiagen) and quantified using the QuantiFluor staining kit (Promega, USA). A second PCR of 7 cycles was then run twice for each sample under similar PCR conditions, with purified PCR products as matrix (7.5 ng of DNA were used for a 25 µl mix of PCR) and dedicated fusion primers ('F479/AdaptorB', 'R888/MID/AdaptorA') integrating the required adaptors, keys and multiplex identifiers at the 5' extremities. All duplicated PCR products were then pooled, purified using a MinElute PCR purification kit (Qiagen), and quantified using the QuantiFluor staining kit (Promega, USA). For all libraries, equal amounts from 30 samples were pooled, and then cleaned to remove excess nucleotides, salts and enzymes using the Agencourt AMPure XP system (Beckman Coulter Genomics). 100 µl of TE buffer (Roche) was used for the elution. Pyrosequencing was then carried out on a GS FLX Titanium (Roche 454 Sequencing System) by Genoscope (Evry, France).

**Bioinformatics sequence analysis.** Bioinformatic analyses were done using the GnS-PIPE pipeline developed by the GenoSol platform (INRA, Dijon, France) (availability: https://doi.org/10.5281/zenodo.1123425) and previously detailed<sup>3,43</sup>. After sequencing, 49,794,516 raw reads were obtained for the 2132 soil samples. After the different preprocessing and filtering steps (detailed previously<sup>3</sup>), 32 634 692 high-quality reads (range of sequencing depth: 48 reads to 49 926 reads by sample) were kept. The number of high-quality reads for each sample was then "rarefied" (i.e. 10 000 high-quality reads for each sample) by random selection to allow efficient comparison of the datasets and avoid biased community comparisons and rarefaction curves. Thus, 1798 soil samples were kept for subsequent analyses, encompassing a total of 17 980 000 reads. A final post-processing step was then applied to this global dataset, as already described regarding for microbial richness across France<sup>3</sup>, to filter the potentially artefactual reads. Briefly, the 17 980 000 reads from all samples were firstly aligned and clustered at 95% of similarity into OTUs. Thereafter, all OTUs that occurred only once in the overall dataset and contained only a single read, were removed. This post-processing step reduced the number of total OTUs from 184 812 to 83 917 (more than 50% lost), but the number of reads only from 17 980 000 to 17 866 981 (less than 1% lost). The number of deleted reads by sample was 62 ± 60 on average (minimum: 10, maximum: 1093). Finally, all kept reads were then compared with the dedicated reference database derived from SILVA to independently determine the composition of each soil community at the *genus* level (procedure, database and programs available online: https://doi.org/10.5281/zenodo.1065438 and https://doi.org/10.5281/zenodo.1064170). The final dataset was composed of 1798 samples with 10 000 sequences per sample distributed among 1355 genera<sup>43</sup>. Unknown sequences represented an average of 11% by sample at the *genus* level. All raw data sets are publicly available in the EBI database system (in the Short Read Archive) under project accession no PRJEB21351.

**Statistical and co-occurrence network analyses.** Comparing the bacterial networks between the four land uses required: i) standardization of the number of soils used to compute the network by land use to avoid a sampling size effect and ii) network replicates which integrated the natural heterogeneity of the soils within each land use. As vineyards provided only 36 soil samples on the territory scale, we standardized the sampling size at 30 sites per network for all land uses, drawn from the respective pool of soils. Thus, the minimum number of combinations (1 947 792) ensured that each network was computed from a unique combination of sites. The number of replicates used to evaluate the range of variation within and between land uses was set at 100. The 30 samples were drawn independently for each of the 100 network replicates. Numerical sampling of the sites was monitored. Among the 1732 sites classified as forest, grassland, crop system or vineyards, 1717 were sampled at least once. This replication procedure ensures that more than 99% of the available sites were sampled, standardized analysis of all four land uses, while being representative of French soils in general.

Then, for each replicate, network computation was based on a contingency matrix which provided the relative abundances of 1355 bacterial genera for the 30 randomly-selected soil samples. The Spearman correlation coefficient for each pair of genera was used as a similarity index to estimate taxa co-occurrence<sup>50</sup>. A correlation was considered as robust and non-random if the p-value was below 0.05 after correction using the False-Discovery Rate method<sup>51</sup>. By choosing this threshold, only correlations with a coefficient above 0.50 were kept. The positive and negative correlations were respectively interpreted as co-occurrence (facilitation) and co-exclusion (antagonistic) relationships. To describe the topology of the resulting networks, a set of metrics were calculated, namely the number of links, the number of positive links, the number of negative links, the ratio between the positive and negative links, the connectance defined as the proportion of potential links which are significantly observed, the average path length and the average degree (for more details, see<sup>11</sup>). Due to the homoscedasticity between land-use and the Gaussian distribution of the model residuals, these metrics were then compared between the four land uses by applying a parametric analysis of variance and a Tukey HSD post-hoc comparison. The degree (defined as the number of edges of each node) of the 20 most connected genera was recorded for all 100 replicates of all four land uses and then statistically compared within each land use by applying a pairwise Kruskal-Wallis rank test corrected by the Bonferroni method. The genera with the highest average degree were then considered as hub genera.

In addition, the bacterial diversity was evaluated for each land use at the taxonomic level of OTUs and genera. Three indices were computed: a- the richness $N_0$ which corresponds to the number of taxa, b- the Shannon entropy index $H = -\sum_i (P_i * \log P_i)$, where $P_i$ is the probability of finding the taxon i and c- the Shannon evenness $E = \exp(H)/N_0$. These indices were compared between the four land uses by applying a parametric analysis of variance and a Tukey HSD post-hoc comparison.

**Data Availability**

All raw data sets are publicly available in the EBI database system (in the Short Read Archive) under project accession no PRJEB21351.
References

1. Dequiedt, S. et al. Biogeographical patterns of soil molecular microbial biomass as influenced by soil characteristics and management. *Glob. Ecol. Biogeogr.* **20**, 641–652 (2011).

2. O'Brien, S. L. et al. Spatial scale drives patterns in soil bacterial diversity. *Environ. Microbiol.* **18**, 2039–2051 (2016).

3. Terrat, S. et al. Mapping and predictive variations of soil bacterial richness across France. *PLoS One* **12**, e0186766 (2017).

4. Ranjard, L. et al. Turnover of soil bacterial diversity driven by wide-scale environmental heterogeneity. *Nat. Commun.* **4**, 1434 (2013).

5. Deng, Y. et al. Network succession reveals the importance of competition in response to emulsified vegetable oil amendment for uranium bioremediation. *Environ. Microbiol.* **18**, 205–218 (2016).

6. Schleuning, M., Fründ, J. & García, D. Predicting ecosystem functions from biodiversity and mutualistic networks: An extension of trait-based concepts to plant-animal interactions. *Ecography (Cop).* **38**, 380–392 (2015).

7. Burkle, L. A., Marlin, J. C. & Knight, T. M. Plant-Pollinator Interactions over 120 Years: Loss of Species, Co-Occurrence and Function. *Science* **(80-).** **339**, 1611–1616 (2013).

8. Van Der Heijden, M. G. A. et al. A widespread plant-fungal-bacterial symbiosis promotes plant biodiversity, plant nutrition and seedling recruitment. *ISME J.* **10**, 389–399 (2016).

9. Massol, F. & Petit, S. Interaction Networks in Agricultural Landscape Mosaics. *Adv. Ecol. Res.* **49**, 291–338 (2013).

10. Faust, K. & Raes, J. Microbial interactions: From networks to models. *Nat. Rev. Microbiol.* **10**, 538–550 (2012).

11. Karimi, B. et al. Microbial diversity and ecological networks as indicators of environmental quality. *Environ. Chem. Lett.* **15**, 265–281 (2017).

12. Fuhrman, J. A. Microbial community structure and its functional implications. *Nature* **459**, 193–199 (2009).

13. Burmølle, M., Ren, D., Bjarnsholt, T. & Sørensen, S. J. Interactions in multispecies biofilms: Do they actually matter? *Trends Microbiol.* **22**, 84–91 (2014).

14. Madsen, J. S., Sørensen, S. J. & Burmølle, M. Bacterial social interactions and the emergence of community-intrinsic properties. *Curr. Opin. Microbiol.* **42**, 104–109 (2018).

15. Shi, S. et al. The interconnected rhizosphere: High network complexity dominates rhizosphere assemblages. *Ecol. Lett.* **19**, 926–936 (2016).

16. Zappelini, C. et al. Diversity and complexity of microbial communities from a chlor-alkali tailings dump. *Soil Biol. Biochem.* **90**, 101–110 (2015).

17. Mandakovic, D. et al. Structure and co-occurrence patterns in microbial communities under acute environmental stress reveal ecological factors fostering resilience. *Sci. Rep.* **8**, 1–12 (2018).

18. Ma, B. et al. Geographic patterns of co-occurrence network topological features for soil microbiota at continental scale in eastern China. *ISME J.* **10**, 1891–1901 (2016).

19. Lupatini, M. et al. Network topology reveals high connectance levels and few key microbial genera within soils. *Front. Environ. Sci.* **2**, 1–11 (2014).

20. Morriën, E. et al. Soil networks become more connected and take up more carbon as nature restoration progresses. *Nat. Commun.* **8** (2017).

21. Constancias, F. et al. Contrasting spatial patterns and ecological attributes of soil bacterial and archaeal taxa across a landscape. *Microbiologyopen* **4**, 518–531 (2015).

22. Horrigue, W. et al. Predictive model of soil molecular microbial biomass. *Ecol. Indic.* **64**, 203–211 (2016).

23. Lentendu, G. et al. Effects of long-term differential fertilization on eukaryotic microbial communities in an arable soil: A multiple barcoding approach. *Mol. Ecol.* **23**, 3341–3355 (2014).

24. Steele, J. A. et al. Marine bacterial, archaeal and protistan association networks reveal ecological linkages. *ISME J.* **5**, 1414–1425 (2011).

25. Karimi, B., Meyer, C., Gilbert, D. & Bernard, N. Air pollution below WHO levels decreases by 40% the links of terrestrial microbial networks. *Environ. Chem. Lett.* **14**, 467–475 (2016).

26. Lienhard, P. et al. Pyrosequencing evidences the impact of cropping on soil bacterial and fungal diversity in Laos tropical grassland. *Agron. Sustain. Dev.* **34**, 525–533 (2014).

27. Cordero, O. X. & Datta, M. S. Microbial interactions and community assembly at microscales. *Curr. Opin. Microbiol.* **31**, 227–234.

28. Pagliai, M., Vignozzi, N. & Pellegrini, S. Soil structure and the effect of management practices. *Soil Tillage Res.* **79**, 131–143 (2004).

29. Tecon, R. & Or, D. Biophysical processes supporting the diversity of microbial life in soil. *FEMS Microbiol. Rev.* **41**, 599–623 (2017).

30. Freilich, M. A., Wieters, E., Broitman, B. R., Marquet, P. A. & Navarrete, S. A. Species co-occurrence networks: Can they reveal trophic and non-trophic interactions in ecological communities? *Ecology* **99**, 690–699 (2018).

31. de Boer, W., Folman, L. B., Summerbell, R. C. & Boddy, L. Living in a fungal world: impact of fungi on soil bacterial niche development. *FEMS Microbiol. Rev.* **29**, 795–811 (2005).

32. Meersmans, J. et al. A high resolution map of French soil organic carbon. *Agron. Sustain. Dev.* **32**, 841–851 (2012).

33. Stoyan, H. et al. Spatial heterogeneity of soil respiration and related properties at the plant scale. *Plant Soil* **222**, 203–214 (2000).

34. Vos, M., Wolf, A. B., Jennings, S. J. & Kowalchuk, G. A. Micro-scale determinants of bacterial diversity in soil. *FEMS Microbiol. Rev.* **37**, 936–954 (2013).

35. Ma, Z. S. The P/N (Positive-to-Negative Links) Ratio in Complex Networks—A Promising In Silico Biomarker for Detecting Changes Occurring in the Human Microbiome. *Microb. Ecol.* **75**, 1063–1073 (2018).

36. Banerjee, S. et al. Network analysis reveals functional redundancy and keystone taxa amongst bacterial and fungal communities during organic matter decomposition in an arable soil. *Soil Biol. Biochem.* **97**, 188–198 (2016).

37. Berry, D. & Widder, S. Deciphering microbial interactions and detecting keystone species with co-occurrence networks. *Front. Microbiol.* **5**, 1–14 (2014).

38. Saito, A., Mitsui, H., Hattori, R., Minamisawa, K. & Hattori, T. Slow-growing and oligotrophic soil bacteria phylogenetically close to Bradyrhizobium japonicum. *FEMS Microbiol. Ecol.* **25**, 277–286 (1998).

39. Siering, P. L. et al. Microbial biogeochemistry of Boiling Springs Lake: A physically dynamic, oligotrophic, low-pH geothermal ecosystem. *Geobiology* **11**, 356–376 (2013).

40. Giovannoni, S. J., Schabtach, E. & Castenholz, R. W. Isosphaera pallida, gen. and comb. nov., a gliding, budding enbacterium from hot springs. *Arch. Microbiol.* **147**, 276–284 (1987).

41. Masson-Boivin, C. & Sachs, J. L. Symbiotic nitrogen fixation by rhizobia — the roots of a success story. *Curr. Opin. Plant Biol.* **44**, 7–15 (2018).

42. Fierer, N., Bradford, M. A. & Jackson, R. B. Toward an Ecological Classification of Soil Bacteria. *Ecology* **88**, 1354–1364 (2007).

43. Karimi, B. et al. Biogeography of soil bacteria and archaea across France. *Sci. Adv.* **4**, 1–14 (2018).

44. Lozupone, C. a & Knight, R. Global patterns in bacterial diversity. *Proc. Natl. Acad. Sci. USA* **104**, 11436–11440 (2007).

45. Terrat, S. et al. Meta-barcoded evaluation of the ISO standard 11063 DNA extraction procedure to characterize soil bacterial and fungal community diversity and composition. *Microb. Biotechnol.* **8**, 131–142 (2015).

46. Amore, R. D. et al. A comprehensive benchmarking study of protocols and sequencing platforms for 16S rRNA community profiling. *BMC Genomics* **17**, 55 (2016).

47. Zhou, J. et al. High-Throughput Metagenomic Technologies for Complex Microbial Community Analysis: Open and Closed Formats. *MBio* **6**, e02288–14 (2015).
48. Shade, A. et al. Fundamentals of microbial community resistance and resilience. *Front. Microbiol.* **3**, 1–19 (2012).

49. Tylianakis, J. M., Laliberté, E., Nielsen, A. & Bascompte, J. Conservation of species interaction networks. *Biol. Conserv.* **143**, 2270–2279 (2010).

50. Barberán, A., Bates, S. T., Casamayor, E. O. & Fierer, N. Using network analysis to explore co-occurrence patterns in soil microbial communities. *ISME J* **6**, 343–351 (2012).

51. Benjamini, Y. & Hochberg, Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. *J. R. Stat. Soc. Ser. B* **57**, 289–300 (1995).

## Acknowledgements

This study was granted by ADEME (French Environment and Energy Management Agency) and by “France Génomique” through involvement of the technical facilities of Genoscope (project number ANR-10-INBS-09-08). In addition, due to the involvement of technical facilities of the ANAEE France infrastructure at the GenoSol platform, it also received a grant from the French state through the National Agency for Research under the program “Investments for the Future” (reference ANR-11-INBS-0001). RMQS soil sampling and physico-chemical analyses were supported by a French Scientific Group of Interest on soils: the “GIS Sol”, involving the French Ministry of Ecology, Sustainable Development and Energy (MEEM), the French Ministry of Agriculture (MAP), the French Institute for Environment (IFEN), the Environment and Energy Management Agency (ADEME), the French Institute for Research and Development (IRD) and the National Institute for Agronomic Research (INRA). We thank all the soil surveyors and technical assistants involved in sampling the sites. Thanks to Christophe Djemiel for his comments on the results. Thanks are also extended to Diana Warwick for her comments on the manuscript.

## Author Contributions

All authors conceptualized the research project. C.J. coordinated the sampling at the territory scale. P.W. and C.C. contributed to the DNA sequencing. S.T. and S.D. performed the bioinformatic analyses. C.J. and D.A. provided environmental dataset. B.K. performed the data analyses and prepared the figures. B.K. and L.R. wrote the original draft. All authors reviewed the manuscript.

## Additional Information

Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-019-40422-y.

**Competing Interests:** The authors declare no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_5.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2019