The experiment was carried out in an agricultural area belonging to the Teaching, Research and Extension Farm (FEPE) of UNESP—Campus of Jaboticabal, starting in November 2019 and ending in March 2020. The climate of the Jaboticabal region is Cwa, according to the Köppen classification, with wet summers and dry winters. The UNESP unit is located in the northwestern region of São Paulo state (21° 15′ 22″ S, 48° 18′ 58″ W, and an altitude of 605 m).

Plots in the field condition of 30 m² (6 m × 5 m) were established, with a spacing of 0.5 m between rows and 0.1 m between plants. The soil of the area was of the latosol type, and sowing was carried out mechanically at a depth of 5 cm. Chemical analysis of the soil was previously carried out to evaluate its acidity and nutrient availability. The nitrogen dose, urea corresponding to the control treatment was applied manually by repetitions, row by row, after the emergence of the plants. Due to the history of having successive crops with soybean inoculated with strains of Bradyrhizobium, there was no need for reinoculation of this inoculant in this area.

The microorganisms used in this study were sourced from diverse sources. Specifically, Bacillus subtilis was obtained from soil and identified through sequencing and was deposited in the bank under the accession number MZ133755. Bacillus aryabhattai CCT4005 and Streptomyces sp. 244 were provided by André Tosello, and Saccharopolyspora spinosa ATCC 49460 was obtained from the American Type Culture Collection (ATCC). The bacterial strains belong to Laboratory of Soil Microbiology (LSM) of FCAV/UNESP in Jaboticabal. Microorganisms were inoculated in suitable culture media for growth. The inoculants B. subtilis and B. aryabhattai were prepared in nutrient broth, Streptomyces sp. in SM medium (1% glucose, 3% casein peptone, 1.0% yeast extract, 0.05% potassium dihydrogen phosphate, and 0.2% magnesium sulfate heptahydrate) and S. spinosa in Jensen’s medium (2% sucrose, 0.1% dipotassium phosphate, 0.05% magnesium sulfate, 0.05% sodium chloride, 0.01% ferrous sulfate, 0.0005% sodium molybdate, 0.2% calcium carbonate and 1.5% agar). All the strains were grown in a BOD chamber for 72 h, after which a final concentration of 1 × 10⁸ colony forming units (CFU) mL⁻¹ was applied. In total, two applications were carried out, the first 15 days after planting (DAP) and the second 30 days after planting. The applications of the inoculants were carried out using 1000 mL ha⁻¹ backpack sprayers via foliar way.

After the physiological maturity of the grains, which were obtained at 120 days after sowing, the soybean was harvested by means of a specialized parcel harvester, and the grains were then packed in paper bags and stored at room temperature until drying.

The dry mass of the aerial part was determined at the beginning of the flowering of the soybean crop. Five randomly chosen plants were collected from each plot during this period, placed in paper bags and dried in an oven with forced air circulation at 65°C for 72 h. After this period, the dry matter was obtained by weighing the material on a semianalytical balance (Lobo et al., 2019).

Five hundred micrograms of dried and ground plant samples were weighed and placed into 50 mL digestion tubes, which were left to decouple at room temperature for 1.5 h. The tubes were then positioned in a digestion block, initially heated to 80°C for 20 min before the temperature was increased to 160°C. The tubes were monitored and removed once the material ascended the tube walls, and most of the HNO₃ evaporated, leaving a clear solution. After cooling, 1.3 mL the concentrated HClO₄ was added to each tube. The tubes were returned to the block, the temperature was increased to 210°C, and digestion was deemed complete when the solution turned colorless with dense white vapor of HClO₄ and H₂O formed above the dissolved material. The tubes were cooled and the contents were diluted to 25 mL with water in a snap-cap glass (Wan et al., 2020).

For phosphorus analysis, 1 mL of the digested sample was transferred to a test tube, to which 4 mL of water and 2 mL of reagent mix (comprising equal parts of 5% ammonium molybdate and 0.25% vanadate) were added. The mixture was allowed to rest for 15 min before measuring the absorbance at 420 nm using a UV VIS spectrophotometer (Meier et al., 2020).

Using a box with holes adapted for the simultaneous counting of 100 grains of the crop, the mass of 1000 grains was determined by weighing eight subsamples with 100 grains each on a semianalytical scale.

At the end of the field experiment, the soybeans were harvested by a specialized plot harvester, stored at room temperature for drying and later weighed on a semianalytical scale to determine the total grain yield. The productivity in kg ha⁻¹ relative to the production per unit of area was also determined.

Data normality and heteroscedasticity were assessed by the Shapiro–Wilk and Bartlett tests, respectively. Analysis of variance (ANOVA) was performed for each experiment to detect differences, and Tukey’s multiple comparison test was used to compare the means at the 5% significance level. PCA was performed after data standardization (to avoid the influence of different units for the means of response variables), and all the statistical analyses were performed with R software (R Core Team, 2023).

Sterilized roots were macerated using a sterile mortar and pestle with liquid nitrogen. The PowerMax Soil DNA Extraction Kit (Mo Bio Laboratories, Carlsbad, CA) was used to extract the genomic DNA according to the Manufacturer’s instructions. The concentration of the extracted DNA was determined by fluorometry (Qubit^™ 3.0, Invitrogen), and purity was estimated by calculating the A260/A280 ratio via spectrophotometry (NanoDrop^™ 1000, Thermo Fisher Scientific) (de Souza et al., 2016). The hypervariable region V4 of the 16S rRNA gene was amplified using the primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). Three forward primers were used for amplification. These were modified by adding degenerate nucleotides (Ns) to the 5′ region to increase the diversity of target sequences (de Souza et al., 2016). PCR was performed in 30 cycles using the HotStarTaq Plus Master Mix Kit (Qiagen) under the following conditions: 94°C for 3 min, followed by 28 cycles at 94°C for 30 s, 53°C for 40 s, and 72°C for 1 min, and a final elongation step at 72°C for 5 min. PNA clamp sequences (PNA Bio) were added to block amplification of the 16S rRNA gene from the ribosomes and mitochondria. The amplification products were analyzed on a 2% agarose gel to determine the amplification success and relative intensity of the bands. Amplicon sequencing was performed on the Illumina MiSeq platform (Lundberg et al., 2013).

The initial data quality was assessed using the “FastQC” program (Version: 0.11.9; Andrews, 2010). Then, the parameters for quality processing were determined through the “USEARCH” program (Version: 11.0.667), where the thresholds for quality pruning (“-fastx_info”) and size (“-fastq_eestats2”), as well as the position of the primers used for amplification of the target region (“-search_oligodb”). The detected primers (515F. 5′-GTGCCAGCMGCCGCGGTAA; 806R 5′-GGACTACHVGGGTWTCTAAT) were then removed by the “Atropos” program (Version: 1.1.31; Didion et al., 2017). The low-quality ends were pruned by the “PRINSEQ-lite” program (Version: 0.20.4; Schmieder and Edwards, 2011), removing bases whose average quality window was less than Q20 (read forward) or Q18 (read reverse) (“-trim_qual_window 3”; “-trim_qual_right 20/18”), as well as complete sequences with a mean quality lower than Q20 (“-min_qual_mean 20”). The processed reads were submitted to the “DADA2” pipeline (Version: 1.22.0; Callahan et al., 2016) through its package for the statistical program “R” (Version: 4.1.2; R Core Team, 2020). The reads were subjected to a quality control process (“filterAndTrim”), where size truncation (“truncLen = 250”) and quality filtering (“maxEE = 2”) were performed based on the values stipulated by the “USEARCH.” Then, the exact amplicon sequence variants (ASVs) were designated for each sample, which was later filtered to retain possible chimeric sequences (“removeBimeraDenovo”) and then taxonomically classified (“assignTaxonomy”) based on the NCBI RefSeq 16S rRNA database supplemented by RDP (Version: 16; Cole et al., 2014). The taxonomic annotations of the ASVs, as well as their per-sample counts, were exported into a “phyloseq” object (R package “phyloseq”; Version: 1.38.0; McMurdie and Holmes, 2013) and transformed into compositional data (function “phyloseq_standardize_otu_abundance” from the R package “metagMisc”; Version: 0.0.4; Tedersoo et al., 2022). The datasets generated for this study were deposited in the NCBI: Sequence Read Archive (SRA), under the accession number PRJNA1086858.

The extent of sequencing coverage per sample was inferred through rarefaction curves obtained by the “amp_rarecurve” function of the R package “ampvis2” (version 2.7.17; Andersen et al., 2004). Richness (Chao1) and diversity (Shannon and Gini-Simpson indices) indices were calculated using the “alpha” function of the R package “microbiome” (Version: 1.10.0). The means of these measurements were compared pairwise, using the Student/Wilcoxon t-test (depending on the adherence of the data to normality and homoscedasticity, measured by the Shapiro–Wilk and Bartlett tests, respectively). For beta diversity, the dissimilarities were calculated using the Bray–Curtis index (“distance” function of the R package “phyloseq”), which was used for principal coordinate analysis (PCoA). The statistical significance of the separation of the evaluated treatments was determined through a PERMANOVA, considering a p-value of 0.05. Finally, the taxa whose abundances significantly changed were identified. For this purpose, the “DESeq2” approach (Version: 1.34.0; Love et al., 2014), which models the data assuming a negative binomial distribution and implements the Wald test to compare the means (p-value adjusted <0.05), was used. The taxonomic profiles and other graphical representations were generated in R using the “ggplot2” package (Version 3.3.5). The networks of each treatment were established based on the classified genera of each treatment. A minimum abundance filter was used to avoid spurious correlations (minimum relative abundance of 0.1%). The relationships between genera were inferred through Pearson’s correlation coefficients. For the calculation of correlations, the tables of counts of ASVs aggregated at the sex level were submitted to the function “corr.test” of the “psych” R package (Version: 2.2.5), considering a minimum correlation threshold of 0.75 (strongly positive or negative relationships) at a confidence level of 95% (p-value <0.05). Correlation values were provided to the function “graph_from_data_frame” of the R package “igraph” (Version: 1.3.1; Csardi, 2005) to construct the network graph. Centrality and topological measures were obtained from the graphs. The numbers of co-occurring genera (nodes), relationships between them (edges), and the formation of modules were evaluated. The number of connections (degree), intermediation ability (betweenness), and tendency to form dense clusters (clustering coefficient) were also estimated. The most important genera (Hubs) were determined from the computation of Kleinberg’s hubbiness score (Kleinberg, 2000).

The microbial communities were mostly dominated by the phylum Proteobacteria (Figure 2A). This phylum represented, on average, 93.75% of the relative abundance. The BS treatment had the greatest relative abundance of Proteobacteria in the samples (average: 95.41%). Sequences derived from other phyla of lower abundance were concentrated in Bacteroidetes (average: 2.82%), Actinobacteria (average: 2.01%), and Firmicutes (average: 1.09%). On average, sequences unclassified at the phylum level represented 0.16% of the relative abundance of the samples.

Among the most abundant genera (Figure 2B), there was a high prevalence of members of the Proteobacteria phylum (21 of the top 25; Figure 2B). Together, the genera of these phyla represented, on average, 88.61% of the samples. There were also genera of the phyla Bacteroidetes (2/25), Actinobacteria (1/25), and Firmicutes (1/25). The massive presence of the genus Bradyrhizobium, which represented, on average, 71.22% of the sequences at the genus level, is remarkable. Bradyrhizobium was the most abundant genus in all the samples, except for sample ST (ST treatment), whose abundance was only 12.66%. Four ASVs of the genus Bradyrhizobium were detected; however, none of them were classified at the species level. The second most abundant genus throughout the samples was Salmonella (average: 1.94%), corresponding to a single ASV whose species was not classified. The third most abundant genus, Massilia (average: 1.59%), presented nine different ASVs, six of which were from unclassified species and three of which were assigned to the species M. agri strain K-3-1, M. phosphatilytica strain 12-OD1, and M. kyonggiensis strain TSA1.

The taxonomic profiles revealed that the BA treatment had greater consistency in the presence of different genera. This can be seen in the high average abundance of genera other than Bradyrhizobium (average 38.77%) and in the higher percentage of less abundant genera (“Others”; average 12.67%) (Figure 2B). In the other treatments, there were occasional samples with less dominance of the genus Bradyrhizobium; however, this was not consistent throughout the entire treatment (Figure 2B).

Despite the noticeable differences observed in the proportions of the taxonomic profiles (Figure 2), alpha diversity measures were similar between treatments since there was no significant difference in pairwise comparisons between the different indices evaluated (Table 4). In the measurement of the richness of ASVs, inferred by the Chao1 index, there was an increase in the values of the treatments inoculated with microorganisms (BS, BA, ST, and SS; Table 1); however, the variance of these values in the samples did not allow us to identify statistical differences when compared with the control. Similarly, for the diversity indices (Shannon and Gini-Simpson), the BA and ST treatments had greater values but were not significantly different from the control treatment (Table 4).

The similarity in sample compositions was also reinforced by beta diversity analysis. The Bray–Curtis distances between the samples did not show any separation between the treatments (PERMANOVA: p-value >0.05), although the composition can explain most of the variation found within them (first two axes = 78.64% of the variation; Figure 3).

Although the prokaryotic communities present in the treatments were narrowly similar, some significant differences were found when prospecting for differentially abundant taxa (DATs). In total, 30 DATs were found, 24 of which were unique. The measures of differences (fold changes), p-values, and relative abundances in the treatments. Differences were concentrated on more specific taxonomic ranks, such as family, genus, and species (Figure 4). Compared with the control treatment (CT), the ST treatment caused the most shifts at 14 DAT (Figure 4C). Of these, 12 were found in greater abundance in the ST samples, while two other species were more prevalent in the CT samples.

The SS treatment had the least number of shifts (4 DAT), mainly reductions in the abundance of taxa present in the CT treatment (Figure 4D). The abundances of the Sphingobacterium mucilaginosum strain THG-SQA8 and the Stenotrophomonas maltophilia strain ATCC-13637 seemed to be most affected by inoculation since their abundances decreased after 4 and 3 of the treatments, respectively.

The microbial community structure was explored through co-occurrence networks (Figure 5 and Table 5). The relationships of taxa at the genus level seemed to be affected by the treatments, as structural shifts were observed (Figure 5). The CT network had intermediate centrality measures in general comparison with other networks.

The BS and BA treatments had greater tendencies to form larger modules, with modularity values of 0.524 and 0.687, respectively (Table 5). Additionally, the BA treatment group presented a greater number of modules (7; Table 5). Notably, this treatment had strong negative correlations (Figure 5 and Table 5). The modularity values of the ST and SS treatments were considerably lower than those of the CT network, at 0.074 and 0.081, respectively. The SS had only two modules.

Another aspect noted in the structuring is the density of the modules. The ST treatment presented a module composed of strongly correlated genera (Figure 5), with a high number of edges and a high mean degree of connections per node (Table 5). This resulted in a high clustering coefficient (0.995; Table 5). The same observations apply—to a lesser extent—to the BS treatment. The BA and SS treatments, on the other hand, had opposite trends, as they had a lower number of edges and mean degree by genus. Thus, their clustering coefficients were lower than those of CT (Table 5). Another consequent effect of the greater number of rarefied modules in the BA and SS treatment networks was the greater dependence on intermediate taxa to maintain the internal links of the modules. This can be verified by the highest means and maximums of betweenness of these treatments (Table 5). Finally, changes are observed in the genera considered hubs since the results differed between treatments (except for the BS treatment, whose taxa are present in the densest submodule tied in this measure) (Table 5). However, considering the same exception, all hub taxa belonged to the phylum Proteobacteria (Figure 5 and Table 5).

The data show statistical significance at a significance level of 5% based on the ANOVA. Figure 6 demonstrates that all bacterial treatments exhibited similar or superior performance in productivity and production compared to the control and S. spinosa (SS), which demonstrated the highest mean productivity (4407.32 kg/ha) in contrast to the control (4201.32 kg/ha). For the productivity parameter, the highest value was observed in the SS treatment compared with the control, whereas no statistically significant difference was found between the other treatments and the control. Regarding nitrogen content, only the ST treatment did not exhibit a statistically significant difference from the control, whereas the other treatments demonstrated higher values than the control.