Skip to main content

Physiological and transcriptome analysis of Magnolia denudata leaf buds during long-term cold acclimation



Magonlia denudata is an important perennial tree species of the Magnoliaceae family, known for its ornamental value, resistance to smoke pollution and wind, role in air purification, and robust cold tolerance. In this study, a high-throughput transcriptome analysis of leaf buds was performed, and gene expression following artificial acclimation 22 °C, 4 °C and 0 °C, was compared by RNA sequencing.


Over 426 million clean reads were produced from three libraries (22 °C, 4 °C and 0 °C). A total of 74,503 non-redundant unigenes were generated, with an average length of 1173.7 bp (N50 = 1548). Based on transcriptional results, 357 and 235 unigenes were identified as being upregulated and downregulated under cold stress conditions, respectively. Differentially expressed genes were annotated using Gene Ontology and the Kyoto Encyclopedia of Genes and Genomes pathway analyses. The transcriptomic analysis focused on carbon metabolism and plant hormone signal transduction associated with cold acclimation. Transcription factors such as those in the basic helix-loop-helix and AP2/ERF families were found to play an important role in M. denudata cold acclimation.


M. denudata exhibits responses to non-freezing cold temperature (4 °C) to increase its cold tolerance. Cold resistance was further strengthened with cold acclimation under freezing conditions (0 °C). Cold tolerance genes, and cold signaling transcriptional pathways, and potential functional key components for the regulation of the cold response were identified in M. denudata. These results provide a basis for further studies, and the verification of key genes involved in cold acclimation responses in M. denudata lays a foundation for developing breeding programs for Magnoliaceae species.


Magnolia denudata, an important species in the Magnoliaceae family, is well-known for its ornamental value, resistance to smoke pollution and wind, and contribution to air purification. As a robust deciduous species thriving under climatic conditions, its origins are attributed to the southeastern region of North America, and is now cultivated in China, in areas south of the Yangtze River and in northern cities such as Beijing and Lanzhou in cold areas [1]. M. denudata can survive low or freezing temperatures conditions and regrow in the spring. This species demonstrates high cold tolerance and can be cultivated in northern China [2]. To understand the mechanisms by which M. denudata adapts to cold stress, it is important to investigate its cold-resistance mechanism.

Cold is considered a major abiotic stress factor that tends to negatively affect growth, development, reproduction, yield, and survival [3]. To reduce the adverse effects of cold stress, temperate -plants have evolved various adaptive mechanisms including cold acclimation a process in which plants improve freezing resistance after exposure to short days and low non-freezing temperatures for a specific period [4]. Cold acclimation is an important strategy for plants to survive cold winters and involves physiological, metabolic, and cellular changes occurring after sensing cold temperatures [5, 6]. Responses include cold signal transmission, alterations in cellular structure and lipid composition, increase in antioxidant levels, overexpression of cold-related transcription factors (TFs) and accumulation of osmoregulatory substances such as soluble sugar and protein [7,8,9,10,11,12]. Cold tolerance, especially freezing tolerance, can be markedly enhanced after cold acclimation [13].

Transcriptional regulation of gene expression plays an important role in plant adaptation and cold stress tolerance [14]; the C-repeat binding factor (CBF)/dehydration responsive element (DRE) plays a crucial role the transcriptional pathways related to cold tolerance. Cold-induced CBF is one of the members of the APETALA2/−ETHYLENE RESPONSE FACTOR (AP2/ERF) gene family, and can induce the expression of cold response (COR) genes by binding to cis-acting elements in their promoters [12]. Many genes belonging to COR, such as dehydrin, encode a wide range of cryoprotective proteins that confer protection to cells against cold-induced damage [15].

RNA sequencing (RNA-seq) has been widely used in the analysis of the tolerance response in plants such as Elymus nutans [16], Camellia oleifera [17], Citrus paradisi [18], and Nicotiana tabacum [19]. Sequencing M. denudata improves our understanding of the physical changes involved in cold tolerance, and provides an in-depth explanation of its remarkable cold-hardiness.

Therefore, it is imperative to further study the gene expression profile of M. denudata under cold acclimation conditions. However, few studies have investigated the metabolic pathways and TFs related to cold resistance mechanisms in this species. This study aims to bridge our understanding of the molecular mechanisms associated with phenotypic and physiological changes during artificial cold acclimation.


Physiological responses of M. denudata to cold acclimation

To determine the ideal time for RNA-seq, we measured the physiological responses of M. denudata besides observing phenotypic changes. We chose five temperatures (22 °C, 12 °C, 4 °C, 0 °C, and − 4 °C) represented as N1 (control), CA1, CA2, CA3, and CA4, respectively, to observe the changes in M. denudata with decreasing temperature.

The median lethal temperature (LT50), obtained by determining the relative conductivity can reflect a plant’s cold resistance level [20]. Thus, we determined the LT50 and found that it was the highest at N1 (10.41 °C), followed by CA1 (7.59 °C), CA2 (− 5.68 °C), CA3 (− 12.83 °C), and CA4 (− 11.50 °C). Based on these results, cold acclimation is the strongest at 0 °C (Fig. 1A).

Fig. 1
figure 1

Physiological changes in leaf buds of M. denudata at different cold conditions. A MDA content, LT50, B soluble protein, soluble sugar content, C proline and ABA content in M. denudata at different temperatures. Lowercase letters indicate significant differences among different temperature treatments at p < 0.05

Malondialdehyde (MDA), as the final product of lipid peroxidation, is an important indicator for measuring oxidative damage under different stresses [21]. In contrast to LT50, the MDA concentration increased with the enhancement of the cold acclimation. At CA1, the MDA concentration of the shoot apex (buds) increased slightly and then increased significantly at CA2, followed by a gradual increase at CA3. The highest MDA concentration (nearly a two-fold increase) was obtained at CA4. The rate of increase in the MDA content of the buds was 23.76, 46.82, 65.95, and 82.51%, respectively, at CA1, CA2, CA3, and CA4. Based on the above values, the increase rate reached a maximum at − 4 °C (Fig. 1A).

Previous studies have shown that cold-resistance and other abiotic stresses in plants are associated with soluble sugar and soluble protein concentration [22]. With the extension of cold acclimation, the concentrations of soluble protein and soluble sugar increased from N1 to CA3, showing nearly a 4-fold and 2-fold increase, respectively, at CA3. However, the concentration decreased at CA4, with a value lower than that at CA2 (Fig. 1B). Hence, the maximum value was observed at CA3.

Similar to the above trend, the contents of proline and abscisic acid (ABA) showed slight increases at CA1, followed by sharp increases at CA2. Finally, the peak appeared at CA3 with a 7-fold and 25-fold increase for proline and abscisic acid, respectively. There was a steep decrease at CA4. Based on the above values, the maximum was observed at CA3 (Fig. 1C).

According to previous studies, 22 °C is the most suitable temperature for the growth of M. denudata; therefore in our study we set the control at 22 °C. In addition, 12 °C, which is a relatively warm temperature, did not show significant changes compared with the control. Based on our physiological data, at CA4, soluble sugar and protein concentrations decreased after 0 °C. At this time, we estimated that M. denudata had exceeded the limit of cold acclimation and suffered from freezing damage. In summary, we selected 22 °C (N1), 4 °C (CA2), and 0 °C (CA3) for the analysis of the transcriptome data to determine the molecular mechanisms underlying the cold acclimation of M. denudata.

Transcriptome sequencing and de novo assembly

For the transcriptome assembly, a cDNA sample was prepared using high-quality RNA extracted from stem apexes (buds) for three libraries (4 °C [CA2], 0 °C [CA3], and control [N1]), and used for sequencing. Approximately 50 million reads were obtained for the control sample (N1), and 49 million raw reads were obtained for the two cold-treated samples (CA2 and CA3), respectively. The Q30 percentage (sequencing error rate < 1%) and the guanine and cytosine (GC) percentage were approximately 95 and 48% in the nine libraries, respectively (Table 1).

Table 1 Sequencing the M. denudata transcriptome in nine leaf buds samples from plants that controlled (N11, N12, N13) or cold-treated (CA21, CA22, CA23 and CA31, CA32, CA33)

The de novo assembly yielded 74,503 unigenes with lengths ranging from 301 to 23,158 bp, and the average length was 1173.7 bp (Table 2). In general, the number of unigenes decreased with the increase of gene length after 400 bp; unigenes with a length of 401–500 bp (8658; 11.62%) accounted the largest proportion, followed by length from 501 bp to 600 bp (8263; 11.09%) (Fig. 2A).

Table 2 Statistics of transcriptome unigenes
Fig. 2
figure 2

Characteristics and annotations of unigenes. A Unigene lengths distribution in M. denudata B Top 10 Species distribution for each unigene based on NR database C GO categorization functional annotation of the unigenes D KEGG annotation of the unigenes E Distribution of unigenes with KOG functional classification annotation

Annotation and classification of M. denudata unigenes

The assembled unigenes were searched against Nr, SwissProt, KEGG, KOG, eggNOG, GO, and Pfam. databases using BLASTx with an E-value < 10− 5. Among these unigenes, there were 31,803 and 19,666 exhibited significant differences (E-value <1e− 5) in the Nr and GO database, respectively. In addition, 22,385, 9403, 18,935,28,452, and 40 unigenes, respectively, were successfully annotated to known proteins in the other five databases (SwissProt, KEGG, KOG, eggNOG, and Pfam) (Table 3). In general, the unigene sequences were most similar to genes from N. nucifera (8235), followed by E. guineensis (4304), P. dactylifera (3331), V. vinifera (3311), and M. acuminata subsp. malaccensis (2327) via BLASTx matches (Fig. 2B).

Table 3 Annotation statistics of M. denudata unigenes

Gene annotation and functional classification of M. denudata unigenes

Gene Ontology (GO) was used to annotate the predicted M. denudata genes. According to the GO analysis, 19,666 unigenes were classified into the following three ontologies: biological process, cellular component, and molecular function. In the biological process category, genes involved in the ‘cellular process’ (12,765) and ‘metabolic process’ (11,346) were well represented. The cellular component category mainly comprised ‘cell’ (14,762), ‘cell parts’ (14,740), and ‘organelles’ (11,868). Within the molecular function category, ‘binding’ (12,358) and ‘catalytic activity’ (10,710) were highly represented, and there was a considerable number of genes pertaining to ‘biological regulation’ (5401), ‘membrane’ (6303), and ‘organelle part’ (6105) categories; however, few pertained to ‘cell killing’ (10), ‘virion’ (8), ‘virion part’ (6), and translation regulator activity (8) (Fig. 2C). The GO analysis identified a substantial number of genes related to biological processes and cellular components following cold acclimation.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a signal pathway database used for analyzing the metabolic pathways, utilities of biological systems and functions of genes. Based on a comparison against the KEGG database, 9403 unigenes (12.62%) were found to be significantly enriched in the database and were assigned to 214 KEGG pathways. The most represented pathways were ‘Environmental Information Processing−−Signal transduction’ (965 members), ‘Metabolism−−Carbohydrate metabolism’ (840 members), and ‘Genetic Information Processing−−Translation’ (783 members) (Fig. 2D).

In total, 18,935 unigenes were assigned to a KOG classification and divided into 25 specific categories (Fig. 2E). The largest group was ‘general function prediction only’ (7750), followed by ‘signal transduction mechanisms’ (1589), ‘posttranslational modification’ and ‘protein turnover, chaperones’ (1530). Only a few unigenes were assigned to the ‘nuclear structure’ (60), and ‘cell motility’ (57) ‘cell motility’ (11) categories.

Gene changes among the different phases of cold acclimation

Differentially expressed genes (DEGs) with a p-value < 0.05 and |log2FC| > 1 met the differential expression threshold. Among the DEGs identified, expression levels of 12,666 and 8824 genes were upregulated, and those of 9363 and 8523 genes were downregulated in the CA2 and CA3 groups, respectively. Between CA3 and CA2, 5355 up-regulated and 6840 down-regulated genes were identified; there were more DEGs in CA2 vs N1 than those observed in the other groups, while the number of DEGs in CA3 vs CA2 was the lowest of the three (Fig. 3A). Compared to N1, expression levels of 7633 (39.7%) genes were induced exclusively by CA2, while for CA3 only 1981 (10.3%) unigenes were specific. A similar trend in DEGs was observed in the repressed unigene groups with the development of cold acclimation; there were 5054 (27.5%) unigenes specifically expressed at CA2, while there were only 2305 (12.6%) unigenes expressed at CA3. Among these DEGs, 357 upregulated and 235 downregulated genes were detected in response to all treatments (Fig. 3B and C).

Fig. 3
figure 3

The expression of the gene changes in M. denudata during different cold acclimation stages. A Changes in gene expression profile during different cold acclimation stages. The numbers of up-regulated (red) and down-regulated (blue) genes between CA2 vs. N1, CA3 vs. N1 and CA3 vs. CA2 are summarized B Venn diagrams of CA3 vs. CA2, CA3 vs. N1 and CA2 vs. N1 differential genes expression C Venn diagrams of CA3 vs. CA2 and CA2 vs. N1 differential genes expression

Overall, in M. denudata, transcript abundance changed significantly from N1 to CA3; the expression of cold response genes could be induced and expressed largely from N1 to CA2, although there were fewer changes in gene expression changes at CA3, these may be crucial to cold acclimation. At the same time, we found that there were less DEGs compared to CA2 vs. N1, while there were several DEGs showing expression in response to cold acclimation at CA3 which indicated the role CA3 played in cold acclimation.

DEGs of assembled M. denudata transcripts under cold acclimation conditions

Compared to N1, GO analysis of DEGs at CA2 demonstrated that genes related to membrane structure and transcription were overexpressed. In biological process, cellular component, and molecular function, ‘transcription, DNA-templated’, ‘integral component of membrane’, and ‘DNA binding transcription factor activity’ were the most enriched GO categories. Compared to N1, for biological process, the major subcategory was ‘DNA integration’. In the cellular component, ‘plasma membrane’, and ‘extracellular region’ were the most representative subcategories, ‘nucleic acid binding’, ‘zinc ion binding’, and ‘aspartic-type endopeptidase activity’ were the top three subcategories in the ‘cellular component’ category in the development of cold acclimation at CA3 (Additional file 1: Table S1).

Compared to N1, KEGG pathway enrichment analysis for DEGs indicated that three pathways, namely – ‘phenylpropanoid biosynthesis (ko00940)’, ‘biosynthesis of amino acids (ko01230)’ and ‘starch and sucrose metabolism (ko00500)’ were significantly enriched, and pathways associated with ‘plant hormone signal transduction (ko04075)’ (Fig. 4A), ‘starch and sucrose metabolism (ko00500)’ and ‘biosynthesis of amino acids (ko01230)’ were significantly decreased at CA2. At CA3 vs. N1, ‘carbon metabolism (ko01200)’, ‘protein processing in endoplasmic reticulum (ko04141)’, ‘starch and sucrose metabolism (ko00500)’, and ‘glycolysis / gluconeogenesis (ko00010)’ pathways were significantly enriched, and the ‘starch and sucrose metabolism (ko00500)’ and ‘plant hormone signal transduction (ko04075)’ (Fig. 4B) pathways were significantly depleted. Compared with CA2, ‘plant hormone signal transduction (ko04075)’ and ‘biosynthesis of amino acids (ko01230)’ pathways were enriched, ‘glycerophospholipid metabolism (ko00564)’, ‘glycerolipid metabolism (ko00561)’, and ‘DNA replication (ko03030)’ pathways were significantly depleted at CA3 (Additional file 2: Table S2).

Fig. 4
figure 4

KEGG map about plant hormone signal transduction compared with N1. Red indicates significantly upregulated; green indicates significantly downregulated; yellow indicates proteins encoded by both up- and downregulated genes; purple indicates no significantly changed. A CA2 vs. N1; B CA3 vs. N1

Response of important TFs and protein kinases (PKs) to cold acclimation

TF-regulated gene expression networks have an important effect on plant growth and development. To explore the TFs in response to cold acclimation in M. denudata, 261 TF genes (103 upregulated and 158 downregulated) were identified. These TFs were mainly concentrated in the basic helix-loop-helix (bHLH), ERF, WRKY, bZIP, GATA, MYB, zinc finger and AP2 families. Among all the genes evaluated, the number of genes for bHLH was the largest, with 52 genes (10 upregulated and 42 downregulated), followed by ERF with 27 genes (17 upregulated and 10 downregulated), WRKY with 22 genes (15 upregulated and 7 downregulated), bZIP with 17 genes (6 upregulated and 11 downregulated), GATA with 13 genes (6 upregulated and 7 downregulated), MYB with 13 genes (4 upregulated and 9 downregulated), zinc finger with 3 genes (1 upregulated and 2 downregulated), and AP2 with 3 genes (1 upregulated and 2 downregulated). The expression levels of the same gene was different for different periods of cold acclimation. The expression of CL6131Contig1 (the GT-3b gene) increased significantly at CA2 and then decreased significantly with the extension of treatment time (Additional file 3: Table S3).

A total of 423 DEGs encoding PKs were identified in the cold stress gene sets. There were 222 DEGs (86 upregulated and 136 downregulated) encoding serine/threonine-PKs (STPKs), including 60 DEGs encoding LRR receptor-like STPKs (18 upregulated and 42 downregulated), 33 DEGs encoding G-type lectin S-receptor-like STPKs (19 upregulated and 14 downregulated) and 7 DEGs encoding CBL-interacting STPK (5 upregulated and 2 downregulated). In addition, there were 16 DEGs (4 upregulated and 12 downregulated) encoding calcium-dependent PKs (CDPKs), 4 DEGs (3 upregulated and 1 downregulated) belonging to the CBL-interacting protein kinases (CIPKs) family, 7 DEGs (6 upregulated and 1 downregulated) belonging to the mitogen-activated protein kinase kinases (MAP 3Ks) family. We also found 91 receptor-like protein kinase (RLK) DEGs (36 upregulated and 49 downregulated) included in the group of STPKs, including 27 leucine-rich repeat RLKs, 14 cysteine-rich RLKs and 7 proline-rich RLKs. The expression of CL37311Contig1 and CL57668Contig1, which are inactive leucine-rich RLK and rust resistance kinase, increased at CA2 and then decreased at CA3 (Additional file 4: Table S4).

Gene expression of biosynthesis and metabolism related to proline and hormones

Proline, as one of the most effective osmotic balance adjustment substances, accumulates rapidly when plants subjected to various abiotic stresses; 54 genes that are involved in the regulation of and relief from osmotic stress were involved in ‘arginine and proline metabolism (ko00330)’ based on the KEGG enrichment analysis (Additional file 7: Fig. S1). Expression levels of four genes (CL8828Contig1, CL12870Contig1, CL38215Contig1, and CL53606Contig1) for delta1-pyrroline-5-carboxylate synthetase (P5CS) were significantly upregulated. The expression of CL8828Contig1, CL12870Contig1, and CL53606Contig1 increased by 5-, 6.2-, and 12.9-fold, respectively, at CA2, and CL38215Contig1 increased by 3.2-fold at CA3 compared with CA2. DREB1, a key TF involved in proline synthesis, is highly expressed during cold acclimation and promotes proline synthesis for cold tolerance.

Expression of transportation inhibitor response 1 (TIR1) gene in the auxin signaling pathway was downregulated in CA2 and CA3 compared to the control. The auxin response factor and small auxin up RNA (SAUR) levels also decreased after CA2.

ABA plays a key role in the development of cold tolerance in plant cells. Since protein phosphatase 2C (PP2C) is a negative regulator of the ABA signaling pathway, pyrabactin resistance like (PYL) abscisic acid receptors bound to ABA can initiate ABA signaling by inhibiting PP2C activity. For example, expression levels of three DEGs (CL32501Contig1, CL32603Contig1, and comp14208_c0_seq1_2) encoding ABA receptor PYL family, were upregulated after CA2.

Ethylene, one of the five major types of plant hormones, plays a vital role in the development of resistance against the effect of cold. The expression levels of ethylene-insentive3 (EIN3)-binding F-box protein1/2 (EBF1/2) and constitutive triple response1 (CTR1), which are located downstream of the ethylene receptor, increased at CA2 (Additional file 5: Table S5).

Validation of RNA-Seq by qRT-PCR

To validate the reliability of RNA-seq, 12 DEGs were selected and tested by qRT-PCR. The genes belonged to different categories including cold-responsive protein, sugar transport, PKs, metabolism, TFs and hormone-related signaling. The trends of genes under different treatments conditions observed using qRT-PCR were similar to those observed using RNA-seq, whereas two DEGs, the gene encoding light-harvesting complex II chlorophyll a/b binding protein 1 (LHCB) (CL1043Contig2) and auxin-responsive protein IAA27 (CL3681Contig1) displayed a slightly inconsistent expression pattern (Fig. 5), confirming the reliability of RNA-seq.

Fig. 5
figure 5

The expression of 12 DEGs from 22 °C(N1), 4 °C (CA2), 0 °C (CA3) and comparison between qRT-PCR and RNA-seq. The relative expression of the qRT-PCR assay is calculated by the 2−ΔΔCT method and are displayed as blue columns; UBQ was used as the reference gene. FPKM values of RNA-seq are shown as red polylines


Signaling mediates cold-stress responses

Cell membranes counter effects of low temperature with membrane rigidification after plants are subjected to cold stress, and subsequently the low-temperature signal is transmitted downstream through a secondary messenger such as Ca2+ or ABA. This triggers a Ca2+ influx, which acts upstream of the genes related to CBFs and COR [7, 12]. In this study, GO enrichment helped identify an important role of the ‘integral component of membrane’ in with the development of tolerance to cold stress, suggesting that membrane components were altered with the enhancement of the membrane fluidity. Under cold stress, the saturation of membrane lipid decreases with the upregulating content of fatty acid desaturase (FAD) which increases membrane fluidity, thereby increasing the cold tolerance [23, 24]. In our study, one FAD5 gene (CL30097Contig1) was identified whose expression downregulated at CA2 and upregulated at CA3 indicating a possible reason for the enhancement of cold resistance at 0 °C. Additionally, 31 DEGs related to calcium ion binding (15 upregulated and 16 downregulated) were identified, and only three of these DEGs were involved in the ABA-activated signaling pathway [25, 26]. In a previous study, Ca2+-dependent protein kinases (CDPKs) were verified as important sensors of Ca2+ flux in plants in response to different abiotic stresses [27], and have been shown to enhance abiotic stress tolerance [28]. In this study, 12 genes related to CDPK were identified. Notably, seven genes were downregulated from N1 to CA2 with the enhancement of cold resistance, indicating that genes related to CDPK were not completely functional at this juncture. Moreover, the expression of CL31785Contig1, a gene encoding CDPK, was first downregulated at CA2 and then upregulated at CA3. For M. denudata, the low temperature at CA2 did not completely stimulate the cold stress response, whereas extremely low temperatures such as those at CA3, could further strengthen cold tolerance. These results support the findings of this study, as soluble sugar content increased significantly at CA2 (cold) and CA3 (extreme cold conditions), and the stage from CA2 to CA3 was deemed for cold acclimation. However, in tomato (Lycopersicon esculentum) [29] and rice (Oryza sativa. L) [30], overexpression of CDPK induced the expression of genes related to the development of cold tolerance under cold stress conditions. This is not consistent with the results of this study, indicating a different function of CDPK in M. denudata.

Light signaling is also closely related to the exhibition of cold stress responses [31, 32]. Phytochromes are proved involving in cold tolerance by regulating CBF/DREB1s expression [33, 34]. Phytochrome interacting factor (PIFs), belonging to the bHLH family [35], can directly bind to phytochromes and transcriptionally regulate the expression level of the downstream genes of PHYA/B [36]. Several studies have revealed that PIF1 [33], PIF3 [37], and PIF4 [32] exert negative effects on cold tolerance. In this study, PIF1 (CL38596Contig1), PIF3 (CL39360Contig1), and PIF4 (CL38753Contig1) were identified, which were found to be downregulated at CA2 and CA3 compared with the control. Notably, expression of CL38753Contig1, which encoded PIF4, was first downregulated and then upregulated at CA2. Overall, temperatures below 4 °C had a negative effect on M. denudata.

Hormone signaling and transport in response to cold acclimation

Plant hormones, especially ABA and auxin, participate in signal transduction pathways during the response to abiotic stress [38]. ABA is an important signaling molecule in the plant adaptive response to cold stress [39]. In response to cold conditions, ABA-responsive element binding factor (ABF) can activate the expression of downstream genes [40]. Five ABF genes (four upregulated and one downregulated) were identified, and this finding was consistent with the increase in ABA content from N1 to CA3. Apart from the genes related to signaling, ABA-IMPORTING TRANSPORTER 1 (AIT1), an ABA importer gene, involved in ABA biosynthesis [41] was also identified. In this study, two AIT1 genes were identified, and expression levels of both were upregulated at CA3 indicating that ABA transport occurred actively under cold acclimation conditions. This finding was consistent with the variation in ABA content.

Studies suggest that auxin may lead to the rapid expression of certain genes related to cold tolerance and auxin metabolism is also part of the plant cold-resistance regulation system [42]. Numerous genes in the ‘auxin-activated signaling pathway’ and ‘response to auxin’ were detected in this study. Previous studies have shown that the protein serine/threonine kinase PINOID (PID) [43], an early auxin-inducible gene, regulates PIN localization on cellular membranes and thus regulates polar auxin transport [43, 44]. Currently, studies on PID genes mainly focus on many plants such as Arabidopsis [45], tobacco [46], cucumber [47], and rice [48], whereas studies on woody perennials have rarely been reported because most researchers select leaves or whole seedlings for investigation in experiments. In this study, shoot apexes (leaf bud) were used as they are considered the main existing organ for their vital role in survival during natural overwintering, and two PID genes (both with downregulated expression levels) were identified. The results confirmed that auxin plays an important role in cold response mediation in M. denudata and will thus be the focus of our future research.

Carbohydrate metabolism in response to cold acclimation

Soluble sugars, constituting an important class of carbohydrates, can mitigate the damage caused by cold stress to plants by decreasing the freezing point of the cell, thus preventing membrane injury, and increasing the stability of the cell structure [49,50,51]. Cold acclimation promotes the development of cold resistance, and plants with robust cold resistance mechanisms tend to have higher soluble sugar content [52, 53]. In this study, soluble sugar content increased with an increase in cold stimulus. The rate of increase from N1 to CA2 was faster than that from CA2 to CA3, which indicated that soluble sugar and soluble protein accumulated more actively from N1 to CA2, to prevent the infliction of damage attributable to cold conditions. Meanwhile at CA3, the contents of soluble sugar and soluble protein remained at higher levels supporting the fact that M. denudata persisted at the cold-acclimation stage at CA3.

In addition to phenotypic data, the KEGG analysis helped identify 307 regulated genes in the ‘carbon metabolism’ pathway, 250 regulated genes in the ‘starch and sucrose metabolism’ pathway, 222 regulated genes in the ‘protein processing in endoplasmic reticulum’ pathway, and 141 regulated genes in the ‘amino sugar and nucleotide sugar metabolism’ pathway (Table 4 and Additional file 8: Fig. S2). This demonstrates that carbohydrate metabolism plays an important role in cold acclimation.

Table 4 Categorization of M. denudata unigenes to KEGG biochemical pathways

Sucrose synthase (SuS) and sucrose phosphate synthase (SPS) are close related to sucrose metabolism and cold tolerance development [54, 55]. Expression of SuS and SPS is induced to exhibit responses to cold stress via regulation of sucrose metabolism [56]. In Brassica oleracea L. [57] and Triticum aestivum L. [58], the expression of SuS is induced by cold stress. In the transcriptome of M. denudata, seven genes (three upregulated and four downregulated) encoded SuS and one gene (upregulated) encoded SPS; most genes were expressed from N1 to CA2. Expression of four genes related to SuS were downregulated by low temperature. These results were not consistent with those reported by previous studies that explained the mechanisms by which SuS and SPS were involved in sucrose metabolism, and elucidated the manner in which sucrose metabolism influenced cold acclimation responses. These discrepancies warrant further study.

Similarly, trehalose and raffinose are important in sugar metabolism during cold acclimation [59, 60], and eight trehalose-phosphate phosphatases (four upregulated and four downregulated) and five raffinose synthases (three upregulated and two downregulated) were identified. These genes did not increase one-sidedly, which indicated that they played a limited role in the development of cold resistance in M. denudata. These synthases were found to mostly function at the cold acclimation stage from N1 to CA2; therefore, sugar metabolism might play a major role in the earlier stage of cold acclimation. Previous studies have identified Sugars Will Eventually Be Exported Transporters (SWEETs) as being involved in sugar transport and output [61,62,63]. In addition to the above-mentioned results, seven SWEET genes we also identified, and of these five SWEET genes were downregulated under cold acclimation conditions, including one SWEET16 gene. SWEET16 can catalyze the transport glucose and sucrose, and its expression is downregulated by low temperature in Arabidopsis and tea (Camellia sinensis). This indicated that under cold acclimation conditions, less amount of sugar was exported and more amount of sugar remained in the plants, and that carbohydrate metabolism contributed to cold tolerance development, a result consistent with those reported by previous studies.

Cold-related TFs expressed in response to cold acclimation

To develop adaptation and resistance mechanisms to different stresses, plants correspondingly form self-suited regulatory systems, in which transcriptional regulation plays an important role. TFs play a key role in the environmental stress response network by binding to cis-acting elements in the promoter region and by regulating the expression of a series of genes [64].

AP2/ERF is a plant-specific gene family including four subfamilies, namely AP2, RAV, DREB and ERF. Among the AP2/ERF family, DREB, and ERF are closely related to exhibition of the plant abiotic stress response [65]. In this study, six CBF DEGs and 27 ERF DEGs were identified, which might serve as candidate genes for further study on the cold response mechanisms in M. denudata. In the DREB1 subgroup, DREB1A/CBF3, DREB1B/CBF1, and DREB1C/CBF2 can be induced under cold stress and the overexpression of DREB1s/CBFs significantly improves tolerance under abiotic stress [66,67,68]. PdCBF2 increased in close correlation with cold acclimation in buds of almond (P. dulcis) [69]. In the six CBF DEGs, expression levels of four genes (two CBF3, one CBF4 and one CBF1) were continuously upregulated during cold acclimation, while expression levels of two CBF1 genes (CL64792Contig1 and CL63168Contig1) were downregulated at CA2 and then upregulated at CA3. In Arabidopsis, overexpression of AtCBF1 induces COR genes and enhances cold tolerance [70]. These results further prove that CA3 is an important stage for the strengthening of cold tolerance. Whether the two CBF1 genes in M. denudata possess specific functions warrants further study.

The bHLH gene family [71] plays an important role in the regulation of plant development and stress response [72]. In this study, 52 bHLH (10 upregulated and 42 downregulated) were identified. In a previous study, ICE1 (inducer of CBF expression), a bHLH transcription factor, could specifically encode a MYC-like bHLH transcriptional activator that bound to the promoters of CBF3 at low temperatures to induce the expression of CBF3 [73, 74]. ICE1 has been shown to positively regulate cold tolerance [74]. In this study, two genes encoding ICE (CL32817Contig1 and CL35094Contig1) did not upregulate after cold acclimation according to a previous study in Arabidopsis [74, 75], while the expression of CBF3 was upregulated. Therefore, understanding the role of ICE1 in the development of cold resistance and the causes of CBF3 upregulation in M. denudata merits further study.

In this study, AP2/ERF and bHLH families were identified as two highly active groups in response to cold acclimation, indicating that M. denudata might improve cold tolerance through signaling pathways and transcription regulated by AP2/ERF and bHLH protein. However, further study is warranted to understand the relationships among TFs to enhance cold resistance mechanisms in M. denudata.

Antioxidation mechanisms in response to cold acclimation

Reactive oxygen species (ROS) are rapidly accumulated when plants counter different abiotic stresses. ROS can act as a signal switch in response to environmental stress but can also cause oxidative damage to plants [76]. In this study, the MDA content increased with the extension of cold acclimation, indicating the occurrence of ROS accumulation. Superoxide dismutase (SOD), ascorbate peroxidase (APX), and catalase (CAT) are crucial enzymes that counter ROS, and play an important role in conferring antioxidant protection to plants [77, 78]. Their significant upregulation can effectively reduce oxidative damage inflicted by ROS and enhance low-temperature tolerance in plants. The expression increased with decreased temperature. Under cold stress conditions, SOD, CAT, and APX genes in M. denudata were significantly and differentially expressed, and the up-regulated genes played an important role. The expression levels of CL27100Contig1 (SOD) and CL38755Contig1 (SOD) were up-regulated, whereas those of CL11817Contig1 (CAT) and CL34986Contig1 (APX) were downregulated. This resulted in no significant difference in the MDA content between CA2 and CA3.


In this study, using de novo sequencing, a M. denudata dataset containing 74,503 unigenes was constructed to observe the dynamic changes of gene expression under different cold acclimation stages. A considerable number of DEGs were identified under cold acclimation conditions, especially those involved in cold signal perception and transmission, hormone regulation, TFs related to cold regulation and adaptation, and those involved in the elimination of ROS in antioxidation mechanisms. Based on the transcriptomic data, a responsive diagram of low temperature and SD in M. denudata has been hypothetically summarized in Fig. 6.

Fig. 6
figure 6

A hypothetical responsive diagram in M. denudata buds during cold acclimation

Transcriptome differences were compared based on different temperatures; comparisons with the control (N1 vs. CA3, N1 vs. CA3, and CA2 vs. CA3) helped identify the genes responsive to cold stress. Among the three pairwise comparisons, there were 357 upregulated and 235 downregulated transcripts between the cold-stress and control samples. Based on annotation analysis, pathways especially proline and hormones-related involve in cold acclimation and cold tolerance enhancement in M. denudata. From N1 to CA2, M. denudata experienced the most important cold acclimation stage, and cold tolerance was further strengthened at CA3. The transcriptome expression profiling of M. denudata may facilitate our understanding of the molecular regulation mechanism related to cold acclimation in M. denudata. This study also provides insights to improve cold resistance in species that are sensitive to cold injury in the Magnoliaceae family.


Plant materials and cold acclimation treatments

M. denudata seeds were obtained from a single tree and after conducting sterilization and seeding, the seedlings of M.denudata were cultured in single pots in growth chambers with a temperature of 25 ± 1.0 °C (day) / 22 ± 1.0 °C (night), a photoperiod of 12-h light/12-h dark and 65% relative humidity in February 2017. The three apical buds of a six-month-old seedling per pot were selected randomly and used for experiments during the whole cold acclimation treatment in September 2017 (autumn). The experiment included five temperatures for analysis. A room temperature of 22 °C in a low-temperature incubator (3 M, USA) served as the control (N1). The samples in other groups of M. denudata buds were treated for seven days sequentially in low-temperatures incubators at the following four different experimental temperatures: low temperature of 12 °C (CA1), 4 °C (CA2), 0 °C (CA3), and − 4 °C (CA4). Two days were left for cooling slowly in temperatures incubator between two groups. Additionally, three groups (CA2, CA3, and N1) of living M. denudata buds were subjected to seven days of cold acclimation treatment for RNA-seq. The long-term cold acclimation (22 °C to − 4 °C) of M. denudata was also compared using physiological and biochemistry experiments (Fig. 7).

Fig. 7
figure 7

Experimental design of cold acclimation in M. denudata. A room temperature of 22 °C was served as the control (N1) and then M. denudata buds were treated at low temperature of 12 °C (CA1), 4 °C (CA2), 0 °C (CA3), and − 4 °C (CA4) for seven sequential days. Two days were left for cooling slowly in temperatures incubator between two groups

Assay of physiological index contents

Physiological index contents of the leaf buds were determined by performing the following methods. The LT50 and proline content were determined according to the method described by Duan et al. and Bates et al. [79, 80]. The ABA content was estimated by performing high-performance liquid chromatography according to the method described by Xiao et al. [81]. Soluble sugar and protein contents were determined using the method described by Yang et al. [82]. The MDA content was assessed according to the thiobarbituric acid method. Each experiment included three biological and technical replicates and each biological replicate included three plants.

RNA isolation

Total RNA isolation from the leaf buds of N1, CA2, and CA3 groups was performed using the mirVana™ miRNA Isolation Kit (Ambion-1561) according to the manufacturer’s instructions. RNA samples were assessed using the Agilent Bioanalyzer 2100 (Agilent) and were quantified using the NanoDrop-2000 (Thermo Scientific). After removing the genomic DNA using DNase I (Takara), the high-quality RNA samples were used to construct the sequencing library. Each experiment included three biological replicates.

cDNA library construction and sequencing

Total RNA was extracted using the mirVana miRNA Isolation Kit (Ambion) following the manufacturer’s instructions. RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent). The libraries were constructed using TruSeq Stranded mRNA LTSample Prep Kit (Illumina) according to the manufacturer’s protocol. Then these libraries were sequenced on the Illumina sequencing platform (HiSeqTM 2500) and 125 bp/150 bp paired-end reads were generated. Each experiment included three biological replicates.

De novo assembly and sequence annotation

Raw data (raw reads) were processed using Trimmomatic to control the quality and remove the connector [83]. Reads containing ploy-N and low quality reads were removed to obtain clean reads. After removing the adaptor and low quality sequences, the clean high-quality reads were assembled into expressed sequence tag clusters (contigs) and assembled into transcripts de novo using Trinity [84] (vesion: trinityrnaseq_r20131110) in paired-end method. The longest transcript was selected as a unigene according to the sequence similarity and length for subsequent analysis.

The function of the unigenes was annotated by conducting alignment of the unigenes with the NCBI non-redundant (NR), SwissProt, and Clusters of orthologous groups for eukaryotic complete genomes (KOG) databases using BLASTx [85] with a threshold E-value of 10− 5.The proteins with the highest hits to the unigenes were used to assign functional annotations. Based on the annotation in SwissProt, Swissprot ID was mapped to the GO term to obtain the GO annotation of the unigene. The unigenes were mapped to the KEGG database [86] to annotate their potential metabolic pathways.

Analysis and functional enrichment of DEGs

FPKM [87] and read counts of each unigene were calculated using bowtie2 [88] and eXpress [89]. Unigenes were standardized by estimateSizeFactors. P-values and foldchange of unigenes based on the analysis of different groups were calculated by nbinomTest from DESeq [90]. The data on DEGs in each treatment group were then filtered based on the threshold p-value < 0.05 and |log2FoldChange| > 2. GO enrichment and KEGG pathway enrichment analysis of DEGs were conducted using R, according to the hypergeometric distribution.

TFs analysis

TFs were predicted based on CDS predicted protein sequences. Plant TFs domains were obtained from HMM search ( and the unigenes were classified according to gene family information.

qRT-PCR analysis of gene expression

To validate the accuracy of the RNA-seq results, 12 DEGs were selected for qRT-PCR validation and UBQ [91] was used as a reference gene for analysis. The purified RNA samples were reverse-transcribed to cDNA using 2 μg of RNA in a 20 μl reaction system with the 5X All-In-One RT MasterMix with AccuRT (ABM, USA). The specific qRT-PCR primers were designed using Beacon Designer 7. qRT-PCR was conducted using the StepOnePlus™ Real-Time PCR System (Applied Biosystems, USA) and the final reaction contained 5 μl TB Green Premix Ex Taq (Tli RNaseH Plus) (2X), 0.2 μl ROX Reference Dye (50X) (Takara, Japan), 0.2 μl upstream primer, 0.2 μl downstream primer, 1.0 μL cDNA, and 3.4 μl ddH2O in a final volume of 10 μl. The primer sequences used in qRT-PCR are listed in Additional file 6 (Table S6). Each sample included three biological and technical replicates.

Statistical analysis

The data were analyzed by one-way ANOVA followed by least significant difference test (p value < 0.05) using SPSS 22.0 software (IBM Corporation, USA). Graphs were constructed by SigmaPlot version 10 (Systat Software, USA) and R Project (R Foundation for Statistical Computing, Austria).

Availability of data and materials

The data sets used and analyzed during the study are available from the in the SRA database of NCBI (



Abscisic Acid




ascorbate peroxidase




Basic leucine zipper




C-repeat binding factor


Calcineurin B-like protein


Calcium-dependent protein kinases


CBL-interacting protein kinases


Differentially expressed genes


The evolutionary genealogy of genes: Nonsupervised Orthologous Groups database


Fragments per kb per million fragments


Gene ontology


inducer of CBF expression


Median Lethal Temperature


Kyoto encyclopedia of genes and genomes


Eukaryotic ortholog group




the NCBI non-redundant


The Protein families database


Real-time quantitative reverse transcription PCR


Reactive oxygen species


Transcription factors


Superoxide dismutase


  1. Xu F, Kirchoff BK. Pollen morphology and ultrastructure of selected species of Magnoliaceae. REV PALAEOBOT PALYNO. 2008;150(1):140–53.

    Article  Google Scholar 

  2. Yang Y, Jia Z, Chen F, Sang Z, Ma L. Comparative analysis of natural cold acclimation and deacclimation of two Magnolia species with different winter hardiness. Acta Physiol Plant. 2015;37(7):129.

    Article  CAS  Google Scholar 

  3. Chinnusamy V, Zhu J, Sunkar R. Gene Regulation During Cold Stress Acclimation in Plants. Methods in molecular biology (Clifton, NJ). 2010;639:39–55.

    Article  CAS  Google Scholar 

  4. Weiser CJ. Cold resistance and injury in Woody plants: knowledge of hardy plant adaptations to freezing stress may help us to reduce winter damage. SCIENCE. 1970;169(3952):1269–78.

    Article  CAS  PubMed  Google Scholar 

  5. Iba K. ACCLIMATIVE RESPONSE TO TEMPERATURE STRESS IN HIGHER PLANTS: approaches of gene engineering for temperature tolerance. Annu Rev Plant Biol. 2002;53(1):225–45.

    Article  CAS  PubMed  Google Scholar 

  6. Atıcı Ö, Nalbantoǧlu B. Antifreeze proteins in higher plants. PHYTOCHEMISTRY. 2003;64(7):1187–96.

    Article  CAS  PubMed  Google Scholar 

  7. Reddy ASN, Ali GS, Celesnik H, Day IS. Coping with stresses: roles of calcium- and calcium/calmodulin-regulated gene expression. Plant Cell. 2012;23(6):2010–32.

    Article  Google Scholar 

  8. Mittler R. Oxidative stress, antioxidants and stress tolerance. Trends Plant Sci. 2002;7(9):405–10.

    Article  CAS  PubMed  Google Scholar 

  9. Gusta LV, Wisniewski M, Nesbitt NT, Gusta ML. The effect of water, sugars, and proteins on the pattern of ice nucleation and propagation in acclimated and nonacclimated canola leaves. Plant Physiol. 2004;135(3):1642–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Van Buskirk HA, Thomashow MF. Arabidopsis transcription factors regulating cold acclimation. PHYSIOL PLANTARUM. 2006;126(1):72–80.

    Article  Google Scholar 

  11. Baek K, Skinner DZ. Alteration of antioxidant enzyme gene expression during cold acclimation of near-isogenic wheat lines. Plant Sci. 2003;165(6):1221–7.

    Article  CAS  Google Scholar 

  12. Chinnusamy V, Zhu J, Zhu J. Cold stress regulation of gene expression in plants. Trends Plant Sci. 2007;12(10):444–51.

    Article  CAS  PubMed  Google Scholar 

  13. Hincha DK, Zuther E. Introduction: plant cold acclimation and freezing tolerance. Methods Mol Biol. 2014;1166:1–6.

    Article  CAS  PubMed  Google Scholar 

  14. Lee B, Henderson DA, Zhu J. The Arabidopsis cold-responsive transcriptome and its regulation by ICE1. Plant Cell. 2005;17(11):3155–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Thomashow MF. PLANT COLD ACCLIMATION: freezing tolerance genes and regulatory mechanisms. Annu Rev Plant Physiol Plant Mol Biol. 1999;50(1):571–99.

    Article  CAS  PubMed  Google Scholar 

  16. Fu J, Miao Y, Shao L, Hu T, Yang P. De novo transcriptome sequencing and gene expression profiling of Elymus nutans under cold stress. BMC Genomics. 2016;17(1):870.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Wu L, Li JA, Li Z, Zhang F, Tan X. Transcriptomic analyses of Camellia oleifera ‘Huaxin’ leaf reveal candidate genes related to long-term cold stress. Int J Mol Sci. 2020;21(3):846.

    Article  CAS  PubMed Central  Google Scholar 

  18. MAUL P, McCOLLUM GT, POPP M, GUY CL, PORAT R. Transcriptome profiling of grapefruit flavedo following exposure to low temperature and conditioning treatments uncovers principal molecular components involved in chilling tolerance and susceptibility. Plant Cell Environ. 2008;31(6):752–68.

    Article  CAS  PubMed  Google Scholar 

  19. Hu R, Zhu X, Xiang S, Zhan Y, Zhu M, Yin H, et al. Comparative transcriptome analysis revealed the genotype specific cold response mechanism in tobacco. BIOCHEM BIOPH RES CO. 2016;469(3):535–41.

    Article  CAS  Google Scholar 

  20. AIRAKI M, LETERRIER M, MATEOS RM, VALDERRAMA R, CHAKI M, BARROSO JB, et al. Metabolism of reactive oxygen species and reactive nitrogen species in pepper (Capsicum annuum L.) plants under low temperature stress. Plant Cell Environ. 2012;35(2):281–95.

    Article  CAS  PubMed  Google Scholar 

  21. Jena K, Kim S, Suh J, Yang C, Kim YG. Identification of cold-tolerant breeding lines by quantitative trait loci associated with cold tolerance in Rice. Crop Sci. 2012;52(2):517–23.

    Article  Google Scholar 

  22. Moellering E, Muthan B, Benning C. Freezing tolerance in plants requires lipid remodeling at the outer chloroplast membrane. SCIENCE. 2010;330(6001):226–8.

    Article  CAS  PubMed  Google Scholar 

  23. Murakami Y, Tsuyama M, Kobayashi Y, Kodama H, Iba K. Trienoic fatty acids and plant tolerance of high temperature. SCIENCE. 2000;287(5452):476–9.

    Article  CAS  PubMed  Google Scholar 

  24. Kodama H, Horiguchi G, Nishiuchi T, Nishimura M, Iba K. Fatty acid desaturation during chilling acclimation is one of the factors involved in conferring low-temperature tolerance to young tobacco leaves. Plant Physiol. 1995;107(4):1177–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Kudla J, Batistič O, Hashimoto K. Calcium Signals: the Lead currency of plant information processing. Plant Cell. 2012;22(3):541–63.

    Article  Google Scholar 

  26. Hetherington AM, Brownlee C. THE GENERATION OF Ca2+ SIGNALS IN PLANTS. Annu Rev Plant Biol. 2004;55(1):401–27.

    Article  CAS  PubMed  Google Scholar 

  27. Boudsocq M, Sheen J. CDPKs in immune and stress signaling. Trends Plant Sci. 2013;18(1):30–40.

    Article  CAS  PubMed  Google Scholar 

  28. Cheng S, Willmann MR, Chen H, Sheen J. Calcium signaling through protein kinases. The Arabidopsis calcium-dependent protein kinase gene family. Plant Physiol. 2002;129(2):469–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Leclercq J, Ranty B, Sanchez-Ballesta M, Li Z, Jones B, Jauneau A, et al. Molecular and biochemical characterization of LeCRK1, a ripening-associated tomato CDPK-related kinase. J Exp Bot. 2005;56(409):25–35.

    Article  CAS  PubMed  Google Scholar 

  30. Yang T, Poovaiah BW. Calcium/calmodulin-mediated signal network in plants. Trends Plant Sci. 2003;8(10):505–12.

    Article  CAS  PubMed  Google Scholar 

  31. Ding Y, Shi Y, Yang S. Molecular regulation of plant responses to environmental temperatures. Mol Plant. 2020;13(4):544–64.

    Article  CAS  PubMed  Google Scholar 

  32. Lee CM, Thomashow MF. Photoperiodic regulation of the C-repeat binding factor (CBF) cold acclimation pathway and freezing tolerance in Arabidopsis thaliana. Proc Natl Acad Sci. 2012;109(37):15054–9.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Franklin KA, Whitelam GC. Light-quality regulation of freezing tolerance in Arabidopsis thaliana. Nat Genet. 2007;39(11):1410–3.

    Article  CAS  PubMed  Google Scholar 

  34. Wang F, Guo Z, Li H, Wang M, Onac E, Zhou J, et al. Phytochrome a and B function antagonistically to regulate cold tolerance via abscisic acid-dependent Jasmonate signaling. Plant Physiol. 2015;170(1):459–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Toledo-Ortiz G, Huq E, Quail PH. The Arabidopsis basic/Helix-loop-Helix transcription factor family. Plant Cell. 2003;15(8):1749–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Castillon A, Shen H, Huq E. Phytochrome interacting factors: central players in phytochrome-mediated light signaling networks. Trends Plant Sci. 2007;12(11):514–21.

    Article  CAS  PubMed  Google Scholar 

  37. Jiang B, Shi Y, Zhang X, Xin X, Qi L, Guo H, et al. PIF3 is a negative regulator of the CBF pathway and freezing tolerance in Arabidopsis. Proc Natl Acad Sci. 2017;114(32):E6695–702.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Peleg Z, Blumwald E. Hormone balance and abiotic stress tolerance in crop plants. Curr Opin Plant Biol. 2011;14(3):290–5.

    Article  CAS  PubMed  Google Scholar 

  39. Sah SK, Reddy KR, Li J. Abscisic acid and abiotic stress tolerance in crop plants. Front Plant Sci. 2016;7:571.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Choi H, Hong J, Ha J, Kang J, Kim SY. ABFs, a family of ABA-responsive element binding factors. J Biol Chem. 2000;275(3):1723–30.

    Article  CAS  PubMed  Google Scholar 

  41. Kanno Y, Hanada A, Chiba Y, Ichikawa T, Nakazawa M, Matsui M, et al. Identification of an abscisic acid transporter by functional screening using the receptor complex as a sensor. Proc Natl Acad Sci. 2012;109(24):9653–8.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Eremina M, Rozhon W, Poppenberger B. Hormonal control of cold stress responses in plants. Cell Mol Life Sci. 2016;73(4):797–810.

    Article  CAS  PubMed  Google Scholar 

  43. Friml J. A PINOID-dependent binary switch in apical-basal PIN polar targeting directs auxin efflux. SCIENCE. 2004;306(5697):862–5.

    Article  CAS  PubMed  Google Scholar 

  44. Okada K, Ueda J, Komaki MK, Bell CJ, Shimura Y. Requirement of the auxin polar transport system in early stages of Arabidopsis floral bud formation. Plant Cell. 1991;3(7):677–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Lin F, Xu D, Jiang Y, Chen H, Fan L, Holm M, et al. Phosphorylation and negative regulation of CONSTITUTIVELY PHOTOMORPHOGENIC 1 by PINOID in Arabidopsis. Proc Natl Acad Sci. 2017;114(25):6617–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Lee SH, Cho H. PINOID positively regulates auxin efflux in Arabidopsis root hair cells and tobacco cells. Plant Cell. 2006;18(7):1604–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Song M, Cheng F, Wang J, Wei Q, Fu W, Yu X, et al. A leaf shape mutant provides insight into PINOID serine/threonine kinase function in cucumber (Cucumis sativus L.). J Integr Plant Biol. 2019;61(9):1000–14.

    Article  CAS  PubMed  Google Scholar 

  48. Wu HM, Xie DJ, Tang ZS, Shi DQ, Yang WC. PINOID regulates floral organ development by modulating auxin transport and interacts with MADS16 in rice. Plant Biotechnol J. 2020;18(8):1778–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Krasensky J, Jonak C. Drought, salt, and temperature stress-induced metabolic rearrangements and regulatory networks. J Exp Bot. 2012;63(4):1593–608.

    Article  CAS  PubMed  Google Scholar 

  50. Kasuga J, Arakawa K, Fujikawa S. High accumulation of soluble sugars in deep supercooling Japanese white birch xylem parenchyma cells. NEW PHYTOL. 2007;174(3):569–79.

    Article  CAS  PubMed  Google Scholar 

  51. Morin X, Ameglio T, Ahas R, Kurz-Besson C, Lanta V, Lebourgeois F, et al. Variation in cold hardiness and carbohydrate concentration from dormancy induction to bud burst among provenances of three European oak species. Tree Physiol. 2007;27(6):817–25.

    Article  CAS  PubMed  Google Scholar 

  52. Lee JH, Yu DJ, Kim SJ, Choi D, Lee HJ. Intraspecies differences in cold hardiness, carbohydrate content and -amylase gene expression of Vaccinium corymbosum during cold acclimation and deacclimation. Tree Physiol. 2012;32(12):1533–40.

    Article  CAS  PubMed  Google Scholar 

  53. Benina M, Obata T, Mehterov N, Ivanov I, Petrov V, Toneva V, et al. Comparative metabolic profiling of Haberlea rhodopensis, Thellungiella halophyla, and Arabidopsis thaliana exposed to low temperature. Front Plant Sci. 2013;4.

  54. Geigenberger P, Stitt M. Sucrose synthase catalyses a readily reversible reaction in vivo in developing potato tubers and other plant tissues. PLANTA. 1993;189(3):329–39.

    Article  CAS  PubMed  Google Scholar 

  55. Ruan Y. Sucrose metabolism: gateway to diverse carbon use and sugar signaling. Annu Rev Plant Biol. 2014;65(1):33–67.

    Article  CAS  PubMed  Google Scholar 

  56. Baud S, Vaultier MN, Rochat C. Structure and expression profile of the sucrose synthase multigene family in Arabidopsis. J Exp Bot. 2004;55(396):397–409.

    Article  CAS  PubMed  Google Scholar 

  57. Sasaki H, Ichimura K, Imada S, Yamaki S. Sucrose synthase and sucrose phosphate synthase, but not acid invertase, are regulated by cold acclimation and deacclimation in cabbage seedlings. J Plant Physiol. 2001;158(7):847–52.

    Article  CAS  Google Scholar 

  58. Crespi MD, Zabaleta EJ, Pontis HG, Salerno GL. Sucrose synthase expression during cold acclimation in wheat. Plant Physiol. 1991;96(3):887–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Jang I, Oh S, Seo J, Choi W, Song SI, Kim CH, et al. Expression of a bifunctional fusion of the Escherichia coli genes for Trehalose-6-phosphate synthase and Trehalose-6-phosphate phosphatase in transgenic Rice plants increases Trehalose accumulation and abiotic stress tolerance without stunting growth. Plant Physiol. 2003;131(2):516–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Zuther E, Büchel K, Hundertmark M, Stitt M, Hincha DK, Heyer AG. The role of raffinose in the cold acclimation response of Arabidopsis thaliana. FEBS Lett. 2004;576(1–2):169–73.

    Article  CAS  PubMed  Google Scholar 

  61. Chen L, Hou B, Lalonde S, Takanaga H, Hartung ML, Qu X, et al. Sugar transporters for intercellular exchange and nutrition of pathogens. NATURE. 2010;468(7323):527–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Xu Y, Tao Y, Cheung LS, Fan C, Chen L, Xu S, et al. Structures of bacterial homologues of SWEET transporters in two distinct conformations. NATURE. 2014;515(7527):448–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Chen L, Cheung LS, Feng L, Tanner W, Frommer WB. Transport of sugars. Annu Rev Biochem. 2015;84(1):865–94.

    Article  CAS  PubMed  Google Scholar 

  64. Kasuga M, Liu Q, Miura S, Yamaguchi-Shinozaki K, Shinozaki K. Improving plant drought, salt, and freezing tolerance by gene transfer of a single stress-inducible transcription factor. Nat Biotechnol. 1999;17(3):287–91.

    Article  CAS  PubMed  Google Scholar 

  65. Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K. AP2/ERF family transcription factors in plant abiotic stress responses. Biochim Biophys Acta. 1819;2011:86–96.

    Google Scholar 

  66. Liu Q, Kasuga M, Sakuma Y, Abe H, Miura S, Yamaguchi-Shinozaki K, et al. Two transcription factors, DREB1 and DREB2, with an EREBP/AP2 DNA binding domain separate two cellular signal transduction pathways in drought- and low-temperature-responsive gene expression, respectively, in Arabidopsis. Plant Cell. 1998;10(8):1391–406.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Stockinger EJ, Gilmour SJ, Thomashow MF. Arabidopsis thaliana CBF1 encodes an AP2 domain-containing transcriptional activator that binds to the C-repeat/DRE, a cis-acting DNA regulatory element that stimulates transcription in response to low temperature and water deficit. Proc Natl Acad Sci. 1997;94(3):1035–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K. AP2/ERF family transcription factors in plant abiotic stress responses. Biochimica et Biophysica Acta (BBA) - Gene Regulatory Mechanisms. 2012;1819(2):86–96.

    Article  CAS  Google Scholar 

  69. Barros P, Gonçalves N, Saibo N, Oliveira M. Cold acclimation and floral development in almond bud break: insights into the regulatory pathways. J Exp Bot. 2012;63(12):4585–96.

    Article  CAS  PubMed  Google Scholar 

  70. Jaglo-Ottosen KR, Gilmour SJ, Zarka DG, Schabenberger O, Thomashow MF. Arabidopsis CBF1 overexpression induces genes and enhances freezing tolerance. SCIENCE. 1998;280(5360):104–6.

    Article  CAS  PubMed  Google Scholar 

  71. Carretero-Paulet L, Galstyan A, Roig-Villanova I, Martínez-García JF, Bilbao-Castro JR, Robertson DL. Genome-wide classification and evolutionary analysis of the bHLH family of transcription factors in Arabidopsis, poplar, Rice, Moss, and algae. Plant Physiol. 2010;153(3):1398–412.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Xu W, Zhang N, Jiao Y, Li R, Xiao D, Wang Z. The grapevine basic helix-loop-helix (bHLH) transcription factor positively modulates CBF-pathway and confers tolerance to cold-stress in Arabidopsis. Mol Biol Rep. 2014;41(8):5329–42.

    Article  CAS  PubMed  Google Scholar 

  73. Kim YS, Lee M, Lee J, Lee H, Park C. The unified ICE–CBF pathway provides a transcriptional feedback control of freezing tolerance during cold acclimation in Arabidopsis. Plant Mol Biol. 2015;89(1–2):187–201.

    Article  CAS  PubMed  Google Scholar 

  74. Chinnusamy V. ICE1: a regulator of cold-induced transcriptome and freezing tolerance in Arabidopsis. GENE DEV. 2003;17(8):1043–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Zarka D, Vogel J, Cook D, Thomashow M. Cold induction of Arabidopsis CBF genes involves multiple ICE (inducer of CBF expression) promoter elements and a cold-regulatory circuit that is desensitized by low temperature. Plant Physiol. 2003;133(2):910–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Das K, Roychoudhury A. Reactive oxygen species (ROS) and response of antioxidants as ROS-scavengers during environmental stress in plants. Frontiers in Environmental Science. 2014;2.

  77. Jalmi SK, Sinha AK. ROS mediated MAPK signaling in abiotic and biotic stress- striking similarities and differences. Front Plant Sci. 2015;6:769.

    Article  PubMed  PubMed Central  Google Scholar 

  78. Hernandez M, Fernandez-Garcia N, Diaz-Vivancos P, Olmos E. A different role for hydrogen peroxide and the antioxidative system under short and long salt stress in Brassica oleracea roots. J Exp Bot. 2010;61(2):521–35.

    Article  CAS  PubMed  Google Scholar 

  79. Duan X, Cai, Yang J, Chen G, Sang M. Fall Ethephon Application Enhances the Freezing Tolerance of Magnolia wufengensis During Overwintering. FORESTS. 2019;10:868.

    Article  Google Scholar 

  80. Bates LS, Waldren RP, Teare ID. Rapid determination of free proline for water-stress studies. Plant Soil. 1973;39(1):205–7.

    Article  CAS  Google Scholar 

  81. Aihua X, Faju C, Zhongkui J, Ziyang S, Zhonglong Z, Lvyi M. Determination of 4 plant hormones in Magnolia wufengensis by gradient elution high performance liquid chromatography (in Chinese). Chinese Journal of Analysis Laboratory. 2019.

  82. Yang Y, Zhang R, Duan X, Hu Z, Shen M, Leng P. Natural cold acclimation of Ligustrum lucidum in response to exogenous application of paclobutrazol in Beijing. Acta Physiol Plant. 2019;41(1):15.

    Article  CAS  Google Scholar 

  83. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. BIOINFORMATICS. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.

    Article  CAS  PubMed  Google Scholar 

  86. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36(Database issue):D480–4.

    Article  CAS  PubMed  Google Scholar 

  87. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  89. Roberts A, Pachter L. Streaming fragment assignment for real-time analysis of sequencing experiments. Nat Methods. 2013;10(1):71–3.

    Article  CAS  PubMed  Google Scholar 

  90. Anders S, Huber W. Differential expression of RNA-Seq data at the gene level -- the DESeq package. EMBL. 2013.

  91. Chang P, Shen Y, Dong B, Bian S, Wang N, Xuan L, et al. Selection of reference gene in Magnolia denudata under salt stress using qRT-PCR (in Chinese). Journal of Agricultural Biotechnology. 2018;26(9):1611–20.

    Google Scholar 

Download references


We are grateful to Forestry Bureau of Wufeng County for providing materials. We also thank OEbiotech Co., Ltd. for their support in RNA sequencing.


Not applicable.

Author information

Authors and Affiliations



KJW analyzed the data and drafted the manuscript. XJD designed and carried out the experiment. KJW and XJD did the bioinformatics analysis. ZLZ and ZYS provided the experimental materials and participated in data analysis. JD helped to draft the manuscript. ZKJ and LYM initiated and supervised the study. Both KJW and XJD contributed to the writing of this article. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Zhongkui Jia or Luyi Ma.

Ethics declarations

Ethics approval and consent to participate

M. denudata seeds were obtained from a single tree from Henan Lvtu Landscape Ecological Construction Co. Ltd. (North China), and plant materials were formally identified by Prof. Luyi Ma and deposited in the laboratory of the School of Landscape Architecture at Beijing University of Agriculture.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1

. Gene Ontology (GO) analysis of the differentially expressed genes (DEGs) for M. denudata buds in response to cold acclimation.

Additional file 2: Table S2

. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed genes (DEGs) for M. denudata buds in response to cold acclimation.

Additional file 3: Table S3

. Differentially expressed transcription factors (TFs) of M. denudata buds in response to cold acclimation.

Additional file 4: Table S4

. Differentially expressed protein kinases (PKs) of M. denudata buds in response to cold acclimation.

Additional file 5: Table S5

. Differentially expressed genes (DEGs) related hormone metabolism and biosynthesis.

Additional file 6: Table S6

. qRT-PCR primer information.

Additional file 7: Fig. S1

. KEGG map about arginine and proline metabolism and arginine and proline metabolism compared with N1.

Additional file 8: Fig. S2

. KEGG map about starch and sucrose metabolism compared with N1.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wu, K., Duan, X., Zhu, Z. et al. Physiological and transcriptome analysis of Magnolia denudata leaf buds during long-term cold acclimation. BMC Plant Biol 21, 460 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: