Combined metagenomic and metabolomic analyses reveal that Bt rice planting alters soil C-N metabolism

Bt rice led to the redistribution of soil nitrogen

To characterize the influence of Bt rice on soil environmental biochemistry, samples were first separated into two portions including soils and surface waters. Bt proteins were not detected in surface waters from all cultivars (Supplemental Table S2). However, Bt protein contents for rhizospheres from all three cultivars and bulk soils ranged between 64.14 and 126.68 pg/g soil (Supplemental Table S3). Bt protein contents in samples from Bt rice grown in IRRI rice nutrient solution reached 850 pg/ml (Supplementary Table S1). We speculated that the vast majority of Bt protein released from Bt plants was bound tightly to soil particles and was thus difficult to isolate, purify, and detect. Total N, NH₄⁺-N, NO₃⁻-N, and NO₂⁻-N contents in T1C-1 rhizospheres were significantly higher than in the Minghui 63 rhizospheres, although the soil pH of T1C-1 rhizospheres was also significantly lower than for Minghui 63 soils (Supplemental Table S3). Interestingly, the total N, NH₄⁺-N, and NO₃⁻-N contents in the Zhonghua11 rhizospheres were significantly higher than in the Minghui 63 rhizospheres, pointing to an apparent impact of genotypic differences from different conventional cultivars on soil nitrogen. No differences in organic matter and total P contents were identified among all soil samples (Supplemental Table S3). In addition, the surface waters of T1C-1 exhibited higher NO₃-N contents than Minghui 63 soils, but lower pH values than Minghui 63 (Supplemental Table 2), consistent with soil results.

Bt rice altered soil microbial communities, but not surface water communities

Soil and surface water samples were collected and analyzed to characterize metagenomic profiles associated with different cultivars. A total of 11,529,157 and 2,880,919 genes were obtained for soil and surface water samples, respectively (Supplementary Table S4). The α diversity indices, Shannon–Wiener index (H’), Simpson index (D), and Evenness (E) were significantly higher in soils than in surface waters, but significant differences were not observed for Richness (R) (Fig. 2A). Except for R, the α diversity indices E, H′, and D were significantly higher in the T1C-1 rhizosphere than in the other samples, suggesting that Bt rice increased soil microbial diversity rather than altering taxonomic compositions. Differences in α diversity indices were not observed among all of the surface water samples (Supplementary Table S5). Principal coordinates analysis (PCoA) (Fig. 2B) based on microbial taxonomic level (genera) and functional classifications (clusters of orthologous groups of proteins, COG) indicated that soil samples from different rice cultivars and bulk soils formed distinct clusters in ordination space. These distinct groupings were not observed for surface water samples, suggesting that Bt rice cultivation altered soil microbial community composition and functions, but these changes did not occur in surface waters. The rhizospheres of T1C-1, Minghui 63, and Zhonghua 11 shared substantial overlap in total genera (Supplementary Fig. S2A). In addition, 40 genera specifically inhabited T1C-1 rhizospheres (Supplementary Fig. S2B). To further identify taxa that were differential between T1C-1 and Minghui 63 soils, the 50 most abundant genera that were differentially abundant for T1C-1 or Minghui 63 were specifically analyzed using a T-test. Among these, 33 were elevated in T1C-1 soils compared with Minghui 63 soils (Supplementary Fig. S3). Thus, the strongest enrichment was observed for taxa in T1C-1 soils, which is consistent with the general increased α diversity indices for T1C-1 communities (Supplementary Table S5).

Network analysis can provide comprehensive insights into microbial community dynamics and help identify co-occurrence patterns at phylogenetic levels. Co-occurrence network analysis was performed based on the 50 most abundant microbial genera in soils and surface waters. Soil networks consisted of 48 nodes (i.e., genera) and 926 edges (Fig. 2C), while surface water networks comprised 39 nodes and 154 edges (Fig. 2D). The average network transitivity was higher in soil (0.7355) than in surface water (0.4588) networks. In contrast, the average shortest path length of soil samples (1.7243) was lower than for surface water samples (3.6505). These properties suggested that the soil microbial community network was significantly larger and more complex than that of the surface waters. The average modularity and negative correlations in microbial co-occurrence networks were evaluated to assess soil microbiome stability across plant cultivar treatments. Soil microbial networks comprised a larger proportion of negative correlations in T1C-1 soils compared to the Minghui 63 network (Supplementary Fig. S4). The presence of negative interactions contributes to the stability of co-oscillations in communities and promotes the stability of networks. Modular structure analysis indicated that nodes in the networks of T1C-1, Minghui 63, Zhonghua11, and bulk soils grouped into three major modules, respectively, indicating that T1C-1 soils did not exhibit decreased network modularity (Supplementary Fig. S5). Candidatus Accumulibacter, Haliangium, Rhodoplanes, Variovorax, Lysobacter, Dongia, Ramlibacter, etc. were identified as keystone taxa defined as highly connected taxa for T1C-1. The most abundant phyla in rhizosphere of different cultivars are shown in Fig. 2E. ANOVA analysis identified four phyla (Actinobacteria, Chloroflexi, Acidobacteria, and Gemmatimonadetes) that were significantly different between different cultivar communities. Actinobacteria and Chloroflexi abundances were highest in T1C-1 rhizospheres compared to soils of other cultivars, while Acidobacteria and Gemmatimonadetes abundances were lowest in T1C-1 rhizospheres (Fig. 2F).

Bt rice cultivation did not significantly alter the abundances of potential probiotic or phytopathogenic microorganisms

The effects of GM planting on potential probiotic or plant pathogen microbial taxa is an important aspect of environmental risk assessment for GM plants. Of the seven probiotic microorganisms that were differentially abundant in the rhizosphere of T1C-1 and Minghui 63, four exhibited enriched abundances and three exhibited decreased abundances in the T1C-1 rhizosphere compared to those of Minghui 63 (Fig. 3). Notably, the three taxa with decreased abundances belonged to the Lactobacillus genus that was also among the enriched microbial taxa. We speculate that the three Lactobacillus with decreased abundances did not decrease the beneficial impacts on plants demonstrated by other Lactobacillus. When considering paired communities of T1C-1 against Zhonghua 11 or T1C-1 against bulk soils, the proportion of increased or depleted probiotic microorganisms followed similar trends, indicating that T1C-1 plants tended to recruit rhizosphere bacteria that could aid in plant health and fitness. In addition, 21 microbial plant pathogens exhibited significantly different abundances in the rhizospheres of T1C-1 and Minghui 63 (Supplementary Fig. S6). Among all of the differentially abundant microorganisms, 10 exhibited increased abundances and 11 exhibited decreased abundances in the rhizosphere of T1C-1 plants relative to those of Minghui 63 plants. In particular, T1C-1 plants significantly inhibited the growth of the dominant phytopathogenic taxa Xanthomonas oryzae compared to Minghui 63 plants. Similar results were observed when comparing the communities of T1C-1 against Zhonghua 11 or T1C-1 against bulk soil communities, suggesting that T1C-1 planting alone altered the identity of the pathogenic microorganisms, but it did not trigger the accumulation of pathogenic microorganisms. Taken together, these results indicated that no adverse effects were apparent on soil microbial communities due to T1C-1 planting.

Associations between nitrogen-transforming microorganisms with variation in different N forms among Bt rice cultivation soils

Bt rice cultivation led to altered distributions of N forms in soils and thus, the effects of Bt rice on nitrogen-transforming microbial taxa were further investigated. A total of 41,543 non-redundant genes were identified that were involved in nitrogen-transforming reactions among the total 11,529,157 genes present in the soil samples. In addition to unclassified microorganisms, the most abundant microbial genera associated with nitrogen transformations were Pedosphaera, Gemmatimonas, Nitrospira, and Sorangium (Fig. 4A). PCoA analysis demonstrated that variation in the abundances of nitrogen-transforming taxa and their associated functions largely separated T1C-1 communities from Minghui 63 communities, consistent with the broad differences observed in soil profiles between these two cultivars (Fig. 4B, C). The relationships between soil properties and nitrogen-transforming microbial genera were also analyzed with RDA (Fig. 4D). The first two RDA components explained 64.14% of the total variation in genera-level composition. Soil pH was most correlated with variation in nitrogen-transforming genera composition, followed closely by NO₃^–-N, NO₂⁻-N and NH₄⁺-N (Supplementary Table S6). A paired sample t test approach was used to identify microbial genera that were significantly enriched or depleted in T1C-1 soils compared to Minghui 63 soils. The 15 most abundant nitrogen-transforming microbial genera with significant differences in abundance between T1C-1 and Minghui 63 soils included Candidatus Rokubacteria, unclassified Acidobacteria, Geobacter, and Anaeromyxobacter, among others (Fig. 4E). In addition, three nitrogen-transforming enzymes exhibited significantly different abundances between T1C-1 and Minghui 63 metagenomes including ferredoxin-nitrate reductase [EC:1.7.7.2], nitrite reductase (cytochrome c-552) [EC: 1.7.2.2], and carbamoyl-phosphate synthase (ammonia) [EC:6.3.4.16] (Fig. 4F).

The associations between soil microbial functions and the contents of different N forms were also evaluated. The levels of the aforementioned enzymes were lower in T1C-1 soils relative to Minghui 63 soils and these were particularly prevalent in the nitrogen metabolism pathway (ko 00910) (Fig. 4G). Considering the higher NH₄-N, NO₃-N, and NO₂-N contents in T1C-1 rhizospheres compared with Minghui 63 rhizospheres (Supplementary Table S3 and Fig. 4G), such relationships are consistent with different N forms acting as enzyme substrates. Thus, these results suggest that the cultivation of Bt rice altered soil nitrogen transformations by regulating the expression of associated genes. Similar results were also observed when comparing the paired samples of T1C-1 against Zhonghua 11 (Supplementary Fig. S7) or Zhonghua 11 against Minghui 63 (Supplementary Fig. S8).

Soil metabolomic profiles were altered by Bt rice cultivation, but not in surface waters

Changes in the metabolites that microorganisms may directly or indirectly interact with in soils were subsequently investigated using untargeted metabolomics. These analyses were used to assess the responses of microbial metabolic processes to Bt rice cultivation and were guided by the metagenomic insights into taxonomic and functional potentials. A total of 6527 and 6998 metabolites were detected in soils and surface waters, respectively. Among of these, a total of 527 and 539 metabolites were identified in soils and surface waters, respectively. PCA analysis of the metabolite profiles indicated that the first and second principal components significantly separated T1C-1 and Minghui 63 (or Zhonghua 11) metabolic profiles, suggesting that Bt rice significantly altered the low molecular weight metabolite profile in soils (Fig. 5A). Similar patterns were not observed in surface water metabolite profiles (Fig. 5A). A total of 988 metabolites were significantly up-regulated and 326 metabolites were significantly down-regulated in T1C-1 rhizospheres compared with Minghui63 rhizospheres (Fig. 5B), further suggesting that T1C-1 planting generally positively impacted soil metabolomes. Moreover, up-regulated metabolites were significantly higher than down-regulated metabolites in plant rhizospheres compared to bulk soils, indicating that plant cultivation increased soil metabolites. Discrimination with VIP values for OPLS-DA model variables led to the identification of 13 representative metabolites among the 200 most abundant metabolites (VIP > 1, p < 0.05) that differed in abundance between T1C-1 soils and others (Supplementary Fig. S9). The 13 representative metabolites comprised lipids, nucleotides, organoheterocyclic compounds, organic acids, and other compounds that enabled clear separation of T1C-1 profiles from others (Fig. 5C). Correlations between the 13 representative metabolites and soil physico-chemical properties were further analyzed (Fig. 5D). Soil Bt protein concentrations were positively correlated to the contents of deoxycytidine, N,N-dimethyl-Safingol and (S)-3,4-dihydroxybutyric acid, but negatively correlated to asteltoxin contents.

Associations between Bt rice-linked microorganisms and metabolites

A multi-omic framework integrating soil microbiomes and metabolomes was constructed to identify robust associations in the variation of differentially abundant microorganisms and metabolites in Bt and conventional rice variety soils. These robust associations contributed to characterize possible mechanistic relationships that occurred to Bt rice. To identify representative bacterial and eukaryotic taxa and metabolites most likely to explain differences among different soil samples, biomarker analysis was conducted with linear discriminant analysis (LDA) effect size (LEfSe) tests (Supplementary Figs. S10–S12). A total of 9017 and 6540 FDR-significant (q < 0.05) associations were identified between representative differentially abundant metabolites and bacterial or eukaryotic genera (Supplementary Fig. S13). A total of 24,000 bacterial and 10,750 eukaryotic associations between metabolites and genera were identified, with 37.6% and 60.8% being statistically significant, respectively. Although many metabolites were associated with one or more genera, they tended to not be mechanistically associated with most genera (and vice versa). In addition, a total of 978 associations with T1C-1 were identified, including 967 positive associations and 11 negative correlations that further contributed to understanding the putative mechanistic relationships that are altered when T1C-1 is grown compared to conventional rice.

To understand the functional consequences of microbial community changes due to T1C-1 planting, we first functionally profiled genes in the metagenomes and then summed their abundances according to Enzyme Commission (EC) annotations. The LEfSe approach was again used, but with abundance data of microbial gene encoding enzymes, revealing 27 representative enzymes that were differentially abundant (FDR-corrected q < 0.05) in T1C-1, Minghui 63, Zhonghua 11, or bulk soil metagenomes (Supplementary Fig. S14). A total of 1,661 FDR-significant (q < 0.05) associations (18.0% of all associations) were observed between representative differentially abundant enzymes and metabolites.

Bt rice re-establishes the carbon metabolism of soil microbial communities

A total of 105 metabolites annotated as standards were differential in T1C-1 soils compared with Minghui 63 soils, including 78 up-regulated and 27 down-regulated metabolites (VIP > 1 and p < 0.05) (Supplementary Fig. S15). The KEGG pathways associated with representative metabolites were primarily involved in lipid metabolism, carbohydrate metabolism, and energy metabolism, among other areas (Supplementary Fig. S16). Combined ‘omics’ approaches revealed that differentially enriched KEGG pathways in T1C-1 metagenomes compared to Minghui 63 metagenomes included galactose metabolism (ko00052), which is a component of carbohydrate metabolism (Fig. 6). In addition, the abundances of raffinose and melibiose were lower in T1C-1 soils compared to Minghui 63 soils based on non-targeted metabolomic profiling (Fig. 6A). Consistently, the abundances of microbial gene encoding enzymes involved in raffinose and melibiose-related metabolic pathways were significantly higher in T1C-1 metagenomes relative to Minghui 63 metagenomes, including alpha-glucosidase [EC: 3.2.1.20], alpha-galactosidase [EC: 3.2.1.22], and beta-fructofuranosidase [EC: 3.2.1.26] (Fig. 6B). These results suggest a re-establishment of carbon metabolism pathways in which raffinose and melibiose involvement might be attributed to the modulation of associated gene expression.

Putative mechanisms underlying the correlation of Bt rice-linked root exudate metabolites with rhizosphere microorganisms

A total of 6236 metabolites were detected in rice root exudates, of which 649 metabolites were structurally identified. PCA analysis indicated clear separation of the T1C-1 metabolite profiles from the Minghui 63 and Zhonghua 11 metabolomes, clearly suggesting altered metabolite profiles among different cultivars (Supplementary Fig. S17A). Significantly up-regulated metabolite fold-changes were higher than the significantly down-regulated metabolites in root exudates of T1C-1 profiles compared against Minghui63 or Zhonghua 11 profiles, further suggesting that T1C-1 planting contributed to enrichment of root exudates (Supplementary Fig. S17B–D). Thirty of the most abundant metabolites that were differentially enriched in T1C-1 profiles compared to Minghui 63 profiles were identified, including 27 up-regulated and three down-regulated metabolites (VIP > 1 and p < 0.05) (Supplementary Fig. S17E). A total of 81 identical metabolites were identified from root exudates and rhizosphere samples, among which 25 metabolites exhibited similar patterns of variation and highly statistically significant correlations (Spearman’s r > 0.7 and p < 0.05) between root exudate and rhizosphere profiles (Supplementary Table S7).

To determine how root exudates of different rice cultivars influence soil microbiomes, Spearman’s correlation analysis was conducted using compounds released as root exudates and the rhizosphere microorganisms that exhibited significant differences in abundances in soils of different plants. Acidobacteria abundances were significantly (p < 0.05) correlated with most root exudate compounds (194 correlations) while Chloroflexi were correlated with the least amount of compounds (27 correlations) (Supplementary Table S8). To identify whether a specific class of root exudate metabolite was associated with soil microbial communities, root exudates were classified as amino acids, carbohydrates, lipids, and organoheterocyclic compounds. Lipid metabolite levels were significantly correlated with the abundances of Acidobacteria, Actinobacteria, Chloroflexi, and Gemmatimonadetes (166 significant correlations; Supplementary Table S9), followed by organoheterocyclic compounds (54), amino acids (46), and carbohydrates (31). Numerous possible associations between T1C-1-representative root exudate metabolites and microbial taxa were observed (Supplementary Fig. S18, Fig. 7A). To experimentally validate the potential for T1C-1-representative metabolites to exert effects on the growth of T1C-1-representative microbial species, we cultured Terrimonas ferruginea that was dominant in soils, in the presence of four metabolites that it was associated with (Fig. 7B). Among the four predicted negative associations, three metabolites (except 20-hydroxy-leukotriene E4) indeed inhibited the growth of Terrimonas ferruginea at the tested concentrations. These results provide direct evidence for the influence of Bt rice-representative root exudates on Bt rice-representative microorganisms.

Source: Ecology - nature.com