The 10-year-old tea cultivar “Longjing 43” (Camellia sinensis (L.) O. Kuntze cv. Longjing 43) grown in a field in Hangzhou, Zhejiang Province, China (N30°18′, E120°10′), was used as plant material. The tea plants were provided thorough pest and fertilizer management. We sampled axillary buds from the middle position of the annual branches on different dates between October 2013 and October 2014. Buds collected from adjacent multiple plants were pooled as one sample, and samples collected from three different locations in the field were used as biological replicates. All sampled buds were frozen immediately in liquid nitrogen after being detached from the branches and then stored at -80°C until RNA isolation. The temperature and sunshine duration of each day from September 2013 to April 2014 are shown in Supplementary Material 1-S1.

RNAs were isolated separately from the samples collected on December 1, 2013, February 14, 2014, March 14, 2014, and June 4, 2014, referring to the CTAB method described by Chang et al. (1993). The appearances of shoot and axillary buds on the above sampling dates are recorded in Supplementary Material 1-S2. The extracted RNAs were treated with DNase I (Invitrogen, Carlsbad, CA, USA) to remove contaminating genomic DNA, and then the RNA quality was verified by 1% gel electrophoresis. High-quality RNAs were treated with a Poly(A) Purist^TM-MAG Magnetic mRNA Purification Kit (Ambion, Life Technologies, Lithuania) following the manufacturer’s protocol to enrich mRNA. The enriched mRNAs were used to prepare cDNA libraries using a TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA, USA) following the manufacturer’s instruction. The fragment length distributions of the prepared cDNA libraries were examined with an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). A total of twelve cDNA libraries were sequenced using a whole lane of an Illumina HiSeq 2000 sequencer after qRT-PCR normalization.

After sequencing data were acquired, the adapter sequences were removed, and low-quality reads were removed by the “Sickle” quality-based-trimming program¹ with default settings. The coverage of repeat reads was decreased using the “Trinity normalize by Kmer coverage r2103-08-14” program with parameters of “K-mer size = 25, maximum coverage = 30, maximum pct of mean for stdev of coverage across read = 100”². Finally, the Trinity (v2.2.0) pipeline was run on Kmer-normalized reads. The de novo assembly transcripts were evaluated by “Compute Contig Statistics” and “BUSCO” in the Discovery Environment of CyVerse³. By BlastX analysis (BLAST 2.2.30+) with the NR database⁴ and protein databases of Populus (Ptrichocarpa_210_v3.0.protein.fa and Ptrichocarpa_210_v3.0.annotation_info.txt) and Arabidopsis (Athaliana_167_TAIR10.protein.fa and Athaliana_167_TAIR10. annotation_info.txt) released as parts of Phytozome v11.0.8⁵, the best hits (with a significant E-value of < 1e^-5) were assigned to the assembled transcripts of tea plant.

Using the de novo assembled transcripts as reference sequences, transcript abundance from different sets of RNA-seq reads were quantified using the RSEM v1.2.11 program (Li and Dewey, 2011). The expression patterns and “PPEE” values of each gene and contig were estimated by EBSeq v1.1.5, a useful R package for differential expression analysis of RNA-seq data (Leng et al., 2015). “PPEE” stands for the posterior probability that a gene/transcript is equally expressed. Meanwhile, EBSeq also provides the normalized mean count value for each gene/transcript under each condition, which can be used to calculate fold changes between conditions. Finally, a total of 15 expression patterns were generated to describe the expression changes between four specific conditions in this study (Supplementary Material 2). Each row refers to a different expression pattern, and each column gives the expression status of a different condition. Two conditions are equally expressed if and only if their statuses are the same. We summarized the DEGs among the comparisons of Para vs. Endo, Endo vs. Eco, and Eco vs. Flush (Supplementary Material 2).

A total of 12 samples used in RNA-seq analysis were collected on four different dates. According to the gene expression pattern analysis results, a total of 16,126 DEGs among four different sampling dates were used for clustering analysis. The log2 values of normalized mean counts were used as fold changes of gene/transcript expression level. Using “heatmap.2” in R packages, we made a heat map to show the expression changes of all DEGs between conditions. In this analysis, all DEGs were clustered by the Pearson method, and all samples were clustered by the Spearman method. This information was used to support the dormancy status assessments: paradormant (Para, June 4, 2014), endodormant (Endo, December 1, 2013), ecodormant (Eco, February 14, 2014), and bud flush (Flush, March 14, 2014). To further verify the dormancy states of bud samples on December 1, 2013, and February 14, 2014, we used the following regrowth method. The 9-year-old potted tea cultivar “Longjing 43” grown in the field at the same location where the bud sampling was conducted was used to identify the dormancy status of axillary buds by moving the plants from the field to a greenhouse (14 h of light at 23°C and 10 h of dark at 20°C) after light pruning. Dormancy identification was performed almost every half month from September 2013 to March 2014.

Gene set enrichment analysis was used to identify gene sets that were overrepresented among the DEGs, especially in the comparisons of Para vs. Endo, Endo vs. Eco, and Eco vs. Flush. GSEA is a useful statistical approach for identifying DEG sets (Subramanian et al., 2005). In the above comparisons, we analyzed the up-regulated genes and down-regulated genes separately; in addition, we analyzed the DEGs without considering the direction of change of the individual genes. The GSEA was performed using the Java application GSEA v2.0.13 (Broad Institute, Cambridge, MA, USA); the parameters were set up following the description by Howe et al. (2015). Gene sets were considered statistically significant at an FDR p-value of 0.05. Three Gene Ontology (GO) datasets including biological process, cellular component and molecular function were analyzed. The GO terms of tea genes were assigned using the annotation information of Arabidopsis through available functional annotation. Moreover, signaling pathway enrichment analysis was also included. For the Pathway Studio dataset, six sets of genes, including expression targets, miRNA targets, binding partners, protein modification targets, proteins regulating cell processes and disease/cell process-associated proteins, were used.

We identified the transcription factor genes involved in DEGs using the Arabidopsis thaliana Plant Transcription Factor Database v4.0 (PlantTFDB 4.0⁶; Jin et al., 2017). The chromatin-associated genes and phytohormone-associated genes were also identified based on the previous functional annotation and the available classification data in an extensive transcriptome study of Populus (Howe et al., 2015). The expression changes (log2 value) of all the identified chromatin-associated and phytohormone-associated genes among the four different developmental phases are shown as heat maps produced by the R program. We estimated the standard deviation (SD) of log2 values of mean count values of samples under four different conditions and then ranked the deviation values. The differentially expressed phytohormone-associated genes with SD > 1.5 are shown in the heat map analysis; moreover, the top 60 differentially expressed transcription factor genes with large SDs are shown in the heat map as well.

Recently, extensive transcriptome changes have been investigated during the natural onset and release of vegetative bud dormancy in Populus (Howe et al., 2015). To compare the results of transcriptome analyses in Populus with this study, we summarized the matching DEGs in both studies by sequence homology. All the DEGs among the paradormancy, endodormancy and ecodormancy phases in this study were analyzed together with the data in Supplemental ‘Data File S1’ provided by Howe et al. (2015). In tea plant, the expression patterns of DEGs during two dormancy transitions of Para/Endo and Endo/Eco were identified following the standards described by Howe et al. (2015). The difference in the significance and direction of expression during transitions was identified based on previous pattern analysis in Supplementary Material 2.

A total of 15 DEGs were selected from different expression patterns for qRT-PCR validation. Five micrograms of the total RNA used in the previous RNA-Seq library construction was used for cDNA synthesis. A DNase I kit (Invitrogen, Carlsbad, CA, USA) and a SuperScript III First-Strand Synthesis System (Invitrogen, Carlsbad, CA, USA) were used for cDNA synthesis following the exact instructions of the kits. The products were diluted 20-fold with ultrapure water and used as templates in qRT-PCR analysis on a ViiA^TM 7 Real-Time PCR System (Thermo Fisher Scientific, USA). The primer information for qRT-PCR is listed in Supplementary Material 3-S1. qRT-PCR reactions were set up with 10 μl of 2 X SYBR^® Select Master Mix (4472908, Applied Biosystems^®, USA), 1 μl of forward primer (10 pM/μl), 1 μl of reverse primer (10 pM/μl), 1 μl of cDNA template and 7 μl of water with the following program: 95°C for 2 min and 50 cycles of 95°C for 10 s, annealing at 60°C for 10 s, and elongation at 72°C for 40 s. A heat dissociation protocol from 55 to 95°C was added to detect the specificity of primers. Each reaction was performed in triplicate, and PTB was chosen as a reference gene (Hao et al., 2014). The expression levels of tested genes were calculated by the delta-delta Ct method. For RNA-Seq analysis, we calculated the log2 value using the average of mean count values from three biological replicates and used it as fold change in expression. The expression levels of each detected transcript under different conditions were normalized to the paradormancy phase. The expression levels detected by RNA-Seq and qRT-PCR are shown in the same bar chart.

In this study, the axillary buds collected on December 1, 2013, February 14, 2014, March 14, 2014, and June 4, 2014, were used for RNA-Seq analysis. The appearances of shoots and axillary buds on the above sample collection dates are recorded in Supplementary Material 1. In June, without pruning, tea plants had already produced long shoots after sprouting in the spring. In practice, regrowth always appears after pruning at this time, therefore the axillary buds collected at June were in paradormancy. On March 14, the ambient temperature rapidly increased (Supplementary Material 1-S1), and the dormant buds entered flush after a long, cold winter. On December 1, the axillary buds had undergone months of dormancy induction in early winter after bud set and should have been in endodormancy. On February 14, the dormant axillary buds had already endured a long period of cold, but the outside temperature was still fairly low. The axillary buds should have been in ecodormancy at that time.

To confirm the dormancy status assignment on collected buds, we performed regrowth experiments (Supplementary Material 1-S3). The time to bud flush by the plants collected on February 14 (11 days) was more than that of buds from pruning on the eve of bud flush (March 7th) to “one bud and two leaves” (8 days). Likewise, significantly more days (about 28 days) to bud flush were required by the plants collected on December 1. Therefore, the axillary buds collected on December 1, 2013, February 14, 2014, March 14th, 2014, and June 4th, 2014, were in likely in endodormancy (Endo), ecodormancy (Eco), bud flush (Flush), and paradormancy (Para), respectively. Clustering analysis of DEGs showed that the samples collected at different dates were classified into different group (Figure 1), indicates that the employed axillary buds were in different developmental phases. This further validated our assumption.

The total RNAs extracted from more than twelve axillary bud samples were used for RNA-Seq analysis. In total, twelve RNA-Seq libraries were constructed and sequenced on an Illumina HiSeq 2000 platform. Approximately 70 GB of raw data were harvested (Submission number in NCBI: SRR5040773 to SRR5040784), and more than 60 GB of clean data were obtained after quality-based trimming. In total, 313,388 transcripts were produced after de novo assembly by Trinity (Table 1). The average length of all transcripts was 849 bp, and N50 (length for which half of the total bases are in contigs of this length or longer) was 1538 bp.

All transcripts were annotated by BlastX using the NR database, Arabidopsis protein database and Populus protein database with 50.60, 38.90, and 37.53% of the transcripts annotated, respectively. The four different conditions produced a total of 15 expression patterns (Supplementary Material 2). Only the transcripts assigned as “pattern 1” showed no difference in expression among the four conditions. In this study, we only focused on the DEGs in the comparisons of Para vs. Endo, Endo vs. Eco, and Eco vs. Flush. The expression patterns involved in the above comparisons are listed in Supplementary Material 2, and the DEGs are listed on the following pages. The number of DEGs involved in these comparisons is shown in Table 2. We used hierarchal clustering of expression patterns to group all transcripts (16,125 in total; pattern 2 to pattern 15) by different conditions (Figure 1). The biological replicates collected at each time point had very close expression levels and clustered together. As expected, the samples collected on different dates were clustered into different branches. All the DEGs among the four different conditions that had similar expression patterns were clustered as well. The DEGs clustered into multiple groups that showed clearly different expression patterns under different conditions. These results are consistent with the previous developmental phase assignment.

To validate RNA-seq, we performed cDNA synthesis using the total RNAs for RNA-seq and used the cDNA as a template for qRT-PCR. The expression levels of 15 selected transcripts were detected under four different conditions. Most selected transcripts were found through RNA-Seq analysis to have different expression patterns. The expression results are shown in Supplementary Material 3-S2. The 15 selected transcripts had similar expression patterns between RNA-Seq and qRT-PCR, suggesting reliable expression data by RNA-Seq.

During the DEG analysis, many gene/transcripts were significantly different at the RNA level. It is difficult to identify the major pathways or metabolic processes responsible for axillary bud development regulation by investigating individual genes. In this study, DEG sets were estimated by the GSEA method, which is a powerful tool for analyzing the expression of large numbers of genes (Subramanian et al., 2005; Howe et al., 2015). Using GSEA, GO term enrichment, and Pathway Studio, signaling pathway analyses were carried out according to gene set enrichment identification. Detailed information is provided in Supplementary Materials 4-S1–S3.

We analyzed the enriched gene sets based on GO category including biological process, cellular component and molecular function. The up-regulated, down-regulated, and up-/down-regulated transcripts were considered separately. More GO terms were down-regulated in Para vs. Endo and Endo vs. Eco than in Eco vs. Flush, with a p-value < 0.001 (Supplementary Materials 4-S1–S3). Without considering the direction of change of the individual transcripts, the top 15 enriched GO terms in the biological process, cellular component and molecular function categories among the three comparisons are shown in Figure 2. We also identified the GO term enrichment analysis while limiting the gene sets according to the direction of gene expression, the results of enrichment analysis according to the up- or down-regulated genes provided additional insights into the changes that occurred in axillary buds (Supplementary Materials 4-S1–S3).

In GO term analysis, the terms ‘histone lysine methylation,’ ‘regulation of DNA replication,’ ‘DNA methylation,’ and ‘DNA replication’ were significantly enriched, especially in comparisons of Para vs. Endo and Eco vs. Flush. This result indicates that the chromatin-associated genes may play important roles in regulating axillary bud dormancy formation and release. In total, we identified 127 chromatin-associated genes with differential expression among different developmental phases (Supplementary Material 5). The expression of these genes at each phase is shown in Figure 3. The gene sets related to chromatin in GSEA are listed in Supplementary Materials 4-S5–S7.

The GO term of ‘sequence-specific DNA binding transcription factor activity’ contained the largest number of entities in the molecular function catalog in GSEA. Therefore, we identified all the transcription factors in DEGs among the four conditions. In total, 798 transcription factors were identified from DEGs (Supplementary Material 6), and we chose the top 60 differentially expressed transcription factors and show their expression levels in Figure 4.

To determine the major signaling pathways involved in regulating dormancy status transition, signaling pathway enrichment analysis was performed (Supplementary Material 4-S4). Phytohormone-associated gene sets were mainly enriched in all three comparisons. In Para vs. Endo, ‘cytokinins (CK) signaling,’ and ‘salicylic acid signaling’ were up-regulated; in contrast, ‘gibberellin (GA) signaling’ and ‘plant growth auxin signaling’ were down-regulated. In Endo vs. Eco, ‘cytokinins signaling,’ ‘plant growth auxin signaling,’ and ‘stress ABA signaling’ were up-regulated, while ‘salicylic acid signaling’ and ‘jasmonic acid, ethylene, and salicylic acid crosstalk signaling’ were significantly down-regulated. In Eco vs. Flush, ‘gibberellin signaling’ was down-regulated. These indicated that phytohormones play important roles in dormancy transition regulation. Furthermore, phytohormone-associated genes were identified in all DEGs among the four different conditions. In total, 38 ABA-associated genes, 35 GA-associated genes, 119 auxin-associated genes, 33 CK-associated genes, 31 brassinosteroid (BR)-associated genes, 41 ET-associated genes, and 31 JA-associated genes were identified. The expression patterns of these hormone-associated genes are listed separately in Supplementary Materials 7-S1–S7.

In extensive transcriptome analysis of vegetative buds in Populus, 1,362 DEGs were identified during the dormancy transitions Para/Endo and Endo/Eco. By functional annotation, 417 DEGs in dormancy transitions from Para to Eco in tea plant were matched to 394 out of 1362 DEGs from the poplar study, in which two or more DEGs in tea plant could be matched to same gene in poplar (Supplementary Material 8). The chromatin-associated genes, phytohormone-associated genes, dormancy-related QTL neighbors and transcription factors were examined in the poplar study. Therefore, we statistically analyzed the number of DEGs involved in the above gene groups during endodormancy onset and release (Table 3). In total, 6 out of 13 chromatin-associated genes were enriched in both the poplar and tea plant studies. IAA and JA were the top two hormone-associated signaling pathways enriched by large numbers of DEGs in both studies. In poplar, more DEGs were located in dormancy-related QTLs of LG3, LG5, LG8a and LG8b. Similarly, relatively more DEGs from tea plant were matched to the dormancy-related QTLs. Many DEGs were identified as transcription factors in both the poplar and tea plant studies. In particular, ‘MYB’ was the most enriched transcription factor in both studies. ‘ERF,’ ‘bHLH,’ ‘C2h2,’ ‘WRKY’ and ‘HD-ZIP’ were enriched in sequence. Although the major enriched terms in the poplar study were also matched by DEGs in the tea plant study, only approximately 60% of terms were matched in total. Additionally, the expression patterns of matching DEGs did not show much consistency between poplar and tea plant.

A total of 12 libraries were constructed in RNA-Seq analysis, and 313,388 transcripts were produced after de novo assembly. This transcriptome set was substantially larger than previous dormancy-related RNA-Seq analysis, where only 42,916 assembled sequences were identified from a two leaves and a bud samples (Paul A. et al., 2014). Additionally, our data set appears to be more complete as our average length and N50 of transcripts were up to 849 bp and 1538 bp, respectively, whereas in the previous study, the average length of sequences was approximately 459 bp (Paul A. et al., 2014). Other recent global transcriptome profile analysis on the same tea cultivar reported a total of 127,094 unigenes with a mean length of 355 and an N50 of 506 were obtained using different tea organs as materials, and later a total of 216,831 transcripts, with an average length of 356 bp and an N50 of 529 bp, were harvested from mature leaves (Shi et al., 2011; Wang et al., 2013). In view of the statistical results of our assembly (Table 1), 95.2% of conserved eukaryotic genes that were identified as complete or partial in our assembly, therefore we obtained high-quality sequencing and assembly. The large number of transcripts produced in this study could be attributed to the deep RNA sequencing and large genome size of tea plant (Tanaka et al., 2006). The obtained transcripts were subsequently annotated using three different protein databases.

Although expression profile analyses of buds in paradormancy, endodormancy and ecodormancy have been repeatedly reported in other species, little attention has been paid to buds in the sprouting stage (Ruttink et al., 2007; Horvath et al., 2008; Bai et al., 2013). The introduction of the bud flush phase here helps to draw a relatively complete development cycle of axillary buds. In a previous study, a total of 5,204 DEGs were identified when comparing gene expression in actively growing shoots and dormant shoots of tea plant (Paul A. et al., 2014). Compared with previous studies, more DEGs were identified in this study, and all the DEGs were employed in the subsequent GSEA analysis.

Aside from transcription factor-mediated expression regulation, DNA methylation, histone modifications, and RNA-based mechanisms are the three main epigenetic mechanisms operating in plant development, and they are particularly prominent in plant dormancy studies (Cooke et al., 2012; Ríos et al., 2014). In this study, chromatin-associated gene set analysis showed that DNA methylation- and histone modification-related terms were enriched in the three comparisons. In particular, the terms associated with methylation maintenance and deacetylation, such as ‘neighbors of maintenance of DNA methylation,’ ‘neighbors of MET1’ and ‘binding partners of HD1,’ were down-regulated in Para- vs. Endo-dormancy. In Arabidopsis, MET1 and CMT3 are major maintenance-type DNA methyltransferases, and DRM2 is the major de novo-type DNA methyltransferase; histone H3-K9 methyltransferase KYP and the catalytically non-active protein DRM3 are their functionally related proteins (Ashapkin et al., 2016). Reversible histone acetylation and deacetylation at the N-termini of histone tails constitute an important histone modification mechanism playing a crucial role in the regulation of gene activity (Roth et al., 2001). Histone acetyltransferases (HATs) and histone deacetylases (HDACs) catalyze histone acetylation and deacetylation reactions, respectively, which play essential roles in regulation of gene expression in plant development and plant responses to environmental stresses (Liu et al., 2016). Santamaría et al. (2009) measured DNA methylation and acetylation levels genome-wide during bud set and bud burst and discovered their general association with bud dormancy. Recently, RNA-Seq analysis in the apple tree further reported that high expression levels of DNA methyltransferases and histone methyltransferases were detected during the dormancy and fruit set stages (Kumar et al., 2016). In Picea glauca after short-day induction, HDAC and CMT-type DNA methyltransferase were down-regulated (El Kayal et al., 2011).

The histone methylation and histone H3 acetylation levels in the promoter or intron regions of DAM, a gene involved in dormancy regulation, are associated with dormancy transition (Horvath et al., 2010; Leida et al., 2012; Saito et al., 2015). These findings suggest a significant association between bud development regulation and epigenetic changes, especially in dormancy transitions, that have been observed in other perennials appears to be conserved in the evergreen broad-leafed tea plant as well. In apple buds, HTA8 and HTA12 together with DAM were identified as major dormancy-related genes (Falavigna et al., 2014). Histone modification was investigated in the DAM homolog PpMADS13-1 from Japanese pear; the loss of histone variant H2A.Z coincided with the down-regulation of PpMADS13-1 (Saito et al., 2015). The functions of CHR32 and chromatin-associated genes still require further study, but they might play important roles in regulating DNA unwinding and transcription, depending on functional annotation. The genes FLT1, HMGB13 and HFO8 were the top three DEGs in group C1 and showed significantly elevated expression in the ecodormancy stage.

This study also presents expression profiling of all identified differentially expressed chromatin-associated genes. Most chromatin-associated DEGs were grouped into group B. In this group, DEGs had obvious low expression levels at the endodormancy and ecodormancy stages. Histone-lysine N-methyltransferase, HTA7, CHR32 and chromatin-associated gene were the most enriched four DEGs in group B. In Populus transcriptomic analysis, the GO molecular function gene set of ‘Histone lysine N-methyltransferase’ was also down-regulated from paradormancy to endodormancy (Howe et al., 2015).

Barakat et al. (2012) identified chilling-responsive microRNAs genome-wide in Prunus persica, and histone lysine N-methyltransferase was listed as an important target of four microRNA families (miR5021, miR164, miR396 and miR2919) in the vegetative to reproductive phase transition of the meristem during winter dormancy. Interestingly, in our gene set analysis (Supplementary Material 4) for the comparison of Para vs. Endo, ‘miRNA targets of MIR164A’ was highly enriched.

The functional identification of the CO/FT module in flowering and seasonal growth cessation regulation in trees was a breakthrough in the study of dormancy (Bohlenius et al., 2006). Methylation of H3K27 and H3K4 in FT chromatin regulates its expression level (Jiang et al., 2008; Jeong et al., 2009). Polycomb Repressive Complex 2 (PRC2) and LIKE-HETEROCHROMATIN PROTEIN (LHP1) play roles in the abovementioned methylation regulation mechanism (Turck et al., 2007; Jiang et al., 2008). Histone mark readers MRG1/2 can even modulate FT expression through interaction with CO (Bu et al., 2014). The HMGB family is one group of chromosomal high-mobility-group (HMG) proteins in plants and contains a single HMG-box DNA-binding domain (Grasser, 2003). Members of the HMGB family act as versatile modulators of chromatin function, including transcriptional regulation (Grasser et al., 2006). Histone H4 (HFO) is an important carrier of chromatin modifications on the nucleosome core (Iglesias and Cerdan, 2016).

The chromatin-associated DEGs that showed specific low expression at the ecodormancy phase were mainly grouped into group A. DNG1, ABHF10 and AGO4 were the top genes with expression variation. Similarly, ABHF4 and DNG1 were important chromatin-associated genes identified in poplar transcriptomic analysis; these genes had the same expression patterns during dormancy transition (Howe et al., 2015). The only two Argonaute (AGO) genes with significant expression differences were clustered into Group A in this study. AGO proteins are key players in RNA silencing and can bind small non-coding RNAs, control protein synthesis, and affect messenger RNA stability (Hutvagner and Simard, 2008). The gene sets related to miRNA targets were enriched in the three comparisons, and the AGO proteins play important roles in miRNA function execution. According to the gene set enrichment and chromatin-associated DEG analyses, epigenetic mechanisms may play important roles in tea plant bud development and dormancy regulation.

Many phytohormone-associated transcriptome factors were identified among comparisons, and the expression targets or neighbors of GAI, ABI4, MYC2, BRI1, EIN2/3, ETR1, and PIN3 were majorly enriched in gene set analysis. ABA is a phytohormone regulating many aspects of plant growth and development, including tolerance of a variety of environment stresses (And and Giraudat, 1998). Changes in GA metabolism and signaling occur during early dormancy induction and release (Cooke et al., 2012). ABI4 is an important regulator in seed dormancy by regulating the biogenesis of ABA and GA in Arabidopsis (Shu et al., 2013). MYC2 and other ABA-related transcription factors are up-regulated after 1–2 weeks of short-day induction (Ruttink et al., 2007). Aside from ABA and GA, the roles of BR in promoting germination have attracted more attention in recent seed dormancy studies (Steber and McCourt, 2001). In this study during the bud flush stage, BRI1-related expression targets were mainly up-regulated, which indicated the possible roles of BR in bud sprouting regulation. Implicating BR signaling in tea plant bud dormancy transitions provides further evidence for a potential role for BR first noted in poplar (Howe et al., 2015). Ethylene metabolism- and signaling-related gene sets (expression targets of EIN2/3, ETR1) were enriched in comparison of Endo vs. Eco. Moreover, 12 ERFs were included in the top 60 differentially expressed transcription factors. Based on the expression profiling of the top 60 transcription factors with different expression, three ERFs had high expression levels in the paradormancy and endodormancy stages, and another 9 ERFs showed high expression levels in the endodormancy and ecodormancy stages. Therefore, ERFs might play important roles at the dormancy initiation and deep dormancy stages. APETALA2-ethylene-responsive element binding protein (AP2-EREBP) was identified as a differentially expressed transcription factor among active and dormant tea shoots (Paul A. et al., 2014). In the poplar study, five of six differentially expressed ERFs had high expression levels in endodormancy; in addition, gene sets associated with EIN2 and EIN4 were expressed at high levels during endodormancy (Howe et al., 2015). The abovementioned expression analysis confirmed the hypothesis about the importance of the ethylene signaling pathway in previous dormancy studies (Rohde and Bhalerao, 2007; Ruttink et al., 2007).

Light signal transduction- and circadian clock-associated gene sets and transcription factors were significantly enriched during dormancy transition. The photoreceptor-associated gene sets of ZTL and PHYB were mainly up-regulated in the comparison of Para vs. Endo. ZTL family proteins together with cryptochromes and phototropins are well-known receptors for blue light and phytochromes, particularly PHYA and PHYB, for red and far-red light (Kami et al., 2010). PIFs are bHLH family transcription factors and can directly interact with PHYs (Duek and Fankhauser, 2005; Leivar and Quail, 2011). The transcription of PIF4/PIF5, both negative regulators of PHYB signaling, is regulated by circadian clock-mediated coincidence mechanisms, while PIF3 is not (Huq and Quail, 2002; Nomoto et al., 2012; Soy et al., 2012). The gene sets of PIF4 and PIL6 were significantly down-regulated during the transition from paradormancy to endodormancy in this study. This coincides with the low expression levels of bHLH family genes at the endodormancy and ecodormancy stages in expression profiling analysis (Figure 4). PIF4 also functions in plant responses to environmental cues, increasing auxin levels to the point necessary to promote elongation growth at the proper temperature (Nomoto et al., 2012). PIF4 is even involved downstream of GA signaling and directly interacts with BZR1, a brassinosteroid-activated transcription factor (de Lucas et al., 2008; Oh et al., 2012).

The gene sets related to GI, CCA1, LHY, ELF3, and ZTL were enriched among the dormancy transitions in this study. These gene sets are major components of the circadian clock, a network of transcription-translation negative feedback loops (Horvath, 2009; Cooke et al., 2012). GI-, ZTL- and LHY-related gene sets were up-regulated when buds went into para- and endodormancy. In leafy spurge, several circadian regulators, such as ELF4, GI, FLAVIN-BINDING, KELCHREPEAT, F BOX 1 (FKF1) and CAA1, were significantly up-regulated during endodormancy and ecodormancy. In poplar, repression of LHY delayed bud set was observed; however, the expression of the clock genes LHY1/2, TOC1, and PRR5 varied from bud dormancy formation to release (Ibáñez et al., 2010). The function of the CO/FT module in flowering regulation has been well described (Pin and Nilsson, 2012). Regulation of CO expression and protein stability by light signaling and circadian clock is a critical mechanism in Arabidopsis flowering regulation (Golembeski and Imaizumi, 2015; Hajdu et al., 2015; Lazaro et al., 2015). However, the recent functional discovery of FT in seasonal growth regulation indicated the multiple roles of the CO/FT module in regulating both flowering and dormancy in perennial plants (Bohlenius et al., 2006; Horvath, 2009; Hsu et al., 2011). CO can interact directly with DELLA protein and histone mark readers MRG1/2 during FT regulation (Bu et al., 2014; Xu et al., 2016). When the ecodormant tea buds sprouted, the gene set ‘binding partners of CO’ was significantly up-regulated; this indicated that the CO/FT module may play important roles in seasonal growth regulation in tea plant as well.

MADS-box genes are possible downstream components of CO/FT signaling that are closely related to bud dormancy regulation (Cooke et al., 2012). Dormancy-associated MADS-box (DAM) genes, which are MIKC^C-type MADS-box genes and homologs to SVP and AGL24 in Arabidopsis, were identified in multiple species (Bielenberg et al., 2008; Horvath et al., 2010; Ubi et al., 2010). In leafy spurge, the expression levels of DAM genes were up-regulated during dormancy formation, especially at the endodormancy stage (Horvath et al., 2010). In peach and Japanese apricot trees, the multiple DAM genes showed different seasonal and photoperiodic expression patterns (Sasaki et al., 2011; Yamane et al., 2011). In an extensive expression profiling analysis of poplar, several identified DAM-like genes surprisingly showed lower expression levels during endodormancy (Howe et al., 2015). In our study, several MADS-box genes showed significantly differential expression under the measured conditions (Figure 4). Except for one MADS-box gene (comp27798_c0_seq2) that showed very high expression at ecodormancy, all other MADS-box genes only had high expression in paradormancy, suggesting possible regulatory roles more in line with similar genes from poplar. Further studies should be performed to determine if these MADS-box genes had similar functions in bud dormancy to known DAMs in other species.

Seven MYB genes were involved in the top 60 differentially expressed transcription factors in this study, and three of them (comp100540_c0_seq1, comp76266_c0_seq2, comp62057_c0_seq1) showed high expression levels in endodormancy and ecodormancy. Similarly, MYB62 and MYB4 showed differential expression during short-day induction or dormancy transition in poplar studies (Ruttink et al., 2007; Howe et al., 2015). Functional identification of MYB genes in Arabidopsis indicated that MYB genes control many aspects of plant secondary metabolism, as well as the identity and fate of plant cells (Stracke et al., 2001). However, functional study of the roles of MYB genes in bud dormancy is required.

In plants, FT, auxin and sugars are typical moveable molecules associated with activity-dormancy transitions (Cooke et al., 2012). Plasmodesmata are important channels for cell-to-cell signaling within the apex (van der Schoot and Rinne, 2011). The dynamic accumulation and hydrolyzation of callose on the plasmodesmata are similar to a switch for cellular communication that is closely associated with the states of bud dormancy (Rinne et al., 2001). 1,3-Beta-glucan synthases and 1,3-beta-glucanases are the key enzymes involved in the abovementioned processes (Rinne et al., 2011; van der Schoot and Rinne, 2011). In this study, the GO term ‘plasmodesma’ was up-regulated in paradormant and endodormant buds and down-regulated during the bud flush period. Moreover, the gene set ‘neighbors of CAL’ was up-regulated in the paradormancy and bud flush periods and down-regulated in the endodormancy and ecodormancy periods. More attention has been paid to the β-1,3-glucanase gene family in recent studies. In Arabidopsis, at least fifty β-1,3-glucanase genes mainly grouped into α, β and γ clades were identified, all of which contain an N-terminal signal peptide and a glycosyl hydrolase family 17 domain (Doxey et al., 2007). Approximately one hundred β-1,3-glucanase genes were identified genome-wide in poplar (Rinne et al., 2016). However, only the members of the α-clade in Arabidopsis and the α-clade and γ-clade members in hybrid aspen can localize to the plasmodesma (Knox and Benitez-Alfonso, 2014; Paul L. K. et al., 2014). In addition, their expression patterns vary in response to chilling and short days (Rinne et al., 2011). Therefore, it is crucial to analyze both the expression data and phylogenetic information when we consider the function of a β-1,3-glucanase gene. GA is the most important signaling hormone in β-1,3-glucanase induction during dormancy release; it acts by promoting the hydrolyzation of callose at sieve plates and plasmodesmata (Rinne et al., 2011, 2016). In tea shoots, the level of active GA is extremely low in the early and deep dormancy periods and dramatically increases prior to dormancy release (Nagar and Kumar, 2000). Our current observations also show that the overwintering buds of tea plant have fairly low substance exchange with adjacent organs in deep dormancy periods (the result is in press). The abovementioned evidence indicates that plasmodesma-mediated cellular communication is an important regulation mechanism in the dormancy transition of tea plant. In addition, the regulation of cell wall modification should be considered. The GO terms and gene sets related to the cellulose synthesis pathway were differentially expressed during dormancy transition. Cellulase 2 and polygalacturonase were also identified as DEGs involved in cell wall modification in dormant tea shoots in a previous study. Although the two glucan polymers cellulose and callose are synthesized at the plasma membrane by cellulose or callose synthase complexes, respectively, much cross-talk occurs between the two processes (Schneider et al., 2016). Additionally, cell wall modification can affect the structure of plasmodesmata, causing rapid change in intercellular communication in response to environmental signals (Knox and Benitez-Alfonso, 2014).

Phytohormones are a group of naturally occurring, organic substances that influence multiple physiological processes, mainly growth, differentiation and development, at low concentrations. The important roles of phytohormones in bud dormancy regulation have been frequently reviewed (Horvath et al., 2003; Cooke et al., 2012; Singh et al., 2017). GA and ABA attract substantial attention in early dormancy studies, and their functions in the dormancy mechanism are better understood. An extensive analysis of metabolite and gene expression dynamics was conducted following the timetable of dormancy initiation in poplar (Ruttink et al., 2007). GA signaling repression-related modulators were dramatically up-regulated immediately after the onset of dormancy induction. GA signaling, light signal transduction, and chromatin remodeling were recognized as major early actions stimulated by dormancy induction. The decrease in active GA levels in the early stage of dormancy formation is closely related to growth cessation (Olsen, 2010; Hoffman, 2011). As discussed above, GA is also an important regulator during dormancy release by directly inducing β-1,3-glucanases, playing key roles in the hydrolyzation of callose around plasmodesmata (Rinne et al., 2011). Flowering uses a similar strategy; the induction of GA by decapitation is also important in regrowth regulation of paradormant axillary buds in hybrid aspen (Rinne et al., 2016) and leafy spurge (Chao et al., 2006). ABA is another important hormone induced by short days in poplar following GA signaling (Ruttink et al., 2007). More than 60 genes involved in ABA signal transduction have been studied, and approximately 146 genes have been classified as ectopic ABI3 targets in poplar. In this study, signaling pathway, gene set and DEG analyses generally showed that the GA signal was repressed in the endodormancy and ecodormancy stages while the ABA signal was stimulated. Interestingly, unlike ABI3 in poplar, the gene set ‘expression targets of ABI4’ was significantly enriched in tea plant. Recently, the roles of GA and ABA signaling in poplar endodormancy maintenance were highlighted by Howe et al. (2015). This evidence indicates that GA and ABA are important hormones involved in bud dormancy initiation and maintenance.

Ethylene biosynthesis and signal transduction are triggered between GA and ABA signal peaks during dormancy induction in poplar (Ruttink et al., 2007). Furthermore, ethylene-associated gene set is up-regulated during endodormancy in poplar (Howe et al., 2015). However, the ethylene-associated gene sets were down-regulated during the transitions to endodormancy and ecodormancy in our study. A previous RNA-seq study also showed that transcriptome factor ERF2 was significantly down-regulated in dormant shoots of tea plant compared to active shoots and identified ERF2 as an important DEG related to leaf senescence (Paul A. et al., 2014). The up-regulation of ERFs was closely related to the phenotype of precocious leaf senescence in Arabidopsis (Koyama et al., 2013). Unlike poplar, the mature leaves of tea plant do not fall in autumn and winter. The differences in expression patterns of ethylene-related genes between tea plant and poplar may be caused by the unique features of deciduous plants.

Brassinosteroids (BRs) and jasmonates (JAs) are also universal hormones in plants. BRs play important roles in endogenous regulation of growth and development and JAs in plant defense as well as growth, seed germination, senescence and abscission (Davies, 2010). In previous bud dormancy studies, little attention has been paid to BRs and JAs (Rohde and Bhalerao, 2007; Cooke et al., 2012; Singh et al., 2017). Recently, BR- and JA-associated DEGs were described in detail in a poplar transcriptome analysis (Howe et al., 2015). Both signaling pathways showed differential modulation during transition to endodormancy. BRs together with auxin metabolism can be affected by temperature (Olsen, 2010). JAs are inhibitory signals for cyclin-dependent kinase activity in both the G1–S-phase and the G2–M-phase transitions (Horvath et al., 2003). Based on GSEA results in this study, the BR signaling pathway was mainly down-regulated at the endodormancy and ecodormancy stages, and the JA pathway was down-regulated at paradormancy and endodormancy but stimulated at the ecodormancy and bud flush stages. This clear differential modulation related to dormancy status transition indicates the roles of BRs and JAs in bud dormancy regulation.

Auxin was the first identified plant hormone and is now well known, playing roles in multiple biological processes in plants (Mockaitis and Estelle, 2008; Vanneste and Friml, 2009). The positive roles of auxin in removing dormancy callose from phloem were validated early on by exogenous application of synthetic auxin in Magnolia kobus and Vitis vinifera (Aloni et al., 1991; Aloni and Peterson, 1997). Subsequently, putative polar auxin transport genes, including PttLAX1-PttLAX3 and PttPIN1-PttPIN3, were identified from hybrid aspen, and their expression reduction related to dormancy formation was attributed to induction by environmental cues at the end of the growing season (Schrader et al., 2003). Horvath et al. (2003) drew a comprehensive pathway model for bud dormancy regulation and highlighted the interactions of auxin with other hormone signals. In this study, the enriched gene set of PIN3, CDK and CDC2 indicated the functional roles of auxin during dormancy transition. In total, 119 auxin-associated DEGs were identified during dormancy transitions; they composed the largest group of hormone-related genes (Supplementary Material 7). Many of the auxin associated DEGs with large expression differences were classified into Group A, which clearly showed lower expression levels in endodormancy and ecodormancy buds compared with paradormancy and bud-flush buds (Figure 5). Generally, the GO term ‘plant growth auxin signaling’ was down-regulated in Para vs. Endo. Similarly, the formation of endodormancy in poplar buds was consistent with the low expression of auxin-associated genes in multiple reports (Horvath et al., 2008; Baba et al., 2011; Howe et al., 2015). Low auxin levels during the deep dormancy stage were observed in tea and other species (Li et al., 2003; Nagar and Sood, 2006). These findings suggest that changes in auxin content and auxin-associated genes are necessary for dormancy transitions in tea plant as has been observed in other systems. In particular, a large proportion of auxin-associated genes with large expression differences were auxin early response factors, such as SAUR and its response proteins. The auxin polar transport-related gene sets were mainly down-regulated in Endo vs. Eco. These findings indicate the greater importance of auxin early response factors and auxin polar transport-related genes in dormancy regulation.

In tea plant, from paradormancy to bud flush, the size of overwintering buds gradually increased (Wang et al., 2014). It is unclear from our results if bud expansion is due to cell division, however, many cell cycle regulators are up-regulated as the buds transitioned from endodormancy to paradormancy, suggesting that there may be cell division occurring in the winter once paradormancy is released (Supplementary Material 2). CK signaling was up-regulated in both endodormant and ecodormant buds, while a couple of particular CK-related gene sets were down-regulated in Para vs. Endo and up-regulated in Endo vs. Eco. Interestingly, compared with other hormones, the CK-associated genes with large expression differences did not show clear expression tendencies. CK likely plays important roles in cell division in dormant buds. However, weak linkage has been constructed between CK and dormancy regulation (Horvath et al., 2003; Cooke et al., 2012). Obviously, dormancy relies on a complex regulation network between hormones in plants, and most physiological processes involve multiple hormones (Vanstraelen and Benková, 2012). For example, cytokinin signaling directly interacts with auxin by regulating PIN-FORMED auxin transporters, and cytokinin-auxin cross-talk may function in outgrowth of paradormant buds (Muller and Leyser, 2011; Simaskova et al., 2015). D-type cyclins and cyclin-dependent kinase involved in cell cycle regulation are regulated collectively by various hormones, including auxin, CK, BR and GA (Horvath et al., 2003). Therefore, the functions of hormones in dormancy regulation should be investigated from a broad viewpoint.

Tea plant is a thermophilic evergreen woody species that only goes into dormancy when exposed to a cold winter (beyond approximately 16° north or south) (Barua, 1969). Poplar is a deciduous woody species that is mainly distributed in temperate and cold temperate zones. As in peach, plum, and pear trees, winter cold/dormancy is necessary for the normal life cycle of poplar, but it is not necessary for tea plant. Its winter dormancy feature may have developed over the course of its slow, progressive northward migration because tea plant originates in the warm southern region of China (Yu, 1986). Moreover, the dormancy of poplar can be induced by short days, and tea plant dormancy is induced by cold and short days (Barua, 1969; Resman et al., 2010; Paul and Kumar, 2011). The mechanism of bud dormancy has been well studied in poplar, and many genes involved in dormancy regulation have been identified (Ruttink et al., 2007; Anderson et al., 2010; Singh et al., 2017). Recently, an extensive transcriptome study discovered the most important dormancy-related genes in poplar following a long time scale (Howe et al., 2015). This gives us a chance to compare the major dormancy-associated genes in tea plant and poplar on a global transcriptomic level. Based on the sequence homology analysis, we annotated the transcripts of tea plant using a poplar protein database and then summarized the DEGs, focusing particularly on chromatin, hormone, QTL and transcription factor classifications in both species (Table last one). Surprisingly, the enriched genes or gene families with a large number of members were almost the same in both species. Epigenetic mechanisms are important regulators in bud dormancy regulation (Cooke et al., 2012; Singh et al., 2017). Approximately half of the chromatin-associated genes identified in poplar also show significantly different expression in tea plant. One example is the AGO4 gene, which shows clear down-regulation during endodormancy in poplar and leafy spurge and can form RNA-induced silencing complexes with short interfering RNAs to mediate DNA methylation and transcriptional gene silencing (Horvath et al., 2008; Matzke and Mosher, 2014; Howe et al., 2015). AGO4 (comp37668_c0_seq1, comp42629_c0_seq1) also showed significant down-regulation during endodormancy in this study. Regarding hormone signaling pathways, auxin-associated genes had the most members in both poplar and tea plant. The auxin-associated DEGs were comprehensively summarized, and their potential roles in dormancy were highlighted in our study and previous studies (Wang et al., 2014; Howe et al., 2015). Moreover, the DEGs located near dormancy-related QTLs had a high match rate with DEGs in tea plant. Although there is no validated dormancy-related QTL in tea plant yet, this result provides us important inspiration for tea plant dormancy research. bHLH, C2H2, ERF, HD-ZIP, MYB, and WRKY are transcription factor families with large numbers of DEGs in both poplar and tea plant. Numerous reports about these transcription factors in current and previous dormancy studies indicate their general and important roles in dormancy regulation (Ruttink et al., 2007; Paul A. et al., 2014; Wang et al., 2014; Howe et al., 2015). However, the differentially expressed MIKC^C-type MADS-box transcription factors, including DAM-like genes in poplar, did not have matching targets in tea plant. This suggests the functional variation of DAM-related genes in tea plant. It is speculated that the molecular pathways involved in tea plant dormancy regulation are consistent with those of poplar to a certain extent. However, the specific functions of individual genes in dormancy regulation of tea plant require further identification. Moreover, the lack of intact transcript information from tea plant limits homology analysis with poplar; therefore, more reliable comparison is still meaningful to improve our understanding of tea plant dormancy when tea plant genome information is available.