Time-coursed transcriptome analysis identifies key expressional regulation in growth cessation and dormancy induced by short days in Paulownia

Maintaining the viability of the apical shoot is critical for continued vertical growth in plants. Terminal shoot of tree species Paulownia cannot regrow in subsequent years. The short day (SD) treatment leads to apical growth cessation and dormancy. To understand the molecular basis of this, we further conducted global RNA-Seq based transcriptomic analysis in apical shoots to check regulation of gene expression. We obtained ~219 million paired-end 125-bp Illumina reads from five time-courses and de novo assembled them to yield 49,054 unigenes. Compared with the untreated control, we identified 1540 differentially expressed genes (DEGs) which were found to involve in 116 metabolic pathways. Expression of 87% of DEGs exhibited switch-on or switch-off pattern, indicating key roles in growth cessation. Most DEGs were enriched in the biological process of gene ontology categories and at later treatment stages. The pathways of auxin and circadian network were most affected and the expression of associated DEGs was characterised. During SD induction, auxin genes IAA, ARF and SAURs were down-regulated and circadian genes including PIF3 and PRR5 were up-regulated. PEPC in photosynthesis was constitutively upregulated, suggesting a still high CO2 concentrating activity; however, the converting CO2 to G3P in the Calvin cycle is low, supported by reduced expression of GAPDH encoding the catalysing enzyme for this step. This indicates a de-coupling point in the carbon fixation. The results help elucidate the molecular mechanisms for SD inducing dormancy and cessation in apical shoots.


Results
Short day induces growth cessation and dormancy in Paulownia. We treated one-year-old young trees propagated from roots of Paulownia cultivar P. tomentosa with short-day cycles of 8 h light and 16 h darkness in a greenhouse during the fastest growing season in August. 0.5 cm shoot tips from eight trees were sampled on Aug 1 st (S1) as a control before the SD treatment, and similar tissues were collected at five-day intervals on Aug 6 (S2), 11 (S3), 16 (S4) and 21 (S5). Obvious growth cessation and dormancy induction in shoot tips were observed after SD treatment (Fig. 1a-e), whereas the shoot tips under a long-day photoperiod were still in a fast-growing condition, the same as at S1. The apical shoots from stage S4 had similar phenotype to that of dormant shoots at the end of September when the growth was stopped and enlarged vessel cells and increased xylem were formed in the cambia zone of apical shoot 34 . We found that the resulting terminal buds did not regrow in the following year. However, after treating one-year-old trees with constant light, the trees exhibited a very late growth cessation www.nature.com/scientificreports www.nature.com/scientificreports/ followed by dormancy and terminal buds can regrow in the next year. The terminal buds of one-year-old trees did not re-grow in the coming year even after being transferred to a tropical area in Guangzhou, China (GPS location 23.157 N, 113.359 E). Taking these results together, we concluded that SD, instead of LD, induced dormancy and death in the terminal shoots of the Paulownia.
De novo assembly of transcripts in Paulownia. Transcriptomes from each stage of the collected samples were examined with RNA-Seq technology. More than 39.7 million paired-end 125-bp Illumina reads were generated from mRNA extracted from samples at each stage, totalling 219,619,968 reads with a quality score Phred Q30 > 92% from the five stages (Supplementary Table S1). All reads were de novo assembled into transcripts with Trinity software (version 201308) 35 , and 49,054 unigenes were obtained. More than 75% of the paired reads were properly mapped back to the unigene assembly. The RNA-Seq reads have an average coverage depth of 220 X and unigenes have average length of 853 bp. The RNA-Seq data and unigene assemblies are publicly available at NCBI (Supplementary Table S1).
functional annotation of transcripts. The transcripts were annotated by similarity searching against the NR databases of NCBI, EggNOG, KEGG and InterPro (E-value ≤ 10-10 for Pfam, E-value ≤ 10-5 for others). A total of 54,880 (55%) transcripts were annotated with at least one hit from the databases (Supplementary  Table S1). A total of 24,184 (24%) transcripts were functionally annotated in the orthologous database EggNOG. Among the orthologs, the largest portion of the transcripts (14,175; 43%) belonged to the class of cellular process and signalling, dominated by functional groups "posttranslational modification, protein turnover, chaperones," "signal transduction mechanisms" and "intracellular trafficking, secretion, and vesicular transport" (Supplementary Table S2). A total of 9150 (28%) transcripts were in the class of information storage and processing, dominated by "transcription. " Of the transcripts, 9475 (29%) belonged to metabolism, dominated by "carbohydrate transport and metabolism" (Supplementary Table S2). Further enrichment analysis by using KOBAS 36 revealed that transcripts in Paulownia were mapped to 255 KEGG pathways. Overview of differential gene expression induced by short day treatment. We measured the genome-wide transcript abundance from our RNA-Seq data in fragments per kilobase length per million reads (FPKM). The results showed that the average abundance of all transcripts display a decreased trend with increasing days of SD treatment, reaching a lower level after 15 days of treatment at stage S4 (p < 0.05, Fig. 1f). Pearson correlation analysis showed that the expression levels were more similar at two successive stages, indicating a gradual change of expression in apical shoots ( Supplementary Fig. S1).
We then investigated the differential expression (DE) of genes in response to the SD treatment in time courses and then focused on the DE compared with the starting point, as a control under LD (16 h light and 8 h darkness) just before the SD treatment. In total, the expression of 1540, 1540 and 1346 genes were changed at least two-fold (p < 0.0001), four-fold (p < 0.0001) and eight-fold (p < 0.0001), respectively, suggesting an expressional regulation of many genes in response to SD treatment (Supplementary Table S3). Most (87.4%) of the DE genes that had two or four-fold change also had an eight-fold change, suggesting that a short-day photoperiod induced a dramatic change in gene expression (Supplementary Table S3, Fig. 2a). Then, we focused on the 1540 DE genes for subsequent analyses. The expression patterns fell into five major clusters, and 318 genes in cluster 4 were constitutively switched off after SD treatment (Fig. 2a,b). Across all SD stages, 491common DE genes, consisting of 188 up-regulated and 303 down-regulated DE genes, were identified compared with LD stage S1 (Fig. 2c). Comparison revealed that more DE genes were differentially regulated with accumulated SD treatment (Fig. 2d). The information of DE gene IDs and their up-or down-regulation at each SD stage is available at Supplementary  Table S4.

Gene ontology (GO) classification and involved pathways of differentially expressed genes.
The annotation of DE genes was retrieved from unigene annotation for functional characterisation and annotation results against database NR, Interpro, EGGNOG, KEGG were available in Supplementary Tables S5-S8, respectively. In total, 860, 1414, 2983 and 2734 of the DE transcripts from stages S2, S3, S4 and S5, respectively, were assigned GO terms (Supplementary Table S6). Across all stages, the biological process was the most enriched GO category, followed by the molecular function category (Fig. 3a). Metabolic process, single-organism process, Figure 2. Expression patterns of differential genes in time-course series to short-day treatment in Paulownia. Genes expression displayed here was at least 2 folds of change and statistically p < 0.0001. S2, S3, S4, and S5 represent the shoot tip samples after 5, 10, 15, and 20 days of short-day treatment while S1 represents for long day condition before treatment. Heatmap (a) shows the gene expression level in color after log2 transformation of FPKM value. Each row represents expression level of one gene at different samples. The tree in the left shows the clusters of expression pattern. Image (b) shows the clustered patterns of gene expression. Venn diagram (c) shows the common and specific differentially up-regulated and down-regulated genes across samples. Image d shows the number of up-and down-regulated genes relative to S1. www.nature.com/scientificreports www.nature.com/scientificreports/ cellular process, localisation, and biological regulation were the most abundant groups in the biological process category. In the molecular function category, binding and catalytic activity were the most abundant groups. In the cellular component category, cell part, macromolecular complex, membrane, membrane part and organelle were the most abundant groups, though fewer genes were involved than in the other two categories (Fig. 3a).
To understand the regulatory networks of DE genes, we conducted pathway mapping against a pathway database KEGG using the KOBAS tool (version 3.0). The results revealed that DE genes were mapped to 209 pathways. Of those, 116 pathways were common in all stages of SD treatment, indicating continued gene regulation in these pathways. Two, eleven and seventeen additional pathways were specific to stage S3, S4 and S5, respectively (Fig. 3b, Supplementary Table S9). This suggested that longer SD treatment affected more pathways. The most enriched DE pathways included ribosome for translation, plant hormone signal transduction, carbon fixation in photosynthetic organisms for energy metabolism, and the circadian rhythm network. Change in ribosome for translation suggested altered protein synthesis. Here, we focused our subsequent analysis on the latter three pathways. No specific pathway was found in early stage S2 after only five days of SD treatment (Fig. 3b).

Expression patterns of DE genes in hormone signal transduction.
We characterised the DEGs involved in plant hormone signal transduction. Further GO network enrichment analysis using Blast2GO (version 4.0.2) revealed that 12 DEGs were under the same GO accession (0009725), which is described as "in response to hormone" (Supplementary Fig. S2). KEGG pathway mapping with KAAS tool against Arabidopsis and Populus identified that these genes encode key proteins functioning in multiple hormone signalling pathways (Table 1), suggesting regulation of hormones under SD. The auxin pathway had the greatest number of enriched DEGs, and we annotated these gene functions based on protein homolog similarity (Table 1). Key genes encoding AUX/IAA, auxin response factor (ARF), and small auxin up RNA (SAUR) family proteins in the auxin pathway were regulated significantly in both RNA-Seq results and RT-qPCR results with good coefficients of 0.9, 1, 0.9, 0.8, 1, 0.9 for IAA_C42198, ARF_C52200, ARF_C52619, SAUR_C39339, SAUR_C43763 and SAUR_C10372, respectively, suggesting a large impact on the auxin pathway (Fig. 4). The gene IAA was down-regulated through all stages with SD treatment and reached the lowest expression level at S4, after 15 days of treatment. Two copies of ARF were both down-regulated from S2 to S4 and then up-regulated at S5. As in other plants 37 , multiple SAURs were identified in the auxin pathway in Paulownia in this study. Two of three copies of SAUR were up-regulated, www.nature.com/scientificreports www.nature.com/scientificreports/ but at different rates, and the third copy displayed alternative up-and down-regulation (Fig. 4). In each of the other hormone pathways, only one DEG was regulated significantly during SD treatment, e.g. CRE1 in the cytokinin pathway, EBF1/2 in the auxin and ethylene pathway, BZR1/2 in the brassinosteroid pathway, JAR1 in the jasmonic acid pathway, GID1 in the gibberellin pathway, and SnRK2 in the abscisic acid pathway. The detailed gene expression patterns of these genes during each stage under SD were compared and are displayed in Fig. 4.

Gene regulation in the circadian rhythmic network.
Gene expression enrichment showed that many differentially expressed genes were in the circadian rhythm network. We examined the expression of these genes and mapped them to the well-documented circadian rhythm network of Arabidopsis and Populus with KAAS at the KEGG website. We identified and annotated six key genes, such as PIF3, PRR5, PKF1, MYB and SAP1 (Table 2), with altered regulation in Paulownia under SD and mapped them onto the KEGG map MAP04712 (Fig. 5a). We validated the relative expression patterns with RT-qPCR which were consistent with RNA-Seq results with the coefficients of 1, 0.9, 1, 0.9, 0.9 and 0.9 for gene SPA1, PIF3, PKF1, PRR5, MYB and FT, respectively. Most of these genes showed up-regulated expression patterns but at different rates, except for MYB and FT (Fig. 5b). The expression of the key transcription factor PIF3, which functions in passing the SD signal from phytochromes to the internal clock loop, was up-regulated, which was consistent with previous reports in Populus under SD 3 (Fig. 5a). CO/FT is responsible for seasonal dormancy in trees 3 , but we did not find significant changes in expression of CO under SD in one-year-old Paulownia. FT expression was very low under both LD and SD in RNA-Seq results, except that an up-regulated level of FT at a late stage, e.g. S4 (Fig. 5b,c), which was also reported previously in Populus under SD 3 .

Gene regulation in carbon fixation of photosynthesis under SD.
To understand the role of DEGs in the enriched carbon fixation pathway, we further characterised these genes based on Arabidopsis and Populus in KEGG. Ten DEGs were mapped to three processes of the carbon fixation pathway of photosynthesis, including the C4-Dicarboxylic acid cycle, crassulacean acid metabolism (CAM), and the Calvin cycle (Fig. 6a). These DEGs encoded ten key enzymes catalysing the reactions of this pathway (Fig. 7, Supplementary Table S10). Most genes were upregulated, and the remaining genes showed dynamic up or down patterns (Fig. 6b). In the CO 2 concentration process known in mesophyll cells, the gene c53559_g1, encoding PEPC EC:4.1.1.31 in the first step of binding CO 2 to oxaloacetic acid, and the gene c48076_g1 encoding malate dehydrogenase (MDH) EC:1.1.1.37 for the conversion of oxaloacetic acid to malate were up-regulated under SD. The expression patterns of MDH and PEPC in both RNA-Seq and RT-qPCR analysis showed a good correlation with coefficients of 0.5 and 0.5, respectively (Fig. 6b,c, Supplementary Table S10), suggesting a higher CO 2 -concentrating activity. All known genes encoding the enzymes in the known process CAM for CO 2 residue delivery were upregulated. In the Calvin cycle, being the hub of metabolite conversion, the expression of six genes was altered. The gene GAPDH coding key enzyme glyceraldehyde 3-phosphate dehydrogenase (EC:1.2.1.12) was down-regulated, which indicated low glyceraldehyde-3-phosphate (G3P) generation from delivered CO 2 (Fig. 6). Expression of the gene encoding Rubisco EC:4.1.1.39 was not changed.

Discussion
Dormancy is an evolutionary adaptation to environmental factors, such as stresses, in perennial plants. It is known that an SD photoperiod can induce bud dormancy in some perennial plants 1 . Different from existing physiological reports on Paulownia 33 , here, we reported that SD, similar to winter days, induced growth cessation and dormancy, which mimics natural apical shoot tip in winter. Our transcriptome analyses identified DE genes and their regulation in pathways in Paulownia shoot tips under SD-induced growth cessation and dormancy induction in a time course during the fast-growing season. To our knowledge, this is a novel characterisation of the dynamics of global transcriptomes in response to SD treatment. www.nature.com/scientificreports www.nature.com/scientificreports/ It had been thought that growth cessation, dormancy induction and bud formation are independent of each other 38,39 . Fifteen days of SD is sufficient to induce growth cessation and dormancy in the apical buds of Paulownia, and the buds could enter certain stages of dormancy. Previous report already shows that shoot tips stopped growth at the end of September under natural season and have enlarged vessel cells and increased xylem in the cambia zone of apical shoot 34 . The phenotype at stage S4 here was similar to that at the end of September. Therefore, the S4 stage should have growth cessation and could be at dormancy stage. The highest number of differentially expressed genes was approximately 20 days in our study under SD, which is comparable to the reported 21 days under SD in grapevine 5 . This indicated that 2-3 weeks may be enough to trigger most of the gene changes, which leads dormancy induction. A previous study showed that the leaf primordia in the apical meristem will grow into young leaves and then will die very soon, and the naked growth cone will gradually die too in the same cultivar during the dormancy processes 40 . Similar changes were observed in our study, which also supports SD-induced dormancy. The terminal shoot tips will die in the end, but we did not find any DE gene related to   www.nature.com/scientificreports www.nature.com/scientificreports/ programmed cell death (PCD) in this study, which may be interpreted by our SD treatment at the fast-growing season or the induced cessation still at an early stage before PCD symptoms appeared.
Genetic and genomic studies on Paulownia lag far behind those on other plants, even though there have been several reports on gene expression relating to drought and hybridization of Paulownia 41,42 . Here, we enriched the transcripts and their gene expression in a time-course series, as well 491 commonly shared DE genes across all stages under SD, which will be the key set of gene regulated by SD induction although stage-specific genes may play roles too. The induced 116 shared pathways provide detailed regulation.
SD induces dormancy in perennial plants by down-regulating CO and FT, and up-regulating PIF3/PIF4 1 . The diurnal circadian rhythm plays crucial roles in the timing of gene expression to adapt to environmental conditions 43 . Our study found regulation of several circadian genes in SD-induced growth cessation and dormancy. The regulation occurred through genes in both the far red/red light and blue light signalling pathways and through interaction with circadian oscillators. These included significant changes (>2-fold) of five genes. PIF3 was increasingly up-regulated in the far red and red signalling pathways. SPA1 and FKF1, which control photomorphogenesis and interact with circadian oscillators in blue light, respectively, were regulated dynamically. PRR5 (pseudo-response regulator 5) is a transcriptional repressor for regulating the clock 44 , and its gene expression was up-regulated at all stages in the RNA-Seq analysis and, to our best knowledge, is the first observation of PRR5 regulation under SD. The protein MYB, known to encode a circadian clock regulator 45 , was down-regulated through stages S2 to S4 and up-regulated at S5 after long SD treatment. We did not observe regulation of CO expression which was expected to be at low level at day time 1 when the samples was collected, but we found that FT was regulated at late stage only. This may have been caused by the very young status of the one-year-old www.nature.com/scientificreports www.nature.com/scientificreports/ Paulownia plants we investigated here. A transcriptome analysis from SD to LD shows the FT is upregulated in Medicago 46 . Therefore, FT gene could be one of the key genes for dormancy induction under SD in apical shoot tips, which is of further interest in gene manipulation during breeding improvement.
Carbon fixation metabolism was also affected the growth cessation and dormancy induction in Paulownia. We found many DEGs involved in three processes of carbon fixation metabolism and genes in C4 model of carbon fixation in Paulownia. However, the Kranz anatomy of C4 plant was not present in Paulownia 47 . This is also reported that the Kranz anatomy is not essential in C4 plants 48 . The presence of expressed genes in both the C4 and CAM cycles suggest that Paulownia has evolved to separate initial CO2 fixation, minimize photorespiration and save water. The CAM is known to fix CO2 in the night which is beneficial to close stomata in the daytime to reduce water loss via transpiration 16,17 . Therefore, Paulownia may fix CO2 via CAM to adapt to hot and dry environments, which is good for the large leaf area of Paulownia. Leaf is the major tissue for photosynthesis while shoot tips should be minor tissue compared with leaves. However, the shoot tips of Paulownia are bigger than other plants and green, which may also conduct strong carbon fixation. In addition, many detected DEGs in photosynthesis suggested that the shoot tips do have photosynthesis and were affected. Whether these DEGs is the cause or the result of growth cessation is worthy of further investigation with the candidate genes identified here. Figure 6. The differential gene regulation in carbon fixation of photosynthesis. Image a shows the metabolism of carbon fixation to which the differentially expressed genes (shaded) were mapped. The pathway was modified from map MAP00710 (version 2013) from KEGG database (http://www.kegg.jp). Image b shows the expression patterns of differential expressed genes encoding corresponding enzymes in pathway before (S1) and after (S2-S5) short day treatment. Image c shows the relative gene expression level detected with RT-qPCR. The pearson coefficient between RNA-Seq result and RT-qPCR results is 0.5 and 0.5 for MDH and PEPC, respectively. (2019) 9:16602 | https://doi.org/10.1038/s41598-019-53283-2 www.nature.com/scientificreports www.nature.com/scientificreports/ Consistently, previous report shows that the carbon efflux was decreased at autumn in deciduous plants 19 . Thus, the photosynthesis is definitely associated with the growth cessation. During SD induction, the expression of the key gene PEPC was up-regulated, indicating an even higher efficiency of CO 2 capture in the shoot tips. However, the gene GAPDH, which is responsible for converting released CO 2 from malate to G3P, was down-regulated in the Calvin cycle. This suggested a high CO 2 input but a low G3P output. Taking these data together, the extra captured carbon in the form of malate may be converted to another compound instead of going to the Calvin cycle under SD. The expression of RUBISCO, the most important gene in carbon fixation, was not changed in our study, although it was reported to be reduced in ABA-induced dormancy in Spirodela 21 . These may be explained by the fact that higher CO 2 -concentrating activity may suppress the decrease of RUBISCO expression in our study. Of course, posttranscriptional and translational modification will also affect the final activity of Rubisco.
Research efforts have investigated the role of hormones in inducing and maintaining dormancy and found that hormones often crosstalk with each other 7,8,12,15 . Microarray expression analysis showed genes in four hormones, including GA, ethylene, auxin and ABA, were regulated during bud dormancy in grape 9 . In our research, seven hormone pathways were significantly affected by SD treatment. The auxin signalling pathway showed the most significant changes because expression of six transcripts encoding AUX/IAA, ARF and SAUR genes was changed. Here, transcripts of one IAA and two copies of ARF were down-regulated and showed a dramatic downward trend at 15 days, which is consistent with those in previous transcriptomic studies during growth cessation and bud dormancy in other plants 6,7 , meaning common regulation in auxin pathways. SAURs, a gene family, are effectors that integrate environmental signals in the auxin pathway for plant growth by inhibiting ATPase and cell expansion, integrating other hormones, such as ABA, GA, Jasmonate, and Brassinosteroids (BR), to regulate plant growth and development 37 . In our results, three copies of SAURs displayed differential expression, meaning that they may have distinct regulation roles during plant development, as had been reported in a previous study 37 . Cytokinin synthesis was reported to be inhibited by auxin 9 . Two transcripts of CRE genes were differentially regulated. Taking these results together, we propose that the auxin signalling pathway has key roles in growth cessation and bud dormancy under SD in Paulownia.
Here our results show that genes in signalling of additional hormones, including BR, GA, JA, ABA and ethylene, were also regulated in growth cessation and bud dormancy. BR functions to promote stem elongation and cell division. Here, the up-regulated BZR1 in SD treatment is a key transcription factor activated by BR. BZR1 can interact with DELLA family proteins, which inhibit the GA signalling pathway 49,50 . Thus, BR signalling and GA signalling crosstalk under SD. JAR1, up-regulated in SD treatment, catalyses the synthesis of JA-amido conjugates, including the most active form, JA-Ile 51,52 . ABA is generally believed to play a positive role in seed and bud dormancy in woody plants 12 , and SnRK2 is known to play an important role in ABA signalling 53 . Therefore, the up-regulated SnRK2 in our study indicated an increase in ABA under SD to promote dormancy. Ethylene counteracts ABA signalling pathway and releases dormancy 15 . Here, the ethylene signalling pathway inhibitor encoding gene EBF1 was up-regulated during SD treatment, which may be required for maintaining high ABA for the dormancy status. Together, many hormone-related pathway genes work together to induce and maintain growth cessation and dormancy under SD in Paulownia.
In summary, our RNA-Seq-based transcriptome analyses clarified that many genes were regulated in multiple pathways and metabolisms in apical shoot tips during SD-induced growth cessation and dormancy induction in Paulownia. The regulation of key genes in the circadian rhythm network, such as PIF3 and PRR5, plays a role in environmental signal integration in dormancy induction. Expression of genes in all known plant hormones, especially IAA, ARF and SAUR in auxin signalling, was altered. Under SD, increased expression of PEPC was not coupled with higher assimilation to G3P in photosynthetic carbon fixation pathway. We summarized the discovery in this study in Fig. 7, where expressional regulation on key DE genes integrates the SD signal and leads to dormancy www.nature.com/scientificreports www.nature.com/scientificreports/ control in Paulownia. Our basic understanding of expressional regulation under SD in Paulownia provided here may help to guide dormancy management via editing key gene's regulation for breeding higher quality timber.

Materials and treatment.
To investigate short day effects on terminal growth of shoot tips, plant clones were generated from the Paulownia cultivar P. tomentosa roots buried in natural soil inside a greenhouse during the previous winter. The plants were grown until the beginning of Aug 2015, and the plants were chosen for this study when they were approximately 50-60 cm high with one healthy apical bud. Forty plants, eight plants for each sampling stage, were treated with short-day light cycles of 8 h light and 16 h darkness, while the natural day light cycle was 16 h light and 8 h darkness during the experiment. Approximately 0.5 cm of the shoot tip from each plant was sampled just before the treatment on Aug 1 st as a control, and then, the same tissues were sampled every five days (Aug 6, 11, 16 and 21). The samples were collected between 9 and 10 o' clock in the morning. The temperature was 27.8 °C, and the humidity was 65-85%. All plants were watered appropriately. The plants were grown in the greenhouse at Henan Agricultural University (GPS location 34.798 N, 113.650 E). Samples were then immediately frozen in liquid nitrogen and store at −80 °C until being processed.
RNA extraction and sequencing. Shoot tips from the same sampling date were pooled and ground into powder with a mortar and pestle in liquid nitrogen. Total RNA was extracted from 80-120 mg powder for each sample using kit (TIANGEN, catalog # DP441) and checked in 1% agarose gel and then Agilent 2100 Bioanalyzer followed by subsequent Illumina library preparation as described previously 54 . Paired-end 125 bp sequencing was performed for each library on a HiSeq. 2500 machine.
Transcript assembly, annotation and expression analysis. Before transcript assembly, raw reads were cleaned. The removal of poor-quality sequences and trimming of adaptor sequences were carried out using cutadapt (v1.8.1). The clean reads were merged and assembled into transcripts using Trinity software (version 201308) 55 with the default parameters except with the parameters-min_kmer_cov 3-min_glue 3. All transcripts were then clustered by cdhit (v4.6) with default parameters 56 . Then, the clustered transcripts were aligned to the Rfam database (v.11; http://www.sanger.ac.uk/Software/Rfm/) to exclude reads corresponding to rRNAs. The software package RSEM was adopted to quantify total expressed transcripts in FPKM (fragments per kilobase of exon model per million mapped reads) 57 . The read count was used for differentially expressed genes analysis with EdgeR package (version 3.14) with settings as at least two-fold change and qvalue < 0.0001 58 . The expression patterns of genes across sampling stages were clustered with mfuzz package (version 2.44) 59 . Unigenes were defined as the longest sequences in the assembly cluster called component in Trinity 55 . Unigenes were then searched against the databases of NR, KEGG and EggNOG with BLAST with E-value threshold 1E-5. Gene function analysis was predicted using InterProScan software (version 5.8-49) 60 , and enrichment analysis was done using WEGO 61 and Blast2GO software with default settings 62 . KEGG mapping and enrichment were analysed using KOBAS with default settings 36 and annotated with the Automatic Annotation Server at the KEGG website (www. genome.jp/kegg/) using default settings against Arabidopsis and Populus species 63,64 . The permission to publish KEGG pathway map images of Circadian rhythm, Carbon fixation in photosynthetic organisms and Plant hormone signal transduction is granted by Kanehisa Laboratories.
Real-time quantitative (Rt-qpcR) analysis. RT-qPCR analysis was conducted for selected DEGs to validate the expressional patterns mined from RNA-Seq analysis. Total RNAs were extracted from shoot tips of triplicate biological plants with a RNAprep Pure kit (TIANGEN, catalog # DP441) and cDNA was synthesized from mRNA with FastQuant RT kit (TIANGEN, catalog # KR106) as described previously 54 . Three independent experiments were conducted. Primers are designed from DEGs with primer3 and provided in supplementary file (Supplementary Table S11). The primers efficiency in RT-qPCR was optimized and efficiency was tested in dilution series of 1, 1/10, 1/1000 and 1/10000 (Supplementary Table S11, Supplementary Fig. S3). RT-qPCR analysis was conducted for selected DEGs using SYBR green master according to previous methods except annealing temperature set to 60 °C here 65 . 18 S and H 2 O was used as reference gene and blank controls. The relative expression level was calculated with 2-ΔΔCt method 66 . The relative expression levels were presented as mean value ± standard error (S.E.) relative to S1. Therefore, the RT-qPCR results show the trend of expression pattern, but not the absolute level. Different letters indicate statistically significant differences to Duncan's multiple range tests with p < 0.05. The coefficient between RNA-Seq result and RT-qPCR result was calculated with Pearson correlation.

Data availability
The RNA-Seq reads and transcriptome assembly have also been deposited at the Short Read Archive and DDBJ/ ENA/GenBank under the Bioproject accessions PRJNA393208 and GFVK00000000, respectively. The version described in this paper is the first version, GFVK01000000.