Figures
Abstract
Plant microRNAs (miRNAs) are small, regulatory non-coding RNAs involved in a wide range of biological processes, from organ development to response to stimuli. In recent years, an increasing number of studies on model plant species have highlighted the evolutionary conservation of a high number of miRNA families and the existence of taxon-specific ones. However, few studies have examined miRNAs in non-model species such as orchids, which are characterized by highly diversified floral structures and pollination strategies. Therefore, we analysed a small RNA library of inflorescence tissue of the Mediterranean orchid Orchis italica to increase the knowledge on miRNAs in a non-model plant species. The high-throughput sequencing and analysis of a small RNA library of inflorescence of O. italica revealed 23 conserved and 161 putative novel miRNA families. Among the putative miRNA targets, experimental validation demonstrated that a DEF-like MADS-box transcript is cleaved by the homolog of miR5179 of O. italica. The presence of conserved miRNA families in the inflorescence of O. italica indicates that the basic developmental flower regulatory mechanisms mediated by miRNAs are maintained through evolution. Because, according to the “orchid code” theory, DEF-like genes exert a key function in the diversification of tepals and lip, the cleavage-mediated inhibitory activity of miR5179 on a OitaDEF-like transcript suggests that, in orchids, miRNAs play an important role in the diversification of the perianth organs.
Citation: Aceto S, Sica M, De Paolo S, D'Argenio V, Cantiello P, Salvatore F, et al. (2014) 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 9(5): e97839. https://doi.org/10.1371/journal.pone.0097839
Editor: M. Lucrecia Alvarez, TGen, United States of America
Received: November 26, 2013; Accepted: April 25, 2014; Published: May 15, 2014
Copyright: © 2014 Aceto et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This work was supported by the 2009 Regione Campania Grant L.R. N5/2002. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
MicroRNAs (miRNAs) are small (18–25 nt) non-coding RNAs that regulate gene expression at post-transcriptional level. They originate from single-stranded precursors (pre-miRNAs) able to self-pair and form hairpin structures. In plants, the Dicer1-like RNAse III enzyme excises the miRNA/miRNA* duplex from the pre-miRNA. Subsequently, the guide strand of the duplex is recruited by an RNA-induced silencing complex (RISC), whose core component is the protein ARGONAUTE1 (AGO1). This complex permits the interaction between the miRNA and its target mRNA, thereby regulating its expression through mRNA cleavage, translational repression or epigenetic modifications [1], [2]. Plant miRNAs play crucial roles in different processes, from leaf and root morphogenesis to floral induction, organ formation, reproduction and stress response [2], [3], [4], [5].
Advances in next-generation sequencing techniques have prompted a plethora of studies on miRNAs in both model and non-model species and have led to the development of specific in silico analysis tools [6], [7], [8], [9]. Recently, the high-throughput sequencing approach has been applied to the study of miRNAs in two orchid species, Phalaenopsis aphrodite and Erycina pusilla, both belonging to the sub-family Epidendroideae (Orchidaceae) [10], [11]. Conserved and novel miRNAs have been identified in both species, and some differences between the two species have been detected.
The Orchidaceae family is one of the largest among flowering plants and is characterized by highly specialized reproductive strategies and extremely diversified flowers. Despite their remarkable morphological diversity, orchid flowers share bilateral symmetry and the organization into three outer tepals, two inner lateral tepals and an inner median tepal (labellum or lip). In addition, the column, the inner whorl, is constituted by male and female reproductive tissues fused together. At the base of the column is located the ovary, whose maturation is triggered by pollination. The pollen grains (pollinia) are located at the top of the column.
The (A)BCDE model describes the function of different functional classes of homeotic genes involved in flower development [12], [13], [14]. All the genes involved in the (A)BCDE model, except the class A gene APETALA2 (AP2), are MADS-box transcriptional factors. Although this model is applicable to most flowering plants, in orchids it shows modifications in the expression domains of the class B GLOBOSA-like (GLO-like) and DEFICIENS-like (DEF-like) genes [15], [16], [17], [18] and of the class C (AGAMOUS-like) and D (SEEDSTICK-like) genes [19], [20]. To date, the only gene belonging to the (A)BCDE model and negatively regulated by a miRNA is AP2, whose transcript is cleaved by miR172 with a conserved mechanism present also in orchids [10], [21].
The aim of this study is to characterize the miRNA population expressed in the floral buds of Orchis italica, a wild Mediterranean orchid species belonging to the sub-family Orchidoideae. In the attempt to shed light on conserved and novel miRNAs expressed in a non-model plant species, we performed a computational analysis of the floral miRNAs of O. italica based on deep sequencing data in order to identify possible new miRNA targets involved in floral development.
Materials and Methods
Deep sequencing of small RNAs of O. italica
Total RNA was extracted from inflorescence of O. italica before anthesis (∼9 mm diameter size) using the Trizol Reagent (Ambion). Ten different florets from a single inflorescence were used as starting material. Although not synchronous, the selected florets were collected from the bottom of the inflorescence, displayed approximately the same size and could be considered in the same developmental stage. After DNase treatment, RNA was quantified using the Nanodrop 2000 c spectrophotometer (ThermoScientific). The quality of the extracted RNA was checked using the Agilent 2100 BioAnalyzer (Agilent). Small RNA library preparation and sequencing were carried out according to the manufacturer instruction (Illumina). The library obtained was sequenced using the MiSeq instrument (Illumina). The raw reads were deposited at NCBI under the accession PRJNA222705.
In silico analysis of the library
Raw reads of inflorescence tissue of O. italica were processed using the plant version of the UEA sRNA workbench [22], [23] to remove the adaptor sequence, the low quality reads and the reads with abundance lower than 5. Filtering options were set to exclude reads shorter than 18 and longer than 35 nucleotides. After removal of tRNA and rRNA contaminants, reads were collapsed and the read count for identical sequences was summed.
Using the filtered dataset of O. italica as query, standalone BLAST was conducted against the known plant mature and hairpin miRNA sequences downloaded from mirBase 20 [24] to search for conserved miRNAs in the inflorescence of O. italica. In addition, short reads were blasted against the miRNAs expressed in the floral bud of the orchid P. aphrodite [10]. Only reads matching at least 18 nt and with less than 3 mismatches were considered positive.
We used the miRDeep-P software [25], which is specific for plant miRNA analysis, to search for new small RNA sequences in the short reads of inflorescence tissue of O. italica and matching the criteria of plant putative miRNAs. The first step of this analysis is based on the bowtie alignment of the cleaned and collapsed short reads to a reference transcriptome. Unfortunately, the O. italica transcriptome is not available and only a few orchid transcriptomes [10], [26], [27], [28] are deposited in public databases as Sequence Read Archive files. As more data are available for the orchid P. aphrodite than for the other orchids (including tissue-specific sRNA libraries), we used this species as reference although it belongs to the distantly-related sub-family Epidendroideae. We assembled the transcriptome of P. aphrodite starting from the deep sequencing files deposited in the Sequence Read Archive under accession code SRA030409. After adapter removal, filtering and collapsing, reads of vegetative, seed and inflorescence tissue of P. aphrodite were separately assembled using the Trinity software [29]. The three assembled tissue-specific transcriptomes were compared to obtain a not-redundant collection of unigenes that was used as reference transcriptome. The annotation of the assembled transcriptome was performed using the FastAnnotator online tool [30]. Under miRDeep-P, the bowtie alignment was performed allowing for a maximum of three mismatches and reads mapping to multiple positions (maximum 15) were retained. The potential miRNA precursors were then selected by setting the maximum length to 250 nt, and their secondary structure was predicted using RNAfold.
Putative target analysis
We used the psRNATarget online tool [31] for the prediction of the putative targets of the short reads matching the plant miRNA criteria. The cleaned short reads of O. italica were analyzed against the newly assembled transcriptome of P. aphrodite setting the search parameters to the default value (maximum expectation 3.0).
Using more stringent search parameters (maximum expectation 0.0), the cleaned reads were analyzed also against different transcripts of genes involved in flower development of O. italica that we have previously isolated (OrcLFY [GenBank: AB088851], OitaAP2 [GenBank: KF152921], OrcPI [GenBank: AB094985], OrcPI2 [GenBank: AB537504], OitaAG [GenBank: JX205496], OitaSTK [GenBank: JX205497]). Within these transcripts, the class B MADS-box genes are represented by OrcPI and OrcPI2, both of which are members of the GLO-like lineage. We decided to isolate also the other class B genes of O. italica belonging to the DEF-like lineage and to verify if they are targets of miRNAs.
One microgram of total RNA extracted from floral buds of O. italica was reverse transcribed using the Advantage RT-PCR kit (Clontech) and an oligo dT primer. The MADS-box degenerate primer MADS_F (Table 1) and a poly-T primer were used to amplify 1 µl of first strand cDNA using the LongAmp Taq PCR Kit (New England Biolabs) [19], [32]. The amplification products were cloned into the pGEM-T Easy vector (Promega); more than 50 clones were sequenced using the plasmid primers T7 and SP6 and sequencing reactions were run on an ABI 310 Automated Sequencer (Applied Biosystems).
BLAST search of the obtained sequences revealed four different DEF-like cDNAs that were deposited as OitaDEF1-4 [DDBJ: AB857726, AB857727, AB857728, AB857729, respectively]. Virtual translation of the coding sequence of the OitaDEF-like genes was aligned to other DEF-like and GLO-like (class B MADS-box proteins) sequences present in GenBank using ClustalW. The Neighbor Joining (NJ) tree of the resulting alignment was constructed using the MEGA5 software [33], performing 1000 bootstrap replicates.
Based on the predicted position of the putative miRNA target site on the OitaDEF-like cDNAs, specific reverse primers were designed and used to verify the presence of the cleavage product performing a modified 5′-RACE experiment [34]. All the DEF-like specific reverse primers annealed downstream of the predicted putative miRNA cleavage site. Using the RLM-RACE GeneRace kit (Invitrogen), the 5′ adaptor was ligated to the 5′-terminus of the RNA extracted from inflorescence tissue of O. italica (500 ng) without any enzymatic treatment to remove the 5′ cap. After reverse transcription, cDNA was amplified using the GeneRace 5′ Primer and the OitaDEF1R1, OitaDEF2R1, OitaDEF3R1 or the OitaDEF4R1 specific reverse primers (Table 1). A second PCR reaction was conducted on 1 µl of the first reaction using the nested adaptor forward primer and the nested OitaDEF1R2, OitaDEF2R2, OitaDEF3R2 or the OitaDEF4R2 specific reverse primers (Table 1). The amplification products were cloned and sequenced as described above.
Tissue expression analysis
Using the protocol described above, we extracted total RNA from leaf, outer and inner tepal, lip and column dissected from the inflorescence of O. italica before anthesis. In this bud stage, the flower organs are already formed and cell division is completed, while cell elongation is still ongoing. We extracted total RNA from ovary tissue before and after manual fertilization (3, 7 and 10 days after pollination). Within the immature ovules of O. italica (before pollination) the megaspore mother cell is undergoing the first meiotic division. At 3 days after pollination, the female gametophyte has completed its development; at 7 days after pollination, fertilization has occurred and the seeds are in the early developmental stage; at 10 days after pollination the seeds are mature (Barone Lumaga, personal communication).
Within the miRNAs in the inflorescence of O. italica, we selected 10 of those most highly expressed and/or with interesting putative targets (e.g. involved in miRNA biogenesis, in flower development, etc), and analyzed their relative expression profile in the different tissues. Eight of the selected miRNAs correspond to known plant miRNAs and two fulfil the plant miRNA criteria but do not correspond to any known miRNAs. The primer sequences of the selected miRNAs are listed in Table 1.
The Poly(T) Adaptor RT-PCR method [35] was used to analyse the expression pattern of the selected miRNAs starting from 350 ng of total RNA from each tissue. A poly-T adaptor was ligated to the 3′-terminus of the miRNAs and subsequently reverse transcribed. The forward primer specific for each selected miRNA was used in combination with the poly-T adaptor reverse primer during the Real Time PCR experiments. The latter were performed on 30 ng of the first strand cDNA from each tissue using the 5.8S RNA as endogenous control gene [35] and the SYBR Green PCR Master Mix (Life Technologies). Reactions were conducted in technical and biological triplicates. The Real Time PCR Miner online tool [36] was used to calculate PCR efficiency (E) and optimal threshold cycle (CT) for each well. The mean relative expression ratio (Rn) and standard deviation of the selected miRNAs in the different tissues was calculated using the 5.8S RNA as the endogenous control gene with the following formula Rn = (1+EmiRNA)-CTmiRNA/(1+E5.8S)-CT 5.8S. The Rn values were multiplied by an arbitrary value (100,000). Differences in relative expression levels of the 10 selected miRNAs in the various tissues were assessed by ANOVA followed by the Tukey HSD post hoc test. Real Time PCR product of several samples was cloned and sequenced to exclude the presence of amplification artefacts.
Real Time-PCR experiments were conducted to analyze the relative expression of the OitaDEF-like genes of O. italica in the same tissues. Total RNA (350 ng) from each tissue was reverse transcribed as described above. The Real Time PCR reactions were conducted using 30 ng of the first strand cDNA from each tissue and the following primer pairs: OitaDEF1F/OitaDEF1R1, OitaDEF2F/OitaDEF2R1, OitaDEF3F/OitaDEF3R1, OitaDEF4F/OitaDEF4R1 (Table 1). The actin OitaAct gene was used as endogenous control [GenBank: AB630020]. Amplification conditions are reported elsewhere [16]. The amplification data were processed as described above.
Results and Discussion
Conserved and novel miRNAs in the inflorescence of O. italica
The high-throughput sequencing of the small RNA library of inflorescence tissue of O. italica produced 4,718,127 total reads. Adaptor trimming and length filtering (18–35 nucleotides) reduced the total reads to 3,019,126. Short sequences representing less than 5 reads, and tRNA and rRNA contaminant sequences were removed thereby resulting in 1,064,237 total and 37,818 distinct reads (Table 2, File S1). Within the specified length range, the number of reads was highest between 21–24 nucleotides, with a peak at 24 nt (Figure 1). This is in agreement with the distribution pattern of small RNAs in other plant species such as Arabidopsis thaliana [37], Caragana intermedia [8], Solanum tuberosum [38], Avicennia marina [7], Oryza sativa [39] and the orchid P. aphrodite [10]. BLAST analysis revealed 175 distinct sequences in the inflorescence tissue that corresponded to 23 known plant miRNAs and 6 star sequences (Figure 2A, Table 3, File S2). Among the former, the most abundant are miR6300, miR166, miR319, miR396, miR393 and miR168, which are also the most heterogeneous, being represented by the highest number of members (Figure 2B).
A) Sum of reads; B) number of members.
The comparison of the conserved miRNAs between the inflorescence tissue of O. italica and P. aphrodite (File S3) revealed that almost all miRNAs are shared in both species, albeit with some differences in their relative abundance. For example, miR166 and miR168 are among the most abundant miRNAs in the floral buds of both O. italica and P. aphrodite, miR319 and miR396 are expressed at a higher level in O. italica than in Phalaenopsis, miR528 is highly expressed in Phalaenopsis whereas it has only five reads in O. italica. In addition, miR6300 is highly expressed in O. italica whereas it is absent from Phalaenopsis. These differences could reflect the different developmental stage of the flower bud used to construct the small RNA libraries in the two orchid species and the size difference of the two libraries [10].
The known miRNAs identified in the inflorescence of O. italica correspond to 16 of the 37 evolutionary and taxonomically conserved miRNA families [40] (Figure 3). Given that 5 of these conserved miRNAs are specific of dicots and 1 is specific of Bryophyta and Lycopodiophyta, approximately 52% (16/31) of the evolutionary conserved miRNAs is expressed in the inflorescence of O. italica. Among them, one (miR528) is specific of monocots. In addition, all the 8 evolutionary conserved miRNA families involved in flower development [41] are present in the inflorescence of O. italica (Figure 3), even though miR164 is present only as star sequence (miR164*). These results suggest that the basic regulatory pathways of flower development mediated by miRNAs are maintained through evolution. In addition to mature miRNAs, we detected 6 miRNAs* (currently annotated in miRBase as -3p) in the inflorescence of O. italica (Table 3). They have almost the same number of members as their respective mature miRNA, although they are differently accumulated. In particular, there are more miR166*, miR159* and miR171* reads than miR166, miR159 and miR171 reads, respectively, whereas the other miRNAs* have fewer reads than their mature miRNAs. In addition, miR164 is present only as star sequence. These differences in the abundance of mature and star miRNAs might be due to a functional role of miRNA* when the number of reads is high, or to miRNA* degradation due to its exclusion from the silencing complex in favour of the mature miRNA [6], [42], [43].
The list of the evolutionary conserved miRNAs is from Cuperus et al. 2011 [40]. Dotted lines group together the evolutionary conserved miRNAs involved in flower development as described in Luo et al. 2013 [41]. Continuous lines group together the evolutionary conserved miRNAs detected in the inflorescence of O. italica. ∧ indicates miRNA specific of dicots; § indicates miRNA specific of Bryophyta and Lycopodiophyta.
Using the newly assembled transcriptome of P. aphrodite as reference, we looked for novel miRNAs using miRDeep-P, and identified 478 distinct sequences of inflorescence that matched the criteria of plant miRNA structures. Among these, 109 correspond to known plant miRNA families present in miRBase and to miRNA sequences of P. aphrodite, and 100 show similarity with known miRNAs although limited to 11–15 nucleotides. The latter putative miRNAs are not abundant (number of reads ranging from 5 to 154) and might represent known miRNAs belonging to families with variable nucleotide sequences. We analysed the remaining 296 distinct reads using the CD-HIT suite [44], and identified 161 clusters with a sequence identity cut-off of 80%. With the exception of 3 sequences of 26, 28 and 31 nt, respectively, their distribution length ranges between 19–24 nt, with a peak at 21 nt (File S3). These putative novel miRNAs should be considered orchid-specific as their prediction was made using the P. aphrodite transcriptome as reference.
Putative targets
Using a MADS-box degenerated primer and a poly-T primer, we identified four different cDNAs belonging to the DEF-like lineage of the class B MADS-box gene family. These cDNAs range from the first codon (the ATG start codon) to the poly-A tail. OitaDEF1 is 804 bp with an ORF of 681 bp; OitaDEF2 is 893 bp with an ORF of 603 bp; OitaDEF3 is 733 bp with an ORF of 672 bp; OitaDEF4 is 818 bp with and ORF of 675 bp. The NJ tree constructed on the amino acid alignment of the sequence of the virtual translation of these cDNAs with other DEF-like and GLO-like proteins present in GenBank shows that they belong to the known orchid DEF-like clades 1–4 (Figure 4). The tree topology confirms that the DEF-like genes duplications predate the origin of the Orchidaceae [45]. All the OitaDEF-like putative proteins have the typical domain organization of the MIKC transcription factors: the highly conserved MADS-box domain, the moderately conserved intervening I-domain and keratin-like K-domain and the variable C-terminal domain [46].
Numbers represent the bootstrap percentages. Taxon labels reflect the GenBank annotation and do not always correspond to the clade name of the orchid DEF-like lineage. Alpinia hainanensis AHMADS8 [AAT99429], AHMADS5 [AAT99427];Asparagus officinalis AODEF [BAC75969], AOGLOA [BAD13495], AOGLOB [BAD13496]; Cymbidium hybrid cultivar MADS1 Cymb [ABG78568]; Dendrobium crumenatum DcOAP3A [AAZ95248], DcOAP3B [AAZ95249]; Dendrobium thyrsiflorum DthyrPI [AAY86363]; Elaeis guineensis EgGLO1 [AAQ13915], EgGLO2 [AAQ03229], EgDEF1 [AAW66883]; Epipactis palustris EpalPI [AAY56590]; Gongora galeata GogalDEF1 [ACR16036], GogalDEF2 [ACR16037], GogalDEF3 [ACR16038], GogalGLO1 [ACR16039]; Habenaria radiata HrDEF [BAH03320], HrGLO1 [BAH03322], HrGLO2 [BAH03321]; Oncidium hybrid cultivar OMADS3 Onc [AAO45824]; Orchis italica OitaDEF1 [BAO00916], OitaDEF2 [BAO00917http://getentry.ddbj.nig.ac.jp/getentry/dad/BAO00916.1/?filetype=html], OitaDEF3 [BAO00918], OitaDEF4 [BAO00919], OrcPI [BAC22579], OrcPI2 [BAI78360]; Phalaenopsis equestris PeMADS2 [AAR26628], PeMADS3 [AAR26629], PeMADS4 [AAR26626], PeMADS5 [AAR26627], PeMADS6 [AAV28175]; Phalaenopsis hybrid cultivar PhPI9 [AAV28175], PhPI10 [AAV28490], PhPI15 [AAV28491]; Phragmipedium longifolium PhlonDEF1 [ACR16044], PhlonDEF2 [ACR16045], PhlonDEF3 [ACR16046], PhlonDEF4 [ACR16047], PhlonGLO1 [ACR16048]; Spiranthes odorata SpodoDEF1 [ACR16049], SpodoDEF2 [ACR16050], SpodoDEF3 [ACR16051], SpodoDEF4 [ACR16052], SpodoGLO1 [ACR16053]; Vanilla planifolia VaplaDEF1 [ACR16054], VaplaDEF2 [ACR16055], VaplaDEF3 [ACR16056], VaplaGLO1 [ACR16057].
The putative targets of the known and novel miRNAs identified in our study were predicted using the psRNATarget online tool. A total of 5349 putative targets were detected against the transcriptome of P. aphrodite using the default settings of psRNATarget (maximum expectation 3.0). These putative targets correspond to 1456 distinct Gene Ontology (GO) annotations involved in a broad variety of biological processes including cellular and metabolic processes, regulatory processes, developmental processes, response to stimuli, etc. Homologs of known miRNA targets were predicted for several conserved miRNAs of O. italica, including AP2 (target of miR172), Dicer1-like (target of miR162), MYB (target of miR159/319), NF-YA (target of miR169), HD-ZIP transcription factors (targets of miR166), etc.
When the psRNATarget analysis was conducted against the coding transcripts of O. italica, we used a more stringent parameter setting (maximum expectation 0.0), obtaining positive results for OitaAP2 and the two class B MADS-box transcripts OitaDEF2 and OitaDEF4. As expected, the O. italica miRNAs that target the OitaAP2 transcript are homologs of miR172 [47]. The two O. italica miRNAs INF_28502 and INF_26041 are predicted to target the OitaDEF2 and OitaDEF4 transcripts in a corresponding position. Both these miRNAs are homologs of miR5179 of Oryza sativa and Brachypodium distachyon whose function is unknown. Recently, an orchid-specific miRNA that putatively cleaves the transcripts of class B MADS-box genes was identified in E. pusilla and P. aphrodite, without any further characterization [11]. Using a modified 5′ RACE protocol, we successfully amplified a fragment of the expected size (129 bp) and sequence representing the cleavage product of miR5179 on the OitaDEF2 transcript of O. italica (Figure 5), thus confirming that this class B MADS-box gene is a target of the homolog of miR5179 in O. italica. Although the in silico analysis predicted also OitaDEF4 as a possible target of the homolog of miR5179 with the same expectation value of OitaDEF2, we did not detected any cleavage product of this miRNA on OitaDEF4. In addition, even though the nucleotide sequence of the putative target site of miR5179 is quite conserved also in the OitaDEF1 and OitaDEF3 transcripts (Figure 5A), the respective cleavage product was not detectable. This result could be related to the nucleotide differences of the regions surrounding the predicted target site on the four DEF-like transcripts, which results in decreased accessibility of the RISC complex on OitaDEF4 as well as on OitaDEF1 and OitaDEF3. Alternatively, it is conceivable that the activity of miR5179 on the transcript of OitaDEF1, OitaDEF3 and OitaDEF4 might be mediated by translational repression or epigenetic modifications, without a direct cleavage.
A) Nucleotide alignment of the putative target site of miR5179 on OitaDEF2 and its surrounding region with the corresponding region on OitaDEF1, OitaDEF3 and OitaDEF4. The nucleotide sequences span from the position 330 to 410, considering the first nucleotide of the ATG start codon as position 1. B) Agarose gel electrophoresis of the modified 5′ RACE experiment on OitaDEF2 run together with the 100 bp ladder (Fermentas). The alignment of the miR5179 and its target site on OitaDEF4 is reported and 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 (see Methods).
Tissue expression
Ten miRNAs of O. italica, 8 matching known plant miRNAs and 2 putative novel (Table 1), were selected to conduct a relative expression analysis within different inflorescence tissues (Figures 6 and 7). As the Real Time PCR experiments were conducted on the different floral whorls, while the raw reads count was performed on the whole inflorescence, differences were detected between the two datasets. The examined miRNAs belong to families with a different number of members (Figure 2B) and the primers used in the Real Time PCR experiments amplify only part of them. This can further explain the differences between the two datasets.
The predicted structure of the pre-miRNA is reported on the left of each graph, where the miRNA and miRNA* sequences are shown in red and pink, respectively. MFE, minimum free energy; Rn, relative expression ratio; Te_out, outer tepal; Te_inn, inner tepal; Co, column; Ov_np, not pollinated ovary; Ov_3dap, Ov_7dap, Ov_10dap, ovary 3, 7 and 10 days after pollination, respectively; Le, leaf. Bars indicate the standard deviation.
A) Expression profile of the OitaDEF1-4 genes and of the homolog of miR5179 in different tissues of O. italica. The predicted structure of the miR5179 pre-miRNA of O. italica is reported at the rigth side of the expression graph, where the miRNA and miRNA* sequences are shown in red and pink, respectively. MFE, minimum free energy; Rn, relative expression ratio; Te_out, outer tepal; Te_inn, inner tepal; Co, column; Ov_np, not pollinated ovary; Ov_3dap, Ov_7dap, Ov_10dap, ovary 3, 7 and 10 days after pollination, respectively; Le, leaf. Bars indicate the standard deviation. B) From left to right: floral tissues of O. italica used in the expression analysis; single floret and inflorescence of O. italica after anthesis.
The expression pattern of the homolog of miR390 in O. italica (Figure 6A) reveals a not uniform expression profile in the examined tissues. miR390 triggers the biogenesis of the trans-acting silencing RNAs (ta-siRNAs) by directly interacting with the TAS3 mRNA [48]. These ta-siRNAs specifically inhibit the auxine response factor mRNAs ARF2, ARF3 and ARF4, which are involved in the responses to the phytohormone auxin through a pathway that is highly conserved among land plants [48].
The miR160 homolog is highly expressed in all floral organs of O. italica (Figure 6B), particularly in the lip and ovary, which is in agreement with its role in the shaping of the inflorescence and of the floral organs, as well as in fertility [49], [50].
As shown in Figure 6C, the expression of the homolog of miR169 of O. italica is higher in the tepals and lip than in the column, ovary or leaf. The targets of miR169 are NF-YA transcription factors that activate the expression of the class C MADS-box gene AGAMOUS (AG) [51]. Consequently, miR169 restricts, in an indirect manner, the expression domain of the AG gene to the inner floral whorls. The expression pattern of miR169 in O. italica is complementary to that reported for OitaAG in the same floral tissues [19], which demonstrates that the role of miR169 is highly conserved among plants.
The expression of the homolog of miR162 of O. italica is relatively low in all the examined tissues, particularly in the lip and leaf (Figure 6D). miR162 targets the transcript of the Dicer1-like gene that is involved in the biogenesis of the miRNAs [52]. The low levels of miR162 could reflect a high level of processing of the pre-miRNAs within the floral tissues of O. italica.
The expression of the homolog of miR168 of O. italica is higher in tepals, lip, column, not pollinated ovary and leaf than in ovary after pollination (Figure 6E). miR168 cleaves the transcripts of the AGO1 gene that encodes a core component of the RISC complex [53]. In addition, miR168 seems to be involved in stress responses and signal transduction [54].
The homolog of miR166 of O. italica was expressed at variable levels within the examined tissues, with the lowest value detected in the column and the highest in ovary 10 days after pollination (Figure 6F). miR166 targets the mRNA of the HD-ZIP III transcription factors and is involved in shoot apical and lateral meristem formation, organ polarity, and vascular development [55], [56].
The homolog of miR396 of O. italica is expressed at various levels in the floral and leaf tissue; it is low in the ovary before pollination and high in outer tepals and leaf (Figure 6G). miR396 is involved in flower and leaf development targeting the transcripts of growth regulating factors [57].
The expression profile of two putative novel orchid-specific miRNAs is shown in Figure 6H and I. The predicted putative targets of IN_27201 and IN_20892 are a transmembrane kinase and a leucine-rich repeat protein kinase, respectively (Table 1). The expression of both the putative novel miRNAs was lower than that of the selected conserved miRNAs. This is in agreement with reports that the expression of the non-conserved miRNAs is lower compared with that of conserved miRNAs [7], [8].
Finally, as shown in Figure 7A, the miR5179 of O. italica is expressed at low levels in the organs of the perianth. Its expression is considerably higher in ovary tissue than in the other examined tissues, and reaches its maximum level at 10 days after pollination. This pattern is opposite to that of miR5179's putative target transcript OitaDEF2, which is highly expressed in outer tepals, less expressed in inner tepals and lip, whereas its expression is very low in the column and ovary. This complementary expression profile fully fits with the possible repressive role of miR5179 on OitaDEF2. Figure 7A also shows the expression profile of the other DEF-like genes of O. italica. OitaDEF1 is expressed mainly in outer tepals, whereas the expression of OitaDEF3 and OitaDEF4 is higher in inner tepals and lip than in the other tissues.
According to the “orchid code” theory, the class B MADS-box genes belonging to the DEF-like group are responsible for the morphological differences among the organs of the orchid perianth (outer and inner tepals and lip) (Figure 7B) [18], [58]. Fine regulation of the relative expression levels of the four DEF-like genes differentiates the tepals (in which the clade 1 and 2 genes play a predominant role) from the lip (in which the clade 3 and 4 genes play a predominant role) [18]. The expression pattern of the DEF-like genes in O. italica is in agreement with the “orchid code” theory, and supports the concept that the co-expression and the relative abundance of the four DEF-like genes (not simply their presence/absence) is responsible for the morphological diversification of the orchid perianth. In addition, our results suggest that a miRNA regulates the activity of a clade 2 DEF-like gene, which implicates miR5179 in the differentiation of the organs of the perianth in orchids.
Conclusions
The analysis of the inflorescence miRNome of the orchid O. italica revealed the presence of conserved and novel miRNAs. The in silico search for the possible miRNA targets showed a conserved miRNA cleavage site within the four OitaDEF-like transcripts, which we experimentally validated for OitaDEF2. This result suggests that miRNAs play an important role in the diversification of the organs of the perianth in orchids through the inhibitory regulation of the clade 2 DEF-like gene, and highlights the potential of next-generation sequencing in the study of non model plant species. Different mechanisms might act to regulate the expression level of the other DEF-like genes, suggesting the existence of lineage-specific regulatory mechanisms contributing to the functional specialization of the DEF-like clades in orchids.
Supporting Information
File S1.
Distinct reads of the small RNA library of inflorescence of O. italica after pre-processing. The file is in fasta format and for each sequence the number of reads is indicated after “_x”.
https://doi.org/10.1371/journal.pone.0097839.s001
(TXT)
File S2.
Results of the BLAST analysis of the cleaned reads of O. italica versus the plant miRNAs annotated in miRBase20.
https://doi.org/10.1371/journal.pone.0097839.s002
(XLSX)
File S3.
Results of the BLAST analysis of the cleaned reads of O. italica versus the miRNAs of floral buds of P. aphrodite.
https://doi.org/10.1371/journal.pone.0097839.s003
(XLSX)
File S4.
Length distribution of the clustered orchid-specific putative miRNAs in the inflorescence of O. italica.
https://doi.org/10.1371/journal.pone.0097839.s004
(TIF)
Acknowledgments
The authors are grateful to Dr. Marinella Salemme for her contribution to the RNA extraction, to Dr. Maria Rosaria Barone Lumaga for the description of the ovary development stages and to Prof. Salvatore Cozzolino for supplying the plant material and performing the manual pollination of ovaries in O. italica. The authors thank Mrs Rosaria Terracciano and Mr Enzo Iacueo for their technical support, and Jean Ann Gilder (Scientific Communication, srl., Naples, Italy).
Author Contributions
Conceived and designed the experiments: SA LG. Performed the experiments: MS SDP VD PC. Analyzed the data: SA. Contributed reagents/materials/analysis tools: SA LG FS. Wrote the paper: SA LG FS.
References
- 1. Axtell MJ (2013) Classification and comparison of small RNAs from plants. Annu Rev Plant Biol 64: 137–159.
- 2. Wu G (2013) Plant microRNAs and development. J Genet Genomics 40: 217–230.
- 3. Chuck G, Candela H, Hake S (2009) Big impacts by small RNAs in plant development. Curr Opin Plant Biol 12: 81–86.
- 4. Chen X (2005) MicroRNA biogenesis and function in plants. FEBS Lett 579: 5923–5931.
- 5. Nag A, Jack T (2010) Sculpting the flower; the role of microRNAs in flower development. Curr Top Dev Biol 91: 349–378.
- 6. Lin Y, Lai Z (2013) Comparative analysis reveals dynamic changes in miRNAs and their targets and expression during somatic embryogenesis in longan (Dimocarpus longan Lour.). PLoS One 8: e60337.
- 7. Khraiwesh B, Pugalenthi G, Fedoroff NV (2013) Identification and analysis of red sea mangrove (Avicennia marina) microRNAs by high-throughput sequencing and their association with stress responses. PLoS One 8: e60774.
- 8.
Zhu J, Li W, Yang W, Qi L, Han S (2013) Identification of microRNAs in Caragana intermedia by high-throughput sequencing and expression analysis of 12 microRNAs and their targets under salt stress. Plant Cell Rep.
- 9. Thiebaut F, Grativol C, Carnavale-Bottino M, Rojas CA, Tanurdzic M, et al. (2012) Computational identification and analysis of novel sugarcane microRNAs. BMC Genomics 13: 290.
- 10. An FM, Hsiao SR, Chan MT (2011) Sequencing-based approaches reveal low ambient temperature-responsive and tissue-specific microRNAs in phalaenopsis orchid. PLoS One 6: e18937.
- 11. Lin CS, Chen JJ, Huang YT, Hsu CT, Lu HC, et al. (2013) Catalog of Erycina pusilla miRNA and categorization of reproductive phase-related miRNAs and their target gene families. Plant Mol Biol 82: 193–204.
- 12. Coen ES, Meyerowitz EM (1991) The War of the Whorls - Genetic Interactions Controlling Flower Development. Nature 353: 31–37.
- 13. Krizek BA, Fletcher JC (2005) Molecular mechanisms of flower development: An armchair guide. Nature Reviews Genetics 6: 688–698.
- 14. Causier B, Schwarz-Sommer Z, Davies B (2010) Floral organ identity: 20 years of ABCs. Semin Cell Dev Biol 21: 73–79.
- 15. Kanno A, Nakada M, Akita Y, Hirai M (2007) Class B gene expression and the modified ABC model in nongrass monocots. ScientificWorldJournal 7: 268–279.
- 16. Salemme M, Sica M, Gaudio L, Aceto S (2011) Expression pattern of two paralogs of the PI/GLO-like locus during Orchis italica (Orchidaceae, Orchidinae) flower development. Dev Genes Evol 221: 241–246.
- 17. Cantone C, Sica M, Gaudio L, Aceto S (2009) The OrcPI locus: Genomic organization, expression pattern, and noncoding regions variability in Orchis italica (Orchidaceae) and related species. Gene 434: 9–15.
- 18. Mondragon-Palomino M, Theissen G (2011) Conserved differential expression of paralogous DEFICIENS- and GLOBOSA-like MADS-box genes in the flowers of Orchidaceae: refining the ‘orchid code’. Plant J 66: 1008–1019.
- 19. Salemme M, Sica M, Gaudio L, Aceto S (2013) The OitaAG and OitaSTK genes of the orchid Orchis italica: a comparative analysis with other C- and D-class MADS-box genes. Mol Biol Rep 40: 3523–3535.
- 20. Aceto S, Gaudio L (2011) The MADS and the Beauty: Genes Involved in the Development of Orchid Flowers. Curr Genomics 12: 342–356.
- 21. Salemme M, Sica M, Iazzetti G, Gaudio L, Aceto S (2013) The AP2-Like Gene OitaAP2 Is Alternatively Spliced and Differentially Expressed in Inflorescence and Vegetative Tissues of the Orchid Orchis italica. PLoS One 8: e77454.
- 22. Moxon S, Schwach F, Dalmay T, Maclean D, Studholme DJ, et al. (2008) A toolkit for analysing large-scale plant small RNA datasets. Bioinformatics 24: 2252–2253.
- 23. Stocks MB, Moxon S, Mapleson D, Woolfenden HC, Mohorianu I, et al. (2012) The UEA sRNA workbench: a suite of tools for analysing and visualizing next generation sequencing microRNA and small RNA datasets. Bioinformatics 28: 2059–2061.
- 24. Kozomara A, Griffiths-Jones S (2011) miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res 39: D152–157.
- 25. Yang X, Li L (2011) miRDeep-P: a computational tool for analyzing the microRNA transcriptome in plants. Bioinformatics 27: 2614–2615.
- 26. Sedeek KE, Qi W, Schauer MA, Gupta AK, Poveda L, et al. (2013) Transcriptome and proteome data reveal candidate genes for pollinator attraction in sexually deceptive orchids. PLoS One 8: e64621.
- 27. Zhang J, Wu K, Zeng S, Teixeira da Silva JA, Zhao X, et al. (2013) Transcriptome analysis of Cymbidium sinense and its application to the identification of genes associated with floral development. BMC Genomics 14: 279.
- 28. Chou ML, Shih MC, Chan MT, Liao SY, Hsu CT, et al. (2013) Global transcriptome analysis and identification of a CONSTANS-like gene family in the orchid Erycina pusilla. Planta 237: 1425–1441.
- 29. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, et al. (2011) Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol 29: 644–652.
- 30. Chen TW, Gan RC, Wu TH, Huang PJ, Lee CY, et al. (2012) FastAnnotator—an efficient transcript annotation web tool. BMC Genomics 13 Suppl 7S9.
- 31. Dai X, Zhao PX (2011) psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res 39: W155–159.
- 32. Cantone C, Gaudio L, Aceto S (2011) The GLO-like locus in orchids: duplication and purifying selection at synonymous sites within Orchidinae (Orchidaceae). Gene 481: 48–55.
- 33. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, et al. (2011) MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol 28: 2731–2739.
- 34. Llave C, Xie Z, Kasschau KD, Carrington JC (2002) Cleavage of Scarecrow-like mRNA targets directed by a class of Arabidopsis miRNA. Science 297: 2053–2056.
- 35. Shi R, Sun YH, Zhang XH, Chiang VL (2012) Poly(T) adaptor RT-PCR. Methods Mol Biol 822: 53–66.
- 36. Zhao S, Fernald RD (2005) Comprehensive algorithm for quantitative real-time polymerase chain reaction. J Comput Biol 12: 1047–1064.
- 37. Fahlgren N, Howell MD, Kasschau KD, Chapman EJ, Sullivan CM, et al. (2007) High-throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of MIRNA genes. PLoS ONE 2: e219.
- 38. Zhang R, Marshall D, Bryan GJ, Hornyik C (2013) Identification and characterization of miRNA transcriptome in potato by high-throughput sequencing. PLoS One 8: e57233.
- 39. Morin RD, Aksay G, Dolgosheina E, Ebhardt HA, Magrini V, et al. (2008) Comparative analysis of the small RNA transcriptomes of Pinus contorta and Oryza sativa. Genome Res 18: 571–584.
- 40. Cuperus JT, Fahlgren N, Carrington JC (2011) Evolution and functional diversification of MIRNA genes. Plant Cell 23: 431–442.
- 41. Luo Y, Guo Z, Li L (2013) Evolutionary conservation of microRNA regulatory programs in plant flower development. Dev Biol 380: 133–144.
- 42. Jones-Rhoades MW, Bartel DP, Bartel B (2006) MicroRNAS and their regulatory roles in plants. Annu Rev Plant Biol 57: 19–53.
- 43. Manavella PA, Koenig D, Rubio-Somoza I, Burbano HA, Becker C, et al. (2013) Tissue-specific silencing of Arabidopsis SU(VAR)3-9 HOMOLOG8 by miR171a. Plant Physiol 161: 805–812.
- 44. Huang Y, Niu B, Gao Y, Fu L, Li W (2010) CD-HIT Suite: a web server for clustering and comparing biological sequences. Bioinformatics 26: 680–682.
- 45. Pan ZJ, Cheng CC, Tsai WC, Chung MC, Chen WH, et al. (2011) The duplicated B-class MADS-box genes display dualistic characters in orchid floral organ identity and growth. Plant Cell Physiol 52: 1515–1531.
- 46. Kaufmann K, Melzer R, Theissen G (2005) MIKC-type MADS-domain proteins: structural modularity, protein interactions and network evolution in land plants. Gene 347: 183–198.
- 47.
Salemme M, Sica M, Iazzetti G, Gaudio L, Aceto S (2013) The AP2-like gene OitaAP2 is alternatively spliced and differentially expressed in inflorescence and vegetative tissues of the orchid Orchis italica. PLoS One.
- 48. Jouannet V, Moreno AB, Elmayan T, Vaucheret H, Crespi MD, et al. (2012) Cytoplasmic Arabidopsis AGO7 accumulates in membrane-associated siRNA bodies and is required for ta-siRNA biogenesis. EMBO J 31: 1704–1713.
- 49. Liu PP, Montgomery TA, Fahlgren N, Kasschau KD, Nonogaki H, et al. (2007) Repression of AUXIN RESPONSE FACTOR10 by microRNA160 is critical for seed germination and post-germination stages. Plant J 52: 133–146.
- 50. Mallory AC, Bartel DP, Bartel B (2005) MicroRNA-directed regulation of Arabidopsis AUXIN RESPONSE FACTOR17 is essential for proper development and modulates expression of early auxin response genes. Plant Cell 17: 1360–1375.
- 51. Hong RL, Hamaguchi L, Busch MA, Weigel D (2003) Regulatory elements of the floral homeotic gene AGAMOUS identified by phylogenetic footprinting and shadowing. Plant Cell 15: 1296–1309.
- 52. Hirsch J, Lefort V, Vankersschaver M, Boualem A, Lucas A, et al. (2006) Characterization of 43 non-protein-coding mRNA genes in Arabidopsis, including the MIR162a-derived transcripts. Plant Physiol 140: 1192–1204.
- 53. Vaucheret H (2009) AGO1 homeostasis involves differential production of 21-nt and 22-nt miR168 species by MIR168a and MIR168b. PLoS One 4: e6442.
- 54. Li W, Cui X, Meng Z, Huang X, Xie Q, et al. (2012) Transcriptional regulation of Arabidopsis MIR168a and argonaute1 homeostasis in abscisic acid and abiotic stress responses. Plant Physiol 158: 1279–1292.
- 55. Miyashima S, Honda M, Hashimoto K, Tatematsu K, Hashimoto T, et al. (2013) A comprehensive expression analysis of the Arabidopsis MICRORNA165/6 gene family during embryogenesis reveals a conserved role in meristem specification and a non-cell-autonomous function. Plant Cell Physiol 54: 375–384.
- 56. Husbands AY, Chitwood DH, Plavskin Y, Timmermans MC (2009) Signals and prepatterns: new insights into organ polarity in plants. Genes Dev 23: 1986–1997.
- 57. Debernardi JM, Rodriguez RE, Mecchia MA, Palatnik JF (2012) Functional specialization of the plant miR396 regulatory network through distinct microRNA-target interactions. PLoS Genet 8: e1002419.
- 58. Mondragon-Palomino M, Theissen G (2009) Why are orchid flowers so diverse? Reduction of evolutionary constraints by paralogues of class B floral homeotic genes. Ann Bot (Lond) 104: 583–594.