Eight-year-old M. wufengensis and M. denudata were cultivated in Jiufeng National Forest Park (Beijing, China; 40°3′25″N, 116°6′39″E). The trees were not clipped or chemically treated before sampling. In 2019–2020, one-year shoots with one apical bud were collected from M. wufengensis on 2 November (19N₁), 23 November (19N₂), 6 December (19D₁), 14 December (19D₂), 21 December (19D₃), 30 December (19D₄), 12 January (20J₁), 18 January (20J₂), and 19 February (20F₁). In 2020–2021, M. wufengensis and M. denudata shoots and buds were collected on 5 October (20O₁), 20 October (20O₂), 5 November (20N₁), 20 November (20N₂), 5 December (20D₁), 20 December (20D₂), 5 January (21J₁), and 20 January (21J₂). 12 trees were divided into three replicates, and all samples were collected from the same 4 trees at each stage. All leaf bud samples of M. wufengensis were stored in liquid nitrogen immediately after collection and then at –80°C until RNA extraction.

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).

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.

Chilling units (CUs) of M. wufengensis and M. denudata leaf buds were calculated based on 0–7.2°C (Weinberger, 1950) and Utah models (Richardson, 1974).

The extraction, purification, and determination of endogenous abscisic acid and gibberellin (GA₁, GA₃ and GA₄) 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₂. 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 for RNA-seq was isolated separately from the 19N₁, 19D₁, and 19D₂ samples using Plant RNA Purification Reagent for plant tissue according to the manufacturer’s instructions (Invitrogen, Carlsbad, CA, United States) and genomic DNA was removed using DNase I (Takara Bio, Shiga, Japan). Then, RNA quality was determined using the 2100 Bioanalyzer Instrument (Agilent Technologies, Santa Clara, CA, United States) and quantified using NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, United States). Only high-quality RNA samples (optical density (OD)260/280 = ∼1.8–2.2, OD260/230 ≥ 2.0, RNA integrity number ≥ 6.5, 28S:18S ≥ 1.0, quantity > 1 μg) were used to construct the sequencing library.

RNA purification, reverse transcription, library construction, and sequencing were performed at Shanghai Majorbio Bio-pharm 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.

The raw paired-end reads were trimmed and quality controlled using SeqPrep¹ and Sickle² with default parameters. Subsequently, clean data from M. wufengensis were used to perform de novo assembly with Trinity³ (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⁴ (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⁵ (Ogata et al., 1999).

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⁶ (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⁷ and KOBAS⁸ (Xie et al., 2011).

Venn diagrams were drawn and trend analysis was performed using Venny 2.1⁹ and Short Time-series Expression Miner software (STEM) (Ernst and Bar-Joseph, 2006), respectively.

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.

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).

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 double-distilled 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.

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).

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₂; 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₂ and ED release between 20N₂ and 20D₁, 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).

GA₃ content was low before 21J₁ and peaked at 21J₂, and a small peak appeared at 20N₂ (before ED release) (Figure 2A) and a similar result was observed regarding the content of GA₁ (Figure 2B). In addition, the content of GA₄ decreased during ED release phase and a peak at 20D₂ was observed (Figure 2B). Moreover, before ED release, the concentration of ABA kept increasing until at 20O₂ and then decreased rapidly to the lowest level at 20N₂ with dormancy release, with a slight increase at 20D₁ and considerable reduction immediately at 20D₂. After ED release, ABA levels sharply increased after 20D₂ and increased further at 21J₂ (Figure 2C). In addition, the ratio of content of ABA/GA₃ increased before ED release and experienced a sharp decrease after ED release (Figure 2D).

Three libraries 19N₁ (ED, control, BBP = 0%), 19D₁ (ED release phase, 0% < BBP < 50%), and 19D₂ (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).

The assembled unigenes were compared against NR, Swiss-Prot, Pfam, COG, GO, and KEGG databases using BLASTX (E-value < 1e^–5). Simultaneously, we found that the number of unigenes successfully annotated to the NR database was the highest (78,142; 41.7%), followed by COG (71,437; 38.1%), KEGG (61,840; 33.0%), Pfam (55,774; 29.8%), Swiss-Prot (55,692; 29.7%), and GO (43,500; 23.2%) (Supplementary Figure 4B). With respect to species, the unigene sequences were most similar to genes from Quercus suber (29,500), Cinnamomum micranthum (11,576), Carpinus fangiana (5,871), Nelumbo nucifera (3,592), and V. vinifera (3,052) using BLASTX matches (Supplementary Figure 4C).

KEGG analysis revealed 61,840 unigenes to be significantly mapped to 139 KEGG pathways and classified into six categories: Metabolism, Genetic Information Processing, Environmental Information Processing, Cellular Processes, Organismal Systems, and Human Diseases. The highest unigenes representation pathways in “Metabolism” were carbohydrate metabolism (4,090) and amino acid metabolism (2,911), those in “Genetic Information Processing” were translation (4,895) and folding, sorting, and degradation (3,342). Signal transduction (823), transport and catabolism (2,138), environmental adaptation (852), and endocrine and metabolic disease (116) were most associated with “Environmental Information Processing,” “Cellular Processes,” “Organismal Systems,” and “Human Diseases,” respectively (Supplementary Figure 4D).

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).

Unigenes with p value < 0.05 or |log2FC| ≥ 1 were defined as DEGs. Among the DEGs, 8,565 and 30,321 were upregulated and 7,675 and 13,672 genes were downregulated at 19D₁ and 19D₂, respectively. Moreover, 28,849 upregulated and 9,715 downregulated unigenes were identified between 19D₂ and 19D₁. The number of DEGs of 19D₂ versus 19N₁ were the highest, followed by 19D₂ versus 19D₁ and 19D₁ versus 19N₁ (Figure 3A). To further explore DEGs related to dormancy release under natural conditions, a Venn diagram was drawn between 19D₁ versus 19N₁, 19D₂ versus 19N₁, and 19D₂ versus 19D₁, 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₁ to 19D₂ were identified (Figure 3C).

Compared to 19N₁, GO analysis of DEGs at 19D₁ 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₂, 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₁.

Compared to 19N₁, KEGG pathway enrichment analysis for DEGs indicated that four pathways – “ribosome (map03010),” “oxidative phosphorylation (map00190),” “plant-pathogen interaction (map04626)” and “plant hormone signal transduction (map04075)” were significantly enriched at 19D₁. At 19D₂ versus 19N₁, “spliceosome (map03040)” followed by “oxidative phosphorylation (map00190),” and “plant hormone signal transduction (map04075)” were significantly enriched. Compared with D₁, “ribosome (map03010),” “oxidative phosphorylation (map00190),” and “spliceosome (map03040),” were significantly enriched at 19D₂ (Supplementary Table 2).

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₁ (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₁, and sharply increased at 19D₁before decreasing to a relatively low level at 19D₂. 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₁, which gradually increased at 19D₁ and part of 19D₂. 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.

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 downregulated), C2C2 with 93 genes (56 upregulated and 37 downregulated), bHLH with 86 genes (45 upregulated and 41 downregulated), bZIP with 80 genes (64 upregulated and 16 downregulated), AP2/ERF with 66 genes (35 upregulated and 31 downregulated), NAC with 55 genes (25 upregulated and 30 downregulated), and WRKY with 48 genes (26 upregulated and 22 downregulated) (Supplementary Table 4).

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).

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₁, 20N₂, 20D₁, 20D₂, 21J₁, and 21J₂ 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₁ and 20N₂ were in ED and ready to release from ED, and 20D₁, 20D₂, 21J₁, and 21J₂ were in ECD. BBP on 20N₁ and 20D₁ 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₂, 21J₁ and 21J₂). 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.

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₂, 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).

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₂ 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₁ and then increased slowly at ECD. The expression level of MwICE-1-2 showed a high expression at 20D₂ but maintained a relatively low expression at other phases (Figure 7).

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₁ 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.

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.

The 0–7.2°C and Utah models are widely used to assess CRs in perennial trees (Yang et al., 2021). As shown in Figure 1B, temperature steadily dropped to 7.2°C relatively late, so no CUs were accumulated before 20N₁ 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₁.

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₃ 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₃ decreases with the length of dormancy (Li et al., 2020). 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₃ was higher than the other two GAs, indicating an important role GA₃ 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₂ and 21J₂, and GA₃ content steadily increased after dormancy release, similar those observed in other species, which suggests a relationship between content of ABA and GA₃ 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 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₃ 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 (Li et al., 2018). 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.

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 (Yang et al., 2021). 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₂, while the temperature was still above 0°C (Supplementary Figure 7). In the coldest two months, which included D₁, D₂, J₁, and J₂, 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.