Identification and characterization of microRNAs in tree peony during chilling induced dormancy release by high-throughput sequencing

Tree peony, one of the most valuable horticultural and medicinal plants in the world, has to go through winter to break dormancy. Growing studies from molecular aspects on dormancy release process have been reported, but inadequate study has been done on miRNA-guided regulation in tree peony. In this study, high-throughput sequencing was employed to identify and characterize miRNAs in three libraries (6 d, 18 d and 24 d chilling treatments). There were 7,122, 10,076 and 9,097 unique miRNA sequences belonging to 52, 87 and 68 miRNA families, respectively. A total of 32 conserved miRNAs and 17 putative novel miRNAs were identified during dormancy release. There were 771 unigenes as potential targets of 62 miRNA families. Total 112 known miRNAs were differentially expressed, of which 55 miRNAs were shared among three libraries and 28 miRNAs were only found in 18 d chilling duration library. The expression patterns of 15 conserved miRNAs were validated and classified into four types by RT-qPCR. Combining with our microarray data under same treatments, five miRNAs (miR156k, miR159a, miR167a, miR169a and miR172a) were inversely correlated to those of their target genes. Our results would provide new molecular basis about dormancy release in tree peony.


Results
Deep sequencing of Paoenia ostii sRNAs. To investigate small RNA expression profiles in Paoenia ostii during physiological dormancy stages based on the results of morphologic observations 11 , three small RNA libraries of flower buds were constructed after 6 d, 18 d and 24 d chilling enduring. For each library, small RNAs were collected, pooled together and sequenced. A total of 19,762,599 reads with lengths of 16 bp to 30 bp were obtained from the three libraries, and average 3.8 million (range: 2.69-5.05 million) clean small RNA reads were acquired from each library after removing adapters and low-quality reads. The average number of unique reads per library was 2.06 million ranging from 1.63 to 2.68 million (Table 1). Most obtained sRNA sequences were 21-24 nt in all of the three libraries, of which 24 nt long sRNAs were the most abundant, accounting for approximately 52.3% on average (Fig. 1).
Clean data were searched against Rfam databases to annotate sRNAs, and known miRNAs were identified by alignment to sequences in miRBase 20.0 with no mismatch. siRNAs, ribosomal RNAs (rRNAs), tRNAs, snRNAs and snoRNAs were annotated by BLASTn to NCBI Genbank database and Rfam database. In order to eliminate the possibility of degraded mRNAs in three libraries, we aligned them through intron/exon alignment with unigenes in tree peony cDNA libraries 11 . The remaining unannotated sRNAs were used to predict novel miRNAs  and potential miRNA seeds ( Table 2). It is noticeable that the miRNAs represented 19.62% of the total sRNA in 6 d chilling treatments, but only 14.01% and 14.78% in 18 d and 24 d chilling treatments, which may as a result of many genes associated with endo-and eco-dormancy release are activated during dormancy release. There are about more 9,000 unique miRNAs at the physiological stage of dormancy release and eco-dormancy than at the status of dormancy, which indicates that the miRNA populations in flower buds after dormancy release are more diversified, as well as biological processes are more complex.

Nucleotide Preference of 19-24 nucleotide small RNAs.
Previous studies have shown that most miRNA sequences start with uridine (U), whereas the majority of 24-nucleotide siRNAs have adenosine (A) as their 5′ first nucleotide in plants [26][27][28][29] . In our result, the same trends were observed among cloned tree peony small RNAs, about 67.1% of small RNAs sequences started with uridine, and all 24-nt siRNAs had start-nucleotide preference with adenosine (Fig. 2b). Besides, we found that about 45.5% of total small RNAs also had a clear preference for adenosine as the last nucleotide, of which all 23-nt had a clear preference for uridine and all 24-nt for adenosine as the last nucleotide (Fig. 2b). In order to investigate whether the last nucleotide preference from 19 nt to 24 nt small RNAs also exist in other model plants like Arabidopsis, we downloaded 427 Arabidopsis small RNA deep sequencing datasets from miRBase database (http://www.mirbase.org/ftp.shtml) and analyzed their nucleotide compositions (Supplementary Figure S1). In all 427 Arabidopsis datasets, strong last nucleotide preference for uridine was observed in 23 and 24 nt small RNAs, indicating that difference of nucleotide preference might exist among species.
Identification of conserved and novel miRNAs in tree peony. To identify conserved miRNAs in tree peony, but its genome is not available, known plant miRNAs (miRNA precursors and mature miRNAs) registered in miRBase 21.0 were used as reference for miRNA mapping. Clean data that aligned to known miRNAs allowing two mismatches and had no less than 5 reads per million (RPM) in at least one library were regarded as conserved miRNAs. In three libraries, total 112 known miRNAs belonging to 99 miRNA families were identified in the three libraries, of which there were 7,122, 10,076 and 9,097 unique miRNA sequences belonging to 52, 87 and 68 miRNA families, respectively (Additional file 1, Table 3). Of which, there was 32 conserved miRNAs (Table 4). In our data, 15 miRNAs sequences were found anchored in the 5p-arm and 17 miRNAs anchored in the 3p-arm (Addition file 1, Table 4). Unexpectedly, one less-conserved miRNA (miR5072) was obtained,  which was previously found only in monocots 30 . Furthermore, there were 55 miRNAs belonging to 46 miRNA families shared in the three libraries (Fig. 3), and the most abundant miRNA identified by sequence homology was PsmiR159 with more than 500,000 actual sequencing reads, accounting for approximately 69% of the total conserved miRNA reads, following by PsmiR5266 with more than 100,000 actual sequencing reads, PsmiR166, PsmiR319, PsmiR1509 with more than 10,000 actual sequencing reads, and miR398 showed the minimum amount (Table 4). At the same time, the frequencies of miRNAs read varied from 8 (PsmiR2111a-5p) to 709,087 (PsmiR159a), which indicate that miRNAs displayed drastically different expression level in tree peony during dormancy release (Fig. 3). After normalization, more than half of the conserved miRNAs were less than 100 times. In addition, the relative abundance of certain members within same miRNA family varied widely (  (Table 4). These results suggest that members showed different expression trends within same miRNA family, probably because their expressions are development-stage specific or either induced or suppressed during dormancy release in tree peony. Based on the miRNA annotation criteria 31 , novel miRNAs could be predicted by mapping the remaining non-annotated sRNAs to Populus genome. In our data, seventeen novel miRNAs were obtained and named as PsmiR1 to PsmiR17 (Table 5). All precursors of novel miRNAs had regular stem-loop structures and the predicted hairpin structures. To investigate the conservation of these 17 novel miRNAs in other plant species including Malus domestica, Physcomitrella patens and Populus trichocarpa, they were used to perform BLAST searches against miRBase databases. PsmiR5, PsmiR7 and PsmiR16 matched genomes of other plant species (Table 5). Reverse transcript PCR (RT-PCR) was performed to validate the expression of some new predicted miRNAs in flower buds after chilling treatments. The primer sequences were listed in Supplementary Table S1. We found five of the 17 predicted miRNAs including PsmiR9, PsmiR3, PsmiR1, PsmiR4 and PsmiR13 expressed in tree peony flower buds (Fig. 4). PsU6 was amplified as a positive control. We found that these novel miRNAs could be detected in flower buds after 18 d chilling treatments.
We found 469 mature miRNAs could be aligned to other species' precursors (precursors data from miRbase r21), mostly in Glycine max (36.7%), Oryza sativa (30.7%), and Arabidopsis thaliana (30.06%). The 17 precursors of novel miRNAs derived from tree peony transcriptome was listed in Table 5, only 1 precursors coded both 5p and 3p side mature miRNAs, others only possessed 1-side mature sequence.
Prediction of miRNA targets in tree peony. Previous study found that plant miRNA target sites mainly situate at opening reading frames (ORFs) 32 . To understand possible biological functions of the identified miRNAs in tree peony, we predicted the miRNA targets using the mRNA transcriptome of Paoenia ostii flower buds as a reference sequence since the genome of Paoenia ostii is not publicly available 10 . A total of 771 unigenes were predicted as potential targets of 62 known miRNA families (Additional file 2), and the majority of target proteins and corresponding annotations were shown in Table 6. Most miRNAs had more than one predicted target proteins, and some of the miRNAs have more than 10. Based on GO annotation, more than half of the predicted target genes were involved in biological process (metabolic process, regulation of transcription, signal transduction, transport and regulation of act polymerization) and molecular function (binding and methyltransferase activity) (Fig. 5). However, there were many conserved miRNAs target genes that had no functional annotation. Novel miRNAs targets were also predicted, but only two of them have been found target relationship with two unigenes ( Table 5).
5′-RNA ligase-mediated (RLM)-RACEs were carried out to confirm the miRNA-guided cleavage and cleavage sites of predicted target transcripts. Squamosa-promoter-binding protein-like (SPL) family genes and APETALA2 (AP2) had been reported that they were the predicted targets of miR156 and miR172, respectively 15,33 . Our results showed that PsAP2 (JI446524) could be cleaved at the site between bases 12 (T) and 13 (C) within the complementary region to PsmiR172a (Fig. 6). PsSPL (JI446831) could be cleaved by PsmiR156 at the site between 10 (C) and 11 (T), which was also identified as a miRNA cleavage site in rose 33 .
Expression profile of Paoenia ostii miRNAs during dormancy release. To identify miRNAs that were responsive during dormancy release, we compared miRNA expression level among three libraries. All the conserved candidate miRNAs with no less than 10 reads in each library were analyzed. Differentially expressed miRNAs that exhibited more than a 2-fold change were selected between each two treatments. There were 112 known miRNAs differentially expressed among three libraries (Additional file 3), and of which 55 miRNAs were shared among three libraries (See Additional file 3-shared miRNAs in three libraries). There existed 28 miRNAs (including miR8126-3p, miR6479, miR2949b, miR5057, miR6144, miR7743-3p, miR6483, miR5083, etc.) only in 18 d chilling duration library, which might play important roles in regulation dormancy release (Additional   (Table 7). There were 11 down-regulated and 43 up-regulated miRNAs from 6 d to 18 d chilling treatment. Among them, PsmiR5072 showed the highest degree of induction (7.7-fold), PsmiR390a-3p, PsmiR5059, PsmiR170-5p, PsmiR166a, PsmiR2109 and PsmiR5077 were also clearly up-regulated (>5-fold). In addition, PsmiR5519 was specifically and significantly induced during dormancy release. To further elucidate the potential regulatory roles of miRNAs in the transition from dormancy to eco-dormancy, we made a comparative analysis of miRNA expression between 18 d chilling and 24 d chilling. There were 31 down-regulated and 24 up-regulated miRNAs, among them PsmiR1509b was the most significantly down-regulated (>35-fold). Similar to the report in poplar 34 , eight known development-related miRNA families were also detected differentially expressed during dormancy release in our data, including PsmiR164, PsmiR396, PsmiR168, PsmiR319, PsmiR171, PsmiR166, PsmiR156 and PsmiR172. These miRNAs mainly acted on cell proliferation (miR164, miR396 and miR319) 28,35 , vascular development (miR166) 36 and miRNA biogenesis (miR168) 37 . PsmiR164 and PsmiR168a were continuously induced from dormancy to eco-dormancy stage, which were different from that in poplar 34 . PsmiR166 was up-regulated during dormancy release and repressed in eco-dormancy stage, similar results were detected in poplar during chilling induced dormancy-release 34 . Members of the PsmiR171 and PsmiR166 families showed the same expression patterns, but distinct differences of expression levels were also observed within other families during the same process. For example, PsmiR159a was up-regulated from 18 d to 24 d chilling treatments, while PsmiR159b-3p was down-regulated (>2-fold). PsmiR168a-3p was repressed during dormancy release (from 6 d chilling to 18 d chilling treatments), but PsmiR168a was continuously induced, indicating that members from same miRNA family might play different roles during this process. In addition, there were 28 miRNAs detected only in 18 d chilling treatment library and 7 miRNAs (PsmiR5227, PsmiR5665, PsmiR1886.1, PsmiR774b-5p, PsmiR4357, PsmiR1217-5p and PsmiR5224b) only in 24 d chilling treatment library (Additional file 4), which might mainly function in the transition from dormancy to dormancy release stage and eco-dormancy stage, respectively.
The sequencing results showed that the abundance of novel miRNAs was relatively less than that of conserved miRNAs (Table 5). Among these 17 predicted novel miRNAs, PsmiR1 was dramatically up-regulated during dormancy release. PsmiR3 and PsmiR4 were sharply reduced from 6 d to 18 d chilling treatments, PsmiR9, PsmiR10 and PsmiR13 showed opposite pattern with that of PsmiR3 and PsmiR4. PsmiR7, PsmiR11, PsmiR12, PsmiR15 and PsmiR17 were detected only in flower buds after 18 d chilling treatments.
Expression Validation of tree peony miRNAs. To confirm the expression patterns of miRNAs, as well as detect their responses to chilling treatments at three physiological stages, the expression of 15 conserved miR-NAs, whose sequencing counts altered significantly after treatment, were analyzed by RT-qPCR. It is showed that the expression levels of miRNAs were a constant change process with the time of treatment. We classified them into four types (Fig. 7). Type a: slowly increased (PsmiR3630, PsmiR390b-5p, PsmiR159a and PsmiR164a); type b: suddenly increased (PsmiR168a and PsmiR5072); type c: first increased and then decreased (PsmiR159b-3p, PsmiR160a, PsmiR166a, PsmiR167a, PsmiR169a, PsmiR319-3p and PsmiR172a); type d: first decreased and then increased (PsmiR156k and PsmiR157a). These results suggest that miRNAs belonging to type c were early stage   response miRNAs, those belonging to type a and type b might accelerate endo-dormancy release. Notably, the expression trends of two PsmiR159 family members, PsmiR159a and PsmiR159b-3p were different, indicating that the expression of miRNAs is a multiform process with the altered time of chilling treatments. Combining with our microarray data under same treatments 11 , five target genes of five miRNAs (miR156k, miR159a, miR167a, miR169a and miR172a) showed inverse expression patterns (Fig. 7e).

Discussion
Tree peony is an important horticultural crop worldwide of great ornamental, medical and economic value. Native to China, tree peony is regarded as "King of flower" and have deep botanical history in Chinese culture. It is crucial to understand the molecular mechanism of dormancy, which is a main obstacle for tree peony forcing culture. Based on morphological changes of Paoenia ostii 'Feng Dan' and global mRNA expression profiling, the physiological status of flower buds receiving less than 18 d chilling treatment are regarded as endo-dormancy, and that receiving more than 18 d are defined as eco-dormancy 11 . miRNAs are paid more and more attention as key regulators of gene activity in animals and plants 26,38,39 . In this study, we adapt high-throughput sequencing technology to identify sRNAs from Paeonia ostii and analyze their response to dormancy release. This work will provide useful information to deepen our understanding of the miRNA regulatory mechanisms during dormancy release.
There are 243, 511 and 207 miRNAs annotated in Arabidopsis, rice and soybean according to the miRBase database, respectively [40][41][42] . In this study, we first completed construction of sRNA libraries (   than others. Similar results were observed from most plants, such as Arabidopsis, rice, tomato 43 , cucumber 44 , apple 45 , peach 46 and rose 33 . However, 21 nt-long sRNAs were the second enriched in this study, which was different with previous reports in Arabidopsis (Arabidopsis thaliana) and rice (Oryza sativa) 30,47,48 . But in poplar, the 21 nt sRNAs are the most abundant 49 . Most of 21 nt sRNAs in our data started with uridine, which was consistent with the observation that ARGONAUTE1 (AGO1) usually harbors miRNAs with a 5′ terminal uridine. 24 nt sRNAs had start-nucleotide preference of adenosine, which was also reported in previous work 26,27,29 . Tree peony endo-dormancy transcriptome database was employed to predict putative targets of tree peony miRNAs. The well-known conserved miRNAs including miR156, miR159, and miR164 have been identified.  Table 6. Majority of the predicted target genes and corresponding annotation of known miRNAs in tree peony. However, nearly half of known miRNAs and three novel miRNAs were not predicted homologous to any proteins in the Genbank nr database, which might because of the incomplete tree peony genome and limited number of transcript data in public database. Hundreds of miRNAs have been surveyed since high-throughput sequencing technology is widely used, but little has been done on identifying and analyzing miRNAs in tree peony and their response during dormancy release. Total 112 known miRNAs belonging to 99 families were identified in the three libraries. Two miRNAs families, PsmiR159 and PsmiR166, were relatively abundant. PsmiR166a increased continuously until dormancy release (6 d -18 d) and had a very low expression level at eco-dormant stage (24 d), which was consistent with recent studies in poplar 34,50 . For the expression level of 17 novel miRNAs, PsmiR1 was continuously up-regulated from dormancy to eco-dormancy stage, PsmiR9, PsmiR10 and PsmiR13 (up-regulation) and PsmiR3 and PsmiR4 (down-regulation) had converse expression pattern at the early stage of dormancy release.

Cold-responsive and auxin-related miRNAs. A continuous chilling accumulation is an effective natural
or artificial way to release dormancy in tree peony. The endo-dormant bud can respond to chilling treatment, which stimulates growth-promoting respond signals including auxin or appropriate outside conditions (such as exogenous GA). In our study, there were 112 conserved miRNAs differentially expressed between 6 d and 18 d chilling library and between 18 d and 24 d chilling library. Among them, PsmiR160a was highly expressed in endo-dormancy release (18 d chilling) and quickly decreased in eco-dormancy (24 d chilling). MiR160 targeted auxin response factor 10/16/17 (ARF10/ARF16/ARF17) and negatively regulated auxin signaling 51,52 . Ding et al.   Table 7. Differentiated expressions of shared miRNAs from tree peony flower bud after different chilling treatments. Note: + and − indicate the induction and repression of miRNA, respectively.
found that miR160 was highly expressed in eco-dormancy (five weeks cold treatment) 34 . In our case, the high expression of PsmiR160 in endo-dormancy release might because of the difference of dormancy mechanism exist in tree peony and poplar. On the other hand, the expression of miR160 in endo-dormancy release may strengthen the auxin signal by inhibition of its target genes, and this hypothesis is also consistent with earlier report that exogenous GA could effectively promote endo-dormancy release 10 . Ding et al. reported other two auxin-related miRNAs, miR390 and miR167, increased during active growth 34 . In our study, the expression of PsmiR390b was steadily increased during the transition from endo-dormancy to eco-dormancy, and PsmiR167a was significantly induced during dormancy release, which suggest that auxin signal pathway participated in the process of dormancy release. We also found that PsmiR168a was continuously up-regulated from endo-dormancy stage to eco-dormancy stage, the same trend was found in poplar 34 . Similarly, miR168, member of the Csn-miR168 family, was found to be a cold-responsive miRNA, which was induced in two tea cultivars after 12 h of cold treatment 22 .
MiR168 regulates its target ARGONAUTE1 (AGO1) to participate in miRNA biogenesis 37 . The high expression level of PsmiR168a would lead to the repression of AGO1, which would cause a reduction in the miRNA expression levels. The up-regulation of PsmiR168a suggested that cold-responsive miRNA participated in the release of endo-dormancy, their inductions were also consistent with its function of miRNA biogenesis 37 .
MiRNA targets. Since the genome of Paoenia ostii is not publicly available, the mRNA transcriptome of Paoenia ostii flower buds 11 were employed as a reference to predict the putative miRNA targets. Based on GO annotation, more than half of the predicted targets in tree peony were involved in binding, catalytic activity, metabolic process and cellular process. For example, miR5141 targets gene encoding ATP synthase, which have been reported to be involved in ATP synthesis and ATP utilization during seed dormancy breaking 53 . In pear, specific enrichment of unigenes was observed for 15 pathways involved in metabolic processes including oxidative phosphorylation 4 . Several other target transcripts, which encode alpha N-terminal protein methyltransferase, Endoglucanase, GTPase-activating protein and F-box domain associated with various biological processes and cellular activities were also detected. For instance, PsmiR395 targets genes encoding enzyme beta galactosidase, which recently have been reported to be involved in cell wall modification during the transition from dormancy to eco-dormancy in onion bulbs 54 . PsmiR171 targets genes coding endoglucanase, which had been shown to be antifreeze proteins during seed germination in sunflower 55 . F-box proteins, the target of PsmiR169, have been identified previously as a key regulator of karrikin signaling and seed dormancy in Arabidopsis 56 . Novel miRNAs target genes were also predicted, but only two of them have not been found target relationship.
MiR169 and miR166 regulated cellular process and biological process by acting on their target genes. MiR166 function mainly in vascular development 36 , and the down-regulation of PsmiR166a at eco-dormant stage might help to increase the expression level of its target gene. In addition, Potkar et al. found that ptrmiR169a and its target gene PtrHAP2-5 showed inverse expression patterns during the dormancy period, which suggests that miR169 mediate attenuation of the target HAP2-5 transcript at this process 24 . Jeyaraj et al. found that CsmiR169 targeted COBRA-like protein encoding gene and regulated cellulose synthase, which suggests that miR169 have possible role in cell cycle and other biological function during the bud development 23 . In our study, PsmiR169a was highly expressed at the early stage of dormancy release and steadily down regulated at eco-dormant stage, and similar results were obtained during vegetative bud dormancy period of aspen 24 .
MiR156 and miR172 regulate and control the juvenile-to adult vegetative transition by targeting transcription factors SQUAMOSA promoter-binding protein-like (SPL) and APETALA2 (AP2) genes in both annual herbs 17,57 and woody perennial plants 58 , showing converse expression patterns and regulatory relationships 57 . It is noteworthy that in our study we found the expression levels of PsmiR156k and PsmiR172a during the transition from endo-dormancy to eco-dormancy also had the converse expression patterns. Similarly, Ding et al. also found miR156 and miR172 showed completely converse expression patterns during the dormancy-active growth transition 34 . Our results showed that target genes (JI446524 and JI446831, putative AP2 and SPL genes) had cleavage sits of PsmiR172a and PsmiR156a, respectively, which suggested that miR156 and miR172 might play an important role during dormancy transition, which need to be further confirmed by experiments.

Materials and Methods
Plant materials. Four-year-old tree peony plants (Paoenia ostii 'Feng Dan') were potted and moved to refrigeration house with temperature 0-4 °C from 5 November to 30 December, 2014 in Qingdao, Shandong, China. The morphologic observation indicated flower buds receiving less than 18 d chilling treatment are in physiological status of endo-dormancy, while those receiving more than 18 d chilling treatment were in eco-dormancy physiological status 11 . Therefore, in order to decrease individuality, more than 5 plants were mixed buds-three buds for each individual were collected after 6 d, 18 d and 24 d chilling requirement fulfilling. Three replicates samples were harvested and immediately frozen in liquid nitrogen and stored at −80 °C until further use.
Small RNA library construction and sequencing. Total RNA from tree peony flower buds after chilling treatments (6 d, 18 d and 24 d) was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instruction and separated on 15% denaturing polyacrylamide gels. The 16-30 nt sRNAs were excised and recovered. The adapters (5′ and 3′) were ligated to the isolated sRNAs, which were sequentially reverse transcribed and amplified by PCR. The purified PCR products were sequenced using an Illumina Genome Analyzer (Illumina, USA) at Beijing Biomarker Technologies, Beijing, China.
Analysis of sequencing data. Raw sequence reads were produced by Illumina Genome Analyzer at Biomarker-Beijing, China and processed into clean full length reads by the Biomarker small RNA pipeline. During this procedure, all low-quality reads including 3′ adapter reads and 5′ adapter contaminants were removed. The remaining high-quality sequences were trimmed of their adapter sequences. Sequences larger than 30 nt and smaller than 16 nt were discarded. All high-quality sequences were considered as significant and further analyzed.
All matched sRNA sequences were categorized into classes including miRNAs, siRNAs, ribosomal RNAs (rRNAs), tRNAs, snRNAs, snoRNAs, repeat-associated sRNAs and degrade tags of extrons of introns, etc. Then, SCientifiC RepoRts | (2018) 8:4537 | DOI:10.1038/s41598-018-22415-5 the clean sequences were annotated by performing BLASTn searches against the Rfam (http://www.Sanger.ac.uk/ Software/Rfam) and NCBI (http://blast.ncbi.nlm.nih.gov/Blast.cgi) databases, the detailed processes were following: the clean data were aligned with known miRNAs (miRNA precursors and mature miRNAs) registered in miRBase 21.0 (http://microran.sanger.ac.uk/sequence/index.html/) because of the difference among species, this process allowed two mismatches and free gaps to get temporary miRNA sequences and count of miRNA families; The highest-expression miRNA for each temporary mature miRNA family were selected to form a temporary miRNA database of Paoenia ostii. Finally, alignment of clean data to temporary miRNA database to identify conserved miRNAs in Paoenia ostii, only those perfect matching (≤two mismatches) were considered as conserved miRNAs. Potential novel miRNAs were identified using criteria that were previously developed for plant miRNA prediction 59 . The unique fold back structures of miRNA precursors can be utilized to predict novel miRNAs using MIREAP program (http://sourceforge.net/projects/mireap/). Potential targets for both known and novel miRNAs were identified on TAPIR and Target Finder based on Paeonia ostii transcriptome sequencing data 10 according to the search algorithm that only three or fewer mismatches and no gap are allowed to be present in the complementarily between miRNAs and their corresponding targets 32 . The biological function category of the predicted targets was annotated against the Universal Protein Resource (http://www.uniprot.org).
Differential expression analysis of miRNA and Reverse Transcriptase quantitative PCR (RT-qPCR) and 5′ RLM-RACE. Differential expression analysis of miRNAs was performed based on sequence reads generated from three libraries after different chilling treatments according to the method described by Ren 49 . In detail, the expression of miRNAs was normalized to obtain the number of miRNAs per million reads [normalized expression = (number of miRNA reads/total number of clean reads) × 1,000,000]. Normalized miRNA reads with values less than one in three libraries were excluded. The remaining miRNAs were used to calculate differences in expression by fold change (normalized miRNA reads in 18 d or 24 d chilling treatment/normalized miRNA reads in 6 d chilling treatment) and significant P-values 60,61 .
To validate miRNA expression, sRNAs were isolated from flower buds after different chilling treatments using an RNAiso for small RNA (TaKaRa, Dalian, China) following the manufacturer's instructions. Then, the sRNA was polyadenylated by poly (A) polymerase, and first-strand cDNA was obtained using SYBR ® Primescript miRNA RT-PCR Kit (TaKaRa, Dalian, China). Briefly, the polyA was added to the 3′ of total RNA, then the RNA was reverse-transcribed with an oligo-dT adaptor. Quantitative PCR was performed in a total volume of 25 μL, containing 2 μL cDNA, 0.4 μM PCR forward primer (1 μL), 0.4 μM Uni-miR RT-qPCR primer (1 μL), 12.5 μL of 2 × SYBR premix Ex Taq ΙI, and 8.5 μL dd H 2 O. The reactions were completed using Roche Light Cycler 480 (Roche, Mannheim, Germany) with the following program: 95 °C for 10s and 40 cycles of 95 °C for 5s, 55 °C for 30s and 72 °C for 15s. The reactions were run in triplicate and the 2 −ΔΔCt relative quantification method was used to calculate the relative changes in gene expression 62 . Small nuclear RNA U6 was used as endogenous reference, primers used in this study were listed in Supplementary Table S1.
To conform whether the predicted targets were cut by miRNAs and cleavage sites, the 5′ RLM-RACE were carried out using the FirstChoice RLM-RACE Kit (Ambion). Specifically, one microgram total RNA was firstly ligated to 5′ RACE oligo adaptor without calf intestine alkaline phosphatase and tobacco acid pyrophosphatase treatments. Then, the ligated RNA was used to synthesize the cDNA. The primers of miR172a target gene (JI446524) (5′-TCGGAGAAATGCTTTGTCCATGGCCAT-3′) and miR156a target gene (JI446831) (5′-TTGCGAGGTTCTGGGTTTGGAG-3′) for 5′ RLM-RACE were designed by Primer premier 5.0 software (Supplementary Table S1). PCR was carried out according to the manufacturer instructions, and the PCR products were purified by 1.0% agarose gel electrophoresis and cloned into the pMD18-T vector (Takara, Dalian, China) for sequencing.