Small RNA and transcriptome deep sequencing proffers insight into floral gene regulation in Rosa cultivars

Background Roses (Rosa sp.), which belong to the family Rosaceae, are the most economically important ornamental plants—making up 30% of the floriculture market. However, given high demand for roses, rose breeding programs are limited in molecular resources which can greatly enhance and speed breeding efforts. A better understanding of important genes that contribute to important floral development and desired phenotypes will lead to improved rose cultivars. For this study, we analyzed rose miRNAs and the rose flower transcriptome in order to generate a database to expound upon current knowledge regarding regulation of important floral characteristics. A rose genetic database will enable comprehensive analysis of gene expression and regulation via miRNA among different Rosa cultivars. Results We produced more than 0.5 million reads from expressed sequences, totalling more than 110 million bp. From these, we generated 35,657, 31,434, 34,725, and 39,722 flower unigenes from Rosa hybrid: ‘Vital’, ‘Maroussia’, and ‘Sympathy’ and Rosa rugosa Thunb. , respectively. The unigenes were assigned functional annotations, domains, metabolic pathways, Gene Ontology (GO) terms, Plant Ontology (PO) terms, and MIPS Functional Catalogue (FunCat) terms. Rose flower transcripts were compared with genes from whole genome sequences of Rosaceae members (apple, strawberry, and peach) and grape. We also produced approximately 40 million small RNA reads from flower tissue for Rosa, representing 267 unique miRNA tags. Among identified miRNAs, 25 of them were novel and 242 of them were conserved miRNAs. Statistical analyses of miRNA profiles revealed both shared and species-specific miRNAs, which presumably effect flower development and phenotypes. Conclusions In this study, we constructed a Rose miRNA and transcriptome database, and we analyzed the miRNAs and transcriptome generated from the flower tissues of four Rosa cultivars. The database provides a comprehensive genetic resource which can be used to better understand rose flower development and to identify candidate genes for important phenotypes.

The development of molecular markers for roses began with the first molecular linkage map covering over 300 markers from Rosa multiflora hybrids [4], and several genetic maps were constructed recently [5,6]. However, the genetic resources for roses are relatively weak compared to those for other Rosaceae families. As of June 2011, approximately 4,834 unigenes were available. These were generated from 9,289 expressed sequence tags (ESTs) for rose in the Genome Database for Rosaceae (GDR) [7,8]. These unigenes cover only 7.6% of apple genes, 13.89% of strawberry genes, and 16.84% of peach genes. Clearly, more abundant transcriptomic resources generated from different roses are needed to allow for the investigation of key traits, including resistance to disease or stress, flower morphology, and scent [9,10].
Transcriptome sequences are often analyzed from both model and non-model plants to monitor whole gene expression. Whole gene expression is useful to identify biotic [11] or abiotic stress related genes [12,13], understand organ development [14,15] and characterize differential traits between closely related species of rose [16]. Next-generation sequencing (NGS) technologies, such as Illumina, SOLEXA, Genome Sequencer FLX system (GS FLX), and ABI SOLiD, allow analysis of the transcriptome because of increased throughput and reduced sequencing cost [17,18]. The GS FLX is considered by many to be the most powerful platform to analyze protein-coding sequence data with strengths, including long reads, good accuracy, and the ability to support ultra-high-throughput analysis [19]. Because of these strengths, GS FLX is often applied to generate transcriptome data (summarized in Table 1 in [20]). Transcriptomic data in ornamental plants like roses expands our knowledge of the genetic control of flower quality. These findings can be applied in the floricultural industry to advance efforts to screen economically important phenotypes [21]. Here, we generated transcriptomes of four Rosa cultivars using GS FLX sequencing to compare floral development and other features.
MicroRNAs (miRNAs) are short (20-24 nt) nonprotein-coding RNAs [22,23], which play essential roles in regulating plant growth and development [22]. miRNA genes, called MIR genes, are transcribed by RNA polymerase II (Pol II), and the miRNA transcripts form self-complementary fold-back structures called primary miRNA (pri-miRNA). Pri-miRNAs are processed to mature RNAs (miRNA/miRNA* duplex) through cleavage by Dicer-like 1 (DCL1) protein [23]. After release into the cytosol, miRNAs bind near-perfectly to their target mRNAs, and the remaining strands (miRNA*) are degraded. miRNAs regulate expression of target mRNA post-transcriptionally through either cleavage of the target mRNA or translational inhibition. miRNAs are often identified by cloning or by bioinformatic alignment to known miRNAs. However, several strategies have been developed to computationally identify miRNAs from deep sequencing data [24][25][26][27]. These algorithms, which were initially designed for vertebrate miRNAs, are now able to predict plant miRNAs considering different features of plant miRNAs [28,29]. There are several databases of annotated and predicted miRNAs for plants that can be used in bioinformatic approaches for miRNA identification, including miRBase [30] and the Plant miRNA Database (PMRD) [31]. These recent developments of high-throughput sequencing techniques and analysis tools, as well as databases for miRNAs, challenge us to identify miRNAs from various plant species.
In this study, we constructed a Rose database of transcriptomes and miRNA sequences and annotations generated from flower tissues of three Rosa hybrid cultivars (Vital, Maroussia, and Sympathy) and R. rugosa Thunb. (Haedang). We included data from the three sequenced Rosaceae genomes (apple, strawberry, and peach) and the grape genome. In the rose database, users can analyze gene contents, gene families, phylogenetics, and miRNA profiles from a publicly available website (http:// 210.218.199.249/rose/). These data can be used to identify candidate miRNAs and target genes that may regulate flower development and hormonal regulation. It can be used to establish a genetic basis of valuable phenotypes of rose flowers. To demonstrate the utility of the database, we describe miRNAs that are highly conserved among Rosa ('Vital' , 'Maroussia' , 'Sympathy' , and 'Haedang'), as well as statistically abundant miRNAs in one rose than others. Several target genes of miRNAs were experimentally validated using 5' RACE, and most of them were transcription factors regulating flower development. Therefore, rose database provides a comprehensive resource to understand flower development and to identify new candidate genes for studying rose phenotypes.

Results
Transcriptome sequence assembly from rose flower libraries  Table 1). The average read length was 219-226 bp for each library. Prior to assembly, we masked low-complexity and poly (A/T) sequences and removed reads shorter than 100 bp using the SeqClean program downloaded from the Dana Farber Cancer Institute (http://compbio.dfci.harvard.edu/tgi/software/) ( Figure 1). After pre-processing, 88-94% of raw data remained to be clustered and assembled into contigs (Table 1).

Comparative analysis and functional annotation of rose flower miRNAs
The scheme used to develop functional annotation for roses is depicted in Figure 1. We downloaded genes for three sequenced Rosaceae genomes (apple, strawberry, and peach) [7] and grape [33] to compare to rose unigenes. We annotated rose unigenes (including contigs, un-assembled singlets, and un-clustered singlets) and other protein-coding genes by implementing BLASTx against non-redundant (NR)-protein and Arabidopsis protein sequences. We found that 51%, 61%, 58%, and 51% of unigenes aligned to known and/or hypothetical proteins for 'Vital' , 'Maroussia' , 'Sympathy' , and 'Haedang' , respectively (Table 2). Although rose unigenes were annotated at a lower level than those of related species, (i.e. 81%, 71%, 91%, and 89% of genes were annotated for apple, strawberry, peach, and grape, respectively), approximately 64-77% of rose contigs were annotated. Using Pfam, we identified 2,262, 2,315, 2,383, and 2,369 distinct domains within the 'Vital' , 'Maroussia' , 'Sympathy' , and 'Haedang' unigenes, respectively (Table 2). These numbers are slightly less than the numbers of domains identified in apple (3,311), strawberry (3,333), peach (3,297), and grape (3,305) genes ( Table 2). The reason why smaller numbers of domains were identified in Rosa is probably due to the fact that the genes which contain unidentified domains were expressed at lower levels in rose flowers or only expressed in non-flower tissues.
To identify functional and evolutional relationships, we analyzed orthologs and paralogs of 294,936 gene/ unigenes from 4 Rosa, apple, strawberry, peach and grape, by implementing orthoMCL [34]. We found that 242,234 genes were clustered into 39,447 gene families, suggesting that 82.13% of total genes or unigenes have similarity among species (Figure 2A). Approximately 23% of the 39,447 gene families were homologous to Arabidopsis gene families in The Arabidopsis Information Resource (TAIR) database. The TAIR provides 996 gene families and 8,331 genes, covering 32.16% of Arabidopsis genes [35]. Gene family comparisons provide the number of conserved and shared genes in plants available in the rose database and the number of species-specific gene families for each species ( Figure 2B). On the whole, 4,881 gene families corresponding to 75,314 genes were conserved in 8 plants analyzed in this study ( Figure 2B), and 592 gene families corresponding to 7,451 genes were common to all Rosaceae species, indicating that conserved functions were specific to Rosaceae ( Figure 2B). In rose, 1,484 gene families with 6,808 genes were common to 4 Rosa and were not conserved in fruitbearing Rosaceae species ( Figure 2B). Additionally, 396, 107, 121, and 442 gene families were specific to Vital, Maroussia, Sympathy, and Haedang, respectively ( Figure 2B). Of the species-specific gene families, only 3-6% was annotated based on homology to Arabidopsis gene families, suggesting these may be the gene families evolved within each rose.

miRNA analysis
We sequenced flower small RNA (sRNA) sequences for roses using the Illumina GAIIx platform and generated an average of 9,958,816 sRNAs per rose (Table 3). After pre-processing, 9,698,841 (94.18 %) unique sRNA sequences, with 18-26 nucleotides (nt) in length, were remained for mapping. The sRNAs were dominantly 24nt, 21nt, and 22nt in length, whose production relied on DCL3, DCL4 and DCL2, respectively (Additional file 1). [36], as previously observed in other plant species [37,38]. In order to distinguish miRNAs from all other sRNAs, we predicted secondary structures of miRNAs by mapping 9.6 million sRNA sequences to the strawberry genome because there is no genome available for rose. 3,090,334 (31.86%) sRNA tags were perfectly mapped to the strawberry genome after discarding sRNAs that are aligned to more than 10 loci in the strawberry genome. We predicted secondary structures using RNAfold, and the sequences used for this prediction were 500 bp marginal sequences from mapped sRNA. We validated miRNAs which meet the criterion  Transcript reads were assembled by TGICL [32] after pre-processing with SeqClean. Assembled sequences were annotated by searching homologs against non-redundant (NR) database and Arabidopsis (ver. 10). Domains and Gene Ontology (GO) were analyzed by Pfam. KEGG metabolic pathways and plant ontology (PO) were assigned by referring the description of top matched Arabidopsis genes. miRNAs were identified by applying modified miRDeep [24,29]. To search for known miRNAs, sRNA sequences were BLAST searched against currently known miRNAs in miRBase [86] databases.

Conserved miRNAs in four roses and Rosaceae
We investigated conservation and/or variation of miRNAs in rose flowers by analyzing the presence or absence of miRNA tags among four Rosa cultivars ( Figure 4A). The miRNA tags denote individual miRNA with sequence redundancies. We found that two miRNA tags are uniquely expressed in Maroussia and Haedang, respectively. Most (297 miRNAs, 90%) of the miRNAs were conserved in the four Rosa cultivars ( Figure 4A). We also analyzed miRNAs conserved in Rosaceae. For this analysis, we retrieved sRNA sequences from Gene Expression Omnibus (GEO) database for strawberry (GSE34813), peach (GSE18764), and apple (GSE36065). miRNAs from these species were re-analyzed according  Distribution of orthologous and paralogous gene families, including Rosa and apple, strawberry, peach and grape, were analyzed by OrthoMCL [34]. A. A total of 294,936 sequences from the eight different organisms were clustered into 242,234 (82.13%) families and singletons using OrthoMCL. B. Numbers of families and genes (in parenthesis) presented in each organism combination are given in the individual sections. 4,881 families (with 75,314 genes) were common to all species analyzed in this study, 592 families (with 7,451 genes) were shared in Rosaceae group. 6,808 were clustered in rose specific group with a total of 1,484 families. The species-specific gene families denote genes clustered within singe organisms.
to our pipeline (Additional file 3). As a result, a total of 123 miRNAs were conserved among Rosaceae. Among them, 13, 7, 12 and 7 miRNAs were uniquely expressed in rose, strawberry, peach and apple, respectively ( Figure 4B). Thirty miRNAs were conserved in all Rosaceae ( Figure 4B). Strawberry has the largest number of miRNAs (17) conserved in Rosa flower, comparing to those in peach (2) and apple (7) ( Figure 4B), which suggest that strawberry is the most closely related species of Rosa.

Regulatory roles of conserved miRNAs in Rosa
Plant miRNAs are negative regulators of gene expression, and thus, play essential roles in developmental patterning [22]. Therefore, we hypothesized that the conserved miRNAs in four Rosa cultivars may have evolutionally conserved functions. To address this hypothesis, we analyzed conserved miRNAs whose target genes were commonly identified in Rosa, and miRNAs with higher reads were representatively shown (Table 4).
With these criteria, we selected fifteen miRNA families, two of which were novel. These miRNAs may have conserved function, but it is necessary to experimentally validate whether these miRNAs actually regulate computationally predicted target genes. For the target validation, we selected 7 miRNA families and 15 target genes whose penalty score were four or less (Additional file 4). In the specific circumstance where the expression of a miRNA target gene is regulated through miRNA-directed cleavage mechanism, 5 0 RNA ligase-mediatedrapid amplification of cDNA ends (5 0 RLM-RACE) is generally employed to confirm such targeting. Therefore, it is necessary to experimentally validate whether computationally predicted targets are actually regulated by miR-NAs by means of miRNA-directed cleavage of these targets because plant miRNAs regulate their target genes mainly by cleaving them [42,43]. Therefore, we isolated total RNAs from four Rosa cultivars and performed 5' RLM-RACE for squamosa-promoter binding protein (SBP), MYB, APETELA2 (AP2), no apical meristem (NAC), F-BOX, and auxin response factor (ARF). As a result, nine target genes of seven miRNAs were successfully validated ( Figure 5). The target gene, SBP transcription factor, was validated to be target of miR156/157 families in three Rosa cultivars (Maroussia, Sympathy and Haedang). We observed that the SBP transcription factor was targeted both by miR156 and miR157 families in 'Sympathy' , whereas other SBP transcription factors were targeted only by miR156 family in Haedang and 'Maroussia' ( Figure 5A). miR156 and miR157 in plants have been grouped in one miRNA family due to their high degree of sequence similarity and their conserved target, the SBP transcription factors [22]. Additionally, for other miRNAs (miR159, miR172, miR164, miR394, and miR160 families), we confirmed miRNA-directed cleavage in one or two Rosa cultivars ( Figure 5B-F). Considering the function of target genes, most of targets validated in 5' RLM-RACE assay were transcription  (105) Novel miRNAs (Unique) 28 (21) 31 (23) 28 (21) 25 (20) 33 (25) Found by homology search with miRBase a 136 137 137 128 137 Total miRNA groups (Major) 78 (40) 82 (42) 80 (40) 75 (41) 84 (53) Conserved miRNA groups (Major) 63 (27) 65 (28) 65 (28) 61 (28) 65 (37) Novel miRNA groups (Major) 15 (13) 17 (14) 15 (12) 14 (13) 19 (16) a The number of miRNAs searched by BLAST against known miRNAs in the miRBase [86]. Flavonoids synthesize diverse secondary metabolites determining flower colour such as antocyanins, flavonols, flavones and proantocyanidins. We aligned Rosa unigenes against Arabidopsis genes involved in flavonoid or carotenoid biosynthetic pathways. Subsequently, we analyzed miRNAs whose target genes are related to those pathways. We identified thirteen miRNA families putatively targeting 27 target genes (Table 5). To analyze miRNAs regulating colour-related pathways, we examined miRNA enrichment in given Rosa with p-values less       (Table 5). The miRNA enrichment means that statistically over-expressed miRNA tags in a given rose cultivar comparing with miRNA tag constructed in another cultivar. Five of the miRNAs (miR171, miR166i, miR159e, miR845, and miR396e) were enriched in the white flowers of Maroussia. Especially, miR396e, which is predicted to target Cytochrome P450 (CYP), had two fold higher read counts than other roses. CYP is an essential gene for synthesizing red colour in flowers, fruit, and epidermal tissues [44,45]. It is suggested that these miRNAs may negatively regulate target genes to prevent accumulation of carotenoids or flavonoids, resulting in white flowers. However, no target genes with penalty scores of four or less were predicted.  The bold characters denote enriched miRNA tags with less than p<0.001 by analyzing Audic's test [88]. Thus, we were not able to confirm the miRNA-directed cleavage of target genes involved in colour determination. However, among the target genes, the SPL and R2R3-MYB transcription factors, both of which are known to negatively regulate flavonoid biosynthesis, were experimentally validated to be targets of miR156 and miR159, respectively [46,47]. In addition, the target gene of miR159 was predicted only in Maroussia (white) and Haedang (pink), which indicates that the colours of the rose flowers, may be tightly regulated via complex mechanism of various miRNAs in nature (Additional file 5).
In addition to miRNA profiles, we also analyzed differential expression of colour-related genes involved in carotenoid and flavonoid biosynthesis ( Table 6). We identified 28 (with 68 reads), 15 (15), 31 (76) and 28 (100) unigenes for carotenoid biosynthesis (Additional file 6) and 31 (96), 12 (12), 40 (232), and 38 (201) Unigenes for flavonoid biosynthesis (Additional file 7) for Vital, Maroussia, Sympathy and Haedang, respectively. Based on these data, the number of colour-related unigenes in Maroussia is lower than other roses and fruit-bearing Rosaceae. Moreover, all colour-related unigenes in Maroussia were singlets due to smaller number of sequence reads, whereas unigenes in other roses were assembled into contigs. In summary, we found that some of the miRNAs were enriched in Maroussia, and the smaller number of sequencing reads of colour-related genes were found in Maroussia, which raises an intriguing possibility of miRNA-directed regulation of colour-related gene in Maroussia.

Discussion
Rose breeders select roses according to particular criteria, which include cold and disease resistance, flower form, recurrent flowering, and to some degree, scent [21]. In spite of the importance of phenotypes for roses, only a few studies have addressed flower features at the molecular level. Prior to this study, only a few thousand rose unigenes had been deposited in the EST database (dbEST) at NCBI (search query: "Rosa"[Organism]). With the aim of providing valuable resources for molecular studies of rose flowers, we constructed the Rose transcriptome database and analyzed miRNA sequences from three R. hybrids (Vital, Maroussia, and Sympathy) and Haedang.
These rose cultivars were selected because of a range of structures (from a single to multiple layers of petals) and colours (from white to dark red) in their flowers (Figure 1). The database includes 35,657, 31,434, 34,725, and 39,722 unigenes for 'Vital' , 'Maroussia' , 'Sympathy' , and 'Haedang' , respectively (Table 1), covering an average of 55.96%, 60.70%, 59.14%, and 61.57% of apple, strawberry, peach, and grape protein-coding genes, respectively. We annotated unigenes up to similar levels of relatedspecies, especially in contigs (Table 2). However, overall numbers of identified domains were slightly lower than those identified in related species (Table 2). These domains may be expressed at lower levels in rose flowers or not expressed in flower tissues. In addition, the Rose database provides additional information for the unigenes, including gene description, domains, GO, PO, metabolic pathway, FunCat, and homologous gene families among roses, Rosaceae families (apple, strawberry, and peach) and grape (Figure 1). The sequence dataset and highquality annotation enabled us to analyse and compare the transcripts which are specialized in rose flowers and related species.
With the help of large-scale sRNA sequencing, a large number of miRNAs have been identified, characterized, and reported [48][49][50][51][52][53][54]. However, these techniques are best applied to plant species whose genomes have been assembled, annotated and published (i.e., Arabidopsis, Rice and Grape [51,52,54]) because it is necessary to predict and construct hairpin precursors of potential miRNAs, using neighbouring genomic sequences of the mapped sRNAs to distinguish high quality miRNAs from other sRNAs [24,28]. Therefore, in the absence of the genome sequences, the strategy to identify miRNAs is limited [55,56].
In this study, we set out to identify highly qualified conserved and novel miRNAs from floriculture plants Rosa × hybrida (Vital, Maroussia, and Sympathy) and Haedang, for which genomes are not available. To predict hairpin precursors for potential miRNAs, we mapped sRNA tags to the strawberry genome, the most closely related species among genome sequenced plants, and identified mature miRNAs using modified miRDeep, a program that employs a probabilistic model of miRNA biogenesis [24,29]. Mapping sRNA to cross-species is possible because secondary structures for pre-miRNAs are highly conserved between species [22]. Identification of miRNAs by mapping of sRNA to related species has strengths in that it i) maximizes the screening of novel miRNAs in plants for which genome is not available, ii) allows identification of the conserved and novel miRNAs between related species, and iii) provides an incremental improvement in accuracy of miRNA identification by utilizing secondary structures. We might have missed some novel miRNAs that are specific to roses; however, these novel miRNAs are less likely to be important in functional studies because they are expressed at low levels and lack target-genes, whereas conserved miRNAs are usually highly expressed [23,57].
With this approach, we identified 122, 125, 125, and 120 miRNAs for 'Vital' , 'Maroussia' , 'Sympathy' , and 'Haedang' , respectively (Table 3). Of these miRNAs, 21, 23, 21, and 20 were novel miRNAs for 'Vital' , 'Maroussia' , 'Sympathy' , and 'Haedang' , respectively (Table 3). These novel miRNAs correspond to 16-18% of the total identified miRNAs. We also identified 101, 102, 104, and 100 conserved miRNAs by aligning sRNA tags to known miRNAs for 'Vital' , 'Maroussia' , 'Sympathy' , and 'Haednag' , respectively (Table 3). This hybrid method maximized the identification of expressed miRNAs for rose flowers, therefore providing many miRNA candidates (Figure 1). In addition, this allows us to understand the conserved and unique regulatory processes occurring in rose flowers of different species and hybrids. miRNA research is advancing from analysis of single tissue or species to comparison of miRNAs between species [50], varieties [58], tissues [48,[51][52][53]59], or different development stages [54]. Today, NGS technologies allow for the comparison of miRNA profiles with statistical analysis of the redundancy of miRNA tags and enables discussion of tissue-or organ-specific miRNAs [51,52,54,58,59]. We first compared miRNA profiles among three R. hybrida cultivars and Haednag to analyze the conservation and variation of miRNA tags in the four roses ( Figure 4). We identified miRNA tags detected in single rose cultivar, leading to specific function to Rosa. We also identified 297 conserved miRNAs in all four roses ( Figure 4A), including 274 known miRNAs and 23 novel miRNAs. Among 274 conserved miRNA, 70 miRNA tags were also identified in strawberry. In addition, 30 conserved miRNAs were also conserved in all Rosaceae ( Figure 4B), suggesting that Rosaceae share many conserved miRNAs between ornamental and fruit-bearing plants. Furthermore, among the conserved miRNAs, seven miRNA families were experimentally confirmed to regulate their target genes by cleavage mechanism using 5' RACE assay ( Figure 5), suggesting that the conserved miRNAs identified in this study are actually functional.
Most of the target genes validated in this study were transcription factors ( Figure 5A-D, F), and their mutant phenotypes were characterized in many model plants (Table 4). Over-expression of miR156 ( Figure 5A) and miR159 ( Figure 5B) induced delayed flowering in Arabidopsis by negatively regulating SPL and MYB family transcription factors genes, respectively [60,61]. The expression of miR167-resistant ARF6 ( Figure 5A) leads to arrested ovule development and indehiscent anthers [62]. miR172 ( Figure 5C) is crucial for development of reproductive organs and for timely termination of floral stem cells by regulating AP2 RNA stability [63]. The expression of miR172-resistant AP2 induces the formation of variable numbers of floral organs with numerous petals and lacking inner whorl organs [63,64]. miR160 regulates development by altering expression of auxininduced genes through ARF families [65,66]. miR164 targets CUC1 and CUC2 transcripts in Arabidopsis and controls leaf margin development. Therefore, we assumed that the conserved miRNAs, which were experimentally validated in this study, may have important roles in floral organ identity or flower developments.
In addition, these eight miRNAs are evolutionary conserved and abundantly expressed miRNAs in roses. According to previous studies, miR156, miR159, and miR160 are evolutionary conserved in all land plants, and miR164, and miR172 are conserved in seed-bearing plants [57]. Evolutionarily conserved miRNAs in plants tend to regulate ancestral transcription factors that specify basic meristem functions, organ polarity and separation, cell division, or hormonal control (reviewed by Garcia [67]). Based on experimental validation of conserved miRNAs and the current discussions [22,23,67], we might expect that novel and un-validated miRNAs identified in this study (Table 4) possibly play important roles such as flower development or hormonal control. Our analysis suggests that the Rose database is a useful tool to search for candidate target transcripts or miRNAs that play roles in flower development in rose and for those have a variety of other specific functions.
Flower colour in most angiosperms is one of the most important targets for plant breeders and many differentcoloured cultivars have been bred using natural mutants or genetically-related species. Flower colours are determined by an accumulation of secondary metabolites such as flavonoids, carotenoids, and betalains [68,69]. We examined miRNA profiles in which target genes were involved in colour-metabolite related biosynthesis pathways to gain insight into the regulation of the white flowered cultivar, Maroussia (Figure 1), which seems to be regulated by miRNAs (Table 5). We hypothesized that miRNA enrichment in Maroussia may negatively regulate the colour-related genes, leading to white colour of Maroussia. Although miRNA enrichments were also observed in other Rosa (Vital, Sympathy, Haedang), the most interesting thing was that five miRNA tags were enriched in Maroussia. Especially, miR356e was expressed up to five times more than other Rosa. Thus, the function of its target, CYP, which involves catalyzing the biosynthesis of flavonoids and cyanidin (red to magenta) and delphinidin (violet to blue) [44,45], would be more negatively regulated in Maroussia, which may possibly lead to lack of colour . Unfortunately, we were not able to validate miRNA-directed cleavage of these targets due to high penalty score of the target genes (Table 5). However, we would rather expect that miRNA-directed regulation of target gene possibly lead to low read counts of target genes (i.e. low level of expression) in the transcriptome library, which makes target identification more challenging.
However, previous studies reported that two transcription factors, SPL and R2R3-MYB, both of which regulate expression of antocyanins-related genes. Moreover, overexpressed miR156 directly prevent the expression of anthocyanin biosynthetic genes (Additional file 5) by targeting SPL9, in Arabidopsis [47]. In this study, we identified nine miR156 members from all Rosa (Additional file 2), and their target genes, SPL transcription factors, were experimentally validated by 5' RACE assay ( Figure 5). The miR159 were among the most frequent in our library (187,579; 271,208; 264,412 and 327,436 for 'R.thunb.' , 'Marcia' , 'Sympathy' , and 'Vital' , respectively) and its sequencing frequencies were 10 to 100 times more than other relatively abundant miRNA families, including miR156, miR157, and miR167 (Table 4). Along with this, it has been previously reported that differential expression of R2R3-MYB gene determine colour patterning in plants that are linked with anthocyanin production [70,71]. Given the fact that the miR159 is very highly expressed in Rosa and miR159-directed cleavage of R2R3-MYB gene is confirmed using 5' RACE, our results raise an intriguing possibility that miRNAs in roses may be involved in pigment synthesis pathway. In addition, miR828, and miR858, which are also involved in R2R3-MYB regulation [72] was predicted target MYB genes in Rosa (Additional file 2).
Based on the miRNA profiles (Table 5, Additional file 5) and the previous research, flower colour seems to be regulated by combinatorial mechanisms. Therefore, it is difficult to elucidate the molecular mechanism or miRNA-directed regulation of colour determinacy of Rosa flowers. Nevertheless, miRNA profiles give potential clues to examine the colour of Rosa flowers. Moreover, there are limited numbers of unigenes currently available. Therefore, we would expect to reveal more direct evidence of miRNA-mediated regulation of colour development in Rosa when genome resources become more available.
In addition to the miRNA profiles, we also compared the expression profiles of colour-related transcripts in Rosa and Rosaceae families ( Table 6). It shows that Maroussia contains less number of unigenes (similar to strawberry or peach) than other Rosa cultivars. The higher number of colour-related gene in apple arose by current whole genome duplication [2]. In addition, the average number of reads per gene in Maroussia was also smaller than other Rosa cultivars because all of them were singlets (Table 6). On the whole, based on expression profiles of miRNA and the smaller number of transcripts involved in colour-metabolite related biosynthesis in Maroussia, it is possible to expect that colour-related genes are systemically repressed in Maroussia.

Conclusions
The Rose database is the first database that allows for the comparison of transcripts and miRNAs for the flowers of Rosa hybrida and Haedang, containing more than 30,000 transcript sequences for each rose and more than 300 conserved and novel miRNAs with high quality annotations. These data and analysis are available at http://210.218.199.249/rose/. We identified several miRNA families that are conserved and highly expressed in four Rosa, including novel miRNAs. We also present preliminary data which raises intriguing possibilities of miRNA-directed regulation in the white colouration of the Maroussia flower. Our study demonstrates that the Rose database may facilitate a comprehensive understanding of rose flower development, and is a valuable genetic resource for the analysis of gene functions and regulatory pathways that determine flower phenotypes.

Rose flower materials and RNA preparation
The flowers of two red roses (Vital and Sympathy) and white one (Maroussia) at various developmental stages (from bud to open flower) were purchased from a rose farm located in Goyang, GyeongGi-Do province, South Korea. Flowers at various developmental stages from vegetatively propagated R. rugosa Thunb were collected from natural habitats.
A modified Suzuki's method [73] was used for total RNA extraction as follows: Rose flowers were ground to a fine powder in liquid nitrogen using a mortar and pestle. The powder was transferred to a tube containing extraction buffer (100 mM Tris-HCl pH 6.8, 25 mM EDTA pH 8.0, 2% (w/v) cetyltrimethylammonium bromide, 1.4 M NaCl, 2% (w/v) polyvinyl polypyrrolidone, 7% (v/v) 2-mercaptoethanol, and 0.5 g/l spermidine trihydrochloride) pre-heated to 65°C in a powder-to-buffer ratio of 1:3 to 1:5. The sample was mixed rapidly and incubated for 5 min at 65°C then cooled to room temperature. An equal volume of chloroform was added to the homogenate followed by vortexing for 1 min. The mixture was centrifuged at 12,000 rpm for 10 min at 4°C and the upper phase was transferred to a new tube. Chloroform extraction was repeated 4-6 times until the middle layer was clear. Then, 0.6 volumes of 10 M LiCl was added to the sample before incubation at −20°C for 1 h and subsequent centrifugation at 10,000 rpm for 20 min at 4°C. Pellets were resuspended with DEPC-treated water and 1/10 volume of 3 M sodium acetate (pH 5.5) and an equal volume of isopropanol was added before centrifugation at 12,000 rpm for 30 min at 4°C. The precipitated total RNA pellet was washed with 75% ethyl alcohol, airdried for 10 min, and dissolved in DEPC-treated water.

Sequencing of cDNA and small RNA
A total of 2 μg of mRNA isolated from total RNA using a PolyATract W isolation system (Promega, Medison, WI, USA) was converted to cDNA using the cDNA Synthesis Kit (Stratagene, La Jolla, CA, USA). Approximately 5 μg of cDNA was sheared by neubulization to produce random 300-800 bp fragments. Oligonucleotide adaptors were ligated to the fragmented cDNA samples. Fragments were denatured to generate single-stranded DNA that was amplified by emulsion PCR for sequencing. Sequencing was performed on a 454 GS FLX system (Roche, Fresno, CA, USA) at NICEM (Seoul National University, Korea; http://nature.snu.ac.kr). The raw data submitted to NCBI under accession no. SRA049095.
Small RNA libraries were constructed as described previously [74] with some modifications. Low molecular weight RNA was isolated from 200 mg of total RNA by PEG 8000/NaCl precipitation. Small RNAs (20-30 nt) were purified from 15% denaturing PAGE gels and ligated first with the 5' RNA adaptor and then with the 3' RNA adaptor provided by Illumina. In each step, the ligated product was PAGE-gel purified. After first-strand synthesis and 18 cycles of PCR amplification, the product was PAGE-gel purified and submitted for sequencing on an Illumina GAIIx at MACROGEN (Seoul, Korea; www.macrogen.com). The raw data and miRNA profiles submitted to NCBI under accession no.GSE39882.

Pre-processing and assembly of rose transcripts
In a pre-processing step, we masked low-complexity and poly (A/T) sequences, and removed reads less than 100 bp using the SeqClean program (downloaded from the Dana Farber Cancer Institute; http://compbio.dfci. harvard.edu/tgi/software/). For contig assembly, we clustered reads into groups using megablast with criteria of more than 30 bp alignments between reads and identities greater than 94%. We then applied CAP3 to assemble the clustered reads into contigs [32].

Transcript annotation and comparative analysis
The overall process and integrated databases are represented in Figure 1. To analyze functional annotation for each transcript, we aligned contig and singlet sequences with a lower expect value of 1e -10 against plant nonredundant (NR) proteins downloaded from the NCBI database [75] and against Arabidopsis proteins downloaded from TAIR (ver. 10; http://www.arabidopsis.org/) [76]. Based on the Arabidopsis annotation, we integrated pathway information referring to KEGG [77], FunCat [78], and plant ontology (PO) [79,80]. Domains were identified with hmmPfam in the InterProScan Package [81]. Gene ontology (GO) was also analyzed based on the hmmPfam annotation. To exploit functionally related genes among related species, we downloaded protein and CDS sequences of the Rosaceae family (apple, strawberry, and peach), from GDR (http://www.rosaceae.org/, [7]), and grape, from Grape Genome Browser (http://www. genoscope.cns.fr/externe/GenomeBrowser/Vitis/, [33]), and analyzed orthologs and paralogs with OrthoMCL [34].

miRNA analysis and miRNA target prediction
The scheme used for identifying miRNAs is presented in Figure 1. We consulted Breakfield et. al. paper for miRNA prediction [37]. To get high-quality small-RNA (sRNA) reads, we removed poor quality reads and adaptor sequences (5'-ATCTCGTATGCCGTCTTCTGCTTG-3') using the FASTX toolkit (downloaded from http:// hannonlab.cshl.edu/fastx_toolkit/download.html, [82]), from raw data (Table 3). We removed sRNA sequences shorter than 18 nt in length or containing ambiguous nucleotides. From cleaned sequences, we further removed sRNAs aligned to the non-coding RNAs (ncRNA) such as rRNAs, tRNAs, snRNAs and snoRNAs. We collected ncRNAs from Rfam [83], the Plant snoRNA Database [84], and 19,350 rRNAs from the NCBI database. tRNAs were predicted with tRNAScan-SE [85] for the strawberry reference genome. We predicted tRNAs from strawberry because it was used to predict secondary structure of miRNA. The remaining sRNAs were aligned to strawberry genomes. We used sRNAs for further analysis with following criterion 1) perfect match, 2) less than 10 locations of the strawberry genome. The secondary structures were constructed with RNAfold (for sRNAs more than 25 copy numbers). From validated sRNA sequence, we predicted proper mature-star hairpin pairing with 500 bp marginal sequence from mapped sRNAs. Among the sRNAs with the same 5' alignments, we selected representative sRNA with maximum read. We selected the mature-star sequences within the 15 bp margin and refolded by RNAfold. We applied a modified miRDeep program to predict Rosa miRNAs [24] and selected valid miRNAs suitable for plant miRNA annotation [39]. We confirmed conserved miRNA by searching homologous miRNAs in miRBase (http://www.mirbase. org/, ver. 17, [86]) and grouped miRNA families. During the sRNA mapping to the strawberry genome, several putative miRNAs were removed because of lower similarity. Adding these miRNAs, we identified conserved miRNAs by aligning unmapped sRNAs to the miRNA sequences in miRBase. Target genes for miRNAs were predicted using the web-based TargetFinder program (http://carringtonlab.org/resources/targetfinder, [87]).

miRNA target validation assays
For miRNA target validation, gene-specific 5 0 RNA ligase-mediated rapid amplification of cDNA ends (5 0 RLM-RACE) was performed using the GeneRacer Kit (Invitrogen). 1.5μg of total RNAs from R. thunb, Marcia, Sympathy and Vital were ligated to 0.25μg of the GeneRacer RNA oligo adapter (5'-CGACUGGAGCAC GAGGACAC UGACAUGGACUGAAGGAGUAGAAA-3'). The combination of oligo(dT) and random hexamers were then used to prime 1 st strand cDNA synthesis in a reverse transcription reaction. The resulting cDNA was PCR-amplified with GeneRacer 5' primer (5'-CGACTG GAGCACGAGGACACTGA-3') and each respective gene-specific primer (shown in Additional file 4). The PCR product was further amplified by nested PCR using GeneRacer 5' nested primer (50-GGACACTGACATG GACTGAAGGAGTA-30) and each respective genespecific primer (shown in Additional file 4). The final PCR product was gel-purified and finally cloned into TA vector (RBC Bioscience) for sequencing.

miRNA profile analysis
Audic's test is statistic analysis which often applied to analyze digital gene expression profiles [88]. We applied this test to estimate miRNA profiles for each Rosa. Read count of each identified miRNA is normalized to the total number of miRNA read counts that are matched to the reference genome or known miRNA in each sample. The expression of miRNAs were represented by the number of reads (or redundancies) of the same miRNA tag (or sequence). We hypothesized that miRNA is likely to be specifically abundant in Rosa, when a miRNA tag has a large number of reads derived from a given Rosa cultivar as compared to another Rosa cultivar. We applied Audic's test to estimate the probability of differential expression for individual miRNA between two pools (a specific Rosa cultivar versus other Rosa cultivars) of reads [88]. Details were denoted in Additional file 8 with the 2 * 2 contingency table and the probability equation for Audic's test [88].