Rhizosphere Bacterial Networks, but Not Diversity, Are Impacted by Pea-Wheat Intercropping

Plant-plant associations, notably cereal-legume intercropping, have been proposed in agroecology to better value resources and thus reduce the use of chemical inputs in agriculture. Wheat-pea intercropping allows to decreasing the use of nitrogen fertilization through ecological processes such as niche complementarity and facilitation. Rhizosphere microbial communities may account for these processes, since they play a major role in biogeochemical cycles and impact plant nutrition. Still, knowledge on the effect of intecropping on the rhizosphere microbiota remains scarce. Especially, it is an open question whether rhizosphere microbial communities in cereal-legume intercropping are the sum or not of the microbiota of each plant species cultivated in sole cropping. In the present study, we assessed the impact of wheat and pea in IC on the diversity and structure of their respective rhizosphere microbiota. For this purpose, several cultivars of wheat and pea were cultivated in sole and intercropping. Roots of wheat and pea were collected separately in intercropping for microbiota analyses to allow deciphering the effect of IC on the bacterial community of each plant species/cultivar tested. Our data confirmed the well-known specificity of the rhizosphere effect and further stress the differentiation of bacterial communities between pea genotypes (Hr and hr). As regards the intercropping effect, diversity and structure of the rhizosphere microbiota were comparable to sole cropping. However, a specific co-occurrence pattern in each crop rhizosphere due to intercropping was revealed through network analysis. Bacterial co-occurrence network of wheat rhizosphere in IC was dominated by OTUs belonging to Alphaproteobacteria, Bacteroidetes and Gammaproteobacteria. We also evidenced a common network found in both rhizosphere under IC, indicating the interaction between the plant species; this common network was dominated by Acidobacteria, Alphaproteobacteria, and Bacteroidetes, with three OTUs belonging to Acidobacteria, Betaproteobacteria and Chloroflexi that were identified as keystone taxa. These findings indicate more complex rhizosphere bacterial networks in intercropping. Possible implications of these conclusions are discussed in relation with the functioning of rhizosphere microbiota in intercropping accounting for its beneficial effects.


INTRODUCTION
Rhizosphere is a dynamic zone of interactions between microorganisms and their host plants (Hiltner, 1903;Hartmann et al., 2008). These interactions can be assimilated to a feedback loop: plants release a significant part of their photosynthates in the form of rhizodeposits, which results in the recruitment of a microbial community best adapted to the rhizosphere environment; rhizosphere microorganisms interact with each other and with the host plant, and impact plant growth, nutrition and health (Philippot et al., 2013). Rhizosphere ecology has received a great deal of attention with major progress made in understanding plant-microorganism interactions (Philippot et al., 2013;Guttman et al., 2014). They have allowed to demonstrate the specificity of the socalled rhizosphere effect at the species (Lemanceau et al., 1995;Grayston et al., 1998;Berg and Smalla, 2009;Lakshmanan et al., 2014;Tkacz et al., 2020) and even at the genotype level, for maize (Peiffer et al., 2013), soybean (Zhong et al., 2019), and medic . The importance of the rhizosphere microbiota in terms of abundance, diversity and beneficial effects for the host plant has led to an holistic vision of the plant and its microbiota, rather than considering plants and microbiota as standalone entities (Hacquard and Schadt, 2015;Vandenkoornhuyse et al., 2015;Theis et al., 2016). Plant growth, development, health and fitness are mediated by plant but also microbial traits, with variations in plant phenotypes directly linked to their rhizosphere microbiota (e.g., biomass: Swenson et al., 2000;flowering time: Panke-Buisse et al., 2015). Thus, the holobiont concept has been recently proposed as encompassing the plant per se and its associated microbiota (Vandenkoornhuyse et al., 2015). Lemanceau et al. (2017) have further proposed the concept of functional core microbiota, in which plants recruit given microbial functional genes whatever the soils in which they are cultivated. Identification of plant and microbial traits involved in positive feedback loops has become a major target for plantbreeding in order to take better advantage of beneficial effects of rhizosphere microbiota (Wei and Jousset, 2017) for decreasing the use of chemical inputs in a more sustainable agriculture (Lemanceau et al., 2015).
Agroecology aims at valuing biotic interactions in agroecosystems in order to reduce the use of chemical inputs. A specific attention is given to crop diversification to promote agriculture sustainability (Altieri, 1999;Wezel et al., 2014;Bedoussac et al., 2015;Lemanceau et al., 2015). A classic strategy for increasing plant diversity in cropping systems is the intercropping (IC) that consists in cultivation of different plant species or cultivars on the same field and at the same time (Willey, 1979). Intercropping is a longstanding and widespread practice in low-input cropping systems throughout the world (Altieri, 1999;Knörzer et al., 2009;Maitra et al., 2021). Intercrop area represents 20-25% of arable land in China,and 17% in India,and up to 83% in Northern Nigeria,and 94% in Malawi (Knörzer et al., 2009). Indeed, in Europe, intercropping systems, such as wheat-pea, encounter different obstacles that contribute to their slow adoption and dissemination (Mamine and Farès, 2020).
Varietal selection is one of the main technical limit identified by authors (Mamine and Farès, 2020).
Intercropping may allow to increasing yields (Bedoussac and Justes, 2010a;Mamine and Farès, 2020), while reducing or even avoiding the use of nitrogen fertilizers when using legumes thanks to their ability to fix atmospheric nitrogen (Pelzer et al., 2012). Indeed, legumes promote the uptake of nutrients (e.g., nitrogen, phosphorus and iron, . . . ;Hinsinger et al., 2011;Zuo and Zhang, 2009;Xue et al., 2016) and grain protein content of the associated cereals (Bedoussac et al., 2015). Additionally, IC allows to: (i) reducing the pressure of weeds, by occupying available ecological niches, and that of pests, through the physical barrier effect (Corre-Hellou et al., 2011), and to (ii) providing mechanical support to the peas by the cereals (Bedoussac et al., 2015).
Nitrogen-fixing bacteria contribute to the above-referred beneficial effects of the associated legumes. More generally, it has been proposed that rhizosphere microbiota may account for the added value of IC. Thus, attempts have been to assess the impact of IC on rhizosphere microbiota. Total microbial communities (Taschen et al., 2017;Gao et al., 2019;Liu et al., 2021;Tang et al., 2021) and specific functional guilds (e.g., ammonia oxidizing bacteria, Song et al., 2007;diazotrophic Proteobacteria, Solanki et al., 2020) from the total root systems of plant genotypes associated in IC have been analyzed. These reports evidenced that rhizosphere microbiota from IC and sole-cropping (SC) differ significantly, these differences being more strongly expressed for bacteria than for fungi (Gong et al., 2019). Changes in bacterial communities were mostly associated with differences in the abundance of specific phyla. These phyla were either increased (e.g., Proteobacteria, Chloroflexi, Gemmatimonatedes, Acidobacteria, Nitrospirae, and Firmicutes in proso millet and mung bean; Rhizobiales, Burkholderiales, Pseudomonadales and Bacillus populations in wheat and alfalfa; Actinobacteria in wheat and pea) or decreased (Actinobacteria in proso millet and mung bean; Sphingomonadales and Xanthomonadales populations in wheat and alfalfa, α-Proteobacteria and Acidobacteria in wheat and pea) (Taschen et al., 2017;Gong et al., 2019;Li et al., 2020). However, considering separately the roots of each plant genotypes cultivated in IC, no consistent conclusion can be drawn. No difference could be detected between microbiota from IC and SC of fababean and wheat , or between rhizobia populations in IC and SC of maize and soybean (Herrmann et al., 2014). In contrast, abundance of ammonia oxidizing bacteria was increased in maize and fababean in IC compared to the respective SC (Song et al., 2007). Conclusion variations between reports could possibly be ascribed to differences in compatibility between plant genotypes cultivated together. Optimization of plant-plant interactions in intercropping by an appropriate choice of plant genotypes and cultivars is a major issue. It has been hypothesized that this choice is also crucial to value beneficial plant-plant-microbe interactions. This is supported by the different responses of the root bacterial community of two sugarcane varieties when intercropped with soybean (Liu et al., 2021). Thus, identifying the appropriate plant partners in IC represent a major issue, this require to test different combinations of plant genotypes/cultivars. Biotic interactions do not only occur between plant-plant and plants-microbes, but also among microorganisms. An increasing attention is given to interaction networks between rhizosphere organisms (van der Heijden and Hartmann, 2016) and the impact of IC on these networks has recently been stressed (Liu et al., 2021;Tang et al., 2021). Thus, co-occurrence networks have been proposed as an additional parameter to characterize microbial communities (Barberán et al., 2012;Berry and Widder, 2014;Morriën et al., 2017;de Vries et al., 2018;Jacquiod et al., 2020). Concerning IC co-occurrence network, Tang et al. (2021) showed, through a shotgun metagenome analysis, that sugarcane and peanut IC increased the abundance of bacterial genes involved in organic matter turnover comparing to SC, without correlating these differences to changes in microbiota diversity. Liu et al. (2021) showed differences between the co-occurrence networks between two sugarcane varieties in IC. None of these studies analyzed the differences in co-occurring network between bacterial taxa in IC and SC.These complex interactions may account for the increased yield of the IC but also of the following crop, indicating a positive legacy effect of multispecies cropping systems .
In the present study, we evaluated the impact of pea-wheat intercropping on rhizosphere microbiota. More specifically, we assessed how this impact may differ according to the pea and wheat cultivars. Biodiversity, structure and network of co-occurrence of bacterial community from roots of plants cultivated separately and in combination were characterized by high throughput sequencing of 16S rRNA genes. Results are discussed in the impact of IC of wheat and pea on the diversity, structure and networks of bacterial communities.

Site Description, Experimental Design, and Sampling
The impact of a given plant species on the rhizosphere bacterial community of the other plant species, cultivated in intercropping (IC), was tested in two independent experiments. Both experiments were performed at the Experimental Unit INRAE-Epoisses, France (47 • 14 11.2 N 5 • 05 56.1 E), the first from October 2016 to July 2017, the second from October 2017 to July 2018 (Table 1 and Supplementary Figure 1). Both followed a spring oat crop ( Table 1). Soil physico-chemical parameters are indicated in Table 1.
In the first experiment, emphasis was given to the impact of wheat on pea bacterial community, by testing the effect of seven winter wheat (Triticum aestivum) cultivars (CF11007, CF14336, Ehogold, Flamenko, Forcali, RE13003, Renan) on the bacterial community of three winter pea (Pisum sativum) cultivars (cv. Fresnel -hr genotype, Geronimo and Spencer -Hr genotypes).
In the second experiment, emphasis was given to the impact of pea on wheat bacterial community, by testing the effect of 11 winter pea cultivars (Aviron, China S-29, Fresnel, Furious, Geronimo, Isard, Isard H3 1.2, Isard ttl, Joker, IVD 304/10, Spencer) on the bacterial community of two winter wheat  Table 2).
In both experiments, each wheat and pea cultivars were cultivated in sole cropping (SC) and in IC.
The IC set up was full mixed of the two plant species on the row IC, as previously described to be the best suitable for cereals and herbaceous legumes in intercropping (Malézieux et al., 2009). Sowing rates varied according to experimental treatments as follows: (i) wheat in SC: 300 grains/m 2 , in IC: 150 grains/m 2 ; (ii) hr peas in SC: 80 grains/m 2 , in IC: 60 grains/m 2 ; (iii) Hr peas in SC and IC: 40 grains/m 2 . Plants did not receive any chemical inputs, or watering. The sowing rate was optimized in order to reach 50% wheat and 50% pea at harvest in IC. In all cases, sowing rate was at 50% for wheat, that of pea differed upon cultivars. It was at 75% and 100% for hr and Hr genotypes, respectively, to take in account the difference of competitiveness of the pea genotypes (Bedoussac and Justes, 2010b).
These sowing rates allowed a plant emergence, expressed as the average ratio between IC/SC, equal to 48% for the wheat and 87% for the hr genotype, and 47,5% for the wheat and 104,5% for the Hr genotypes in the first experiment, and 62% for the wheat, 89% for the pea hr genotype, and 78% for the Hr genotypes in the second experiment. Treatments were replicated in three blocks, each encompassing 31 plots (1.5 m × 8 m) in the first experiment and 35 in the second.
Ten root systems were randomly sampled per plot to a depth of 20 cm. In intercropping cultures, only wheat and pea root systems in close contact were sampled and their roots were further carefully separated on site. Samplings were performed at an early flowering stage for peas, which was reached 15 days earlier in hr than in Hr pea genotypes. Wheat roots were sampled at both these dates corresponding to heading stage. Bare soil was collected in three uncultivated plots integrated into blocks in each experiment and at each sampling date.
Root systems and bare soils were kept cold and transferred immediately to the laboratory. Rhizosphere soils were taken from the root systems as described by Offre et al. (2007). Samples of rhizosphere and bare soils were lyophilized at −80 • C and stored at −20 • C.

Molecular Characterization of Bacterial Communities
One hundred fifty-six samples of rhizosphere soil were analyzed from the first experiment (84 from the wheat rhizosphere and 72 from the pea rhizosphere) and 177 from the second (78 from the wheat rhizosphere, 99 from the pea rhizosphere), for a total of 333 rhizosphere samples. Moreover, a total of 12 bare soil samples were also analyzed as controls, for a total of 345 samples.
DNA was extracted from soil samples (1 g dry weight) according to ISO standard 11063 (Petric et al., 2011). The library for MiSeq sequencing was generated through two PCR steps according to Berry et al. (2011). The first step consisted in amplifying all the taxa present in the samples. The bacterial 16S rRNA gene V3-V4 hypervariable region was amplified using primers Pro341F (5 -CCTACGGGAGGCAGCAG-3 ) and Pro805R (5 -CCTACGGGNBGCASCAG-3 ) (Takahashi et al., 2014). The PCR1 mix was prepared by adding 10X Advantage 2 PCR Buffer (Ozyme, Saint-Cyr-l'École, France), 10 mM each dNTP Mix (Thermo Fisher Scientific, Waltham, MA, United States), 10 µM of Pro805R and Pro341F each (Eurogentec, Liège, Belgium), 1.5U of 50X Advantage 2 Polymerase Mix (Ozyme, Saint-Cyr-l'École, France), 10 ng of DNA to be amplified and water up to 25 µl of final volume for each tube. Thermal cycling conditions were 2 min at 94 • C, followed by 30 cycles of 30 s at 94 • C, 30 s at 55 • C, and 1 min at 72 • C, with a final extension at 72 • C for 1 min. Duplicate first step PCR (PCR1) products were pooled and then used as template for the second step PCR (PCR2). PCR2 amplification added multiplexing index-sequences to the overhang adapters using a unique multiplex primer pair combination for each sample according to Illumina guidelines. The conditions for PCR2 were the same as for PCR1, except for the number of cycles (8 cycles instead of 30).
Both PCR steps were performed on the Applied Biosystem 9700 thermal cycler (Applied Biosystem, Foster City, CA, United States). The PCRs were checked by electrophoresis (1.5% agarose, TAE1X, 100V). Two technical replicates of each PCR (1 and 2) were made, the products were then pooled and purified by Agencourt AMPure XP magnetic beads (Beckman Coulter, Brea, CA, United States), according to the supplier's recommendations. The amplified DNA was finally quantified at StepOnePlus Real Time PCR Systems (Thermo Fisher Scientific, Waltham, MA, United States). The equimolar mixture of the samples was prepared before being sent for sequencing (GenoScreen, Lille, France). Sequencing was performed using 300-bp paired-end sequencing chemistry on the Illumina MiSeq platform (Illumina, San Diego, CA, United States). Raw paired-end reads were then demultiplexed and assembled per sample, with the Illumina MiSeq Reporter software (version 3.1).

Bioinformatic Analyses
Sequence data were analyzed using an in-house developed Python notebook piping together different bioinformatics tools (available upon request). Briefly, quality checks of the 16S rRNA sequences were conducted using the QIIME pipeline (Caporaso et al., 2010b) and short sequences were removed (<400 bp). Reference-based and de novo chimera detection as well as clustering in Operational taxonomic Units (OTUs) were performed using VSEARCH (Rognes et al., 2016) and RDP representative set of 16S rRNA sequences as the reference database. The identity thresholds were set at 97%. Representative sequences for each OTU were aligned using PyNAST (Caporaso et al., 2010a) and a 16S rRNA phylogenetic tree was constructed using FastTree (Price et al., 2010). Taxonomy was assigned using UCLUST (Edgar, 2010) and the SILVA database (SILVA SSU 138 update release; Quast et al., 2012).
The raw sequences for this study have been deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB42023 1 .
After cleaning, a total of 5 824 759 sequences (a mean of 16 316 sequences for each sample) were kept for OTU picking. 10 549 OTUs were delineated, and 5 713 OTUs were considered for further analysis after rarefying using the "rarefy_even_depth" function in phyloseq package.

Statistical Analyses
Statistical analyses were conducted using R statistical software version 4.0.2 (R Development Core Team, 2014). The α-diversity of the bacterial communities was assessed by quantifying the number of OTUs per sample (richness), and by calculating the Shannon (both richness and evenness) and Simpson (evenness) indexes. β-diversity was investigated through Principal Coordinate Analysis (PCoA) of unweighted UniFrac distances. Both α and β diversity were calculated using phyloseq package (McMurdie and Holmes, 2013). Effects of plant species, cultivars, and associated plant species and cultivars on bacterial communities were tested using one-and two-way Permutational MANOVA analyses (PERMANOVA; this statistical test allows a direct additive partitioning of the variation for complex models), with 999 permutations and, when necessary, applying the 'strata' correction in order to restrict permutations between the 2 years of culture. Significant tests were further carried on individual pair-wise comparisons between experimental treatments, as described by Anderson (2001). PERMANOVA analyses were run using the Vegan package (Oksanen et al., 2019). OTUs explaining differences between treatments were identified by differential OTU abundance analysis. This was achieved by fitting a generalized linear model with a negative binomial distribution to normalized values for each of the OTUs and testing for differential abundance using a likelihood ratio test under DESeq package (Anders and Huber, 2010), after performing the extension DESeq with phyloseq (McMurdie and Holmes, 2014).
Wilcoxon-Mann-Whitney and Kruskal-Wallis tests were applied to identify possible significant effects on α-diversity of plant species and cultivars in monocropping and in intercropping. Kruskal-Wallis tests were followed by pair-wise comparisons with Dunn test. Wilcoxon-Mann-Whitney and Kruskal-Wallis tests were run under the dplyr package (Wickham et al., 2021) and Dunn test under the dunn.test package (Dinno, 2015).

Network Analyses
The interactions between coexisting OTUs in rhizosphere and in bare soils were further analyzed through co-occurrence network. As suggested by Berry and Widder (2014), data of the two experiments were pooled for constructing the corresponding matrices in order to increase the size of the sample dataset and thus obtain a more stable co-occurrence network and accurate correlation estimation. Before outputting the five matrices (bare soil, SC and IC wheat, SC and IC pea), only OTUs represented in at least 50% of the samples of the entire dataset were kept 1 https://www.ebi.ac.uk/ena/browser/view/PRJEB42023 for the co-occurrence network computation. Thus, OTUs only occurring in one of the two experiments were eliminated during this initial trimming, in order to keep only the common to the two experiments. The five correlation matrices amongst OTUs were calculated using Poisson Log Normal models (PLN, Chiquet et al., 2020). The models were validated by using the Bayesian Information Criterion (BIC, only r2-values provided here) and the significance of partial correlations was evaluated by a resampling of each matrix (n = 30) to test the robustness of the networks, using the Stability Approach to Regularization Selection (StARS) method (Liu et al., 2010). StARS method was developed for high dimensional problems and is based on random subsamples (30 iterations in our study, as stated before) and the construction of an highly stable graph from subsamples, in order to evaluate the robustness of the network along the path of solutions (Liu et al., 2010). Moreover, StARS method also allowed reducing possible biais in OTUs relative abundances between the two experiments, through the use of partial correlations. Hereafter we elaborated a network approach based on edge arithmetic (Jacquiod et al., 2020) to identify OTU correlations that were specific of intercropping (Supplementary Figure 2). Briefly, in the matrix of the pea and the wheat intercropping, we systematically removed all correlation interferences that were attributed to (i) the bare soil, (ii) the pea monocropping, and (iii) the wheat monocropping (Supplementary Figure 2). This resulted in two trimmed networks of OTUs showing specific links in the pea and in the wheat intercropping rhizosphere, respectively. Then, we intersected these two networks in order to retain the unique fractions of wheat, the unique fraction of pea, and the common conserved links found in both rhizospheres only in intercropping context. Complexity of networks was investigated by means of the degree index, the node betweenness and the edge betweenness (Newman, 2003).

Effect of Plant Species and Cultivars Grown in SC on the Rhizosphere Bacterial Communities
The effects of the plant species and cultivars on the rhizosphere bacterial communities were analyzed on either the separate or the pooled dataset of the two experiments.
Microbial α-diversity was significantly higher in wheat than in pea rhizosphere in the first experiment including seven wheat cultivars (Wilcoxon-Mann-Whitney, p = 7.46e-04; p = 3.3e-04; p = 1.27e-03 for Observed, Shannon and Simpson diversity indices, respectively), but not in the second only including two wheat cultivars (Supplementary Figure 3A). In the pooled dataset, microbial α-diversity was significantly higher in wheat than in pea rhizosphere (Wilcoxon-Mann-Whitney, p = 3.65e-06; p = 1.08e-05; p = 6.52e-05 for Observed, Shannon and Simpson diversity indices, respectively) (Supplementary Figure 3B). β-diversity was significantly different between wheat and pea rhizosphere in both experiments (separate datasets, PERMANOVA analysis; F-model = 2.61, p = 0.001, based on 999 permutations; pooled datasets, PERMANOVA analysis; F-model = 3.11, p = 0.001, based on 999 permutations) (Supplementary Figure 3C). Differences between the species in the pooled data sets were explained by five OTUs of Proteobacteria, three Alphaproteobacteria (two Rhizobiales and one Sphingomonadales orders) and two Betaproteobacteria (Burkholderiales order) that were preferentially associated with pea (Supplementary Data Sheet 1).
Within wheat species, no significant differences between cultivars were detected in microbial αand β-diversity in the first experiment (Supplementary Figures 4A,B). Within pea species, α-diversity was significantly higher in Hr than in hr cultivars in the second experiment (Supplementary Figures 5A,B; Wilcoxon-Mann-Whitney, p = 9.8e-04; p = 8.59e-04; p = 5.67e-03 for Observed, Shannon and Simpson diversity indices, respectively). β-diversity also differed significantly between

Compared Effects of SC and IC on the Diversity of the Rhizosphere Bacterial Communities
α-diversity of wheat bacterial communities did not differ between SC and IC with pea in any of the two experiments (Figure 1A), nor in the pooled dataset ( Figure 1B). Similarly, β-diversity of wheat bacterial communities did not differ in SC and IC ( Figure 1C).
αand β-diversity of pea bacterial communities did not either differ in SC and in IC in any of the experiments (Figures 1A,C), although significant differences were detected in the pooled dataset, in which the richness of the bacterial communities associated to IC pea plants was higher than the one associated to SC pea plants ( Figure 1B, Kruskal-Wallis, p = 2.35e-09, p = 5.26e-08, p = 3.05e-06, for Observed, Shannon and Simpson, respectively). In order to better explore the impact of the IC on the β-diversity of wheat and pea bacterial communities, a PERMANOVA with two covariate (plant species -wheat or peaand culture -SC or IC-) has further been performed. In all cases (separate and pooled experiments), this test confirmed that only plant species had a significant impact on β-diversity (PERMANOVA analysis for plant species factor; p = 0.001, based on 999 permutations).
No differences in αand β-diversity was either detected between intercropping and monocropping when testing a data subset only including cultivars which were shared in the two experiments (i.e., wheat: Ehogold and Fresnel, and pea: Fresnel, Geronimo and Spencer).

Effect of Mono-and Intercropping on the Co-occurrence Network of the Rhizosphere Bacterial Communities
Co-occurrence networks were produced from the pooled dataset, as a high number (>25) of samples is required to obtain a stable network with accurate correlation estimation, as recommended by Berry and Widder (2014).
The BIC R 2 (0.98 for monocropped and intercropped wheat, 0.99 for monocropped and intercropped pea, and 0.97 for bare soil) clearly showed that the PLN models fitted the dataset, resulting in accurate correlation matrices.
After removing the edges observed in the bare soil and in the SC rhizospheres (Supplementary Figure 2), we obtained cleaned intercropping networks featuring edges only present in the rhizosphere of each IC plant rhizosphere.
The resulting cleaned intercropped wheat microbial network consisted of 573 nodes (OTUs) and 1673 edges (1462 positive and 211 negative edges; mean degree or node connectivity 5.8). The mean node and edge betweenness centrality were 851.4 and 389.6 respectively. The resulting cleaned intercropped pea microbial network consisted of 451 nodes (OTUs) and 1189 edges (1112 positive and 77 negative edges; mean degree or node connectivity 5.3). The mean node and edge betweenness centrality were 716.4 and 357.1 respectively. Pea microbial network showed a higher positive to negative edge ratio in comparison to wheat network (14.4 vs. 6.9 respectively, Supplementary Figure 6A). Mean degree, node and edge betweenness were significantly higher in wheat than pea network (Wilcoxon-Mann-Whitney, p < 0.05; Supplementary Figures 6B-D). 50% of the OTUs in wheat and pea networks were affiliated to Acidobacteria, Bacteroidetes and Alphaproteobacteria phyla (Supplementary Figure 6E).
We then applied a network intersection (Supplementary Figure 2) between the pea and wheat networks under IC to specifically identify: (i) the unique fraction of the pea network that was only seen in the pea under IC; (ii) the unique fraction of the wheat network that was only seen in the wheat under IC; and (iii) the common network that was shared amongst both rhizosphere under IC.
The network of intercropped wheat was characterized by a dominance of Alphaproteobacteria, Bacteroidetes and Gammaproteobacteria OTUs (Figure 2A); but that intercropped pea did not show any significant taxonomic dominance ( Figure 2B).
Regarding the common network fraction shared amongst both pea and wheat rhizosphere under IC, a clear organization in modules was observed, with three keystone OTUs that had a very strong degree compared to the others. Networks belonging to intercropped wheat and pea shared three main modules including OTUs assigned to Acidobacteria, Alphaproteobacteria and Bacteroidetes phyla (Figure 3). Three keystone OTUs were further identified (Figure 3): OTU-496 belonging to Acidobacteria, OTU-152 belonging to Betaproteobacteria (order Burkholderiales, family Alcaligenaceae) and OTU-233 belonging to Chloroflexi (order Thermomicrobia).

DISCUSSION
An increasing attention is given to intercropping in agroecology to better value resources and to decrease the use of synthetic inputs (i.e., fertilizers, pesticides). In wheat-pea intercropping, reports indicated the promotion of nitrogen nutrition of wheat (Ghaley et al., 2005;Bedoussac and Justes, 2010a;Guiducci et al., 2018) and suggested a promotion of iron nutrition for pea (Zuo and Zhang, 2009), and of phosphorus nutrition for both plant species (Li et al., 2003Hinsinger et al., 2011). The possible contribution of rhizosphere microbiota of plants to these beneficial effects mostly remains to be untapped. This requires to disentangle the complex interactions between the plants grown in association and their rhizosphere microbiota. A major issue is to determine whether the rhizosphere microbiota of the plant species cultivated together differ or not from that of the plant species cultivated separately, and how the IC impact would vary upon the cultivars chosen to be cultivated together. For this purpose, we have compared the rhizosphere bacterial communities of wheat and pea when cultivated in intercropping and in sole cropping, and have tested different cultivar combinations. Bacterial communities were characterized on the basis of their biodiversity, structure and of co-occurrence network of 16S rRNA genes.
In sole cropping, we confirmed the specificity of the rhizosphere effect, which was reported for long (Lemanceau et al., 1995;Grayston et al., 1998;Berg and Smalla, 2009;Lakshmanan et al., 2014;Tkacz et al., 2020). Indeed, bacterial communities from wheat and pea differed significantly (Supplementary  Figure 3). Both richness and evenness of bacterial OTUs were higher in wheat than in pea rhizosphere, this greater number of bacterial OTUs is in agreement with previous reports (Turner et al., 2013;Taschen et al., 2017;Cordero et al., 2020). Differences in richness and evenness were at least partly explained by a higher representation of Proteobacteria in pea rhizosphere, especially Comamonadaceae, Sphingomonadaceae, but also, as expected, Bradyrhizobiaceae and Rhizobiaceae (Supplementary Data Sheet 1), known to be leguminous-associated bacteria (Wielbo, 2012;Chaudhari et al., 2020). Differences were less clear-cut at the cultivar level. In wheat, no difference could be detected between the rhizosphere microbiota of the tested cultivars (Supplementary Figure 4). This is in agreement with Corneo et al. (2016), Simonin et al. (2020). However, Kavamura et al. (2020) identified differences of richness between tall and semi-dwarf cultivars, which contrast with the lack of differences recorded in the present study between the four tall cultivars (RE 13 088, CF 14 336, RE 13 003 and Ehogold) and the three semi-dwarf (Flamenko, Forcali, and Renan) tested. In pea, we detected a significant higher richness of the bacterial FIGURE 3 | Co-occurring bacterial network of common OTUs belonging to intercropped wheat and pea rhizosphere. Each network node (individual circle) represents an OTU. Network edges are represented as straight lines connecting the nodes and indicate significant co-occurrences based on partial correlation obtained from a Poisson Log Normal model (r > | 0.03|, p < 0.05, n = 30 iterations); green for positives and red for negative co-occurrence. OTU-496, OTU-152 and OTU-233 have been identified as keystone OTUs.
Frontiers in Microbiology | www.frontiersin.org communities in Hr than in hr genotypes (Supplementary Figure 5). This observation is consistent with the physiological differences between the genotypes known for their differential sensitivity to photoperiod, that leads to a later floral initiation (Murfet, 1973;Alcalde et al., 1999), flowering and maturity in Hr genotypes. Still, to our knowledge, this is the first report pinpointing significant discriminating effect of Hr/hr genotypes on their rhizosphere bacterial communities.
In intercropping, biodiversity and structure of rhizosphere bacterial community of plant species/cultivars did not differ significantly from sole cropping in any of the two experiments (Figure 1). This is consistent with previous reports made on wheat-fababean IC  and on sugarcane varieties IC (Liu et al., 2021). Differences between intercropping and sole cropping were only detected in pea, when pooling data from the two experiments. Then, bacterial community richness appeared to be higher in intercropping than in sole cropping. Our observation is consistent with previous research on intercropping involving other plant species .
However, in overall, our results on biodiversity and structure are not in favor of a differential effect of IC compared to SC on wheat and pea bacterial communities. This would suggest that both plant species have similar impact on their bacteria independently of their neighboring plant. This is in agreement with Tkacz et al. (2020), who showed that bacterial community is influenced more by the root fraction than by the soil or plant species. However, the lack of differences between IC and SC could also be ascribed to the characterization methods of the bacterial community. Indeed, they provided information on the taxonomic composition and diversity, but not on the interactions between microbial groups or on their functions. Pivato et al. (2017) previously reported that despite their low impact on the total bacterial community, combination of plant species had a significant effect on the functional bacterial community mediating nitrification. Thus, additional analyses were performed to compare the co-occurrence networks in IC and SC. More specifically; we searched for the possible existence of specific cooccurrence links amongst rhizosphere OTUs that would only be recorded in IC. After filtering for all potential interference sources in our data (Supplementary Figure 2), we identified three networks whose edges were only observed in IC, only found in wheat cultivated in IC, another specific to pea cultivated in IC and finally a common network for wheat and pea cultivated in IC.
The specific wheat IC co-occurrence network was characterized by a dominance of OTUs belonging to Alphaproteobacteria, Bacteroidetes and Gammaproteobacteria (Figure 2). Alphaproteobacteria (e.g., Rhizobiales order) are known to be well represented in wheat rhizosphere (Bartoli et al., 2020), both in wild and domesticated cultivars, and thus to show a high heritability . Higher abundance of Bacteroidetes and Proteobacteria (e.g., Alphaproteobacteria, Gammaproteobacteria) was reported in tall than in semi-dwarf wheat cultivars (Kavamura et al., 2020). Connector bacterial OTUs belonging to Gammaproteobacteria, Alphaproteobacteria and Bacteroidetes were previously shown to be more represented in wheat rhizosphere than in bulk soil (Fan et al., 2018). The rationale for the further increase of these nodes recorded here in intercropping remained to be investigated. In the specific pea IC co-occurrence network, no dominant OTUs nodes were observed. Larger and stronger networks in bacterial communities were also described in pea when combining two cultivars (Horner et al., 2019). Lastly, wheat and pea IC specific networks had three common main modules with dominant OTUs belonging to Acidobacteria, Alphaproteobacteria and Bacteroidetes phyla. Three additional keystone OTUs were identified to be shared in networks of wheat and pea cultivated in intercropping: belonging to Acidobacteria (class RB41, order Ellin6075), Betaproteobacteria (order Burkholderiales, family Alcaligenaceae), Chloroflexi (order Thermomicrobia). Acidobacteria was described to be a key taxa in microbiota network associated with wheat (Kavamura et al., 2020), and among this taxa, order Ellin6075 to be part of the core microbiota of Brassica napus rhizosphere (Taye et al., 2020). A lot of attention has recently been dedicated to Acidobacteria in rhizosphere ecology (da Rocha et al., 2013;Kielak et al., 2016a,b;Kalam et al., 2020) with populations of Acidobacteria enriched in the rhizosphere (da Rocha et al., 2013;Kielak et al., 2016b). Genomic and metagenomic analyses allowed Kielak et al. (2016a) to predict a range of activities (e.g., the ability to attach roots thanks to exopolysaccharide production, promotion of plant iron uptake, indole-3-acetic acid production) in Acidobacteria with populations beneficial to the host-plant (Afzal et al., 2019). Acidobacteria have frequently been described as cooccurring with Proteobacteria, however it is not yet clear if this co-occurence stems from overlapping niches and/or from metabolic interactions. Alcaligenaceae have been identified as being associated with soil suppressiveness to soilborne diseases (Chapelle et al., 2016;Gómez Expósito et al., 2017). Abundance of Chloroflexi, known for their ability to oxidize nitrites, varies in wheat rhizosphere upon N addition (Ma et al., 2020). Since legumes may transfer ammonium from nodules to the surrounding soil and plants (Zhang et al., 2017), wheat intercropped with pea may benefit from an increased content in ammonium that would promote nitrifiers. Indeed, amoA genes appear to be more represented in maize-peanut intercropping than in maize monocropping. Among Chloroflexi, Thermomicrobia are dominant taxa in wheat rhizosphere (Latif et al., 2020) and key taxa in microbial network associated with tall wheat cultivars (Kavamura et al., 2020).
In conclusions, the present study shows that bacterial communities associated with wheat and pea differ between IC and SC, despite the lack of significant differences of their biodiversity and structure. Among the key taxa of specific of IC networks, some could be candidate promoting plant growth, nutrition and health. Our data also point out more complex networks within bacterial communities in the IC rhizosphere of wheat and pea, whereas their biodiversity and structure were not impacted. Co-occurring networks of plant microbiome were described to be more structured and complex in rhizosphere than in bare soil (van der Heijden and Hartmann, 2016). How this increased complexity may account for the beneficial effects of the intercropping on the plant growth and nutrition remains to explored. Still, recent studies clearly showed an increased function expression in belowground communities when the networks of co-occurence between populations were more complex, despite the lack of biodiversity variations (Morriën et al., 2017). The possible enhancement of functionalities in more complex microbial networks could be assumed to be related to modified activities of the populations when closely interacting. This hypothesis is currently being tested with the help of transcriptomic approaches in synthetic bacterial communities.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
BP, JB, and PL conceived the study. TG, NM, CL, and JB designed the field experiments. TG and JM performed the field experiments. BP and FD performed the sampling and were involved in the experiments in molecular biology. BP and AS performed bioinformatics and statistical analysis. BP and SJ conceived the analysis of the co-occurrence networks. SJ performed the analysis of the co-occurrence network and revised the manuscript. BP and PL wrote and revised the manuscript. All authors contributed to the article and approved the submitted version.