De novo transcriptome sequencing and gene expression profiling of Magnolia wufengensis in response to cold stress

Background Magnolia wufengensis is a new species of Magnolia L. and has considerable ornamental and economic value due to its unique characteristics. However, because of its characteristic of poor low temperature resistance, M. wufengensis is hardly popularization and application in the north of China. Furthermore, the mechanisms of gene regulation and signaling pathways involved in the cold-stress response remained unclear in this species. In order to solve the above-mentioned problems, we performed de novo transcriptome assembly and compared the gene expression under the natural (25 °C) and cold (4 °C) conditions for M. wufengensis seedlings. Results More than 46 million high-quality clean reads were produced from six samples (RNA was extracted from the leaves) and were used for performing de novo transcriptome assembly. A total of 59,764 non-redundant unigenes with an average length of 899 bp (N50 = 1,110) were generated. Among these unigenes, 31,038 unigenes exhibited significant sequence similarity to known genes, as determined by BLASTx searches (E-value ≤1.0E-05) against the Nr, SwissProt, String, GO, KEGG, and Cluster of COG databases. Based on a comparative transcriptome analysis, 3,910 unigenes were significantly differentially expressed (false discovery rate [FDR] < 0.05 and |log2FC (CT/CK)| ≥ 1) in the cold-treated samples, and 2,616 and 1,294 unigenes were up- and down-regulated by cold stress, respectively. Analysis of the expression patterns of 16 differentially expressed genes (DEGs) by quantitative real-time RT-PCR (qRT-PCR) confirmed the accuracy of the RNA-Seq results. Gene Ontology and KEGG pathway functional enrichment analyses allowed us to better understand these differentially expressed unigenes. The most significant transcriptomic changes observed under cold stress were related to plant hormone and signal transduction pathways, primary and secondary metabolism, and photosynthesis. In addition, 113 transcription factors, including members of the AP2-EREBP, bHLH, WRKY, MYB, NAC, HSF, and bZIP families, were identified as cold responsive. Conclusion We generated a genome-wide transcript profile of M. wufengensis and a de novo-assembled transcriptome that can be used to analyze genes involved in biological processes. In this study, we provide the first report of transcriptome sequencing of cold-stressed M. wufengensis. Our findings provide important clues not only for understanding the molecular mechanisms of cold stress in plants but also for introducing cold hardiness into M. wufengensis. Electronic supplementary material The online version of this article (10.1186/s12870-019-1933-5) contains supplementary material, which is available to authorized users.


Background
Magnolia is an important landscape tree species around the world. Magnolia, which belongs to the subfamily Magnolioideae of the family Magnoliaceae, is a large genus of approximately 210 flowering plant species. Magnolia wufengensis (Fig. 1) is a new species of Magnolia exhibiting great ornamental and economic value, with the unique tepal characteristics of evenly red color both inside and outside, flaps of a single color and color uniformity [1]. M. wufengensis was first discovered in Wufeng, Hubei Province, southern China [2,3], and has since been introduced to many cities as an excellent landscaping tree. However, compared with most Magnolia species, M. wufengensis is more sensitive to chilling injury. This sensitivity to chilling injury leads to unsustainability and low returns regarding M. wufengensis seedling production and majorly limits its market value in northern China [4]. To ensure normal growth in cold regions, novel M. wufengensis cultivars with improved cold tolerance are needed.
Cold stress is the main form of abiotic stress affecting plant health worldwide and is a major environmental factor limiting plant growth, development, yield and geographical distribution [5]. To overcome this barrier, plants have developed a series of response mechanisms that initiate a battery of events to regulate gene expression, and then physiochemical modifications to enhance their cold resistance [6]. These physiochemical modifications include signaling cascade, modification of the membrane lipid composition, synthesis of cryoprotectant molecules, and an increase in the scavenging activity of reactive oxygen species (ROS) [7][8][9][10].
The changes in gene expression triggered by environmental stress are major molecular response mechanisms of plants adapting to environmental challenges [11]. Numerous studies have demonstrated that cold stress leads to global changes in gene expression, and these changes are highly important for transcriptional regulation in plant adaptation to cold stress [12][13][14][15]. The changes in the expression of many cold-responsive genes involve both regulatory genes and structural genes [16]. Changes in gene expression lead to the accumulation of many protective proteins, such as heat shock proteins (HSPs) [17], cold-regulated (COR) proteins, dehydrins, cryoprotective proteins [6], and various metabolites [18]. Although numerous cold-responsive genes are continually identified and annotated in arrays of plant species, the number of these genes is still insufficient to elucidate the complete cold-responsive mechanism. Furthermore, genes may be involved in different pathways in different species, and different species may exhibit different coldresponsive genes. Therefore, we must continuously explore cold-responsive regulatory pathways and genes in different species.
Genome-wide transcriptome analysis is an efficient method for detecting and elucidating the breadth of molecular mechanisms involved in physiological processes and thus increasing the efficiency of identifying genes of interest [19]. RNA-Seq technology provides robust methods for genome-wide transcriptome analyses and is increasingly used for obtaining a general overview of RNA transcript profiles in both model and non-model species. Recently, cold-responsive transcriptomes have been studied in many crops and other economically important plants, such as Prunus dulcis Mill. [20], Beta vulgaris L. [21], Poncirus trifoliata (L.) Raf. [19], Camellia sinensis [10], and Eucalyptus dunnii [22]. However, no studies using such analytical approaches to obtain a comprehensive overview of the genes involved in cold-responsive molecular mechanisms in M. wufengensis have been performed thus far.
In order to obtain original insight into the molecular mechanisms underlying the cold response of M. wufengensis at the transcriptomic level, we performed transcriptome sequencing and gene expression profiling in normal temperature-treated and cold-treated M. wufengensis using the Illumina sequencing technique. As a result, we identified numerous simple sequence repeats (SSRs) and candidate genes that take part in plant development, primary and secondary metabolism, and hormone and signal transduction pathways. This study provides the first report of the M. wufengensis response to lowtemperature stress. A comprehensive reading of the gene expression changes in leaves under chilling stress can help us understand the cold-stress response pathway at the molecular level and identify efficient ways to enhance cold tolerance in M. wufengensis. In addition, the results of this research provide a basis for better understanding the cold responses of many plants because M. wufengensis is a basal angiosperm.

Results
Response of M. wufengensis to chilling stress M. wufengensis seeds were collected from Wufeng County in Hubei Province, China and germinated in moistened sand. Under normal conditions, 2-month-old M. wufengensis plants exhibit fully expanded mature leaves (Fig. 2a). Following the application of lowtemperature treatment at 4°C, the M. wufengensis plants displayed visible morphological changes. No significant changes were initially observed in the plants for 4 h (Fig. 2b). After the plants were subjected to the cold stress for 8 h, the leaves began to show a slight downward trend (Fig. 2c), and after12 h of the low-temperature treatment, the leaves all faced downward and had softened (Fig. 2d). The whole plants began to show evident phenotypic damage as the duration of the cold treatment was prolonged. After 24 h of the low-temperature treatment, the immature stems exhibited loss of strength, and more severely wilted leaves were observed; furthermore, the edges of the leaves began to dry (Fig. 2e). After 48 h of the chilling stress at 4°C, the leaves had completely withered (Fig. 2f ). These results suggest that M. wufengensis is likely a cold-sensitive plant.  In order to determine the ideal time points for transcriptome sequencing, in addition to observing morphological changes, we measured the physiological response of the plants and the changes in the expression of key genes. Malondialdehyde (MDA), which is a final product of lipid peroxidation, is an important physiological index to measure oxidative damage under abiotic and biotic stress [23]; therefore, we tested the MDA content of the leaves during this experiment. In contrast to the control group, the MDA concentration increased as the duration of the cold treatment increased (Fig. 3a) and was highest (showing a nearly 81% increase) after 12 h of the cold-stress treatment, suggesting that the greatest change in gene expression occurred at this time point. The relative eletrolyte leakage also showed similar changes (Fig. 3c). Previous studies have shown that the tolerance of plants to cold and other stresses is closely correlated with the proline and soluble sugar concentration [24]. With extension of the treatment time, the concentration of proline and soluble sugar increased, showing a nearly fivefold and twofold changes after 12 h; although the proline and soluble sugar concentration slightly decreased later, a high concentration was still maintained ( Fig. 3b and f ). In addition, we also examined the maximum quantum yields of PSII and chloro-phyII content and found that these decreased as the duration of the cold treatment increased (Fig. 3d and e), with the maximum decline observed at 12 h. The C-repeat binding factor (CBF)/dehydration responsive element is regarded as one of the most effective pathways in the Fig. 3 Physiological, biochemical and key gene expression changes in cold-treated Magnolia wufengensis. a MDA contents. b Proline contents. c Relative electrolyte leakage changes. d The maximum quantum yields of PSII. e ChlorophyII content. f Soluble sugar content. g Expression of CBF in M. wufengensis leaves exposed to 4°C for 0 to 48 h. h Expression of SOD in M. wufengensis leaves exposed to 4°C for 0 to 48 h chilling response of plants [25]. Similar to proline, CBF expression reached its highest level at 12 h (Fig. 3g). Superoxide dismutase (SOD) can effectively remove harmful reactive oxygen in plants, which is an important indicator of plant cold resistance [26]. Thus, we examined SOD expression and found that it decreased as the duration of the cold treatment increased (Fig. 3h), with the maximum decline observed at 12 h. Based on these results, we chose the samples that had been cold treated for 12 h for transcriptome sequencing.

Sequencing and de novo transcriptome assembly
Total RNA was extracted from the leaves of 2month-old M. wufengensis from both the control and cold-treated groups (4°C, 12 h). RNA sequencing of 3 biological replicates from the control and coldtreated samples, designated CK1, CK2, CK3, CT1, CT2 and CT3, was then performed on the Illumina HiSeq 4000 platform. The clean reads were obtained after discarding adaptor sequences, low-quality reads (Q-value < 20) and the reads containing more than 10% ambiguous 'N' bases. For the control samples, 86,803,854, 65,817,120, and 61,169,918 clean reads with more than 92% Q30 bases were acquired (Table 1). For the cold-treated samples, 97,721,634, 83,927,960, and 66,717,556 clean reads with more than 93% Q30 bases were acquired ( Table 1). The GC content was approximately 48% in the six samples. The biological replicates produced comparable data ( Table 1).
The filtered clean reads from the six samples were subjected to assembly using the Trinity program [27]. The preliminary assembly generated 94,502 transcripts with an N50 of 1,198 bp. Among these transcripts, the longest transcript was 12,096 bp, and the average transcript length was 916 bp (Table 2). After isoforms were considered, a total of 59,764 unigenes (N50 value =1, 110) were generated. The mean size of these unigenes was 899 bp, and their length ranged from 251 bp to 12, 096 bp ( Table 2). As shown in Fig. 4a, most of the unigenes (55,974, 93.66%) were < 2000 bp in length. In general, the number of unigenes decreased as gene length increased, with lengths between 601 and 800 accounting for the largest proportion (20.41%, 12,198), indicating that the assembly showed high quality.

Annotation and classification of Magnolia wufengensis unigenes
To identify the putative function of the M. wufengensis unigenes, several complementary methods were utilized. The assembled unigenes were searched against the NCBI non-redundant (Nr), SwissProt, String, GO, KEGG, and COG databases using the BLASTX algorithm, with an E-value of less than 1.0 × 10 − 5 . In total, a known protein was found in at least one of the abovementioned databases for 31,038 (51.88%) unigenes (Table 3). Among these unigenes, 30,846 exhibited significant hits (E-value <1e − 5 ) in the Nr database and 20,180 in the GO database. In the other four databases (i.e., SwissProt, KEGG, COG and String), 13,469, 12,180, 7,518 and 13,469 unigenes, respectively, were successfully aligned to known proteins, while 28,756 (48.12%) unigenes did not match any sequences in the six databases. Compared with the other five databases, the Nr database presented the highest proportion of annotations, with 30,846 (99.38%) unigenes being annotated in the Nr database. Overall, the unigene sequences exhibited the most similar BLASTx matches to gene sequences from Vitis vinifera (7,320), followed by Phoenix dactylifera (5,929), Theobroma cacao (1,796), Populus trichocarpa (965) and Jatropha curcas (952) (Fig. 4b).
To identify the active biological pathways in M. wufengensis, pathway annotations of the unigenes was performed using the KEGG pathway tool. The KEGGannotated unigenes (12,180) were distributed to 341 KEGG pathways (Additional file 1: Table S1); among these pathways, the 'plant−pathogen interaction' , 'carbon metabolism' , 'biosynthesis of amino acids' , 'ribosome' and 'starch and sucrose metabolism' pathways were the most abundant (Fig. 4e).
Differentially expressed genes in cold-treated Magnolia wufengensis plants Using the de novo-assembled transcriptome as a reference, the genes expressed in the cold-treated and control groups were identified. We verified consistency among the samples by performing correlation analyses based on FPKM (fragments per kilobase of gene per million mapped reads) values. The high similarity (Fig. 4f, r > 0.91) among the three biological replicates from the control and cold-treated samples demonstrated that the RNA-Seq results were consistent.
To obtain a comprehensive understanding of the transcript expression related to resistance to chilling stress in M. wufengensis, we identified the genes that were differentially expressed between the control M. wufengensis plants and those grown at cold temperatures. Compared with the control plants, 3,910 significantly differentially expressed genes (DEGs) (false discovery rate [FDR] < 0.05 and |log 2 FC (CT/CK)| ≥ 1) were identified in the cold-treated group using edgeR software. Among these DEGs, the number of up-regulated genes was 2,616, and the number of down-regulated genes was 1,294. Functional annotation of these DEGs was performed using the six databases, and approximately 56.75% (2,218) of the DEGs were successfully annotated (Table 3 and Additional file 2:  Table S2).

GO enrichment analysis of DEGs
To obtain insight into the functional categories of the DEGs induced by cold stress, a GO enrichment analysis was performed using Goatools (Fisher exact test, P-value ≤0.05). The most enriched GO category among these DEGs was 'signal transduction' , followed by 'hydrolase activity' , 'lipid metabolic process' , 'oxidoreductase activity' , 'carbohydrate metabolic process' , 'oxidation-reduction process' , 'response to stimulus' , 'single-organism metabolic process' , 'membrane' , and 'chloroplast' (Table 4). Thus, the cold resistance process in M. wufengensis is complex.

Pathway enrichment analysis of DEGs
The pathways that displayed significant changes (P-value ≤0.05) in response to cold stress were identified using the KEGG database. A total of 16 KEGG pathways were significantly enriched (Table 5), among which the 'plant hormone signal transduction' , 'plant-pathogen interaction' , 'starch and sucrose metabolism' , 'mitogen-activated protein kinase (MAPK) signaling pathway' , and 'fatty acid metabolism' pathways were the most highly represented. The 'plant hormone signal transduction' pathway (ko04075) exhibited the most DEGs, suggesting that plant hormones play significant roles in resistance to cold stress in M. wufengensis. The second largest number of DEGs were in the 'plant-pathogen interaction' category (ko04626), indicating that plants are vulnerable to pathogenic bacteria during exposure to cold stress. The 'MAPK signaling pathway' (ko04016) exhibited the fourth largest number of DEGs, indicating that during cold stimulation, the expression of internal genes in plants is regulated by various signaling substances for adaptation to the cold environment. 'Starch and sucrose metabolism' (ko00500) and 'fatty acid metabolism' (ko01212) were the most highly represented, which was unsurprising because cold causes changes in the metabolism of plants to improve their cold tolerance [29].

Validation of RNA-Seq-based DEGs in cold-treated Magnolia wufengensis plants by qRT-PCR
To verify the accuracy and reproducibility of the RNA-Seq data, we analyzed the transcript abundance of the 16 stochastically selected DEGs using qRT-PCR, including 12 up-regulated genes and 4 down-regulated genes in the unigene dataset (Additional file 2: Table  S2). All 12 up-regulated genes were significantly induced by the cold treatment, although their expression changed as the duration of the cold treatment increased (Fig. 5a). Moreover, all four down-regulated   genes were correspondingly repressed by the cold treatment, which was in accord with the RNA-Seq results (Fig. 5b). Therefore, the reliability of our RNA-Seq data was confirmed by the consistency between the qRT-PCR results and RNA-Seq analyses.

Identification of plant hormone and signal transductionrelated unigenes
Signal transduction in plants was evidently influenced by cold stress, as verified by both the GO and KEGG enrichment analyses of the DEGs performed in this study.
Previous studies have shown that kinases play an essential role in signal transduction by relaying signals via protein phosphorylation [30]. In our dataset, more than 170 genes were predicted to encode protein kinases with varying expression levels in M. wufengensis plants under cold-treatment conditions (Additional file 3: Table S3). Genes encoding receptor-like protein kinases (RLKs) accounted for the largest proportion of these genes. More than 67 RLK genes were included in the group of serine/threonine-protein kinases; 21 in the leucine-rich repeat receptor-like protein kinases; eight in the cysteine-rich receptor-like protein kinases; eight in the proline-rich receptor-like protein kinases; and two in the wall-associated receptor kinases (c64483_g1, c66494_g2). Furthermore, the expression of BRASSINOSTEROID-RE-SPONSIVE 1 (BRH1, c41432_g1) decreased, while the expression of an annotated brassinosteroid LRR receptor kinase (c2960_g1) increased under the cold stress. Mitogen-activated protein kinase (MAPK) cascades play important roles in signal transduction induced by abiotic stress [31]. Among the identified DEGs, five MAPKKKs (c64011_g2, c26325_g2, c58822_g1, c20119_g1, c26325_ g1) and one MAPKK (c62421_g1) were up-regulated by the cold treatment, while two MAPKKs (c37197_g1, c46922_g1) were down-regulated. In addition, according to the KEGG analysis of the DEGs, the MAPK signaling pathway (ko04016) was induced under the cold treatment conditions, and the expression of several kinases related to the MAPK signaling pathway was also up-regulated. Calcium, a metal ion, is a secondary messenger in plant cells and plays important roles in many signaling pathways [32,33]. In plants, oscillations in calcium levels can be monitored by calcium sensor proteins, including calmodulin-like proteins (CML), calmodulin-binding protein (CBL) and calcium-dependent protein kinases (CDPKs), which relay signals and trigger downstream responses [34]. In this study, we identified one gene (c55887_g3) encoding CDPKs, four genes (c66505_g1, c42914_g1, c48548_g1, and c65757_g1) encoding CML, eight genes (c59040_g1, c60069_g1, c62530_ g1, c50847_g2, c59040_g2, c50263_g1, c58464_g1, and c46181_g2) encoding CBL and three genes (c53516_g1, c61236_g1, and c61236_g2) encoding CBL-interacting protein kinases (CIPKs) that were up-regulated in the cold-treated group compared with the control group (Additional file 3: Table S3).
According to previous research, reactive oxygen species (ROS) can produce signals in plants triggering cellular changes for resisting changes in the external environment [35]. Plants have an antioxidant system that protects against poisonous oxygen. We identified 36 DEGs encoding enzymes associated with ROS scavenging (Additional file 3: Table S3), including ascorbate peroxidase (APX), peroxidase (POD), glutathione Stransferase (GST), glutathione peroxidase (GPX), polyphenol oxidase (PPO), glutaredoxin, thioredoxin, and polyamine oxidase. Among these antioxidant enzymes, the most important enzymes for ROS scavenging are POD and thioredoxin. ROS promote cellular signaling via an ROS-induced MAPK signaling pathway [31]. In the present study, we identified DEGs that were significantly enriched in the 'MAPK signaling pathway' (ko04016 , Table 5).
Furthermore, the expression of genes involved in hormone biosynthesis or signal transmission pathways was substantially altered following the cold exposure (Table 6). According to the KEGG analysis of the DEGs, we found that the cold treatment could induce carotenoid biosynthesis (ko00906), which is related to ABA biosynthesis in plants. For instance, compared with the control group, the expression of c52591_g1 was up-regulated 3.93-fold in the cold-treated plants; additionally, the gene predicted to be encoded by c52591_g1 is 9-cis-epoxycarotenoid dioxygenase (NCED), an important enzyme in ABA biosynthesis. Interestingly, the expression of a unigene (c62542_g3) encoding a xanthoxin dehydrogenase, which is an important protein in ABA biosynthesis, was down-regulated by the chilling treatment. Furthermore, two unigenes (c48091_g1and c61792_g2) encoding an ABA 8′-hydroxylase, which is an important enzyme in ABA catabolism, were down-regulated by chilling treatment. This apparent regulation of genes involved in ABA biosynthesis and catabolism suggests the presence of dynamic and multifaceted mechanisms controlling ABA content in response to stress. The changes in the expression of genes involved in the ABA signaling pathway were similar to those related to ABA biosynthesis. For instance, three genes encoding regulatory components of the ABA receptor (RCARs; c60506_ g1, c7217_g1, and c38913_g1) were up-regulated, and one RCAR gene (c1769_g1) was down-regulated in response to cold stress. In contrast, one unigene (c62935_g2) encoding PP2C (protein phosphatase 2C), which is a negative regulator of ABA, was down-regulated at 12 h, and three unigenes (c58656_g3, c56328_g1, and c63212_g3) encoding SnRKs (SNF1-related protein kinase) were up-regulated.
The expression of the following four genes, which are associated with the biosynthesis of the defense-type hormone jasmonic acid (JA), was up-regulated: one gene encoding long chain acyl-CoA synthetase (c60616_g1), one gene encoding jasmonate O-methyltransferase (c60697_ g4) and two genes encoding 4-coumarate--CoA ligase-like (c43502_g1 and c41372_g1). Furthermore, two genes (c59607_g4 and c57501_g1) encoding JAZ and two genes (c56341_g2 and c42250_g1) encoding MYC2, which take part in JA signal transduction, were up-regulated under cold-treatment conditions. Auxins generally play an important role in plant apical dominance and lateral bud differentiation [36]. Interestingly, we observed significant changes in the 16 genes involved in the auxin signaling pathway, including one gene encoding an auxin transporter-like protein (AUX), six auxin-induced protein IAA6 (AUX/IAA) genes, four auxin response factors (ARF) and two auxin-responsive GH3-like protein (GH3) genes, all of which were up-regulated in the cold-treated plants. In addition, three unigenes encoding SAUR were identified, including two upregulated unigenes and one down-regulated unigene.

Identification of transcription factors (TFs) involved in the response to cold stress
Transcription factors (TFs) have important functions in regulating upstream cold signal transduction, which can Notably, all the bZIP family members were found to be up-regulated in the present study.

Expression of genes involved in metabolism and biosynthesis
When plants are exposed to cold stress, many of their metabolic processes may exhibit changes. Based on the GO and KEGG pathway enrichment analyses, many DEGs are associated with metabolism and biosynthesis. We found that 28 DEGs (Table 7) were enriched in the 'starch and sucrose metabolism' (ko00500) pathway, including 21 up-regulated DEGs and seven down-regulated DEGs. Under cold-treatment conditions, sucrose synthase, pectin methylesterase, glucosidase and pectinesterase inhibitor expression was induced, while the expression of starch synthase and polygalacturonase decreased. Additionally, cold stress increased the abundance of a pectin lyase-like superfamily protein transcript by 6.41-fold. The 'lipid metabolism process' (GO:0006629) was significantly induced by cold stress. A total of 70 genes associated with lipid metabolism were identified, most of which were related to oxylipin, wax, cutin, suberin, sterol, lipid, and fatty acid biosynthetic and metabolic processes (Additional file 5: Table S5). The genes involved in wax biosynthesis were all evidently up-regulated; for example, the expression level of the gene encoding the ADHESION OF CALYX EDGES protein was 8-fold higher in the cold-treated plants than that in the control plants. A notable, approximately 9-fold increase in the gene encoding GDSL esterase/lipase LTL1 (c51965_g1) was observed after cold treatment.
Furthermore, according to the KEGG analysis, secondary metabolic changes occurred after cold treatment, including changes in the 'phenylpropanoid biosynthesis' (ko00940) and 'flavonoid biosynthesis' pathways (ko00941). In total, 31 DEGs were enriched in the two pathways (Table 8). In the flavonoid biosynthesis pathway, 9 DEGs were identified, and only one gene (c67286_g1) encoding flavonol synthase was down-regulated. In the phenylpropanoid biosynthesis pathway, most up-regulated genes were peroxidases. In addition, dehydrogenase, hydroxylase, methyltransferase and glucosidase were up-regulated. Notably, one caffeoyl-CoA O-methyltransferase (CCoAMT; c61212_ g1) was up-regulated by 6.21-fold under cold-treatment conditions, and this gene is involved in both pathways.
Amino acids are also important for the cold response in plants, as indicated by the 'biosynthesis of amino acids' pathway (ko01230). Twenty DEGs were identified in this pathway (Table 9), and most DEGs showed an up-regulated trend. For example, cold treatment increased the expression of glutamate synthase (c62914_g2 and c65821_ g6), serine hydroxymethyltransferase (c71333_g1 and c77854_g1) and aspartate kinase (c48954_g1). However, certain genes, such as tryptophan synthase (c40985_g1) and cysteine synthase (c64613_g2), were down-regulated.

Cold-responsive genes related to photosynthesis
According to our results, in addition to DEGs related to metabolism, biosynthesis, phytohormones and signal transduction, we also identified many DEGs enriched in the 'photosynthesis' (ko00195) and 'chloroplast' (GO: 0009507) categories. Thus, photosynthesis might be an important process in the cold response of M. wufengensis. The chloroplast is a plant cell-specific organelle that supplies most of the plant energy via photosynthesis. Previous studies have shown that cold stress decreases the rate of photosynthesis by destroying the function of many proteins involved in photosynthesis [38]. We identified 46 genes (Additional file 6: Table S6) encoding components of photosystem I (PSI), PSII and chloroplasts that were significantly decreased under cold-treatment conditions in the M. wufengensis plants. Among these genes, seven genes encoding proteins localized to the thylakoids were identified. For example, LOX2S (linoleate 13S-lipoxygenase 2-1), AOS (allene oxide synthase), and a hypothetical protein (c48930_g1) were all down-regulated under cold-stress conditions. Additionally, psbQ-like protein 3 (PSb; c42229_g1), chlorophyll A/B binding protein (c66781_g2 and c51541_g1), photosystem I P700 chlorophyll an apoprotein (c61894_g1) and phosphoenolpyruvate carboxylase (c51792_g1) all decreased following cold exposure.

Discussion
In this study, we described the transcriptional response caused by cold treatment in M. wufengensis plants. A comprehensive transcriptomic analysis was performed in both cryogenically treated (4°C) and untreated plants through Illumina-based 2 × 150 bp paired-end read sequencing. The Illumina a next-generation sequencing platform for transcriptome assembly generates relatively shorter reads but shows a much higher transcriptome coverage at a lower expense than does other platforms [39]. More than 6 Gb of high-quality clean reads were obtained, which were de novo assembled into 59,764 unigenes with an N50 of 1,110 bp, and the average length of the assembled unigenes was 899 bp, indicating that the assembly was of high quality and likely included many full-length cDNAs. Homologous genes were found for 51.88% of the assembled unigenes in at least one of the public databases that was searched, and 51.57% of the sequences presented a homolog that was determined with a high probability score in one of the following four species: Vitis vinifera (23.58%), Phoenix dactylifera (19.10%), Theobroma cacao (5.78%) and Populus tricho carpa (3.11%). Altogether, these results showed that our M. wufengensis SRA dataset represents an important transcriptome resource for gene discovery and functional analyses. The cold resistance of plants can be improved by inducing or inhibiting gene expression under cold stress [40,41]. In this study, 3,910 DEGs were identified, and 45.6% of these DEGs did not exhibit homologous genes in the examined public databases; these genes may either be unique to M. wufengensis or represent homologs of genes that have yet to be identified in previous studies in other plants.

Signaling mediates cold-stress responses
Under stress, plants can trigger the expression of genes involved in multiple signal transduction pathways and further activate downstream regulatory pathways that are associated with physiological adaptation. First, cold stress affects membrane rigidification in plant cells, which is considered to be the primary event that triggers downstream cold stress responses in plants [42]. In our results, we found that 119 genes associated with membrane components were altered (Additional file 8: Table  S8), and 91 of them were up-regulated, suggesting enhanced membrane fluidity of M. wufengensis under cold stress. The effect of cold on the membrane fluidity and the structural dynamics of membrane proteins initiates cold signaling by affecting ion and metabolite transport [43]. Sixty-one transmembrane transport-related DEGs, including genes associated with organic phosphonate, organic anion, organic acid, amino acid, sulfate and ion transmembrane transport, was identified in cold-treated M. wufengensis seedlings (Additional file 9: Table S9). Especially, the release of Ca 2+ ions from the apoplast, the endomembrane systems, and the organelles has a strong impact on initiating cold acclimation [44]. Intracellular Ca 2+ sensors, such as CIPKs and calciumdependent protein kinases (CDPKs), respond to rising cytosolic Ca 2+ -levels [45]. A total of 16 up-regulated Ca 2+ signaling-related genes were identified under coldtreatment conditions, which is consistent with the expected role of Ca 2+ in early cold signal transduction. Ca 2+ induces intracellular biochemical reactions through the sensing and transduction of downstream CDPKs, thereby regulating the plant response to various abiotic stress signals [46]. It is well known CBL proteins interact with a group of serine/threonine protein kinases known as CIPKs [47]. The CBL-CIPK pathway regulates plant responses to various environmental stresses, including cold, drought and salinity stress [48]. CML proteins usually contain four elongation factor (EF)-hand domains to facilitate Ca 2+ binding [49]. One CDPKs, four CMLs, eight CBLs and three CIPKs were identified as up-regulated in this study, suggesting that the Ca 2+ -mediated signaling pathway has important functions in the M. wufengensis response to cold stress. Besides, DEGs contain more than 170 RLKs (Additional file 3: Table S3), which may play an important role in the perception of cold in M. wufengensis. Most RLKs contain an extracellular domain, a single transmembrane motif, and a cytoplasmic kinase domain. It is believed that the extracellular domains are involved in cold recognition that subsequently leads to activation of the cytoplasmic kinase domain [50]. ROS are normal products of plant cellular metabolism and frequently accumulate when plants suffer from environmental stress. The excessive production of ROS can cause oxidative damage to plant cells, but ROS are also well-known secondary messengers in various cellular processes in plants, including tolerance to abiotic stresses [35,51,52]. The inhibition of photosynthesis is another important consequence of ROS production [53]. In the present study, many DEGs were enriched in the 'photosynthesis' (ko00195) and 'chloroplast' (GO:0009507) categories. Interestingly, the expression of genes in these pathways decreased after cold treatment, which indicated that the photosynthetic capacity of M. wufengensis was inhibited, in turn causing the accumulation of ROS. The delicate balance between ROS production and scavenging determines whether ROS act as damaging or signaling molecules. Excess ROS can be scavenged through an efficient antioxidant system including both nonenzymic and enzymic antioxidants [54]. Genes encoding APX, POD, GST, GPX and PPO were identified in the plants (Additional  file 3: Table S3), and previous studies have shown that these genes play important roles in scavenging ROS under abiotic stress [54,55]. Furthermore, many DEGs were enriched in the carotenoids biosynthesis (ko00906) pathway, studies also have shown that up-regulation of nonenzymic antioxidants, including glutaredoxin, thioredoxin, polyamine and carotenoids, can enhance stress tolerance [56]. The maintenance of a higher antioxidant capacity to scavenge harmful ROS has been shown to be associated with enhanced plant tolerance to environmental stresses [57]. Moderate amount of ROS in M. wufengensis mainly play a signaling role, helping to activate the cellular mechanisms associated with stress tolerance [58]. In addition, ROS signals are transduced to downstream targets through the MAPK cascade [59], and a phosphorylation reaction between the upstream receptor and the downstream target is performed by the MAPKs [60]. We observed that the transcript levels of one MAPKK (c62421_g1) and five MAPKKK (c64011_ g2, c26325_g2, c58822_g1, c20119_g1, and c26325_g1) were up-regulated by cold treatment. Various components of the ROS signaling transduction pathway may be involved in the regulation of cold-stress responses, and our dataset revealed several potential candidates.

Phytohormone signals under cold treatment conditions
Phytohormones play key roles and coordinate various signal transduction pathways during the abiotic-stress response [61]. In this study, according to our transcriptome data, the expression of genes involved in signal transduction pathways associated with several plant hormones was altered after cold treatment for 12 h, suggesting that hormone pathways play critical roles in the M. wufengensis response to cold stress. ABA is considered a necessary messenger in the plant adaptive response to abiotic stresses. In response to environmental stresses, the level of endogenous ABA rapidly increases, which in turn activates specific signaling pathways and modifies gene expression [62,63]. The plant response to cold stress has been suggested to be ABA-independent, although contradictory evidence has emerged [64,65]. NCED, which is a rate-limiting enzyme, is a crucial enzyme in ABA synthesis [66,67]. Accordingly, we found that an NCED gene was up-regulated by 3.93-fold after cold treatment, suggesting that the synthesis of ABA might be elevated in M. wufengensis under cold conditions. Furthermore, we detected increased expression of three RCARs and SnRKs and decreased expression of one PP2C, which are key components in the ABA signaling pathway, as follows: RCARs are proposed to be ABA receptors; the PP2C protein acts as a negative regulator; and SnRKs are positive regulators [68,69]. The presence of ABA leads to the formation of RCAR-PP2C complex, thus inhibiting the activity of PP2C [70,71] and activating SnRKs to act on membrane proteins, ion channels and transcription factors and promote expression of ABA-responsive genes [72]. Thus, ABA-mediated signaling is likely involved in the cold response in M. wufengensis.
Ethylene also plays a crucial role in the adaptation and survival of plants under various stress conditions or is the ultimate effector of other hormone-regulated responses [73]. Ethylene signaling has been implicated in the cold-stress response [74,75], but various studies have shown both positive and negative effects. For instance, many reports have demonstrated the positive role of ethylene in cold stress in rice [76], tomato [74], winter rye [75] and A. thaliana [77]. However, Shi et al. [78] reported that ethylene signaling has a negative effect on cold tolerance by inhibiting the expression of CBF and type-A ARR genes. These contradictory observations may arise from differences in the species studied or the growth of the tissues employed in these studies [77]. In our study, the expression of an ACC oxidase gene (ACO, c57132_g1) was found to be repressed following cold exposure. ACO is linked to the synthesis of ethylene, and its down-regulation may therefore reduce ethylene production. We also identified 8 genes related to ethylenemeditated signaling that showed a trend of down-regulation. Thus, we hypothesize that ethylene may play a negative role in the cold response of M. wufengensis.
The auxins are a group of plant hormones that regulate plant metabolism and are thought to promote plant growth [79]. Interestingly, we found significant changes in the 16 genes involved in the auxin signaling pathway. Although knowledge regarding the role of auxins in abiotic stress is limited, studies are increasingly investigating the role of auxins in cold resistance. One study showed that cold stress primarily targets the transport of intracellular auxin and cannot completely turn off all the auxin protein transport pathways [80]. SIZ1, which is a key regulator of the cold signaling pathway, stabilizes ICE1 under cold conditions by inhibiting its polyubiquitination [81] and has been found to have a negative effect on the phosphate-starvation-induced remodeling of root architecture by controlling auxin patterning [82]. Therefore, whether and how auxin-regulated signaling affects the cold responses in M. wufengensis remains to be fully established.
In addition to ABA, ethylene and auxin, other hormones may also affect the plant response to cold stress. The effects of the growth-related hormone GA on abiotic stress and its interaction with other hormones under environmental stress have become increasingly clear [83]. For example, GID1 is a GA receptor that determines GA signaling transduction [84]. The expression of a unigene (c47742_g1) encoding GID1 was observed to be down-regulated in this study, suggesting that GA signaling transduction may be inhibited by cold stress. Thus, cold-stressed plants may conserve energy by slowing their growth in response to adverse conditions. Altogether, different hormonal signals may affect the expression and regulation of many downstream genes.

Transcription factors involved in the cold-stress response
Transcription factors (TFs) play a critical role in plant development and stress tolerance [40]. These proteins regulate the expression of inducible genes by interacting with cis-acting elements of the promoter region of many stress-inducible genes [85]. Furthermore, TFs highlight the function of diverse transcriptional mechanisms in the stress signaling pathways and are associated with the regulation of various abiotic stress signals and the transcription of various genes [86]. Numerous TF families, including the AP2/EREBP, bHLH, WRKY, MYB, NAC and MYC families, are thought to orchestrate the signals that are transduced when plants are exposed to diverse environmental stresses [87]. According to our results, 113 predicted TFs were induced by cold treatment. The largest group of cold-inducible TFs belonged to the AP2-EREBP, and bHLH family, which was composed of 20 members, respectively. Among these TFs, two unigenes (c64671_g1and c47639_g1) were annotated as CBF/DREB genes that have important functions in plant cold tolerance [16]. Therefore, the CBF pathway is conserved in the M. wufengensis response to cold stress. Many studies have demonstrated that MYB TFs are involved in various cold-stress regulatory pathways, especially in the regulation of CBF by the ABA-independent signaling pathway [25]. The NAC family may be involved in the regulation of ABA-dependent signaling pathway genes under diverse stresses, including cold stress [88]. Other TFs, including the bHLH, WRKY, HSF, bZIP, MYB-related, and GRAS families, were also found to be induced by cold stress. Studies have confirmed that members of these families can significantly affect the cold resistance of plants [5], but further study will be required to determine whether these TFs function independently or synergistically to enhance cold resistance in M. wufengensis.

Cold-stress signaling induces cellular protection processes
The complex interactions among diverse signaling pathways and multiple transcription factors collectively affect the physiological responses to cold stress, including the regulation of membranes, cellular redox status, and primary and secondary metabolism [6,25]. M. wufengensis cells rapidly accumulate compatible solutes and osmoprotectants, including sugars, amino acids, fatty acid et al. These osmolytes help to re-establish cellular osmotic balance by increasing the water potential of cells and thereby play a key role in the protection and stabilization of cellular organelles, proteins and membranes [89]. Some studies have confirmed that there is a strong correlation between sugar concentration and cold resistance, and sugar accumulation acts synergistically with cold-inducible gene products [90]. Fructan, sucrose, fructose, and raffinose likely stabilize the membrane system through osmotic regulation and act as ROS scavengers to protect plants from oxidative damage in low-temperature environments [14]. During cold acclimation, sugar level increases under artificial conditions, which could be due to improved transport mediated by sucrose synthase or invertase [91]. These observations are consistent with our results, as greater induction of sugar biosynthesis-related genes was observed in the plants under cold stress. Notably, two genes encoding trehalose synthase (c62697_g1and c64134_g3) were up-regulated, which is consistent with resports by Li et al. [92], who found that plants can counteract abiotic stresses by accumulating trehalose for oxidative detoxification. We also found that the expression of starch synthase and polygalacturonase decreased, likely because most starch breaks down into soluble sugars under cold stress.
Altered lipid biosynthesis, biomembrane rearrangement, and specific fatty acid changes reduce the impairment of cell membranes by cold stress. A group of lipid metabolism-related genes, including genes associated with sterol, cutin, oxylipin, wax, suberin, and fatty acid biosynthesis and metabolism, was identified in coldtreated M. wufengensis plants. Wax may play critical roles in regulating the cold resistance of M. wufengensis, as evidenced by the high abundance of transcripts of genes related to its biosynthesis. Previous studies have demonstrated that wax biosynthesis and deposition are involved in regulating abiotic stresses tolerance [12]. A gene encoding GDSL esterase/lipase, which contributes to the protection of plants against abiotic stresses [93], was found to be strongly induced under cold stress. This result provides powerful evidence indicating that lipid metabolism in M. wufengensis may be crucial in the adaption to cold stress.
A positive correlation has been shown between the accumulation of amino acids and cold tolerance, and amino acids including proline, glutamic acid, glutamine, alanine, aspartic acid, serine and tryptophan have been reported to accumulate in response to low temperatures [94,95]. Additionally, cold stress can induce osmotic stress, and amino acids act as important osmolytes [96]. In our study, we identified 20 amino acid biosynthesisrelated genes that contribute to the protection of M.
wufengensis plants against cold injury by regulating ion transport and stomatal opening and affecting enzyme activity, gene expression, and redox homeostasis [96]. Thus, the synthesis of amino acids might be promoted to help plants combat cold stress.
In addition to activating the primary mechanisms of sugar accumulation, amino acid biosynthesis and lipid metabolism, secondary mechanisms are induced in plants leading to the biosynthesis of other secondary metabolites, most of which are derived from the phenylpropanoid biosynthesis and flavonoid biosynthesis pathways. Flavonoids prevent oxidative damage under abiotic stress by removing excess ROS [97]. Phenolic compounds, which are the main products of phenylpropanoid biosynthesis, also play various physiological roles in abiotic stress [98]. In this study, 9 flavonoid biosynthesisrelated genes were detected, which is similar to the findings reported by Li et al. [41], who demonstrated that the accumulation of flavonoids protects petunias from chilling injury. Among related enzymes in the phenylpropanoid biosynthetic pathway, phenylalanine ammonia-lyase is one of the most relevant [99]. However, we identified only one gene encoding this particular enzyme, whereas most of the identified genes encoded peroxidases, which oxidize phenols and show increased activity in response to various types of stress [100]. Most flavonoids are superior to wellknown antioxidants, such as ascorbate, but various flavonoid oxidation processes could also be involved in ROS scavenging [101]. Therefore, peroxidase might significantly modulate these phenomena [102].

Photosynthesis is inhibited by cold stress
Photosynthesis is the basic physiological process underlying plant growth and development, but low temperatures can significantly affect this process [103]. In our research, we found that photosynthesis was inhibited by cold stress, which inhibited the expression of genes involved in chlorophyll biosynthesis and destroyed the function of the thylakoids. The decline in photosynthesis under cold stress is mainly due to cold stress induced inhibition of the activity of key enzymes associated with photosynthesis [38]. Two photosynthetic enzyme-related genes, encoding linoleate 13S-lipoxygenase 2-1 (LOX2S) and allene oxide synthase (AOS), were observed to be down-regulated, and PSII was inhibited, as indicated by the low abundance of transcripts of chlorophyll A/B binding protein and psbQ-like protein (Psb) genes. In Elymus nutans, cold stress can inhibit the expression of photosynthesis-related genes [104]. Furthermore, the inhibition of photosynthesis reduced assimilation on the photosynthetic electron transport chain, which cause redox imblances, support generation of ROS [53], and damage the photosystems [105]. The inhibition of photosynthesis may lead to the accumulation of ROS, resulting in ROS damage, which is likely a factor in the poor cold tolerance of M. wufengensis.

Conclusions
M. wufengensis is a new Magnolia species with important economic and ornamental value, but it is highly sensitive to low temperatures. In the present study, we comprehensively described the transcriptional responses of M. wufengensis leaves to cold stress. Overall, 3,910 DEGs were obtained using RNA-Seq data; these genes included 2,616 up-regulated DEGs and 1,294 down-regulated DEGs. The identified genes were classified by performing functional annotation through BLASTx, GO, and KEGG analyses. Based on the transcriptomic data, a hypothetical model of the response to cold stress in M. wufengensis is provided in Fig. 6. Cold stimulation may first be sensed through cold-induced membrane rigidification, most likely via transmembrane proteins, which in turn leads to cold signal transduction, including Ca 2+ , phytohormone and ROS signaling, ultimately inducing changes in cellular metabolism. All signaling pathways act on protein kinases that can switch on various downstream transcription factors (TFs) family proteins. The activation of TFs triggers the expression of a cascade of downstream cold-responsive genes, whose transcription further regulates the homeostasis of cellular metabolism, such as lipid metabolism or secondary metabolism. Furthermore, excess ROS can be scavenged via photosynthesis and cellular redox processes, which further affects ROS signaling. The model established by comparing the differences in gene expression between cold-treated and control plants may promote further studies addressing the molecular mechanisms underlying the cold response in plants. In addition, a large number of the identified genes (approximately 45.65% of all 3,910 DEGs) encode proteins of unknown functions. Studies on these genes might reveal new mechanisms by which M. wufengensis responds to cold stress.

Plant materials and stress treatments
The M. wufengensis seeds used in this study were collected from Wufeng County in Hubei Province, China. Before the experiment, the seeds should be cleaned and stored at 4°C in paper bags. First, 0.1% (w/v) sodium hypochlorite were used to sterilize the M. wufengensis seeds for 3 min. The seeds were then washed 5 times with sterile distilled water and germinated in moistened sand for 20 d at 25°C. Next, the seedlings showing uniform growth were transferred to a solid substrate with a mixture of nutritive soil, perlite and vermiculite in a ratio of 5:1:1(v/v), with one seedling per pot. The seedlings were kept in a growth chamber under a 16/8 h (day/night) photoperiod at 25°C and 65% relative humidity. The cold treatment was performed in a climate chamber under the same photoperiod indicated above at 4°C. The control group plants were planted under a normal temperature (25°C), and other conditions were unchanged. When the seedlings were 60 days old, they were exposed to low temperatures. Leaves were collected from the plants after cold treatment for 0 h (control), 4, 8, 12, 24, and 48 h, then immediately frozen in liquid nitrogen and stored at − 80°C for further analyses. Collection was performed from more than 3 plants at each sampling time, and each sample collection was repeated 3 times to obtain biological replicates.

RNA extraction
The total RNA was extracted from the leaves of both the cold-treated and control groups using an EASYspin plant RNA Extraction Kit (Aidlab China) following the manufacturer's instructions. High-quality RNA samples (OD260/280 = 1.8~2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, 28S:18S ≥ 1.0, > 10 μg), assessed using a 2100 Bioanalyzer (Agilent) and quantified with an ND-2000 (NanoDrop Technologies), were employed to construct the sequencing library after the removal of genomic DNA using DNase I (TaKara). Each experiment included three biological replicates.

RNA-Seq library construction and Illumina HiSeq4000 sequencing
The RNA-Seq library was constructed using a TruSeq™ RNA sample preparation kit from Illumina (San Diego, CA). In brief, polyA mRNA was isolated from 5 μg of total RNA using oligo (dT) beads and then fragmented into small pieces with fragmentation buffer. Using these small fragments as templates, double-stranded cDNAs were synthesized with a SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA) and random hexamer primers (Illumina). The double-stranded cDNAs were then end-repaired and phosphorylated, and a polyA tail was added following Illumina's library construction protocol. Next, target fragments of 200-300 bp were selected via agarose gel electrophoresis and amplified using Phusion DNA polymerase (NEB) over 15 PCR cycles, and the products were used to construct the RNA-Seq library. The final library was quantified by TBS380 fluorometer, and sequencing was carried out on the Illumina HiSeq 4000 platform (2 × 150 bp read length). Each experiment included three biological replicates.

De novo assembly and sequence annotation
The adaptor reads, low-quality reads (Q-value < 20) and the reads containing more than 10% ambiguous 'N' bases were deemed poor-quality reads and discarded, after which the remaining were considered clean reads. Following preprocessing, the clean reads from the control and cold-treated (4°C for 12 h) 60-day-old plant samples were employed for de novo assembly using Trinity (http://trinityrnaseq.github.io) [106]. The contigs were clustered and further assembled according to paired-end information. The non-terminal extended sequences of contigs identified by Trinity were assembled, and the resulting sequences were defined as unigenes. All the assembled unigenes were searched against public protein databases using BLASTX with an E-value of less than 1.0 × 10 − 5 ; the databases included the NCBI nr (NCBI non-redundant protein database), String, SwissProt, COG (Clusters of Orthologous Groups) and KEGG (Kyoto Encyclopedia of Genes and Genomes) databases. The Gene Ontology (GO) annotations of the assembled unigenes were obtained using the BLAST2GO (https://www.blast2go.com/) [28] program with nr annotation. Metabolic pathway analysis of the assembled unigenes was performed according to the KEGG (http:// www.genome.jp/kegg/) database [107].

Analysis and functional enrichment of DEGs
The FPKM (fragments per kilobase of exon per million mapped fragments) method was used to determine the expression level of each transcript. Gene and isoform abundances were calculated with RSEM (http://dewey lab.biostat.wisc.edu/rsem/) [108]. Differential expression was analyzed based on the count of the expression level of each transcript in the control and coldtreated samples using EdgeR (Empirical analysis of Digital Gene Expression in R, http://www.bioconduc tor.org/packages/2.12/bioc/html/edgeR.html) software [109]. A false discovery rate [FDR] < 0.05 and |log 2 FC (CT/CK)| ≥ 1 were used as criteria for identifying significant differences in expression. GO enrichment analyses of differentially expressed genes (DEGs) were performed using Goatools (https://github.com/tanghai bao/Goatools) and KOBAS (http://kobas.cbi.pku.edu. cn/expression.php) [110] with GO annotation. Next, we employed the same method for KEGG pathway functional enrichment analysis of DEGs. The DEGs were considered to be significantly enriched when a Bonferroni-corrected P-value ≤0.05 was obtained.

qRT-PCR analysis of gene expression
Total RNA was extracted from the leaves of 60-dayold plants treated at 4°C for 0 (control), 4, 8, 12, 24, and 48 h using an EASYspin plant RNA Extraction Kit (Aidlab China). The RNA was then employed to synthesize first-strand cDNAs using an oligo (dT) 18 primer, after the removal of genomic DNA using DNase I (TaKaRa, Japan). qRT-PCR was performed with the StepOnePlus™ Real-Time PCR System (ABI) using SYBR premix Ex Taq (Takara, Japan). Each experiment included three biological replicates. qRT-PCR was performed under the following conditions: 95°C for 3 min; forty cycles of 95°C for 15 s, and 58°C for 30 s. The DNA primers are listed in Additional file 7: Table S7. The relative expression of each gene was calculated according to the comparative cycle threshold (ct) method, and M. wufengensis ACTION was employed as a reference.

Physiological analyses of cold-treated M. wufengensis plants
The physiological changes in M. wufengensis under the applied cold treatment conditions were analyzed before transcriptome sequencing. The MDA content was determined according to the thiobarbituric acid (TBA) reaction, as proposed by Dhindsa et al. [111] with minor modifications. The proline concentrations were measured using the sulfosalicylic acid-acid ninhydrin method proposed by Bates et al. [112] with slight modifications. Control and cold-treated leaves were also sampled to measure relative eletrolyte leakage [113], soluble sugar concentration [114], maximum quantum yields of PSII [115] and chlorophyII content [116]. Each experiment included three biological replicates.