Genome-Wide Identification of miRNAs Responsive to Drought in Peach (Prunus persica) by High-Throughput Deep Sequencing

Peach (Prunus persica L.) is one of the most important worldwide fresh fruits. Since fruit growth largely depends on adequate water supply, drought stress is considered as the most important abiotic stress limiting fleshy fruit production and quality in peach. Plant responses to drought stress are regulated both at transcriptional and post-transcriptional level. As post-transcriptional gene regulators, miRNAs (miRNAs) are small (19–25 nucleotides in length), endogenous, non-coding RNAs. Recent studies indicate that miRNAs are involved in plant responses to drought. Therefore, Illumina deep sequencing technology was used for genome-wide identification of miRNAs and their expression profile in response to drought in peach. In this study, four sRNA libraries were constructed from leaf control (LC), leaf stress (LS), root control (RC) and root stress (RS) samples. We identified a total of 531, 471, 535 and 487 known mature miRNAs in LC, LS, RC and RS libraries, respectively. The expression level of 262 (104 up-regulated, 158 down-regulated) of the 453 miRNAs changed significantly in leaf tissue, whereas 368 (221 up-regulated, 147 down-regulated) of the 465 miRNAs had expression levels that changed significantly in root tissue upon drought stress. Additionally, a total of 197, 221, 238 and 265 novel miRNA precursor candidates were identified from LC, LS, RC and RS libraries, respectively. Target transcripts (137 for LC, 133 for LS, 148 for RC and 153 for RS) generated significant Gene Ontology (GO) terms related to DNA binding and catalytic activites. Genome-wide miRNA expression analysis of peach by deep sequencing approach helped to expand our understanding of miRNA function in response to drought stress in peach and Rosaceae. A set of differentially expressed miRNAs could pave the way for developing new strategies to alleviate the adverse effects of drought stress on plant growth and development.


Introduction
Peach (Prunus persica L.) is considered to be one of the most widely grown and economically important stone fruit species in the Rosaceae, comprising more than 3,000 species in approximately 110 genera distributed worldwide [1]. In 2010, it was estaminated that world annual production of peaches and nectarines exceeded 19 million metric tons according to FAO statistics (FAOSTAT, http://faostat.fao.org). In addition to its ecological and high economic importance, peach is also emerging as a model tree species for both comparative genomic studies, evolutionary studies and plant development research owing to its small genome size of 300 Mb (just about twice comparing with Arabidopsis thaliana) and the relatively short reproductive time [2,3]. Peach has a haploid chromosome number of 8 [4] and the first draft of peach genome (peach v1.0, obtained from ''Lovell'' haploid) was unraveled by the International Peach Genome Initiative (IPGI), which are available from the Genome Database for Rosaceae (http://www.rosaceae. org/peach/genome). The genus Prunus, which includes peach, nectarine, apricot, weet and sour cherry, have stone fruits with fleshy mesocarp, but the growth and development of these fruits, especially large-fruited species like peach, are seriously affected by drought. For peach, drought stress is one of the major abiotic stresses limiting fruit production and quality during the 4-6 week period before harvesting [5,6].
As a major abiotic factor, drought can be described as basically the water deficiency or insufficient access to water and it has adverse effects on the growth of plants and crop production. However, plants growing in drought stress respond to dehydration and have to develop a variety of mechanisms at morphological and molecular level in order to remain alive. Deciphering the physiological processes and molecular genetic mechanisms involved in drought resistance has certainly made a significant progress in understanding of complex biological response of plants at the molecular and organism levels against the drought. The expression profile of protein-coding genes could highly fluctuate in response to drought at the transcriptional and post-transcriptional levels of miRNA [7,8,9]. The aberrant expresison of genes may be regulated by an newly discovered small RNAs, termed micro-RNAs (miRNAs).
Considerable amount of plant miRNAs have been identified by computational and/or experimental methods and these miRNAs have been deposited in the latest release of miRBase v18 (release 18.0 November 2011, only experimentally validated [23]) and PMRD (Plant microRNA database, both experimental and computational [24]) since the first miRNAs were discovered in plants in 2002 [25,26]. At present, there are 4053 hairpin entries pertaining to 52 plant species in the miRBase v18.0. The identification of any miRNAs has great importance for subsequent research such as miRNA function, nature, target prediction and biogenesis. In recent years, the innovative strategies and practical methodologies have been developed for determining miRNA expression. The main approaches of experimental methodologies can be summarized as follows: direct-cloning [27][28][29], stem-loop qRT-PCR [10,[30][31][32], next-generation sequencing technology [33][34][35][36][37] and hybridization-based detection, such as northern blotting [38], in situ detection [39,40] and microarray [31,41,42]. Although each of these methods has their own particular advantages and disadvantages, the next-generation sequencing technologies play an increasingly prominent role in discovering novel miRNAs [36,43] and measuring quantitively expression levels of low-abundant [44] and species or tissue specific miRNAs [35,45] comparing to other genome-wide transcriptome analysis methods, such as miRNA-microarray [46,47]. In addition to experimental approach, computational approach has also been a preferred method because of its low cost, high efficiency and speed, prior to experimental validation. The efficacy and power of computational approach come from its major characteristic features: (i) high evolutionary conservation from mossess to eudicots in plants (comparative genomics) [48,49], (ii) hairpinshaped stem-loop secondary structure [50,51], (iii) high minimal folding free energy index [48].
Many recent studies have revealed that plant miRNAs have pivotal roles in plant response to abiotic stresses, including drought [9,52,53], salt [54,55,56], cold [33,57], oxidative stress [58][59][60] and UV-B radiation [61]. The miRNAs whose expression level is significantly altered in drought condition compared with normal conditions have been well-reviewed in recent works [7,8]. It has been indicated that a certain number of miRNAs involve in response to drought stress by altering the gene expression [52,62,63]. As for Prunus species, there is no comprehensive data about the expression profiles of drought responsive miRNAs. The aim of the present study is to determine the expression profile of drougt stress-responsive miRNAs in peach. Thus miRNA deep sequencing by Illumina HiSeq 2000 were applied not only for simultaneous evaluation of drought responsive miRNAs' expressions, but also for providing comprehensive information about P. persica miRNA transcriptome on genome-wide scale. Stem-loop real time qRT-PCR (ST-RT PCR) was also employed to futher validate the expression level of a set of miRNAs identified during deep sequencing. Additionally, the identification and characterization of P. persica miRNAs and their target genes were established by using computational methods combined with experimental validation.

Results
We used the Illumina Solexa sequencing platform to investigate the genome-wide identification and expression profiles of miRNAs in peach, particularly for the drought-responsive miRNAs. Four small RNA libraries were constructed by the use of total RNAs isolated from control leaf (LC), drought-stressed leaf (LS), control root (RC), and drought-stressed root (RS) tissues. Small-RNA sequencing yielded a total of 53,878,885 high-quality raw sequence. Total high-quality raw reads in each of LC, LS, RC and RS libraries are 15,499,314, 12,473,137, 12,703,130 and 13,203 (Table 1 and Table S1). In order to get a bigpicture view of sequence distribution of all sRNA reads, all clean reads were mapped against the peach genome database at Genbank (http://www.ncbi.nlm.nih.gov/genome/388), Rfam (http://rfam. sanger.ac.uk/) and miRBase v18.0 (http://www.mirbase.org/), and therefore, they are classified into seven annotation categories: noncoding RNAs (tRNA, rRNA, snRNA and snoRNA), miRNA, exonsense, exon-antisense, intron-sense, intron-antisense, and unknown sRNAs (Table 2, Figure S1). As shown in Table 2, the highest abundance of unique conserved and potential non-conserved miRNAs reads was found in root-drought stress library and leaf control library, whereas most of the total miRNA reads were found in leaf-stress library. The length distribution of unique sRNA reads revealed that the majority of reads from each library were 18-25 nt in length, of which the class of 24 nt was the most abundant group accounted for average ,50% of total reads for each library and it was followed by the group of 21 nt class ( Figure 1). Although these small RNAs unevenly distributed in four groups according to their length, small RNAs in control and drought-exposed group for leaf and root represent similar distribution within each of their own group ( Figure 1).

Identification of known miRNAs in peach
In order to identify known (both conserved and species-specfic) miRNAs from control and drought-exposed root and leaf tissues of peach, small RNA sequences generated from each library were independently aligned with currently known and experimentally validated mature miRNAs deposited in miRBase v18.0, including 4,014 viridiplantae miRNAs belonging to 52 plant species. After homology search, a total of 531, 471, 535 and 487 miRNAs were identified from LC, LS, RC and RS libraries, respectively (Table  S2). These miRNAs belong to 43 evolutionary conserved miRNA families ( Table 3), suggesting that miRNA-mediated biological process are also present in peach as found in other plant species. However, some miRNAs, such as miR416, miR437, miR441 and miR529, were not detected in both leaf and root samples, suggesting that these miRNAs may be tissue-specific expression. The expression levels varied from miRNAs to miRNAs from one copy to more than one million of copies based on the deep sequencing (Table S2). A majority of miRNAs were detected with more than 50 copies; such as a total of 272 miRNAs for LC (51,22%), 229 miRNAs for LS (48,61%), 225 miRNAs for RC (46,20%) and 269 miRNAs for RS (50,28%) were sequenced more than 50 times. As previously reported, evolutionary conserved miRNAs have generally high expression abundances when compared with non-conserved miRNAs. Among the conserved miRNAs, total reads of miR535, miR157, miR166, miR156 and miR408 accounted for vast majority of total miRNAs; LC (82,49%), LS (89,09%), RC (75,54%) and RS (55,02%). Of these, miR535 was the most abundant miRNA in both control and drought-exposed libraries ( Table S2).

Identification of novel miRNAs in peach
After obtaining known miRNAs in peach, the remaining sequences of four libraries, which are classified as ''unannotated'' (excluding known miRNAs and Rfam matching other non-coding RNAs), were taken into consideration to discover novel and potential peach-specific miRNA candidates. To accomplish this, these small RNA sequences were aligned with the P. persica genome to identify genomic regions potentially harbouring potential pre-miRNA sequences whose hairpin-like structures are widely used for distinguishing miRNAs from other small noncoding RNAs. The minimum of free energy (MFE) of the secondary structures was also considered to be another criteria for prediction of potential pre-miRNAs. After aligning these unannotated sequences to the genome, we obtained a total of 197, 221, 238 and 265 novel miRNA precursor candidates for LC, LS, RC and RS libraries, respectively (Table S4) and some of these novel miRNA candidates with characteristic features are listed in Table 4. All novel miRNA prediction were carried out according to the default parameters of MIREAP (MicroRNA Discovery By Deep Sequencing) software developed by BGI. In agreement with previously reported results, the uracil nucleotide is dominant in the first position of 59 end for majority of these newly determined putative novel miRNAs. The first nucleotide bias analysis showed that uracil was the most frequently used first nucleotide in miRNAs of P. persica; with 10,528 uracil nucleotides (47%) for LC, 12,834 uracil nucleotides (43%) for LS, 21,571 uracil nucleotides (63%) and 13,014 uracil nucleotides (34%) for RS library (Table  S3). Our sequence analysis for all libraries showed that the putative pre-miRNAs of each library greatly varied from 70 to 365 nucleoties in length. With the usage of software mFold, these pre-miRNA sequences were applied to predict the characteristic stemloop secondary structure of pre-miRNA and their locations were also determined in the genomic loci (Tables S4 and S5). Some of the stem-loop secondary structures of predictive pre-miRNAs of P. persica determined via mFold can be found in Figure 2. We also calculated the minimum folding free energies of putative peach miRNA precursors for each libraries; ranging from 218,8 to 2157,40 kcal/mol with an average of 253,10 kcal/mol for LC, from 218,3 to 2171,23 kcal/mol with an average of 255,28 kcal/mol for LS, from 218,32 to 2157,4 kcal/mol with an average of 250,99 kcal/mol for RC and from 218,11 to 2181,01 kcal/mol with an average of 250,34 kcal/mol for RS (Table S4). In contrast with the common or evolutionarily conserved miRNAs, the predicted novel miRNAs are often expressed at a very low level as reported before. One possible explanation for this result was that many plant miRNAs are evolutionarily conserved and approximately one hundred miRNA  Table 2. Classification of small RNA sequences from control and drought stress libraries.
Genome-wide expression patterns of drought-responsive miRNAs identified in peach 262 and 368 miRNAs were observed with more than 2 fold change response to drought treatment in peach leaf and root, respectively (Table S6; Figure 3). As reported in previous studies [9,42,53,[66][67][68][69][70][71][72], a series of miRNAs, including miR156, miR159, miR160, miR165/miR166, miR167, miR168, miR169, miR170/ miR171, miR390, miR393, miR395, miR396, miR397, miR398, and miR408 are considered to be drought-responsive miRNAs ( Figure 4). Our analysis also revealed that miR165 and miR167 were commonly down-regulated in both leaf and root, whereas miR156 was slightly up-regulated in root and leaf after stress treatment. The expression level of miR159, miR169, miR393, miR397, miR398 and miR393 were only decreased in root under drought stress while the miR395 were only down-regulated in leaf in response to drought ( Figure 4). The miR160 were solely upregulated in root, whereas there were no changes in the expression level of miR168, miR390 in leaf and root tissues (Table 5). Since our results largely agree with previous studies (well reviewed in [7,8]), it has been said that peach has its own specific miRNA expression profile under drought-stress.

Target prediction and function analysis
The putative miRNA targets in peach were predicted using BlastN search (v2.2.22) against EST and cDNA sequences in P. persica genome annotation database (http://www.plantgdb.org/ XGDB/phplib/download.php?GDB = Pe) based on the rules described in the section of methods. Based on this strategy, a total of 672, 1293, 1194 and 1719 putative miRNA targets were obatined for LC, LS, RC and RS libraries, respectively. Although the number of targets for each miRNA varied considerably among the libraries ranging from 65 to 1 for LC, from 457 to 1 for LS, from 433 to 1 for RC and from 563 to 1 for RS, most of miRNAs in each library have only one predicted target transcripts; LC (30%), LS (30%), RC (32%) and RS (37%) ( Table S5). The sequence alignment between miRNA and its predicted target genes are also found in Tables S4 and S5. For comprehensive annotation, all putative target transcripts in each library were analyzed by Gene Ontology (GO) terms with the aid of Blast2GO program with default parameters. Among the peach miRNA targets identified, a total of 571 target transcripts (137 for LC, 133 for LS, 148 for RC and 153 for RS) generated significant GO terms for further analysis. Then, transcripts representing genes with a known function were categorized by biological process, cellular component and molecular function according to the ontological definitions of the GO terms. The putative target transcripts of miRNAs in the biological process category were   16 terms for RC and 6 terms for RS), transporter activity (5 terms for LC, 2 terms for LS, 4 for RC and 6 terms for RS), structural molecule activity (2 terms for RC and 1 term for RS), enzyme regulatory activity (2 terms for LC, 1 term for LS and RS and 3 terms for RC) and transcription regulator activity (1 term for RC and 2 terms for RS) ( Figure 5A). As shown in Figure 5, most of the miRNA target genes were assigned to the binding category whose present sequences appear to be involved in nucleic acid binding, protein binding and ion binding. Because these sequences encode transcription factors, this is in accord with previously reported notion explained as a large proportion of miRNA targets encode transcription factors [7,8,14]. In the biological-processes category, total number of miRNA target sequences in each library fell into multiple classes, however, it is notable that some of these transcripts in the biological process category were related to stress response process, for example; response to stress (13 terms), oxidation-reduction  process (9 terms) and response to abiotic stimulus (8 terms) for LC, oxidation-reduction process (12 terms), response to abiotic stimulus (7 terms), small molecule metabolic process (10 terms) and cellular catabolic process (6 terms) for LS, response to salt stress (7 terms), small molecule metabolic process (7 terms) and oxidation-reduction process (12 terms) for RC and oxidationreduction process (6 terms), cellular response to stress (7 terms), response to salt stress (5 terms) and catabolic process (5 terms) for RS ( Figure 5B). Based on KEGG biochemical pathway analysis, a total of 351, 530, 563 and 597 target transcripts involved in different cellular pathways were determined for LC (124 pathways), LS (216 pathways), RC (211 pathways) and RS (226 pathways) libraries, respectively (  3.60% for RC and 3.29% for RS). It is interesting to note that most of target genes were associated with plant hormone signal transduction because the environmental stress factors such as drought affect plant hormone balance and these biotic environmental stress factors are regulated by transcription factors, which are potential targets of most plant miRNAs.

qRT-PCR validation of P. persica miRNAs and target transcripts
We applied stem loop quantitative real-time RT-PCR (qRT-PCR) for further experimental verification of the presence of some conserved miRNAs and comparison of the expression pattern of these miRNAs with deep sequencing. Analysis of seven droughtresponsive miRNAs by qRT-PCR show that the expression level of miR156 and miR168 were high in leaves and roots under drought stress in comparison to control samples while the expression of miR164 and miR395 was down-regulated in root and leaf tissues of drought-stressed samples. The expression of miR169 was induced in leaf but inhibited in root tissues after drought treatment, whereas the expression level of miR171 was induced in root but inhibited in root tissues under drought. As for miR166, although its expression was down-regulated in root tissues in response to drought, the expression level of miR166 were not changed between control and drought stressed leaves of peach ( Figure 6). The relative expression profile of miR156, miR164, miR168, miR171 and miR395 using qRT-PCR had a good correlation with deep sequencing. However, there is a discrepancy between the results obtained from deep sequencing and qRT-PCR experiment. Deep sequencing results showed that the expression of miR169 was down-regulated in leaf tissue while qPCR experiment revealed that its expression was up-regulated in leaf tissue after treatment. Rather than experimental methods, duration of drought probably caused expression level differences.
qRT-PCR was also used for detection and quantification of predicted targets of six drougt-responsive miRNAs (miR156, miR164, miR166, miR169, miR171 and miR395). Our results revealed a negative correlation between the levels of miRNAs and those of their target messages ( Figures 6A and 6B). Thus, downregulation of miRNA could lead to increased expression of its target gene. For instance, the decreased expression of miR164 and miR395 promoted the expression of their targets genes in both root and leaf tissues. Conversely, drought-induced up-regulation of miR156 led to down-regulation of its target gene.

Discussion
Among various abiotic stresses, drought is considered to be one of the most detrimental factors to agriculture and adversely influence crop productivity and quality due to its high scale of impact and wide distribution [73]. As land plants are sessile organisms, they cannot escape from unfavourable environmental stress conditions surrounding them. Thus, land plants have to develop various mechanisms at the physiological and molecular    levels in order to cope with stress. Recently, miRNAs have turned out to be new players in plant tolerance to environmental stresses like drought, cold, heat and high salinity [74] and much effort has been devoted to understanding their role in the responses of drought stress in various plants including; Hordeum vulgare [9], Medicago truncatula [66], Populus euphratica [68], Vigna unguiculata [69] and Oryza sativa [75]. In the present study, deep sequencing technology was used for quantitative determination of genomewide miRNA expression patterns of P. persica in response to drought. A total of 535 known miRNAs were detected in peach, although 126, 256, 293, 329, 197, 157 and 126 known miRNAs were identified in M. truncatula [66], P. euphratica [68], V. unguiculata [69], A. hypogaea [76], G. max [77], P. aphrodite [78] and V. amurensis [79], respectively. Hence, it can be possible to deduce that the small RNA repertoire of peach is relatively richer than other plant species. By comparing the expression level of individual peach miRNAs in drought-stressed tissues to control, the expression level of 262 (104 up-regulated, 158 down-regulated) of the 453 miRNAs significantly changed in leaf tissue, whereas 368 (221 up-regulated, 147 down-regulated) of the 465 miRNAs had expression levels that significantly changed in root tissue (Table S6). Among these miRNAs, drought-responsive miRNAs (Table 5, Figure 2) were differentially expressed and showed fluctuations in their expression in both peach leaf and root. The expression of miR165 and miR167 was found to be significantly down-regulated in leaf and especially root under drought stress, whereas miR156 was slightly but not significantly up-regulated in drought-stressed tissues. It should be noted that expression of miR165/166 was induced in leaf but inhibited in root tissues of H. vulgare after dehydration stress [9] while miR166 was down-regulated in both leaf and root tissues of P. persica. Similarly, miR171 expression was up-regulated in leaves of barley while transcript level of miR171 were only increased in the root, but not in the leaf of peach in response to drought. Although miR167 was significantly down-regulated in leaf and root libraries of peach, its up-regulation was observed in the Arabidopsis thaliana [42] and P. euphratica [68]. Some miRNAs in different plant species display different expression patterns in response to drought; for example, miR168 is up-regulated and down-regulated in A.thaliana [42] and O.sativa [80], respectively while its expression level was not significantly affected by drought. As another drought-responsive miRNA, miR395 was downregulated and its expression was restricted to root tissue and was not detected in leaves of peaces ( Table 5). As the expression levels of miR159, miR396 and miR397 were down-regulated in peach after treatment, this finding is inconsistent with previous reports suggesting that expression of these miRNAs was up-regulated in both A. thaliana and P. euphratica [42,68]. Although the measured expression level of the miR166 did not change in leaf tissue of peach under drought stress, the qPCR results indicate that its expression level decreased in root tissue after treatment ( Figure 6). This result is consistent with the previous finding that the miR166 was downregulated in roots of barley after drought stress [71]. Because the function of miRNAs appears to be in gene regulation by targeting specific mRNAs for degradation or translation inhibition, identifying the potential target transcripts of miRNAs is crucially important for understanding miRNAmediated processes such as drought tolerance in plants. Therefore, target prediction analyses were particularly conducted for drought-responsive miRNAs, mentioned above, whose targets generally encode transcription factors and tranporters. Among them, miR159 was up-regulated in response to water limitation and was confirmed to target MYB transcription factors (myb33 and myb101) in Arabidopsis under drought stress in response to ABA accumulation [81]. However, in contrast to Arabidopsis, miR159 was down-regulated in rice [80] and peach root tissue. Although some members of MYB-family transcription factors were found in peach transcriptome libraries (at GDR; Genome database for Rosaceae), we could not determine the Myb transcription factors as targets for miR159 and this result may be consistent with previous findings that miR159 target was not related to MYB in tomato [82]. Another miRNA, miR160 is known to target three Auxin Response Factors (ARF 10, ARF 16 and ARF 17) in Arabidopsis [83] and it was reported that ARFs are transcription factors binding to auxin-responsive promoter elements to induce or repress auxin-regulated transcription [84] during the plant development such as root development and branching. After drought treatment, miR160 was up-regulated in peach roots and this miRNA (miR160a and miR160b) was also up-regulated in drought-tolerant cowpea cultivar in response to drought [69]. Thus, the upregulation of miR160 could be important in drought responses among different plant species. It was revealed that the miR169 family members were associated with drought response and high salt stress [52,55]. The miR169 targets a gene family encoding the alpha subunit of CCAATbinding NFY transcription factors (NFYA) requiring for adoptation to drought stress [8]. Two members, expression level of miR169a and miR169c were substantially down-regulated in Arabidopsis via ABA-dependent way [52] and also showed that level of miR169 was decreased under drought stress in M. truncatula by using high-throughput sequencing and qRT-PCR methods [66]. In addition, expression level of miR169 increased in two cowpea genotypes [69]. These results show good concordance with our findings that miR169 was down-regulated in peach after treatment. However, miR169g was up-regulated in rice during drought stress because the miR169g promoter contains two putative DRE (dehydration-responsive) cis-elements, causing the upregulation in response to drought and cold [53]. As an abiotic stress, drought disturbs the balance between ROS production and ROS elimination and thus leads to ROS accumulation in plant cells, which damages nucleic acids, oxidizes proteins and causes lipid peroxidation [85]. The detoxification of ROS radicals were carried out by Superoxide dismutases (SODs). The miR398 regulates the expression level of two Cu/Zn superoxide dismutases (cytosolic CSD1 and chloroplastic CSD2) under drought stress and level of miR398 was down-regulated in both M. trunctula [66] and maize [70], whereas its up-regulation was found in Triticum dicoccoides after 8-h stress [71]. In this study, we found that the level of miR389 was down-regulated in peach after stress treatment. A recently published paper [62] showed that miR408 are upregulated in response to water deficit in M. truncatula by targeting plantacyanin, and this miRNA was also induced in leaf tissue of H. vulgare under dehydration stress [9]. Contrary to M. truncatula and H. vulgare, expression level of miR408 decreased in peach and this result is consistent with previous finding [67] where they also detected the induction in expression of miR408 upon drought stress in O.sativa. However, it is necessary to specify that plantacyanin, putative target of miR408, was not found in peach transcriptome library during the computational target prediction process, therefore experimental methods such as RLM-RACE or degradome sequencing may be used for accurate determination of miR408 target.

Conclusions
In the present study, we genome-widely identified miRNAs and their expression pattern of drought-responsive miRNAs in roots and leaves of P. persica by using high-throughput sequencing. The expression level of 262 (104 up-regulated, 158 down-regulated) of the 453 miRNAs for LC/LS and 368 (221 up-regulated, 147 down-regulated) of the 465 miRNAs for RC/RS changed significantly in response to drought. Among them, drought responsive miRNAs (miR156, miR164, miR166, miR168, miR169, miR171, and miR395) were detected and their expression levels were measured by qRT-PCR. Our research also represents the first concerted effort to determine the large-scale small RNA datasets of peach. After sequencing, we identified a total of 531, 471, 535 and 487 known miRNAs and a total of 197, 221, 238 and 265 novel miRNA candidates for LC, LS, RC and RS tissues, respectively. Most of putative target transcripts for these miRNAs in biological process were related to nucleic acid binding (transcription factors) and catalytic activity. These results will greatly contribute to the understanding of post transcriptional gene regulation response to drought stress in peach.

Materials and Methods
Plant materials and stress treatment P. persica cultivar Francoise plants obtained from in vitro culture clones were grown in plastic pots for one month [86]. Droughtstress treatments were applied to the plants with similar stem lenght and leaf area for one week by withholding the water until leaves-wrinkle in greenhouse conditions as +24/18uC day/night 16/8 h light/dark. Then, root and leaf tissue samples were collected and immediately frozen in liquid nitrogen and stored at 280uC until RNA isolation.

Total RNA Isolation, small RNA library construction and sequencing
Total RNAs from leaves and roots of control samples and plants exposed to drought stress used in this study was isolated using an TriPure Isolation Reagent (Roche) according to the manufacturer's protocol. The quality and quantity of purified RNA were assessed by using a Nanodrop ND-2000c spectrophotometer (Nanodrop Technologies, Wilmington, DE, USA). Then, DNase treatment was carried out as described before [9] and all samples were stored at 220uC for miRNA quantification. For each control and stress treatment samples, equal amount of total RNA was pooled from three biological replicates to generate enough RNA (approximately 1000 ng) for deep sequencing. P. persica small RNA library construction, cluster generation and deep sequencing were carried out by the BGI (Beijing Genomics Institute, Hong Kong). Briefly, the isolated total RNAs of each sample was resolved on a denaturing 15% polyacrylamide gel for size selection and these small RNAs (#30 bases) were ligated to a pair of Solexa adaptors at the 59 and 39 ends using T4 RNA ligase. After ligation and purification, adapter-ligated small RNAs were reverse transcribed and 15-cycles amplified (cDNA RT-PCR) with a pair of adapter complementary primers in order to produce sequencing libraries. Then, PCR products were purified and were directly sequenced using Illumina HiSeq 2000 instrument according to manufacturer's recommendation (BGI, Shenzhen, China).

Bioinformatics analysis of sequencing data and novel miRNA prediction
After the sequencing reactions were complete, the high-quality small RNA reads ranging from 18 to 30 nucleotides were obtained from raw data analysis pipeline including; removing the low quality tags and trimming adaptor sequences so as to identify conserved and novel miRNAs in P. persica. Then, small RNA reads were used as queries to search against the Rfam family database and NCBI Genbank database to remove non-coding RNAs such as rRNA, tRNA, snRNA, snoRNA; the remaining sequences were searched against the miRBase database v18.0 with up to two mismatches to idnetify ''conserved'' mature miRNA orthologs. small RNAs not mapped to any miRNAs in miRBase were subsequently analyzed for potenial novel miRNAs by the program MIREAP (developed by BGI) with default parameters for mapping the peach genome and obtaining all candidate precursors with hairpin-like structures of novel miRNA candidates. Secondary structures of novel miRNAs were also checked using Mfold 3.2 [87].

Identifying miRNAs respoisve to drought treatment
In order to identify drought-responsive miRNAs and determine their genome-wide expression changes in response to drought, we first compared the gene expression patterns of miRNAs in control and the drought-treated leaf and root tissues in peach. Towards this purpose, we considered the following criteria; (i) adjusted pvalue should be less than 0.01 (p-value,0.01) in at least one data set, (ii) fold change or log 2 ratio of normalized counts between drought and control libraries was greater than 1 or less than 21 in one of the libraries. The frequency of miRNA read counts was normalized as transcripts per million (TPM) and normalization of miRNA expression levels between control and drought-stresses samples was carried out based on the following formula: Normalization formula : Actual miRNA count = Total count of clean reads |10 6 Actual read counts and normalized counts for each miRNA in each library are provided (Table S6). Afterwards, the fold-change between treatment and control and P-value were calculated from the normalized expression using the formula shown below: Fold-change formula : Fold change~log 2 treatment=control ð Þ P À value formula : For miR156, miR164, miR166, miR168, miR169, miR171 and miR395, the miRNA stem-loop reverse transcription reaction was performed in a volume of 10 mL containing 2, 20, and 200 ng of total RNA samples of leaf and root samples (1 mL), 0.5 mL 10 mM dNTP mix, 1 mL stem-loop RT primer (1 mM) and 7.5 mL nuclease free water. All those components of the reaction were mixed singly for the different dilutions of total RNA stem-loop RT primer cDNA syntheses and incubated for 5 min at 65uC and then put on ice for 2 min. Then, 4 mL first strand buffer (56), 2 mL 0.1 M DTT, 0.1 mL RNAseOUT (40 units/mL) and 0.25 mL SuperScript III (200 U/mL; Invitrogen) were added onto each tube and RT reaction temperature program set as follows: 30 min at 16uC, 60 cycles at 30uC for 30 s, 42uC for 30 s, and 50uC for 1 s. The RT reactions were terminated at 85uC for 5 min. During cDNA synthesis for miRNA quantification, we also generated additional reaction tubes including all components without RT primer (no-RT) and RNA template (no-RNA) as control reactions.
To experimentally validate predicted P. persica miRNAs and to measure and compare the expression levels of the miRNAs in leaf and root tissues on different stress conditions, qRT-PCR was carried out via FastStart SYBR Green Master mix (Roche) on The LightCycler H 480 II Real-Time PCR (Roche). By using previously synthesized 2 mL RT stem-looped cDNA products, quantitative PCR reactions were performed as followed; 10 mL 26 PCR Master mix, 1 mL forward (10 pmol), 1 mL reverse (10 pmol) primers, 0.3 mL (30 nM) reference dye and 10.7 mL nuclease-free water were mixed. With specifically designed forward primers for each individual miRNAs, the universal reverse primer (59-GTGCAGGGTCCGAGGT-39) [30] (Table S8) was designed for all the quantifications. Specified qRT-PCR thermal setup was adjusted as follows: 95uC for 15 min, followed by 40 cycles of 95uC for 5 s, 56uC for 10 s and 72uC for 30 s. All PCR products were denatured at 95uC and cooled to 65uC and the fluorescence signals were accumulated consistently from 65uC to 95uC as the temperature increased at 0.2uC per second. The reactions were repeated at least three times for credible statistical analysis.

Target transcript validation
In order to validate and detect expression level of predicted miRNA target genes which are related to drought stress, qRT-PCR was performed with a number of gene-specific primers. The target transcripts of Ppe-mir156, Ppe-mir166, Ppe-mir168, Ppe-mir169, Ppe-mir171, and Ppe-mir395 were obtained using psRNATarget (user-submitted transcripts and miRNA option) and BlastN algoritms. Specific PCR primers were designed using Primer3Plus software and also, primer dimers and hairpin formations were checked with the Autodimer program (Table  S8). At first, complementary DNA (cDNA) was generated from 1500 ng RNA using Superscript III First-Strand Synthesis System (reverse transcriptase, Invitrogen) according to manufacturer's instructions. In brief, the qPCR was performed in a 96-well plate instrument (LightCyclerH 480 Instrument II) and in 20 ml reactions that contained 1 ml of this cDNA, 300 nM each of specific forward and reverse primer, and FastStart SYBR Green I master mix (Roche). Each sample was run in biological and technical triplicates for each gene and relative quantity of these target transcrips calculated based on the housekeeping gene 18s rRNA (forward: GTGACGGGTGACGGAGAATT/reverse: GA-CACTAATGCGCCCGGTAT) as a normalizer. The qRT-PCR conditions were as follows: preheating for 10 min at 95uC before 40 cycles of 95uC for 30 s, 55uC for 1 min followed by 72uC for 10 min. To eliminate false-positive results, the melting curves of the gained real-time PCR data were analysed for each run and the data of the fluorescence signal were obtained from 55uC to 95uC as the temperature increased at 0.5uC per second.

Target Sequence annotation, Gene Ontology (GO) classification and KEGG pathway mapping
Because the majority of plant miRNAs have perfect or nearperfect complementarity with their target sites, the computational methods for finding the putative targets of miRNAs is the preferred way for prediction of conserved and novel peach miRNAs. Therefore, putative mature miRNA sequences were used as query to search against the Prunus persica EST database and high quality cDNA sequence by using BlastN search (http://www. plantgdb.org/XGDB/phplib/download.php?GDB = Pe). Alignments between each miRNA and its putative miRNA target(s) should meet certain criteria as follows; (i) No more than four mismatches between miRNA and its target (G-U bases count as 0.5 mismatches), (ii) No more than two adjacent mismatches in the miRNA/target duplex, (iii) No adjacent mismatches in in positions 2-12 of the miRNA/target duplex (59 of miRNA), (iv) No mismatches in positions 10-11 of miRNA/target duplex, (v) No more than 2.5 mismatches in positions 1-12 of the of the miRNA/ target duplex (59 of miRNA) as noted by Allen [88] and Schwab [11]. The functional annotation and categorization of identified putative miRNA targets were determinated using the Blast2Go (B2G) software suite v2.3.1 with the default parameters (http:// www.blast2go.com/b2ghome) [89]. Beside, these putative miRNA target sequences were used as query against the KEGG (Kyoto Encyclopedia of Genes and Genomes) database using the KEGG automatic annotation server [90] in order to reveal their biological function in various cellular metabolic pathways. With the aid of KAAS annotation tool, an ortology number (KO) in database was assigned to the genes within KEGG Genes database based on the sequence similarity comparisons and then, the KO numbers associated with the corresponding unique KEGG gene were used for mapping one of the KEGG's reference metabolic pathways.