Transcriptomic Analysis Reveals the Positive Role of Abscisic Acid in Endodormancy Maintenance of Leaf Buds of Magnolia wufengensis

Magnolia wufengensis (Magnoliaceae) is a deciduous landscape species, known for its ornamental value with uniquely shaped and coloured tepals. The species has been introduced to many cities in south China, but low temperatures limit the expansion of this species in cold regions. Bud dormancy is critical for plants to survive in cold environments during the winter. In this study, we performed transcriptomic analysis of leaf buds using RNA sequencing and compared their gene expression during endodormancy, endodormancy release, and ecodormancy. A total of 187,406 unigenes were generated with an average length of 621.82 bp (N50 = 895 bp). In the transcriptomic analysis, differentially expressed genes involved in metabolism and signal transduction of hormones especially abscisic acid (ABA) were substantially annotated during dormancy transition. Our results showed that ABA at a concentration of 100 μM promoted dormancy maintenance in buds of M. wufengensis. Furthermore, the expression of genes related to ABA biosynthesis, catabolism, and signalling pathway was analysed by qPCR. We found that the expression of MwCYP707A-1-2 was consistent with ABA content and the dormancy transition phase, indicating that MwCYP707A-1-2 played a role in endodormancy release. In addition, the upregulation of MwCBF1 during dormancy release highlighted the enhancement of cold resistance. This study provides new insights into the cold tolerance of M. wufengensis in the winter from bud dormancy based on RNA-sequencing and offers fundamental data for further research on breeding improvement of M. wufengensis.


INTRODUCTION
Owing to the instability in global climate, many perennial plants have suffered from abnormal weather conditions, including extreme temperatures in winter. Bud dormancy is the temporary suspension of visible growth in plant buds and represents a protective strategy for perennial plants to survive unfavourable climatic changes during winter (Rohde and Bhalerao, 2007). Bud dormancy is traditionally categorised into three phases: paradormancy (PD), inhabited by substances generated from another part of the plant; endodormancy (ED), controlled by internal factors; and ecodormancy (ECD), regulated by the external environment (Lang et al., 1987;Considine and Considine, 2016). Plants cannot resume growth in a favourable environment until ED release (Rohde and Bhalerao, 2007). To break ED, plants need to fulfil chilling requirements (CRs) after accumulating sufficient chilling hours (Arora et al., 2003), as insufficient cold accumulation may delay dormancy release, influence flower morphology, and even impair growth and production (Atkinson et al., 2013). Therefore, it is necessary to evaluate bud dormancy status and assess CRs in perennial trees. Three models are mainly used to calculate CRs in woody perennials: 0-7.2 • C model (Weinberger, 1950), Utah model (Richardson, 1974), and dynamic model (Fishman et al., 1987a,b).
Temperature and photoperiod are important environmental signals controlling the seasonal dormancy cycle in perennials (Anderson et al., 2010;Maurya and Bhalerao, 2017). Short photoperiods induce bud formation, bud dormancy induction and apical meristem cessation of shoots (Weiser, 1970;Singh et al., 2017). Moreover, dormancy release requires sufficient chilling accumulation in winter as low temperature mostly regulate dormancy release and bud break (Heide, 2008;Yamane et al., 2011). According to the two different dormancy-related environmental factors, plants can be classified into three types: temperature-sensitive, photoperiod-sensitive, and temperatureand photoperiod-sensitive (Bai et al., 2016).
Phytohormones are a crucial factor influencing bud dormancy in perennials, and endogenous hormones and their balance regulate the induction of and release from dormancy (Sonnewald and Sonnewald, 2014;Liu and Sherif, 2019). Some conventional hormones such as gibberellin (GA), abscisic acid (ABA), and auxin (IAA) participate in the dormancy cycle (Horvath et al., 2003;Vimont et al., 2021). In general, an increase in the ABA content accompanied by a decrease in the GA 3 and IAA content is observed during the dormant induction phase, whereas the opposite trend is observed during dormancy release in plants (Liu and Sherif, 2019). High levels of IAA and GA 3 accelerate dormancy release (Rinne et al., 2011;Zhuang et al., 2013) whereas ABA maintains dormancy Tylewicz et al., 2018). In addition, exogenous ABA application results in a delay in bud break in, for example, Pyrus pyrifolia (pear) , Vitis vinifera (grape) (Zheng et al., 2015), and Betula pendula (birch) (Rinne et al., 1994).
The role of ABA in dormancy has been widely studied at physiological and molecular levels and evidence has indicated that ABA biosynthesis, catabolism, and signalling pathway are involved in the regulation of bud dormancy .
A rate-limiting enzyme involved in ABA biosynthesis, 9-cisepoxycarotenoid dioxygenase (NCED), has been indicated to control dormancy at the transcriptional level (Zheng et al., 2015;Li et al., 2018). During catabolism, ABA is degraded by ABA 8'-hydroxylase, which is encoded by cytochrome P450 CYP707A, and the relationship between CYP707A and ABA content has been widely investigated (Cutler and Krochko, 1999;Saito et al., 2004). The ABA signalling pathway consists of two groups of ABA regulators: Protein Phosphatase 2c (PP2Cs) and SNF1-Related Protein Kinase 2 (SnRK2s). Besides, the ABA receptors were identified as Pyrabactin Resistance (PYRs), Pyrabactin Resistance-Like (PYLs), and Regulatory Component of ABA Receptor (RCARs) (Hubbard et al., 2010). ABA binds to PYR/PYL/RCARs and forms PP2C complexes, which inhibit the activity of PP2Cs. PP2Cs can suppress SNF1-related protein kinase 2 (SnRK2s) function via dephosphorylation, which negatively affects ABA signalling, allowing SnRK2s to activate the downstream ABRE-binding factor (AREB/ABF) transcription factors (TFs) (Umezawa et al., 2009;Soon et al., 2011). Several studies have shown that genes related to ABA signalling are involved in dormancy regulation . In Hybrid Aspen, short days induce high levels of ABA which suppresses PICKLE (PKL) to induce the expression of SVP-like (SVL), which is an orthologue of short vegetative phase (SVP) and then SVL induces callose synthase 1 (CALS1) expression to promote the establishment of dormancy Singh et al., 2019).
Bud dormancy is an important overwintering process, and many studies have shown that bud dormancy is associated with winter cold resistance at the molecular level. C-repeat binding factor (CBF) belongs to the APETALA2/-ETHYLENE RESPONSE FACTOR (AP2/ERF) gene family, regulates many genes related to cold response and tolerance and can be induced by inducer of CBF expression (ICE) (Chinnusamy et al., 2007;Guo et al., 2018). Dormancy-associated MADS-box (DAM)/SVP/SVL genes are known to control bud dormancy in many species Yang et al., 2018;Gao et al., 2021). Therefore, the relationship between DAM and CBF links between bud dormancy and cold resistance. PmCBFs are known to bind to the promoter of PmDAM6 and activate the expression of PmCBFs in P. mume (Zhao et al., 2018a,b). In P. pyrifolia, the expression of PpDAMs is directly induced by the accumulation of CBF by binding to CRT/DRE motifs (Niu et al., 2015;Saito et al., 2015). Li et al. (2019) reported that low temperature induces PpCBF1-PpDAM2 regulon to function during ED transition . Thus, bud dormancy may be associated with cold tolerance during the winter.
Magnoliaceae plants have high ornamental value and are widely cultivated globally. Magnolia denudata, as a common species in north China, has been widely cultivated for its prominent cold tolerance (Yang et al., 2015). Magnolia wufengensis (Supplementary Figure 1), a new species of Magnoliaceae, was discovered growing in Wufeng County, Hubei Province, People's Republic of China (Ma et al., 2006). As a deciduous landscape species with uniquely shaped and colours of the tepals, M. wufengensis has been introduced to many cities in south China for its rich biological characteristics and will have a place in global horticultural plants (Shi et al., 2021). However, this is difficult in north China where temperatures can be extremely low, because M. wufengensis is more sensitive to the cold and with a deeper dormancy level than other Magnoliaceae species such as M. denudata (Yang et al., 2015;Deng et al., 2019;Duan et al., 2019). Bud dormancy is an important biological process that helps plants survive cold temperature in winter. RNA sequencing (RNA-seq) has been recently used to study bud dormancy in many species such as pear (P. pyrifolia) (Bai et al., 2013), tea (Camellia sinensis) (Hao et al., 2017), sweet cherry (Prunus avium L.) (Vimont et al., 2019), and wintersweet (Chimonanthus praecox) . In this study, using RNA-seq, we aimed to explore: (i) the cycle period between ED and ECD and the effects of different meteorological factors on dormancy release of M. wufengensis, (ii) which key genes and pathways were involved in regulation of different dormancy phases, (iii) the role that hormones, especially ABA, play in endodormancy maintenance, and (iv) the relationship between cold tolerance and different phases of dormancy. This study will provide a foundation for improving cold resistance and thus allowing normal growth in winter and expanding the northern boundary of M. wufengensis cultivation.

Evaluation of Bud Dormancy Status
The bud dormancy status of leaf buds was evaluated as previously described (Liu et al., 2012) with some modifications. In 2019-2020 and 2020-2021, 1-year-old shoots with one apical bud about 10 cm long (n = 7) were sampled and inserted in wet flower mud in a box full of water and allowed to grow in the climate chamber at 25 ± 1.0 • C during the day and 22 ± 1.0 • C during the night, with a photoperiod of 14 h light/10 h dark and 60% relative humidity. Twenty-one leaf buds grown in three boxes of flower mud were divided into three biological replicates. The water was changed and the base of the shoots was cut every 3-4 days. Dormancy status was determined by the bud break percentage (BBP) after 32 days. We defined the beginning of bud break when the as green leaf tips were enclosing visible leaves. If the buds were at more than 50% bud break after 32 days, then the buds were considered to be released from ED (Yooyongwech et al., 2009).

Acquisition of Meteorological Data and Chilling Units in Magnolia wufengensis and Magnolia denudata During Dormancy
Maximum (T max ), minimum (T min ), and average temperatures (T avg ) were recorded every 15 min using the weather station (WeatherHawk, Campbell Scientific, UT, United States) located in Jiufeng National Forest Park.

Measurements of Phytohormones Contents
The extraction, purification, and determination of endogenous abscisic acid and gibberellin (GA 1 , GA 3 and GA 4 ) were performed using an enzyme-linked immunosorbent assay (ELISA) according to manufacturer's instructions. The fresh samples (1 g bud) were homogenised in liquid nitrogen and extracted in pre-cold 80% (v/v) methanol with butylated hydroxytoluene (BHT) (1 mmol/L) and kept at 4 • C overnight. The samples were centrifuged for 15 min at 5,000 rpm (4 • C). Afterward, the extracts were passed through a C18 Sep-Pak Cartridge (Waters, Milford, MA, United States) and dried with N 2 . Then the residues were dissolved in 0.01 mol L −1 PBS (pH 7.4) to determine the levels of ABA and GAs content. Calculations of plant hormones by ELISA followed the protocol described in Zhao et al. (2006). The ELISA kits used for the assay were purchased from Saipei Biotechnology Co., Ltd. (Wuhan, China). Each experiment contained three biological and technical replicates.
RNA purification, reverse transcription, library construction, and sequencing were performed at Shanghai Majorbio Biopharm Biotechnology Co., Ltd. (Shanghai, China) according to the manufacturer's instructions (Illumina, San Diego, CA, United States). The RNA-seq transcriptome libraries of M. wufengensis were prepared using TruSeq RNA Sample Prep Kit (Illumina). Poly(A) mRNA was purified from total RNA using oligo-dT-attached magnetic beads (Invitrogen) and then fragmented using the fragmentation buffer. Using these short fragments as templates, double-stranded complementary DNA (cDNA) was synthesised using SuperScript Double-Stranded cDNA Synthesis Kit (Invitrogen) with random hexamer primers (Illumina). Subsequently, the synthesised cDNA was subjected to end-repair, phosphorylation, and "A" base addition according to Illumina's library construction protocol. Libraries were selected for size using cDNA target fragments of 200-300 bp on 2% Low Range Ultra Agarose (Bio-Rad) followed by PCR amplification using Phusion DNA polymerase (New England Biolabs, Boston, MA, United States) for 15 PCR cycles. After quantification using TBS380, two RNA-seq libraries were sequenced in a single lane on NovaSeq 6000 Sequencing System (Illumina) for 2 × 150 bp paired-end reads. Each experiment included three biological replicates.

De novo Assembly and Sequence Annotation
The raw paired-end reads were trimmed and quality controlled using SeqPrep 1 and Sickle 2 with default parameters. Subsequently, clean data from M. wufengensis were used to perform de novo assembly with Trinity 3 (Grabherr et al., 2011). All the assembled transcripts were searched against the National Center for Biotechnology Information protein non-redundant (NR), Clusters of Orthologous Genes (COG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases using BLASTX to identify the proteins that had the highest sequence similarity with the given transcripts to retrieve their function annotations and typical cut-off E-values were set as less than 1.0 × 10 −5 . Blast2GO 4 (Conesa et al., 2005) programme was used to obtain gene ontology (GO) annotations of uniquely assembled transcripts for describing their biological processes, molecular functions, and cellular components. Metabolic pathway analysis was performed using KEGG 5 (Ogata et al., 1999).

Differential Expression Analysis and Functional Enrichment
To identify differentially expressed genes (DEGs) between two different samples, the expression level of each transcript was calculated according to the transcripts per million reads method. RSEM 6 (Li and Dewey, 2011) was used to quantify gene abundance. DEG analysis was performed using DESeq2 (Love et al., 2014) and EdgeR (Robinson et al., 2009) with DEGs | log2FC| > 1 and Q value ≤ 0.05 (DESeq2 or EdgeR) considered to be significant. In addition, functional enrichment analysis using GO and KEGG were performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways with a Bonferroni-corrected p value ≤ 0.05 compared to the whole-transcriptome background. GO functional enrichment and KEGG pathway analyses were carried out using Goatools 7 and KOBAS 8 (Xie et al., 2011).
Venn diagrams were drawn and trend analysis was performed using Venny 2.1 9 and Short Time-series Expression Miner software (STEM) (Ernst and Bar-Joseph, 2006), respectively.

Exogenous Abscisic Acid Treatment
For ABA treatment, nine shoots of M. wufengensis were collected from the 12 trees from November 2020 to January 2021 and sprayed with 100, 200, or 300 µM ABA (Aidlab, Beijing, China) and 0.2% ethanol (mock treatment) at approximately 13:00 for three consecutive days. Each treatment was executed with three biological and technical replicates. Buds were collected at 7, 14, 21, and 28 days after ABA treatment and stored immediately in liquid nitrogen and then at -80 • C. Other shoots that grew in the chamber environment mentioned above for 32 days were used to measure BBP.

Identification and Validation of Cold-Related Genes by Cold Acclimation
To ensure DEGs identified by RNA-seq involve in cold tolerance of M. wufengensis, a cold acclimation experiment was conducted to valid their functions. Three apical buds were selected randomly and used for experiments during the whole cold acclimation treatment in September 2017 (autumn). The experiment included three temperatures for analysis. A room temperature of 22 • C in a low-temperature incubator (3M, United States) served as the control. The samples in other groups of M. wufengensis buds were treated for seven days sequentially in low-temperatures incubators at the following two different experimental temperatures: low temperature of 12 and 4 • C. Two days were left for cooling slowly in temperatures incubator between two groups (Supplementary Figure 2).

Quantitative PCR Analysis of Gene Expression
RNA of M. wufengensis buds for quantitative PCR (qPCR) was extracted using the HiPure HP Plant RNA Mini Kit (Magen, Shanghai, China) according to the manufacturer's instructions and genomic DNA was removed using DNase I. cDNA used for qPCR was reverse transcribed from 2 µg of purified RNA in a 20 µL reaction volume based on the manufacturer's instructions (G592, Applied Biological Materials, Richmond, BC, Canada). The qPCR primers were designed using Beacon Designer 7 (PREMIER Biosoft International, Palo Alto, CA, United States) and passed the specificity test. qPCR was carried out on a StepOnePlus Real-Time PCR System (Applied Biosystems, United States) with a final reaction volume of 10 µL containing 5 µL TB Green Premix Ex Taq (Tli RNaseH Plus; Takara Bio) (2X), 0.2 µL each of ROX Reference Dye (50X) (Takara Bio), upstream primer, and downstream primer, 1.0 µL cDNA, and 3.4 µL doubledistilled water. MwACTIN was used for a reference gene for analysis. The primer sequences used in qPCR are listed in Supplementary Table 1. Each sample included three biological and technical replicates.

Statistical Analysis
The study was conducted with a completely randomised design. The data were analysed using one-way analysis of variance followed by least significant difference test and p value <0.05 was considered significant. Graphs were constructed using SigmaPlot version 10 (Systat Software, San Jose, CA, United States) and R Project (R Foundation for Statistical Computing, Vienna, Austria). All data were analysed using SPSS Statistics version 20 (IBM, Armonk, NY, United States).

Dormancy Status and Chilling Requirement of Buds in Magnolia wufengensis and Magnolia denudata During Natural Overwintering
To study the relationships between bud dormancy and cold tolerance, it is imperative to define the status of bud dormancy. As is shown in Figure 1A, no bud breaking was observed in M. wufengensis and M. denudata on 20O 2 ; however, BBP increased with progress in chilling accumulation mainly in November and December. In M. wufengensis, the apical leaf buds were determined in the ED phase before 20N 2 and ED release between 20N 2 and 20D 1 , when the number of CUs reached 62-214 CUs and 480-548 CUs based on the 0-7.2 • C and Utah models, respectively. In addition, the ED release occurred between 5 December and 20 December in M. denudata, which was later than that in M. wufengensis with CUs reaching 214-294 CUs and 548-589.5 CUs based on the 0-7.2 • C and Utah models, respectively ( Figure 1B).

Phytohormone Concentration During
Bud Dormancy of Magnolia wufengensis GA 3 content was low before 21J 1 and peaked at 21J 2 , and a small peak appeared at 20N 2 (before ED release) (Figure 2A) and a similar result was observed regarding the content of GA 1 (Figure 2B). In addition, the content of GA 4 decreased during ED release phase and a peak at 20D 2 was observed ( Figure 2B). Moreover, before ED release, the concentration of ABA kept increasing until at 20O 2 and then decreased rapidly to the lowest level at 20N 2 with dormancy release, with a slight increase at 20D 1 and considerable reduction immediately at 20D 2 . After ED release, ABA levels sharply increased after 20D 2 and increased further at 21J 2 ( Figure 2C). In addition, the ratio of content of ABA/GA 3 increased before ED release and experienced a sharp decrease after ED release ( Figure 2D). Transcriptome Sequencing, de novo Assembly, and Annotation of Magnolia wufengensis Unigenes During Bud Dormancy Three libraries 19N 1 (ED, control, BBP = 0%), 19D 1 (ED release phase, 0% < BBP < 50%), and 19D 2 (ECD, BBP > 50%) (Supplementary Figure 3) were constructed from cDNA obtained from more than three apical buds and sequenced on the NovaSeq 6000 platform. Approximately 59.75 GB of clean reads was obtained after quality control, and Q30 percentage and guanine and cytosine content (GC) percentage were more than 92.14 and 46.92% in the nine samples, respectively ( Table 1). The de novo assembly using Trinity yielded 187,406 unigenes ranging from 201 bp to 14,669 bp with an average length of 621.82 bp, and N50 of 895 bp ( Table 2). In general, the number of unigenes decreased with the increase in gene length, and the largest proportion of unigenes was between 200 bp and 500 bp (123,562, 66%), followed by 501 bp to 1,000 bp (35,770, 19%), and 1,001 bp to 1,500 bp (11,895, 6%) (Supplementary Figure 4A).
Based on GO analysis, 43,500 unigenes were classified into three main categories: biological process, cellular component, and molecular function. Biological process was mainly comprised of proteins involved in cellular process (29,572), metabolic process (27,167), and biological regulation (8,340). The most represented cellular components were cell part (29,063), membrane part (22,285), and organelle (16,681). In addition, we found a high number of unigenes involved in binding (36,757), catalytic activity (35,199), and transporter activity (5,244) in molecular function (Supplementary Figure 4E).     (Figure 3A). To further explore DEGs related to dormancy release under natural conditions, a Venn diagram was drawn between 19D 1 versus 19N 1 , 19D 2 versus 19N 1 , and 19D 2 versus 19D 1 , and 4,286 DEGs were found to intersect all three groups ( Figure 3B). To distinguish the changing patterns in gene expression, gene expression profile clustering was performed. From this, 4,286 genes were assigned to 16 different profiles by STEM and six profiles that were significantly enriched from 19N 1 to 19D 2 were identified ( Figure 3C). Compared to 19N 1 , GO analysis of DEGs at 19D 1 demonstrated that genes related to membrane structure and transcription were overexpressed. Biological process, cellular component, and molecular function, "single-organism transport, " "cellular component, " and "oxidoreductase activity" were the most enriched GO categories. At 19D 2 , in biological process, the major subcategories were "metabolic process" and "single-organism process." In cellular component, "cellular component, " "cell part, " and "intracellular part" were the most representative subcategories, and "oxidoreductase activity, " "transporter activity, " and "RNA binding" were the top three subcategories compared to 19N 1 .

Hormone Signal Transduction Related Genes Were Expressed During Bud Dormancy Transition
Based on KEGG annotation, DEGs related to phytohormones play an important role in dormancy transition. A total of 51 DEGs related to plant hormone signal transduction (Supplementary Table 3) divided into three main gene clusters as shown in Figure 4. Cluster A (13 genes) was highly expressed at the 19N 1 (ED) stage, and showed low expression levels during dormancy release, which indicated that the genes in this cluster may be involved in breaking ED. Among the DEGs in cluster A, two genes involved in ABA signalling, PP2C and ABF, and five auxin-related genes, including three IAA, and one each of GH3 and AUX, were found, which indicates that ABA and auxin signalling were activated during dormancy release. Cluster B contained 14 genes that showed low expression level at N 1 , and sharply increased at 19D 1 before decreasing to a relatively low level at 19D 2 . Among them, DEGs associated with auxin (IAA, SAUR, ARF, and T1R1), ethylene (EBF and ETR), and brassinosteroid (TCH4) regulation were most abundant in this expression profile cluster. Cluster C, with the largest number of DEGs (24 genes), exhibited a low expression level at 19N 1 , which gradually increased at 19D 1 and part of 19D 2 . Among them, DEGs that responded to jasmonic acid (COI1, JAZ, and JAR1), cytokinin (ARR-A), auxin (SAUR, ARF, IAA and AUX1), and brassinosteroid (TCH4 and BZR1) exhibited a similar expression pattern to that in cluster C.

Transcription Factors Were Active During Endodormancy Release
A total of 1,001 TF genes (552 upregulated and 449 downregulated) that were active during dormancy transition were identified. These TFs were mainly concentrated in MYB_superfamily, C2H2, C2C2, bHLH, bZIP, AP2/ERF, NAC, and WRKY. Among all the evaluated genes, the number of genes of the MYB_superfamily was the highest with 131 genes (59 upregulated and 72 downregulated), followed by C2H2 with 103 genes (90 upregulated and 13  Validation of RNA-Seq Results Using Quantitative PCR Ten DEGs were randomly selected to demonstrate the reliability of RNA-seq using qPCR. The trends of genes during different dormancy phases using qPCR were consistent with the RNA-seq results, indicating favourable reliability of RNA-seq (Supplementary Figure 5).

Abscisic Acid at a Concentration of 100 µM Promoted Endodormancy Maintenance in Magnolia wufengensis
Based on the transcriptional analysis, we analysed and inferred that hormone metabolism, signal transduction, and especially ABA may play important roles in dormancy transition. Therefore, to further figure out the function of ABA with respect to ED release, three different exogenous concentrations of ABA (100, 200, and 300 µM) were applied to apical buds on 20N 1 , 20N 2 , 20D 1 , 20D 2 , 21J 1 , and 21J 2 in M. wufengensis, and their BBP with exogenous ABA and mock treatments (0.2% ethanol) was compared after 32 days. Based on the evaluated dormancy status, 20N 1 and 20N 2 were in ED and ready to release from ED, and 20D 1 , 20D 2 , 21J 1 , and 21J 2 were in ECD. BBP on 20N 1 and 20D 1 significantly decreased after treatment with 100 µM ABA (Figure 5), whereas ABA at concentration had almost no effect on the germination rate in ECD (20D 2 , 21J 1 and 21J 2 ). This implies that during ED or the ED release phase, ABA at a concentration of 100 µM played a positive role in ED maintenance and was therefore selected for further investigation.

Identification and Expression of Abscisic Acid-Related Genes During Dormancy Transition Under Natural Environment
To further study the molecular mechanism of action of ABA on ED release, one NCED, three CYP707A, two PYL, five PP2C, one SNRK2, and three ABI DEGs related to ABA synthesis, metabolism, and signalling were identified using RNA-seq. Among these genes, the expression of MwNCED-3 was downregulated and showed a low expression before dormancy release, and then increased during dormancy release ( Figure 6A). The expression of MwCYP707A-1-2 declined toward ED release and steadily increased after ED release, which is consistent with the content of ABA, whereas almost no expression of MwCYP707A-1-1 and MwCYP707A-2 was observed during the entire dormancy release phase ( Figure 6B). Moreover, the expression patterns of genes related to ABA signalling were also determined. The expression levels of MwPYL-1/3 genes decreased before dormancy release and increased rapidly during dormancy release, and then decreased steadily at ECD. The same expression pattern was observed in MwPP2C-6. On the contrary, MwSNRK2-10 slowly increased before dormancy release, peaked on 20N 2 , decreased rapidly during dormancy release, and increased thereafter at ECD. In addition, MwPP2C-24 and MwABI-5 increased before dormancy release and dropped with ECD development (Figure 6C).

Expression of C-Repeat Binding Factor and Inducer of C-Repeat Binding Factor Expression Genes During Dormancy Transition Under Natural Conditions
To understand the relationship between dormancy transition and cold tolerance in winter, we identified one CBF1 and two ICE1 genes in our transcriptome data and measured their expression patterns during the natural ED process in M. wufengensis. To confirm their functions on cold resistance, we conducted an experiment under cold acclimation and found that MwCBF-1 and MwICE-1-1 were induced under cold stress, indicating that the two genes were associated with cold tolerance (Supplementary  Figure 6). Among these genes, the expression of MwCBF-1 was upregulated to 20N 2 and showed a high expression before dormancy release, and then decreased but still kept a high expression after ED release. The expression of MwICE-1-1 decreased to N 1 and then increased slowly at ECD. The expression level of MwICE-1-2 showed a high expression at 20D 2 but maintained a relatively low expression at other phases (Figure 7).

Expression Analysis of Genes Related to Abscisic Acid and Cold Tolerance in Response to Exogenous Abscisic Acid Treatment
The significantly decreased BBP indicated that ABA promoted ED maintenance. To further determine the function of ABA in the maintenance of ED, the responses of buds collected at 20D 1 were compared between mock and ABA treatment group. Genes related to dormancy transition under natural conditions were focussed on (Figure 8). The expression of MwCBF-1 induced by ABA was considerably upregulated compared to that in the mock treatment and a similar increase was observed in MwPYL-1 and MwABI-5. Fluctuations were observed in the expression of MwNCED-3 and MwPYL-3 during the whole treatment time. In addition, a slight decreasing trend can be seen in the expression of the two PP2C genes in the previous 21 days, while there was a significant increase in MwPP2C-6 on the 28th day after treatment which was the opposite to the trend observed for MwSnRK-2-10. The expression of MwCYP707A-1-2 was similar to that of MwPP2C-6, showing a high expression at the 7th and 28th days and maintaining a relatively low expression at the 14th and 21st days after ABA treatment.

Expression Analysis of Genes Related to Bud Dormancy in Response to Exogenous Abscisic Acid Treatment
In addition to analysing the expression pattern of ABA and cold-related genes, several genes such as D-type cyclin (CYCD), PKL and CALS1 involving in bud dormancy have also been identified. Among the genes, one CYCD3, one PKL and two CALS1 genes were differently expressed during dormancy transition. To further study whether a similar model induced by ABA exists in M. wufengensis, we measured the expression of the dormancy-related genes under ABA treatment  ( Figure 9). The expression of MwCYCD-3 increased in the previous 21 days and was considerably promoted by ABA and a similar trend exists in MwCALS-1-1 in the first 21 days. In addition, the expression of MwPKL was depressed by ABA throughout the treatment time except the 21st day. During dormancy release, the expression MwCYCD3 and MwCALS-1-1/2 were promoted by ABA, but MwPKL expression was depressed by ABA.

Chilling Requirement and Environmental Factors Affected Bud Dormancy
The 0-7.2 • C and Utah models are widely used to assess CRs in perennial trees . As shown in Figure 1B, temperature steadily dropped to 7.2 • C relatively late, so no CUs were accumulated before 20N 1 under the 0-7.2 • C model. The Utah model is often used for cold regions such as north China, as it takes into consideration the accumulation effect of different temperatures (Erez et al., 1990). Our results were thus consistent with those of previous studies. Based on the Utah model, M. wufengensis underwent lesser CR to break dormancy than M. denudata. In addition, we found that the lower the temperature, the higher the BBP. Based on our results, the buds accumulated enough CUs to break ED before 20D 1 .

Transcriptome Data Revealed Phytohormones Involved in Bud Dormancy Release
Bud dormancy is associated with phytohormones in many species (Cooke et al., 2012). GA and ABA regulate dormancy induction and release (Yue et al., 2018;Yang et al., 2019). In general, ABA levels increase with the establishment of dormancy (PD to ED) and decrease during dormancy release (ED to ECD) while GA 3 exhibits opposite trends during dormancy transition. For example, in wintersweet, the content of ABA increases with the chilling accumulation from PD to ED and then decreases after dormancy release, and that of GA 3 decreases with the length of dormancy . Similar results have been observed in many other perennial species such as pear (P. pyrifolia) (Tuan et al., 2017;Ito et al., 2021), peach (Prunus persica) (Wang et al., 2016), grape (V. vinifera) (Zheng et al., 2015), and leafy spurge (Euphorbia esula) (Chao et al., 2017). In addition, besides changes of GAs content, we also found that the content of GA 3 was higher than the other two GAs, indicating an important role GA 3 plays in dormancy transition of M. wufengensis.
Many DEGs were significantly annotated to hormones during dormancy transition based on KEGG annotation. Therefore, the effect of hormones on dormancy was focussed on. In M. wufengensis, the ABA content decreased during dormancy release but increased rapidly from 20D 2 and 21J 2 , and GA 3 content steadily increased after dormancy release, similar those observed in other species, which suggests a relationship between content of ABA and GA 3 and the depth of ED. Based on our records, the average temperature during this stage was -4.4 • C. Therefore, we hypothesised that a protective strategy to cope FIGURE 9 | Relative gene expression levels of specific differentially expressed genes (DEGs) involved in bud dormancy under ABA treatment during 2020-2021 in Magnolia wufengensis. Each experiment was performed with three biological replicates. Each bar represents the mean ± SEM of three biological replicates.
with low temperatures may exist in M. wufengensis. In addition, a homeostatic network of various hormones is thought to be at the centre of dormancy transition (Shu et al., 2013;Zhang et al., 2018). In C. praecox, the ratio of ABA/GA 3 increases with dormancy breaking (induced by chilling) and decreases during dormancy release in P. mume (Wen et al., 2016;Li et al., 2020). Thus, these results are consistent with those of other studies.
Studies of exogenous hormone application have shown that exogenous ABA application can effectively promote dormancy establishment and maintain ED (Zheng et al., 2015;Li and Dami, 2016). Based on previous studies, inhibitory effect of exogenous ABA does not always work all the time and depends on status of buds. In pear (Pyrus fauriei), ABA inhibitory effect on bud break could be affected by chilling accumulation (Tamura et al., 2002) and in grape (V. vinifera), ABA application to buds which are released from dormancy did not reduce BBP (Zheng et al., 2015). Similar results were found in the present study as 100 µM ABA delayed dormancy release before ECD, whereas 200 and 300 µM ABA could not effectively inhibit BBP. On the one hand, inhibition was not dependent on concentration, which is not consistent with the results observed in grape whose BBP is more efficiently inhibited by high concentration of ABA (Zheng et al., 2015). On the other hand, high concentration of ABA could not suppress bud break in M. wufengensis during dormancy release, which is different from results in pear whose BBP with 100, 200, and 300 µM ABA was similar . Compared to 100 µM whose concentration is relatively low, we infer that 200 or 300 µM ABA may damage and then stimulus the defence system of buds in M. wufengensis, so higher concentration of ABA treatment has not depressed bud break percentage. These results suggest that a non-ABA-regulated controlling mechanism of dormancy may exist in M. wufengensis, and the mechanism needs further study. However, as our experiment was conducted in a climate chamber after buds were cut from trees, whether similar results can be observed in field experiments remains to be investigated.
In addition to physiological data about hormones, we further analysed the data at the molecular level. The content of ABA in plants is not dominated by a single factor but by a balance of biosynthesis and metabolism, and function through various signalling pathways. Therefore, we further analysed the changes in ABA content during dormancy by evaluating ABA biosynthesis, metabolism, and signalling. NCED and CYP707A are two of the main genes involved in ABA biosynthesis and metabolism, respectively. Overexpression of NCED promotes seed dormancy and delays germination in tobacco (Nicotiana plumbaginifolia) (Qin and Zeevaart, 2002). Li et al. (2018) found that PpNCED-2 and PpNCED-3 are highly expressed during ED and decrease rapidly during dormancy release. In this transcriptome data, we identified one NCED DEG and named the gene MwNCED-3. The expression pattern was similar to PpNCED-3 and consistent with the decrease in ABA content during ED released.
Furthermore, CYP707A is known to be highly expressed during dormancy release in peach and pear (Wang et al., 2016;Tuan et al., 2017). Recent studies have shown that CYP707A is involved in dormancy release regulation. In potato, ABA content changes are correlated with NCED and CYP707A gene families and associated with dormancy release in tubers (Destefano-Beltrán et al., 2006). In the present study, the expression of MwCYP707A-1-2 was downregulated during dormancy release and immediately upregulated during ECD ( Figure 6B). This result is not consistent with those of the abovementioned research in peach or pear. Based on these results, we infer that MwCYP707A-1-2 may be a key gene involved in ED release, and the role of MwCYP707A-1-2 in bud dormancy of M. wufengensis needs to be further studied.
ABA signal transduction is also involved in dormancy regulation and bud dormancy. In pear, Bai et al. (2013) found that PpPP2Cs were upregulated while PpSNRK2s were downregulated after dormancy release. Similar results were observed by Li et al. (2018) who found the expression of PpPYLs, PpSNRK2s, and PpABIs to be upregulated from PD to ED, but the expression of PpPP2Cs was low during ED and increased with the decrease of ABA content during dormancy release. The expression of VvPP2Cs higher at PD than at ED (Tuan et al., 2017). From our results, genes associated with ABA signalling showed different expression patterns during dormancy transition. MwPP2C-6/8/51, MwPYL-1/3, and MwSNRK2-10 showed a relatively high expression in ED and ED release, while the expression of MwPP2C-16/24 and MwABI-5/2-2 in ECD was higher than that in ED. However, the consistency between the trends of their expression and dormancy transition was not significant and their functions need to be further studied.
In addition, ABA could not only influence the expression of ABA-related genes, but also involve in regulation of many dormancy-related genes. In the present study, these dormancy-related genes such as CYCD, PKL and CALS1 were differentially affected by ABA: MwCYCD-3 and MwCALS-1-1/2 were promoted by ABA and MwPKL was depressed by ABA which indicated that a similar ABA-centred model of H. Aspen exists in dormancy release M. wufengensis. In addition, functions of the dormancy-related genes need to be further study.

Cold Tolerance and Bud Dormancy in Winter
Bud dormancy, an important biological process for plants to survive the winter, is strongly related to cold hardiness enhancement in winter. Li et al. (2003) found that cold tolerance is enhanced before dormancy development in silver birch (B. pendula). In the winter, perennial plants tend to enhance cold tolerance through natural cold acclimation to survive under long-term cold conditions (Uemura et al., 1995;Lee and Thomashow, 2012). The CBF-dependent signalling pathway is an important cold signalling pathway in plants. CBFs/DREBs bind to cis-elements of cold resistance genes and activate their expression, thus improving cold resistance in plants (Park et al., 2015;Ding et al., 2020). ICE1, a positive regulator of cold response, can activate the expression of CBF (Chinnusamy, 2003;Lee et al., 2005). DAM/SVP can regulate the bud dormancy cycle in perennial trees (Falavigna et al., 2019). Recently, the relationship between CBF and DAM co-regulating dormancy has been widely studied in many species . In P. mume, PmDAM6 and PmCBFs mainly respond to chilling temperature (below 20 • C) and freezing cold (0 • C), respectively. Additionally, eight PmCBFs were upregulated under the stimulus of a cold signal, which then induced the expression of all six DAM genes during dormancy development (Zhang et al., 2018;Zhao et al., 2018a). In the Populus hybrid, PtCBF1 and PtDAM1 induction was found to be related to ED development (Boldizsár et al., 2021). In pear, Niu et al. (2015) found that PpCBF can induce the expression of PpDAM and PpDAM and subsequently inhibit PpFT, which then stimulates growth cessation and promote dormancy maintenance. In the present study, we identified one DREB1B/CBF1 gene (TRINITY_DN8378_c0_g1) and two ICE1 genes (TRINITY_DN20323_c1_g1 and TRINITY_DN20323_c0_g2). MwCBF-1 was highly expressed during overwintering, so we suspected that MwCBF-1 was a positive regulator for cold resistance in M. wufengensis. Before dormancy release, MwCBF-1 achieved a peak at 20N 2 , while the temperature was still above 0 • C (Supplementary Figure 7). In the coldest two months, which included D 1 , D 2 , J 1 , and J 2 , MwCBF-1 decreased but maintained a high expression to cope with cold stress when dormancy was released (Figure 8). Furthermore, the expression of MwCBF-1 can be efficiently induced by exogenous ABA application. We suspected that M. wufengensis can efficiently enhance cold resistance during the dormancy release phase to survive the winter, so it is important to maintain ED and extend the time of the dormancy release phase. Above all, a hypothesis of molecular model for ABA and its biosynthesis, metabolism and signalling pathway; cold tolerance and acclimation; and dormancy during overwintering was proposed (Figure 10). In addition, the role that ECD plays in cold tolerance enhancement during overwintering needs to be further study.

CONCLUSION
Overall, this study provides fundamental insights into the bud dormancy cycle and CR in two Magnoliaceae plants and we hypothesised that M. wufengensis and M. denudata are both sensitive to low temperature and short day based on meteorological data. The content of ABA and GA 3 , and the ABA/GA 3 ratio significantly changed during dormancy release and ECD. A M. wufengensis dataset containing 187,406 unigenes was constructed to observe the dynamic changes in gene expression under different dormancy phases using RNA-seq. Comparison with ED led to the identification of 16,240 and 43,993 DEGs during ED release (19D 1 ) and ECD (19D 2 ), respectively. Among the DEGs, many key genes and metabolic pathways, especially those of plant hormones, were identified using KEGG and GO analyses. Based on heatmap analysis of plant hormone transduction, we found that auxin-and ABA-related genes showed high expression in ED. Thus, auxin and ABA may regulate dormancy transition in M. wufengensis. Application of 100 µM of exogenous ABA before dormancy release could effectively maintain dormancy. Seventeen DEGs involved in ABA biosynthesis, metabolism, and signal transduction were identified based on RNA-seq data. We conducted qPCR on the 17 ABA-related DEGs and found that MwCYP707A-1-2 may be involved in dormancy regulation. Besides, MwCBF-1 was highly expressed during dormancy release, suggesting a relationship between cold tolerance and bud dormancy. Thus, our findings shed light on the mechanism underlying dormancy release and further our understanding of overwintering from bud dormancy in M. wufengensis.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: National Center for Biotechnology Information (NCBI) BioProject database under accession number PRJNA695868 (https://www.ncbi.nlm.nih. gov/bioproject/PRJNA695868).

AUTHOR CONTRIBUTIONS
KW performed most of the experiment, analysed the data, and drafted the manuscript. KW and XD designed the experiment, did the bioinformatics analysis, and contributed to the writing of this article. ZZ and ZS provided the experimental materials and participated in data analysis. YZ involved in conducting experiments. HL analysed the data. ZJ and LM initiated and supervised the study. All authors contributed to the article and approved the submitted version.

FUNDING
This research was supported by lateral research of Beijing Forestry University (662004284).