Integrative metabolomics and transcriptomics

Full text

Turn on search term navigation

Introduction

Maidong [Ophiopogon japonicus (L.f.) Ker-Gawl], a perennial herbaceous plant in the Liliaceae family [1]. This plant serves as both a medicinal and food source [2–5], primarily composed of flavonoids [6,7], steroidal saponins [8,9], polysaccharides [10,11]and terpenoids [12]. Maidong is primarily found in Sichuan and Zhejiang [13], of which Sichuan Maidong (abbreviated as Chuanmaidong, CMD), which has the advantages of shorter planting cycle, higher yield and cost-effectiveness, is the main variety of Maidong [14]. Currently, although an increasing number of studies have explored the pharmacological and nutritional attributes of CMD, such as alleviation of diseases such as colitis and lung cancer [15,16], and as functional food ingredients [17], there is limited in-depth research on its agronomic properties, in particular on the mechanisms of its responses to adapt to abiotic stresses.

Climate change has exacerbated global water scarcity, making drought stress one of the most detrimental abiotic factors impacting crop production and food security, which is primarily due to agriculture’s heavy reliance on water resources [18,19]. Drought stress disrupts normal plant growth, disturbs metabolic balance, reduces water use efficiency, impairs plant productivity, and poses a significant threat to sustainable crop production [20]. Many agricultural areas are affected by drought, resulting in crop yield losses that can reach up to 50% or even higher [21]. Consequently, uncovering distinctive mechanisms of drought tolerance in plants and screening high quality varieties that can withstand drought stress are urgent issues.

To investigate the mechanisms of CMD adaptation to drought and the ability of different CMD species to adapt to drought, we integrated transcriptomic and metabolomic analyses to identify differentially expressed genes (DEGs), differentially expressed metabolites (DEMs) and differentially expressed transcription factors (DETFs) that adapt CMD to different levels of drought stress, as well as the differences in response to drought stress in different cultivars. Like many other plants, CMD undergoes a complex adaptive process when exposed to drought stress, which may significantly affect the transformation of its secondary metabolites [22,23]. Thus, we also explored the flavonoid regulatory pathways in the drought-adaptation process of CMD. Overall, this study provides a theoretical basis for the selection of optimal drought-resistant CMD varieties and a scientific reference for elucidating the drought tolerance mechanism of CMD.

Materials and methods

Plant material and growing conditions

Erect Plants (EP) and Creeping Plants (CP) of CMD varieties were obtained from Sichuan Province, China (31°30’N, 104°95’E and 31°32’N, 104°89’E). These specimens were taxonomically confirmed as Ophiopogon japonicus (L. f.) Ker-Gawl. of the Liliaceae family by researcher Dr. Tiezhu Chen. The seedlings were planted in plastic pots and cultivated in the greenhouses at the Sichuan Academy of Chinese Medicine Sciences (30°37’N, 104°40’E). The pots had a caliber of 24 cm, a base diameter of 19.5 cm, and a height of 26.5 cm. The substrate was formulated with an organic matter content of 1.80 g/kg, total N 2.50 g/kg, total P 1.50 g/kg, pH 6.8, and a maximum water holding capacity of the substrate of 25.87%. During the growth and development of the two CMD germplasm resources, the amount of watering was controlled, 16 h/day of light was given and other cultivation and management modes were kept consistent. Seedlings were acclimatized for 28 days, and then watering was controlled. Then, the seedlings of similar height and growth potential were randomly divided into control (CK), moderate drought (D1) and severe drought (D2) groups, with six pots in each group. In the CK group, standard watering practices were maintained, delivering 1598 mL per week, one seventh a day. The standard weekly watering amount is an average of the period of maximum water need (April-June) for CMD, and watering is done once a day. Meanwhile, the D1 group received reduced watering of 799 mL per week, and the D2 group was subjected to even more stringent watering, receiving 559.3 mL per week. After 21 days of drought stress experiment, the plants were dug up and the roots were washed with pure water, then the fresh tubers were quickly cut and placed in liquid nitrogen in preparation for the next experiment. The three samples from each group (CK, D1, D2) of both EP and CP were randomly selected for comprehensive analysis through transcriptomic and metabolomic profiling.

RNA preparation and RNA-seq

Total RNA was extracted from EP and CP fresh tubers exposed to different drought stress levels, and then RNA libraries were constructed [24,25]. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer(5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Then 3 μL USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer [26]. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. Subsequently, sequencing was performed using the Illumina NovaSeq6000 sequencing platform and transcriptome assembly was performed using Trinity [27]. The default setting for min_kmer_cov was applied, while all other parameters were maintained at their default values.

Transcriptome data analysis

Gene functionality was annotated by referencing several databases, including NR (NCBI non-redundant protein sequences) [28], Pfam (Protein family) [29], KOG/COG/eggNOG (Clusters of Orthologous Groups of proteins) [30], Swiss-Prot (a manually annotated and reviewed protein sequence database) [31], KEGG (Kyoto Encyclopedia of Genes and Genomes) [32], and GO (Gene Ontology) [33]. DEGs among the CK, D1, and D2 groups were identified and analyzed using the DESeq R package (1.10.1), we used the following screening criteria: |fold change| > 2 and FDR < 0.01. The DEGs were subsequently subjected to functional clustering analysis based on the KEGG pathway and GO categories [34].

Metabolomic profile using liquid chromatography-mass spectrometry

The extraction of total metabolites from fresh tuber of CMD was performed as follows: 50 mg of the sample was weighed, and subsequently, 1000 μL of an extraction solution containing an internal standard (composed of methanol: acetonitrile: water in a 2:2:1, v/v ratio, with an internal standard concentration of 20 mg/L) was added. This mixture underwent vortex mixing for 30 seconds, followed by the addition of steel beads. The sample was then processed using a 45Hz grinder for 10 minutes, with an additional 10-minute ultrasonication step performed in an ice water bath. Afterward, the sample was left to stand at 20°C for 1 hour and subsequently centrifuged at 4°C and 12000 rpm for 15 minutes. Following centrifugation, 500 μL of supernatant was subjected to vacuum drying, followed by the addition of 160 μL of an extraction solution (acetonitrile: water, 1:1, v/v). A 30-second vortex mixing step was performed, followed by 10 minutes of sonication in an ice-water bath. Finally, the sample was centrifuged at 4°C, 12000rpm for 15 minutes.

The metabolites (supernatant, 120 μL) were separated using an LC−MS system (Acquity I-Class PLUS ultra-high performance liquid tandem Waters Xevo G2-XS QTOF high-resolution mass spectrometer). The liquid chromatograph column, Acquity UPLC HSS T3 column 1.8 μm, 2.1 × 100 mm), was purchased from Waters. The metabolite detection ofthe samples was conducted under positive and negative ion modes. The parameters of the ESI ion source were as follows: capillary voltage: 2000 V (positive ion mode) or −1500 V (negative ion mode); cone voltage: 30 V; ion source temperature: 150°C; desolvent gas temperature: 500°C; backflush gas flow rate: 50 L/h; and desolventizing gas flow rate: 800 L/h.

Metabolomics data analyses

The raw data obtained through MassLynx V4.2 underwent further processing steps, including peak extraction and alignment, utilizing the Progenesis QI software [35]. To assess the reproducibility of the samples, principal component analysis (PCA) was performed. Metabolites were annotated with reference to KEGG, the Human Metabolome Database (HMDB) [36], and lipid profiles. The multiplicity of differences was calculated and compared based on grouping information, and significant differences (p-values) for each compound were determined using the t-test. Additionally, we constructed an OPLSDA model [37] using the R language software package ropls, with model reliability verified through 200 alignment tests. Multiple cross-validation was performed to calculate the VIP value of the model. DEMs were identified using a combination of multiplicity of difference, p-value, and VIP value (FC > 1, p-value < 0.05, VIP > 1). Statistical analysis of metabolome data was performed on the BMKCloud platform (Biotechnology Co., Ltd., Beijing, China).

Conjoint analysis of transcriptomic and metabolomic data

In the correlation analysis, the relative abundance of genes and the relative abundance of metabolites were used to calculate the correlation coefficients. The key genes and metabolites in the important pathways for drought resistance in CMD were screened separately and characterized using Pearson correlation analysis. The expression of some key enzyme genes in flavonoid metabolism was correlated with the accumulation of flavonoids, and the correlation coefficient, r, indicated the degree of correlation between the two. r closer to 1 indicated a higher degree of positive correlation between the two, while r closer to -1 indicated a higher degree of negative correlation between the two.

Co-expression network analysis

Co-expression analysis was performed using DETFs detected in EP and CP and DEGs related to flavonoid biosynthesis pathway. The Pearson correlation coefficient (PCC) between DETFs and DEGs was calculated. TF-gene pairs with PCC ≥ 0.90 or ≤ -0.90 were selected for analysis, and the co-expression network was constructed and visualized using Cytoscape software.

Real-time quantitative PCR validation

Seven genes were randomly selected for qRT—PCR analysis. Specific primers were designed using Primer Premier 5.0 software, details are shown in S1 File. qRT—PCR was carried out in 20 μL reaction, 2×ChamQ SYBR Color qPCR Master Mix 10 μL, 2 μL each of upstream and downstream primers, 4 μL of template cDNA, and ddH₂O was added to 20 μL. The reaction parameters were as follows: 30 s at 95°C, 40 cycles of 15s at 95°C, and 30s at 60°C. The relative expression of mRNA was calculated using 2^-ΔΔCt, and all experiments were repeated three times [38].

Statistical analysis

Excel 2019 software was used for statistical analysis. qRT-PCR validated bar graphs show the mean ± SD of three independent experiments.

Results

Analysis of the differentially expressed genes (DEGs)

This study yielded a total of 112.15 GB of clean data from the 18 samples, with each sample demonstrating a Q30 base percentage exceeding 92.83%, indicating high data quality. Following the assembly process, we obtained a set of 37,247 Unigenes. Of these, 32,525 Unigenes received annotations, and notably, 22,167 Unigenes exceeded the 1 kb length threshold, further affirming the robustness of our transcriptome dataset. The resulting data are visually presented in Fig 1. Notably, the total count of DEGs showed an increasing trend in both EP and CP with higher levels of drought stress post-treatment. However, the extent of this increase varied between the two varieties. Specifically, DEGs in EP showed a substantial 130% increase, with the number of up-regulated genes rising from 449 to 1003 and down-regulated genes from 742 to 1747. In contrast, CP displayed a relatively stable count of up-regulated genes, while the number of down-regulated genes increased from 2338 to 3110. These findings highlight that drought stress treatment induced notable transcriptome alterations in CMD tissues, with the total number of DEGs positively correlated with the severity of drought.

[Figure omitted. See PDF.]

Red bars indicate an increase, blue bars indicate a decrease.

Differential gene function analysis

To assess the biological functions of the differential genes under drought stress treatments, we separately annotated the D1 and D2 treatment groups’ differential genes using the GO database, categorizing them into biological processes, cellular components and molecular functions. Notably, as shown in Fig 2, the top three major enriched categories were consistent in both D1 and D2, encompassing metabolic processes, cellular processes and single-organism processes in the biological process category, membrane-related components in the cellular component category, and catalytic activities, binding and transporter activities in the molecular function category. However, differences were also observed between D1 and D2 as they exhibited distinct enriched categories. D1 was characterized by enrichments related to locomotion, extracellular region parts and nutrient reservoir activities, while D2 did not exhibit these specific enrichments. Conversely, D2 displayed enrichments associated with growth, other organism parts, supramolecular complexes, transcription factor activities, protein binding and toxin activities, which were not observed in D1.

[Figure omitted. See PDF.]

(a) GO enrichment in the D1 treatment group, (b) GO enrichment in the D2 treatment group.

The KEGG enrichment analysis results show distinct pathway enrichments in response to drought stress treatments (Fig 3). Specifically, in the D1 group, we found enrichment of differential genes in pathways related to the biosynthesis of amino acids, phenylpropanoid biosynthesis and glycolysis/gluconeogenesis (carbon metabolism). On the other hand, in the D2 group, differential genes exhibited enrichment in pathways associated with plant-pathogen interaction, starch and sucrose metabolism, and phenylpropanoid biosynthesis. Significantly, both the D1 and D2 groups displayed substantial enrichment in the phenylpropanoid biosynthesis pathway, which serves as a crucial source for flavonoid biosynthesis—the primary active ingredient in CMD. Furthermore, the D2 group exhibited a significant enrichment in the flavonoid biosynthesis pathway. These findings strongly suggest that flavonoids represent a pivotal pathway in the adaptation of CMD to drought stress.

[Figure omitted. See PDF.]

(a) KEGG enrichment in the D1 treatment group, (b) KEGG enrichment in the D2 treatment group.

Metabolomic analysis

To investigate the inherent material distinctions between EP and CP in their responses to drought, we used metabolomics for the identification of differential metabolites. Three samples were taken from each treatment group for metabolomic analysis. PCA demonstrated a degree of data consistency among replicate samples within each group (S2 File), effectively showcasing the clear segregation of metabolites between the different CMD types. Furthermore, the OPLS-DA results yielded robust Q2Y values of 0.97 and 0.99, respectively (Fig 4). These high Q2Y values indicate the reliability and consistency of the metabolomics data. Metabolite identification in this study was based on metabolite data detected in both positive and negative ion modes, with a total of 23,642 peaks detected, of which 4,240 metabolites were annotated.

[Figure omitted. See PDF.]

Blue triangles indicate EPs, red circles indicate CPs. (a) OPLS-DA analy-sis in the D1 treatment group, (b) OPLS-DA analysis in the D2 treatment group.

The abundance of metabolites differed between EP and CP acclimatized to drought.

To investigate the differences in drought adaptation between EP and CP, we quantitatively analyzed and compared the measured metabolites. In the D1 treatment group, we identified 2,124 DEMs, with 1,103 up-regulated and 1,021 down-regulated metabolites. In the D2 treatment group, we found 2,691 DEMs, consisting of 1,274 up-regulated and 1,417 down-regulated metabolites (S3 File). Notably, as drought severity increased, both EP and CP exhibited an increase in the total number of DEMs, up-regulated DEMs, and down-regulated DEMs. Specifically, the total DEMs increased by 366 and 489, up-regulated DEMs by 238 and 253, and down-regulated DEMs by 128 and 236, respectively, in EP and CP. The clustered heat map of differential metabolites, as shown in Fig 5, indicates a high level of reproducibility among samples within each group, reinforcing the robustness of our findings.

[Figure omitted. See PDF.]

(a)D1 treatment group, (b)D2 treatment group.

EP and CP have different metabolic patterns of flavonoids during drought acclimation.

Many bioactive compounds are synthesized through phenylpropanoid metabolism, among which flavonoids are a crucial branch of this pathway. During drought adaptation, different CMD varieties accumulated substantial quantities of flavonoid secondary metabolites (Fig 6). Our metabolomic analysis identified several DEMs within the flavonoid biosynthesis pathway, including demethoxycurcumin, catechin, peonidin derivatives and Apiforol, among others (Fig 7). Notably, the accumulation patterns of flavonoids in EP and CP differed significantly. For instance, in EP, acacetin levels were found to be increased, while naringin derivatives (specifically naringin dihydrochalcone) were increased in EP but decreased in CP. Similarly, Leucocyanidin and kaempferol increased in CP but decreased in EP. These distinct patterns highlight the contrasting responses of EP and CP to drought stress, particularly concerning the accumulation of flavonoid metabolites.

[Figure omitted. See PDF.]

(a)D1 treatment group, (b)D2 treatment group. Sections boxed in red indicate flavonoid-related metabolic pathways.

[Figure omitted. See PDF.]

The expression profile of each gene was determined based on log2(FC) in the RNA-seq data.

Expression of flavonoid synthesis-related genes in two CMD varieties

To investigate the flavonoid accumulation mechanisms during CMD’s adaptation to drought, we constructed a gene expression heat map focusing on enzymes in the flavonoid biosynthesis pathway based on literature findings and RNA-seq data. Fig 8 illustrates our key observations, whereby in CP, D1 treatment led to an overall upregulation of PAL, 4CL, HCT and PKS gene expression. Consequently, this upregulation indirectly elevated the levels of trans-cinnamic acid, caffeic acid and curcumin derivatives (specifically, demethoxycurcumin and curcumin monoglucoside) (Fig 7). Additionally, p-coumaroyl-CoA levels, a direct precursor in flavonoid biosynthesis, were increased. Conversely, in EP, the D2 treatment resulted in overall increased CHS gene expression, indirectly leading to elevated levels of naringenin chalcone derivatives, specifically naringin dihydrochalcone. These findings shed light on how changes in gene expression contribute to the modulation of flavonoid accumulation in CMD during drought adaptation.

[Figure omitted. See PDF.]

The level of transcription of the enzyme is represented by a heat map. The genes encoding enzymes include phenylalanine ammonia-lyase (PAL), cinnamate 4-hydroxylase (C4H), 4-coumarate-CoA ligase (4CL), p-coumaroyl quinate/shikimate 3’-hydroxylase (C3’H), caffeoyl-CoA-O-methyltransferase (CCoAOMT), polyketide synthase (PKS), caffeoyl shikimate esterase (CSE), shikimate O-hydroxycinnamoyltransferase (HCT), chalcone synthase (CHS), chalcone isomerase (CHI), flavone synthase (FNS), caffeic acid-O-methyltransferase (COMT), flavone 3-hydroxylase (F3H), flavonol synthase (FLS), dihydroflavonol reductase (DFR), flavonoid 3’-hydroxylase (F3’H), leucoananthcyanidin reductase (LAR), anthocyanidin synthase (ANS), flavonoid 3’, 5’methyltransferase (FAOMT).

Correlation analysis between key DEGs expression and accumulation of flavonoids

Pearson correlation analysis was conducted to examine the relationship between the expression levels of DEGs and the contents of compounds within the flavonoid biosynthetic pathway (Fig 9). The results revealed significant correlations between the compounds in the synthesis pathway and the associated DEGs. Specifically, CHI (DN9092901) displayed a negative correlation with naringenin and apigenin (r < -0.83) while exhibiting a positive correlation with trans-cinnamic acid (r > 0.84). Consequently, the expression of CHI in CP was higher than that in EP, but the content of naringenin and apigenin was lower in CP than in EP. The expression of CSE (DN3397101) and HCT (DN10274801) was positively correlated (r < -0.83) with naringin dihydrochalcone, and the content of apigenin exhibited a positive correlation (r > 0.81) with these genes. Additionally, the accumulation of acacetin was positively correlated (r > 0.93) with HCT (DN885801 and DN11868801). Furthermore, COMT (DN878101) showed a positive correlation (r > 0.82) with the content of catechin and peonidin derivatives, while FLS (DN493101) displayed a negative correlation with these compounds (r < -0.84). These correlation findings provide valuable insights into the interplay between gene expression and the accumulation of specific flavonoid compounds in CMD during drought adaptation.

[Figure omitted. See PDF.]

Significant correlations are marked with * (p < 0.05) and ** (p < 0.01). Red color indicates positive correlation and green color indicates negative correlation.

Gene regulatory networks of flavonoids

Structural genes involved in flavonoid metabolism are typically under the regulatory control of transcription factors (TFs). Herein, a variety of TFs were identified among the DEGs in drought-stressed CMD (Fig 10), including AP2/ERF-ERF, AUX/IAA, GRAS, MYB, NAC, C2H2, C3H, TRAF, WRKY, B3, HMG, HSF, LOB and Tify. The number and types of TFs varied considerably across different drought levels. In the D1 treatment, the top 5 TFs were AP2/ERF-ERF, AUX/IAA, GRAS, MYB and NAC, while under D2 treatment, the top 5 TFs were AP2/ERF-ERF, WRKY, C2H2, C3H and NAC. To explore the relationships between these TFs and key enzyme genes, Pearson correlation analysis was conducted based on the LOG2FC of the transcripts of the DETFs and the key enzyme genes. We selected correlations with Pearson correlation coefficients |r| ≥ 0.9 and p < 0.05 for visualization purposes (Fig 11). The results demonstrated that DEGs encoding 4CL, PAL, PKS, CHS, HCT, CCoAOMT and CHI exhibited co-expression with various DETFs. The expression of these genes was regulated by numerous TFs, including WRKY, TIFY, MYB, C2H2, AUX/IAA, NAC and HSF. Overall, these findings elucidate the complex regulatory networks involving TFs and key enzyme genes in CMD’s response to drought stress, shedding light on the intricacies of flavonoid metabolism regulation.

[Figure omitted. See PDF.]

(a) D1 treatment group, (b) D2 treatment group.

[Figure omitted. See PDF.]

Co-expression networks connecting structural genes in flavonoid biosynthesis (circles) with TFs (diamonds).

Validation by qRT-PCR

To verify the accuracy of DEGs, seven differentially expressed genes (DN4609301(4CL), DN966021(CHS), DN1133803(HCT), DN2448301(F3H), DN10614501(PKS), DN878101(FAOMT), DN3397101(CSE)) were randomly selected from flavonoid metabolism (Fig 8) for qRT—PCR validation. As shown in Fig 12, there was no significant difference between the fold expression of CP and EP in qRT—PCR and RNA -seq, which highly confirmed the credibility of the RNA -seq results.

[Figure omitted. See PDF.]

Red bars indicate multiplicity of cp/ep expression in RNA-Seq; blue bars indicate multiplicity of cp/ep expression in qRT-PCR. Data are mean ± standard deviation (SD) (three biological replicates and three technical repeats).

Discussion

CMD is rich in nutrients and medicinal components, and is the raw material of many health products and proprietary Chinese medicines, which are in great demand in the market. Currently, the main focus is on the pharmacological effects of CMD, and few people have studied its abiotic resistance. Here, we comprehensively analyzed the transcriptomics and metabolomics of CMD under drought stress to reveal the differences in the adaptability of different CMD varieties to drought, particularly in terms of flavonoid biosynthesis and oxidase regulation. In the following, we discuss the functional categorization of gene expression and metabolites with varying abundance, which could serve as a valuable resource for enhancing our understanding on the intricate metabolic variations occurring during drought acclimation in different CMD varieties.

Metabolic regulation plays a pivotal role in enabling plants to uphold cellular osmotic potential during drought stress, with plant metabolic genes and metabolites participating in numerous metabolic pathways to adapt to such stressors [39]. For instance, previous research has revealed significant alterations in metabolic pathways, including the TCA cycle and glycolytic pathway, in response to drought stress in different soybean varieties, with notable transcriptional changes occurring in numerous genes [40]. Sprenger et al. similarly observed strong drought-induced responses in transcripts across different potato varieties, with more pronounced transcriptional changes detected in the sensitive variety [41]. In CMD, as drought progressed from moderate to severe, DEGs increased by 130% and DEMs by 366 in EP. In contrast, CP exhibited a 25% increase in DEGs and a 489 increase in DEMs under the same drought conditions. From this, it can be seen that drought stress increased the number of DEGs and DEMs, and the extent of the increase was strongly correlated with the severity of drought and the type of CMD, which is similar to the results of the above references. Our findings suggest that CP efficiently mobilizes a smaller number of genes to regulate a larger array of metabolites, indicating a more effective response to drought compared to EP and underscoring CP’s superior drought tolerance compared to EP.

Flavonoids, a category of polyphenolic compounds renowned for their antioxidant properties, have been implicated in the regulation of plant tolerance to diverse abiotic stressors. Recent investigations, such as the work conducted by Zhao et al. through transcriptome and metabolome analyses, suggest that phenylpropanoid/flavonoid metabolism in soybeans could serve as a pivotal metabolic pathway for drought tolerance in crops [42]. Similarly, Gai et al. demonstrated that elevating the levels of flavonoids, anthocyanins, flavonols and isoflavonoids could enhance drought tolerance in tea trees [43]. In regard to white clover, the accumulation of quercetin and kaempferol in response to water deficit has been linked to increased drought tolerance [44]. Furthermore, studies on rice have identified 4-hydroxycinnamic acid and ferulic acid as critical metabolites for drought tolerance, with DEGs involved in the biosynthesis of these metabolites emerging as promising candidate genes for enhancing drought tolerance [45]. We observed a large enrichment of DEGs in the phenylpropanoid biosynthetic pathway and flavonoid biosynthetic pathway under drought conditions. In particular, most metabolites in the flavonoid biosynthesis pathway in CP showed an upward trend and increased with increasing drought, including trans-cinnamic acid, caffeic acid and curcumin derivatives. This suggests that these substances may be key metabolites for CP adaptation to drought environments, and DEGs involved in the synthesis of these compounds may be able to contribute to the improvement of CMD drought tolerance. In addition, these metabolites have various biological activities [46–48], we can try to subject CMD to moderate drought stress to increase the content of its active ingredients to improve its medicinal and commercial value.

Transcription factors play pivotal roles in plants’ responses to drought stress. Existing literature highlights several examples of their involvement in enhancing drought tolerance. For instance, SbWRKY30 has been shown to regulate the drought stress response gene SbRD19 in sorghum, thus enhancing drought tolerance by increasing the capacity to scavenge reactive oxygen species (ROS) and protect plant cells [49]. In walnuts, the expression of TIFYs has been linked to the regulation of drought tolerance [50]. C2H2-type zinc finger proteins have demonstrated the ability to enhance drought tolerance in transgenic plants by elevating the levels of osmoregulatory substances [51]. MYB proteins, a crucial family of transcription factors, have been associated with improved drought tolerance by regulating flavonoid synthesis and cuticle formation, impacting leaf permeability [52]. Overexpression of Aux/IAA proteins has been found to modulate the ABA signaling pathway, enhance resistance to oxidative stress, and improve drought tolerance in Brachypodium [53]. Additionally, the NAC transcription factor has been shown to increase leaves’ water content, reduce hydrogen peroxide levels, and act as a regulator of drought tolerance in tomatoes [54]. HSF transcription factors promote flavonoid accumulation, ROS scavenging, abscisic acid-induced stomatal closure, and the expression of abscisic acid signaling-related genes under drought conditions, contributing to enhanced plant survival [55]. The co-expression network we constructed provides insights into the crucial roles of TFs in regulating gene expression in CMD. Numerous types of TFs have been reported to be intricately involved in the regulation of key genes responsible for flavonoid synthesis, including a range of families, such as WRKY, TIFY, C2H2, MYB, AUX/IAA, NAC, HSF, among others. Specially, the transcription factors WRKY and TIFY are prominent in the co-expression network (Fig 11). Metabolome and transcriptome co-expression analysis of two different disease resistance genotypes of P. edulis identified PeWRKY30 as a key TF co-expressed with flavonoid accumulation in yellow fruit P. edulis [56]. TIFY transcription factors in Tartary Buckwheat and walnut respond to drought stress and a variety of other abiotic stresses [50–57]. Those suggests that WRKY and TIFY transcription factors may also play a central role in CMD adaptation to drought stress. Notably, a single TF can control the synthesis of multiple proteins associated with flavonoids in CMD, and conversely, each protein can be influenced by several TFs. These intricate interactions illustrate the complexity of the regulatory mechanisms associated with flavonoid synthesis in CMD, urging the need for further experimental investigations to fully elucidate these relationships.

Conclusion

This study represents comprehensive exploration of transcriptomic and metabolomic responses in CMD under varying drought conditions. Drought stress induced notable changes in both gene expression and metabolite levels in CMD. We observed a notable enrichment of DEGs in the phenylpropanoid biosynthesis pathway, accompanied by an accumulation of flavonoid-related metabolites, particularly in creeping CMD. Furthermore, our investigation highlights the pivotal role of transcription factors, including WRKY, TIFY, C2H2, MYB, AUX/IAA, NAC and HSF, in regulating flavonoid synthesis in CMD. Interestingly, creeping CMDs probably displayed superior drought tolerance compared to erect CMDs, as evidenced by a smaller increase in DEGs and DEMs with escalating drought levels and a more pronounced surge in flavonoid synthesis pathway metabolites. Collectively, these findings highlight the differences in drought responses between creeping and erect CMD varieties, offering insights into potential strategies for enhancing drought resilience in CMD plant species.

Supporting information

S1 File. List of qRT-PCR primers in this study.

https://doi.org/10.1371/journal.pone.0313580.s001

(DOCX)

S2 File. PCA analysis of metabolites identified in D1 and D2.

https://doi.org/10.1371/journal.pone.0313580.s002

(DOCX)

S3 File. Number of DEMs after D1 and D2 drought treatments.

https://doi.org/10.1371/journal.pone.0313580.s003

(DOCX)

S4 File. Transcriptome quality control data, assembly data and functional annotation.

https://doi.org/10.1371/journal.pone.0313580.s004

(XLSX)

S1 Raw data.

https://doi.org/10.1371/journal.pone.0313580.s005

(ZIP)

References

1. 1. Zhou Y-F, Wang L-L, Chen L-C, Liu T-B, Sha R-Y, Mao J-W. Enrichment and Separation of Steroidal Saponins from the Fibrous Roots of Ophiopogon japonicus Using Macroporous Adsorption Resins. RSC Adv. 2019, 9 (12), 6689–6698.

* View Article

* Google Scholar

2. 2. Wang L, Jiang X-L, Zhang W-M, Li F, Khan A-A, Liu X, et al. Homo-Aro-Cholestane, Furostane and Spirostane Saponins from the Tubers of Ophiopogon japonicus. Phytochemistry. 2017, 136, 125–132.

* View Article

* Google Scholar

3. 3. Fan S, Zhang J, Xiao Q, Liu P, Zhang Y, Yao E, et al. Cardioprotective Effect of the Polysaccharide from Ophiopogon japonicus on Isoproterenol-Induced Myocardial Ischemia in Rats. Int J Biol Macromol. 2020, 147, 233–240.

* View Article

* Google Scholar

4. 4. Hou J, Zhong L, Liu J, Liu F, Xia C. Interaction of the Main Active Components in Shengmai Formula Mediated by Organic Anion Transporter 1 (OAT1). J Ethnopharmacol. 2022, 296, 115515. pmid:35777609

* View Article

* PubMed/NCBI

* Google Scholar

5. 5. Li Z, Wu Y-Y, Yu B-X. Methylophiopogonanone A, an Ophiopogon Homoisoflavonoid, Alleviates High-Fat Diet-Induced Hyperlipidemia: Assessment of Its Potential Mechanism. Braz J Med Biol Res. 2020, 53 (3), e9201. pmid:32130294

* View Article

* PubMed/NCBI

* Google Scholar

6. 6. Li N, Zhang J-Y, Zeng K-W, Zhang L, Che Y-Y, Tu P-F. Anti-Inflammatory Homoisoflavonoids from the Tuberous Roots of Ophiopogon japonicus. Fitoterapia. 2012, 83 (6), 1042–1045.

* View Article

* Google Scholar

7. 7. Yan X, Ye T, Hu X, Zhao P, Wang X. 58-F, a Flavanone from Ophiopogon japonicus, Prevents Hepatocyte Death by Decreasing Lysosomal Membrane Permeability. Sci Rep. 2016, 6, 27875.

* View Article

* Google Scholar

8. 8. Li W, Ji L, Tian J, Tang W, Shan X, Zhao P, et al. Ophiopogonin D Alleviates Diabetic Myocardial Injuries by Regulating Mitochondrial Dynamics. J Ethnopharmacol. 2021, 271, 113853. pmid:33485986

* View Article

* PubMed/NCBI

* Google Scholar

9. 9. Liu B, Li B, Chen G, Pan Y, Zhou D, Li N. Spirostane Saponins with a Rearranged A/B Ring System Isolated from the Rhizomes of Ophiopogon japonicus. Phytochemistry. 2022, 193, 112975.

* View Article

* Google Scholar

10. 10. Fang J, Wang X, Lu M, He X, Yang X. Recent Advances in Polysaccharides from Ophiopogon japonicus and Liriope Spicata Var. Prolifera. Int J Biol Macromol. 2018, 114, 1257–1266.

* View Article

* Google Scholar

11. 11. Zhang L, Wang Y, Wu F, Wang X, Feng Y, Wang Y. MDG, an Ophiopogon japonicus Polysaccharide, Inhibits Non-Alcoholic Fatty Liver Disease by Regulating the Abundance of Akkermansia Muciniphila. Int J Biol Macromol. 2022, 196, 23–34.

* View Article

* Google Scholar

12. 12. Wu Y, Bi S-X, Huang Z, Qi J, Yu B-Y. Novel Steroidal Saponins with Cytotoxic Activities from the Roots of Ophiopogon japonicus (L. f.) Ker-Gawl. RSC Adv. 2018, 8 (5), 2498–2505.

* View Article

* Google Scholar

13. 13. Chen M-H, Leong F, Gao S-J, Chen X, Lu J-J, Lin L-G, et al. Comparison of Ophiopogon japonicus and Liriope Spicata Var. Prolifera from Different Origins Based on Multi-Component Quantification and Anticancer Activity. Molecules. 2023, 28 (3), 1045.

* View Article

* Google Scholar

14. 14. Chen Z, Zhu B, Peng X, Li S, Zhao J. Quality Evaluation of Ophiopogon japonicus from Two Authentic Geographical Origins in China Based on Physicochemical and Pharmacological Properties of Their Polysaccharides. Biomolecules. 2022, 12 (10), 1491.

* View Article

* Google Scholar

15. 15. Liu L, Lyu J, Yang L, Gao Y, Zhao B. Using Pharmacokinetic-Pharmacodynamic Modeling to Study the Main Active Substances of the Anticancer Effect in Mice from Panax ginseng-Ophiopogon japonicus. Molecules. 2024, 29 (2), 334.

* View Article

* Google Scholar

16. 16. Zhang T, Guo Z, Cheng X, Xia R, Lai S, Liao L. Protective Properties of Ophiopogonin D in DSS-Induced Colitis: Insights into Microbiota Modulation. Inflammopharmacology. 2024. pmid:39039348

* View Article

* PubMed/NCBI

* Google Scholar

17. 17. Lin WL, Su WW, Cai XY, Luo LK, Li PB, Wang YG. Fermentation effects of oligosaccharides of Radix Ophiopogonis on alloxan-induced diabetes in mice. Int J Biol Macromol. 2011, 49(2):194–200. pmid:21549746

* View Article

* PubMed/NCBI

* Google Scholar

18. 18. D’Ippólito S, Rey-Burusco MF, Feingold SE, Guevara MG. Role of Proteases in the Response of Plants to Drought. Plant Physiol Biochem. 2021, 168, 1–9. pmid:34607206

* View Article

* PubMed/NCBI

* Google Scholar

19. 19. Ozturk M, Turkyilmaz Unal B, García-Caparrós P, Khursheed A, Gul A, Hasanuzzaman M. Osmoregulation and Its Actions during the Drought Stress in Plants. Physiol Plant. 2021, 172 (2), 1321–1335. pmid:33280137

* View Article

* PubMed/NCBI

* Google Scholar

20. 20. Verma S, Negi NP, Pareek S, Mudgal G, Kumar D. Auxin Response Factors in Plant Adaptation to Drought and Salinity Stress. Physiol Plant. 2022, 174 (3), e13714. pmid:35560231

* View Article

* PubMed/NCBI

* Google Scholar

21. 21. Razi K, Muneer S. Drought Stress-Induced Physiological Mechanisms, Signaling Pathways and Molecular Response of Chloroplasts in Common Vegetable Crops. Critical Reviews in Biotechnology. 2021, 41 (5), 669–691. pmid:33525946

* View Article

* PubMed/NCBI

* Google Scholar

22. 22. Dwivedi AK, Singh V, Anwar K, Pareek A, Jain M. Integrated Transcriptome, Proteome and Metabolome Analyses Revealed Secondary Metabolites and Auxiliary Carbohydrate Metabolism Augmenting Drought Tolerance in Rice. Plant Physiology and Biochemistry. 2023, 201, 107849. pmid:37393858

* View Article

* PubMed/NCBI

* Google Scholar

23. 23. Yang B, He S, Liu Y, Liu B, Ju Y, Kang D, et al. Transcriptomics Integrated with Metabolomics Reveals the Effect of Regulated Deficit Irrigation on Anthocyanin Biosynthesis in Cabernet Sauvignon Grape Berries. Food Chemistry. 2020, 314, 126170. pmid:31978717

* View Article

* PubMed/NCBI

* Google Scholar

24. 24. Qin Y, Wang F, Lu C, Wang F, Wen Y, Liu Y, et al. First Report of Tobacco Mild Green Mosaic Virus Infecting Rehmannia Glutinosa in China. Plant Dis. 2022. pmid:35286130

* View Article

* PubMed/NCBI

* Google Scholar

25. 25. Panailidou P, Lotos L, Sassalou C-L, Gagiano E, Pietersen G, Katis N, et al. First Report of Grapevine Virus H (GVH) in Grapevine in Greece. Plant Dis. 2021. pmid:33728963

* View Article

* PubMed/NCBI

* Google Scholar

26. 26. Li R-J, Gao X, Li L-M, Liu X-L, Wang Z-Y, Lü S. De Novo Assembly and Characterization of the Fruit Transcriptome of Idesia Polycarpa Reveals Candidate Genes for Lipid Biosynthesis. Front. Plant Sci. 2016, 7. pmid:27375655

* View Article

* PubMed/NCBI

* Google Scholar

27. 27. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-Length Transcriptome Assembly from RNA-Seq Data without a Reference Genome. Nat Biotechnol. 2011, 29 (7), 644–652. pmid:21572440

* View Article

* PubMed/NCBI

* Google Scholar

28. 28. Benson DA, Cavanaugh M, Clark K, Karsch-Mizrachi I, Lipman DJ, Ostell J, et al. GenBank. Nucleic Acids Res. 2013, 41 (Database issue), D36–42. pmid:23193287

* View Article

* PubMed/NCBI

* Google Scholar

29. 29. Mistry J, Chuguransky S, Williams L, Qureshi M, Salazar GA, Sonnhammer ELL, Tosatto SCE, Paladin L, Raj S, Richardson LJ, Finn RD, Bateman A. Pfam: The Protein Families Database in 2021. Nucleic Acids Res. 2021, 49 (D1), D412–D419. pmid:33125078

* View Article

* PubMed/NCBI

* Google Scholar

30. 30. Zhang Y, Liu X, Li Y, Liu X, Ma H, Qu S, et al. Basic Characteristics of Flower Transcriptome Data and Derived Novel EST-SSR Markers of Luculia Yunnanensis, an Endangered Species Endemic to Yunnan, Southwestern China. Plants (Basel). 2022, 11 (9), 1204. pmid:35567204

* View Article

* PubMed/NCBI

* Google Scholar

31. 31. Chen C, Xu L, Xia W, Qi C, Liu J, Hua B, et al. Metabolomics and Transcription Profiling of Pumpkin Fruit Reveals Enhanced Bioactive Flavonoids and Coumarins in Giant Pumpkin (Cucurbita maxima). J Agric Food Chem. 2023, 71 (27), 10459–10469.

* View Article

* Google Scholar

32. 32. Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000, 28 (1), 27–30. pmid:10592173

* View Article

* PubMed/NCBI

* Google Scholar

33. 33. The Gene Ontology Consortium. The Gene Ontology Resource: 20 Years and Still GOing Strong. Nucleic Acids Res. 2019, 47 (D1), D330–D338. pmid:30395331

* View Article

* PubMed/NCBI

* Google Scholar

34. 34. Zhang Z, Long Y, Yin X, Wang W, Li W, Chen T, et al. Metabolome and Transcriptome Analysis of Sulfur-Induced Kiwifruit Stem Laccase Gene Involved in Syringyl Lignin Synthesis against Bacterial Canker. J Agric Food Chem. 2023. pmid:37651104

* View Article

* PubMed/NCBI

* Google Scholar

35. 35. Buchfink B, Xie C, Huson DH. Fast and Sensitive Protein Alignment Using DIAMOND. Nat Methods. 2015, 12 (1), 59–60. pmid:25402007

* View Article

* PubMed/NCBI

* Google Scholar

36. 36. Zhou F, Xu D, Xiong S, Chen C, Liu C, Jiang A. Transcriptomic and Metabolomic Profiling Reveal the Mechanism Underlying the Inhibition of Wound Healing by Ascorbic Acid in Fresh-Cut Potato. Food Chemistry. 2023, 410, 135444. pmid:36641908

* View Article

* PubMed/NCBI

* Google Scholar

37. 37. Aru V, Pisano MB, Savorani F, Engelsen SB, Cosentino S, Cesare Marincola F. Metabolomics Analysis of Shucked Mussels’ Freshness. Food Chem. 2016, 205, 58–65. pmid:27006214

* View Article

* PubMed/NCBI

* Google Scholar

38. 38. Liu Q, Zhang Y, Dong X, Zheng L, Zhou Y, Gao F. Integrated Metabolomics and Transcriptomics Analysis Reveals That the Change of Apoplast Metabolites Contributes to Adaptation to Winter Freezing Stress in Euonymus Japonicus. Plant Physiology and Biochemistry. 2023, 202, 107924. pmid:37541019

* View Article

* PubMed/NCBI

* Google Scholar

39. 39. Kumar M, Kumar Patel M, Kumar N, Bajpai AB, Siddique KHM. Metabolomics and Molecular Approaches Reveal Drought Stress Tolerance in Plants. Int J Mol Sci. 2021, 22 (17), 9108. pmid:34502020

* View Article

* PubMed/NCBI

* Google Scholar

40. 40. Wang X, Song S, Wang X, Liu J, Dong S. Transcriptomic and Metabolomic Analysis of Seedling-Stage Soybean Responses to PEG-Simulated Drought Stress. IJMS. 2022, 23 (12), 6869. pmid:35743316

* View Article

* PubMed/NCBI

* Google Scholar

41. 41. Sprenger H, Kurowsky C, Horn R, Erban A, Seddig S, Rudack K, et al. The Drought Response of Potato Reference Cultivars with Contrasting Tolerance. Plant Cell & Environment. 2016, 39 (11), 2370–2389. pmid:27341794

* View Article

* PubMed/NCBI

* Google Scholar

42. 42. Zhao B, Zhang S, Yang W, Li B, Lan C, Zhang J, et al. Multi-omic Dissection of the Drought Resistance Traits of Soybean Landrace LX. Plant Cell & Environment. 2021, 44 (5), 1379–1398. pmid:33554357

* View Article

* PubMed/NCBI

* Google Scholar

43. 43. Gai Z, Wang Y, Ding Y, Qian W, Qiu C, Xie H, et al. Exogenous Abscisic Acid Induces the Lipid and Flavonoid Metabolism of Tea Plants under Drought Stress. Sci Rep. 2020, 10 (1), 12275. pmid:32704005

* View Article

* PubMed/NCBI

* Google Scholar

44. 44. Lee B-R, Kim K-Y, Jung W-J, Avice J-C, Ourry A, Kim T-H. Peroxidases and Lignification in Relation to the Intensity of Water-Deficit Stress in White Clover (Trifolium repens L.). Journal of Experimental Botany. 2007, 58 (6), 1271–1279.

* View Article

* Google Scholar

45. 45. Ma X, Xia H, Liu Y, Wei H, Zheng X, Song C, et al. Transcriptomic and Metabolomic Studies Disclose Key Metabolism Pathways Contributing to Well-Maintained Photosynthesis under the Drought and the Consequent Drought-Tolerance in Rice. Front. Plant Sci. 2016, 7. pmid:28066455

* View Article

* PubMed/NCBI

* Google Scholar

46. 46. Kang C, Kim J, Ju S, Cho H, Kim HY, Yoon I-S, et al. Colon-Targeted Trans-Cinnamic Acid Ameliorates Rat Colitis by Activating GPR109A. Pharmaceutics. 2022, 15 (1), 41. pmid:36678670

* View Article

* PubMed/NCBI

* Google Scholar

47. 47. Bai X, Li S, Liu X, An H, Kang X, Guo S. Caffeic Acid, an Active Ingredient in Coffee, Combines with DOX for Multitarget Combination Therapy of Lung Cancer. J Agric Food Chem. 2022, 70 (27), 8326–8337. pmid:35772797

* View Article

* PubMed/NCBI

* Google Scholar

48. 48. Shabbir U, Rubab M, Tyagi A, Oh D-H. Curcumin and Its Derivatives as Theranostic Agents in Alzheimer’s Disease: The Implication of Nanotechnology. Int J Mol Sci. 2020, 22 (1), 196. pmid:33375513

* View Article

* PubMed/NCBI

* Google Scholar

49. 49. Yang Z, Chi X, Guo F, Jin X, Luo H, Hawar A, et al. SbWRKY30 Enhances the Drought Tolerance of Plants and Regulates a Drought Stress-Responsive Gene, SbRD19, in Sorghum. J Plant Physiol. 2020, 246–247, 153142. pmid:32112957

* View Article

* PubMed/NCBI

* Google Scholar

50. 50. Liu X, Yu F, Yang G, Liu X, Peng S. Identification of TIFY Gene Family in Walnut and Analysis of Its Expression under Abiotic Stresses. BMC Genomics. 2022, 23 (1), 190. pmid:35255828

* View Article

* PubMed/NCBI

* Google Scholar

51. 51. Han G, Lu C, Guo J, Qiao Z, Sui N, Qiu N, et al. C2H2 Zinc Finger Proteins: Master Regulators of Abiotic Stress Responses in Plants. Front Plant Sci. 2020, 11, 115. pmid:32153617

* View Article

* PubMed/NCBI

* Google Scholar

52. 52. Wang X, Niu Y, Zheng Y. Multiple Functions of MYB Transcription Factors in Abiotic Stress Responses. Int J Mol Sci. 2021, 22 (11), 6125. pmid:34200125

* View Article

* PubMed/NCBI

* Google Scholar

53. 53. Su P, Sui C, Li J, Wan K, Sun H, Wang S, et al. The Aux/IAA Protein TaIAA15-1A Confers Drought Tolerance in Brachypodium by Regulating Abscisic Acid Signal Pathway. Plant Cell Rep. 2023, 42 (2), 385–394. pmid:36566287

* View Article

* PubMed/NCBI

* Google Scholar

54. 54. Thirumalaikumar VP, Devkar V, Mehterov N, Ali S, Ozgur R, Turkan I, et al. NAC Transcription Factor JUNGBRUNNEN1 Enhances Drought Tolerance in Tomato. Plant Biotechnol J. 2018, 16 (2), 354–366. pmid:28640975

* View Article

* PubMed/NCBI

* Google Scholar

55. 55. Wang N, Liu W, Yu L, Guo Z, Chen Z, Jiang S, et al. HEAT SHOCK FACTOR A8a Modulates Flavonoid Synthesis and Drought Tolerance. Plant Physiol. 2020, 184 (3), 1273–1290. pmid:32958560

* View Article

* PubMed/NCBI

* Google Scholar

56. 56. Ma F, Zhou H, Yang H, Huang D, Xing W, Wu B, et al. WRKY Transcription Factors in Passion Fruit Analysis Reveals Key PeWRKYs Involved in Abiotic Stress and Flavonoid Biosynthesis. Int J Biol Macromol. 2024, 256 (Pt 1), 128063. pmid:37963507

* View Article

* PubMed/NCBI

* Google Scholar

57. 57. Zhao Z, Meng G, Zamin I, Wei T, Ma D, An L, et al. Genome-Wide Identification and Functional Analysis of the TIFY Family Genes in Response to Abiotic Stresses and Hormone Treatments in Tartary Buckwheat (Fagopyrum tataricum). Int J Mol Sci. 2023, 24 (13), 10916.

* View Article

* Google Scholar

Citation: Cheng T, Lin J, Zhou X, Wang H, Zhou X, Huang X, et al. (2025) Integrative metabolomics and transcriptomics profiling reveals differential expression of flavonoid synthesis in Ophiopogon japonicus (L. f.) Ker-Gawl. in adaptation to drought. PLoS ONE 20(1): e0313580. https://doi.org/10.1371/journal.pone.0313580

About the Authors:

Tingting Cheng

Roles: Writing – original draft