Three distinct soybean varieties, Zhongliandou (Z1510), Hejiao (H0269), and Maiyu (M579) were used in this study. Z1510 and M579 are both genetically modified soybeans. The Z1510 transformed from ZhongHuang 6,106 with the gat and eps(g2) gene, conferring glyphosate tolerance. M579 derived from DBN9004, contains epsps and pat gene, providing glyphosate and glufosinate-ammonium tolerance. In contrast, H0269 is a conventional soybean variety that does not possess these modifications and relies on traditional weed control methods.

Three distinct soybean varieties cultivated at the experimental base of the Xing’an League Agriculture and Animal Husbandry Science, Institute in Inner Mongolia, a region of significant importance for soybean production in China. The site (45°43′-46°18’N, 121°53′-122°52′E) has a temperate continental climate with chernozem and chestnut soils (pH 8.03, organic matter 34.76 g/kg, total N 3.437 g/kg, total K 2.22%, alkaline N 252.54 mg/kg, available P 20.78 mg/kg, available K 88.75 mg/kg). A randomized complete block design with three replications was employed. Planting utilized a double-row ridge system (60 cm ridge spacing, 0.6 cm plant spacing) at 280,000 plants/ha. Z1510 and M579 received 2,250 mL/ha of 41% glyphosate isopropylamine salt, while H0269 received conventional herbicides. Base fertilizer was applied at >450 kg/ha.

Samples were collected at seedling and flowering stages using a five-point sampling method. Rhizosphere soil was obtained from roots at 10–30 cm depth, sieved (2 mm), and stored at −80°C. DNA was extracted using HiPure Universal DNA Kit. The V3 and V4 region of 16S rRNA and ITS1 region were amplified and sequenced on a Qsep-400 platform (Persoon et al., 2017). Sequences were processed using QIIME 2 (version 2021.4) and DADA2 to obtain amplicon sequence variants (ASVs; Kesim et al., 2023). Taxonomic annotation used the SILVA 138 database (Bars-Cortina et al., 2023). Alpha diversity was assessed in QIIME 2. Beta diversity was evaluated using principal coordinates analysis (PCoA) based on Bray-Curtis distances. Linear discriminant analysis effect size (LEfSe) identified differentially abundant taxa (LDA > 3.5). Functional classification used FUNGuild (version 1.0; Djemiel et al., 2017).

At maturity on September 28, yield components were measured following national standards. Plants per m² were counted in 2 m² areas, while single plant grain weight (g) was calculated by dividing total grain weight by the number of plants in the same 2 m² area. The 100-grain weight (g) was determined as an average of 3 replicates, and moisture content (%) was measured using a moisture meter with 3 replicates. Finally, the theoretical yield (kg/ha) was calculated using the formula: plants/m² × 10,000 × single plant grain weight× 10⁻⁵× 0.9. This comprehensive approach provides a standardized method for estimating crop yields based on key components measured in sample areas, with multiple replicates used for some measurements to enhance accuracy and reliability.

The 16S rRNA gene and ITS amplicon sequencing data were processed using the SILVA 138/16S rRNA database and Unite 8.0 database for taxonomic assignment. Quality control measures were implemented to ensure data integrity, including the detection and removal of chimeric sequences using the UCHIME (v1.9.1) software. Alpha diversity indices were calculated using the MOTHUR (v1.30.2) software to assess the richness and diversity of microbial communities across samples. To investigate the overall community composition and assess the similarities or dissimilarities between samples, beta diversity analysis was performed using principal coordinate analysis (PCoA) implemented in the QIIME (v1.9.1) software package (Caporaso et al., 2010). The resulting PCoA plots were generated using the R packages “vegan,” “ecodist,” and “ggplot2” in R (v4.1.3), enabling the visual representation and interpretation of beta diversity patterns (Schloss et al., 2009; Cole et al., 2014; Chen and Boutros, 2011). Correlations between microbial community composition (top 20 dominant bacterial taxa) and yield traits were examined using Pearson’s correlation coefficient. Data were subjected to a one-way analysis of variance with factors of treatments and expressed as means +SD. Comparisons between any two groups were performed by unpaired Student’s t-tests.

The 16S rDNA and ITS sequencing analyses revealed significant microbial diversity in rhizosphere soil samples. Using the 341F-806R primer pair, 404,470 denoised and non-chimeric ASVs were obtained from rhizosphere soil samples, averaging 67,411 ASVs per replicate (Supplementary Table S1). The ITSIF-ITSIR primer pair yielded 807,404 ASVs from surrounding soil samples, with an average of 134,567 ASVs per replicate (Supplementary Table S2). The higher ASV count in rhizosphere soil samples suggested potentially greater microbial diversity.

The alpha diversity of rhizosphere microbial communities in M579, Z1510, and H0269 soybean varieties revealed distinct patterns for bacterial and fungal populations (Table 1). At both seeding and flowering stages, bacterial abundance consistently followed the order M579 > Z1510 > H0269, with M579 exhibiting significantly higher Chao1 and ACE indices compared to H0269. M579 maintained 5.8% higher Shannon diversity than H0269 across growth stages. Notably, bacterial diversity indices remained stable within each soybean variety, suggesting that genetic factors have a stronger impact on bacterial community composition than developmental changes. In contrast, fungal communities showed more conserved traits, with 123 shared ASVs across all samples and no significant diversity differences among varieties or growth stages. However, fungal ASV abundance showed diverged between seedling and flowering stages. M579 and Z1510 decreased significantly, while H0269 remained stable. Shannon indices showed no significant differences among varieties. These findings revealed the complex dynamics of soybean rhizosphere microbiomes, showcasing distinct bacterial and fungal responses to plant genotype and developmental stage.

Beta diversity analysis using PCoA revealed significant variations in rhizosphere microbial community composition across soybean growth stages and genotypes. The PCoA1 axis accounted for 18.55 and 33.78% of the total variance in bacterial and fungal community composition, respectively. At seedling stage, Z1510 and H0269 showed similar bacterial and fungal communities, distinct from M579. At flowering stage, bacterial communities of Z1510 and M579 converged, while H0269 diverged (Figure 1A). However, Fungal communities at flowering stage showed less distinct clustering among varieties, but M579 displayed high within-group variability (Figure 1B). These findings underscored the dynamic nature of rhizosphere microbiomes, influenced by both plant genotype and developmental stage.

Analysis of soybean microbiomes across growth stages revealed significant variations in bacterial community structure. At the bacterial phylum level, the dominant groups identified were Proteobacteria, Acidobacteriota, and Actinobacteria. During the seedling stage, Proteobacteria exhibited the highest abundance in the H0269 variety, while M579 showed the lowest levels. Conversely, Acidobacteriota was most abundant in M579, with H0269 and Z1510 displaying similar abundances. Notably, the abundance of Actinobacteria remained relatively consistent across all three soybean varieties. The flowering period saw increased Proteobacteria and Acidobacteriota across all varieties, while Actinobacteria declined (Figure 2A). At the genus level, the dominant genera during both seedling and flowering periods were Sphingomonas, Gemmatimonas, RB41, Vicinamibacteraceae, and WD2101_soil_group (Figure 2B).

Fungal communities were dominated by Ascomycota, Basidiomycota, and Zygomycota phyla. At the seedling stage, Ascomycota was most abundant in H0269, followed by Z1510 and M579, while Basidiomycota peaked in M579. These abundance patterns persisted into the flowering stage, with Ascomycota increased in H0269 and Z1510 but decreased in M579 (Figure 2C). At the genus level, Fusarium and Guehomyces dominated, with Guehomyces showing higher abundance in M579 during the seedling stage. Conversely, Fusarium became more prevalent in H0269 during flowering (Figure 2D).

The LEfSe analysis of rhizosphere bacteria across three soybean varieties revealed significant differences in functional bacteria recruitment between seedling and flowering stages. At seedling stage, a total of 89 bacterial species were recruited, with M579 recruiting the highest number (48 species), followed by Z1510 (21 species), and H0269 (20 species). M579 was enriched with taxa such as P_Acidobacteria, c_Vicinamibacteria, and P_Chloroflexi. H0269 had o-Xanthomonaadales, f_Rhodanobaceae, and d_Rhodanobacte; and Z1510 had o-Gemmatimonadales and related taxa. During flowering stage, the total number of bacterial species decreased to 51, with H0269 recruiting 25 species, and both Z1510 and M579 recruiting 13 species (Figures 3A,B). H0269 showed enrichment in c_Alphasproteobacteria and c_Actinobacteria. Z1510 in o_Burkholderiales and f_Comamonadaceae, and M579 in f_Oxalobacteraceae and related taxa. From seedling to flowering stage, M579 experienced a substantial decrease of 33 bacterial species, Z1510 decreased by 8, while H0269 increased by 5.

Similarly, LEfSe analysis of rhizosphere fungi indicated significant differences in functional fungi recruitment between growth stages. In the seedling stage, a total of 67 fungal species were recruited: M579 recruited 36 species, Z1510 recruited 18 species, and H0269 recruited 13 species. During the flowering stage, this number dropped to 25 fungal species overall, with H0269 recruiting 12 species, M579 recruiting 8 species, and Z1510 recruiting 5 species. From seedling to flowering stages, M579 saw a decrease of 28 fungal species, Z1510 decreased by 13, and H0269 decreased by 6. Fungal taxa with LDA scores greater than 5 included f_Nutriaceae in H0269, o_Hypocrales and s_Clonostachys_rasea_f_catenulata in Z1510, and d_Ratobasidium and s_Ratobasidium-ramicola in M579 (Figures 3C,D).

Yield trait analysis across three soybean varieties revealed a consistent pattern, with Z1510 > M579 > H0269 for theoretical yield and grain number per plant. Z1510 showed superior yield indicators but lower 100-kernel weight and moisture content compared to M579 and H0269 (Table 2), suggesting a compensatory mechanism of increased grain number offsetting reduced grain weight.

The study examined relationships between rhizosphere microbes and five yield indicators: plant density, grains per plant, 100-grain weight, moisture content, and theoretical yield. The results revealed 15 bacterial and 13 fungal species exhibited significant associations with yield components. Among the bacteria, taxa such as s_wastewater metagenome and s_Taibaiella sp. demonstrated positive correlations with plant density but negative associations with grain weight. Conversely, species like S_Pedobacter panaciterrae showed positive correlations with theoretical yield, grain weight, and water content (Figure 4A). Fungal species, including s_Ascomycota sp. FL 2010c, negatively correlated with grain number and yield per plant while positively associating with water content. Conversely, S_Guehomyces pullulans exhibited positive correlations with yield components but negative associations with water content (Figure 4B).