Inter-Kingdom Networks of Canola Microbiome Reveal Bradyrhizobium as Keystone Species and Underline the Importance of Bulk Soil in Microbial Studies to Enhance Canola Production.

The subterranean microbiota of plants are of great importance for plant growth and health, as root-associated microbes can perform crucial ecological functions. As the microbial environment of roots is extremely diverse, identifying keystone microorganisms in plant roots, rhizosphere and bulk soil is a necessary step towards understanding the network of inuence within the microbial community associated with roots and enhancing its benecial elements. To target these hot spots of microbial interaction, we used inter-kingdom network analysis on the canola growth phase of a long-term cropping system diversication experiment conducted at four locations in the Canadian prairies. Our aims were: to verify whether bacterial and fungal communities of canola roots, rhizosphere and bulk soil are related and inuenced by diversication of the crop rotation system; to determine whether there are common or specic core fungi and bacteria in the roots, rhizosphere, and bulk soil under canola grown in different environments and with different levels of cropping system diversication; and to identify hub taxa at the inter-kingdom level that could play an important ecological role in the microbiota of canola. Our results showed that fungi were inuenced by crop diversication but not by bacteria. We found no core microbiota in canola roots but identied three core fungi in the rhizosphere, one core mycobiota in the bulk soil and one core bacteria shared by the rhizosphere and bulk soil. We identied two bacterial and one fungal hub taxa in the inter-kingdom networks of the canola rhizosphere, and one bacterial and two fungal hub taxa in the bulk soil. Among these inter-kingdom hub taxa, Bradyrhizobium sp. and Mortierella sp. are particularly inuential on the microbial community and the plant. To our knowledge, this is the rst inter-kingdom network analysis utilized to identify hot spots of interaction in canola microbial communities. cropping systems, with all rotation phases present each year. Our study used the canola phases of three crop rotation systems. The canola grown was the glufosinate-tolerant variety L241C. The three crop diversication treatments used in this study were: (1) monoculture of canola, (2) wheat-canola, and (3) pea-barley-canola (Table 1). These treatments were applied in a randomized complete block design with four blocks at each of four experiment sites.


Introduction
Plant subterranean microbiota have often been described as a "black box" [1][2][3][4]-a term referring to an inherent complexity of inner workings that makes a system di cult to grasp in its entirety. We can de ne plant subterranean microbiota as composed of different biotopes, notably: the root interior, the rhizosphere, and the bulk soil. The root and rhizosphere microbiota are highly in uenced by the plant. Rhizodeposits recruit root symbionts and shape the microbial community of the rhizosphere [5][6][7]. Microbial communities shaped by the plant in the rhizosphere and root interior can protect it against pathogens and enhance its growth and production [6,[8][9][10][11]. Bulk soil microbiota, while important to plant health, are less in uenced by the plant due to their distance from the roots. However, the bulk soil microbiota are the inoculum from which the plant Table 1 Selected crop diversi cation treatments from a long-term experiment established in 2008 at three different sites in the Canadian prairies [39]. The rotation phases examined in this study in 2018 are underlined.
Cropping . Crops were grown according to best management practices, as described in Harker et al., (2015). With the exception of the Lethbridge site, the growing season at all sites was characterized by more frequent rain in July just before sampling ( Figure S4).
Root, rhizosphere, and bulk soil samples were collected at the end of the canola owering period, at full ower after 50% of owers on the raceme had opened. This occurred in the fourth week of July 2018. Three to four plants within each plot were randomly selected and uprooted with a shovel. The shoots were removed, and roots were placed in plastic bags and brought to the laboratory on ice in a cooler. About 5 g of rhizosphere soil per plot was collected by gently brushing the roots. The brushed roots were then gently washed with sterilized distilled water. The bulk soil was sampled with a 2 cm diameter soil probe at a 7 cm depth, exactly in between two seed rows. The samples were kept at 4°C before being shipped on ice to the laboratory in Québec City, Québec, where they were preserved at -80°C until DNA extraction.
More details on site description, experimental design and sampling methods are provided in Floc

ASV determination and bioinformatic pipeline
Bioinformatics were used in QIIME2 environment version 2021.4 [40]. The bioinformatic pipeline used for the processing of our ITS and 16S rRNA gene sequences was DADA2 v1.18.0 [41]. First, we used Cutadapt 3.4 to remove the primer part of the ITS and 16S rRNA gene sequences with "minimum-length" at 50 and "p-error-rate" at 0.1. Then, we excluded the sequences with less than 220 bp with the command "-p-trunc-len", as the base quality of the sequences tended to diminish below that threshold in our data. Next, the amplicon sequence variant (ASV) table was calculated, and chimeras eliminated, resulting in a sequence length ranging from 250 to 253 nucleotides. ASVs were then identi ed using the naïve Bayesian classi er method on the databases SILVA and RDP, and the identities of ASVs of interest were veri ed manually using BLAST on the NCBI nt database. With the taxonomic resolution of the ITS and 16S RNA gene, it is generally not possible to identify a bacterium or fungus at the species level. Thus, the identi cations at species level presented here must be considered with caution despite the fact that they perfectly match (100% similarity and coverage) the NCBI reference sequences.

Data processing and statistical analyses
To assess the variation induced in the canola rhizosphere by crop diversi cation systems, the datasets were ltered from their rare ASVs using QIIME2 with "--p-min-frequency" set to 17 and "--p-min-samples" set to 1.
The effect of crop diversi cation on bacterial and fungal community structure was assessed by permutational multivariate analysis of variance (PERMANOVA) [42], considering 16 blocks (four blocks per each of the four sites), using the function "adonis" of the vegan package v2.5.7 [43] in R v4.1.0, and the set of relative abundance data.
We then aimed to identify the universal fungal and bacterial components of the core microbiota and hub taxa in each biotope of the canola subterranean microbiome. We de ned the core microbiota of a biotope as the set of microbial ASVs that were present in the microbiota of a particular biotope at all sites and plots. To assess the interactions among bacterial and fungal taxa in the microbiota, we created a co-occurrence inter-kingdom network using the package Spiec-Easi v1.1.0 in R 4.1.0 [21]. The analysis incorporated the roots, rhizosphere, and bulk soil fungal and bacterial community. The input data consisted of the raw abundance matrixes of the ITS and 16S ASVs. We rst ltered the datasets to remove ASVs with a frequency lower than 20%. The Spiec-Easi run was conducted with the algorithm "mb" with the lambda min ratio set at 10 − 2 and 50 repetitions. We then imported the networks into Cytoscape 3.8.2 [44] for plotting and used the "organic" layout to draw the networks. Edges were de ned as co-occurrences or mutual exclusion based on the positive or negative values of inverse covariance linking the nodes. Betweenness centrality, de ned as the fraction of the shortest path between all other nodes in the network containing the given node, and degree score, highlighted central nodes and provided information about network architecture. A score of betweenness centrality and degree of connectivity greater than 95% of the network taxa suggested participation in multipartite interactions in the community and allowed us to ag the highly connected taxa as hub taxa.

Results
Bioinformatic yield and taxonomic pro le of the microbiota of canola roots, rhizosphere, and bulk soil Sequencing yielded 10,118,613 ITS and 10,273,048 16S reads. Our bio-informatic pipeline retained 6,934,809 (68.03%) non-chimeric ITS-reads and 7,157,985 (68.21%) non-chimeric 16S reads after ltering, trimming, denoising and merging. These reads were respectively assigned to 2140 ITS amplicon sequence variants (ASV) and 9814 16S ASV. Rarefaction curves indicated that read abundances reached saturation for all samples ( Figure S1).
Principal coordinate analysis (PCoA) showed a site effect on both ITS and 16S communities ( Figure S2 and S3). PERMANOVA thus conducted by site reported a signi cant effect of cropping system diversi cation on fungal communities in certain sites and microbiota compartments, but no signi cant effect of diversi cation on bacterial communities was detected in any compartment of the subterranean microbiota of canola (  Total 11 1 DF, Degree of Freedom "*" = p < 0.05. "**" = p < 0.01 and "***" = p < 0.001, Level of signi cance.
Shared microbiota and Core-microbiota of canola roots, rhizosphere, and bulk soil.
The canola roots, rhizosphere, and bulk soil biotopes shared 318 fungal ASV. These ASV represent 59%, 24% and 21% of the fungal communities of the roots, rhizosphere, and bulk soil. Thirty-three fungal ASV representing respectively 5.6% and 2.5% of the root and rhizosphere fungal communities were shared in these ecospheres, but absent in bulk soil. Forty-two fungal ASV found only in the bulk soil and roots made up 7.2% and 2.7% of these fungal communities, respectively. Five hundred and fty-two fungal ASV representing 42% and 36% of the rhizosphere and bulk soil fungal communities, respectively, were absent in the roots (Fig. 3).
The 2405 bacterial ASV that were shared between the three compartments of the canola subterranean microbiota represented 69.4%, 29.9% and 29.2% of the bacterial communities of the roots, rhizosphere, and bulk soil, respectively. The 453 ASV shared exclusively between canola roots and rhizosphere soil represented 13% and 5.6% of the root and rhizosphere bacterial communities. The 268 ASV shared only between bulk soil and roots represented 7.7% and 3.2% of the roots and bulk soil bacterial communities, respectively. There were 4353 ASV shared exclusively between the rhizosphere and bulk soil. These shared ASV represented 54.2% and 53% of the rhizosphere and bulk soil bacterial communities, respectively (Fig. 4).
We were unable to detect a core microbiota in roots. However, three fungal amplicon sequence variants (FASV) were always present in the rhizosphere of canola, at all sites and regardless of crop rotation speci cations. This fungal core microbiota was composed of FASV1 (Tricocladium sp.), FASV4 (Fusarium sp.) and FASV7 (Cryptococcus sp.). Only one fungal ASV, FASV2 (Fusarium sp.), was present in the bulk soil sample of all plots. Only one bacterial amplicon sequence variant (BASV) was present in the canola rhizosphere sample of all plots. This BASV, BASV46 (Marmoricola sp.), was also the only BASV present in the bulk soil sample of all plots (Table 3). Inter-kingdom network analysis of the subterranean canola microbiota.
Network analysis detected no potential inter-kingdom interaction in canola roots; thus, no network was drawn.

Legends of the gures
In the bulk soil, the inter-kingdom co-occurrence network had the highest complexity, with 96 ASV and 149 edges (Fig. 6a). All three modules of the network were organized around 3 inter-kingdom hub taxa ( Table 4).

Discussion
Our results support the existence of stable core bacterial and fungal components of the canola rhizosphere and in the bulk soil of canola-producing elds. Inter-kingdom network analysis revealed hub taxa in rhizosphere and bulk soil, but no signi cant potential interaction was found in canola roots. Hub taxa BASV45 Canola subterranean microbiota compartmentation and response to crop diversi cation Crop diversi cation is known to increase canola yield by preventing the accumulation of pests and pathogens in soil [27,28,39,45]. We have previously shown the insigni cance of the crop diversi cation effect on the bacterial communities of the canola rhizosphere in the Canadian prairies [28]. In this paper, our results con rm that bacterial communities are insensitive to diversi cation of cropping systems not only in the canola rhizosphere, but also in its roots and bulk soil ( Table 2). This could be explained by the fact that bacteria are more in uenced by soil physical properties and weather conditions than by crop rotations [46].
A contrario, the composition of the fungal communities in canola roots, rhizosphere and bulk soil appeared to be sensitive to cropping systems diversi cation. The sensitivity of the rhizosphere fungal community was previously reported [27,36]. Geographic location and soil physical properties appear to affect fungal community composition as the crop rotation effect was site-dependant (Table 2).
Crop rotation systems are known to in uence microbial community structure in a wide range of crops [47][48][49].
Crop rotations can be used to modulate biological N-xation and to modify soil structure, with feedback on microbial communities [50][51][52]. Changing the fungal microbiota of canola with crop rotation systems could thus be useful for suppressing pathogens or enhancing canola nutrient uptake. The shaping of canola microbial communities must also take into account the environmental conditions, as fungal communities are subject to substantial geographical variations [53][54][55].
The dominance of Olpidiomycota in the taxonomic pro les of our fungal communities was similar to that previously reported in canola [27,28,36,56]. The dominance of Olpidiomycota in the fungal microbiota of the roots and rhizosphere is mainly due to the abundance of Olpidium brassicae [27,57]. This particularity of a single dominant species spread across the roots and rhizosphere is usually linked to situations of pathogen infestations [58-60] or symbiosis [61] in other crops. The predominance of O. brassicae was reported to have little in uence on canola yield [27,57].
In the soil, compartments like roots, rhizosphere and bulk soil are known to show a concentration gradient of plant chemical compounds that attract microbes [5,7,62]. Thus, the colonization of these three different niches share certain similarities in terms of the composition of their microbial communities. Cordero et al.
(2019) reported 77% similarity between bacterial communities in the roots and rhizosphere of canola. In our case, canola subterranean microbiota were similar in terms of the proportion of bacterial ASV shared between the different compartments, with a signi cant percentage of shared community between the bulk soil and rhizosphere (Figs. 3 and 4). The proportion of bacterial species shared between the root interior and rhizosphere of canola (82.4%) was similar to the proportion reported by Cordero et al., (2019). To the best of our knowledge, this is the rst report of the percentage of mycobiota shared between canola roots, rhizosphere, and bulk soil communities. We found the same trend in fungi as in bacteria. The fact that a signi cant part (30 60%) of the microbial communities are shared between these ecological niches could re ect their physical proximity. The fact that between the three biotopes, bacteria are shared more than fungi could be explained by the presence of lamentous fungi that allow bacteria to migrate following their hyphae, thus allowing a wide range of bacteria to navigate between the compartments.
Canola core-microbes in roots, rhizosphere, and bulk soil.
The canola root interior was devoid of any core fungi or bacteria. This lack of detection could be attributed to the fact that canola roots produce glucosinolates that are toxic to microbial life and lead to limited root colonization by fungi and bacteria [32,36,64]. Olpidium brassicae, an obligatory endophyte previously reported as a core fungi in canola rhizosphere in the Canadian prairies [27,36], was present and dominant in canola roots (Fig. 1A) but not enough to be agged as a core fungus according to our threshold of 99%. This difference between our results and results previously published can be explained by the difference in bioinformatics methods used. In this study, DADA2 allowed us to form ASVs (amplicon sequence variants) that are much more discriminating than the 97% identity threshold used in USEARCH to form OTUs. We can thus identify different genetic variants of the same species as different ASVs in DADA2 (Callahan et al., 2016).
This permits more precise inference of the structure and ecology of microbial communities. Olpidium brassicae shows signi cant genotypic variation and similarity with other pathogenic species known to affect canola, such as O. virulentus [57]. Its signi cant genetic variability and abundance in the roots and rhizosphere of canola suggest that O. brassicae should be the target of population genetics studies in the near future.
In the rhizosphere, three fungi and one BASV were detected as core microbes: FASV1 (Trichocladium sp.), FASV4 (Fusarium sp.), FASV7 (Cryptococcus sp.) and BASV46 (Marmoricola sp.). These three fungi were previously reported by [27] to be part of the canola rhizosphere fungal core microbiota. Paired with the observations reported in the previous article, the ndings of the present study illustrate the stability over time of the concept of core microbiota in the canola rhizosphere, reinforcing the need for long term studies with recurrent samplings.
In canola bulk soil, a core fungus and core bacterium were found: FASV2 (Fusarium sp.) and BASV46 (Marmoricola sp.), respectively. The latter was also present as a core bacterium in the canola rhizosphere. Regarding the former, fusaria are well known commensalists and pathogenic fungi, widely abundant in agricultural soils [65-67]. As no core microbiota has ever been identi ed in canola bulk soil, these results should be taken with caution. Core fungi and bacterium are subject to variation in presence and abundance over time and depending on weather conditions [27,28].
We were able to identify hub taxa at the intra-and inter-kingdom levels that are known as inter-kingdom hub taxa ( Table 4). Each of these could be of importance for canola production and manipulation of canola subterranean microbial communities.
BASV45 (Bradyrhizobium sp.) is a hub taxa that has been linked to other inter-kingdom hub taxa of the canola rhizosphere. Bradyrhizobium is a nitrogen-xing bacterial genus known to nodulate Fabaceae such as soybeans (Glycine max), cowpeas (Vigna unguiculata), Bambara groundnuts (Vigna subterranea) and chickpeas (Cicer arietinum) [68][69][70]. Furthermore, these bacteria demonstrate other ecological functions as plant growth promoting rhizobacteria (PGPR) through hormone secretions and antagonism in non-legume plants [10,71]. [72,73], reported Bradyrhizobium-induced nodular structures on canola roots. Thus, our detection of BASV45 (Bradyrhizobium sp.) as an inter-kingdom hub-taxon in the canola rhizosphere highlights this taxon as a potential PGPR for canola production, and as an agent for community manipulation in canola microbial networks. Indeed, high connectivity microbes are potentially bene cial for the plant, particularly in the rhizosphere [30,74]. Given how this taxon interacts with other plant species, ASV45 could be an important actor in the canola microbiome.
In the plant rhizosphere, microbes compete for space. One of the mechanisms used to compete with other microbes is the production of anti-microbial compounds [75][76][77]. Such is the case for the second most connected inter-kingdom hub taxon, BASV134 (Pseudonocardia sp.  [27,36,92]. FASV21, which was negatively linked to BASV45 (Bradyrhizobium sp.), and FASV151 (Exophiala sp.), is also a potentially bene cial organism for canola production. FASV151 (Exophiala sp.) was reported by [27] to be strongly positively linked to canola yield in the Canadian Prairies. Exophiala is a genus of fungi belonging to the dark septate endophytes (DSE) group, which hosts a broad range of plant growth promoting fungi [93,94].
This mutual exclusion between potentially mutualistic organisms could be explained by these microbes competing for similar ecological niches and may feed on similar rhizodeposits [95][96][97].
Bulk soil is an important part of plant microbiota, but very few studies of canola-related microbiota have taken bulk soil into consideration [45,[98][99][100]. Despite the fact that bulk soil is important to plant health and rhizosphere microbiota, most studies of canola-related microbiota have been restricted to root and rhizosphere microbiota. Bulk soil microbial communities emit volatile compounds that can enhance plant growth, protect against pathogens and even in uence plant root architecture [13,[15][16][17]. These microbial volatile compounds also have an in uence on the structure of rhizosphere microbial communities [14,101]. Compared to the rhizosphere microbial network, the bulk soil network showed a higher overall connectivity and a higher number of potential interactions between fungi and bacteria (Fig. 6). The bulk soil network had three modules, each centered on a hub taxon. In canola bulk soil, we were able to identify three inter-kingdom hub taxa, two fungi and one bacterium: FASV8 (Corynascella sp.), FASV114 (Mortierella sp.) and BASV69 (Bacterium sp.). None of these taxa were linked with each other. This nding, coupled with the modularity of the networks, suggests a strong ecological differentiation between the canola bulk soil hub taxa.
FASV114 (Mortierella sp.) was the most connected hub taxon in the network, with a dominance of mutual exclusions. It was negatively linked with Pseudonocardia and with FASV129 (Davidiella sp.). The latter genus, the teleomorph form of Cladosporidium, is known to count among its ranks numerous plant pathogens [102,103]. As we discussed previously, Mortierella can act as a PGPF. Its negative links with potential pathogens reinforce the need to investigate the impact this genus could have on canola production and health, as it is present as hub taxa in both the bulk soil and rhizosphere soil.
The potential ecological role of the two other hub taxa, FASV8 (Corynascella sp.) and BASV69 (unknown bacterium sp.), is rather di cult to attribute, as FASV8 (Corynascella sp.) is only reported in donkey dung from Iraq [104] and BASV69 is unknown. FASV8 was positively linked with DSE FASV151 (Exophiala sp.), which is known to be of interest for canola cropping. The position of these two hub taxa in bulk soil suggests their potential importance as microbes in canola production and thus the need for subsequent studies de ning their ecological functions.
The fungal microbiota in canola roots, rhizosphere and bulk soil respond to cropping system diversi cation, but show different responses depending on geography and weather conditions. The microbiota in canola roots, rhizosphere and bulk soil demonstrate ecological interconnectivity and recruitment, as a signi cant part of the microbial communities of these biomes is shared. This rst inter-kingdom network analysis of canola rhizosphere and bulk soil microbiota allowed identi cation of two particular microbes of interest for canola production: Bradyrhizobium sp. and Mortierella sp. The latter is a hub taxon in canola rhizosphere and roots and is linked to Exophila sp., a taxon previously described as associated with canola and positively correlated with high canola yield. The hub taxa matching Bradyrhizobium and Mortierella at the genus level could be of potential interest for bioengineering of canola subterranean microbiota and enhancing canola production.  Venn diagram of the ASV of the bacterial community shared between root, rhizosphere and bulk soil, taking all sites in account.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.