The Influence of Transgenic Insect-Resistance and Herbicide-Tolerance Soybean KM2208-23 on the Rhizosphere Micro-Biome

1. Introduction

The global cultivation area of genetically modified (GM) crops has consistently expanded, with soybean (Glycine max L.) being one of the most widely adopted transgenic species [1]. These modifications primarily target agronomic traits such as insect resistance (e.g., via Bacillus thuringiensis (Bt) Cry and Vip proteins) and herbicide tolerance (e.g., to glyphosate or glufosinate), aiming to enhance crop yield and reduce management costs [2]. Notably, stacked-trait varieties, combining multiple transgenes, are increasingly prevalent to address complex field challenges [3]. However, the potential environmental impacts of these transgenic plants, particularly on non-target soil ecosystems, remain a critical component of comprehensive biosafety assessment. Among these concerns, the effects on the structure and function of rhizosphere microbial communities warrant in-depth investigation, given their indispensable roles in soil fertility, nutrient cycling, plant health, and overall ecosystem stability [4].

The rhizosphere, a narrow zone of soil directly influenced by root exudates and activities, harbors an immensely diverse and dynamic microbial consortium [5]. Its community composition is highly sensitive to alterations in plant physiology, biochemistry, and root architecture [6]. Consequently, the introduction of transgenes may indirectly modify the rhizosphere microbiome through changes in root exudation patterns, litter quality, or subtle shifts in plant growth and development [7]. While numerous studies have examined the impact of single-trait GM crops (especially Bt crops) on soil microorganisms, conclusions have been variable, often reporting minor, transient, or plant genotype-dependent effects rather than consistent adverse outcomes attributed solely to the transgene [8,9]. Recent meta-analyses and high-throughput sequencing studies suggest that factors such as plant growth stage, soil type, and agricultural practices often exert stronger influences on microbial communities than the transgenic trait itself [10]. Nevertheless, research on multi-gene stacked traits, particularly those combining insect resistance and novel herbicide tolerance mechanisms, is comparatively limited. The potential synergistic or novel interactive effects of multiple transgenes on the soil microbiome are not yet fully understood, representing a significant knowledge gap.

Most existing temporal studies have focused on limited growth stages, potentially missing critical shifts during key physiological transitions. A comprehensive, stage-resolved analysis-from early seedling establishment through vegetative growth to reproductive maturity and senescence-is essential to capture the full dynamic response of microbial communities to transgenic plants. Furthermore, simultaneous characterization of both bacterial (via 16S rRNA gene sequencing) and fungal (via Internal Transcribed Spacer (ITS) sequencing) communities provides a more holistic view of the rhizosphere ecosystem, as these groups respond differently to plant cues and play distinct ecological roles.

This study aims to explore the potential impact of this complex trait transgenic soybean on the microbial community in the rhizosphere. We compared the dynamic changes of bacterial (16S rRNA) and fungal (ITS) communities in the rhizosphere of the transgenic lines and their non-transgenic isogenic lines before sowing and at five key growth stages (V3, R3, R5, and R8). We hypothesized that, given the targeting nature of the exogenous genes and the dominant influence of plant development and environmental factors, any differences between the transgenic and non-transgenic soybean rhizosphere microbial communities would be minor, transient, and within the natural variation range caused by plant growth stages. The results of this study will provide important scientific basis for the ecological safety assessment of transgenic crops with complex traits.

2. Materials and Methods

2.1. Genetically Modified Soybean

The transgenic transgenic insect-resistance and herbicide-tolerance soybean KM2208-23, which contained the cry1Ac, vip3Aa19, mOsPPO2 and pat genes, was selected as the experimental treatment material, while the non-genetically modified corn QYZ014 served as the control. The experience was performed at No. 53, Qinglong River Road, Hongshiya Sub-district, Huangdao District, Qingdao City, Shandong Province, at the test base of Qingyuan Compound Co., Ltd (36°04′37″N，120°06′13″E). For each soybean variety, 3 replicate plots were set up. In each plot, the 5-point diagonal sampling method was adopted. Two soybeans were taken from each point. The 10 root samples from these 5 points were combined into one sample and placed in a sampling bag. Two samples were taken from each plot, totaling 6 replicates for each variety in each period. Rhizosphere soil of GMO and CK soybean at five growth stages including BQ (pre-sowing), V3 (trefoil stage), R3 (beginning pod), R5 (be-ginning seed) and R8 (full maturity) were collected. In brief, the whole root system was pulled out from the field and remove the loosely attached soil on the roots by shaking. Rhizosphere soil were finally obtained by rinsing the roots with sterile water and high-speed centrifugation of the suspension (6000 g, 20 min).

2.2. Determination of Physical and Chemical Properties of Rhizosphere Soil

The total kalium (TK), available kalium (AK), effective nitrogen (EN), pH, organic matter (OM), total phosphorus (TP), available phosphorus (AP) and total nitrogen (TN) of the rhizosphere soil were detected. In brief, the total nitrogen in the soil, with the participation of the accelerator, is digested using concentrated sulfuric acid. It is converted into ammonium nitrogen. It is alkalinized with sodium hydroxide, and the ammonia vaporized by heating is absorbed with boric acid. It is titrated with an acid standard solution to determine the total nitrogen content of the soil. Under high-temperature conditions, the phosphorus-containing minerals and organic phosphorus compounds in the soil react with sulfuric acid and perchloric acid, causing them to completely decompose and all converting into orthophosphate and entering the solution. Then, the molybdenum-sulfur antimony colorimetric method is used for determination. Under heating conditions, the soil organic carbon is oxidized using an excess of potassium dichromate-sulfuric acid solution. The excess potassium dichromate is titrated with a standard solution of ferrous sulfate. The amount of organic carbon is calculated based on the consumed amount of potassium dichromate using the oxidation correction coefficient, and then multiplied by the constant 1.724 to obtain the soil organic matter content. In the diffusion dish, the soil is hydrolyzed with 1.8 mol·L^-1 NaOH to convert the readily hydrolyzable nitrogen (potential available nitrogen) into NH₃ through alkaline hydrolysis. After NH₃ diffuses, it is absorbed by H₃BO₃. The NH₃ in the H₃BO₃ absorption solution is then titrated with standard acid, from which the content of alkali-hydrolyzed nitrogen in the soil can be calculated. The acidic forest soil samples were extracted with ammonium fluoride and hydrochloric acid. By taking advantage of the ability of F^- to complex Fe³⁺ and Al³⁺ in acidic solutions, the relatively active iron phosphate and aluminum salts in these soils were gradually activated and released. The phosphorus in the leaching solution was determined by the molybdenum-sulfur antimony colorimetric method. Fluorine decomposes silicate minerals by reacting with silicon to form fluorosilicon. It can be heated and volatilized under the presence of strong acids. Perchloric acid is a very strong oxidizing agent under high-temperature conditions. It can decompose organic matter in the soil and effectively remove excess HF in the samples.

2.3. DNA Extraction and Amplicon Sequencing

Extract the DNA from 0.5g rhizosphere soil using the Fast DNA SPIN Kit for Soil (MP Bio-medicals, Solon, USA). DNA quality testing was performed using Qubit® dsDNA HS Assay and NanoDrop. The bacterial 16S rRNA gene was amplified using primers 341F (5’-CCTAYGGGRBGCASCAG-3’)/805R (‘-GACTACHVGGGTATCTAATCC-3’). The fungal ITS region was amplified using primers ITS1F (5’-CTTGGTCATTTAGAGGAAGTAA-3’)/ITS2R (5’-GCTGCGTTCTTCATCGATGC-3’). The amplification reaction contains high-fidelity enzymes (such as KAPA HiFi HotStart ReadyMix), primers (10 μM each 0.2 μL), and template DNA. The amplification procedure is 95℃ 3 min, denaturation: 95℃ 30 s, annealing: 55℃ 30 s (for bacteria) or 50℃ 30 s (for fungi), extension: 72℃ 30s (30 cycles), final extension: 72 °C 5 min. The library was conducted using the Nextera XT Index Kit and quantify the library concentration using Qubit®. Library fragment sizes were confirmed using a Gilent 2100 Bioanalyzer or LabChip GX. After library construction, sequencing was performed using the Illumina MiSeq (2×250/300bp, suitable for high throughput) platform.

2.4. Bioinformatic Analysis of Amplicon Sequencing Data

The analysis of bacterial 16S rRNA genes and fungal ITS sequences was conducted using QIIME version 1.9 [11] and USEARCH version 10 [12]. Initially, FastQC version 0.11.5 [13] was employed to assess the quality of the reads. Subsequently, Trimmomatic version 0.39 [14] was utilized to trim paired reads with quality scores below Q30. The processed sequences were then clustered, with those exhibiting a similarity greater than 97% being classified as belonging to the same operational taxonomic unit (OTU). Taxonomic classification of the sequences was performed using the SILVA version 138.1 [15] and UNITE version 8.2 databases [16] to distinguish between bacterial and fungal sequences. For α diversity estimation, the OTU tables for bacteria and fungi were standardized using the normalize_table.py script in QIIME, following the cumulative sum scaling (CSS) method. β diversity analysis was conducted using the beta_diversity.py script, also employing the CSS normalization approach. Operational taxonomic units present in all samples were identified as core taxonomic units. Principal coordinate analysis (PCoA) was performed using the R software environment, version 4.1.0 [17], employing the Vegan package version 2.6.4 [18] and the Tidyverse package version 2.0.0 [19]. PICRUSt2 version 2.6.2 [20] was used to perform functional prediction of 16S rRNA gene sequences in the Kyoto Encyclopedia of Genes and Genomes (KEGG) functional database. Pairwise T-test was performed for different groups, with the P threshold set at 0.05 (P < 0.05 indicate significance). FunGuild version 1.0 [21] was used to perform the ecological function prediction of fungal community. The biomarker bacteria were identified using LEfSe (logarithmic LDA score > 3.0, P < 0.05).

3. Results

3.1. The Influence of GMO on Physical and Chemical Properties of Soil

To investigate the influence of GMO on soil, rhizosphere soil of soybeans from five different growth stages were collected to detect the physical and chemical properties. The five periods were BQ (before sowing), V3 (trefoil stage), R3 (beginning pod), R5 (beginning seed) and R8 (full maturity). Results showed that the EN, OM and TN indexes showed the significant differences between GMO and CK groups at V3 and R3 period (Table 1, Wilcoxon rank sum test, P < 0.05). The pH index showed the between-group variance at the later stage of soybean growth (R5 and R8). The TK and AP indexes showed the differences at BQ and R3 stage, respectively. In general, genetically modified soybean has a greater impact on the physical and chemical properties of the rhizosphere soil at V3 and R3 stage.

3.2. The Basic Information of Amplicon Sequencing Data

At the BQ, V3, R3, R5 and R8 stage, rhizosphere soil of soybean from both the GMO and CK groups were selected for further analysis. We characterized the bacterial and fungal community compositions by sequencing the small ribosomal subunit (16S rRNA) gene fragments and internally transcribed spacer (ITS) sequences, followed by clustering into operational taxonomic units (OTUs) at 97% identity. From 60 samples, we obtained 4,792,307 high-quality 16S rRNA reads and 5,574,493 ITS reads, with individual sample reads ranging from 42,015 to 128,946 and an average of 86,390 reads per sample, representing 3,284 bacterial and 821 fungal OTUs. The dilution curve, constructed based on the Shannon index, demonstrated a tendency to plateau, indicating that species richness in this environment does not significantly increase with additional sequencing, thereby satisfying the requirements for subsequent analysis (Fig. A1).

3.3. The Influence of GMO on Microbial Alpha Diversity of Soybean

Comparative analyses of Shannon Diversity Index (SDI) were conducted between GMO and CK groups for both bacterial and fungal communities. The SDI for bacterial communities did not exhibit statistically significant differences between GMO and CK groups across the five stages, as determined by the Wilcoxon rank sum test (P > 0.05; Figure 1A). In contrast, the SDI for fungal communities was observed to be higher at R5 stage in the GMO groups compared to the CK group (P < 0.05; Figure 1B). In addition, the influence of growth stage on microbial alpha diversity were also investigated. As Figure 2A and 2B shown, no significant difference of the bacterial SDI was found across the five stages for GMO or CK groups (Tukey’s multiple comparisons test P > 0.05). On the contrary, the SDI of fungal community significantly increased from V3 to R3 stages in both GMO and CK groups (Tukey’s multiple comparisons test P < 0.05, Figure 2C and 2D).

3.4. The Influence of GMO on Microbial Beta Diversity of Soybean

Principal coordinates analysis (PCoA) utilizing Bray–Curtis dissimilarity didn’t distinguish GEO and CK groups both in bacterial communities (Figure A2A, Adonis R²=0.022, P=0.097) and fungal communities (Figure A2B, Adonis R²=0.019, P=0.24) without considering the growth stage. Examining separately the different stage, significant compositional differences in bacterial communities were observed between GMO and CK groups at the V3 (Figure 3A, Adonis R²=0.179, P=0.006) and R3 (Adonis R²=0.136, P=0.004) stages along the first principal component. As for fungal community, significant compositional differences were found at V3 (Figure 3B, Adonis R²=0.537, P=0.002), R5 (Adonis R²=0.149, P=0.006) and R8 (Adonis R²=0.138, P=0.043) stages. These results of PCoA indicated that GMO influenced the beta diversity of both bacterial and fungal community at V3 stage.

To investigate whether growth stage affects the composition of microbiome communities, PCoA was conducted across the five stages. For bacteria, samples in each growth stages significantly separated with others for both GMO and CK groups (Figure 4A and 4B). Fungal communities showed a similar pattern in PCoA results. All these PCoA findings suggest that growth stages deeply influence the composition of microbiome communities, and GMO influence these at special growth stages (Figure 4C and 4D).

3.5. The Influence of GMO on Microbial Community Structure of Soybean

To characterize the core microbiota associated with the GMO group, we conducted a comprehensive analysis of the taxonomic compositions of their bacterial and fungal communities. Analysis of microbial communities at the phylum level revealed similar compositions between the GMO and CK groups (Figure A3 and Figure A4). The predominant bacterial phyla (Figure A3A) identified were Proteobacteria (29.81%), Acidobacteriota (18.33%), Actinomycetota (18.32%), and Gemmatimonadota (11.78%). A parallel analysis of fungal communities (Fig. A4A) indicated that Ascomycota (68.97%) and Fungi_phy_Incertae_sedis (6.30%) were the dominant phyla. At the genus level, the principal bacterial taxa included Sphingomonas (5.43%), MND1 (2.84%), Streptomyces (0.49%), and Pseudomonas (0.18%) (Fig. A3B). The dominant fungal genera were Acrophialophora (6.05%), Schizothecium (5.47%), Alternaria (1.45%), and Albifimbria (1.21%) (Fig. A4B). According to the Wilcoxon rank sum test, five bacterial dominant phyla, including Pseudomonadota, Acidobacteriota, Thermoproteota, Chloroflexota and Myxococcota showed significant differences between the GMO and CK groups at V3 stage (P < 0.05). Similar results were also found in fungal community, in which Mortierellomycota, Mucoromycota, Kickxellomycota and Fungi_phy_Incertae_sedis showed the differences between the two groups at V3 stages. It’s worth noting that more fungal genus was found sensitive to GMO than bacteria at V3 stage, that Acrophialophora, Schizothecium, Ascomycota_gen_Incertae_sedis, Chytridiomycota_gen_Incertae_sedis and Alternaria were influenced by GMO. This differential analysis of the bacterial and fungal communities suggested that soybean at V3 growth stage were more sensitive to the GMO treatment.

3.6. The Influence of GMO on Microbial Biomarkers of Soybean

Line Discriminant Analysis (LDA) Effect Size (LEfSe) is an analytical technique that integrates the non-parametric Kruskal-Wallis and Wilcoxon rank sum tests with the effect size derived from LDA. This method is employed to identify biomarkers exhibiting statistically significant differences across various groups. Considering that the bacterial and fungal community of soybean at V3 stage were sensitive to GMO, LDA was performed to find the biomarkers related to GMO at V3 stage. As Figure 5A shown, Lysobacter, Nitrospira, Novosphingobium and Croceibacterium could serve as the bacterial genus-level biomarkers for GMO group. As for fungal community, Ascobolus, Plectosphaerella, Albifimbria, Lectera, Parazalerion, Verticillium and Bisifusarium were the biomarkers for GMO group (Figure 5B).

3.7. The Influence of GMO on Microbial Function of Soybean

To explore the functional disparities among microbial groups, we employed PICRUSt2 to predict 16S rRNA gene sequences within the KEGG database. Subsequently, we conducted pairwise significance tests between samples using the Wilcoxon rank sum test. Given that rhizosphere microbiome of soybean at V3 stage, our analysis concentrated on the functional differences between the GMO and CK groups at this juncture. As illustrated in Figure 6A, pathways at level 2 related to metabolism of cofactors and vitamins, energy metabolism, cell motility, membrane transport, signal transduction, development and regeneration, transport and catabolism and energy metabolism were enriched in the GMO group (Figure 6A). Conversely, pathways at level 3 associated with one carbon pool by folate, glutathione metabolism, bacterial secretion system, biotin metabolism, flagellar assembly, pyruvate metabolism, butanoate metabolism, base excision repair and folate biosynthesis were more prevalent in the GMO group (Figure 6B). In addition, the role that fungi play in the ecosystem was predicted based on fungal taxa using FunGuild. Results showed that, at V3 stage of soybean, endophyte-plant pathogen and wood saprotroph were enriched in GMO group (Figure A5).

4. Discussion

This study comprehensively evaluated the temporal impacts of a novel stacked-trait (insect-resistant and herbicide-tolerant) transgenic soybean on the physicochemical properties and microbial communities of the rhizosphere soil. Our multi-stage analysis revealed a nuanced pattern: while the transgenic soybean did induce detectable shifts in specific soil properties and microbiome composition, these effects were largely confined to specific growth stages, most prominently the V3 stages. Furthermore, the influence of the plant’s developmental stage itself exerted a dominant and consistent effect on microbial community assembly, a force that far outweighed the impact of the transgene over the entire growth cycle. These findings align with the growing consensus that the effects of GM crops on soil ecosystems are often subtle, transient, and context-dependent rather than profound and persistent [22,23].

4.1. Stage-Specific Modulation of Rhizosphere Environment and Microbiome

The most pronounced differences between GM and non-GM soybean rhizospheres occurred during the V3 and R3 stages. Significant alterations in soil available nitrogen (EN), organic matter (OM), and total nitrogen (TN) at these stages suggest a potential shift in root exudation patterns or nitrogen uptake efficiency mediated by the transgenic traits. Root exudates are primary drivers of rhizosphere chemistry and microbial recruitment [24]. Modifications in plant physiology due to the expression of cry1Ac, vip3Aa19, mOsPPO2, and pat genes may lead to qualitative or quantitative changes in exudate profiles, particularly during active vegetative growth and reproductive onset, thereby altering nutrient dynamics. This is supported by the concomitant significant shifts in bacterial and fungal community beta-diversity specifically at V3. The identification of specific bacterial (e.g., Lysobacter, Nitrospira) and fungal (e.g., Albifimbria, Verticillium) genera as biomarkers for the GM rhizosphere at V3 further underscores this critical stage. The sensitivity of the V3 stage microbiome is notable; it may represent a key window where the plant microbiome is being established and is most receptive to plant-mediated changes.

The lack of persistent divergence in alpha-diversity (Shannon index) for bacteria and its only sporadic increase for fungi at R5 indicates that the transgene did not fundamentally reduce or enhance microbial species richness in a lasting manner. This is a crucial finding for ecological risk assessment, as sustained loss of diversity could impair ecosystem resilience. The overwhelming separation of samples by growth stage in PCoA for both treatments powerfully demonstrates that plant phenology is the principal factor structuring the rhizosphere microbiome, consistent with studies on conventional and other GM crops [25]. The transgene effect, while statistically significant at certain points, appears as a secondary modulation within this dominant developmental trajectory.

4.2. Functional Implications of Microbial Community Shifts

The predicted profiling (PICRUSt2) at the sensitive V3 stage revealed intriguing functional differences. The enrichment of pathways related to energy metabolism (e.g., pyruvate metabolism, oxidative phosphorylation), membrane transport, and bacterial motility (flagellar assembly) in the GM rhizosphere suggests a potentially more metabolically active or copiotrophic bacterial consortium. This could be a response to altered carbon substrate availability from GM root exudates. Notably, the enrichment of “bacterial secretion system” pathways might indicate shifts in inter-microbial or plant-microbe interactions, possibly linked to competition or symbiosis [26]. For fungi, the FunGuild prediction indicated a higher relative abundance of endophyte-plant pathogen guilds in the GM rhizosphere at V3. It does not necessarily imply increased disease risk but may reflect a shift in the competitive balance within the fungal community [27]. Some endophytic plant pathogenic fungi are highly responsive to root exudate chemistry [28]. The observed change could be a transient ecological rearrangement rather than a proliferation of pathogens. Long-term field studies are needed to assess the agronomic significance of such compositional shifts.

4.3. Ecological Significance and Biosafety Perspective

From an ecological safety standpoint, our results support the hypothesis that the impact of this stacked-trait transgenic soybean on the rhizosphere microbiome is limited and stage-specific. The effects were most detectable during peak vegetative/reproductive activity when root exudation is high, and they diminished or became non-detectable at later stages (R5, R8). This “transient perturbation” model is consistent with numerous studies on single-trait Bt crops, where differences often converge with those of conventional crops at maturity [29]. The fact that the core microbial phyla (e.g., Proteobacteria, Acidobacteriota, Ascomycota) remained dominant and unchanged in both treatments reinforces the stability of the overall community structure.

The significant changes in key soil nutrients (EN, TN) at V3 stage, coupled with the microbial community shifts, suggest that the transgenic plant interacts differently with the soil nutrient cycle during these phases. However, these biochemical differences also did not persist, indicating a resilient soil system that buffers short-term changes. The ecological relevance of these temporary shifts is likely minimal, as they fall within the range of variation caused by standard agricultural practices, plant genotype, and environmental fluctuations [30].

5. Conclusions

In conclusion, our temporal study demonstrates that the cultivation of the investigated stacked-trait transgenic soybean induces significant but transient and stage-specific alterations in the rhizosphere environment, primarily during the V3 growth stages. These alterations manifest as changes in soil nitrogen parameters, shifts in microbial community composition (particularly for fungi), and modulations in predicted microbial functional potential. Crucially, the overarching influence of soybean developmental stage remains the primary driver of microbial community succession, and no sustained negative impact on microbial diversity was observed. The effects documented here appear to represent minor ecological perturbations rather than fundamental, lasting dysbiosis. These findings con-tribute to a more nuanced understanding of the interactions between complex transgenic plants and soil ecosystems, supporting the view that their environmental impact is con-text-dependent and should be evaluated across the entire crop life cycle.