Integrated transcriptomic and metabolomic analyses provide new insights into alkaline stress tolerance in Gossypium hirsutum

Abstract

Introduction: Cotton, one of the most important economic crops worldwide, has long been bred mainly for improvements in yield and quality, with relatively little focus on salt-alkali resistance.

Methods: In this study, transcriptomic and metabolomic sequencing were performed on Gossypium hirsutum exposed to alkaline stress for different durations.

Results: The results of sample clustering, principal component analysis (PCA), and the number of differentially expressed genes (DEGs) revealed that 12 hours and 24 hours were the periods during which upland cotton presented the strongest response to salt stress, with flavonoid biosynthesis and alpha-linolenic acid metabolism playing significant roles during this time. A total of 6,610 DEGs were identified via comparison to the 0 h time point, including 579 transcription factors (TFs) that were significantly enriched in pathways such as flavonoid biosynthesis, the cell cycle, the cytochrome P450 pathway, phenylalanine metabolism, phototransduction, and alpha-linolenic acid metabolism. Through ultrahigh-performance liquid chromatography-MS (UPLC-MS), 4,225 metabolites were identified, and 1,684 differentially accumulated metabolites (DAMs) were identified by comparison to the levels at 0 h. A joint analysis of RNA-seq and metabolomic data revealed that the flavonoid biosynthesis and alpha-linolenic acid metabolism pathways play key roles in the response of G. hirsutum to alkaline stress, and the key genes in these pathways were identified. The weighted gene correlation network analysis (WGCNA) revealed 15 candidate genes associated with alkali tolerance in cotton, including 4 TFs and 4 genes related to flavonoid and anthocyanin biosynthesis.

Conclusion: In conclusion, our study provides a theoretical foundation for understanding the molecular mechanisms underlying alkali tolerance in cotton and offers new gene resources for future research.

Keywords: Gossypium hirsutum; RNA-seq; alkaline stress; candidate genes; metabolome.

Figure 1

Correlation analysis and PCA of the RNA-seq data from 18 (G) hirsutum samples under alkaline stress. (a) Correlation analysis of the samples; a value from 0.4 to 1 represents the magnitude of the correlation coefficient between samples, with 0.4 indicating the lowest correlation coefficient and 1 indicating the highest correlation coefficient between samples. (b) PCA of the samples; each point represents a sample, with different colors used to identify processing at different times.

Figure 2

Numbers of DEGs and enrichment analysis results compared with those at 0 h of alkaline stress. (a) Numbers of upregulated and downregulated DEGs identified at each time point compared with 0 h of stress. (b) Venn diagram of unique and common DEGs identified at each time point compared with 0 h of alkaline stress. (c) KEGG enrichment analysis of DEGs identified at various time points compared with 0 h of alkaline stress. (d) GO enrichment analysis of all DEGs compared with the genes detected at 0 h of alkaline stress.

Figure 3

Clustering and enrichment analyses of all DEGs. (a) Line plot of the expression patterns of all the DEGs based on the clustering analysis. The red numbers represent the quantity of DEGs and differentially expressed TFs identified in each cluster. (b) KEGG pathway enrichment analysis of each category of DEGs. The intensity of the color represents the magnitude of the p value, specifically the value of -log₁₀ (p value). A darker color indicates a smaller p value, which corresponds to a larger value of -log₁₀ (p value), whereas a lighter color indicates a larger p value, equating to a smaller value of -log₁₀ (p value).

Figure 4

DAMs and enrichment compared with the metabolites detected at 0 h of stress. (a) Percentages of classified metabolites, (b) numbers of upregulated and downregulated DAMs at each time point compared with that at 0 h of stress, (c) Venn diagram showing unique and shared DAMs at each time point compared with those at 0 h of stress, and (d) KEGG enrichment analysis of all DAMs compared with the metabolites detected at 0 h of stress.

Figure 5

Clustering and classification of all the DAMs. (a) Line plot showing the pattern of changes in the contents of all the DAMs based on the clustering analysis; the red numbers represent the quantity of DAMs identified in each cluster. (b) Classification of all the DAMs into different categories.

Figure 6

Results of the combined RNA-seq and metabolomic analyses, analysis of changes in DAMs and DEGs in the flavonoid biosynthesis pathway, and correlation analysis. (a) KEGG pathway annotations of DAMs and DEGs. (b) Changes in the expression of DEGs in the flavonoid biosynthesis pathway were quantified via standardized scoring, with values standardized to range from -2 to 2. (c) Heatmap of changes in the levels of DAMs in the flavonoid biosynthesis pathway. Standardized scoring was used, with values standardized to range from -2 to 2. (d) Correlation network of flavonoid biosynthesis pathway metabolites and genes; the red line represents a significant positive correlation, and the green line represents a significant negative correlation.

Figure 7

Changes in DAMs and DEGs in the JA biosynthesis pathway and correlation analysis. (a) Changes in the expression of DEGs in the JA biosynthesis pathway were quantified via standardized scoring, with values standardized to range from -2 to 2. (b) Heatmap of changes in the levels of DAMs in the JA biosynthesis pathway; standardized scoring was used, with values standardized to range from -2 to 2. (c) Correlation network between metabolites and genes in the JA biosynthesis pathway; the red line represents a significant positive correlation, and the green line represents a significant negative correlation.

Figure 8

Proportions and expression patterns of differentially expressed TFs. (a) Proportional area chart of differentially expressed TFs. (b) Heatmap of the differential expression patterns of TFs, with the TFs exhibiting the greatest fold changes shown on the right.

Figure 9

WGCNA and analysis of the candidate gene expression patterns. (a) WGCNA clustering dendrogram, with different colors representing different modules. (b) Correlation and significance analyses between modules and alkaline stress in cotton at different time points. The thickness of the lines represents the magnitude of the correlation coefficient between the modules and periods, with gray indicating a p value greater than or equal to 0.05 and green indicating a p value less than 0.05. (c) Gene network diagrams for the pink, green, brown, magenta, and blue modules. (d) Analysis of the expression patterns of the 15 candidate genes in cotton under alkaline stress. The error bars represent the average values ± SDs from three replicates (*P<0.05 and **P<0.01).