Deep sequencing on a genome-wide scale reveals diverse stage-specific microRNAs in cambium during dormancy-release induced by chilling in poplar

Background Trees in temperate zones show periodicity by alternating active and dormant states to adapt to environmental conditions. Although phytohormones and transcriptional regulation were found to be involved in growth cessation and dormancy transition, little is known about the mechanisms of the dormancy-active growth transition, especially dormancy maintenance and release. Small RNAs are a group of short non-coding RNAs regulating gene expressions at the post-transcriptional level during plant development and the responses to environmental stress. No report on the expression profiling of small RNAs in the cambial meristem during the dormancy-active growth transition has been reported to date. Results Three small RNA libraries from the cambium of poplar, representing endodormancy induced by short day conditions, ecodormancy induced by chilling and active growth induced by long day conditions, respectively, were generated and sequenced by Illumina high-throughput sequencing technology. This yielded 123 known microRNAs (miRNAs) with significant expression changes, which included developmental-, phytohormone- and stress-related miRNAs. Interestingly, miR156 and miR172 showed opposite expression patterns in the cambial dormancy-active growth transition. Additionally, miR160, which is involved in the auxin signaling pathway, was expressed specifically during endodormancy release by chilling. Furthermore, 275 novel miRNAs expressed in the cambial zone were identified, and 34 of them had high detection frequencies and unique expression patterns. Finally, the target genes of these novel miRNAs were predicted and some were validated experimentally by 5′RACE. Conclusions Our results provided a comprehensive analysis of small RNAs in the cambial meristem during dormancy-release at the genome-wide level and novel evidence of miRNAs involved in the regulation of this biological process. Electronic supplementary material The online version of this article (doi:10.1186/s12870-014-0267-6) contains supplementary material, which is available to authorized users.


Background
Trees in temperate zones show periodicity by alternating active and dormant states to adapt to natural conditions, such as light, temperature and drought. Growth arrest is the first step of plant dormancy, followed by the dormant state, which can be divided into two stages: ecodormancy (or quiescence) and endodormancy (or rest). In ecodormancy, plants can restore active growth upon exposure to growth-promoting conditions, while plants in endodormancy cannot [1].
Photoperiod has been known to govern the growth cessation of many trees in temperate climates. Plants sense changes in the photoperiod through the leaves and send a graft-transmitted message to the terminal, so that the terminal initiates dormancy [1,2]. The identification of poplar FLOWERING LOCUS T (FT) and CONSTANS (CO) as mediators of growth cessation induced by the short day (SD) photoperiod was a significant breakthrough in the study of dormancy transition regulations [3]. Like-AP1 (LAP1), a poplar ortholog of the Arabidopsis gene APETALA1 (AP1), mediates the photoperiodic control of seasonal growth cessation downstream of CO/FT [4]. Another environmental factor that controls dormancy transition is temperature. Low temperature plays an important role in inducing growth cessation and dormancy [5,6]. However, a continuous chilling must occur to release endodormancy and switch to ecodormancy, and then warm temperatures in the spring subsequently reinitiate growth [6].
Environmental factors are thought to regulate the precise annual cycle's time course by modulating phytohormone levels or altering the sensitivity of the cells to phytohormones. So far, gibberellins and auxin are widely recognized as the most important phytohormones involved in the dormancy transition [7][8][9][10][11]. Applications of exogenous gibberellins could cause the dormant poplar buds to sprout without chilling, and it is shown that a low temperature could alter the expression of key regulators in the gibberellin signal pathway [11,12]. Auxin has crucial roles in cambial cell division, which makes it very important in the dormancy-active growth transition [13,14]. A recent study shows that the induction of cambial growth cessation and dormancy involves changes in auxin responses rather than auxin content [7]. Another phytohormone that may participate in the dormancyactive growth regulation is abscisic acid (ABA), which peaks in poplar apical buds after growth cessation and before bud set [15][16][17].
So far, the studies on mechanisms of the cambial dormancy-active growth cycle have mainly focused on hormonal [9,10,12,15] and transcriptional regulation [16][17][18]. The switch between plant dormancy and active growth is a complex biological phenomenon that involves a large number of genes and many metabolic processes, as well as the interactions of a variety of hormones. Multiple levels of control networks are involved in such complex biological events, in addition to transcriptional and protein regulation.
Small RNAs (sRNAs), short (~21 nt) non-coding RNAs, are important regulators of gene expression at the posttranscriptional level during plant development and response to environmental stress [19]. sRNAs, in particular microRNAs (miRNAs), have been studied extensively in poplar, including genome-wide profiling of sRNAs and miRNAs [20,21] and stress responses to drought [22,23], salt [24], cold [25] and pathogens [26]. In addition, some miRNAs have been found to be of great importance in tree development. For instance, miR166 is reported to be involved in vascular tissue development [27,28] and may be related to the cambial active period [29]. MiR156 and miR172, which is well studied in Arabidopsis, appear not only to control flowering and the timing of sensitivity in response to vernalization, but also vegetative phase changes in trees [30][31][32][33][34][35]. Although comprehensive work has been done to describe miRNAs in trees during various cellular processes, there is no report on the expression profiling of miRNAs in the cambial meristem during the dormancy-active growth transition and little is known about the regulation of miRNAs in the process. In this paper, we present a deep sequencing profile on a genome-wide scale that reveals stage-specific miRNAs in the cambial zone during this process. Millions of sRNA reads were obtained, and after further analysis, we found 123 known miRNAs, including developmental-, phytohormone-and stress-related miRNAs, which showed significant expression-level changes during dormancy-release by chilling. Furthermore, 275 novel miRNAs expressed in the cambium zone were identified, and 34 of them had high detection frequencies and unique expression patterns. The target genes of these novel miRNAs were predicted and some of them were validated. Our results revealed the expression changes of miRNAs in cambium dormancy-release by chilling in poplar, and provided evidence of miRNA involvement in the regulation of the dormancy-active growth transition of trees.

Dormancy-active growth transition induced by photoperiod and chilling in poplar
The induction of dormancy and resumption of growth in poplar were constructed according to Espinosa-Ruiz et al. [36] with some modifications. After 8 weeks of the short day (SD) treatment of 8 h light/16 h dark, the tree growth was arrested. Dormant apical buds formed (Figure 1a, c) and the layers of cambial cells (Figure 1b, d) decreased from 6-8 to 1-2 ( Figure 1h). Although trees were transferred to the long day (LD) condition of 16 h light/8 h dark at this time, they would not resume growth, indicating their endodormant state. To release endodormancy, the 8-week SD-treated trees were exposed to chilling temperatures of 4°C. Only trees exposed to a chilling treatment for at least 4 weeks could resume growth, which was shown in bud burst, and cambial cell division and differentiation, when they were transferred to LD conditions at 25°C for 3 weeks (Figure 1e, f, g, i). The results indicated that the endodormant state was released and that the trees had shifted to the ecodormant state after 5 weeks of the chilling treatment. Then, the active growth state was induced by 3 weeks of the LD condition at room temperature.

Deep sequencing of sRNAs in cambium during dormancy-release in poplar
To investigate the miRNA component of sRNAs and the changes of miRNAs in cambial meristem during dormancy-release in poplar, three sRNA libraries from the cambium of poplar, representing endodormancy with an 8-week SD treatment (SD8), ecodormancy with a 5-week chilling treatment (C5) and active growth under a 3-week LD condition (LD3) after chilling, respectively, were generated and sequenced by Illumina high-throughput sequencing technology. Raw read totals of 16,688,990, 21,379,082 and 15,942,869 from SD8, C5 and LD3, respectively, were acquired. After removal of low-quality sequences, adapter sequences, polyA sequences, sequences smaller than 18 nucleotides and other artifacts, we obtained 16,339,437, 20,887,480 and 15,649,238 high-quality 18 to 30 nt sRNA clean reads in SD8, C5 and LD3, respectively, for further analysis (Table 1).
Among the 18 to 30 nt sRNA clean reads from sequencing, the majority of them (65%) were in the range of 20 to 24 nt in length, with sequences of 21 nt or 24 nt representing the most abundant classes in each library ( Figure 2). The major component of the sRNAs in SD8 and C5 was 21 nt long; however, the proportion of 24 nt sRNAs peaked in LD3 ( Figure 2). sRNA libraries generated by sequencing were complex in composition, including miRNAs, siRNAs, rRNAs, tRNAs, small nuclear RNAs (snRNAs) and small nucleolar RNAs (snoRNAs). To annotate the sRNAs, we first mapped the sRNAs of 18 to 30 nt to the Populus trichocarpa genome (www.phytozome.net) using SOAP software (http://soap. genomics.org.cn), and then characterized each kind of sRNA by aligning them to certain databases. Known miRNAs were identified by alignment to sequences in miRBase 20.0 with no mismatches. Meanwhile, the Rfam9.1, NCBI and Gen-Bank databases were employed to annotate the other kinds of sRNAs, including scRNAs, rRNAs, tRNAs, snRNAs and snoRNAs. The repeats that represented the sRNAs positioned at repeat loci were identified using Tag2repeat software. In addition, there were possibly degraded species of mRNAs in the sRNA libraries, which were determined through intron/exon alignment. The remaining unannotated sRNAs were candidates for predicting novel miRNAs and potential miRNA seeds edit. As a result, 573,822, 618,526 and 844,787 unique sRNAs in SD8, C5 and LD3, respectively, were mapped perfectly to the genome, and the proportions for each kind of sRNA were listed in Table 2. Interestingly, the miRNAs represented 22.68% and 24.92% of the total sRNA reads in SD8 and C5, respectively, but only 13.45% in LD3. There were~200 more unique miR-NAs in LD3 than in both endodormancy and ecodormancy ( Table 2), indicating that the miRNA population in active cambium was more diversified, which may be due to the complex cellular processes associated with active growth.

Identification and expression profiles of known miRNAs in cambial meristem during dormancy-release in poplar
Known miRNAs in the cambium of poplar were annotated by alignment to the sequences in the available poplar miRNA database. As a result, we identified 182 mature miRNA, two miRNA-5p, two miRNA-3p and 183 pre-miRNAs in SD8, 176 mature miRNA, two miRNA-5p, two miRNA-3p and 177 pre-miRNAs in C5, and 175 mature miRNA, two miRNA-5p, two miRNA-3p and 176 pre-miRNAs in LD3. All the mature miRNAs identified belonged to 33 conserved and non-conserved miRNA families, of which 123 known miRNAs in 26 miRNA families showed significant expression-level changes during this process (Additional file 1: Table S1).
To elucidate the potential regulatory roles of these miR-NAs in the dormancy-active growth transition, we analyzed the miRNAs with unique expression patterns during the process, which were mainly involved in plant development and stress response, as well as the plant hormone signal pathway.
In our dataset, eight differentially expressed known development-related miRNA families were detected, including miR164, miR396, miR168, miR319, miR171, miR166, miR156 and miR172 ( Figure 3). These miRNAs  functioned in cell proliferation (miR164, miR396 and miR319) [37][38][39][40], vascular development (miR166) [41] and miRNA biogenesis (miR168) [42]. Most of these development-related miRNAs were enriched in the active growth stage. miR319 increased dramatically from SD8 to C5 and continued at a high expression level from C5 to LD3, suggesting that the expression of miR319 may be affected by the chilling treatment. However, miR164, miR168 and miR396 showed no obvious, or only slight, changes from SD8 to C5, but increased in LD3. Intriguingly, unlike other development-related miRNAs, miR166 was enriched in SD8 and C5, and was nearly undetectable in LD3. The members of the miR171 family showed different expression patterns; some of them were highly expressed, while some were repressed during the active growth, indicating that members from one family could have distinct functions in this process.  Figure 3 The fold change of development-related known miRNAs during dormancy-release in poplar. The y-axis represented the fold change (log2 value) of normalized miRNA counts. SD8 was arbitrarily set to be the control of C5 and LD3. SD8: short day treatment for 8 weeks; C5: chilling treatment for 5 weeks; LD3: long day treatment for 3 weeks. miR156 and miR172 are well known for controlling the meristem cell fate transition in maize [30][31][32][33], Arabidopsis [34] and the vegetative phase change in trees [35]. In our study, 10 miRNA members of the miR156/157 family and six miRNA members of the miR172 family were identified. Intriguingly, miR156 was highly expressed in SD8 and C5, and then decreased in LD3, while miR172 had the opposite expression pattern (Figure 3), which showed a similar expression pattern during the vegetative phase change in trees [35].
To investigate miRNAs involved in the process through the plant hormone pathway, the dynamic expression levels of hormone-related miRNAs were analyzed. Auxin signaling-related miR160, miR167 and miR390 had distinct differential expression patterns in the dormancyactive growth transition ( Figure 4). The expression of miR160 peaked in C5, which was the phase sensitive to auxin treatment in dormancy. The unique enrichment in ecodormancy suggested miR160 had an important role in the transition from endodormancy to ecodormancy. Unlike miR160, miR167 and miR390 maintained low expression levels from SD8 to C5, and then increased dramatically from C5 to LD3 (Figure 4), indicating that miR167 and miR390 may function in the auxin pathway during active growth.
Several members of the miR169 family, whose target genes participated in ABA resistance [43,44], were identified in the cambium during the dormancy-active growth transition. The expression levels of all the miR169 members stayed basically unchanged during SD8 and C5, while in LD3, they displayed opposite trends ( Figure 4). These findings showed that members of the miR169 family had different functions in this process.
The miR159 family repressed the conserved GAMYBlike genes that have been implicated in gibberellin (GA) signaling in anthers and germinating seeds [45]. We found that miR159 was highly expressed in C5 and kept rising in LD3 (Figure 4). The GA signal had already been proven to be a key factor during dormancy-release in poplar [11,12]. The expression change of miR159 raised the possibility that it may be involved in this mechanism through the GA signal pathway.
Lu et al. identified 68 stress tolerance-related miRNAs in poplar [46]. Among them, miR472, miR475, miR477, miR1444 and miR1446 were found to show differential expression levels in this study ( Figure 5). The abundance of miR1444 greatly dropped, while those of the others changed slightly from SD8 to C5. All five of these miR-NAs had lower expression levels in LD3 ( Figure 5).

Identification and expression profiles of novel miRNAs
The Mireap software was employed to screen novel miR-NAs from candidates by exploring not only the secondary hairpin structure, but also the Dicer cleavage site and the minimum folding free energy (MFE). According to the analyses, more than 80% of the candidate novel miRNAs in SD8 and C5 began with a 5′uridine, which was a conserved feature of miRNAs recognized by the ARGO-NAUTE 1 (AGO1) protein [47]. However, in LD3, this ratio was reduced to~50%. To ensure the authenticity of novel miRNAs, several conditions must be satisfied: the lengths of the mature candidate novel miRNAs varied from 20 to 23 nucleotides, the number of reads was greater than five, and all the unique sequences were identified in at least one library. As a result, we obtained 128, 110 and 147 unique sequences in SD8, C5 and LD3, respectively. After removing the redundancy, 275 novel miR-NAs were identified in the three libraries. The average MFE value of novel miRNAs in each library was −55.37 ± 21.50 kcal/mol in SD8, −50.39 ± 19.03 kcal/mol in C5 and Figure 4 The fold change of phytohormone-related known miRNAs during dormancy-release in poplar. The y-axis represented the fold change (log2 value) of normalized miRNA counts. SD8 was arbitrarily set to be the control of C5 and LD3. SD8: short day treatment for 8 weeks; C5: chilling treatment for 5 weeks; LD3: long day treatment for 3 weeks. −53.75 ± 21.72 kcal/mol in LD3. Most of these novel miR-NAs were only expressed at a specific stage, and only a few of them were expressed in two of three libraries. We found 55 specific unique sequences in SD8, 50 in C5 and 92 in LD3 (Additional file 2: Table S2). Most of these novel miRNAs had low detection frequencies in all three libraries. Here, we listed the miRNAs whose detection frequencies were greater than 20 in at least one library and that had marked expression-level changes (Table 3).
To validate the predicted novel miRNAs and confirm the expression profiles determined by Illumina highthroughput sequencing technology, we performed quantitative real-time PCR (qRT-PCR) on a subset of six miRNAs sequences, including two conserved and four novel miRNAs from SD8, C5 and LD3 ( Figure 6). Most of the expression patterns were in agreement with our sequencing data, while a few miRNAs did not show the same expression trends. For example, the expression level of A-m0126_3p in C5 was measured to be higher by qRT-PCR than by the sequencing reads, which may be caused by a lack of sequence depth.

Prediction of novel miRNA targets and RACE validation
A web-based miRNA target prediction program was employed to hunt for potential miRNA target genes. A total of 763 unigene sequences were predicted to be the targets of 119 novel miRNAs in SD8, 942 unigene sequences to be the targets of 107 novel miRNAs in C5 and 833 unigene sequences to be the targets of 126 novel miRNAs in LD3. The number of predicted targets varied from 1 to 34 per miRNA and most had three to seven targets. To focus on the biological processes, we predicted the targets of novel miRNAs that were specifically expressed in one phase or that had expression changes during the phase transition (Table 4). Although the target genes of some of these miRNAs showed distinct functions, a portion of them predicted the target as a single gene or members of a gene family. Many of these targets were involved in energy metabolism and solute transport, representing the dramatic metabolism changes between dormancy and active growth. Several targets were annotated as NBS resistance protein and leucine rich repeat protein, suggesting the direct response to adverse environment. Additional, some novel miRNAs targeted cell signaling-related genes, which could lead to expression change of these genes in the annual cycle.
To validate the cleavage events of novel and known miRNAs, a modified RNA ligase-mediated rapid amplification of cDNA ends (RLM-RACE) experiment was performed to verify the miRNA-guided mRNA cleavage events. We tested four novel miRNAs and two known miRNAs to verify their ability to cleave their targets. All of the cleavage sites were located between 10 and 11 nucleotides relative to the 5′end of the complementary miRNA sequence, which was the characterized cleavage site of almost all of the known miRNAs (Figure 7). The RACE products of miR156 and miR172 were cloned and sequenced. Their alignments to the poplar genome showed that the targets of miR156 and miR172 were homologs of the DNA-binding transcription factors SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) and APETALA2 (AP2), respectively, and the targets of other novel miRNAs were identical with the computer prediction in Table 4.

Discussion
Plant miRNAs have a wide range of regulatory functions in many biological and metabolic processes, including developmental regulation, cell differentiation, signal transduction, growth control, and biotic and abiotic stresses [40]. Although an increasing number of poplar miRNAs have been identified in tissues or under certain environmental conditions [48], and some of them have been well characterized to involve various developmental process [35,49], little is known about the roles of miRNAs in the cambium dormancy regulation in trees. We have presented here a comprehensive analysis of sRNAs in the dormancy-active growth transition at the genome-wide level, which revealed dynamic features of sRNA populations in the annual growth cycle and expression patterns of miRNAs involved in this process. In addition, a set of novel miRNAs with notable expression pattern changes was identified. Together, these results provide novel insights into the regulatory mechanism of the dormancyactive growth transition mediated by miRNAs.

Deep sequencing reveals a diverse set of sRNAs in the cambium of poplar
Using high-throughput sequencing technology, we obtained more than 3 million unique sRNAs reads from three cambium samples during the dormancy-active growth transition in poplar. Although sRNAs are complex in composition, the large majority are 21 nt and 24 nt in plants [19], and the proportion of miRNAs varies between different species and upon environmental conditions [22,50,51]. The 24 nt sRNAs were mainly composed of siRNAs associated with repeats and transposons [52]. In our case, the sRNA length distribution patterns diverged during the dormancy-active growth transition. In dormancy, including endodormancy and ecodormancy, the 21 nt long sRNAs constituted the most abundant class, while in active growth the 24 nt long sRNAs constituted the most abundant class. We determined the size distributions in previous studies in poplar, and found that the 21 nt sRNAs were the major component in leaves and vegetable buds [20,48], while the xylem tissue has a major peak at 24 nt [48], which was in agreement with our data during active growth. We also found that the proportion of total miRNAs in dormancy, including both endodormancy and ecodormancy, was greater than that in active growth, which confirmed the induction of the 21 nt miRNA by dormancy. The increase in 24 nt sRNAs during active growth suggested that the 24 nt sRNAs, which would be mainly siRNAs known to guide DNA methylation and heterochromatin formation [53], may participate in the regulation of cambium activity, including cell division, cell differentiation and phytohormone regulation. The reversal of 21 nt and 24 nt sRNA abundance in the dormancy-active growth transition also indicated that these two kinds of sRNAs may play different roles during the annual growth cycle.

Unique expression patterns of miRNAs in dormancyrelease in poplar
Hundreds of miRNAs have been surveyed in poplar since next generation sequencing technology has become widely used, but little is known about miRNAs in tree dormancy regulation, especially the transition between dormancy and active growth. In this study, we found a series of miR-NAs that might be involved in this process. for instance, most of the developmental-related miRNAs, especially those involved in meristem activity or cell proliferation, presented specific expression patterns. In Arabidopsis, increasing evidence shows that miR164, miR319, miR396 and their targets form a miRNA regulatory network to regulate cell proliferation, leaf development and meristem activity [40,54,55]. Intriguingly, these three miRNAs showed similar expression pattern during the endodormancy release process. The increasing expression levels of these three miR-NAs in active growth suggested that a miRNA network regulating cell proliferation also existed in active cambium cells. Considering that both cambium and the leaf primordium are capable of cell division, the high level of these three miRNAs in active growth is quite reasonable, and the cessation of cell division in dormancy may cause the low abundance of these three miRNAs. Another miRNA that is crucial for vascular development is miR166, which regulates the class III HOMEODOMAIN-LEUCINE ZIPPER family of transcription factors. The relationship between miR166 and its target is essential for leaf abaxial/adaxial polarity establishment [40,56]. Unlike other developmental-related miRNAs, miR166 was more abundant in dormancy, and  Plant invertase/pectin methylesterase inhibitor had a very low expression level in active growth, which was in agreement with a previous study in poplar [29]. These results indicated that miR166 was down-regulated in active growth to increase the expression level of its target gene, which had an important role in vascular development. In addition, the original functions of these miRNAs were mainly found in the shoot apical meristem or leaves, thus the existence of these miRNAs in cambium suggested they may share the same regulatory mechanism in different tissues.
Interestingly, miR168 was found to be up-regulated in active cambium. The target of miR168 was AGO1, which was the key regulator of miRNA biogenesis [42]. The high expression level of miR168 leads to the repression of AGO1, which causes a reduction in the miRNA expression level. As expected, the total miRNAs in active growth decreased, indicating that miR168 was involved in the miRNA biogenesis as a feedback regulator in the cambium stage transition. miR156 and miR172 target DNA-binding transcription factors SPL and AP2 genes, respectively, which control the juvenile-to-adult vegetative transition both in annual herbs [57,58] and woody perennial plants [35]. They show converse expression patterns and regulatory relationships during the phase transition [57]. Surprisingly, the expression levels of miR156 and miR172 also had opposite expression patterns in our study. The 5′RACE results confirmed the cleavage events in miR156 and miR172, suggesting that the two miRNAs were functional during the dormancy-active growth transition. These similar expression patterns suggested the complementary of miR156 and miR172 might play an important role in this process, which needed to be experimentally confirmed.

Auxin-related miRNAs may participate in the regulation of endodormancy release
A continuous chilling is the only natural way to release endodormancy and transition to ecodormancy. The main difference between endodormancy and ecodormancy is that ecodormant trees have the ability to respond to growth-promoting signals, such as auxin or appropriate outside conditions. In other words, the chilling triggers the ability of the tree to respond to auxin. In our data, most miRNAs had significant expression-level changes between dormancy and active growth, but only a few of them had changes between endodormancy and ecodormancy. Among them, miR160, whose target was AUXIN RESPONSE FACTOR 10/16/17 (ARF10/ARF16/ARF17) [59][60][61], was highly expressed in ecodormancy. Lu et al. studied the cold-responsive miRNAs in poplar by microarray analysis, and showed that ptc-miR160a-g were strongly induced by cold treatment for 12 h, and that this induction disappeared after 16 h of treatment [46]. In our case, the chilling lasted 5 weeks, so this high abundance of miR160 in ecodormancy may not be due to cold tolerance, but ecodormancy itself. In Arabidopsis, miR160's target genes are negative regulators of auxin signaling [62,63], so it is possible that the highly expressed miR160 may enhance the auxin signal by repressing its targets. The other two auxin-related miRNAs, miR390 and miR167, increased dramatically during active growth, suggesting that they are involved in the auxin signal pathway during active growth. The results showed that the miRNAs mediating the auxin signal pathway had complex regulatory roles in the cambium dormancy phase transition.

Novel miRNAs and their putative targets during the dormancy-release in poplar
Hundreds of novel miRNAs as well as their targets were identified in this study. Some of them may have important roles in the dormancy-activity transition. For instance, 0A-m0062_3p which had a high expression level in active growth targeted a histone acetyltransferase gene. Considering the close link between histone acetylation and gene activation [64], this result suggested that some changes of gene expression in endodormancy could be caused by the regulation of the histone-modifying enzyme by miRNAs. Additional, one putative target of an active growth-specific novel miRNA, A-m0165_5p, was annotated as a homolog of DOMAINS REARRANGED METHYLTRANSFERASE 2, which catalyzes de novo methylation and is responsible for RNA-directed DNA methylation in Arabidopsis [65]. Together, these findings might suggest a miRNA-guided regulation of epigenetic modifications in dormancy and active growth. Among dormancy-highly-expression novel miRNAs, 0A-m0077_5p, which pairs with an ubiquitinconjugating enzyme and 0A-m0149_5p, which matches a plant invertase/pectin methylesterase inhibitor, were also detected, suggesting different protein levels and cell wall components in annual cycle. Compared with the higher numbers in active growth and endodormancy, only four novel miRNAs with high expression levels in ecodormancy are listed in Table 4. Among them, 5A-m0010_5p targeted a homolog of the No Apical Meristem (NAM) protein, suggesting that this novel miRNA may participate in meristem activity. Interestingly, the target of 5A-m0104_3p with the highest expression in ecodormancy was RACE validated and annotated as a homolog of 6phosphogluconate dehydrogenase, which is an oxidative carboxylase that catalyzes the decarboxylating reduction of 6-phosphogluconate into ribulose 5-phosphate in the presence of NADP. Since a previous study showed that the activity of 6-phosphogluconate dehydrogenase underwent a significant change in poplar xylem between winter and summer [66], the high expression level of 5A-m0104 in ecodormancy may contribute to the physiological change during the transition from endodormancy to active growth. These findings raised the possibility of a regulatory role for miRNAs in metabolism and cell signaling, as well as epigenetic changes between dormancy and active growth.

Conclusions
In summary, a genome-wide sRNA profile of the cambial meristem was performed to present the miRNAs involved in the cambial dormancy-active growth transition. As a result, 123 known miRNAs, comprising 26 miRNAs families with obvious expression changes, were obtained, which included developmental-, phytohormone-, stress-and physiological-related miRNAs. In addition, 275 novel miR-NAs expressed in the cambium were identified, and 34 of them displayed unique expression patterns during the dormancy-active growth transition process. The relative expression levels of four novel and two known miRNAs were also confirmed by qRT-PCR. We predicted the target genes of these novel miRNAs and experimentally validated some of them using 5′RACE. This revealed not only important known miRNAs, which may contribute to the regulation of the dormancy-active growth transition, but also novel miRNAs and their possible target genes, which would provide new insights into the regulatory mechanisms of the process in trees.

Plant material
Poplar (Populus alba × Populus glandulosa cv. "84 k") plantlets were cultured on 1/2 Murashige and Skoog (MS) media with 20 g/L sucrose at 25°C under a LD photoperiod for a month. The induction of dormancy and resumption of growth in poplar was constructed as previously described with some modifications [36]: The plantlets were transferred into a greenhouse at 25°C under a LD photoperiod for at least 6 months until approximately 1 meter in height. Then, the healthy trees were transferred to a growth chamber for SD treatment, the apical buds were observed and the cambial cell layers of the 10 th internodes were counted in anatomical sections taken weekly during the process to evaluate the dormant state. After 8 weeks of the SD treatment, the dormant trees were treated with chilling (4°C) in an illumination incubator under the same conditions for 5 weeks. The chilling-treated trees were transferred to the LD condition at 25°C and the sprouting of dormant buds was inspected every week. Then, the trees released from dormancy grew in the LD condition at 25°C for another 3 weeks to achieve the active growth stage. At least 12 trees were used in every step of the process.
For the anatomical sections, the stems of the 10 th internodes were fixed in formalin/acetic acid/alcohol, dehydrated in a gradient of ethanol solutions, and finally embedded into Spurr's resin according to the manufacturer's description. Sectioning was performed using a Leica microtome, stained with 0.1% (w/v) toluidine blue O (Sigma, St. Louis, MO, USA), and observed under a Zeiss Axioskop 2 Plus microscope equipped with a computer-assisted digital camera.

sRNA sequencing and bioinformatic analysis
The cambial zones of 12 plants were sampled in SD8, C5 and LD3. The Illumina sequencing of sRNAs was performed following a previously published protocol [67]. The plant materials from the cambial zones were carefully scraped and ground in liquid nitrogen, and total RNAs were immediately extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA), and then separated on a 15% denaturing polyacrylamide gel. The 18-26 nt long sRNAs were excised and recovered. 5′ and 3′adapters were ligated to the isolated sRNAs, which were sequentially reverse transcribed and amplified by PCR. The purified PCR products were sequenced using a Solexa 1G Genetic Analyzer (Illumina, USA) at the Beijing Genomics Institute (BGI), Shenzhen, China.
The sequenced raw data were transferred into clean reads after removing the contaminants, low-quantity reads and adapters, which were used for the size distribution. SOAP software [68] was employed to map the clean reads to the poplar genome (http://www.phytozome. net/poplar). GenBank (http://www.ncbi.nlm.nih.gov/ genbank/), NCBI (http://www.ncbi.nlm.nih.gov) and Rfm (http://www.sanger.ac.uk/Software/Rfm/) databases were used to annotate the rRNA, scRNA, snoRNA, snRNA and tRNA in the sRNA library [69]. The clean reads mapped to the exons and introns of mRNA in the poplar genome were annotated as degraded sequences. Known miRNAs were identified by alignment to sequence in miRBase 18.0 (http://www.mirbase.org/) with no mismatches [70].
The remaining unannotated reads were used to predict novel miRNAs through the following methods: the MIR-EAP (http://sourceforge.net/projects/mireap/) software was employed to predict potentially novel miRNAs. First, the secondary structures of the sRNA precursors were predicted by the RNAfold web server [71] (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi) with default parameters. Additionally, we used the essential criteria for screening the miRNA candidates according to Meyer et al. [72].

Target prediction of novel miRNAs
The miRNA target candidates obtained through the MIR-EAP software were used to predict their target genes. The sequences of novel miRNAs were submitted to the psRNA-Target server (http://bioinfo3.noble.org/psRNATarget/), which contains plant miRNAs to screen target genes from Phytozome v. 9.1 (http://www.phytozome.net/poplar/) with criteria described previously [73,74]. The predicted target genes were annotated by searching the Phytozome v. 9.1 and NCBI databases.

miRNA expression analysis
To calculate the relative miRNA expression level and determine if there was a significant expression-level change, we used the log2-ratio and Scatter plot to compare the expression levels of miRNAs expressed in the three libraries based on previously established methods [75,76]. First, samples were normalized to 1 million, regardless of the total number of miRNAs in each sample. After normalization, if the expression level of a miRNA was 0, then it was revised to 0.01; if the expression level of a miRNA gene in all the libraries was less than 1, then this miRNA was removed because its expression level was too low. The normalized reads were used to calculate the fold change and p-value as follows: Fold change = log2 (the normalized treatment reads/ the normalized control reads), and where N1 is the total number of reads in the control sequencing library (SD8), N2 is the total number of reads in the treatment sequencing library, x is the number of reads for a miRNA in the control library and y is the number of reads for a miRNA in the treatment library. In this study, we used endodormancy (SD8) as the control. All calculations were performed on a BGI Bio-Cloud Computing platform (http://cloud.genomics.org.cn).
qRT-PCR of miRNA expression sRNAs were isolated from the cambial zone materials of 12 plants in LD3, SD8 and C5 using a miRNApure Mini Kit (CW Biotech, Beijing, China) following the manufacturer's instructions. Then, the sRNA was polyadenylated by poly (A) polymerase, and first-strand cDNA was obtained from polyadenylated sRNAs using the miRNA cDNA Kit (CW Biotech, Beijing, China) following the manufacturer's instructions. qRT-PCR was carried out as described: the SYBR Premix Ex Taq™ kit (TaKaRa Bio Inc., Japan) and an ABI 7500 Fast Real-time PCR machine (Applied Biosystems, Foster City, CA, USA) were used to complete the amplification, and the reaction procedure was set up according to the manufacturer's protocol. Three replicates were performed for each sample with 5.8S rRNA as an internal reference [46], and we used the 2 -ΔΔCT relative quantification method to calculate relative changes in gene expression [77]. All the primers are listed in Additional file 3: Table S3.

miRNA-mediated cleavage of mRNA
To identify cleavage sites in the target mRNAs, a modified RLM-RACE was performed using a GeneRacer Kit (Invitrogen, Carlsbad, CA, USA). All the steps followed the manufacturer's description, except that the calf intestinal phosphatase treatment was omitted to maintain the cleaved transcripts. All the primers are listed in Additional file 3: Table S3.

Availability of supporting data
The data sets supporting the results of this article are included within the article and its additional files.