Skip to main content

Characterization and comparative profiling of the small RNA transcriptomes in two phases of flowering in Cymbidium ensifolium

Abstract

Background

Cymbidium ensifolium is one of the most important ornamental flowers in China, with an elegant shape, beautiful appearance, and a fragrant aroma. Its unique flower shape has long attracted scientists. MicroRNAs (miRNAs) are critical regulators in plant development and physiology, including floral development. However, to date, few studies have examined miRNAs in C. ensifolium.

Results

In this study, we employed Solexa technology to sequence four small RNA libraries from two flowering phases to identify miRNAs related to floral development. We identified 48 mature conserved miRNA and 71 precursors. These conserved miRNA belonged to 20 families. We also identified 45 novel miRNA which includes 21 putative novel miRNAs*, and 28 hairpin forming precursors. Two trans-acting small interfering RNAs (ta-siRNAs) were identified, one of which was homologous to TAS3a1. TAS3a1 belongs to the TAS3 family, which has been previously reported to target auxin response factors (ARF) and be involved in plant growth and floral development. Moreover, we built a C. ensifolium transctriptome database to identify genes targeted by miRNA, which resulted in 790 transcriptomic target unigenes. The target unigenes were annotated with information from the non-redundant (Nr), gene ontology database (GO), eukaryotic orthologous groups (KOGs) and Kyoto encyclopedia of genes and genomes (KEGG) database. The unigenes included MADS-box transcription factors targeted by miR156, miR172 and miR5179, and various hormone responding factors targeted by miR159. The MADS-box transcription factors are well known to determine the identity of flower organs and hormone responding factors involved in floral development. In expression analysis, three novel and four conserved miRNA were differentially expressed between two phases of flowering. The results were confirmed by RNA-seq and qRT-PCR. The differential expression of two miRNA, miR160 and miR396, targeted ARFs and growth regulating factor (GRF), respectively. However, most of these small RNA were clustered in the uncharacterized group, which suggests there may be many novel small non-coding RNAs yet to be discovered.

Conclusion

Our study provides a diverse set of miRNAs related to cymbidium floral development and serves as a useful resource for investigating miRNA-mediated regulatory mechanisms of floral development.

Background

Cymbidium ensifolium, is a diploid plant with 40 chromosomes (2n) and an estimated haploid genome size of 4,000 Mb [1]. C. ensifolium, is under subgenus Jensoa of genus Cymbidium in the orchid family (Orchidaceae). C. ensifolium has a long juvenile phase before flowering, generally five to six years, which delays breeding programs [2]. The floral morphology is one of the most important factors determining its ornamental value. Thus, it will be necessary to understand the genetic mechanisms underlying floral development in C. ensifolium to change flowering time and modify floral traits. Although much effort has been devoted to the cloning and identification of key genes involved in floral development and flowering of Cymbidium species [3, 4], the role of non-coding RNAs in C. ensifolium floral development is poorly understood.

In general, small, non-coding RNAs are grouped into two major classes: microRNAs (miRNAs) and short-interfering RNAs (siRNAs). Both small RNA classes have the same chemical composition and mechanism of action. However, siRNAs and miRNAs can be distinguished by their origin, evolutionary conservation and the types of genes that they silence [5]. MiRNAs are endogenous non-coding RNAs of ~21 nucleotides (nt) in length, derived from single-stranded RNA hairpin precursors cleaved by a double-stranded-specific RNase (Dicer) in animals and Dicer-like1 (DCL1) in plants [6]. After DCL1 trims the hairpin precursor, the miRNA/miRNA* duplex is released, and most of miRNA* sequences are quickly degraded [7, 8]. However, some miRNA* occur in a large abundance, such as two conserved miRNA families, miR171 and miR396, which were identified in Rosa hybrida [9] and Brassica napus [10]. There also exist non-conserved or species-specific small interfering RNA in plants, which are the trans-acting small interfering RNAs (ta-siRNAs). Ta-siRNAs are similar to miRNAs and act in trans to regulate messenger RNAs (mRNAs) [11]. Ta-siRNAs were generated from TAS gene transcripts, which are cleaved by a miRNA, resulting in the production of small RNA fragments 21 nt in length that are in phase with the miRNA cleavage site [12].

MiRNAs play an essential role in plant flowering time and floral organ identity. In Arabidopsis, three miRNA families, miR156/ 157, miR159 and miR172, have been shown to be involved in flowering time control under special environmental conditions. The overexpression of miR156/157 can delay flowering significantly [13]. MiR159 can regulate floral transition in short-day photoperiods, by repressing the floral meristem-identity gene [14]. MiR172 can promote flowering by integrating ambient temperature signals into the flowering genetic network [15]. Similarly, in rice and maize, miR156 and miR172 have been shown to control the conversion of spikelet meristems to floral meristems to ensure the initiation of floral organ primordia [16, 17]. Moreover, miR156 and miR172 are also involved in floral organ formation. MiR172 controls the inner whorl organ formation [18], whereas the miR164 gene is involved in establishment of boundaries between components within floral organs, and regulating the resulting floral organs sizes [19]. Recently, there are a few studies on miRNA in other orchid genera, i.e. phalaenopsis [20], orchis [21], and erycina [22] but not in cymbidium.

Three major approaches have been used for identifying miRNAs in plants: forward genetics, bioinformatic prediction and direct cloning and sequencing [23]. Only a few miRNAs have been identified by forward genetic studies [24, 25] and bioinformatic prediction for species-specific miRNAs is difficult, especially for species without a well-defined genome or a reference genome. The development of high-throughput sequencing technologies has greatly improved the success of detecting miRNA’s through direct cloning and sequencing. The Solexa platform, which is capable of yielding a high number of reads up to 35 bp in size [26] is suitable for sequencing miRNA, and identifying low-abundance and tissue-specific miRNA [27].

In our study, we deep sequenced the small RNA transcriptome in C. ensifolium using Solexa sequencing. The short C. ensifolium sequences were used to predict conserved and novel miRNAs. In order to explore the potential genes targeted by miRNA, we also created a large C. ensifolium mRNA transcriptome for comparison. The objective of the present study was to discover miRNAs and their corresponding target genes involved with C. ensifolium floral development.

Methods

Sample collection and preparation

Native cultivars of C. ensifolium “Tiegusu” with light green flowers were grown in a greenhouse under natural light with temperatures between 23 °C and 28 °C (Hangzhou, China). “Tiegusu” is one of the most widely known commercial cultivars in China. Flower buds were collected from two phases of floral development. In phase 1 (PS1) the bud dimension is < =0.5 cm; in phase 4 (PS4) bud dimension is between 2 cm and 3 cm in length (Fig. 1). Each tissue sample consisted of a mixture of five to eight flower buds. The two biological replicates for phase 1 and phase 2 resulted in four libraries.

Fig. 1
figure1

Developmental phases of the flower bud in C. ensifolium. Phase 1 (PS1), bud length < =0.5 cm, is the primary phase; Phase 4 (PS4), with bud length between 2.0 cm to 3.0 cm, is the final phase before blooming

RNA isolation and quality control

Total RNA was isolated using a Trizol kit (Invitrogen, USA). Samples were run on a 1 % agarose gel to check for RNA degradation and contamination. RNA purity was checked using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library construction and small RNA sequencing

A total amount of 3 μg total RNA per sample was used as input material to construct small RNA libraries. Sequencing libraries were generated using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA.) following manufacturer’s recommendations and index codes were added to attribute sequences to each sample. Briefly, NEB 3’ SR adaptors were ligated to the 3' end of small RNA. After the 3' ligation reaction, the SR RT primer was hybridized to the 3' SR adaptor transforming the single-stranded DNA adaptor into a double-stranded DNA molecule. A 5´end adapter was then ligated to the 5´end and first strand cDNA was synthesized using M-MuLV Reverse Transcriptase (RNase H). PCR amplification was performed using LongAmp Taq 2X Master Mix, SR Primer for illumina and an index (X) primer. PCR products were purified on a 8 % polyacrylamide gel (100 V, 80 min). DNA fragments corresponding to 140 ~ 160 bp (the length of small noncoding RNA plus the 3' and 5' adaptors) were recovered and dissolved in 8 μL elution buffer. Library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. All of these processes were performed at the Beijing Novo Gene Genomics Institute, China.

The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500/2000 platform and 50 bp single-end reads were generated. The data were uploaded to NCBI (http://www.ncbi.nlm.nih.gov/) for public use (Accession: SRR2037124).

Data filtering and read mapping

Reads with Poly N’s, 5’ adapter contaminants, low quality, or lacking the 3’ adapter or inset tags were removed. The Q20, Q30, and GC-content of the raw data were calculated and read lengths within a specific size range from 17 to 30 nt were chosen for analysis.

To remove tags originating from protein-coding genes, repeat sequences, ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear ribonucleic acid (snRNA), and small nucleolar RNA (snoRNA), small RNA tags were mapped to RepeatMasker, Rfam, and related databases of the model species, i.e. Oryza and Arabidopsis.

The small RNA tags were mapped to C. ensifolium transctriptome sequence [28] using Bowtie [29] and analyzed for their expression and distribution on the reference genome without mismatches present.

Conserved miRNA alignment

Mapped small RNA tags were used to look for conserved miRNA. miRBase20.0 was used as a reference, modified software mirdeep2 [30] and srna-tools-cli [31] were used to identify potential miRNA and draw the secondary structures. Novo Custom scripts (Beijing Novo Gene Genomics Institute, China) were used to obtain the miRNA counts as well as base bias on the first position of identified miRNA with certain length and on each position of all identified miRNA respectively.

Novel miRNA prediction

The characteristics of hairpin structure of miRNA precursor can be used to predict novel miRNA. The available software miREvo [32] and mirdeep2 [30] were integrated to predict novel miRNA through exploring the secondary structure by identifying the Dicer cleavage site and the minimum free energy of the small RNA tags unannotated in the former steps. Ta-siRNA was identified by comparing with ta-siRNA of Arabidopsis and rice using BLAST, or predicted by the software “UEA sRNA tools” [31]. Analysis of miRNA counts and base bias was performed as conserved miRNA.

Small RNA annotation summary

During the alignment and annotation process, some small RNA tags may be mapped to more than one category. To make every unique small RNA map to only one annotation, we used the following priority rule: conserved miRNA > rRNA > tRNA > snRNA > snoRNA > repeat > novel miRNA > ta-siRNA.

MiRNA family analysis, target gene prediction and annotation

In our analysis pipeline, conserved miRNA used miFam.dat (http://www.mirbase.org/ftp.shtml) to look for families, and novel miRNA precursors were submitted to Rfam (http://rfam.sanger.ac.uk/search/) to identify Rfam families. Target gene prediction of miRNA was performed by psRobot_tar in psRobot for plants [33]. The 5’-end of the cleavage product of C. ensifolium DEF3 (CenDEF3) (JQ326260.1) was determined by a modified 5’-RACE using the RLM-RACE GeneRacerTM kit (Invitrogen, Carlsbad, CA, USA). The amplification was performed according to the manufacturer’s instructions. After reverse transcription, the 5’-end of the cleavage product was amplified using the GeneRace 5’ primer and gene specific reverse primers CenDEF3R 5’-TCCCGCAGTAAGTTCCTGTGGGTTTC-3’. The amplified product was sequenced.

Target genes were compared with protein databases, including the non-redundant database (http://www.ncbi.nlm.nih.gov/), using BLASTX with a significance cut-off E-value of 1e-5. Gene Ontology (GO) enrichment analysis was performed on the target genes of differentially expressed miRNAs using GOseq R package [34]. GO terms with corrected P-value less than 0.05 were considered significantly enriched by target genes. GO terms with significance were summarized into major categories using GOSlim Viewer with the Plant GOSlim set. KEGG pathway analyses of the target genes were performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.ad.jp/kegg/) database. We used KOBAS [35] software to test the statistical enrichment of the target unigenes in KEGG pathways, with the threshold of corrected P value <0.05.

Expression analysis of miRNA

miRNA expression levels were estimated by TPM (transcript per million) using the following criteria [36]. Normalized expression = mapped readcount/Total reads*1000000. Differential expression analysis of two flowering phases was performed using the DESeq R package (1.8.3). MiRNAs that had change ratios of more than 2 or less than 0.5 (Fold change Log2 > 1 or < −1) and P < 0.01 was set as the threshold for significant differential expression by default.

To validate the presence and expression of the identified miRNAs, miRNAs with differential expression and eight random miRNA were selected for stem-loop Real-time PCR (RT-qPCR) as described previously [37]. The total RNA was reverse-transcribed using miRNA specific stem-loop primers (Additional file 1). The reverse transcription reaction was performed with M-MLV (Takara, China), and the reverse-transcribed products were used as the template for RT-qPCR with gene-specific primers (Additional file 1). All reactions were assayed in three biological and technical replications, and performed in an ABI PRISM 7900HT (Applied Biosystems, USA) using Platinum SYBR Green qPCR SuperMix-UDG (Invitrogen, USA). PCR condition consisted of: pre-denaturation and hot start Taq activation at 95 °C for 5 min, then 40 cycles of 95 °C for 15 sec, and 60 °C for 30 sec. The U6 snRNA was used for normalization. The relative expression level of miRNA was calculated according to the method of Livak and Schmittgen [38]. In RT-qPCR, the significant threshold was set at the changes in ratios more than 2 or less than 0.5 (Fold change Log2 > 1 or < −1) and P < 0.01.

Results

sRNA profile in C. ensifolium floral buds

To identify novel miRNAs and investigate the effects of miRNAs during floral development in C. ensifolium, we constructed four small RNA libraries of the floral buds at phase 1 (PS1A, PS1B) and phase 4 (PS4A, PS4B) respectively. Deep sequencing of the small RNA libraries generated the mean number of raw reads, 10,372,041 and 9,105,582, for PS1 and PS4, respectively (Table 1). After removing adapter sequences and discarding low-quality reads, the mean number of clean reads, 10,131,741 and 8,911,915, for PS1 and PS4, respectively, remained for analysis. Finally, the reads ranging from 17 to 30 nt were selected for subsequent analysis. The selection resulted in the mean number of reads, 7,004,064 and 7,231,254, for PS1 and PS4, respectively (Table 2).

Table 1 Summary of small RNA in four libraries from two flowering phases
Table 2 Comparison between small RNA and the Cymbidium transcriptome

The size distribution of small RNAs is summarized in Fig. 2. Overall sRNA ranging from 21 to 24 nt was common, and the 24 nt RNA was the most abundant. The 22-nt sRNAs was the second most abundant in PS1, while the 21 nt sRNAs were the second most abundant in PS4. Of these, almost half of the reads in two PS1 libraries, and one third of the reads in two PS4 libraries could be mapped to the C. ensifolium transcriptome (Table 2). However, a large percentage of small RNAs could not be mapped since the complete C. ensifolium genome sequence is not available.

Fig. 2
figure2

Length distribution of small RNA in four libraries derived from two C. ensifolium flowering phases. A1 and A2: two PS1 libraries from flowering phase I, D1 and D2: two PS1 libraries from flowering phase II. The unigenes are grouped from 17 nt to 30 nt with each column representing the number of unigenes of that specific length

The mapped small RNA sequences were clustered into several RNA classes, such as conserved miRNAs, rRNAs, tRNAs, snRNAs, snoRNAs, repeat and novel miRNA, TAS, and an uncharacterized group (Table 3). Most of these small RNA (89.45 %, group reads/the total reads*100 %) were clustered in the uncharacterized group. The second most abundant group (6.15 %) was rRNA. A search against miRbase20.0 identified 0.41 % conserved miRNAs, and 1.30 % novel miRNAs. The abundance of so many unknown sRNAs suggests that there may be many novel small non-coding RNAs yet to be discovered.

Table 3 The category of all small RNA in four libraries from two flowering phases

Conserved miRNA, their nucleotide bias and family designations in C. ensifolium

To identify conserved miRNAs in C. ensifolium, all mapped small RNA reads in the transcriptome [28] were searched against known miRNAs in the database. The result showed that 47,743 total reads representing 667 unique sRNAs matched to conserved miRNAs (Table 4). In these libraries, two PS1 libraries had the mean of 11,339.00 (161.50 unique) reads, and two PS4 libraries had the mean of 12,532.50 (177.00 unique) reads. Of these unique sRNA, 48 mature miRNAs and 71 precursors were identified. All of precursors are able to adopt hairpin structures resembling the fold-back structure of a miRNA precursor.

Table 4 Summary of the conserved miRNA

Analysis of nucleotide bias was performed on conserved miRNA. The analysis showed that uracil (U) appeared at the first three positions with a high frequency (Fig. 3a). At the 1st, 2nd and 3rd positions, the percentage of U averaged 90.99 %, 57.91 % and 53.73 % respectively. In contrast, guanine (G) or cytosine (C), were seldom present in the first two nucleotides. G and C were the most common nucleotides found in the last four base pair positions of miRNA. C occupied very high percentages, varying from 59.64 % to 90.19 %. By contrast, adenine (A) and uracil (U) were seldom seen near the end.

Fig. 3
figure3

Examination of nucleotide bias within conserved (a) and novel miRNA (b) reads. Height of bar is proportional to the frequency of the corresponding base at the given position from 1st to 22nd

A total of 20 miRNA families were identified among these conserved miRNAs, and have been shown to be conserved across 66 plant species (Additional files 2 and 3). In C. ensifolium, the biggest family was miR166 with 10 members, followed by miR156 with 7 members, and miR171_1 and miR159 both with 6 members. Moreover, miR156/157 was the most popular family, and found in 51 plant species, followed by miR159/319, miR396, and miR166, which are found across 42 plant species. The miR171_1, and miR160 families were identified as well and occur in 40 and 39 different plant species, respectively. The abundance of conserved miRNA varied greatly. The majority of miRNAs had less than 10,000 transcripts per million (TPM) in our study, and some rare miRNAs had less than 100 TPM. Three miRNA, cen-miR159a.1, cen-miR166a-3p and cen-miR166m were very prevalent. Cen-miR159a.1 was the most abundant, followed by cen-miR166a-3p and cen-miR166m, and their TPM averaged 118,849.48, 60,627.64 and 60,627.64, respectively (Additional file 4).

Novel miRNAs and their nucleotide bias in C. ensifolium

In addition to the conserved miRNA, we identified novel miRNAs in C. ensifolium floral buds. 149,612 novel miRNA reads were identified from four flower bud libraries, representing 1,627 unique miRNA (Table 5). The number of reads for these predicted novel miRNAs was more than that observed for conserved miRNAs. These unique sRNAs were derived from 45 novel miRNAs including 21 novel miRNA*. These novel miRNA were involved in 28 hairpin miRNA precursors. Precursors formed characteristic stem-loop structures (Fig. 4). The corresponding miRNA* for these novel families were predicted, providing more evidence that they are indeed novel miRNAs. In general, the miRNA had much more reads than miRNA*. For instance, the novel-miR01 and novel-miR44 were most abundant among novel miRNA families, whereas the corresponding miRNA* were relatively deficient. Interestingly, the novel-miR05* was the most common among miRNA*, while the counterpart was barely detected. The expression of novel miRNA was also tissue specific. For example, two novel miRNA, i.e. novel-miR16 and novel-miR20*, were found only at PS4 (Additional file 4).

Table 5 The summary of the novel miRNA
Fig. 4
figure4

Predicted fold-back structures for the novel miRNAs in C. ensifolium. Mature miRNA sequence is shown in red

Novel miRNA was analyzed for nucleotide bias at each position. Uracil (U) appeared at the beginning of the miRNA with a high frequency (Fig. 3b). At the 1st and 2nd positions, the percentage of U averaged 77.86 % and 91.50 % respectively. Conversely, guanine (G) and cytosine (C) were seldom used as the first two nucleotides. Cytosine frequently occurred near the end, and varied from 62.14 % to 94.18 % from the 18th to 21st position. Conversely, adenine (A) and uracil (U) were seldom at the end. The obvious difference with conserved miRNA was that G and C account for almost half of the nucleotide occurrences at the 3rd and the last position.

Trans-acting siRNAs (ta-siRNAs) are a class of small RNAs that regulate plant development. We searched against TAS database and predicted Ta-siRNAs in C. ensifolium. The search resulted in two ta-siRNA genes, one of which is homologous to TAS3a1.

Target prediction and functional analysis of miRNA in C. ensifolium floral buds

Due to the lack of genomic information available for C. ensifolium, we constructed a transcriptome database for predicting the miRNA targets [28]. A total of 47 conserved and one novel miRNAs targeted 790 potential unigenes, with an average of 16.81 per unigenes/miRNA (Additional file 5). Among the miRNA, cen-miR156k targeted the most unigenes (110), followed by cen-miR396a-3p (88), and cen-miR156a (75). cen-miR160e-5p and cen-miR166b-5p targeted two sequences, while cen-miR5538 targeted only one. Most of miRNA targeted multiple sites, indicating that a single miRNA can regulate more than one unigene. Conversely one unigene was also targeted by several miRNAs. A total of 162 unigenes could be regulated by more than one miRNA, two of which were targeted by 8 conserved miRNA.

The 790 miRNA target unigenes were searched against the non-redundant database, 656 of which had hits that exceeded the E-value threshold (Additional file 5). Among 656 targets, 62 were involved in floral development/flowering time. Examples include squamosa promoter-binding-like protein (SPB) (cen-miR156a, cen-miR529b, and cen-miR156k), floral homeotic protein APETALA 2 (cen-miR172a and cen-miR172c), and AP3/MADS5/DEFICIENS-like MADS-box transcription factor (cen-miR156a and cen-miR5179). The cleavage product of cen-miR5179 on C. ensifolium CenDEF3 transcript were successfully amplified, and miR5179 mediated cleavage was confirmed by modified 5’ RLM-RACE (Fig. 5). Additionally, 17 targets were involved in various hormone signaling pathways, such as auxin response factor, (targeted by cen-miR171a, cen-miR167d-5p, cen-miR167a-5p, cen-miR160e-5p, and cen-miR160a-5p), and ethylene responsive element binding protein 1 (targeted by cen-miR529b). However, 134 targets did not match any homologue in database. A large number of targets were only annotated with vague information, e.g. 95 targets were “unnamed protein”, and 93 were not predicted.

Fig. 5
figure5

Cleavage of cen-miR5179 on CenDEF3 in C. ensifolium. The corresponding region from the position 360 to 380, considering the first nucleotide of the ATG start codon as position (+1).The arrows indicate the position of the cleavage and the number of clones corresponding to each site as deduced by the cloning and sequencing of the obtained fragment

Differential expression of miRNA between two C. ensifolium flowering phases and their target genes

To identify differential expression of miRNAs between two flowering phases, normalized miRNA levels were compared between PS1 and PS4. MiRNAs that had fold change log2 > 1 or < −1, and P < 0.01, were considered to be differentially expressed. Five novel and seven conserved miRNA sequences were differentially expressed. Based on the results of high-throughput sequencing, we performed qRT-PCR on these miRNA. In our qRT-PCR experiment, seven (58.33 %) of these miRNA exhibited significantly different expression (Fig. 6; Additional file 6). Among the seven miRNA, three miRNA were significantly up-regulated at PS4, while the other four miRNA were down-regulated at PS4. The expression patterns were similar to those detected by the RNA-Seq data (Pearson correlation coefficients = 0.948, P < 0.01) (Additional file 7). The results demonstrated the reliability of the mRNA-seq results. Thus, our experiments confirmed that the seven miRNAs were differentially expressed during two flowering phases.

Fig. 6
figure6

Expression levels of four conserved and three novel miRNAs between two flowering phases with significance in qRT-PCR. Based on RNA-seq results, 12 significantly different expressed miRNA were exposed to qRT-PCR assay for validation, which resulted in seven validated miRNA. Each experiment consisted of three biological and technical reps. The relative expression levels of the selected genes were calculated using the 2–ΔΔCT method, and the significant threshold was set at 0.01 (**: <0.01). The U6 gene was used as a Control. Error bars represent the standard deviation of the mean expression values. At PS4, novel-miR11, novel-miR31, novel-miR37 and cen-miR396a-5p were down-regulated, while cen-miR160a-5p, cen-miR394 and cen-miR535-5p were up-regulated

The annotation of the target unigene of differentially expressed miRNAs was conducted based on GO enrichment and KEGG analyses. GO categories were assigned for target unigene, and 50 putative targets were assigned into 527 GO terms (Additional file 8). A total of 164 significant GO terms (Corrected P value < 0.05) were categorized as “biological process”, “cellular component” and “molecular function” (Fig. 7). Within the biological process category, the unigenes were classified into 20 groups, of which the most represented GO terms were related to “multicellular organismal development”, followed by “metabolic process” and “cellular process”. Based on molecular function, the unigenes were classified into three groups, of which the top one was “RNA binding”. Based on cellular components, the unigenes were classified into six groups, of which all were similarly represented. KEGG pathway analysis was performed for biological interpretation of these target unigenes of differentially expressed miRNAs. A total of nine biochemical pathways were involved with the miRNAs target unigenes (Additional file 9). On the threshold of corrected P-value <0.05, two were significantly enriched pathways involved with six unigenes, i.e. “Arginine and proline metabolism” and “Cysteine and methionine metabolism” (Fig. 8).

Fig. 7
figure7

Gene ontology of unigenes targeted by differently expressed miRNAs. The significantly enriched GO terms were categorized according to cellular component, molecular function and biological process

Fig. 8
figure8

The nine enriched KEGG pathways of the target genes. The genes were targeted by differentially expressed miRNAs between floral development phases in C. ensifolium. The enrichment factor is on the x-axis and the y-axis contains the pathway names. The size of each point represents the number of genes enriched in a particular pathway. A larger enrichment factor value and lower Q-values indicates a greater degree of enrichment

Discussion

Distribution and nucleotide bias of sRNA in C. ensifolium floral buds

The sRNAs ranging from 21 to 24 nt were common in this study. These sRNA represent the typical plant mature miRNAs [39]. The majority of 24 nt sRNA is congruent with previous reports on Arabidopsis [40], Oryza [41], Arachis [42], Medicago [39], Populus [43], Gossypium [44, 45], Phalaenopsis [20] and Orchis [21]. However, the size distribution differed from wheat and conifer small RNA obtained through 454 high-throughput sequencing [46, 47] and Chinese yew small RNA obtained through Solexa technology [48]. In conifers, 24 nt miRNAs comprise only a small percentage of sRNA, possibly due to the lack of DCL3, the enzyme that cleaves the precursor into mature 24 nt RNAs in angiosperms [46, 49]. In the current study the main difference between PS1 and PS4 was seen within the 21 and 22 nt sRNAs.

In analysis of nucleotide bias, uridine (U) dominated the first position of the novel and conserved miRNAs in this study. The result was consistent with previous studies [50]. The 5′ terminal nucleotide depends on an miRNAs biogenesis mechanism, for example, AGO1 (a protein from Dicer) harbors miRNAs with a 5’ terminal U [50]. The difference between conserved and novel miRNA, was that G/C accounted for almost 60 % at the 3rd base from the 5’ end in novel miRNA and almost 80 % at the last positions in conserved miRNA. The region from 2 to 8 is the “seed region” of an miRNA, which binds with the target gene(s) during gene regulation [51]. The difference in this region suggested conserved and novel miRNA might target different mRNA or different positions on the same mRNA. miRNAs targeting different positions is correlated with different silencing activity, for example, lower GC content was correlated with higher activity [50].

Specific expression of conserved miRNA in C. ensifolium flower buds

The expression of miRNA families differs across species (Additional files 2 and 3). In C. ensifolium, cen-miR159a.1 was the most abundant, followed by cen-miR166a-3p and cen-miR166m. In previous studies on other plant species tae-miR169b in wheat, osa-miR169 in rice, zma-miR156a in maize [52] and mtr-miR156a in Medicago truncatula [53] are the most frequently observed miRNAs, while miR156 in rice and wheat exhibits low abundance [47]. This may suggest a tissue-specific or species-specific expression profile for miRNAs. The highly expressed miRNAs may be involved in regulating essential functions. For example, miR159 is associated with flowering in short days and male sterility (Palatnik et al. 2007; Achard et al. 2004; Quesada et al. 2005). MiR166 targets genes that regulate shoot apical meristem and floral development [54].

In C. ensifolium, the largest miRNA family size was miR166 which consisted of 8 members, followed by miR171 and miR396 both with 5 members. The miRNA families with a large number of members may be indicative of their function in floral development. Different family members also displayed drastically different expression levels (Additional files 2 and 3). For example, in PS1 members of the miR159 family had TPM ranging from 271.6 (cen-miR159d) to 105939.6 (cen-miR159a.1). This was also the case for some other miRNA families, such as miR166 (0 ~ 30090.4 tpm) and miR396 (23.9 ~ 6146.7 tpm). It is possible that the most abundant member of a miRNA family correlates with its importance for performing key regulatory functions at a particular flowering phase.

Identification of novel miRNAs and ta-siRNA in C. ensifolium flower buds

We identified 24 novel miRNA and 21 novel miRNA* which corresponded to 28 precursor sRNA. Some members were tissue-specific miRNA, and most of these members were miRNA*, such as novel-miR02*, novel-miR10*, novel-miR20*, and novel-miR24*. The expression of these miRNA* was much lower than others which was consistent with previous reports [10, 45, 52]. As miRNA* strands are degraded rapidly during the biogenesis of mature miRNAs, it is possible that the miRNA* strand accumulates at a lower level [10].

The novel C. ensifolium miRNAs displayed much higher expression levels compared to conserved families. Due to the limited genomic information available for C. ensifolium, many of the highly expressed miRNA’s identified in this study are being reported for the first time. These miRNA may represent only a small portion of novel miRNA families involved with C. ensifolium floral development and the sequence information provided will be critical for future research.

Trans-acting siRNAs (ta-siRNAs) are a class of small RNAs that regulate plant development [55]. Ta-siRNAs are miRNA-triggered secondary siRNAs from miRNA-targeted noncoding transcripts. We identified two ta-siRNAs, one ta-siRNA gene was homologous to TAS3a1 and the other is a novel TAS. TAS3 is the center of an autoregulatory network involving miR390, which targets auxin response factors 2 (ARF2), ARF3, and ARF4 [56, 57]. Mi390 assists in cleaving the TAS3 transcript. TAS3 inhibits ARF2/3/4 expression, and ARF4 down regulates miR390 accumulation. In contrast, ARF3 up-regulates miR390 in response to auxin [58]. These target ARF genes, are conserved in angiosperms and play critical roles in plant development [59, 60].

Genes targeted by miRNA and their function during C. ensifolium flowering

Analysis of miRNA targets revealed 47 conserved miRNAs, one new miRNA, and 790 potential targets. Most of these miRNA targeted more than one gene, suggesting these miRNA are involved in several biological processes. Similarly, a single gene was also targeted by several miRNAs. This result could be due to miRNAs mapping to the same cDNA at different sites, and cleaving the mRNA into different-sized fragments [16].

Of these target genes, 62 belonged to MADS-box genes in the expanded ABCDE model of floral development, which determines identity of flower organs [28]. As reported in other plant species, miR156 can act on squamosa promoter-binding (SPB) gene [61, 62]. The SPB promotes flowering by activating miR172 and other MADS box genes [62]. Overexpression of miR156 decreases the SPB’s expression, which results in delayed flowering [17], whereas overexpression of miR172 in Arabidopsis accelerates flowering [24, 63]. Also, miR172 was involved in determining the identity of floral organs by regulating the expression of its target AP2 genes [18, 24]. However, miR156 and miR172 appear to have opposite effects on flowering and regulating floral development in coordination with one another [45, 64]. The conserved mechanism of miR156 and miR172 which mediates transition from the vegetative to the reproductive phase is also reported in orchids Phalaenopsis aphrodite [20], Erycina pusilla [22] and Orchis italica [21]. In our study, cen-miR156 and cen-miR172 were identified during both flowering phases, suggesting that they may have roles in flowering and determining identity of floral organs in C. ensifolium. Another MADS-box gene, DEF-like transcription factor, plays a role in the diversification of tepals and lip. Here, the DEF-like gene 3 was targeted by cen-miR5179, and the cleavage was confirmed by RLM-Race. The similar result is also reported in orchid O. italica [21]. The regulation of miR5179 on DEF-like transcript factor suggests that miR5179 exerts an important function in the diversification of the perianth [21].

MiR159 targeted several MYB genes that act as floral developmental regulators [16, 65]. For example, miR159 degrades its target MYB33 in the process of gibberellin (GA)-induced flowering in Arabidopsis, resulting in non-expression of LEAFY, as well as delayed flowering in short days and male sterility [14, 66, 67]. MYB33 may mediate GA signaling by binding to the promoter of the floral meristem-identity gene, LEAFY [68]. In addition, MYB21, MYB24, and MYB57 are all DELLA repressible GA-response genes that mediate stamen filament growth, and regulate stamen maturation through jasmonate [69].

Xyloglucan endotransglucosylase (XTH) expression was targeted by cen-miR159f. The xyloglucan plays an important structural role in the primary and secondary walls [70]. In plant cell walls, the availability of xyloglucan and its oligosaccharides determine tissue tension and is involved in the control of plant growth and development [71, 72]. In C. ensifolium, XTH may play an important role in petal morphogenesis [73]. In Arabidopsis, XTH9 expression is very high in meristematic tissues [74]; AtXTH28 is specifically involved in the growth of stamen filaments, and is required for successful self-pollination in certain flowers [75].

In addition, 21 target genes were involved in various hormone signaling pathways. As shown in our results, three members of cen-miR167 family were all expressed at low levels during both flowering phases and targeted auxin response factor. In Arabidopsis, miR167 targets two members (ARF6 and ARF8) of auxin response factor (ARF) family, which regulate gynoecium and stamen development in immature flowers [76, 77]. The antagonistic effects of the two hormones, auxin and ethylene, control abscission of flowers [78]. In low concentrations, auxin can delay the senescence of flowers, and the biogenesis of ethylene [79]. Notably, in previous studies on cymbidium, flower abscission is not sensitive to ethylene, although ethylene-sensitive flower abscission is very common in other plant species [80].

Additionally, 134 target genes did not match any homologs in the database, and many of these genes were annotated with only vague information i.e. unnamed protein product, predicted protein etc. This may be due to the limited genomic information available for Cymbidium.

Differential expression of miRNAs during C. ensifolium floral development

MiRNAs that regulate flowering phenology during floral development have been previously reported [18, 81]. In this study, seven miRNA showed significantly different expression between two flowering phases. This was confirmed by RNA-seq and RT-qPCR. The results indicate these miRNA play different roles during floral development. Among the seven miRNA, three of them, cen-miR394, cen-miR160a-5p, and cen-miR535-5p, were significantly up-regulated at PS4, while the other four miRNA, novel-miR37, novel-miR11, novel-miR31 and cen-miR396a-5p, were down-regulated at PS4. In GO analysis, the genes targeted by differently expressed miRNA were significantly enriched within “multicellular organismal development” and “RNA binding” in term of “biological process” and “molecular function”, respectively. The results suggested these target genes could be floral development genes and transcription factors. In previous studies, miR394 and its target, an F-box gene referred to as Leaf Curling Responsiveness (LCR), were reported to be involved in the regulation of leaf curling-related morphological development [82]. In Arabidopsis, miR396 can regulate the expression of Growth Regulating Factor (GRF) genes, which are known to be involved in the control of cell proliferation during leaf and root development [83]. MiR160 negatively regulates ARF10, ARF16, and ARF17 [84, 85]. A loss-of-function mutant floral organs in carpel (foc) exhibits some defects of flowering, such as the formation of irregularly shaped flowers and floral organs inside siliques, as well as reduced fertility [84, 85]. In C. ensifolium, miR396a, miR394, and miR160 may also be involved in floral cell proliferation and morphological development.

Conclusion

This study established a miRNA database for C. ensifolium, which includes the identification of conserved and novel miRNAs. All of these miRNAs were reported for the first time in cymbidium. When the miRNA was compared to the C. ensifolium transcriptome, we investigated a number of transcriptomic unigenes targeted by miRNAs. Some unigenes included MADS-box genes, which are known to determine the identity of flower organs, and were targeted by miR156, miR172 and miR5179. Other unigenes affect various hormones signaling pathways involved in floral development, and were targeted by miR159. Expression analysis of RNA-seq and qRT-PCR identified novel and conserved miRNA that were differentially expressed during flowering phases. They also targeted flowering related genes, such as ARF and growth regulating factor (GRF). These results suggest that these miRNAs may play a crucial role in cymbidium flowering, through regulating the expression level of flowering related genes. Identification of these miRNAs and their targets provides new insight into gene regulation of cymbidium floral development. More novel miRNAs will be identified and their biological function will be revealed as the genomic sequences become available and transformation is established in cymbidium.

References

  1. 1.

    Leitch IJ, Kahandawala I, Suda J, Hanson L, Ingrouille MJ, Chase MW, et al. Genome size diversity in orchids: consequences and evolution. Annals of botany. 2009;104(3):469–81.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  2. 2.

    Li X, Xiang L, Wang Y, Luo J, Wu C, Sun C, et al. Genetic diversity, population structure, pollen morphology and cross-compatibility among Chinese Cymbidiums. Plant Breed. 2014;133(1):145–52.

    CAS  Article  Google Scholar 

  3. 3.

    Wang S, Lee P, Lee Y, Hsiao Y, Chen Y, Pan Z, et al. Duplicated C-Class MADS-Box genes reveal distinct roles in gynostemium development in Cymbidium ensifolium (Orchidaceae). Plant Cell Physiol. 2011;52:563–77.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Xiang L, Li X, Qin D, Guo F, Wu C, Miao L, et al. Functional analysis of FLOWERING LOCUS T orthologs from spring orchid (Cymbidium goeringii Rchb. f.) that regulates the vegetative to reproductive transition. Plant Physiol Biochem. 2012;58:98–105.

    CAS  Article  PubMed  Google Scholar 

  5. 5.

    Bartel B, Bartel DP. MicroRNAs: at the root of plant development. Plant Physiol. 2003;132(2):709–17.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  6. 6.

    Tang G, Reinhart BJ, Bartel DP, Zamore PD. A biochemical framework for RNA silencing in plants. Genes Dev. 2003;17:49–63.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  7. 7.

    Chen X. Small RNAs and their roles in plant development. Annu Rev Cell Dev Biol. 2009;25:21–44.

    Article  PubMed  Google Scholar 

  8. 8.

    Voinnet O. Origin, biogenesis, and activity of plant microRNAs. Cell. 2009;136:669–87.

    CAS  Article  PubMed  Google Scholar 

  9. 9.

    Pei H, Ma N, Chen J, Zheng Y, Tian J, Li J, et al. Integrative analysis of miRNA and mRNA profiles in response to ethylene in rose petals during flower opening. PLoS ONE. 2013;8(5):e64290.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  10. 10.

    Zhao YT, Wang M, Fu SX, Yang WC, Qi CK, Wang XJ. Small RNA profiling in two Brassica napus cultivars identifies microRNAs with oil production- and development-correlated expression and new small RNA classes. Plant Physiol. 2012;158(2):813–23.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  11. 11.

    Vazquez F, Vaucheret H, Rajagopalan R, Lepers C, Gasciolli V, Mallory AC, et al. Endogenous trans-acting siRNAs regulate the accumulation of Arabidopsis mRNAs. Mol Cell. 2004;16(1):69–79.

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    Allen E, Xie Z, Gustafson AM, Carrington JC. microRNAdirected phasing during trans-acting siRNA biogenesis in plants. Cell. 2005;121:207–21.

    CAS  Article  PubMed  Google Scholar 

  13. 13.

    Schwab R, Palatnik JF, Riester M, Schommer C, Schmid M, Weigel D. Specific effects of microRNAs on the plant transcriptome. Dev Cell. 2005;8(4):517–27.

    CAS  Article  PubMed  Google Scholar 

  14. 14.

    Palatnik JF, Wollmann H, Schommer C, Schwab R, Boisbouvier J. Sequence and expression differences underlie functional specialization of Arabidopsis microRNAs miR159 and miR319. Dev Cell. 2007;13:115–25.

    CAS  Article  PubMed  Google Scholar 

  15. 15.

    Jung JH, Seo PJ, Ahn JH, Park CM. Arabidopsis RNA-binding protein FCA regulates microRNA172 processing in thermosensory flowering. J Biol Chem. 2012;287(19):16007–16.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  16. 16.

    Chen Z, Li F, Yang S, Dong Y, Yuan Q, Wang F, et al. Identification and functional analysis of flowering related microRNAs in common wild rice (Oryza rufipogon Griff.). PLoS ONE. 2013;8(12):e82844.

    PubMed Central  Article  PubMed  Google Scholar 

  17. 17.

    Chuck G, Meeley R, Irish E, Sakai H, Hake S. The maize tasselseed4 microRNA controls sex determination and meristem cell fate by targeting Tasselseed6/indeterminate spikelet1. Nat Genet. 2007;39:1517–21.

    CAS  Article  PubMed  Google Scholar 

  18. 18.

    Chen X. A microRNA as a translational repressor of APETALA2 in Arabidopsis flower development. Science. 2004;303(5666):2022–5.

    CAS  Article  PubMed  Google Scholar 

  19. 19.

    Luo Y, Guo Z, Li L. Evolutionary conservation of microRNA regulatory programs in plant flower development. Dev Biol. 2013;380(2):133–44.

    CAS  Article  PubMed  Google Scholar 

  20. 20.

    An FM, Hsiao SR, Chan MT. Sequencing-based approaches reveal low ambient temperature-responsive and tissue-specific microRNAs in phalaenopsis orchid. PLoS ONE. 2011;6(5):e18937.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  21. 21.

    Aceto S, Sica M, De Paolo S, D'Argenio V, Cantiello P, Salvatore F, et al. The analysis of the inflorescence miRNome of the orchid Orchis italica reveals a DEF-like MADS-box gene as a new miRNA target. PLoS ONE. 2014;9(5):e97839.

    PubMed Central  Article  PubMed  Google Scholar 

  22. 22.

    Lin CS, Chen JJ, Huang YT, Hsu CT, Lu HC, Chou ML, et al. Catalog of Erycina pusilla miRNA and categorization of reproductive phase-related miRNAs and their target gene families. Plant Mol Biol. 2013;82(1–2):193–204.

    CAS  Article  PubMed  Google Scholar 

  23. 23.

    Liu B, Li J, Cairns MJ. Identifying miRNAs, targets and functions. Brief Bioinform. 2014;15(1):1–19.

    PubMed Central  Article  PubMed  Google Scholar 

  24. 24.

    Aukerman MJ, Sakai H. Regulation of flowering time and floral organ identity by a MicroRNA and its APETALA2-like target genes. Plant Cell. 2003;15:2730–41.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  25. 25.

    Williams L, Grigg SP, Xie M, Christensen S, Fletcher JC. Regulation of Arabidopsis shoot apical meristem and lateral organ formation by microRNA miR166g and its AtHD-ZIP target genes. Development. 2005;132:3657–68.

    CAS  Article  PubMed  Google Scholar 

  26. 26.

    Szittya G, Moxon S, Santos DM, Jing R, Fevereiro MP. Highthroughput sequencing of Medicago truncatula short RNAs identifies eight new miRNA families. BMC Genomics. 2008;9:593.

    PubMed Central  Article  PubMed  Google Scholar 

  27. 27.

    Song C, Wang C, Zhang C, Korir NK, Yu H, Ma Z. Deep sequencing discovery of novel and conserved microRNAs in trifoliate orange (Citrus trifoliata). BMC Genomics. 2010;11:431.

    PubMed Central  Article  PubMed  Google Scholar 

  28. 28.

    Li X, Luo J, Yan T, Xiang L, Jin F, Qin D, et al. Deep sequencing-based analysis of the Cymbidium ensifolium floral transcriptome. PLoS ONE. 2013;8(12):e85480.

    PubMed Central  Article  PubMed  Google Scholar 

  29. 29.

    Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.

    PubMed Central  Article  PubMed  Google Scholar 

  30. 30.

    Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52.

    PubMed Central  Article  PubMed  Google Scholar 

  31. 31.

    Moxon S, Schwach F, Dalmay T, Maclean D, Studholme DJ, Moulton V. A toolkit for analysing large-scale plant small RNA datasets. Bioinformatics. 2008;24(19):2252–3.

    CAS  Article  PubMed  Google Scholar 

  32. 32.

    Wen M, Shen Y, Shi S, Tang T. miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics. 2012;13:140.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  33. 33.

    Wu HJ, Ma YK, Chen T, Wang M, Wang XJ. PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res. 2012;40(Web Server issue):W22–28.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  34. 34.

    Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.

    PubMed Central  Article  PubMed  Google Scholar 

  35. 35.

    Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93.

    CAS  Article  PubMed  Google Scholar 

  36. 36.

    Zhou L, Chen J, Li Z, Li X, Hu X, Huang Y, et al. Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS ONE. 2010;5(12):e15224.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  37. 37.

    Chen C, Ridzon DA, Broomer AJ, Zhou Z, Lee DH, Nguyen JT, et al. Real-time quantification of microRNAs by stem-loop RT-PCR. Nucleic Acids Res. 2005;33(20):e179.

    PubMed Central  Article  PubMed  Google Scholar 

  38. 38.

    Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) Method. Methods. 2001;25:402–8.

    CAS  Article  PubMed  Google Scholar 

  39. 39.

    Wang T, Chen L, Zhao M, Tian Q, Zhang WH. Identification of drought-responsive microRNAs in Medicago truncatula bygenome-wide high-throughput sequencing. BMC Genomics. 2011;12:367.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  40. 40.

    Hsieh LC, Lin SI, Shih AC, Chen JW, Lin WY, Tseng CY, et al. Uncovering small RNA-mediated responses to phosphate deficiency in Arabidopsis by deep sequencing. Plant Physiol. 2009;151:2120–32.

    PubMed Central  Article  PubMed  Google Scholar 

  41. 41.

    Wei LQ, Yan LF, Wang T. Deep sequencing on genome-wide scale reveals the unique composition and expression patterns of microRNAs in developing pollen of Oryza sativa. Genome Biol. 2011;12(6):R53.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  42. 42.

    Zhao CZ, Xia H, Frazier T, Yao YY, Bi YP, Li AQ, et al. Deep sequencing identifies novel and conserved microRNAs in peanuts (Arachis hypogaea L.). BMC Plant Biol. 2010;10:3.

    PubMed Central  Article  PubMed  Google Scholar 

  43. 43.

    Klevebring D, Street NR, Fahlgren N, Kasschau KD, Carrington JC, Lundeberg J, et al. Genome-wide profiling of Populus small RNAs. BMC Genomics. 2009;10:620.

    PubMed Central  Article  PubMed  Google Scholar 

  44. 44.

    Yang X, Wang L, Yuan D, Lindsey K, Zhang X. Small RNA and degradome sequencing reveal complex miRNA regulation during cotton somatic embryogenesis. J Exp Bot. 2013;64(6):1521–36.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  45. 45.

    Wang ZJ, Huang JQ, Huang YJ, Li Z, Zheng BS. Discovery and profiling of novel and conserved microRNAs during flower development in Carya cathayensis via deep sequencing. Planta. 2012;236:2.

    Google Scholar 

  46. 46.

    Morin RD, Aksay G, Dolgosheina E, Ebhardt HA, Magrini V, Mardis ER, et al. Comparative analysis of the small RNA transcriptomes of Pinus contorta and Oryza sativa. Genome Res. 2008;18(4):571–84.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  47. 47.

    Yao YY, Guo GG, Ni ZF, Sunkar R, Du JK, Zhu JK, et al. Cloning and characterization of microRNAs from wheat (Triticum aestivum L.). Genome Biol. 2007;8:6.

    Article  Google Scholar 

  48. 48.

    Qiu D, Pan X, Wilson IW, Ketchum REB, Li F, Liu M, et al. High throughput sequencing technology reveals that the taxoid elicitor methyl jasmonate regulates microRNA expression in Chinese yew (Taxus chinensis). Gene. 2009;436(1–2):37–44.

    CAS  Article  PubMed  Google Scholar 

  49. 49.

    Dolgosheina EV, Morin RD, Aksay G, Sahinalp SC, Magrini V, Mardis ER, et al. Conifers have a unique small RNA silencing signature. RNA. 2008;14(8):1508–15.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  50. 50.

    Mi S, Cai T, Hu Y, Chen Y, Hodges E, Ni F, et al. Sorting of small RNAs into Arabidopsis argonaute complexes is directed by the 5' terminal nucleotide. Cell. 2008;133:116–27.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  51. 51.

    Zhang B, Stellwag EJ, Pan X. Large-scale genome analysis reveals unique features of microRNAs. Gene. 2009;443(1–2):100–9.

    CAS  Article  PubMed  Google Scholar 

  52. 52.

    Wang L, Liu H, Li D, Chen H. Identification and characterization of maize microRNAs involved in the very early stage of seed germination. BMC Genomics. 2011;12:154.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  53. 53.

    Reinhart BJ, Weinstein EG, Rhoades MW, Bartel B, Bartel DP. MicroRNAs in plants. Genes Dev. 2002;16(13):1616–26.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  54. 54.

    Jung JH, Park CM. MIR166/165 genes exhibit dynamic expression patterns in regulating shoot apical meristem and floral development in Arabidopsis. Planta. 2007;225:1327–8.

    CAS  Article  PubMed  Google Scholar 

  55. 55.

    Chen HM, Li YH, Wu SH. Bioinformatic prediction and experimental validation of a microRNA-directed tandem trans-acting siRNA cascade in Arabidopsis. Proc Natl Acad Sci U S A. 2007;104(9):3318–23.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  56. 56.

    Yoon EK, Yang JH, Lim J, Kim SH, Kim SK, Lee WS. Auxin regulation of the microRNA390-dependent transacting small interfering RNA pathway in Arabidopsis lateral root development. Nucleic Acids Res. 2010;38(4):1382–91.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  57. 57.

    Marin E, Jouannet V, Herz A, Lokerse AS, Weijers D, Vaucheret H, et al. miR390, Arabidopsis TAS3 tasiRNAs, and their auxin response factor targets define an autoregulatory network quantitatively regulating lateral root growth. Plant Cell. 2010;22:1104–17.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  58. 58.

    Luo Q, Mittal A, Jia F, Christopher DR. An autoregulatory feedback loop involving PAP1 and TAS4 in response to sugars in Arabidopsis. Plant Mol Biol. 2012;80:117–29.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  59. 59.

    Felippes FF, Weigel D. Triggering the formation of tasiRNAs in Arabidopsis thaliana: the role of microRNA miR173. EMBO Rep. 2009;10(3):264–70.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  60. 60.

    Xie Z, Allen E, Wilken A, Carrington JC. DICER-LIKE 4 functions in trans-acting small interfering RNA biogenesis and vegetative phase change in Arabidopsis thaliana. Proc Natl Acad Sci USA. 2005;102:12984–9.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  61. 61.

    Lu S, Sun YH, Shi R, Clark C, Li L, Chiang VL. Novel and mechanical stress-responsive microRNAs in Populus trichocarpa that are absent from Arabidopsis. Plant Cell. 2005;17:2186–203.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  62. 62.

    Yu S, Galvão VC, Zhang YC, Horrer D, Zhang TQ, Hao YH, et al. Gibberellin regulates the Arabidopsis floral transition through miR156-targeted SQUAMOSA promoter binding-like transcription factors. Plant Cell. 2012;24(8):3320–32.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  63. 63.

    Jung JH, Seo YH, Seo PJ, Reyes JL, Yun J, Chua NH, et al. The GIGANTEA-regulated microRNA172 mediates photoperiodic flowering independent of CONSTANS in Arabidopsis. Plant Cell. 2007;19:2736–48.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  64. 64.

    Wu G, Park MY, Conway SR, Wang JW, Weigel D, Poethig RS. The sequential action of miR156 and miR172 regulates developmental timing in Arabidopsis. Cell. 2009;138:750–9.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  65. 65.

    Peng J. Gibberellin and jasmonate crosstalk during stamen development. J Integr Plant Biol. 2009;51(12):1064–70.

    CAS  Article  PubMed  Google Scholar 

  66. 66.

    Achard P, Herr A, Baulcombe DC, Harberd NP. Modulation of floral development by a gibberellin-regulated microRNA. Development. 2004;131:3357–65.

    CAS  Article  PubMed  Google Scholar 

  67. 67.

    Quesada V, Dean C, Simpson GG. Regulated RNA processing in the control of Arabidopsis flowering. Int J Dev Biol. 2005;49:773–80.

    CAS  Article  PubMed  Google Scholar 

  68. 68.

    Gocal GF, Sheldon CC, Gubler F, Moritz T, Bagnall DJ, MacMillan CP, et al. GAMYB-like genes, flowering, and gibberellin signaling in Arabidopsis. Plant Physiol. 2001;127(4):1682–93.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  69. 69.

    Cheng H, Song S, Xiao L, Soo HM, Cheng Z, Xie D, et al. Gibberellin acts through jasmonate to control the expression of MYB21, MYB24, and MYB57 to promote stamen filament growth in Arabidopsis. PLoS Genet. 2009;5:e1000440.

    PubMed Central  Article  PubMed  Google Scholar 

  70. 70.

    Saladie M, Rose JKC, Cosgrove DJ, Catala C. Characterization of a new xyloglucan endotransglucosylase/hydrolase (XTH) from ripening tomato fruit and implications for the diverse modes of enzymic action. The Plant Journal. 2006;47(2):282–95.

    CAS  Article  PubMed  Google Scholar 

  71. 71.

    Takeda T, Furuta Y, Awano T, Mizuno K, Mitsuishi Y, Hayashi T. Suppression and acceleration of cell elongation by integration of xyloglucans in pea stem segments. Proc Natl Acad Sci U S A. 2002;99(13):9055–60.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  72. 72.

    Alonso-Simon A, Neumetzler L, Garcia-Angulo P, Encina AE, Acebes JL, Alvarez JM, et al. Plasticity of xyloglucan composition in bean (Phaseolus vulgaris)-cultured cells during habituation and dehabituation to lethal concentrations of dichlobenil. Mol Plant. 2010;3(3):603–9.

    CAS  Article  PubMed  Google Scholar 

  73. 73.

    Li X, Xu W, Chowdhury M, Jin F. Comparative Proteomic Analysis of Labellum and Inner Lateral Petals in Cymbidium ensifolium Flowers. INT J MOL SCI. 2014;15:19877–97.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  74. 74.

    Hyodo H, Yamakawa S, Takeda Y, Tsuduki M, Yokota A, Nishitani K, et al. Active gene expression of a xyloglucan endotransglucosylase/hydrolase gene, XTH9, in inflorescence apices is related to cell elongation in Arabidopsis thaliana. Plant Mol Biol. 2003;52(2):473–82.

    CAS  Article  PubMed  Google Scholar 

  75. 75.

    Harada T, Torii Y, Morita S, Onodera R, Hara Y, Yokoyama R, et al. Cloning, characterization, and expression of xyloglucan endotransglucosylase/hydrolase and expansin genes associated with petal growth and development during carnation flower opening. J Exp Bot. 2011;62(2):815–23.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  76. 76.

    Ru P, Xu L, Ma H, Huang H. Plant fertility defects induced by the enhanced expression of microRNA167. Cell Res. 2006;16:457–65.

    CAS  Article  PubMed  Google Scholar 

  77. 77.

    Wu MF, Tian Q, Reed JW. Arabidopsis microRNA167 controls patterns of ARF6 and ARF8 expression, and regulates both female and male reproduction. Development. 2006;133:4211–8.

    CAS  Article  PubMed  Google Scholar 

  78. 78.

    Stepanova AN, Yun J, Likhacheva AV, Alonso JM. Multilevel interactions between ethylene and auxin in Arabidopsis roots. Plant Cell. 2007;19(7):2169–85.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  79. 79.

    McSteen P, Malcomber S, Skirpan A, Lunde C, Wu X, Kellogg E, et al. Barren inflorescence2 Encodes a co-ortholog of the PINOID serine/threonine kinase and is required for organogenesis during inflorescence and vegetative development in maize. Plant Physiol. 2007;144(2):1000–11.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  80. 80.

    van Doorn WG. Effect of ethylene on flower abscission: a survey. Annals of botany. 2002;89(6):689–93.

    Article  PubMed  Google Scholar 

  81. 81.

    Wang JW, Czech B, Weigel D. miR156-regulated SPL transcription factors define an endogenous flowering pathway in Arabidopsis thaliana. Cell. 2009;138(4):738–49.

    CAS  Article  PubMed  Google Scholar 

  82. 82.

    Song JB, Huang SQ, Dalmay T, Yang ZM: Regulation of leaf morphology by microrna394 and its target LEAF CURLING RESPONSIVENESS. Plant Cell Physiol 2012.

  83. 83.

    Hewezi T, Maier TR, Nettleton D, Baum TJ. The Arabidopsis microRNA396-GRF1/GRF3 regulatory module acts as a developmental regulator in the reprogramming of root cells during cyst nematode infection. Plant Physiol. 2012;159(1):321–35.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  84. 84.

    Mallory AC, Bartel DP, Bartel B. MicroRNA-directed regulation of Arabidopsis AUXIN RESPONSE FACTOR 17 is essential for proper development and modulates expression of early auxin response genes. Plant Cell. 2005;17:1360–75.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  85. 85.

    Liu X, Huang J, Wang Y, Khanna K, Xie Z, Owen HA, et al. The role of floral organs in carpels, an Arabdopsis loss-of-function mutation in micro-RNA160a, in organogenesis and the mechanism regulating its expression. Plant J. 2010;62:416–28.

    Article  PubMed  Google Scholar 

Download references

Acknowledgments

This research was supported by the National Basic Research Program funded by the Nature Science Foundation of China (No.31201648), the Postdoctoral Science Foundation of China (No. 2012 M521203), the Special Postdoctoral Science Foundation of China (No. 2013 T60607), and the Foundation for Selected Postdoctoral project of Zhejiang (Bsh1201032), the Qianjiang talents project (No. 2013R10081).

Author information

Affiliations

Authors

Corresponding authors

Correspondence to Xiaobai Li or Dianxing Wu.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

XL and DW conceived of the study, XL, FJ, LJ, and AJ participated in the experiment design and coordination, and manuscript draft. XL, FJ, LJ, and XM performed the sequence alignment and annotation. AJ, XS and XM carried out the statistical analysis. All authors read and approved the final manuscript.

Additional files

Additional file 1:

Real-time quantitative PCR primer sequences for candidate miRNAs.

Additional file 2:

The conserved miRNA family among all of the species.

Additional file 3:

The graphic display of conserved miRNA family among plant species.

Additional file 4:

The expression of the novel miRNA in four libraries from two flowering phases.

Additional file 5:

Annotation of unigenes targeted by miRNA.

Additional file 6:

The comparison of seven differentially expressed miRNA between RT-qPCR and RNA-seq.

Additional file 7:

Coefficient analysis of fold change data between RT-qPCR and RNA-seq. Based on RNA-seq results, seven significant differentially expressed miRNA were subject to RT-qPCR assay. Data indicating relative transcript level from RT-qPCR are means of three replicates, and data from RNA-seq are means of two replicates. Scatterplots were generated by the log2 expression ratios from qPCR Log2(PS4/PS1) (x-axis) and RNA-seq Log2(PS4/PS1) (y-axis). Pearson correlation coefficient (r) was 0.948 (P < 0.01).

Additional file 8:

Gene ontology (GO) enrichment of the target genes. The genes were targeted by differently expressed miRNA between two floral development phases in C. ensifolium.

Additional file 9:

The KEGG pathway enrichment of the target genes. The genes were targeted by differently expressed miRNA between two floral development phases in C. ensifolium.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Li, X., Jin, F., Jin, L. et al. Characterization and comparative profiling of the small RNA transcriptomes in two phases of flowering in Cymbidium ensifolium . BMC Genomics 16, 622 (2015). https://doi.org/10.1186/s12864-015-1764-1

Download citation

Keywords

  • Gene Ontology
  • Floral Organ
  • miRNA Family
  • Floral Development
  • Auxin Response Factor