Soil samples were collected from the artificial infection background of ScALDH21 transgenic and non-transgenic cotton at the Manasi experiment station of the Xinjiang Academy of Agriculture Science. The Latin square design divided the cropping area into 9 plots where transgenic cotton lines 92, 96 (described as L1, L2, respectively) and non-transgenic cotton line 779 (described as NT) were planted, and each line had 9 repetitions. Since 2013, the same cotton lines have been continuously planted in the experimental plots used in this study. Bulk soil samples from the 5 centimeter periphery and 5-10 centimeter depth of cotton plants were collected after a growth season from the center of plots and stored at -80°C.

The soil chemical properties were detected in the Central Laboratory of Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences, based on the methods of forest soil nitrogen (LY/T 1228-2015), phosphorus in forest soil (LY/T1232-2015), potassium in forest soil (LY/T1234-2015), soil organic matter (NY/T1121.6-2006 Soil Testing, part 6), and soil pH (NY/T1121.6-2006 Soil testing, part 2).

Total soil genomic DNA was extracted with the FastDNA SPIN Kit for Soil (MP Biomedicals, Inc., CA). DNA concentration and purity level were monitored on a 1% agarose gel and a Nanodrop 2000 spectrophotometer, respectively. Based on the concentration, DNA was diluted to 1 ng µL^-1 using sterile water. Specific primer pairs (515F: 5’-GTGCCAGCMGCCGCGGTAA-3’ and 806R: 5’-GGACTACHVGGGTWTCTAAT-3’) for the V4 region of 16S rRNA were used to identify bacterial diversity, and ITS (ITS5-1737F: 5’-GGAAGTAAAAGTCGTAACAAGG-3’and ITS2-2043R: 5’-GCTGCGTTCTTCATCGATGC-3’) was used to evaluate fungus diversity.

All PCR reactions were carried out in 30 µL reactions with 15 µL of Phusion® High-Fidelity PCR Master Mix (New England Biolabs). Thermal cycling consisted of initial denaturation at 98°C for 1 min, followed by 30 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s, and elongation at 72°C for 30 s. Finally, elongation was at 72°C for 7 min. The PCR products were purified with the GeneJETTM Gel Extraction Kit (Thermo Fisher Scientific). Sequencing libraries were generated using the Ion Plus Fragment Library Kit (Thermo Fisher Scientific) following the manufacturer’s instructions. The library’s quality was assessed on the Qubit@2.0 fluorometer (Thermo Fisher Scientific). Finally, the library was sequenced on an Ion S5TM XL platform, and 400/600 bp single-end reads were generated (Novogene company).

During analysis the barcodes and primers were removed from the raw sequencing data. After removing the barcode and primer sequences, the sample reads were merged using the FLASH software (V1.2.11, http://ccb.jhu.edu/software/FLASH/) (Magoc and Salzberg, 2011) to obtain Raw Tags. Then, quality control was performed on the Raw Tags using the fastp software to obtain high-quality Clean Tags. Finally, the Clean Tags with a database using the Usearch software to detect chimeric sequences and remove them (Haas et al., 2011), resulting in the final set of effective data known as Effective Tags.

The Effective Tags obtained above, was denoised with the DADA2 module using QIIME2 software (Li et al., 2020), and sequences with and abundance below 5 were filtered to obtain the final Amplicon Sequence Variants (ASVs) and feature table. Next, the classify-sklearn module in QIIME2 software was used align the obtained ASVs with a database in order to obtain species information for each ASV (Callahan et al., 2016). For the 16S regions, the Silva138.1 database (https://www.arb-silva.de/) was used and for the ITS region, used the UNITEv8.2 database (https://unite.ut.ee/).

Principal coordinate analysis (PCoA) and principal component analysis (PCA) analyses were performed using the vegan package in the R programming language (McMurdie and Holmes, 2013; Bro and Smilde, 2014). The vegan package provides specific functions, cmdscale and rda, that can be used to calculate PCoA and PCA. Non-metric multidimensional scaling (NMDS) analysis can be performed using the metaMDS function in the vegan package, which is capable of generating NMDS analysis based on the Bray-Curtis distance.

Alpha diversity was assessed using five indices, including Chao1, Shannon, Simpson, ACE, and PD whole tree. QIIME2 was utilized for index calculation, while R software (version 2.15.3) was employed for dilution curve visualization and analysis of differences among groups based on the alpha diversity index. We utilized the Tukey HSD and Duncan methods of ANOVA to analyze the significance of differences between groups of the alpha diversity index.

Beta diversity was evaluated to determine differences in species complexity among samples. QIIME2 was employed to calculate the Unifrac distance and to construct a UPGMA sample clustering tree. Using R software (version 2.15.3). AMOVA, based on weighted Unifrac distance, was used to test for significant differences between the groups.

We conducted redundancy analysis to demonstrate the relationship between soil chemical properties and OTUs. The relative abundance data of fungal and bacterial communities at the OTU level were considered species data, while soil chemical data were viewed as soil environmental variables. The CANOCO5 software (Smilauer and Jan, 2014) was used to perform RDA analysis to evaluate the correlation between soil microbial communities and soil chemical factors. To assess differences in soil chemical data, ANOVA’s Duncan method was employed for analysis of variance and multiple comparisons. The CANOCO5 software was utilized to display the impact of environmental factors on microbial community structure. Given the nonlinear relationship between variables, we computed correlation coefficients between the two variables using the Spearman correlation analysis method, employing the psych package in the R programming language.

The FunGuild software was utilized to predict fungal functional profiles based on ITS sequencing data using the default database and settings. FAPROTAX (Louca et al., 2016) was used to predict bacterial functional profiles based on 16S rRNA sequencing data using the default database and settings.

To determine chemical properties, soil samples were collected from 5 centimeters around the periphery and 5-10 centimeters deep. We measured eight soil chemical properties in plots planted with two transgenic lines (L1, L2) and non-transgenic cotton (NT). The analysis revealed that the levels of available nitrogen (N) ( Figure 1A ), available phosphorus (P) ( Figure 1B ), and available kalium (K) ( Figure 1C ) in the independent transformed line L2 were slightly higher compared to the transgenic line L1 and the non-transformed line NT ( Figures 1A–C ). The available N ( Figure 1A ) and available P ( Figure 1B ) content in each cotton line did not exhibit significant changes, while a significant difference was observed for available K ( Figure 1C ) between L2 and NT. The levels of soil organic matter in L1 (16.29 g/kg) and L2 (15.67 g/kg) were slightly higher than NT, although this was not significant ( Figure 1D ). Likewise, no significant difference was observed in total N and total K in the soil among the different cotton lines ( Figures 1E, G ). However, it should be noted that the total P content in L1 and L2 was significantly higher than that of NT ( Figure 1F ). The soil pH values of all cotton lines were found to be alkaline, exceeding 8, and there was no significant difference among the three groups ( Figure 1H ).

After filtering out low-quality and short sequence reads, a total of 81073.85 and 83378.48 raw reads were obtained from fungal and bacterial sequencing, respectively. Next, after the initial quality control process, on average 80165.93 and 78516.26 clean reads were obtained, respectively. The average length of fungal ITS was 227.3 bp and that of bacterial 16S rRNA was 417 bp. The sequencing quality was high, with a sequencing error rate less than 1% in clean reads. The operational transcriptional unit (OTU) clustering was performed with 97% consistency, and fungal sequencing yielded 3566 OTUs ( Table S1 ), and bacterial sequencing yielded 7481 OTUs ( Table S2 ). The rationality of total OTUs of each soil sample was detected by rarefaction curves. With the increase in sequencing depth, OTU numbers of fungi ( Figure 2A ) slowly increased, and the curve tended to be smooth. The rarefaction curves of L1, L2, and NT samples were consistent with this trend, and the sequencing depth of fungi was reasonable. The bacterial ( Figure 2B ) rarefaction curve was basically flat but still not saturated, and a small number of bacterial species may not be detected. The number of OTUs increases with the sample size, as depicted in the species accumulation boxplot ( Figures 2C, D ). The curve reaches stability as the number of OTUs saturates and fewer new OTUs are added to each sample, signifying the sequencing depth’s accuracy in representing fungal and bacterial composition. We used rank abundance curves to compare the differences in species composition between transgenic and non-transgenic cotton microbiomes. The results showed that, in different groups, the curve shapes of fungal communities were roughly similar but deviated at the same abundance threshold ( Figure 2E ), while the curve shapes of bacterial communities were also similar but did not deviate at the same abundance threshold ( Figure 2F ). In addition, the inflection points of the curves were analyzed, and they all appeared to have similar abundance values. The analysis of rank abundance curves indicated that the species composition and distribution patterns of root-associated microbiomes in transgenic and non-transgenic cotton were consistent.

The principal coordinate analysis (PCoA) was executed to explore the dissimilarity among microbial communities obtained from different soils. Variation in fungi was explained by 20.43% and 8.89% by PCoA1 and PCoA2, respectively ( Figure 3A ), while variation in bacteria was explained by 15.63% and 9.47% by PCoA1 and PCoA2, respectively ( Figure 3B ). PCoA analysis demonstrated the sample clustering of L1, L2, and NT, where soil samples from different groups clustered together, indicating a similarity in the microbial community structure between the different groups. We carried out principal component analysis (PCA) to evaluate the structural variances of fungal and bacterial communities in the different samples. PC1 and PC2 explained 9.74% and 5.6% of the total variation in fungi, respectively ( Figure 3C ), while PC1 and PC2 explained 5.6% and 4.65% of the total variation in bacteria, respectively ( Figure 3D ). PCA1 and PCA2 scatter plots showed no significant separation between the different groups, suggesting that no significant variation in fungal and bacterial community composition was observed between the samples. Non-metric multi-dimensional scaling (NMDS) was used to analyze the variability between fungal and bacterial communities at different ecological niches. MDS1 and MDS2 scatter plots indicated no significant variances in the fungal and bacterial communities of different groups and showed some similarity in their composition and relative abundance ( Figures 3E, F ).

In this study, we examined the alpha and beta diversity indices for fungi and bacteria in soil samples from transgenic (L1 and L2) and non-transgenic (NT) cotton. These indices reflect the abundance and diversity of fungi and bacteria in the soil. The coverage in all three samples exceeded 99.9%, indicating near-perfect coverage ( Figure 4A ). We used Shannon’s diversity index (Shannon), the Chao1 index (Chao 1), and the abundance Coverage-Based Estimator (ACE) indices, which are widely used in ecology to assess fungal and bacterial abundance. Our findings show that the Shannon and Chao 1 indices did not exhibit significant variation between the transgenic and non-transgenic cotton soils ( Figures 4B, C ). The ACE fungal level was slightly higher than L1 in NT but not different from L2 ( Figure 4D ). Beta diversity analysis revealed variability between L2 and NT and L1 in the fungal community, resulting in reduced diversity, while L2 exhibited reduced bacterial diversity, but overall, the magnitude of change was not significant ( Figure S1 ). In conclusion, our findings suggest that there were no substantial changes in the diversity of soil microorganisms between transgenic and non-transgenic cotton samples.

The statistical analysis revealed that the L1, L2, and NT groups had 2422, 2431, and 2982 OTUs for fungal species, respectively, and 1777 OTUs were common to all groups ( Figure 5A ). In the case of bacterial OTUs, the statistical analysis revealed that L1, L2, and NT had 6726, 6682, and 6844 OTUs, respectively, and 5737 OTUs were common to all groups ( Figure 5B ). A statistical analysis was conducted at the phylum level for fungi and bacteria in the three groups. The results showed that the three groups had similar dominant species at phylum level in both fungi and bacteria. The fungal phyla Ascomycota, Mortierellomycota, Basidiomycota, Kickxellomycota, and Mucoromycota were widely distributed across all groups ( Figure 5C ). Similarly, Proteobacteria, Actinobacteria, Acidobacteria, Chloroflexi, and Bacteroidetes were relatively dominant among bacteria ( Figure 5D ). To conduct further investigation on the phylogenetic relationships and abundance at the genus level, representative sequences for the most abundant genera were aligned through multiple sequences, subsequently constructing a phylogenetic tree. The analysis of the bacterial and fungal genera in three groups, showed that the dominant genera were similar across all three groups, and all from the dominant phylum-level category. Specifically, Mortierella, Cephalotrichum, Fusarium, Alternari, and Penicillium were the main soil fungal genera ( Figure 5E ), while Arthrobacter, Acidobacteria, Sphingomonas, Pseudomonas, Bacillus, and Skermanella were the dominant soil bacterial genera ( Figure 5F ).

To uncover variations in dominant species among the three sample groups at different taxonomic levels (phylum, genus, and species), we identified the top 10 most abundant species for each level. These species were used to create a ternary plot, which allows for clear visualization of the differences in dominant species across the three sample groups at each taxonomic level. Olpidiomycota and Blastocladiomycota were found to be dominant in NT, whereas Glomeromycota and Kickxellomycota were more abundant in L1 and L2 ( Figure 6A ). In contrast, there were no significant differences in bacterial phyla abundance among the three groups ( Figure 6B ). At the fungal genus level, Neocamarosporium and Lecanicillium were found to be dominant in NT, while Leptosphaeria, Gibberella, and Chaetomium were found to be dominated in L1 and L2 ( Figure 6C ). The bacterial genera Lachnospiraceae and Lactobacillus showed dominance in relative abundance in NT, whereas Intestinimonas dominated in relative abundance in L1 and L2 ( Figure 6D ). Notably, there were no significant differences in Bacillus abundance among the three groups ( Figure 6D ). Fungal species including Mortierella alpina, Fusarium delphinoides, and Cephalotrichum nanum were identified as dominant among the three groups. Furthermore, Neocamarosporium chichastianum was found to have relative dominance in the group categorized as NT, while Chaetomium longicolleum, Albifimbria verrucaria, Gibberella intricans, and Leptosphaeria sclerotioides were found to have relative dominance L1 and L2 groups ( Figure S2A ). In addition, the bacterial communities were dominated by Pseudomonas frederiksbergensis shared among all three groups. However, Lachnospiraceae bacterium A4 showed dominant relative abundance in the NT group. Similarly, L1 and L2 groups, the Lachnospiraceae bacterium DW17, Pedobacter duraquae, and Stenotrophomonas chelatiphaga were also found to be dominant ( Figure S2B ).

We investigated the correlation between microbial and environmental factors by conducting a redundancy analysis (RDA) of bacteria and fungi on a genus level. The RDA results indicated that 29.76% of the correlation was identified on axis 1 and 23.23% on axis 2 for the fungal genus level. The analysis indicates that available K had a prominent role in explaining axis 1, while pH had a crucial role in explaining axis 2. Among them, available K showed the highest impact on Penicillium, whereas pH had the most significant effect on Alternaria ( Figure 7A ). Our findings also revealed that the RDA analysis carried out at the bacterial genus level explained 28.92% on axis 1 and 19.24% on axis 2. The analysis identified that available P had the highest contribution in Axis 1, whereas total P had the greatest contribution to Axis 2. Available P had a more profound impact on Acidobacteria, whereas total P significantly impacted Skermanella ( Figure 7B ). Spearman correlation analysis and coefficients were employed in this study to determine the correlation between environmental factors and microbial species richness (alpha diversity). Spearman rank correlation was used to examine the inter-variation relationship between environmental factors and species. The purpose of this analysis was to obtain a correlation and determine significant differed values between the two variables. The abundance of fungal genera was significantly influenced by various environmental factors. Specifically, Dactylonectria was significantly affected by available P, organic matter, total N, and total P, whereas Conocybe abundance was significantly affected by available K. Leptosphaeria abundance was significantly influenced by total K, while Gibberella abundance was significantly affected by total P. In the case of Lecanicillium, available P had a significant impact on its abundance, while the abundance of Solicoccozyma was significantly influenced by pH ( Figure 7C ). At the bacterial genus level, total K had a significant impact on Solirubrobacter abundance, while Rhizobiaceae abundance was significantly affected by available P. Available P and pH had a significant impact on Lactobacillus abundance, whereas these environmental factors had no significant effect on Bacillus ( Figure 7D ).

We used the FUNGuild functional prediction software for ITS, and its prediction of guilds allowed us to study fungal functions from an additional ecological perspective. The results showed that the software could not predict the majority of fungal functions. The most abundant predicted functions were plant pathogen, soil saprotroph, and wood saprotroph, while animal pathogen was also prevalent ( Figure 8A ). We predicted the effects of bacterial biochemical cycling processes for different environmental samples using the FAPROTAX functional classification database based on species information. The results indicated that chemoheterotrophy had the highest abundance, and fermentation, nitrate reduction, nitrogen respiration, and denitrification also accounted for a significant proportion of them ( Figure 8B ). In conclusion, our functional predictions of soil microbial communities were quite similar between transgenic and non-transgenic cotton, with no significant differences in microbial community functions observed.