Integrating transcriptome and metabolome analyses of the response to cold stress in pumpkin (Cucurbita maxima)

Abstract

Cucurbita maxima belong to the genus Cucurbita and are of nutritional and economic importance. Physiological activity, transcriptome, and metabolome analyses of leaf samples from the C. maxima inbreding line IL7 treated at 5 °C and 25 °C were performed. Cold stress resulted in a significant increase in the malondialdehyde content, relative electrical conductivity, soluble protein, sugar content, and catalase activity. A total of 5,553 differentially expressed genes were identified, of which 2,871 were up-regulated and 2,682 down-regulated. In addition, the transcription of differentially expressed genes in the plant hormone signal transduction pathway and transcription factor families of AP2/ERF, bHLH, WRKY, MYB, and HSF was activated. Moreover, 114 differentially expressed metabolites were identified by gas chromatography time-of-flight mass spectrometry, particularly through the analysis of carboxylic acids and derivatives, and organooxygen compounds. The demonstration of a series of potential metabolites and corresponding genes highlighted a comprehensive regulatory mechanism. These findings will provide novel insights into the molecular mechanisms associated with the response to cold stress in C. maxima.

About the Authors:

Fengmei Li

Roles Writing – original draft

Affiliation: Shandong Provincial Key Laboratory of Biochemical Engineering, College of Marine Science and Biological Engineering, Qingdao University of Science & Technology, Qingdao, Shandong, China

Xiuping Lu

Roles Data curation, Software

Affiliation: Shandong Provincial Key Laboratory of Biochemical Engineering, College of Marine Science and Biological Engineering, Qingdao University of Science & Technology, Qingdao, Shandong, China

Pengfei Duan

Roles Software, Visualization

Affiliation: Shandong Provincial Key Laboratory of Biochemical Engineering, College of Marine Science and Biological Engineering, Qingdao University of Science & Technology, Qingdao, Shandong, China

Yanjiao Liang

Roles Formal analysis

Affiliation: Shandong Provincial Key Laboratory of Biochemical Engineering, College of Marine Science and Biological Engineering, Qingdao University of Science & Technology, Qingdao, Shandong, China

Jian Cui

Roles Funding acquisition, Methodology

* E-mail: [email protected]

Affiliation: Qingdao Institute of Agricultural Science Research, Qingdao, Shandong, China

ORCID logo https://orcid.org/0000-0001-8961-7354

Introduction

Low temperature is one of the major environmental factors limiting agricultural productivity and the geographic distribution of many plant species [1]. During early spring, low temperatures are not conducive to plant survival and development, leading to a decline in fruit yield and quality. Some plant species have developed a remarkable ability to adapt to low temperatures [2]. Through intensive research, a variety of mechanisms employed to respond to temperature have been identified, but the way plants sense low temperature remains unknown [3].

After plant exposure to cold stress, a series of physiological and biochemical changes at the molecular and cellular levels are induced to adapt and survive under cold environments [4,5]. These changes include increased accumulation of antioxidant enzymes and non-enzymatic molecules which may neutralize or counteract the harmful effects of reactive oxygen specie (ROS) and protect the plants from cold stress [6]. In the process of cold stress response, a series of osmoregulatory metabolites including soluble sugar (e.g., sucrose, glucose and galactose) and low molecular weight compounds (such as proline, glycine betaine), these substances can help alleviate osmotic stress, expansion, water uptake and metabolic activity of plant cells [7,8]. To date, considerable research has been conducted on the physiological and biochemical mechanisms of plant cold resistance. For example, physiological analyses of Santalum album L. leaves treated at 4 °C for 48 h demonstrated that cold stress induces the accumulation of malondialdehyde (MDA), proline, and soluble sugar (SS), while increasing the levels of superoxide dismutase (SOD) and peroxidase (POD) activity [9]. Ma et al. analyzed the physiological response of Brassica rapa L. held at 4 °C for 24 h to reveal that the POD and catalase (CAT) activity, SS content, and MDA content of the cultivars Longyou-7 and Lenox increased under cold stress [10]. These results indicated that different plants exhibit different physiological responses to cold stress.

Plant transcriptomics and metabolomics are now used to explore the genes and metabolites that enable cold survival. These resources are greatly expanding our understanding of important molecular processes. The study of global gene expression using RNA-Seq can help researchers to unveil the regulatory genes and key pathways that may induce cold resistance. Hence, an increasing number of cold sensitive genes and transcription factors have been identified and analyzed in plants under cold stress [9,11–14]. WRKY, NAC, MYB, AP2/ERF, and bZIP were the most abundant transcription factor families under cold stress. These transcription factor families are known for role in ABA-induced signaling pathways under low temperature [15,16].

Moreover, multiple pathways have been recognized to contribute directly to the cold stress response in plants. Most pathways such as flavonoid biosynthesis, phenylpropanoid biosynthesis, plant hormone signal transduction, cutin, suberine and wax biosynthesis, pentose and glucuronate interconversions, phenylalanine metabolism and starch and sucrose metabolism were significantly enriched KEGG in rubber tree (Hevea brasiliensis) after cold stress [17]. The plant hormone signal transduction pathway has been found to provide tolerance to cold stress in several species including Arabidopsis thaliana [12] and Strawberry (Fragaria × ananassa) [13,14]. Recently, the best approach has been to integrate transcriptomics and metabolomics methods to reveal the complex reaction process that occurs under low temperature stress and provide more information about this response. In recent years, this integrative method has become a useful tool for producing a large amount of molecular information about global cellular events and has been used to study peach, tobacco, seaweed, tomato, and grape, etc [11,13,16–18]. Although several molecular mechanisms related to cold stress have been revealed in these plants, little is known about the response of pumpkin, Cucurbita maxita, to cold stress.

In agricultural production, C.maxita is an important horticultural crop [19] and it is worldwide-cultureted crop because of its remarkable nutritious and economic value. Moreover, C. maxita is often used as important rootstocks for grafting watermelon, cucumber, and melon, it can increase tolerance to cold stress [20]. Although seven cold stress-related glutathione transferase (GST) genes were explored in the cold-tolerant line (C. maxima) and are putative candidates for use in breeding cold-tolerant crops varieties [21], the molecular studies of pumpkin are still in their infancy. Studies of differences in transcriptional and metabolic levels of pumpkin under cold stress can improve our understanding of cold stress responses mechanism.

Materials and methods

Plant material, cold treatments, and physiological analyses

The C. maxima inbred line IL7, provided by the Qingdao Academy of Agricultural Sciences of Shandong Province, China, was used in this study. IL7 has been purified by more than eight consecutive generations of self-pollination.

Seeds were sown in a universal horticultural peat-based soil mixture. Plants were grown for 2 weeks in cultivation chambers at 25/19 °C (day/night) and 80% relative humidity; the digital illumination reached 24,000 Lux. These plants were grown inside cultivation chambers (Ningbo Safe Instrument Co., Ltd., Ningbo, China) set at 25 °C (control). After 3 weeks of growth, corresponding to the 3–4 leaf stage, cold treatment at 5 °C (16 h of light/8 h of dark) was applied and visible disease symptoms were confirmed after 24 h. Leaves were collected 24 h after the onset of the cold stress treatment and prepared in triplicate.

Physiological analyses

The experiments were conducted in a completely randomized block design with three replications so that each data point was the average of three independent samples. Membrane permeability, expressed by relative electric conductivity (REC), was determined by a modified method [22]. The Soil Plant Analysis Development (SPAD) value, the index of total chlorophyll, was detected using a portable chlorophyll meter (SPAD-502 Plus, Minolta Camera Co. Ltd., Japan) [23]. The determination of the soluble sugar (SS) and soluble protein contents, SOD activity, and MDA content was carried out using a kit according to the instructions provided by the Nanjing Jiancheng Biology Engineering Institute. The POD and CAT activities were detected using a kit according to the instructions provided by the Suzhou Comin Biotechnology Co., Ltd.

Total RNA extraction, library preparation, and Illumina sequencing for transcriptome analysis

Total RNA was extracted from young C. maxima leaves of the control (untreated) and cold-stressed (treated) (5 °C, 24 h) plants using the mirVana miRNA Isolation Kit (Ambion). Three biological replicates were used and the RNA was stored at -70 °C until the RNA-seq library preparation was complete. Real biological variability can be better simulated and statistical accuracy improved by using biological replicates. The integrity of the total RNA from the C. maxima leaves was analyzed using an Agilent 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA). We pooled the high-quality reads from the samples to perform de novo transcriptome assembly. Samples with an RNA integrity number ≥ 7 were subjected to subsequent analyses. A TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, USA) was used to construct the libraries. The libraries were sequenced using an Illumina HiSeqTM 2500 sequencing platform, and 150 bp paired-end reads were generated.

Identification and annotation of DEGs for transcriptome analysis

Raw reads were removed with poly-N and low-quality reads were removed to obtain clean reads. The RNA-Seq samples were then mapped to the C. maxima genome (https://www.ncbi.nlm.nih.gov/assembly/GCF_002738345.1/) as the reference genome using Tophat [24]. HTSeq software (Heidelberg, Germany) was used to obtain the read counts of each gene [25] and Cufflinks [26] software computed the gene expression fragments per kilobase of transcript per million mapped reads (FPKM) values, and the read counts for each gene were obtained from the HTSeq-count [25]. The linear association between the gene expression levels of the sequenced samples was evaluated using Pearson’s correlation coefficient. The distances between the samples were calculated using clustering methods to perform principal component analyses of the gene expression levels.

For quantification at the transcript level, the FPKM and read count value of each transcript were calculated using bowtie2 [27] and eXpress. The DEGs were identified using the estimate Size Factors and negative binomial (nbinom) tests in the DESeq function of R package [28]. The software enabled us to obtain the number of reads of a gene in each sample, normalize the data using DESeq, and use the nbinom test function to calculate the p-value and fold change values. RNA-seq means from the same gene in two samples were considered to indicate statistically differential expression when the p-value was <0.05 and fold change was >2. A hierarchical cluster analysis of DEGs was performed to explore the transcript expression pattern. Both a gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEGs were respectively performed using R, based on the hypergeometric distribution. The results were reassembled using the Tophat-Cufflinks platform [26]. The results from the Cuff compare software were compared with the gene annotation information from the reference to obtain sequenced gene structure extensions and identify new transcripts. The alternative splicing analysis of differentially regulated transcripts, isoforms, or exons was performed using Asprofile [28].

Verification of differential expression using real-time quantitative reverse transcription-PCR

The total RNA was extracted from the leaves of pumpkin seedlings treated with cold stress for 0 or 24 h as described above. Quantification was performed using a two-step reaction process involving reverse transcription (RT) and PCR. In the first step, 0.5 μg RNA and 2 μL of 4 × gDNA wiper Mix were added to nuclease-free H2O to give a final volume of 8 μL. The reaction was carried out on a GeneAmp® PCR System 9700 (Applied Biosystems, USA) for 2 min at 42 °C. Then, 2 μl of 5 × HiScript II Q RT SuperMix II were added. The reaction was then carried out for a further 15 min at 50 °C, followed by 5 s at 85 °C. The reaction mixture was then diluted in nuclease-free H2O and stored at -20 °C. Real-time PCR was performed using a LightCycler® 480II real-time PCR instrument (Roche, Switzerland). The reactions were conducted in a mixture volume of 10 μl that included 5 μl of 2 × ChamQ SYBR qPCR Master Mix. There were three biological replicates of each sample. The PCR thermocycling program consisted of an initial incubation at 95 °C for 30 s, followed by 40 cycles at 95 °C for 10 s, and 60 °C for 30 s. The β-actin gene was used as an internal control, and gene specific primers were designed using Primer Premier 5.0 software. The primers are presented in S1 Table. After normalizing the mRNA expression level, it was calculated using the 2-ΔΔCt method [29].

Metabolites

Ten samples from young leaves of the control (untreated) and cold treatment (5 °C, 24 h) groups were analyzed using a gas chromatography time-of-flight mass spectrometry system. The freeze-dried samples were crushed at 45 Hz for 4 min using a mixer mill. A 50 mg aliquot of sample was added to a 1.5 mL Eppendorf tube, and the sample was extracted by adding 1mL of 75% methanol and placing the tube on a shaker overnight at 4 °C. Then, the samples were centrifuged for 15 min (12,500 rpm, 4 °C). The resulting supernatants were transferred to 2 mL sample vials and stored at -80 °C until the ultra-high performance liquid chromatography (UHPLC)-tandem mass spectrometry (MS/MS) analysis. A quality control sample was prepared by mixing equal aliquots of the supernatants from all samples.

LC-MS analysis

The UHPLC separation was performed using a Waters ACQUITY UPLC HSS T3 column (100 × 2.1 mm, 1.8 μm). We used 0.1% formic acid in H2O as mobile phase A and acetonitrile as mobile phase B. During the LC/MS experiment, the MS/MS spectra were acquired by a Triple TOF 6600 mass spectrometer in the information-dependent acquisition (IDA) mode.

In the IDA mode, the data acquisition software (Analyst TF 1.7, AB Sciex) collects and triggers the acquisition of MS/MS spectra based on MS data, depending on preselected criteria. For each cycle, 12 precursor ions with intensities greater than 100 were selected. The collision energy (CE) was 30 V, and the product ion accumulation time of each 12 MS/MS events was 50 ms.

Data were collected in four segments according to the mass ranges 50–300, 290–600, 590–900, and 890–1500. We used an AB Sciex QTrap 6500 mass spectrometer to develop the assay.

Results

Physiological changes of IL7 in response to cold stress

The REC, SS content, soluble protein content, POD, and CAT activity in cold-treated IL7 showed a significant increase at 24 h compared with the control, and there was a significant decrease in the SPAD index of cold-treated IL7 relative to the control (Fig 1A–1D, 1F and 1G). However, there was no significant change in SOD activity and MDA content in IL7 (Fig 1E and 1H). In particular, the SS content in the leaves of C. maxima increased 1.88-fold in fighting the cold stress. These results strongly suggested that the cold-stressed IL7 seedlings actively mounted a stress response after they were exposed to 5 °C.

[Figure omitted. See PDF.]

Fig 1. Physiological changes of IL7 under control and cold conditions.

A, REC; B, SPAD index; C, Soluble sugar content; D, Soluble protein content; E, SOD; F, POD, G, CAT; H, MDA. Error bars denote standard error of the mean. Bars with the same small letters do not statistically differ by the Duncan’s multiple range test at P < 0.05.

https://doi.org/10.1371/journal.pone.0249108.g001

Transcriptome analysis of IL7 under cold stress

To further elucidate the molecular basis for cold responses in IL7, a transcriptomic analysis was performed using RNA-seq with a paired-end sequencing strategy on an Illumina HiSeq 2500 platform. Six cDNA libraries representing the control (Sample_IL7C1, Sample_IL7C2, Sample_IL7C3) and, low temperature-treated (Sample_IL7T1, Sample_IL7T2, Sample_IL7T3) treatments samples were constructed. A The full-scale sequencing analysis from of six samples is shown in Table 1. A total of 47.35 to 48.06 million (M) uniquely mapped clean reads were used for the subsequent analysis analyses. A total of 41.53 Gb of clean data were obtained, with an average of 6.92 Gb per library. The percentages of Q30 base was were greater than 94%, and the average GC content was 45.81%. The sample_IL7s were mapped to the C. maxima genome assembly. The mapped reads in six sample_IL7s were relatively high in the range of 98.40%–98.92%.

[Figure omitted. See PDF.]

Table 1. Summary of sequencing data of IL7.

https://doi.org/10.1371/journal.pone.0249108.t001

Identification of DEGs in response to cold stress

We identified genes expressed in the control and cold-treated leaves of IL7 and found the FPKM value of most genes to be >1 (S1 Fig). The high similarity between the biological replicates of both the control and cold-treated samples indicated that the RNA-seq results were consistent. A correlation analysis of samples also supported consistency (r > 0.99) (S2 Fig). In order to obtain reliable gene expression profiles, reads with a log2 |fold change| > 1 and FKPM values > 10 were selected to annotate the DEGs. A total of 5,553 DEGs with 2,871 up-regulated and 2,682 down-regulated genes between the control and cold treated IL7 samples were detected. In addition, 208 up-regulated genes and 44 down-regulated genes had a log2 |fold change| > 4.

The transcript abundance of all DEGs after cold stress treatment is shown by a hierarchical heat map (Fig 2). The results also showed a series of changes in gene expression in IL7 under cold stress.

[Figure omitted. See PDF.]

Fig 2. Heat map of IL7T vs IL7C (P < 0.05, log2 |Fold Change| > 1).

The deepening color in the heat map indicates an increase in the number of gene transcripts at that growth point.

https://doi.org/10.1371/journal.pone.0249108.g002

Analysis of transcription factors

As critical regulators, transcription factors (TFs) play key roles in the regulation of gene expression under cold stress in plants. Our RNA-seq results showed that 354 DEGs (221 up-regulated and 133 down-regulated genes) were identified and assigned to various TFs. Among these DEGs, 15 TF genes (two bHLHs, seven ERFs, four WRKYs, and two MYBs) were significantly up-regulated by more than 20-fold during cold treatment. Interestingly, most of the predicted bHLH93 (LOC111476397) and ERF1B (LOC111499095) genes showed dramatic levels of activation, as they were up-regulated by 138- and 109-fold, respectively. A total of 252 DEGs were comprised of AP2/ERF, bHLH, MYB, WRKY, HSF, and GATA TFs (Table 2). These genes were the most frequently identified as being up or down-regulated, implying that they might regulate resistance to cold stress. Notably, the large number of AP2/EFR family members indicated that different members may play specific roles in the response to cold stress. These results suggested that these TF families could be involved in the regulation of cold stress.

[Figure omitted. See PDF.]

Table 2. Major transcription factor in DEGs.

https://doi.org/10.1371/journal.pone.0249108.t002

GO and KEGG enrichment results

To investigate the functions of the DEGs detected in the IL7 response to cold stress, we analyzed the DEGs by GO term enrichment. GO analysis showed that the DEGs were categorized into functional groups consisting of 1,963 GO terms in biological processes, 405 GO terms in cellular components, and 1,186 GO terms in molecular function subcategories in the comparison between treated and control IL7 samples. The biological process category represented 23 functional subcategories, out of which metabolic process, response to stimulus, and regulation of biological process were predominantly enriched under cold stress (Fig 3). The cellular component and molecular function categories comprise 13 and 14 functional categories, respectively (Fig 3). In the cellular component category, cell, cell part, and organelle were the most highly represented groups. In the molecular function category, binding, catalytic activity, and nucleic acid binding transcription factor activity were the three major groups associated with the C. maxima unigenes (Fig 3).

[Figure omitted. See PDF.]

Fig 3. Functional classification of C. maxima unigenes by using Gene Ontology (GO).

https://doi.org/10.1371/journal.pone.0249108.g003

KEGG pathway classification assigned 1,529 unigenes were enriched 120 pathways that are involved in metabolism, genetic information, organismal systems, environmental information processing, cellular processes, and human diseases (Fig 4). The most abundant term was signal transduction, followed by the carbohydrate metabolism and amino acid metabolism. The top 20 pathways are shown in Fig 5. The first prominent pathways of significance were plant hormone signal transduction (Fig 5) phenylpropanoid biosynthesis and phenylalanine metabolism, suggesting that these pathways may be closely related to cold stress in IL7.

[Figure omitted. See PDF.]

Fig 4. Functional classification of unigenes by using KEGG terms.

https://doi.org/10.1371/journal.pone.0249108.g004

[Figure omitted. See PDF.]

Fig 5. The top 20 KEGG pathways.

https://doi.org/10.1371/journal.pone.0249108.g005

To verify the RNA-seq data, 15 genes were selected for qPCR validation, including DREB3, ERF1B, and GH3.6, which have already been reported to be related to cold responses. The results of the RNA-seq and qPCR analyses were similar (Fig 6). The data demonstrated that there was a strong positive correlation (R2 = 0.963) between the RNA-seq data and qPCR results.

[Figure omitted. See PDF.]

Fig 6. The quantitative RT-PCR results from IL7.

Relative quantification was obtained through 2-ΔΔCT method using Actin reference gene. Data represent the average from three biological replicates and the error bars indicate the standard deviation (± SD).

https://doi.org/10.1371/journal.pone.0249108.g006

Differential expression of plant hormone signal transduction

Hormones are signaling molecules that regulate gene expression under cold stress. A total of 148 DEGs were associated with auxin (AUX/indole-3-acetic acid) (IAA), cytokinin (CK), gibberellin (GA), abscisic acid (ABA), ethylene (ET), brassinosteroid (BR), jasmonic acid (JA), and salicylic acid (SA)-mediated pathways (S2 Table).

In detail, 53 genes related to auxin were differentially expressed, including 20 AUX/indole-3-acetic acids (Aux/IAA), 21 small Auxin-Up RNAs (SAURs), six AUX1s, four Gretchen Hagen 3s (GH3s), and one auxin response factor (ARF) involved in auxin-induced signal transduction pathways. Most of the gene expression were down-regulated, indicating that the auxin signaling pathway was inhibited. Thirteen DEGs were significantly enriched in cytokinin-mediated pathway; one of the five up-regulated A-ARR genes displayed the highest level of expression with a fold change of 39.15. Six DEGs (five up-regulated and one down-regulated) were involved in GA-mediated pathways; the GIDI and GELLA genes showed prominent up-regulation with 16.34-fold and 7.17-fold changes, respectively. Up-regulated expression of these genes suggested cytokinin-mediated pathway was activted Twenty-three DEGs were related to ABA hormone-mediated signaling; the PP2C gene (LOC111467353) showed up-regulation of 16.1-fold, suggesting that it might play a vital role against cold stress. Moreover, 11 up-regulated and two down-regulated DEGs were found to be regulated via ET. An ERF1/2 gene associated with the ET-signaling pathway was the most highly expressed and exhibited a 109.0-fold change in expression. Fifteen DEGs were associated with BR-mediated pathways (four up-regulated and 11 down-regulated); TCH4, which exhibited a 35.83-fold change in expression was highly induced under cold stress. Furthermore, ten DEGs (nine up-regulated and one down-regulated) were found to be linked with JA-mediated pathways; JAR1, JAZ, and MYC2 were induced by more than three-fold. Three NPR1 genes (two up-regulated and one down-regulated) and five TGA genes (three up-regulated and two down-regulated) were found to be linked with SA-mediated pathways; one NPR1 gene exhibited a 3.02-fold change in expression and one TGA gene displaying a 3.11-fold change in expression were highly expressed in IL7 under cold stress.

Differential expression of genes related to phenylpropanoid biosynthesis

The phenylpropanoid biosynthesis pathway is activated under abiotic stresses including drought, high/low temperature, heavy metal toxicity, salinity, and ultraviolet radiation; thus, this pathway has the potential to scavenge excess ROS and maintain the cellular redox status [30,31].

In the present study, KEGG enrichment indicated the presence of 59 DEGs (31 up-regulated and 18 down-regulated) and demonstrated that phenylpropanoid biosynthesis (Fig 7) played important roles in the resistance of C. maxima to cold stress. The fist enzyme in the phenylpropanoid biosynthesis pathway is phenylalanine ammonia lyase (PAL), which catalyzes the conversion of L-Phe into cinnamic acid [32]. In all of the plants studied, PAL is present in multiple isoforms; thus, this enzyme must belong to a multi-gene family [33]. Likewise, nine PAL isoform genes were up-regulated in C. maxima under cold stress. The trans-cinnamate 4-monooxygenase-like (C4H) enzyme, is a member of the cytochrome P450 monooxygenase superfamily and is also the second key enzyme in the phenylpropanoid biosynthesis pathway, which catalyzes the hydroxylation of trans-cinnamic acid to p-coumaric acid [34]. In our study, only two C4Hs of C. maxima were up-regulated with a log2 |fold change| = 7.8 and log2 |fold change| = 4.1 under cold treatment. The third enzyme in the phenylpropanoid biosynthesis pathway is 4-coumarate CoA ligase, which may play a central role in regulating the entry of hydroxycinnamic acid into the subsequent biosynthetic pathway. Under cold treatment in C. maxima, three 4-coumarate CoA ligases were up-regulated and four down-regulated. In addition, cinnamoyl-CoA reductase is the most important and final step in the synthesis of cinnamaldehyde, and two cinnamoyl-CoA reductases were more highly expressed in the leaves of C. maxima under cold stress.

[Figure omitted. See PDF.]

Fig 7. Phenylpropanoid biosynthesis pathway.

https://doi.org/10.1371/journal.pone.0249108.g007

Detection and profile analysis of metabolites under cold stress

To understand the metabolome changes associated with IL7 under cold stress, an untargeted metabolite analysis was performed on young leaves using UHPLC quantitative time-of-flight MS. The samples consisted of two groups: IL7C (control) and IL7T (treatment). As seen from the results of the orthogonal partial least squares discriminant analysis score map, the differences in the two groups of samples were very significant and all were within the 95% confidence interval (S3 and S4 Figs). A total of 114 metabolites showed differential accumulation (using VIP ≥ 1 and T-test P < 0.05,) in IL7, with 34 increased and 80 decreased metabolites, including 20 different classes of substances (S2 Table). Among increased metabolites, carboxylic acids and derivatives, organooxygen compounds, fatty acyls, and cinnamic acids and derivatives were themoste prominent in C. maxima (Fig 8). More increase of carboxylic acids and derivatives, such as L-phenylalanine and L-tyrosine, were important members of amino acids metabolism (S2 Table).

[Figure omitted. See PDF.]

Fig 8. Differentially accumulated metabolites in the C. maxima after cold stress mainly belonged to four classes.

https://doi.org/10.1371/journal.pone.0249108.g008

Discussion

The response to cold stress in plants is a process by which plants adapt to the environment through gene expression and metabolism [35–38]. In this study, the physiological changes in C. maxima after cold treatment revealed considerable variation related to cold stress. Significant variations in gene expression and metabolite levels after cold stress were detected using RNA-Seq and LC-MS analysis, respectively.

Transcription factor in response to cold stress

The transcriptional regulation of plant genes is an important step in the response to cold stress [39–42]. Thus, it is important to elucidate the complex regulatory mechanisms behind gene expression in plants to reveal the molecular basis of plant tolerance to low temperatures [43–46]. Members of several TF families, including AP2/ERF [47–49], WRKY [50,51], and MYB [14,15,52,53] play crucial roles in diverse stress responses. The AP2/ERF family, which is divided into the AP2/ERF, RAV and DREB subfamilies, plays important roles in multiple stress responses [54,55]. Since the first AP2/ERF TF was cloned from A. Thaliana [56], the biological functions of some AP2/ERF genes have been identified. For instance, overexpression of the rice OsDRE1F gene can increase tolerance to salt, drought, and low temperatures in both Arabidopsis and rice [57], and ectopic overexpression of a citrus ERF1 gene enhances cold tolerance in tobacco [58], while the birch BpERF13 gene demonstrates a positive regulation of cold stress by upregulating a CBF gene and reducing the number of reactive oxygen species (ROS) [59]. In the present study, AP2/ERF TFs were represented by 53 up-regulated and 18 down-regulated genes under cold stress. These results are in conformity with the previous reports on low temperature stress regulation due to presence of AP2/ERF family [60–62]. WRKY TFs represent a valuable family in the resistance to abiotic stresses such as cold, heat, drought, and salt. For example, a total of 17 members of this family were induced by chilling in japonica rice [60]. And a set of WRKY genes were also found to be chilling-regulated in indica rice, of which 12 were commonly induced, 22 were exclusively up-regulated in cold susceptible variety and 18 were exclusively up-regulated in cold tolerant variety genotype [63]. Moreover, the roles of VaWRKY12 from Vitis amurensis in cold stress tolerance were confirmed through its ectopic overexpression in Arabidopsis and grapevine [64]. A VbWRKY32 in the seedling stage of Verbena bonariensis upregulates the transcriptional level of cold response genes, which increases the antioxidant enzyme activities and the contents of osmotic adjustment substance, thereby improving the survival ability under cold stress [65]. Taken together, 33 WRKY DEGs of C. maxima may contribute to the further study of cold tolerance mechanisms. In addition, members of the MYB family have been shown to be pivotal factors in regulating responses to abiotic stress. Some MYB TFs are involved in both cold and chilling stresses. For example, the expression of LcMYB4 in Leymus chinensis was rapidly induced by cold treatment, and positively modulates chilling and freezing tolerance in Arabidopsis [66]. However, VcMYB4a of Vaccinium corymbosum transcript levels gradually decreased during the cold treatment [67]. Moreover, SlMYB102 in tomato enhances antioxidant enzyme activity and promotes the expression of cold resistance under low temperature stress conditions [52]. Similar to previous studies, the expressions of MYB genes in C. maxima were induced with the upregulated and downregulated pattern during cold stress. In the future, the mechanism of cold tolerance of these genes needs to be further studied. Recently, some studies have shown that AP2/ERF, WRKY, and MYB TFs are also involved in the expression of genes related to the phenylpropanoid pathway [68–70]. In our results, major TF families such as AP2/ERF, MYB, WRKY, and MYB were represented by a large number of DEGs and therefore have a significant influence on the resistance to cold stress in C. maxima. Thus, these TFs are likely to be involved in complex mechanisms associated with temperature regulation.

Plant hormone signaling involved in cold stress

Plant hormones play important roles in mediating development processes and defense response to various biotic and abiotic stresses [71,72]. Many genes indeed displayed altered expression levels in several hormone-mediated branches, such as ABA, IAA, CK, GA, ET, BR, JA, and SA. Obviously, most of the genes in auxin signaling pathway in C. maxima leaves exposed to cold stress were inhibited, which generally are obligated to cell elongation, division, cycle, and growth [13,73]. Auxin is an endogenous small molecule with an incredibly large impact on plant growth and and response to stress in concert with other hormonal pathways [74–76]. The Auxin response factor (ARF), Gretchen Hagen 3 (GH3), and Auxin/indole-3-acetic acid (Aux/IAA), small Auxin-Up RNA (SAUR) gene families, are key components of the auxin signaling pathway, and function as regulators of plant growth and response to stress [77,78]. More recently, auxin signal in strawberry (Fragaria×ananassa) was severely blocked after low temperature treatment, indicating that auxin signal plays a very important role in cold stress [13]. This conclusion is consistent with our finding in C. maxima. CK is strongly cross-talks with auxins in plant developmental processes [2]. In Arabidopsis, CKs initiate signalling utilising Arabidopsis histidine kinase (AHK)2, AHK3 and AHK4, which act as CK receptors. Upon CK binding, the receptors transfer the CK signal via Arabidopsis histidine phosphotransfer proteins (AHPs) to nuclear localised Arabidopsis response regulators (ARRs), TFs that regulate expression of CK target genes [79]. In Arabidopsis, there are 23 ARR proteins, and transcript levels of A-type ARR5, ARR7 and ARR15 were induced in response to cold in a AHK-dependent manner [80,81]. Some A-type and B-type ARRs in C. maxima were induced expression, sugessting that CK signaling was activeted after cold stress. Quite on the contrary, as opposed to GAs, BRs are thought to positively control cold stress responses, because it has been shown that BR application increases cold tolerance of many plant species, including chilling-sensitive crops such as maize and Cucumis sativus (cucumber) [81–84]. Previous studies showed that GA and JA can cooperatively regulate diverse aspects of plant growth, development, and defense through DELLA and JAZ proteins [85,86], and JA prioritizes defense over growth by interfering with GA signaling cascade [48]. ABA, a central regulator of cold stress signalling with emerging roles in the CBF-dependent pathway [2]. Cold treatment of pumpkin activated hormone signaling pathway, which was confirmed by some differential expression of genes.

Carboxylic acids and derivatives metabolites in response to cold stress

In addition, the associated accumulation of specific amino acids has been correlated with increased tolerance to stressful environmental conditions [87]. Phenylalanine and tyrosine have different functions in stress resistance and are central to plant metabolism [88]. The increase in amino acids in A. thaliana, tomato, and other plants under during stress acclimation indicate that the accumulation of amino acids is related to the degradation of proteins [13,89,90]. It has been reported that the exogenous application of phenylalanine can promote chilling tolerance in tomato by activating ROS-scavenging systems; an increase in the abundance of the phenylalanine is also induced by low temperature [91]. L-phenylalanine (L-phe) is produced from phenylpyruvic acid-mediated amination by aspartate aminotransferase [92]. Detection of metabolism showed that the accumulation of L-Phe and tyrosine provided important materials to the phenylpropanoid biosynthesis pathway. Moreover, the pathway is also involved in cold tolerance in winter rapeseed and grape [10,18]. The current results indicated that the amino acid metabolism in C. maxima was induced and helped in the resistance to a cold environment.

Conclusions

Cold stress has an important influence on plant growth and development, leading to physiological, biochemical, metabolic, and molecular changes. In the present study, we comprehensively analyzed the changes in the transcriptome and metabolome from C. maxima leaves during cold stress. The results revealed that cold stress broadly activated the pathways for plant hormone signal transduction, phenylpropanoid biosynthesis. Important genes and metabolic processes in these pathways during cold stress could provide valuable information for further functional characterization in future studies.

Supporting information

S1 Fig. The FPKM expression distribution map of IL7T vs IL7C.

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

(TIF)

S2 Fig. Correlation coefficient heat map between samples of IL7.

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

(TIF)

S3 Fig. The Score scatter plot of OPLS-DA for group IL7T vs IL7C.

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

(TIF)

S4 Fig. Permutation test of OPLS-DA model for group IL7T vs IL7C.

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

(TIF)

S1 Table. Information of gene specific primers for qRT-PCR.

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

(DOCX)

S2 Table. DEGs involving in plant hormone signal transduction.

https://doi.org/10.1371/journal.pone.0249108.s006

(DOCX)

Acknowledgments

The authors are grateful to other laboratory members for technical advices and helpful discussions.

Citation: Li F, Lu X, Duan P, Liang Y, Cui J (2021) Integrating transcriptome and metabolome analyses of the response to cold stress in pumpkin (Cucurbita maxima). PLoS ONE 16(5): e0249108. https://doi.org/10.1371/journal.pone.0249108

References

1. John R, Anjum N A, Sopory S K, Akram N A, Ashrafet M. Some key physiological and molecular processes of cold acclimation. Bio Plantarum. 2016;60:603–618.

2. Eremina M, Rozhon W, Poppenberger B. Hormonal control of cold stress responses in plants. Cell Mol Life Sci. 2016;73:797–810. pmid:26598281

3. Knight M R, Knight H. Low-temperature perception leading to gene expression and cold tolerance in higher plants. New Phytol. 2012;195:737–751. pmid:22816520

4. Bohnert HJ, Nelson DE, Jensen RG. Adaptations to Environmental Stresses. Plant Cell. 1995; 7(7):1099–1111. pmid:12242400

5. Browse J, Xin Z. Temperature sensing and cold acclimation. Curr Opin Plant Biol. 2001; 4(3):241–6. pmid:11312135

6. Hussain S, Khan F, Hussain HA, Nie L. Physiological and biochemical mechanisms of seed priming-induced chilling tolerance in rice cultivars. Front Plant Sci. 2016;7:116. pmid:26904078

7. Verslues PE, Agarwal M, Katiyar-Agarwal S, Zhu J, Zhu JK. Methods and concepts in quantifying resistance to drought, salt and freezing, abiotic stresses that affect plant water status. Plant J. 2006;45(4):523–539. pmid:16441347

8. Li YY, Wang XW, Ban QY, Zhu XX, Jiang CJ, Wei CL, et al. Comparative transcriptomic analysis reveals gene expression associated with cold adaptation in the tea plant Camellia sinensis. BMC genomics. 2019;20(1):624. pmid:31366321

9. Zhang X, Teixeira da Silva JA, Niu M, Li M, He C, Zhao J, et al. Physiological and transcriptomic analyses reveal a response mechanism to cold stress in Santalum album L. leaves. Sci Rep. 2017;7:42165. pmid:28169358

10. Ma L, Coulter JA, Liu L, Zhao Y, Chang Y, Pu Y, et al. Transcriptome analysis reveals key cold-stress-responsive genes in winter rapeseed (Brassica rapa L.). Int J Mol Sci. 2019;20:1071. pmid:30832221

11. Jin J, Zhang H, Zhang J, Liu P, Chen X, Li Z, et al. Integrated transcriptomics and metabolomics analysis to characterize cold stress responses in Nicotiana tabacum. BMC Genomics. 2017;18:496. pmid:28662642

12. Sharma R, Singh G, Bhattacharya S, Singh A. Comparative transcriptome meta-analysis of Arabidopsis thaliana under drought and cold stress. PLoS One. 2018;13:e0203266. pmid:30192796

13. Zhang Y, Zhang Y, Lin Y, Luo Y, Wang X, Chen Q, et al. A transcriptomic analysis reveals diverse regulatory networks that respond to cold stress in strawberry (Fragaria×ananassa). Int J Genomics. 2019;2019:7106092. pmid:31467865

14. Cao M, Wang D, Mao Y, Kong F, Bi G, Xing Q, et al. Integrating transcriptomics and metabolomics to characterize the regulation of EPA biosynthesis in response to cold stress in seaweed Bangia fuscopurpurea. PLoS One. 2017;12:e0186986. pmid:29240755

15. Zhou MQ, Shen C, Wu LH, Tang KX, Lin J. CBF-dependent signaling pathway: a key responder to low temperature stress in plants. Crit Rev Biotechnol. 2011;31:186–192. pmid:20919819

16. Agarwal PK, Jha B. Transcription factors in plant and ABA dependent and independent abiotic stress signaling. Biol Plant. 2010; 54(2):201–12.

17. Wang K, Yin XR, Zhang B, Grierson D, Xu CJ, Chen KS. Transcriptomic and metabolic analyses provide new insights into chilling injury in peach fruit. Plant Cell Environ. 2017;40:1531–1551. pmid:28337785

18. Maoz I, De Rosso M, Kaplunov T, Vedova AD, Sela N, Flamini R, et al. Metabolomic and transcriptomic changes underlying cold and anaerobic stresses after storage of table grapes. Sci, Rep. 2019;9:2917. pmid:30814549

19. Paris HS, Yonash N, Portnoy V, Mozes-Daube N, Tzuri G, Katzir N. Assessment of genetic relationships in Cucurbita pepo (Cucurbitaceae) using DNA markers. Theor Appl Genet, 2003;106: 971–978. pmid:12671744

20. Schwarza D, Rouphael Y, Colla G, Venemad JH. Grafting as a tool to improve tolerance of vegetables to abiotic stresses: Thermal stress, water stress and organic pollutants. Sci Hortic. 2010;127:162–171.

21. Abdul Kayum M, Nath UK, Park JI, Biswas MK, Choi EK, et al. Genome-wide identification, characterization, and expression profiling of glutathione s-transferase (GST) family in pumpkin reveals likely role in cold-stress tolerance. Genes (Basel). 2018;9:84 pmid:29439434

22. Mao L, Wang G, Zhu C, Panga H. Involvement of phospholipase D and lipoxygenase in response to chilling stress in postharvest cucumber fruits. Plant Sci. 2007;172:400–405

23. Jiang C, Johkan M, Hohjo M, Tsukagoshi S, Maruo T. Acorrelation analysis on chlorophyll content and SPAD value intomato leaves. Hort Research. 2017;71:37–42.

24. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:R364. pmid:23618408

25. Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–169. pmid:25260700

26. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7:562–578. pmid:22383036

27. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–359. pmid:22388286

28. Florea L, Song L, Salzberg SL. Thousands of exon skipping events differentiate among splicing patterns in sixteen human tissues. F1000Res. 2013;2:188 pmid:24555089

29. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods (San Diego, Calif.). 2001;25:402–408. pmid:11846609

30. Dong CJ, Li L, Cao N, Shang QM, Zhang ZG. Roles of phenylalanine ammonia-lyase in low temperature tolerance in cucumber seedlings. Chinese Journal of Applied Ecology. 2015;26:2041–2049. pmid:26710630

31. Sharma A, Shahzad B, Rehman A, Bhardwaj R, Landi M, Zheng B. Response of phenylpropanoid pathway and the role of polyphenols in plants under abiotic Stress. Molecules. 2019;24:2452. pmid:31277395

32. Zhang X, Liu CJ. Multifaceted regulations of gateway enzyme phenylalanine ammonia-lyase in the biosynthesis of phenylpropanoids. Mol Plant. 2015;8(1):17–27. pmid:25578269

33. Cochrane FC, Davin LB, Lewis NG. The Arabidopsis phenylalanine ammonia lyase gene family: kinetic characterization of the four PAL isoforms. Phytochemistry. 2004;65:1557–1564. pmid:15276452

34. Chen A, Chai Y, Li J, Chen L. Molecular cloning of two genes encoding cinnamate 4-hydroxylase (C4H) from oilseed rape (Brassica napus). J Biochem Mol Biol. 2007;40:247–260. pmid:17394776

35. Dong H, Li W, Tang W, Li Z, Zhang D, Niu Y. Yield, quality and leaf senescence of cotton grown at varying planting dates and plant densities in the Yellow River Valley of China. Field Crops Res. 2006; 98:106–115.

36. Wong C, Ling YS, Wee JLS, Mujahid A, Müller M. A comparative UHPLC-Q/TOF-MS-based eco-metabolomics approach reveals temperature adaptation of four Nepenthes species. Sci Rep. 2020;14;10(1):21861. pmid:33318532

37. Ye T, Wang Y, Feng YQ, Chan Z. Physiological and metabolomic responses of bermudagrass (Cynodon dactylon) to alkali stress. Physiol Plant. 2021;171(1):22–33. pmid:32909624

38. Wang J, Guo J, Zhang Y, Yan X. Integrated transcriptomic and metabolomic analyses of yellow horn (Xanthoceras sorbifolia) in response to cold stress. PLoS One. 2020;24;15(7):e0236588. pmid:32706804

39. Chinnusamy V, Zhu J, Zhu JK. Cold stress regulation of gene expression in plants. Trends in Plant Science. 2007, 12(10): 444–451. pmid:17855156

40. Singh K, Foley RC, Oñate-Sánchez L. Transcription factors in plant defense and stress responses. Curr Opin Plant Biol. 2002;5(5):430–436. pmid:12183182

41. Liu Q, Sun C, Han J, Li L, Wang K, Wang Y, et al. Identification, characterization and functional differentiation of the NAC gene family and its roles in response to cold stress in ginseng, Panax ginseng C.A. Meyer. PLoS One. 2020;11;15(6):e0234423. pmid:32525906

42. Park S, Lee CM, Doherty CJ, Gilmour SJ, Kim Y, Thomashow MF. Regulation of the Arabidopsis CBF regulon by a complex low-temperature regulatory network. Plant J. 2015;82(2):193–207. pmid:25736223

43. Liu Y, Dang P, Liu L, He C. Cold acclimation by the CBF-COR pathway in a changing climate: Lessons from Arabidopsis thaliana. Plant Cell Rep. 2019;38(5):511–519. pmid:30652229

44. Shi Y, Ding Y, Yang S. Molecular Regulation of CBF Signaling in Cold Acclimation. Trends Plant Sci. 2018;23(7):623–637. pmid:29735429

45. Barrero-Gil J, Salinas J. Gene Regulatory Networks Mediating Cold Acclimation: The CBF Pathway. Adv Exp Med Biol. 2018;1081:3–22. pmid:30288701

46. Kidokoro S, Yoneda K, Takasaki H, Takahashi F, Shinozaki K, Yamaguchi-Shinozaki K. Different Cold-Signaling Pathways Function in the Responses to Rapid and Gradual Decreases in Temperature. Plant Cell. 2017;29(4):760–774. pmid:28351986

47. Chen J, Zhou Y, Zhang Q, Liu Q, Li L, Sun C, et al. Structural variation, functional differentiation and expression characteristics of the AP2/ERF gene family and its response to cold stress and methyl jasmonate in Panax ginseng C.A. Meyer. PLoS One. 2020;15(3):e0226055. pmid:32176699 eCollection 2020.

48. Yang Y, Al-Baidhani HHJ, Harris J, Riboni M, Li Y, Mazonka I, et al. DREB/CBF expression in wheat and barley using the stress-inducible promoters of HD-Zip I genes: impact on plant development, stress tolerance and yield. Plant Biotechnol J. 2020;18(3):829–844. pmid:31487424

49. Song Y, Cui H, Shi Y, Xue J, Ji C, Zhang C, et al. Genome-wide identification and functional characterization of the Camelina sativa WRKY gene family in response to abiotic stress. BMC Genomics. 2020;21(1):786. pmid:33176698

50. Li Y, Zhang L, Zhu P, Cao Q, Sun J, Li Z, et al. Genome-wide identification, characterisation and functional evaluation of WRKY genes in the sweet potato wild ancestor Ipomoea trifida (H.B.K.) G. Don. under abiotic stresses. BMC Genet. 2019;20(1):90. pmid:31795942

51. Ullah A, Ul Qamar MT, Nisar M, Hazrat A, Rahim G, Khan AH, et al. Characterization of a novel cotton MYB gene, GhMYB108-like responsive to abiotic stresses. Mol Biol Rep. 2020;47(3):1573–1581. pmid:31933260

52. Wang M, Hao J, Chen X, Zhang X. SlMYB102 expression enhances low-temperature stress resistance in tomato plants. PeerJ. 2020;8:e10059. pmid:33083130

53. Sun HJ, Luo ML, Zhou X, Zhou Q, Sun YY, Ge WY, et al. PuMYB21/PuMYB54 coordinate to activate PuPLDbeta1 transcription during peel browning of cold-stored "Nanguo" pears. Hortic Res. 2020 Sep 1;7:136. pmid:32922808

54. Vazquez-Hernandez M, Romero I, Escribano MI, Merodio C, Sanchez-Ballesta MT. Deciphering the Role of CBF/DREB Transcription Factors and Dehydrins in Maintaining the Quality of Table Grapes cv. Autumn Royal Treated with High CO2 Levels and Stored at 0 °C. Front Plant Sci. 2017;8:1591. pmid:28970842

55. Jofuku KD, den Boer BG, Van Montagu M, Okamuro JK. Control of Arabidopsis flower and seed development by the homeotic gene APETALA2. Plant Cell. 1994;6:1211–1225 pmid:7919989

56. Wang Q, Guan Y, Wu Y, Chen H, Chen F, Chu C. Overexpression of a rice OsDREB1F gene increases salt, drought, and low temperature tolerance in both Arabidopsis and rice. Plant Mol Biol. 2008;67:589–602. pmid:18470484

57. Ma Y, Zhang L, Zhang J, Chen Jiao, Wu T, Zhu S, et al. Expressing a Citrus ortholog of Arabidopsis ERF1 enhanced cold-tolerance in tobacco. Sci. Hortic. 2014;174:65–76.

58. Lv K, Li J, Zhao K, Chen S, Nie J, Zhang W, et al. Overexpression of an AP2/ERF family gene, BpERF13, in birch enhances cold tolerance through upregulating CBF genes and mitigating reactive oxygen species. Plant Sci. 2020;292:110375. pmid:32005381

59. Yun KY, Park MR, Mohanty B, Herath V, Xu F, Mauleon R, et al. Transcriptional regulatory network triggered by oxidative signals configures the early response mechanisms of japonica rice to chilling stress. BMC Plant Biol. 2010;10:16. pmid:20100339

60. Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K. AP2/ERF family transcription factors in plant abiotic stress responses. Biochim Biophys Acta. 2012;1819(2):86–96. 1819;2012:86–96. pmid:21867785

61. Dubouzet JG, Sakuma Y, Ito Y, Kasuga M, Dubouzet EG, Miura S, et al. OsDREB genes in rice, Oryzasativa L., encode transcription activators that function in drought-, high-salt- and cold-responsive gene expression. Plant J. 2003;33:751–63. pmid:12609047

62. Pradhan SK, Pandit E, Nayak DK, Behera L, Mohapatra T. Genes, pathways and transcription factors involved in seedling stage chilling stress tolerance in indica rice through RNA-Seq analysis. BMC Plant Biol. 2019;19(1):352. pmid:31412781

63. Zhang L, Zhao T, Sun X, Wang Y, Du C, Zhu Z, et al. Overexpression of VaWRKY12, a transcription factor from Vitis amurensis with increased nuclear localization under low temperature, enhances cold tolerance of plants. Plant Mol Biol. 2019;100(1–2):95–110. pmid:31011887

64. Wang MQ, Huang QX, Lin P, Zeng QH, Li Y, Liu QL, et al. The Overexpression of a Transcription Factor Gene VbWRKY32 Enhances the Cold Tolerance in Verbena bonariensis. Front Plant Sci. 2020;29;10:1746. pmid:32063911 eCollection 2019.

65. Li X, Jia J, Zhao P, Guo X, Chen S, Qi D, et al. LcMYB4, an unknown function transcription factor gene from sheepgrass, as a positive regulator of chilling and freezing tolerance in transgenic Arabidopsis. BMC Plant Biol. 2020;20(1):238. pmid:32460695

66. Zhang CY, Liu HC, Zhang XS, Guo QX, Bian SM, Wang JY, et al. VcMYB4a, an R2R3-MYB transcription factor from Vaccinium corymbosum, negatively regulates salt, drought, and temperature stress. Gene. 2020;757:144935. pmid:32653482

67. Zhang CY, Liu HC, Zhang XS, Guo QX, Bian SM, Wang JY, et al. VcMYB4a, an R2R3-MYB transcription factor from Vaccinium corymbosum, negatively regulates salt, drought, and temperature stress. Gene. 2020;757:144935. pmid:32653482

68. Meng J, Wang B, He G, Wang Y, Tang X, Wang S, et al. Metabolomics integrated with transcriptomics reveals redirection of the phenylpropanoids metabolic flux in Ginkgo biloba. J. Agric. Food Chem. 2019;67:3284–3291. pmid:30802049

69. Wang Z, Cui Y, Vainstein A, Chen S, Ma H. Regulation of Fig (Ficus carica L.) Fruit Color: Metabolomic and Transcriptomic Analyses of the Flavonoid Biosynthetic Pathway. Front Plant Sci. 2017;8:1990. pmid:29209349

70. Peleg Z. and Blumwald E. Hormone balance and abiotic stress tolerance in crop plants.Current Opinion in Plant Biology, 2011;14(3);290–295. pmid:21377404

71. Verma V. Ravindran P. and Kumar P. P. Plant hormonemediated regulation of stress responses. Plant hormonemediated regulation of stress responses. 2016;16(1)86. pmid:27079791

72. Chaiwanon J, Wang W, Zhu JY, Oh E, Wang ZY. Information integration and communication in plant growth regulation.Cell. 2016;164(6)1257–1268. pmid:26967291

73. Rahman A. Auxin: a regulator of cold stress response. Auxin: a regulator of cold stress response. 2013;147(1);28–35. pmid:22435366

74. Paterlini A. Uncharted routes: exploring the relevance of auxin movement via plasmodesmata. Biol Open. 2020;9(11):bio055541. pmid:33184092

75. Shibasaki K, Uemura M, Tsurumi S, Rahman A. Auxin response in Arabidopsis under cold stress: underlying molecular mechanisms.The Plant Cell. 2009;21(12);3823–3838. pmid:20040541

76. Chapman EJ, Estelle M. Mechanism of auxin-regulated gene expression in plants. Annu. Rev. Genet. 2009; 43:265–285. pmid:19686081

77. Tiwari S, Hagen BG, Guilfoyle TJ. The roles of auxin response factor domains in auxin-responsive transcription.Plant Cell, 2003; 15:533–543. pmid:12566590

78. Kumar MN, Verslues PE. Stress physiology functions of the Arabidopsis histidine kinase cytokinin receptors. Physiol Plant. 2015;154(3):369–80. pmid:25263537

79. Jeon J, Kim NY, Kim S, Kang NY, Novák O, Ku Sj, et al. A subset of cytokinin two-component signaling system plays a role in cold temperature stress response in Arabidopsis. J. Biol. Chem. 2010; 285:23371–23386. pmid:20463025

80. Shi Y, Tian S, Hou L, Huang X, Zhang X, Guo H,et al. Ethylene signaling negatively regulates freezing tolerance by repressing expression of CBF and type-A ARR genes in Arabidopsis. Plant Cell. 2012;24:2578–2595. pmid:22706288

81. Xia XJ, Wang YJ, Zhou YH, Tao Y, Mao WH, Shi K,et al. Reactive oxygen species are involved in brassinosteroid-induced stress tolerance in cucumber. Plant Physiol. 2009;150:801–814. pmid:19386805

82. Singh I, Kumar U, Singh SK, Gupta C, Singh M, Kushwaha SR. Physiological and biochemical effect of 24-epibrassinoslide on cold tolerance in maize seedlings. Physiol Mol Biol. Plants. 2012;18:229–236. pmid:23814437

83. Jiang YP, Huang LF, Cheng F, Zhou YH, Xia XJ, Mao WH, et al. Brassinosteroids accelerate recovery of photosynthetic apparatus from cold stress by balancing the electron partitioning, carboxylation and redox homeostasis in cucumber. Physiol Plant.2013;148:133–145. pmid:22998725

84. Kagale S, Divi UK, Krochko JE, Keller WA, Krishna P. Brassinosteroid confers tolerance in Arabidopsis thaliana and Brassica napus to a range of abiotic stresses. Planta.2007;225:353–364. pmid:16906434

85. Qi T, Huang H, Wu D, Yan J, Qi Y, Song S, Xie D. Arabidopsis DELLA and JAZ proteins bind the WD-repeat/bHLH/MYB complex to modulate gibberellin and jasmonate signaling synergy.The Plant Cell.2014;26, 3:1118–1133. pmid:24659329

86. Yang D.-L., Yao J., Mei C.-S. et al. Plant hormone jasmonate prioritizes defense over growth by interfering with gibberellin signaling cascade. Plant hormone jasmonate prioritizes defense over growth by interfering with gibberellin signaling cascade. Proc Natl Acad Sci USA.2012;109(19):E1192–E1200. pmid:22529386

87. Hildebrandt TM. Synthesis versus degradation: directions of amino acid metabolism during Arabidopsis abiotic stress response. Plant Mol Biol. 2018;98:121–135. pmid:30143990

88. Tzin V, Galili G. New insights into the shikimate and aromatic amino acids biosynthesis pathways in plants. Mol Plant. 2010; 3:956–972. pmid:20817774

89. Usadel B, Bläsing OE, Gibon Y, Poree F, Höhne M, Günter M, et al. Multilevel genomic analysis of the response of transcripts, enzyme activities and metabolites in Arabidopsis rosettes to a progressive decrease of temperature in the non-freezing range. Plant Cell Environ. 2008;31:518–547. pmid:18088337

90. Gonzalez C, Zanor MI, Ré MD, Otaiza S, Asis R, Valle EM, et al. Chilling tolerance of Micro-Tom fruit involves changes in the primary metabolite levels and in the stress response. Postharvest Biol. Tec. 2019;148:58–67.

91. Aghdam MS, Moradi M, Razavi F, Rabiei V. Exogenous phenylalanine application promotes chilling tolerance in tomato fruits during cold storage by ensuring supply of NADPH for activation of ROS scavenging systems. Sci. Hortic. 2019;246:818–825

92. Chao Y, Lo T, Luo N. Selective production of L-aspartic acid and L-phenylalanine by coupling reactions of aspartase and aminotransferase in Escherichia coli. Enzyme Microb Technol. 2000;27:9–25. pmid:10862897