To explore the effects of co-inoculating Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 on soybean, no-inoculation (100% chemical fertilizer), one-inoculation (75% chemical fertilizer + Bradyrhizobium japonicum 5038), and two inoculation (75% chemical fertilizer + Bradyrhizobium japonicum 5038 + Bacillus aryabhattai MB35-5) treatments were conducted in black and cinnamon soils with a total of six samples (B0, B1, B2, C0, C1, and C2) marked in this study. The doses of chemical fertilizers were 100 kg N/ha, 100 kg P₂O₅/ha, and 100 kg K₂O/ha for no-inoculation treatment, and 75 kg N/ha, 75 kg P₂O₅/ha, and 75 kg K₂O/ha for one-inoculation and two inoculation treatments. The chemical fertilizer used in this study was composed of urea (CON₂H₄), Ca(H₂PO₄)₂⋅H₂O and K₂SO₄.

The Bradyrhizobium japonicum 5038 (CGMCC16747) and Bacillus aryabhattai MB35-5 used in this study were screened in our laboratory. Bradyrhizobium japonicum 5038, a nitrogen-fixing strain that can tolerate water stress, was screened from 160 purified strains of Bradyrhizobium that had previously been isolated from the root nodules of different Asiatic soybean cultivars. This research was published by Zhu et al. (2022). Bacillus aryabhattai MB35-5 was screened from 204 purified strains of bacteria, which had previously been isolated from soils in different rice and wheat fields, that were mostly located in 16 Chinese cities. This strain is resistant to salt, alkali, acid and anaerobic conditions. It increased the root dry matter of corn and maize by 24.02 and 16.21%, respectively, compared with the control Silicate bacteria 3016. These results have already been used to obtain a Chinese invention patent (CN112772678B).

Single colonies of Bradyrhizobium japonicum 5038 were cultured in 100 ml liquid yeast extract mannitol broth (YMB) containing 10.0 g mannitol, 0.2 g MgSO4⋅7H₂O, 0.1 g NaCl, 0.5 g K₂HPO₄, 0.5 g yeast extract, 1000 ml H₂O (pH = 7.0). Single colonies of B. aryabhattai MB35-5 were cultured in 100 ml liquid R2A containing 0.5 g yeast extract, 0.5 g tryptone, 0.5 g glucose, 0.3g K₂HPO₄, 5.0 g MgO₃Si, 1000 ml H₂O (pH = 7.0). The cells were incubated at 30°C with shaking at 180 rpm until reached to the value of OD₆₀₀ for 0.6–0.8. The cells were incubated at 30°C with shaking at 180 rpm until they reached OD600 levels of 0.6–0.8. For the one-inoculation treatment, 3 ml of active culture of Br. japonicum 5038 was mixed with 1 kg of soybean seeds. When both strains were co-inoculated simultaneously, 1.5 ml of each activated bacterial liquid was mixed with 1 kg of soybean seeds. All the inoculated seeds were air-dried 2–3 h in the shade before planting.

Black soil samples were collected from a trial site located in Songyuan, Jilin Province (N 45°13′, E 124°82′), whereas cinnamon soil samples were collected from a trial site located in Shijiazhuang, Hebei Province (N 38°05′, E 114°50′). These sites are two experimental stations belonging to China’s agricultural research system. Soybean cultivars were Jiyu47 and Jidou12 for black soil and cinnamon soil, respectively. Soil samples were collected in September 7 and August 22 of 2020 at soybean early pod stage. Rhizosphere soil at a depth of 10–20 cm (plow layer) was collected, gently shaking off the soil which loosely adhered to the roots, retaining the soil tightly adhered to roots, collecting the soil by brush, and mixing to obtain one composite soil sample per plot. Then, we randomly selected five samples from each plot for the different treatments, and a total of 30 samples were collected. After the collection, all soil samples were stored in an icebox and transferred to a laboratory for examination. Each sample was separated into two parts: one part was stored at −80°C for biological analysis and the other part was air-dried in the shade for physical and chemical analyses.

The soil nitrate-nitrogen (NN) and ammonium-nitrogen (AN) contents were determined using samples of fresh soil, whereas the other soil properties were analyzed after drying soil samples at 25 ± 5°C and passing them through a 1.0 mm mesh. The soil pH was measured using a pH meter with soil to water ratio of 1:1 (Li et al., 2013). The NN and AN contents were determined via flow injection analysis (Mulvaney, 1996), whereas the soil total nitrogen (TN) and organic carbon (TC) levels were determined according to the procedure presented by Strickland and Sollins (1987). The alkali-hydrolyzable nitrogen (HN) was measured by an alkaline extraction method (alkaline potassium permanganate method). The available phosphorus (AP) was quantified using the NaHCO₃^– extraction method, whereas the available potassium (AK) was quantified using the NH₄OAc extraction method and a flame photometer (FP640, Shanghai, China) (Mulvaney, 1996). In addition, the soil metallic elements (Na, Ca, and Mg) were measured according to Acero et al. (2006).

When the soil sample was collected, at the same time, five completely plant roots of each treatment was collected and transported to the lab. Then, we selected the relatively full and suitable size nodules for counting and then weighted the dry weight of the nodules at 80°C until we reached a constant weight. The increasing numbers and dry matter of nodules that enable the inoculation of different PGPR species are illustrated in Supplementary Figure S1.

We extracted total DNA from rhizosphere soil with a Power Soil DNA Isolation Kit (MOBIO Laboratories Inc., Carlsbad, CA, United States) according to the manufacturer’s instructions. The DNA quality and concentration (A260/A280) were estimated using a NanoDrop ND-1000 UV-Vis Spectrophotometer (Thermo Scientific, Rockwood, TN, United States). Each DNA sample was stored at −80°C until further analysis. Subsequently, bacteria were analyzed by sequencing the V3-V4 hypervariable region of the 16S ribosomal RNA (rRNA) gene. The V3-V4 region was amplified using universal primers, 338F and 806R (338F: 5′-ACTCCTACGGGAGGCAGCA-3′ and 806R: 5′-GGACTACHVGGGTWTCTAAT-3′) (Xu et al., 2016). The PCR mixtures contain 5 × TransStart FastPfu buffer 4 μL, 2.5 mM dNTPs 2 μL, forward primer (5 μM) 0.8 μL, reverse primer (5 μM) 0.8 μL, TransStart FastPfu DNA Polymerase 0.4 μL, BSA 0.2 μL, template DNA 10 ng, and finally ddH₂O up to 20 μL. The hypervariable regions of the fungal 18S rRNA gene were amplified with the 817F and 1196R primers (817F: 5′-TTAGCATGGAATAATRRAATAGGA-3′ and 1196R: 5′-TCTGGACCTGGTGAGTTTCC-3′) using a thermocycler polymerase chain reaction (PCR) system (GeneAmp 9700, ABI, California, United States) (Rousk et al., 2010). The PCR mixtures contain 10 × buffer 2 μL, 2.5 mM dNTPs 2 μL, forward primer (5 μM) 0.8 μL, reverse primer (5 μM) 0.8 μL, rTaq Polymerase 0.2 μL, BSA 0.2 μL, template DNA 10 ng, and finally ddH₂O up to 20 μL. The PCR amplification of 16S rRNA gene was performed as follows: initial denaturation at 95°C for 3 min, followed by 30 cycles of denaturing at 95°C for 30 s, annealing at 55°C for 30 s and extension at 72°C for 45 s, and single extension at 72°C for 10 min, and end at 4°C.

The PCR products were utilized to purify and generate the amplicon libraries; subsequently, the obtained data were analyzed using QIIME pipeline version 1.8.0. We excluded low-quality reads in an initial quality-filtering step. Chimeric sequences were then identified using UCHIME. Operational taxonomic units (OTUs) were defined via clustering at 97% similarity. The rarefaction and diversity indices were calculated after clustering the OTUs. The taxonomy of each 16S rRNA gene sequence was analyzed by the ribosomal database project classifier¹ against the silva138/16s_bacteria 16S rRNA database and Unite8.0/its_fungi as references with a confidence threshold of 70% (Amato et al., 2013). Data were analyzed using a free online platform, Majorbio I-Sanger Cloud Platform². The raw data were deposited in the short reads archive database of the National Center for Biotechnology Information (SRA accession: PRJNA788980, PRJNA788360).

The network analysis was designed based on all the 30 samples at the phylum and class levels of bacteria and fungi. To explore all the pairwise associations, correlation scores (Spearman’s correlation), OTUs with relative abundances less than 0.005% or in less than five soil samples were discarded. The relationship among OTUs was examined by Pearson’s correlation, Spearman’s correlation, Bray–Curtis dissimilarity, and Kullback–Leibler dissimilarity using R (version 3.6.3). The P-values of the five methods were integrated using the Brown method, and the correlations | P| > 0.7 were retained for the downstream procedure. The resulting correlations were imported into the Gephi platform and then visualized by the Fruchterman–Reingold algorithms.

To determine the potential importance of stochastic processes on community assembly, we adopted a neutral community model (NCM) to predict the relationship between OTU detection frequencies and their relative abundance across the wider metacommunity, which were performed using R (version 3.6.3) and the program presented by Chen et al. (2019). The assembly processes of bacterial and fungal communities were evaluated by calculating the nearest taxon index and beta nearest taxon index (betaNTI) using the “ses.mntd” function in a picante package (Stegen et al., 2012; Xun et al., 2019). When | betaNTI| < 2, the contribution was considered a stochastic process, and when | betaNTI| > 2 the shifts in community composition were deterministic processes.

Soil physicochemical property data were analyzed using Microsoft Excel 2010 pro and SPSS version 20. A Mantel test was adopted to determine the significance of the relationship between soil physicochemical properties and community structure using QIIME. A random forest model was used to identify the bacterial and fungal taxa on the order level that accurately predict the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 by adopting the random forest package (importance = TRUE, proximity = TRUE, ntree = 1000). We further adopted the variation partitioning analysis (VPA) with adjusted R² coefficients based on the redundancy analysis (RDA) and partial Mantel test to quantify the relative effects of environmental and spatial factors in shaping community composition (R version 3.3.1).

An analysis of similarities indicated the inoculation treatments that produced substantial effects on the diversity of bacterial (R² = 0.76, p < 0.01) and fungal (R² = 0.78, p < 0.01) communities relative to the OTU level (Figure 1). The Chao and Shannon indices of the three treatments were all higher in cinnamon soil than that in black soil, respectively. In black soil, the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 increased the Shannon indices of bacteria comparing with that of the Bradyrhizobium japonicum 5038 mono-inoculated treatment. In cinnamon soil, the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 decreased the Chao indices of fungi comparing with that in Bradyrhizobium japonicum 5038 mono-inoculated treatment (Supplementary Figure S3).

A total of 7592 rarefied OTUs from 40 phyla were identified for bacteria. The top 10 most abundant phyla were Actinobacteria (37.37, 15.84%), Proteobacteria (21.03, 25.30%), Acidobacteria (12.68, 26.82%), Chloroflexi (14.31, 7.14%), Bacteroidetes (2.34, 6.71%), Gemmatimonadetes (3.15, 3.53%), Myxococcota (1.77, 4.07%), Firmicutes (2.61, 1.52%), Patescibacteria (1.13, 1.59%), and Methylomirabilota (0.59, 1.49%) in the black and cinnamon soils, respectively (Figure 2). To identify the most important predictors for shaping bacterial communities, 13 soil physicochemical properties were measured (Supplementary Table S1). The Mantel tests demonstrated that the bacterial communities were significantly (p < 0.05) correlated with TC, HN, and Ca for black soil and TC, TN, NN, AN, Ca, and pH for the cinnamon soil (Supplementary Table S2). Among these variables, HN (Mantel r = 0.36, p < 0.05) and NN (Mantel r = 0.58, p < 0.05) exhibited the strongest correlation with the bacterial community structure for black and cinnamon soil, respectively. Principal coordinates analysis (PCoA) results indicated that the bacterial communities were divided along the soil type and number of bacteria inoculation in axes 1 and 2, which could explain 72.21 and 4.68% of the total variations, respectively (Figure 1). The relative abundance of most of the top 10 phyla was significantly (p < 0.05) correlated with the main soil physicochemical properties. Among the soil physicochemical properties, pH and AK exhibit the significant correlation (p < 0.05) with more bacterial communities in black soil (Supplementary Table S4); TN and Mg exhibit the significant correlation (p < 0.05) with more bacterial communities in the cinnamon soil (Supplementary Table S5). Random forest analysis was adopted to identify biomarkers at the family level that are strongly associated with the PGPR species (Supplementary Figure S4). We determined that several of them were enriched by the increase of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5, including Pezizaceae, Aspergillaceae, Ascodesmidaceae, and Rhizophlyctidaceae.

A total of 2405 rarefied OTUs from 16 phyla were identified for fungi. The top three most abundant phyla were Ascomycota (85.38, 74.54%), Basidiomycota (8.23, 10.54%), and Mortierellomycota (3.97, 8.63%) in the black and cinnamon soils, respectively (Figure 2). Moreover, the same phenomenon emerged in the two soil types, such that with the increase in the species of inoculated bacteria, Ascomycota first decreased and then increased, whereas Basidiomycota increased first before decreasing. In other words, the inoculation with Bradyrhizobium japonicum 5038 reduced Ascomycota and increased Basidiomycota, and then the addition of Bacillus increased Ascomycota and decreased Basidiomycota. Mantel tests demonstrated that the fungal communities were significantly (p < 0.05) correlated with TC, NN, AN, AP, TK, Na, Ga, and Mg for the black soil and with TC, TN, NN, AP, AK, Na, Ca, and pH for the cinnamon soil (Supplementary Table S2). Among these variables, NN (Mantel r = 0.48, p < 0.05; r = 0.58, p < 0.05) exhibited the strongest correlation with the fungal community structure for both black and cinnamon soils. PCoA results showed that the fungal communities were divided along the soil type and number of inoculation of bacteria in axes 1 and 2, which could explain 58.60 and 7.82% of the total variations, respectively (Figure 1). The relative abundance of most of the top three phyla was significantly (p < 0.05) correlated with the main soil physicochemical properties. Among the soil physicochemical properties, HN and TN are more significantly correlated with the microbial communities in the black soil (Supplementary Table S4); pH and Na have a more significant correlation with the microbial communities in the cinnamon soil (Supplementary Table S5). Rhizophlyctidales, Auriculariales, Diversisporales, Capnodiales, and other fungi in the order level were analyzed as biomarkers that were strongly associated with the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 (Supplementary Figure S4).

In the co-occurrence network constructed by all 30 soil samples, 124 nodes and 1231 edges were identified for soil bacteria, whereas 80 nodes and 830 edges were identified for soil fungi (Figures 3, 4). In total, 12 sub-network analyses were also performed, as illustrated in Figures 3, 4 for bacteria and fungi, respectively, which showed that as the Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 inoculated, the nodes and edges increase. In other words, the inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 could increase the complexity of the co-occurrence pattern for both bacteria and fungi. The rich microorganisms in the sub-network were different in the two soil types, and the complexity of the co-occurrence pattern in the cinnamon soil was higher than that in the black soil, which is especially reflected in the increase in the amount of Proteobacteria and Dothideomycetes.

In total, 21 and 3 key species of bacteria and fungi within the network were identified (Supplementary Figure S5), which included Gymnoascus, Cyphellophora, Azospirillales, and several other microorganisms (Supplementary Table S6). Moreover, we identified several key species as biomarkers, Cyphellophora (R² = 0.484, P = 0.002) in the black soil and Solirubrobacteraceae (R² = 0.379, P = 0.009) in the cinnamon soil, which was exhibited a strong negative association with PGPR numbers, and Leucoagaricus (R² = 0.425, P = 0.005) in the cinnamon soil, which exhibited a strong positive correlation with PGPR species (Supplementary Figures S6, S7).

The relationship between network characteristics and the physical and chemical properties of soil is illustrated in Supplementary Figure S8. HN, TN, and TK were significantly positive with the number of nodes in 12 sub-networks, whereas the Ca and Mg contents were also important factors with a significant negative relationship with the number of bacteria and fungi nodes, respectively. In addition, environmental factors via the degree of bacteria network and average path length of the fungi network affected the co-occurrence pattern.

The NCM successfully estimated a large fraction of the relationship between the occurrence frequency of OTUs and their relative abundance variations (Figure 5), with 73, 71, and 72% of explained bacterial community variance for non-inoculation, one-inoculation, and two co-inoculation treatments, respectively. The obtained results indicated that stochastic processes were crucial in shaping the bacterial community assembly, and the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 can increase the stochastic processes for driving the microbial community assembly. The contribution of the stochastic and deterministic processes with the inoculation of the Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 was also in agreement with these results (Figure 6). Although the relationships between the occurrence frequency of OTUs and their relative abundance variations were 35, 30, and 37%, we acknowledge that the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5increased the stochastic processes for the fungi community assembly. The contribution of the deterministic processes in fungi was higher than those in bacteria, which also indicates that the deterministic process plays a more significant role in driving the fungal community assembly (Figure 6). The Nm-value was lower for bacterial taxa with the one-inoculation treatment of Bradyrhizobium japonicum 5038 (Nm = 26 170) and two co-inoculation treatment of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 (Nm = 26 858) than the non-inoculation treatment (Nm = 27 825). This phenomenon was more significant for the fungi community, and the Nm-value decreased to the one-inoculation treatment of Bradyrhizobium japonicum 5038 (Nm = 3961) and two co-inoculation treatments of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 (Nm = 3514) from non-inoculation (Nm = 5902). These results indicate that the species diffusion decreased after inoculating PGPR.

Recently, the co-occurrence pattern in ecology has been widely adopted to explore the ecological interactions among microbial communities across various ecosystems (Williams et al., 2014; Fan et al., 2021; Gao et al., 2021). In this study, the obtained results indicated that the network of the microbial community in two soil types after inoculation of Bradyrhizobium japonicum 5038 became complex, and some functional microorganisms increased, e.g., the nitrogen-fixing Azospirillales (Supplementary Table S6). Several previous studies also presented similar results, which indicated that in addition to simulating the proliferation of potential beneficial microbes, Bradyrhizobium japonicum 5038 inoculation also increased connections in rhizobacterial networks and altered the hub taxa (Zhong et al., 2019; Shang et al., 2021). In another study, the results indicated that compared to grass treatments, legumes increased soil bacterial diversity and the abundance of nitrogen-fixing Nitrospirae groups (Zhou et al., 2017b). For fungi, Bradyrhizobium japonicum 5038 inoculation in the soybean also facilitated an increase in the number of connections between fungi, including the hub fungi shift to Phaeosphaeria, Sporobolomyces, Septoria, Edenia, and Leptospora (Xu et al., 2020). Moreover, by building mutants with a defective Nod factor, previous research confirmed that the symbiosis compromised NoeI mutant-inoculation could reduce microbial diversity and co-occurrence interactions, and decrease the abundance of beneficial microbes (Liu et al., 2020). This could be because the inoculated Bradyrhizobium japonicum 5038 forms a more symbiotic system with the host plant, which was changed by altering plant root exudates such as flavonoid (White et al., 2015), Catechin (Mommer et al., 2016), and some volatile organic compounds to change the soil physical and chemical environment and influence the prevalence of rhizosphere microorganisms (Vishwakarma et al., 2020).

It is widely acknowledged that Bacillus produces highly resistant spores that can survive for prolonged periods in the soil, which are very promising as potential inoculants in agriculture (Kumar et al., 2016; Ferreira et al., 2020; Oleńska et al., 2020). In this study, the nodes and edges of the co-occurrence networks increased with the PGPR species, thus indicating that the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 leads to networks with more complexity than the single inoculated Bradyrhizobium japonicum 5038 (Figures 3, 4). Previous studies indicated with an increase in the diversity of microorganisms, the interaction relationship becomes stronger, which increases plant tolerance to a range of other abiotic stresses such as heavy metals contamination, water stress (flooding and drought), salinity, and cold (Gerhardt et al., 2009; Ju et al., 2020). After introducing Bacillus, pathogenic Fusarium certainly decreased to 0.89% from 1.24% in the individually inoculated Bradyrhizobium japonicum 5038 treatment in our study. Several studies verified that Bacillus spp. could produce 1-aminocyclopropane-1-carboxylate (ACC) to inhibit the stress-related hormone production in the plants and increase biotic stress resistance (Saleem et al., 2007; Hayat et al., 2010). Given the important role of environmental factors in driving microbial communities, the increase in the number of PGPR species can specifically provide microbial populations with more interactions owing to niche differentiation. In addition, keystone species are highly connected with other species within the networks, and they potentially exert a considerable influence on the entire microbial community (Gao et al., 2021). The primary reason for such effect is that different co-inoculated PGPR species can secrete different enzyme systems, such as iron transporter IRT1 (Zhang et al., 2009) and cell wall degrading enzyme (Bell et al., 2019), to increase the nutrient utilization rate. Increasing the supply of nutrients will increase the ecological niche of microorganisms, reduce competition, and allow more microorganisms to participate in the biochemical process. From the random forest analysis, our study also verified that Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 are the main microorganisms that alter the community structure.

Soil is a natural medium for cultivating soil microorganisms, and changes in its physical and chemical properties will trigger changes in the structure and interaction of the microbial community. The physical and chemical properties of different soil types will vary significantly, including the microbial communities (Gowda et al., 2021; Zhu et al., 2021). In the black soil, HN and Ca are the main influencing factors of the bacterial community, and NN and Na are the main influencing factors of the fungal community; in the cinnamon soil, the effects of TC, and NN are more significant in determining the microbial community (Supplementary Tables S2, S3). Interestingly, the same phenomenon emerged after co-inoculating Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 in the two soil types. In other words, the complexity of the network increased after inoculating the PGPR, and as the species of PGPR increased, the interaction between microorganisms also increased. Hence, the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 has universal and positive effects on the soil microbial interaction. Collectively, the inoculation of beneficial microorganisms plays a crucial role in maintaining soil ecological health and increasing its resistance. Compound rhizobia agent positively influences plant growth and represents a promising sustainable trend to increase plant production.

In this study, we examined bacterial and fungal community assembly in agroecosystems across two major soybean production areas throughout north China. Elucidating the relationship between community assembly and species coexistence is fundamental and crucial for understanding ecosystem diversity and functioning. Previous reports indicated that microbial co-occurrence associations tended to be higher when communities were primarily driven by the dispersal limitation relative to species sorting (Jiao et al., 2020). Our findings are in agreement with the finding that the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 increased the diversity of bacterial and fungal communities and the complexity of the co-occurrence pattern, at the same time, species diffusion in the community assembly process decreased after inoculating the PGPR species.

The deterministic processes refer to local environmental factors and biotic interactions that are responsible for shaping the microbial community, whereas stochastic processes consider birth, death, speciation, and immigration, which play a key role in shaping the microbial community assembly (Sloan et al., 2006; Zhou and Ning, 2017). In our study, the changes in the stochastic and deterministic processes in the bacterial and fungal communities differed within the inoculation of PGPR. The main assembly processes were stochastic processes for bacteria in the two ecological models, whereas the deterministic process kept increasing for fungi. These results could be attributed to the diversity of the microbial community. Xun et al. (2019) demonstrated that deterministic and stochastic assembly processes are dominant in low- and high-diversity communities, respectively. The bacterial communities exhibit more diversity (Chao index of 4002-4401) than fungi communities (Chao index of 604–814); thus, stochastic processes are dominant for bacteria. Although a previous study demonstrated that deterministic processes play a dominant role in bacterial community assembly in soybean fields (Zhang et al., 2018), the assembly processes of soybean rhizobacterial communities differed in different soybean cultivars and soil types (Liu et al., 2019; Zhong et al., 2019). The inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5s also exerts less impact on the assembly process, because bacteria have more species in sufficient ecological niches. Hence, although the interaction of bacteria increased with the PGPR species (Figure 3), the contribution of the biotic factor was limited. The environmental factor might be the main driving factor for the changes in the bacterial community. In addition, these speculations are well supported by the VPA, which indicated that the environmental factor could explain 86% of the changes in the bacterial community (Supplementary Figure S9).

The contributions of the deterministic assembly process were higher for fungi than bacteria. The VPA demonstrated that the environment factor explains 48% of the changes in the fungal communities, and the 52% residuals, probably owing to the increasing interaction following the co-inoculation of PGPR (Supplementary Figure S9). The inoculation of Bradyrhizobium japonicum 5038, due to the increase in the deterministic assembly process (75.60%) (Figure 6), indicated that the environmental factor could exert a more significant effect; thus, the resistance of the fungal community may decrease. When Bacillus was added, the contribution of the deterministic assembly processes decreased to 33.30%, whereas the contribution of the stochastic processes increased to 66.70% (Figure 6). Bacillus further increased the diversity of fungi, could replenish the ecological niche in time, and competitively prevent the invasion of pathogenic microorganisms (Ongena and Jacques, 2008; Cawoy et al., 2011). Furthermore, the added Bacillus further limited the species diffusion and restricted the direction of movement of pathogenic microorganisms. The shift of balance between the deterministic and stochastic processes was also due to different ecological responses of special taxa (e.g., habitat specialists) to the environmental changes under a specific situation (Leibold et al., 2017; Gad et al., 2020). For example, the abundance of Cyphellophora with a significant relationship with the PGPR species decreased after the inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 in black soil (Supplementary Figure S6). Previous research showed that several Cyphellophora species have been associated with potential pathogens (Decock et al., 2003). Consequently, the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 could form a protective microdomain in the soybean rhizosphere to defense against soil-borne fungal pathogens, especially Fusarium (wilt disease) (Fu et al., 2017; Mendes et al., 2018).

The VPA (Supplementary Table S7) and Mantel test (Supplementary Tables S2, S3) all demonstrated that nitrogen content was a vital environment factor for altering the microbial community, and by altering HN and AN contents, the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 could change the microbial community and further drive the assembly process in the black and cinnamon soil, respectively. This consistency may be attributed to the increase in the abundance of Azospirillales with the co-inoculation of Bradyrhizobium japonicum 5038 and Bacillus aryabhattai MB35-5 (Supplementary Table S6). According to the above analysis, although the change in the microbial diversity is not significant, the co-inoculation of Rhizobium and Bacillus can substantially influence the interaction and community construction of soil microorganisms to a great extent via changes in the nitrogen content.