- Open Access
The response of dominant and rare taxa for fungal diversity within different root environments to the cultivation of Bt and conventional cotton varieties
Microbiomevolume 6, Article number: 184 (2018)
Bacillus thuringiensis (Bt) crops have been cultivated at a large scale over the past several decades, which have raised concern about unintended effects on natural environments. Microbial communities typically contain numerous rare taxa that make up the majority of community populations. However, the response of dominant and rare taxa for fungal diversity to the different root environments of Bt plants remains unclear.
We quantified fungal population sizes and community composition via quantitative PCR of ITS genes and 18S rRNA gene sequencing of, respectively, that were associated with Bt and conventional cotton variety rhizosphere soils from different plant growth stages. qPCR analyses indicated that fungal abundances reached their peak at the seedling stage and that the taproots and lateral root rhizospheres of the Bt cotton SGK321 were significantly different. However, no significant differences in population sizes were detected between the same root zones from Bt and the conventional cotton varieties. The overall patterns of fungal genera abundances followed that of the dominant genera, whereas overall patterns of fungal genera richness followed those of the rare genera. These results suggest that the dominant and rare taxa play different roles in the maintenance of rhizosphere microhabitat ecosystems. Cluster analyses indicated a separation of fungal communities based on the lateral roots or taproots from the three cotton varieties at the seedling stage, suggesting that root microhabitats had marked effects on fungal community composition. Redundancy analyses indicated that pH was more correlated to soil fungal community composition than Bt protein content.
In conclusion, these results indicate that dominant and rare fungal taxa differentially contribute to community dynamics in different root microhabitats of both Bt and conventional cotton varieties. Moreover, these results showed that the rhizosphere fungal community of Bt cotton did not respond significantly to the presence of Bt protein when compared to the two conventional cotton varieties.
The Bt gene originating from Bacillus thuringiensis is widely used to confer pest resistance in plants . Cultivation of these Bacillus thuringiensis (Bt) protein-containing crops has increased dramatically in recent years, with 75.4 million hectares of Bt crops planted in 2016, representing over 41% of all GM crops that were planted. Moreover, planting of Bt cotton globally accounted for 22.3 million hectares in 2016, and four countries each grew more than 1.0 million hectares . The large-scale adoption of Bt crops can be attributed to the increased crop productivity and the reduced need for chemical pesticides, thereby resulting in a reduced environmental footprint for agriculture . GM Bt crops produce insecticidal recombinant Cry1Ac protein, typically also in their root tissue. Such recombinant products may thus potentially enter the rhizosphere as an additional nutrient source for the soil microbial community. Consequently, one of the major potential environmental risks associated with Bt plants is their effect on soils and the potential to inhibit non-target organisms , such as fungi. Several studies have demonstrated that Bt proteins expressed by Bt crops enter soils through root exudates and then bind rapidly onto surface-active components, thereby becoming less accessible to microbial degradation but still retaining their insecticidal activity [5, 6]. Consequently, Bt proteins could accumulate in soils and may influence biological and chemical processes in addition to microbial community composition . Fungi living in the rhizosphere can impact health, nutrition, and productivity of plants in agricultural ecosystems [8, 9], and understanding how fungal communities are affected by Bt crop cultivation is an essential aspect of elucidating soil biological processes at work in the rhizosphere by exposure to Bt proteins produced by Bt cotton. Li et al.  observed significant seasonal variation in the abundances and diversity of culturally viable fungi, but no significant differences could be attributable to the long-term cultivation of Bt cotton in the abundances of fungal populations or diversity. However, the communities identified appeared to be less taxonomically rich, and changes in the relative abundance of species were easy to overlook based on plate counts of cultivable organisms .
Natural microbial communities typically contain only a few dominant species, but a very large number of less abundant, rare taxa . These low-abundance organisms can be representatives of novel microbial lineages [13, 14] and play crucial roles in biogeochemical cycles and metabolic fluxes [15, 16]. Rare taxa, in particular, have been reported to play important roles in the rhizosphere communities of Bt and conventional maize varieties . For example, Hua et al.  found that dominant and rare taxa exhibited different ecological roles in acid mine drainage communities. Qi et al.  assessed fungal population dynamics in soils planted with Bt cotton and its conventional parental line using 18S and ITS rDNA gene FLX-pyrosequencing. However, little information is available regarding the response of dominant and rare fungal taxa to root exudates, including Bt protein. Studies on the importance of dominant and rare taxa for fungal diversity in the rhizosphere of GM crops expressing Bt proteins revealed by Illumina MiSeq sequencing contribute to their general environmental risk assessment, which is required to evaluate their safety and sustainable use in agriculture.
The root system plays a large role in the acquisition of resources in natural heterogeneous soil environments. Moreover, it has been proposed that root branching and root system architecture (RSA) play a significant role in determining the quantitative and qualitative composition of exudates, thereby modulating shifts in the diversity of soil microbial communities . Yang and Crowley  observed that rhizosphere bacterial communities are substantially different in different root zones and that rhizosphere communities may be altered by changes in the root exudate composition. The RSA of cotton is typically composed of a primary (tap) root and lateral roots. A few studies have investigated the effects of consecutive cultivation of Bt cotton on soil microbe-mediated enzymatic properties and microbial biomass [21, 22]. However, no studies have investigated these effects on different root zones nor the effects of Cry proteins released in Bt cotton root exudates on fungal community structure in rhizospheres. Hence, investigating the effects of microhabitats within different root zones would be valuable to assess and precisely monitor the impacts of genetically modified plants on agro-ecosystems, and especially for dicotyledonous species.
The aim of this study was to characterize the dynamics of rhizosphere fungal community composition associated with Bt cotton SGK321 at different growth stages. Specifically, we assessed the response of dominant and rare fungal taxa to Bt cotton cultivation within different root environments (taproots and lateral roots) using qPCR of ITS gene abundances and high-resolution sequencing of 18S rRNA genes using Illumina MiSeq sequencing. To appropriately contextualize the potential differences caused by genetic modification, controls included two other cotton varieties that were cultivated in the same field and subjected to the same analyses as the lateral roots and taproots of Bt cotton. A specific evaluation was conducted of the possible links between soil Bt protein content and shifts in fungal community composition.
Experimental design and soil sampling
The transgenic Bt cotton variety SGK321 that produces Bt protein to inhibit bollworms was developed by the Shijiazhuang Academy of Agricultural and Forestry Sciences (China) and approved for commercial use by the Chinese government in 2002. A conventional parental cotton line of SY321 was used for comparison, in addition to one conventional variety XLZ13 that is widely cultivated. Seeds of SGK321 and SY321 were obtained from the Institute of Cotton Research of CAAS (Anyang, China), and those of XLZ13 from Economical Crop Research Institute of Xinjiang Academy of Agricultural Sciences (Urumqi, China). Plants of the three cotton varieties were grown from seeds in adjoining fields under the same environmental conditions and field management strategies. Each field comprised an area of 30 m × 30 m, with three replicates of each treatment arranged in a randomized block experimental design. Ten meters at both ends of each field was not used to avoid marginal effects. The survey was conducted in fields of the Dezhou Academy of Agricultural Sciences (Shandong, China). The field study did not involve endangered or protected species, and the specific location of this study was 37° 21′ N, 116° 20′ E. The fore-crop of the experimental field was a conventional wheat variety, and the wheat straw and litter were cleaned out before planting the cotton plants. The soil was classified as a sandy loam soil with 11.23 ± 0.89 g/kg organic matter, 1.72 ± 0.08 g/kg total nitrogen (N), 68.24 ± 5.78 mg/kg alk-hydr N, 17.54 ± 1.76 mg/kg available P, 109.78 ± 8.78 mg/kg available K, and a pH of 7.22 ± 0.35 (soil to water ratio of 1:2.5). The field experiment was conducted on March 30, 2016, according to the conventional local agriculture operations. In contrast to the two control varieties SY321 or XLZ13, SGK321showed a higher resistance against the cotton worm, which can especially prevent the normal growth and development of cotton worm. Moreover, the lint yield of SGK321 before frost was found to be higher than that of SY321 (or XLZ13), while its growth period is shorter than that of the two control varieties. The taproots and lateral root rhizospheres and bulk soils were sampled four times, corresponding to the major growth stages of cotton: seedling (July), budding (August), flowering (September), and bolling (October). After weeds and litter were removed from the soil surface, two types of soil samples were harvested. The bulk soil was excised from at least 4 mm away from any roots, and soils within 2 mm of the root surface were considered rhizosphere soils . Plants were gently removed from the soils, and rhizospheres were collected by gently shaking the roots to dislodge small adhering soil clumps . Rhizosphere samples were collected from five randomly selected plants in each plot and mixed to form a composite sample. Sampling was conducted in triplicate on each sampling day. The samples were placed on ice in a cooler and transported on the same day to the laboratory, where they were passed through a 2.0-mm sieve. After sieving, a portion of the soil sample was stored at − 80 °C for DNA extraction, and the remaining portion was stored at 4 °C for soil property and Bt protein content analysis.
Analysis of soil physicochemical properties and soil microbial DNA extraction
To assess Bt protein content, up to 0.5 g of soil was added to a 2-mL microcentrifuge tube. One milliliter of extraction buffer (0.1 M Na2CO3, 0.1 M NaHCO3, 5.0 mM EDTA, 50 mM Na4P2O7.10H2O, 0.1%Triton X-100, pH 10) was then added to fill the tube, and the sample was mixed for 30 min at room temperature, followed by centrifugation for 5 min at 4 °C at 12,000g. The supernatant was collected, and Bt protein content was measured using a commercial ELISA kit (QuantiPlate™ Kit for Cry1Ab/Cry1Ac, EnviroLogix Inc., Portland, Maine, USA) according to the manufacturer’s protocol. Soil pH was determined using a pH meter (FE20-FiveEasy™ pH, Mettler Toledo, Germany) at a ratio of 1:2.5 (weight/volume) soil to distilled water. Soil microbial DNA was extracted as previously described by Li et al. .
Quantification of fungal communities
We conducted quantitative PCR to determine population sizes of the rhizosphere and bulk soil fungal communities using the universal fungal ITS gene primers NSI1 and 58A2R . Standard curves were generated using tenfold serial dilutions of a plasmid containing the ITS gene of Aspergillus sp. isolated from detritusphere. Each PCR reaction contained 5 μl × 2 SYBR qPCR Premix Ex Taq™ (TaKaRa Biotechnology (Dalian) Co., Ltd.), 0.25 μl each of 10 mM forward and reverse primers, 0.5 μl of deionized distilled water, and 4 μl of either standard or extracted soil DNA. PCR was performed on a StepOnePlus™ Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) using a program, including initial denaturation at 95 °C for 3 min, followed by 40 cycles at 94 °C for 30 s, 53 °C for 30 s, and 72 °C for 45 s, with a final extension at 72 °C for 5 min. The copy number of fungal ITS rRNA genes was calculated using a regression equation to convert the cycle threshold (Ct) value to a known number of copies in the standards. All of the qPCR reactions were run in triplicate for each soil sample. The average fungal PCR efficiency was 93.3% with an R2 of the standard curves of 0.998.
Fungal 18S rRNA gene Illumina MiSeq sequencing and sequence processing
A total of 72 independent samples that were obtained from the taproot and lateral root rhizospheres from the three cotton varieties (SGK321, SY321, and XLZ13) were selected for Illumina MiSeq sequencing of 18S rRNA genes (Additional file 1: Table S1). Fungal 18S rRNA genes were amplified using the primers SSU0817F (5′- TTA GCA TGG AAT AAT RRA ATA GGA -3′) and SSU1196R (5′- TCT GGA CCT GGT GAG TTT CC -3′) on a PCR thermal cycler system (GeneAmp 9700, ABI, USA). PCR reactions were conducted using the following program: 3 min of denaturation at 95 °C, followed by 35 cycles of 95 °C for 30 s; annealing at 53 °C for 30 s; and elongation at 72 °C for 45 s, followed by a final extension at 72 °C for 10 min. PCR reactions were performed in triplicate in 20-μl reaction mixtures containing 4 μl of × 5 FastPfu buffer, 2 μl of 2.5 mM dNTPs, 0.8 μl of each primer (5 μM), 0.2 μl of BSA, 0.4 μl of FastPfu polymerase, 10 ng of template DNA, and deionized distilled water for a total volume of 20 μl. The resultant PCR products were extracted from a 2% agarose gel and further purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA), followed by quantification using a QuantiFluor™-ST (Promega, USA) system according to the manufacturer’s protocols. Purified amplicons were pooled in equimolar concentrations and paired-end sequenced (2 × 300 bp) on the Illumina MiSeq platform (Illumina, San Diego, USA) according to the standard protocols at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China). A total of 12 independent bulk soil samples which were included in this study as the control group were also selected for Illumina MiSeq sequencing. All of the fungal 18S rRNA gene sequence data of rhizosphere and bulk soil analyzed in this study have been deposited into the NCBI Sequence Read Archive (SRA) database under the accession numbers SRP126911 and SRP154528, respectively.
Raw FASTQ files were de-multiplexed, quality filtered, and analyzed with QIIME Pipeline Version 1.7.0 (http://qiime.org/tutorials/tutorial.html) , and the paired reads were joined with fast length adjustment of short reads (FLASH) software  using the following criteria: (i) reads were truncated at any region receiving an average quality score < 20 over a 50-bp sliding window, (ii) reads containing > 2 nucleotide mismatches to primers were discarded along with reads containing ambiguous bases, and (iii) sequences with > 10 bp overlap were merged according to their overlapping sequence. Quality-filtered reads were then clustered into operational taxonomic units (OTUs) at a 97% similarity cutoff using the UPARSE algorithm (version 7.1 http://drive5.com/uparse/). Additionally, chimeric sequences were identified and removed using UCHIME . The taxonomy of each 18S rRNA gene sequence was assessed using the RDP Classifier algorithm (http://rdp.cme.msu.edu) and the Silva (123/18) 18S rRNA database with a confidence threshold of 70%.
Alpha diversity indices (Good’s coverage, Shannon diversity index, and Simpson diversity index) were calculated in QIIME (version 1.7.0) . Statistically significant differences in soil pH, Bt protein content, fungal alpha diversity, ITS rRNA gene copies, and relative abundances of different fungal taxa among samples at each sampling time were examined using one-way analysis of variance (ANOVA) tests in SPSS 19.0 (SPSS Institute, Inc., 2010). Significant differences were considered as P < 0.05. Dominant, rare, and responsive fungal taxa at different taxonomic ranks were identified according to the methods described by Dohrmann et al. . An analysis of similarities (ANOSIM) was performed using QIIME 1.8.0 software  to determine whether different group samples had significantly different microbial communities. Principal component analysis (PCA) based on Bray-Curtis similarities was performed using the R software package (https://www.r-project.org/). Redundancy analysis (RDA) was used to determine the environmental variables that most correlated to fungal community compositional differences, and the results were used to construct a soil property matrix for variation partitioning analysis in R v.2.8.1 using the “vegan” package (v.1.15-1) . Cluster analysis of fungal communities was conducted based on Bray-Curtis dissimilarities using the “picante” and “vegan” packages in the R environment (R Development Core Team, 2006).
Fungal population sizes in cotton rhizospheres and bulk soils
The abundances of fungi in the 72 rhizosphere and 12 bulk soil samples ranged from 1.90 × 107 to 10.71 × 107 ITS rRNA gene copy numbers per gram dry soil (Fig. 1), with peak sizes at the seedling stage. Copy numbers in SGK321 and SY321 (or XLZ13) were not significantly different at the budding, flowering, and bolling stages. Notably, the population size of rhizosphere fungi associated with the taproot of SGK321 at the seedling stage was significantly higher compared to the lateral roots, suggesting that the root microhabitat was influential. There were no significant differences between fungal community sizes between rhizospheres of SGK321 and SY321 from the same root zones, suggesting that fungal abundance was not affected by the production of Bt proteins in the root tissue of SGK321. Moreover, the fungal population sizes of rhizosphere soil within the three cotton varieties were significantly higher than those of bulk soil at seedling, budding, and flowering stages, but no significant difference was noted between rhizosphere and bulk soils at bolling stage.
Taxonomic assignment of the 18S rRNA gene sequences of rhizosphere fungi
A total of 1907,500 high-quality 18S rRNA gene sequences were generated by Illumina MiSeq sequencing of the 72 samples, with 16,635–35,015 sequences obtained per sample (Additional file 1: Table S2A, B, C, D). When clustered at the 97% nucleotide similarity level, 304 different phylotypes (OTUs) were present among the soils, and richness was significantly higher at the budding and flowering stages (Fig. 2a). OTU richness was significantly different between the taproot and lateral root rhizospheres of SGK321 at the seedling stage, but no significant differences were observed when comparing the same root zones between SGK321 and SY321. Moreover, there were no significant differences in OTU richness between different varieties and root zones at other stages. The 26 most abundant OTUs (> 0.5% relative abundance) comprised sequences from all of the 72 samples (Additional file 1: Table S3A, B, C, D), and 25 of these OTUs were present in all of the samples. Rarefaction curves of all soil samples reached their asymptotes, indicating that the data generated in this study was enough for the analysis of the diversity of fungi in the rhizosphere (Additional file 1: Figure S1). A total of 40 fungal phyla and 127 genera were identified throughout the growth stages, with the predominant phyla being Ascomycota (83.31%), Basidiomycota (7.81%), Zygomycota (4.40%), unclassified Fungi (1.96%), and unclassified Eukaryota (0.69%), among others. The predominant genera were unclassified Pezizales (27.28%), unclassified Sordariales (11.42%), Fusarium (9.29%), unclassified Hypocreales (5.90%), and unclassified Ascomycota (5.49%), among others (Fig. 3a, b). The six most abundant phyla (each accounting for > 0.05% of fungal sequences) represented 98.71% of the total sequence data, while the remaining 34 rare phyla comprised only 1.29% of the total sequence data. Likewise, the 21 most abundant genera represented 94.90% of the total sequence data, while the remaining 106 rare genera constituted only 5.10% of the total sequence data. Thus, the trend of few dominant taxa to the group with high sequence load and most of the rare taxa to the group with low sequence load could be found at different taxonomic ranks, and the relative sequence of rare taxa increased gradually at the lower hierarchical ranks (Fig. 4).
Diversity of rhizosphere fungi
An ANOSIM test (R2 = 0.8092, P = 0.001) (Additional file 1: Figure S2) showed that the difference between the groups was significantly greater than that within the group, suggesting the results obtained from this study were reliable. The library coverage, C, based on OTU richness was > 99.93% (Fig. 2b), suggesting that the vast majority of all taxonomic units were detected in 72 samples, including in the sample SGK321_T41, that had the lowest number of sequences of any sample (16,635; 0.96%). The alpha diversity indices comprising the Shannon index H′ and Simpson index D of the fungal communities significantly differed among rhizosphere communities at the different growth stages of cotton (Fig. 2c, d). Specifically, H′ was significantly higher at the budding and flowering stages, while D was significantly higher at the seedling and bolling stages. However, the H′ and D values were not significantly different at any of the sampling time points when comparing the taproot and lateral root rhizospheres of SGK321, SY321, and XLZ13. Thus, these results indicate that the variation in diversity of fungal communities was determined by cotton growth stages but not by varieties, including Bt cotton SGK321.
Effect of Bt cotton cultivation on dominant and rare fungal taxa
The richness and abundance of all of the genera, the dominant genera, and the rare genera subsets of the fungal communities changed significantly at the different growth stages (Table 1). The richness of all of the genera in the lateral root and taproot soil communities of SGK321 was significantly different between the seedling and bolling stages. However, comparisons of SGK321 and SY321 indicated no differences at the same stage for the same root zones, and similar results were observed for rare genera. Notably, the abundances of all of the genera, the dominant genera, and the rare genera were not different throughout the growth stages for the three cotton varieties. When compared with the richness and abundance of SY321 lateral root of different growth stages, the richness and abundance changes of all of the genera, the dominant genera, and the rare genera of fungal communities of different cotton varieties and root zones are shown in Table 2. The richness of rare genera and all of the genera correlated similarly to the cultivation of Bt cotton, suggesting that changes in fungal community richness were mainly due to the effects arising from the variety and the root zones on the richness of rare genera. In addition, similar progressions were observed for the abundance distribution of all of the genera and the dominant genera at different growth stages, indicating that the changes in fungal community abundances at different stages were mainly derived from the effects arising from variety type and the root zones on the abundances of dominant genera. Consequently, rare genus abundances were more responsive to external interference when compared to the dominant genera. In order to assess whether distinct taxa were preferred among the dominant genera, the class-level taxonomic assignments of the 21 dominant genera, 106 rare genera, and 120 responsive genera were compared with that of all of the 127 genera (Table 3). The proportion of contributing genera (> 7 genera of one class) among all of the genera was similar to that of the responsive genera for most classes. The dominant genera were generally preferred to associate with the Sordariomycetes and Agaricomycetes, while the rare genera were generally preferred to associate with the Intramacronucleata. In addition, the proportion of contributing genera belonging to the Dothideomycetes was in accordance with that of the dominant genera.
Cluster analysis of rhizosphere fungal communities of different roots environment
Cluster analyses of the relative abundances of all of the genera indicated that the lateral root and taproot rhizospheres of the conventional variety XLZ13 clustered together. These results highlighted the distinct fungal community compositions of XLZ13 compared to rhizospheres from the GM variety (SGK321) or the parental cotton variety (SY321) (Fig. 5). Separation of fungal communities from lateral roots or taproots was observed at the seedling stage of the three varieties, suggesting an impact of the root environment on fungal community composition. The lack of similar separation at other growth stages suggests that fungal communities were also affected by other factors, such as climate or seasonality, for example. Rhizosphere samples from SGK321 and SY321 were highly similar, and communities from all of the growth stages clustered closely. Thus, Bt cotton, such as SGK321, did not have specific effects on fungal community composition. When cluster analyses of the dominant, rare, and responsive genera subsets were conducted (Additional file 1: Figures S3-S5), similar clustering patterns were observed between the dominant and all of the genera subsets, suggesting that the dominant genera primarily contributed to the overall composition dynamics of the communities. In contrast, differences based on the rare genera subset resulted in clustering that was highly distinct from the clustering based on all of the genera. The fungal communities of the soil samples based on the responsive genera clustered roughly into six groups corresponding to different treatments (cotton varieties and root zones). A total of 120 genera significantly responded in its abundance to one variety in comparison with its abundance in all other varieties or to the root environment at different growth stages.
Relationship between rhizosphere fungal community composition and soil pH and Bt protein
The relationship between fungal community composition and soil pH and Bt protein content was analyzed by RDA (Fig. 6a). The first two axes of the RDA accounted for 23.96% of the total variance in fungal community composition, with the first axis accounting for 20.55% of the variance. The distribution pattern of the samples from the seedling and bolling stages was discrete, while that at the budding and flowering stages clustered together. Soil pH was more correlated to community composition than soil Bt protein content, suggesting more important effects arising from soil pH. This conclusion was also supported by a Mantel test, indicating that soil pH (r = 0.700, P = 0.001) was more correlated to fungal community composition than Bt protein content (r = 0.417, P = 0.001). We also used Spearman’s rank order correlations to investigate the response of specific fungal genera over time to soil pH and Bt protein content (Fig. 6b). The abundances of unclassified Zygomycota (r = 0.804, P = 0.001), unclassified Fungi (r = 0.764, P = 0), and unclassified Eukaryota (r = 0.784, P = 00) were positively correlated with soil pH. In contrast, unclassified Nectriaceae (r = − 0.605, P = 0.001) abundances were negatively correlated with soil pH. Lastly, the abundances of no rank Agaricomycetes (r = 0.655, P = 0) and unclassified Agaricomycetes (r = 0.605, P = 0) were positively correlated with soil Bt protein content.
Differences in the variation of dominant and rare fungal taxa in relation to the cultivation of Bt and conventional cotton varieties
Microorganisms present in the rhizosphere are affected by root exudates and play important roles in the growth and ecological fitness of their plant host. Rhizosphere microorganisms are considered to be an important component of soil ecological system. Consequently, the microbial community and their associated dynamic processes are the foundations of ecosystems. However, the importance of rare, low-abundance microorganisms is largely unknown, although it is hypothesized that they contribute to the community stability by acting as a reservoir that can rapidly respond to environmental changes . Fungi are ubiquitous microorganisms that play important roles in soil ecosystems as major decomposers of organic matter, and they release nutrients during nutrient cycling that stimulate plant growth . Thus, it is critical to assess the contribution of dominant and rare fungi toward root exudates including Bt protein in light of the diversity and processes of fungal communities. To explore this issue, we employed deep sequencing of uncultured fungal communities to investigate the response of less abundant, rare fungal community members to the cultivation of Bt cotton plants.
Here, the community structure of the dominant genera obtained the similar clustering pattern with all genera throughout all the growth stages, which is similar to the variation of the bacterial community . In contrast, the communities represented by the rare genera clustered differently, and this was consistent with the community clustering due to the responsive genera. Previous researches have investigated microbial indicators for the purpose of evaluating correlations between microbes and for identifying possible sources of specific environment variation [33,34,35]. However, the majority of these researches lacked sampling power to comprehensively identify the indicators and assess the dynamic nature of indicators due to the use of low-throughput sequencing. We used high-throughput sequencing on the Illumina platform to better understand the rare fungal taxa, which are more sensitive to disturbances. Only 20 of the 120 responsive genera identified in this study were dominant genera, whereas 100 were rare. The detection of responsive genera associated with Bt cotton cultivation provides targets to identify the fungal taxa potentially responsible for degrading and utilizing Bt protein in the rhizosphere. Furthermore, our analyses revealed responsive fungal taxa that may be vulnerable to environmental change, with cotton plant varieties and/or root zones acting as proxies for environmental differences. These results suggest that some of these taxa may be indicators of unmeasured environmental changes and thus provide targets to assess biological drivers that are responsive to subtle physical, chemical, or biological changes. Finally, we found that the response of dominant genera abundances to differences in cotton varieties and/or root zones was consistent with that of all of the genera at different growth stages, while the abundances of rare genera displayed distinctive clustering trends. These results indicate that changes of abundance in the fungal community were caused by the effects of aforementioned treatments on the dominant genera. Indeed, the abundances of rare genera were markedly susceptible to environmental changes, which are similar to the observations from Shade et al. . Compared to the remarkable change of abundances throughout the growth stages, the richness of dominant genera changed only slightly. Further, the variation in richness values for the rare and all of the genera subsets was similar, suggesting that changes in community richness were mainly due to the effects of treatments on the rare genera. Although this study indicated the importance of dominant and rare taxa for microbial diversity in the rhizosphere of Bt and conventional cotton varieties, the metabolic potentials and ecological roles of rhizosphere dominant and rare microbes to the consecutive cultivation of Bt cotton revealed by metagenomics and metatranscriptomics need to be further analyzed in a subsequent experiment.
Differences in fungal community structure with different root zones and possible factors regulating fungal distribution in roots
The plant rhizosphere is a dynamic environment in which many factors may affect the structure and species composition of microbial communities that colonize the roots . These factors can include host plant species  or cultivars  grown in the same soil, and even different root zones of the same plant . Here, we sampled lateral roots and taproots at all of the sampling times, and we collected soils in triplicate in order to minimize sample variation but reflect overall differences associated with root zones. The abundances of ITS rRNA gene copies were clearly affected by (i) cotton species, (ii) growth stages, and (iii) root zones. In particular, the lateral root of SGK321 exhibited lower fungal abundances than taproots. Notably, Zhang et al.  reported that the rhizosphere had a lower fungal diversity (observed OTUs and Chao1) than bulk soil and a distinct fungal community structure in the rhizosphere compared with bulk soil. In this study, during the periods of seedling, budding, and flowering stages, fungal abundance of rhizosphere soil revealed higher than that of bulk soil, but the fungal abundance decreased in rhizosphere soil at bolling stage, without a significant difference with that in bulk soil. These results are consistent with those of Marchand et al. , who found that different microbial population were obtained in rhizosphere and bulk soils during the vegetative growth phase, but there is no significant difference of diversity of culturable microflora between rhizosphere and the bulk soil during the reproductive growth stage. Consequently, in the present study, we speculated that the remarkable decrease of root exudates due to the decline of root function at the bolling stage is unable to alter soil fungal abundance significantly. Moreover, PCA at OTU level (Additional file 1: Figure S6) showed that the fungal community structure in rhizosphere soil had significant differences compared with the bulk soil at the different growth stages. In this study, cluster analysis indicated that the conventional variety, XLZ13, had a different fungal community composition compared to the Bt cotton variety, SGK321, or its conventional parental variety, SY321, throughout the various growth stages. These results suggest that the variation introduced by different conventional cotton varieties might be greater than those between transgenic and conventional cotton varieties. Moreover, lateral roots or taproots of the three cotton varieties grouped separately at the seedling stage, indicating that the rhizosphere fungal community was highly influenced by root microhabitats during this growth period. Importantly, previous studies centered on Bt cotton [22, 40,41,42], and it is indispensable to assess and monitor the impacts of Bt cotton cultivation on agroecosystems according to the case-by-case principle of biosafety assessment. Thus, to thoroughly assess the effects of genetically modified dicotyledons (cotton, soybean, rapeseed, among others) on soil microbial communities, particular root zones should be considered carefully, because they play significant roles in determining the quantitative and qualitative composition of exudates. Badri and Vivanco  suggested that root-secreted chemicals mediate multi-partite interactions in the rhizosphere and that root exudates initiate and modulate dialog between roots and soil microorganisms. Thus, we hypothesized that wide variation in the quantity and quality of root exudates, including Bt protein content, between lateral roots and taproots might affect the dynamics of microbial communities, including nutrient cycling. Consequently, the effect of organic compounds produced in the lateral roots and taproots on selecting microbial community members warrants further investigation. It should also be noted that different fungal populations between the taproots and lateral roots may play specific roles in the different root zones and in the overall ecological fitness of their plant host.
Cultivation of Bt cotton alters Bt protein content and soil pH associated with shifts in soil fungal communities
Recent research has demonstrated that Bt proteins encoded by cry genes from Bacillus thuringiensis are released in root exudates from Bt plants [43, 44]. However, it is not yet clear whether Bt proteins in root exudates significantly influence soil microbial communities. In this study, the Bt protein content from SGK321 was significantly higher than that of SY321 and XLZ13 when considering the same root zones. Notably, this study observed differences of Bt protein from the lateral roots and taproots of Bt cotton at the seedling and budding stages for the first time. These differences in Bt protein content may be due to the differences in Bt protein content within the cotton taproot and lateral root tissues. Hence, we measured Bt protein content in different root tissues zones. However, the values were not significantly different between the taproot and lateral root tissues (Additional file 1: Figure S7). Parker et al.  described that root hairs comprised as much as 77% of the total root surface area of cultivated crops, forming the major point of contact between the plant and rhizosphere. Accordingly, the enlarged interface between cotton and rhizosphere due to the presence of more root hairs in the lateral roots may result in marked alterations of rhizosphere Bt protein contents in different root location. Although RDA indicated that the soil fungal communities were significantly correlated to Bt protein content (range, 54.43 to 303.94 pg/g soil; r = 0.417, P = 0.001) through all of the growth stages, the correlation was not significant when analyzed for individual growth stages. This finding implies that the shifts in the composition of soil fungal communities were due to the effects of the cotton growth stage rather than the variation in Bt protein content. Soil pH is one of the most influential chemical factors affecting soil microbial communities, as it affects numerous abiotic factors, including carbon and nutrient availability . To explain the pH variability within the rhizosphere between SGK321 and SY321 (or XLZ13) at the seedling and budding stages, we cultivated the varieties in sterile Hoagland’s nutrient solution. The total amount of organic acids in root exudates of SGK321 was dramatically lower than that of SY21 and XLZ13 (Additional file 1: Table S4). Thus, the higher rhizosphere pH value of SGK321 might be related to pronounced organic acid decreases in the root exudates, in comparison to the control varieties. Further, soil pH was the strongest predictor of soil fungal community composition, followed closely by Bt protein content, suggesting that Bt protein introduced by root exudates did not impact fungal community composition more strongly than soil pH. Notably, several studies have shown that Bt protein expressed by Bt crops does not persist in the soil after long-term cultivation of Bt crops [10, 47, 48]. Moreover, environmental factors and cotton growth traits are specific at different growth stages, which affect the diversity and dynamics of rhizosphere microbial communities differently. Accordingly, we deduced that these findings obtained from rhizosphere and bulk soil samples at four different growth stages of cotton in this study were probably consistent with the research of the longer-term cultivation of Bt cotton.
In summary, we monitored the diversity and dynamics of the rhizosphere fungal community from lateral roots and taproots of Bt plants for the first time using qPCR and high-throughput sequencing approaches. The dominant and rare fungal community taxa differentially responded to the cultivation of Bt cotton. Furthermore, fungal community structure was substantially affected by different root zones at the seedling stage, while the Bt cotton cultivar had no specific effect on fungal diversity in comparison to the conventional cotton varieties. Our results indicate that soil pH more strongly affected fungal community structure compared to soil Bt protein content. However, further investigations are needed to assess the variation in the organic compound composition of root exudates that were observed between the lateral roots and taproots of the Bt cotton variety and the consequent interactions between root exudates and fungal communities.
Benedict J H, Ring D R. Transgenic crops expressing Bt proteins: current status, challenges and outlook. In: Koul O, Dhaliwal GS. Transgenic crop protection: concepts and strategies. Oxford University Press and IBH Publishing Co, New Delhi, 2004.
Clive J. Global status of commercialized biotech/GM crops:2016. ISAAA Brief 52. Ithaca: ISAAA; 2016. p. 4–8.
Michael SK, Jason MW, Steven LL, James AB, John LV, Bruce GH. The food and environmental safety of Bt crops. Front Plant Sci. 2015;6:283.
Dohrmann AB, Küting M, Jünemann S, Jaenicke S, Schlüter A, Tebbe CC. Importance of rare taxa for bacterial diversity in the rhizosphere of Bt- and conventional maize. ISME J. 2013;7:37e49.
Vettoria C, Paffettib D, Saxenac D, Stotzkyc G, Giannini R. Persistence of toxins and cells of Bacillus thuringiensis subsp. kurstaki introduced in sprays to Sardinia soils. Soil Biol Biochem. 2003;35(12):1635–42.
Valldor P, Graff RM, Martens R, Tebbe CC. Fate of the insecticidal Cry1Ab protein of GM crops in two agricultural soils as revealed by 14C-tracer studies. Appl Microbiol Biotechnol. 2015;99:7333–41.
Singh AK, Singh M, Dubey SK. Rhizospheric fungal community structure of a Bt brinjal and a near isogenic variety. J Appl Microbiol. 2014;117:750–65.
Berendsen RL, Pieterse CM, Bakker PA. The rhizosphere microbiome and plant health. Trends Plant Sci. 2012;17:478–86.
Philippot L, Raaijmakers JM, Lemanceau P, van der Putten WH. Going back to the roots: the microbial ecology of the rhizosphere. Nat Rev Microbiol. 2013;11:789–99.
Li XG, Liu B, Cui JJ, Liu DD, Ding S, Gilna B, Luo JY, Fang ZX, Cao W, Han ZM. No evidence of persistent effects of continuously planted transgenic insect-resistant cotton on soil microorganisms. Plant Soil. 2011;339:247–57.
Epstein SS. The phenomenon of microbial uncultivability. Curr Opin Microbiol. 2013;16:636–42.
Sogin ML, Morrison HG, Huber JA, Mark Welch D, Huse SM, Neal PR, Arrieta JM, Tourna M, Freitag TE, Nicol GW, Prosser JI. Growth, activity and temperature responses of ammonia-oxidizing archaea and bacteria in soil microcosms. Environ Microbiol. 2008;10:57–64.
Castelle CJ, Hug LA, Wrighton KC, Thomas BC, Williams KH, Wu DY, Tringe SG, Singer SW, Eisen JA, Banfield JF. Extraordinary phylogenetic diversity and metabolic versatility in aquifer sediment. Nat Commun. 2013;4:1–10.
Kantor RS, Wrighton KC, Handley KM, Sharon I, Hug LA, Castelle CJ, Thomas BC, Banfield JF.Small genomes and sparse metabolisms of sediment-associated bacteria from four candidate phyla. mBio. 2013;4:e00708–13.
Musat N, Halm H, Winterholler B, Hoppe P, Peduzzi S, Hillion F, Horreard F, Amann R, Jørgensen BB, Kuypers MM. A single-cell view on the ecophysiology of anaerobic phototrophic bacteria. Proc Natl Acad Sci U S A. 2008;105(46):17861–6.
Wrighton KC, Thomas BC, Sharon I, Miller CS, Castelle CJ, VerBerkmoes NC, Wilkins MJ, Hettich RL, Lipton MS, Williams KH, Long PE, Banfield JF. Fermentation, hydrogen, and sulfur metabolism in multiple uncultivated bacterial phyla. Science. 2012;337(6102):1661–5.
Hua ZS, Han YJ, Chen LX, Liu J, Hu M, Li SJ, Kuang JL, Chain PS, Huang LN, Shu WS. Ecological roles of dominant and rare prokaryotes in acid mine drainage revealed by metagenomics and metatranscriptomics. The ISME Journal. 2014;9(6):1–15.
Qi XM, Liu B, Song QX, Zou BJ, Bu Y, Wu HP, Ding L, Zhou GH. Assessing fungal population in soil planted with Cry1Ac and CPTI transgenic cotton and ITS conventional parental line using 18S and ITS rDNA sequences over four seasons. Front Plant Sci. 2016;7:1023.
Badri DV, Vivanco JM. Regulation and function of root exudates. Plant, Cell and Environment. 2009;32:666–81.
Yang CH, Crowley DE. Rhizosphere microbial community structure in relation to root location and plant iron nutritional status. Appl Environ Microbiol. 2000;66(1):345–51.
Chen ZH, Wei K, Chen LJ, Wu ZJ, Luo JY, Cui JJ. Effects of the consecutive cultivation and periodic residue incorporation of Bacillus thuringiensis (Bt) cotton on soil microbe-mediated enzymatic properties agriculture. Ecosystems and Environment. 2017;239:154–60.
Luo JY, Zhang S, Zhu XZ, Lu LM, Wang CY, Li CH, Cui JJ, Zhou ZG. Effects of soil salinity on rhizosphere soil microbes in transgenic Bt cotton fields. J Integr Agric. 2017;16(7):1624–33.
DeAngelis KM, Brodie EL, DeSantis TZ, Andersen GL, Lindow SE, Firestone MK. Selective progressive response of soil microbial community to wild oat roots. ISME J. 2009;3:168–78.
Li P, Ye SF, Liu H, Pan AH, Ming F, Tang XM. Cultivation of drought-tolerant and insect-resistant rice affects soil bacterial, but not fungal, abundances and community structures. Frontier in microbiology. 2018;9:1390.
Li P, Dong JY, Yang SF, Bai L, Wang JB, Wu GG, Wu X, Yao QH, Tang XM. Impact of beta-carotene transgenic rice with four synthetic genes on rhizosphere enzyme activities and bacterial communities at different growth stages. Eur J Soil Biol. 2014;65:40–6.
Mitchell JI, Zuccaro A. Sequences, the environment and fungi. Mycologist. 2006;20:62–74.
Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, Fierer N, Pena AG, Goodrich JK, Gordon JI, Huttley GA, Kelley ST, Knights D, Koenig JE, Ley RE, Lozupone CA, McDonald D, Muegge BD, Pirrung M, Reeder J, Sevinsky JR, Tumbaugh PJ, Walters WA, Widmann J, Yatsunenko T, Zaneveld J, Knight R. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7:335e336.
Magoč T, Salzberg SL. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics. 2011;27:2957e2963.
Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. UCHIME improves sensitivity and speed of chimera detection. Bioinformatics. 2011;27:2194e2200.
RDC Team. R: a language and environment for statistical computing In: Computing RFfS (ed) Vienna, Austria, 2006. https://www.r-project.org/.
Shade A, Jones SE, Caporaso JG, Handelsman J, Knight R, Fierer N, Gilbert JA. Conditionally rare taxa disproportionately contribute to temporal changes in microbial diversity. Contribute to Temporal Changes in Microbial Diversity mBio. 2014;5(4):e01371.
Lodge JD. Nutrient cycling by fungi in wet tropical forests. In: Isaac S, Franckland JC, Watling R, Whalley AJS, editors. Aspects of tropical mycology. Cambridge: Cambridge University Press; 1995. p. 37–57.
Mendes LW, Kuramae EE, Navarrete AA, van Veen JA, Tsai SM. Taxonomical and functional microbial community selection in soybean rhizosphere. ISME J. 2014;8(8):1577–87.
Shelton JL, Akob DM, McIntosh JC, Fierer N, Spear JR, Warwick PD, McCray JE. Environmental drivers of differences in microbial community structure in crude oil reservoirs across a methanogenic gradient. Front Microbiol. 2016;28(7):1535.
Ahmed W, Staley C, Kaiser T, Sadowsky MJ, Kozak S, Beale D, Simpson S. Decay of sewage-associated bacterial communities in fresh and marine environmental waters and sediment. Appl Microbiol Biotechnol. 2018;102(16):7159–70. https://doi.org/10.1007/s00253-018-9112-4.
Dohrmann AB, Tebbe CC. Effect of elevated tropospheric ozone on the structure of bacterial communities inhabiting the rhizosphere of herbaceous plants native to Germany. Appl Environ Microbiol. 2005;71:7750–8.
Buée M, De Boer W, Martin FM, van Overbeek L, Jurkevitch E. The rhizosphere zoo: an overview of plant-associated communities of microorganisms, including phages, bacteria, archaea, and fungi, and of some of their structuring factors. Plant Soil. 2009;321:189–212.
Watt M, Hugenholtz P, White R, Vinall K. Numbers and locations of native bacteria on field-grown wheat roots quantified by fluorescence in situ hybridization (FISH). Environ Microbiol. 2006;8:871–84.
Zhang KP, Jonathan M, Adams JM, Shi Y, Yang T, Sun RB, He D, Ni YY, Chu HY. Environment and geographic distance differ in relative importance for determining fungal community of rhizosphere and bulk soil. Environ Microbiol. 2017;19(9):3649–59.
Marchand AL, Piutti S, Lagacherie B, Soulas G. Atrazine mineralization in bulk soil and maize rhizosphere. Biol Fertil Soils. 2002;35(4):288–92.
Chen ZH, Wei K, Chen LJ, Wu ZJ, Luo JY, Cui JJ. Effects of the consecutive cultivation and periodic residue incorporation of Bacillus thuringiensis (Bt) cotton on soil microbe-mediated enzymatic properties. Agric Ecosyst Environ. 2017;239:154–60.
Li P, Li YC, Shi JL, Yu ZB, Pan AH, Tang XM, Ming F. Impact of transgenic Cry1Ac + CpTI cotton on diversity and dynamics of rhizosphere bacterial community of different root environments. Sci Total Environ. 2018;637-638:233–43.
Rui YK, Yi GX, Zhao J, Wang BM, Li ZH, Zhai ZX, He ZP, Li Q. Changes of Bt toxin in the rhizosphere of transgenic Bt cotton and its influence on soil functional bacteria. World J Microbiol Biotechnol. 2005;21:1279–84.
Knox OGG, Gupta VVSR, Nehl DB, Stiller WN. Constitutive expression of Cry proteins in roots and border cells of transgenic cotton. Euphytica. 2007;154:83–90.
Parker JS, Cavell AC, Dolan L, Roberts K, Grierson CS. Genetic interactions during root hair morphogenesis in Arabidopsis. Plant Cell. 2000;12:1961–74.
Qin H, Wang HL, Strong PJ, Li YC, Xu QF, Wu QF. Rapid soil fungal community response to intensive management in a bamboo forest developed from rice paddies. Soil Biology Biochemistry. 2014;68:177e184.
Dubelman S, Ayden BR, Bader BM, Brown CR, Jiang C, Vlachos D. Cry1Ab protein does not persist in soil after 3 years of sustained Bt corn use. Environ Entomol. 2005;34:915–21.
Lee ZL, Bu NS, Cui J, Chen XP, Xiao MQ, Wang F, Song ZP, Fang CM. Effects of long-term cultivation of transgenic Bt rice (Kefeng-6) on soil microbial functioning and C cycling. Sci Rep. 2017;7:4647.
We would like to thank Christoph C. Tebbe and Anja B. Dohrmann (Institute of Biodiversity, Johann Heinrich von Thünen-Institute (vTI), Federal Research Institute for Rural Areas, Forestry and Fisheries, Braunschweig, Germany) for the excellent analytical method assistance. We thank Xinjiu Dong (Cash Crop Research Institute of Xinjiang Academy of Agricultural Sciences, China) for supplying the seeds of XLZ13. We thank Yuji Jiang (The Institute of Soil Science, Chinese Academy of Sciences) for the analysis of soil properties. We gratefully acknowledge Jun Chen, Wen Zhai (Yuan Long Ping High-Tech Agriculture Co., Ltd., China), and Yongjin Qiao (Crop Breeding and Cultivation Research Institute, Shanghai Academy of Agricultural Sciences) for the statistical analysis.
This work was financially supported by the National Natural Science Foundation of China (No. 31500461), the SAAS Program for Excellent Research Team No. 2017 (B-07) and Shanghai leading talent plan fund.
Availability of data and materials
The datasets generated and analyzed during this study have been deposited in the GenBank short-read archive SRA accessions SRP126911 and SRP154528.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Terminology of rhizosphere samples included in this study. Table S2. Number of sequences obtained per soil sample. Seedling (A), Budding (B), Flowering (C) and Bolling (D) stages. Table S3. Most abundant OTUs (> 0.5% relative abundance) shared by sequences of at least 97% sequence similarity. Seedling (A), Budding (B), Flowering (C) and Bolling (D) stages. Table S4. Types and quantities of the secreted organic acids in the culture solution (μg/d) (n = 3). Figure S1. Rarefaction curves of all of the rhizosphere samples based on OTUs. Figure S2. Analysis of similarities (ANOSIM) for fungal communities of the rhizosphere samples based on OTUs. Figure S3. Hierarchical cluster analysis of the dominant fungal genera from the rhizosphere of different cotton varieties and root environments at seedling (A), budding (B), flowering (C), and bolling (D) stages, based on the Hellinger distances of microbial communities. Figure S4. Hierarchical cluster analysis of the rare fungal genera from the rhizosphere of different cotton varieties and root environments at seedling (A), budding (B), flowering (C), and bolling (D) stages, based on the Hellinger distances of microbial communities. Figure S5. Hierarchical cluster analysis of the responsive fungal genera from the rhizosphere of different cotton varieties and root environments at seedling (A), budding (B), flowering (C), and bolling (D) stages, based on the Hellinger distances of microbial communities. Figure S6. PCA analysis of fungal community at seedling (A), budding (B), flowering (C), and bolling (D) stages based on OTU level. BS indicates bulk soil control. Figure S7. Bt protein contents of the different root tissues collected at different growth stages of the cotton varieties. (DOC 5941 kb)