Metabolomic and transcriptomic analyses of drought resistance mechanisms in sorghum varieties

Drought treatment and screening of drought-tolerant materials

Plump and uniformly sized sorghum seeds were selected for the filter paper germination method. They were first washed three times with sterilized distilled water, then soaked in 1% sodium hypochlorite solution for 20 min, and rinsed 3–5 times with sterile water. The seed germination test was conducted in a constant temperature and humidity chamber at 25 °C, 60% relative humidity, 12 h of light and darkness each, and continuous culture for 72 h. Each culture dish contained 40 sterilized sorghum seeds on filter paper, with the addition of 15 mL of 20% PEG solution by mass. An equivalent amount of distilled water served as a control, and each treatment was replicated five times. Photos were taken at 12, 24, 36, 48, 60, and 72 h to measure shoot and root lengths. Samples collected at 72 h were promptly frozen in liquid nitrogen and stored at −80 °C for subsequent transcriptome and metabolome sequencing.

Transcriptome sequencing

Total RNA from the plant was extracted following the manufacturer’s protocol for the RNAprep Pure Plant Kit (Tiangen, Beijing, China). The RNA concentration, purity, and integrity were assessed with a NanoDrop 2000 and an Agilent Bioanalyzer 2100 system to verify the quality of the RNA prior to further procedures. One microgram of high-quality RNA served as the starting material for sample preparation, which we carried out using a Hieff NGS Ultima Kit for mRNA library construction. The mRNA was purified, cDNA was synthesized, and the sequencing library was prepared following the manufacturers’ recommendations, including ligation of adaptors and PCR to produce the final library, all of which were validated using a bioanalyzer. Sequencing was conducted on an Illumina NovaSeq platform, creating an extensive set of 150 bp paired-end reads. The Fastp software (version 0.23.4) was utilized to eliminate adapter sequences, filter out low-quality and N sequences exceeding a 5% ratio, resulting in clean reads for subsequent analysis (Chen et al., 2018). Subsequently, HISAT2 was employed to align the clean reads to the reference genome (Sorghum_bicolor.v3.1.1, https://phytozome-next.jgi.doe.gov/) (Kim et al., 2019). Gene function annotations were enriched and sourced from various databases, and the quantification of expression levels was accurately conducted using the fragments per kilobase of exon model per million mapped fragments (FPKM) metric.

RNA-seq differential expression analysis

DESeq2 software (version: 1.46.0) was employed to identify differentially expressed genes (DEGs) based on gene count values in each sample, using |fold change| ≥ 2 and FDR < 0.01 as screening criteria (Love, Huber & Anders, 2014). Subsequently, for all DEGs, the ClusterProfiler package (version: 4.0) in R was utilized to conduct GO and KEGG enrichment analyses via the hypergeometric test method (Wu et al., 2021). Transcription factors (TFs) were predicted and annotated utilizing the Plant Transcription Factor Database (Plant TFDB, https://planttfdb.gao-lab.org/).

Metabolite extraction and detection

The samples were freeze-dried under vacuum using a SCIENTZ-100F lyophilization dryer (LCD Display Freeze Dryer; SCIENTZ). Following this, the dried samples were ground into a powder using an MM 400 grinding machine (30 Hz, 1.5 min; Retsch). The powdered samples were dissolved in a 70% methanol extract (1.2 mL), refrigerated overnight at 4 °C, with gentle swirling performed six times during the incubation to enhance the extraction efficiency. Subsequently, each sample was centrifuged at 12,000 revolutions per minute for 10 min, and the upper liquid fraction was filtered through a 0.22-mm membrane filter (0.22-mm pore size).

Ultra performance liquid chromatography (UPLC) conditions: The Waters Xevo G2-XS QTOF high-resolution mass spectrometer was utilized to gather primary and secondary mass spectrometry data in MSe mode using MassLynx V4.2 acquisition software (Waters). In each data acquisition cycle, dual-channel data acquisition was conducted with low and high collision energy settings. The low collision energy was maintained at 2 V, while the high collision energy ranged from 10 to 40 V. The scanning frequency for a mass spectrum was set at 0.2 s. The electrospray ionization (ESI) ion source parameters were configured as follows: capillary voltage at 2,000 V (positive ion mode) or −1,500 V (negative ion mode), cone voltage at 30 V, ion source temperature at 150 °C, desolvent gas temperature at 500 °C, backflush gas flow rate at 50 L/h, and desolventizing gas flow rate at 800 L/h.

Metabolomics data analysis

The filtered samples were stored in containers for further analysis using ultra-performance liquid chromatography and tandem mass spectrometry. Metabolome profiling was carried out employing a broadly targeted metabolomics approach utilizing the MWDB database developed by Wuhan Mettware Biotechnology (http://www.metware.cn/). The MWDB-assigned metabolites were quantified utilizing triple-level quadrupole mass spectra acquired in multiple-reaction monitoring mode. The MWDB is a database for metabolomics analysis that can help researchers quickly and accurately analyze metabolites in different plant species; in the metabolomics analysis of sorghum, the MWDB can provide rich metabolite data. Prior to data analysis, quality control checks were carried out to verify data reliability. Principal component analysis (PCA) was utilized to investigate variabilities both between and within groups. Differential accumulated metabolites (DAMs) were assessed using orthogonal partial least-squares discriminant analysis (OPLS-DA) as described by Kang et al. (2022). Metabolites with a variable importance in projection (VIP) score of ≥1 and a fold change (FC) of ≥2 (or ≤0.5) were identified and classified as DAMs.

Phenotypic evaluation of GL98 and GL220 under drought stress

To evaluate the drought resistance of sorghum varieties GL98 and GL220, we used 20% PEG to simulate drought and detected the phenotypes under stress conditions (Fig. 1A). Without PEG treatment, the roots and shoots of GL98 and GL220 were able to grow normally. The shoot lengths of GL98 and GL220 were similar, but the root length of GL98 was longer than that of GL220 (Fig. 1B). After drought stress, the lengths of the roots and shoots of GL98 and GL220 decreased, but as the treatment time increased, the growth rates of the roots and shoots of GL220 exceeded those of GL98 at 72 h (Fig. 1C). These findings indicate that GL220 has better drought resistance than does GL98 and grows more vigorously under drought stress conditions. To analyze the mechanism of the difference between GL220 and GL98 under drought stress conditions, three replicate transcriptome and five replicate metabolomics sequencing experiments were performed on roots treated with 20% PEG for 72 h and roots grown normally for 72 h.

RNA-seq analysis

Transcriptome sequencing analysis was performed on 12 samples of roots of GL220 and GL98 under normal growth conditions and 20% PEG treatment for 72 h. After the raw data were filtered, 75.05 Gb of clean data were obtained. The proportion of Q30 bases exceeded 94.52%, and the guanine (G) and cytosine (C) content exceeded 53.73% (Table S1). First, the Pearson correlation coefficient between samples was calculated based on the fragments per kilobase of exon model per million mapped fragments (FPKM) value of each gene, and the correlation coefficient of the same sample exceeded 0.88 (Fig. 2A). Principal component analysis (PCA) clustered the replicates of the same sample together, whereas the samples between the treatment and the control were separated (Fig. 2B). These findings indicate that the RNA-seq data were reliable and reproducible and were suitable for further analysis.

RNA-seq differential expression analysis

Through differential expression analysis, we identified a total of 6,344 differentially expressed genes (DEGs) (Figs. 3A and 3B). A total of 4,384 DEGs were identified between GL98 under drought stress and the control at 72 h, of which 2,682 were upregulated and 1,702 were downregulated, including 2,220 unique DEGs (Figs. 3A and 3B). A total of 3,453 DEGs were identified between the GL220 drought stress treatment group and the control group at 72 h, of which 2,435 were upregulated and 1,018 were downregulated, including 1,306 unique DEGs. A total of 724 DEGs were identified between the GL98 and GL220 controls, 319 of which were upregulated and 405 were downregulated, including 256 unique DEGs. A total of 343 DEGs were identified in GL98 and GL220 under drought stress for 72 h, of which 191 were upregulated and 152 were downregulated, including 86 unique DEGs.

To further determine the functions of the 6,344 DEGs, these DEGs were annotated via GO and KEGG. GO enrichment analysis revealed significant annotations related to secondary metabolite biosynthetic processes, the regulation of vesicle-mediated transport, lignin metabolic processes, lignin biosynthetic processes, phenylpropanoid biosynthetic processes, carbohydrate metabolic processes, flavonoid metabolic processes, starch metabolic processes, polysaccharide catabolic processes and the regulation of hormone levels (Fig. 3C). KEGG enrichment analysis revealed significant annotations related to galactose metabolism, transporters, carbohydrate metabolism, phenylpropanoid biosynthesis, thiamine metabolism, starch and sucrose metabolism, amino acid metabolism, fructose and mannose metabolism, flavonoid biosynthesis, cytochrome P450, carbon fixation in photosynthetic organisms and phenylalanine metabolism pathways (Fig. 3D).

DEG clustering analysis

Genes with similar expression patterns usually have similar functions. The K-means clustering algorithm was used to cluster all the DEGs, and ultimately, four clusters were obtained (Fig. 4). The expression of Cluster 1 genes increased after drought stress in GL98 and slightly increased between the GL220 treatment and the control. A total of 437 DEGs were significantly enriched in the membrane trafficking and transporter pathways. The expression of Cluster 2 genes increased after drought stress in both GL98 and GL220. The analysis revealed 3,810 DEGs that were significantly enriched in the phenylpropanoid biosynthesis and flavonoid biosynthesis pathways. The expression of Cluster 3 genes decreased after drought stress in GL98 and between the GL220 treatment and the control. There were 836 DEGs that were significantly enriched in the starch and sucrose metabolism and carbohydrate metabolism pathways. The expression of Cluster 4 genes decreased after drought stress in both GL98 and GL220. The analysis revealed 1,261 DEGs that were significantly enriched in the thiamine metabolism and fructose and mannose metabolism pathways.

TF analysis

A total of 309 differentially expressed TFs were identified among the 6,344 DEGs, with AP2/ERF and MYB accounting for the greatest proportions, with 28 and 27, respectively (Fig. 5A). Hierarchical clustering was used to identify two clusters for the 309 differentially expressed TFs (Fig. 5B). Cluster 1 contained 159 TFs, and their expression levels increased in both GL98 and GL220 after drought stress. Four top TFs were identified, namely, Sobic.003G017601 (bHLH), Sobic.002G225100 (bZIP), Sobic.002G189000 (HD-ZIP) and Sobic.004G101400 (HSF). Cluster 2 contained 150 TFs, and their expression levels decreased in both GL98 and GL220 after drought stress. Four top TFs were identified, namely, Sobic.006G080500 (bHLH), Sobic.001G526800 (MADS), Sobic.003G078100 (AP2/ERF) and Sobic.003G102200 (bHLH).

Global analysis of metabolomics data

A total of 20 samples from the four groups were detected and analyzed by UPLC‒MS, which revealed 3,913 metabolites (Table S2). For example, there were 138 flavonoids, 47 steroids and 25 steroids related to plant stress resistance. These metabolites play key roles in the development of plant stress resistance. PCA revealed that the five biological replicates of the same treatment were clustered together, indicating that the metabolite extraction and detection techniques had good reproducibility and high reliability (Fig. 6A). After treatment, drought stress had an impact on sorghum metabolites, indicating that these metabolites changed dynamically during the response of sorghum to drought stress. The 3,913 metabolites identified were classified and annotated, and the results revealed that the metabolites were mainly divided into 10 categories: lipids and lipid-like molecules accounted for 31.20%; organic acids and derivatives accounted for 17.55%; phenylpropanoids and polyketides accounted for 12.01%; organoheterocyclic compounds accounted for 11.68%; organic oxygen compounds accounted for 10.33%; benzenoids accounted for 7.77%; nucleosides, nucleotides, and analogs accounted for 5.22%; alkaloids and derivatives accounted for 2.17%; organic nitrogen compounds accounted for 1.58%; and lignans, neolignans and related compounds accounted for 0.49% (Fig. 6B).

DAM analysis

DAMs between samples were identified by differential expression analysis, which revealed a total of 1,942 DAMs, including five common DAMs (Figs. 7A and 7B). A total of 1,363 DAMs were identified in the GL98 drought stress and control groups at 72 h, of which 347 were upregulated and 1,016 were downregulated, including 342 unique DAMs (Figs. 7A and 7B). A total of 1,411 DAMs were identified in the GL220 drought stress and control groups after 72 h, of which 389 were upregulated and 1,022 were downregulated, including 370 unique DAMs. A total of 152 DAMs were identified between GL98 and GL220 controls, of which 68 were upregulated and 84 were downregulated, including 49 unique DAMs. A total of 232 DAMs were identified in GL98 and GL220 under drought stress for 72 h, of which 87 were upregulated and 146 were downregulated, including 71 unique DAMs.

KEGG enrichment analysis revealed that 1,942 DAMs were significantly enriched in linoleic acid metabolism, steroid biosynthesis, alpha-linolenic acid metabolism, anthocyanin biosynthesis, carotenoid biosynthesis, starch and sucrose metabolism, biosynthesis of unsaturated fatty acids, biosynthesis of various plant secondary metabolites, fructose and mannose metabolism, flavonoid biosynthesis, phenylpropanoid biosynthesis, and arginine and proline metabolism pathways (Fig. 8A). All the DAMs were clustered via the K-means clustering algorithm, and four clusters were ultimately obtained (Fig. 8B). The expression of Cluster 1, which included 1,030 DAMs, decreased after drought stress in both GL98 and GL220. The expression of Cluster 2, which included 338 DAMs, decreased after drought stress in GL98, with a smaller decrease in expression detected between the GL220 treatment and the control. In both GL98 and GL220, Cluster 3, which included 426 DAMs, presented increased expression after drought stress. In GL220, Cluster 2, which included 148 DAMs, presented increased expression after drought stress and a smaller increase in expression between the GL98 treatment and the control.

Combined transcriptomics and metabolomics analyses of the drought tolerance mechanism of sorghum

Analysis of the DAMs indicated significant enrichment in pathways associated with metabolites such as phenylpropanoid biosynthesis, flavonoid biosynthesis, and anthocyanin biosynthesis post-treatment, highlighting the crucial roles of flavonoids in enhancing plant drought resistance. First, the expression patterns of 20 DEGs related to flavonoid biosynthesis were examined. Following drought stress, the majority of DEGs associated with flavonoid biosynthesis exhibited increased expression levels, with GL98 showing higher expression levels compared to GL220 (Fig. 9A). Thirteen different flavonoid metabolites were detected, but the contents of most of them decreased after drought stress. Only the contents of 5-O-caffeoylshikimic acid, pinostrobin, luteoforol, bracteatin 6-O-glucoside, and p-coumaroyl quinic acid increased after drought stress (Fig. 9B). Among them, the contents of liquorol, bracteatin 6-O-glucoside and p-coumaroyl quinic acid increased in both GL98 and GL220 after drought stress. After drought stress, the content of 5-O-caffeoylshikimic acid increased in GL98 and decreased in GL220. The pinostrobin content remained unchanged in GL98 but increased in GL220. The Pearson correlation coefficients between the DEGs of the flavonoid pathway and DAMs were computed. Genes and metabolites with absolute correlation coefficients exceeding 0.9 and p values below 0.05 were selected for visualization (Fig. 9C). A total of 17 flavonoid biosynthesis genes were significantly correlated with 10 flavonoid metabolites, of which 17 genes were significantly negatively correlated with seven flavonoid metabolites and 15 genes were significantly positively correlated with three flavonoid metabolites. Among them, Sobic.007G058600 was significantly correlated with the metabolites of the 10 most abundant flavonoid pathways, and Sobic.010G052200 was significantly negatively correlated with only Apiforol.