Integrative analysis of the microbiome and metabolome in understanding the causes of sugarcane bitterness

Soil microbial diversity and community structure

In total, 1,334,381 reads were obtained for the bacterial 16S rRNA genes by high-throughput sequencing. After screening these gene sequences with strict criteria (described in “Materials and methods”), 1,061,916 valid sequences were obtained, accounting for 79.6% of the raw reads. Figure 1A shows that the observed richness, Chao1, and Shannon index in the SS (sweet sugarcane) group supported significantly higher richness (P < 0.01) and diversity (P < 0.05) than in the BS (bitter sugarcane) group. It suggested that BS and SS differed regarding several metrics considered in alpha-diversity. Beta-diversity analysis of the two groups showed that the microbial community structure of BS was significantly different from that of SS (P < 0.05), i.e., a distinct separation was observed (Fig. 1B, Fig. S1).

A phylum-level analysis demonstrated the compositional differences of soil microbial communities, showing that all samples were dominated by Proteobacteria, Acidobacteria, Bacteroidetes, Actinobacteria, Firmicutes, Chloroflexi, Gemmatimonadetes, Latescibacteria, and Verrucomicrobia (Fig. 1C). With respect to the top 10 phyla, the soil bacterial community of the BS group had greater abundances of Bacteroidetes, Firmicutes, and Chloroflexi (P < 0.05), compared with lower abundances of Actinobacteria, Latescibacteria, and Verrucomicrobia relative to the SS group (P < 0.05). In the analysis of the microbiota composition at the family level, Rhodanobacteraceae, Xanthomonadaceae, and Ktedonobacteraceae were highly enriched in BS compared to SS (P < 0.05), while Nitrosomonadaceae and Gemmatimonadaceae were highly enriched in SS versus BS (P < 0.05) (Fig. 1D). In order to distinguish the predominant taxa, LEfSe [linear discriminant analysis (LDA) integrated with effect size] was used to generate a cladogram to identify the specific bacteria associated with the BS group (Fig. 2A,B). We identified 81 distinctive genera from the soil bacterial community of the two groups (Table S1). Some bacteria, including Gammaproteobacteria, Acidobacteriia, Ktedonobacteraceae, Muribaculaceae, and Stenotrophomonas, were all significantly overrepresented (all LDA scores (log₁₀) > 4) in the bitter sugarcane group of BS, whereas Deltaproteobacteria, Verrucomicrobiae, Thermoleophilia, and Nitrosomonadaceae were the most abundant microbiota in the control group (SS; LDA scores (log₁₀) > 4). Among them, the top 18 distinct genera (P < 0.05, relative abundance ≥ 0.01%) were selected as key discriminants (Fig. 2C, Table S2). In particular, the genus Chujaibacter in Rhodanobacteraceae and Stenotrophomonas in Xanthomonadaceae were highly abundant in BS versus SS (P < 0.05).

At the OTU level, Venn diagrams were constructed to better represent the results of shared richness between microbial samples from two sites, which, according to a detailed analysis, represented the number of OTUs shared between the two groups (Fig. S2). The two soil groups had a high number of shared sequences (3173 OTUs, 63.1%), revealing the presence of a strong core microbiota composed of Proteobacteria (27.4%), Actinobacteria (9.9%), Acidobacteria (9.7%), and Bacteroidetes (8.9%) (Fig. S2, Table S3). Bacterial functional profiles were inferred from 16S rRNA gene data using PICRUSt. Differences in the microbial function of the soil microbiome between the two groups were compared (Fig. S3). Compared to the SS group, the BS group exhibited lower metabolism of carbon fixation, citrate cycle, and chlorophyll synthesis (P < 0.05).

In addition, the relative abundance of soil nitrogen-fixing bacteria in the two groups was compared. Figure S4 shows that Nitrospira japonica was significantly enriched in SS compared with BS (P < 0.05), indicating that the soil in SS could have a higher efficiency of nitrogen conversion and metabolism versus BS. The bacterial nitrogen-fixing-associated gene NifH was further detected and showed a greater abundance in SS compared to BS (Fig. S4).

Metabolomics analysis of the sugarcane stalks

A multivariate analysis method, PLS-DA, was used to identify the key compounds (detected by LC–MS) responsible for the differentiation. Score plots based on the first two components showed that the two groups were well separated (Fig. 3A). In total, we detected 1245 metabolites in stalk tissues based on LC–MS analyses. Among these metabolites, 247 (19.8%) were differentially abundant between sugarcane stalks collected from site A (BS) and from site C (SS) (Fig. 3B, P < 0.05; Table S4). These metabolites, including amino acids, organic acids, sugars, and sugar alcohols, fatty acids, lipids, benzenoids, phenylpropanoids, and polyketides, are involved in multiple biochemical processes in sugarcane stalks.

To identify the compounds responsible for the difference, the VIP combined with the univariate statistical results (Table S4) was used to select variables with the most significant contribution to the discrimination of the sugarcane compounds of the two groups. Statistically significant differences were observed between the BS and SS groups shown in the volcano map (Fig. 3C). Compared to SS, 112 metabolites were significantly increased in BS, and 135 metabolites were significantly decreased (VIP > 1, P < 0.05). Among them, lamioside (62.94-fold increase), primisulfuron (57.93-fold increase), 3,5-dimethylpyrazin-2-ol (40.66-fold increase), sitaxentan (33.53-fold increase), crolibulin (26.50-fold increase), 2-formylpyridine (7.09-fold increase), (z)-desulfoglucotropeolin (6.45-fold increase), and asparagine (3.23-fold increase) were produced in greater concentrations in BS, while leucoside (67.59-fold decrease), isovitexin 2″-O-arabinoside (29.84-fold decrease), vaccarin (25.65-fold decrease), kuromanin (21.22-fold decrease), and robinin (20.61-fold decrease) were greatly reduced in BS relative to SS (VIP > 3, P < 0.05).

Through KEGG pathway annotation, 194 metabolites of these 247 distinct metabolites were mapped onto 13 different KEGG metabolic pathways, including amino acid metabolism, metabolism of cofactors and vitamins, metabolism of terpenoids and polyketides, and biosynthesis of other secondary metabolites. However, the KEGG-enrichment scatterplot showed that there were only two distinguishing pathways including flavonoid biosynthesis, flavone and flavonol biosynthesis which were significantly enriched in BS compared to SS (P < 0.05; Fig. 3D). Therefore, 22 distinct metabolites were selected from these two KEGG pathways for correlation analysis, including vitexin, naringin, apigenin, hesperetin, rhoifolin, and chlorogenic acid (Table 1).

Correlation between soil properties, the soil microbiome, and the sugarcane metabolome

Analysis showed that the soil from the BS group only had a significantly lower available nitrogen content and significantly higher potassium ion content compared to soil from the SS group (P < 0.05; Table 2). The contents of total carbon, PO₄-P, and other elements did not show dramatic changes. Canonical correlation analysis (CCA) based on the genus-specific abundances (Table S2) from each soil sample was used to determine how the various soil microbial communities correlated with soil edaphic factors. The bacterial community of BS was positively correlated with P and K, but negatively correlated with C, N, pH, and Na⁺ (Fig. 4A). On the contrary, the correlation result was opposite for the bacterial community of SS. Among the 18 distinct bacterial genera associated with soil properties in Fig. 4A, we found that nitrogen had a significantly positive correlation (P < 0.05) with Reyranella (r > 0.99), Dongia (r > 0.94), Gaiella (r > 0.92), Candidatus_Udaeobacter (r > 0.91), and unidentified_Acidobacteria (r > 0.88), and had a significantly negative correlation (P < 0.05) with Ralstonia (r < − 0.83). On the contrary, potassium was found negatively correlated (P < 0.05) with Candidatus_Udaeobacter (r < − 0.82) and Reyranella (r < − 0.81), and positively correlated (P < 0.05) with Ralstonia (r > 0.84) (Table S5).

To determine the impact of soil chemicals on the sugarcane metabolome, CCA was performed and heatmaps (Fig. 4B, Fig. S5) constructed using the Pearson’s correlation coefficient to determine the correlations between the distinct sugarcane metabolites and soil edaphic factors. The correlation analysis revealed that only certain soil elements, including nitrogen and potassium in particular, had significant correlations with the distinct metabolites (Fig. 4B). Nitrogen was significantly positively associated (P < 0.05) with hesperetin (r > 0.86), luteolin (r > 0.82), and neodiosmin (r > 0.84), and significantly negatively correlated (P < 0.05) with dl-arginine (r < − 0.94), 7-hydroxycoumarine (r < − 0.81) and l-pyroglutamic acid (r < − 0.81). Potassium was positively correlated (P < 0.05) with chlorogenic acid (r > 0.86) and dl-arginine (r > 0.83) (Table S6).

Soil microbiome strongly linked to the sugarcane metabolites

To explore the functional correlation between the soil microbiome changes and metabolite perturbations, a Pearson’s correlation matrix was generated by calculating the Pearson’s correlation coefficient between the different bacterial abundances (18 genera) and distinct metabolites (22 compounds mainly involved in the flavonoid biosynthesis pathway, VIP > 1.0) of the BS and SS groups (Table S7). Significant associations were selected (P < 0.05), and their biological networks were visualized using Cytoscape v3.7.1, based upon which clear correlations were identified between the perturbed soil microbiome and altered metabolite profiles (Fig. 5). As expected, the composition of the soil microbiome was widely associated with the sugarcane metabolic content. At a 5% false discovery rate (FDR), we found 108 associations between sugarcane metabolites and microbial species (Fig. 5A). In particular, 88.9% of the species were associated with 86.4% of the sugarcane metabolites. We observed 50.9% positive and 49.1% negative associations with microbial species and metabolites, respectively. On average, each metabolite was associated with five species. Six microbial species played a major metabolic role and were independently associated with 90.9% of the sugarcane metabolites: unidentified_Acidobacteria (12 metabolites), unidentified_Nitrospiraceae (11 metabolites), unidentified_ Acidimicrobiia (11 metabolites), Dongia (10 metabolites), Galiella (9 metabolites), Reyranella (8 metabolites), and Anaeromyxobacter (8 metabolites).

Among the distinct metabolites, apigenin, biochanin A, hesperidin, luteolin, and trigonelline had ten correlations on average with different microbial species. Examining these metabolites, we noted that apigenin had the most associations including five positive correlations and eight negative correlations. The down-regulation of apigenin (0.28-fold) was highly positively correlated (P < 0.05) with the accumulation of unidentified_Nitrospiraceae (r > 0.95), unidentified_Acidobacteria (r > 0.94), Anaeromyxobacter (r > 0.87), Dongia (r > 0.84), and Gaiella (r > 0.83) and was negatively correlated with the accumulation of Acidibacter (r < − 0.93), Conexibacter (r < − 0.92), and the remaining six genera (Fig. 5B). Chlorogenic acid (15.26-fold change, the most up-regulated one among the 22 distinct metabolites) was negatively associated (P < 0.05) only with Dongia (r < − 0.81). Kuromanin (0.05-fold change, the most down-regulated metabolite) was positively associated (P < 0.05) with unidentified_Nitrospiraceae (r > 0.91) and unidentified_Acidobacteria (r > 0.90) and was negatively associated (P < 0.05) with Conexibacter (r < − 0.92) and Acidibacter (r < − 0.86). Thus, the soil microbiome strongly linked to the sugarcane metabolites may contribute to the different flavors of the two groups.

Source: Ecology - nature.com