MicroRNAs (miRNAs) are small noncoding RNAs that regulate gene expression by posttranscriptional silencing mechanisms in plants and animals and control development, disease resistance, and stress responses in plants (Casadevall et al., 2013; Li et al., 2010; Wu et al., 2006; Zhu and Helliwell, 2011). Plant miRNA genes are transcribed by RNA polymerase II as primary miRNAs that are processed to precursor miRNAs (pre-miRNAs) via Dicer-like 1 (DCL1) in the nucleus. Pre-miRNAs are exported to the cytoplasm where DCL1 cleaves the pre-miRNAs near its stem-loop, giving rise to the characteristic 21–23 nucleotide (nt) long double-strand RNAs with two nucleotides 3′ overhangs. These miRNA duplexes are incorporated into the Argonaute (AGO) protein to form the effector complex, which is referred to as the RNA-induced silencing complex (RISC). The two strands of the miRNA duplex are then separated in AGO, with one strand removed and the other one retained in RISC for target recognition and silencing. In plants, most miRNAs exhibit near-perfect complementarity to their targets, often resulting in mRNA cleavage (Rogers and Chen, 2013). Because of the small size, low abundance, and tissue- or developmental stage-specific expression patterns of miRNAs, coupled with stress-induced dynamic changes in their expression, experimentally validating miRNAs has been challenging. Since the advent of next-generation sequencing, high-throughput sequencing data have become available, and computational approaches of processing large-scale data have led to the discovery of a large number of miRNA sequences (Hwang et al., 2013; Kim et al., 2012).

Hibiscus syriacus, also known as the rose of Sharon, is the national flower of South Korea and belongs to the member of the Malvaceae family such as cotton and cacao. H. syriacus has significance in gardening, medicine, and cultural perspective in Korea. Meanwhile, Mugunghwa, which means “flowering forever” and is its Korean name, is derived from its unique characteristic that it blooms and falls away in several months. A recent study provided the draft genome of H. syriacus and suggested that the flowering-related phenotype might be associated with polyploidization (Kim et al., 2017). Despite its interesting features, the miRNA transcriptome of H. syriacus remains unknown. Therefore, identifying miRNAs and their targets in H. syriacus could help to understand their biological functions.

Taking advantage of the recent genome sequencing in H. syriacus (Kim et al., 2017), we employed high-throughput sequencing to identify conserved and novel miRNAs in H. syriacus via a computational approach. We successfully identified 93 conserved and 36 novel miRNAs from four different tissues. We also predicted their corresponding target genes and further validated several targets using the 5′ RNA ligase-mediated rapid amplification of cDNA ends (RLM-RACE) analysis. One of the predicted targets of conserved miRNAs, miR396, is SHORT VEGETATIVE PHASE (SVP), which is involved in flower initiation and is duplicated in H. syriacus (Kim et al., 2017; Yang et al., 2015). In addition, we showed that miR477 directs the cleavage of terpene synthase, which has functions in isoprenoid synthesis. Altogether, our study lays the foundation for understanding the functions of miRNAs in H. syriacus and provides the beneficial data regarding molecular breeding.

Construction of small RNA library

Total RNAs were extracted from H. syriacus tissue samples (leaf, root, flower, and ovary) using TRI Reagent solution (MRC), according to the manufacturer’s protocol. Both the flower and ovary were pretreated with CTAB buffer [2% CTAB, 2% polyvinylpyrrolidone K 30, 100 mM Tris (pH 8.0), 25 mM EDTA, 2 M NaCl, and 2% β-mercaptoethanol] followed by TRI Reagent extraction, as previously described (Peng et al., 2014). A small RNA sequencing library was constructed using the Small RNA Sample Prep Kit (Illumina), according to the manufacturer’s instructions. As previously described (Hwang et al., 2013), 10 μg of total RNAs was resolved using 15% urea-PAGE, and size fractions of small RNAs between 18 and 30 nt were prepared. The prepared small RNAs were ligated to the 3′ and 5′ adapter, and cDNA was synthesized using RT-PCR. cDNAs were amplified using PCR and submitted for Illumina/Solexa sequencing. The GEO accession number for our series is GSE99329.

Identification of miRNA and target prediction

The pipeline of miRNA identification was written in Python language, and the details of the pipeline are shown in Supplementary Fig. S1. Poor-quality small RNA reads and adapter sequences were eliminated using the FASTX toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html), and the reads whose length was between 18 and 26 nt were selected for further analysis. High-quality small RNAs (reads of ≥ 50) were mapped to the reference genome of H. syriacus (Kim et al., 2017) using Bowtie-1.1.2 (Langmead et al., 2009), and the perfectly mapped sequences were selected (≤ 42 times). The filtered reads were analyzed to identify the putative miRNAs in the modified algorithm (Hwang et al., 2013; Meyers et al., 2008). After identification, the miRNA candidates were classified to either the conserved miRNAs or novel miRNAs through a BLASTN search (Camacho et al., 2009) with miRBase release 21.0 (http://www.mirbase.org/) (Kozomara and Griffiths-Jones, 2014) and filtered with several criteria (Supplementary Fig. S1) (Hwang et al., 2013; Trotta, 2014). The reads of miRNAs from each library were normalized to read per 10 million (RP10M). The potential target genes of miRNAs were predicted using psRNAtarget with a penalty score of ≤ 5 (Dai and Zhao, 2011).

miRNA target validation assays

For miRNA target validation, gene-specific 5′ RLM-RACE was performed as previously described (Hwang et al., 2013). In brief, 5 μg of pooled total RNA obtained from leaf, root, flower, and ovary tissues were ligated to 5′ the GeneRacer RNA adapter (Invitrogen). Reverse transcription producing cDNAs were synthesized with both Oligo (dT) primers and random hexamers. Primary PCR was performed with the GeneRacer 5′ primer and 3′ reverse primers containing each respective gene-specific sequence (Supplementary Table S1). Before cloning into a TA cloning vector (RBC) for sequencing, nested PCR was performed using the GeneRacer 5′ nested primer and each respective gene-specific nested reverse primer (Supplementary Table S1).

Small RNA analysis of H. syriacus

Small RNA libraries from leaf, root, flower and ovary were constructed and processed using a computation pipeline, as described in Supplementary Fig. S1. From this, 20–27 million raw reads were obtained from each library (Table 1). After filtering low-quality reads and those with no adapter sequences, 16–25 million reads remained, representing 10–14 million unique sequences (Table 1). These small RNA reads were mapped to the H. syriacus genome for small RNA identification. Small RNA sequences of 21 and 24 nt in length were dominant in all four libraries (Fig. 1), which is mostly consistent with those previously reported on other plant species (Hwang et al., 2013; Kim et al., 2012; Pantaleo et al., 2010). For the 21-nt reads, the 5′-end U was dominant in all four libraries, followed by A. In contrast, the 5′-end A was the major form in 24-nt reads, followed by U, representing the canonical small RNA distribution (Fig. 1). For example, Arabidopsis thaliana miRNAs are processed by DCL1 and loaded into AGO1; a strong bias toward a 5′-terminal U is observed. The 24-nt small RNAs are processed by DCL3 and preferred by AGO4 and have sequences with 5′-terminal A (Voinnet, 2009).

Identification of miRNAs in H. syriacus

Aligned small RNA sequences were analyzed to identify candidate miRNAs on the basis of the strict criteria for annotation of plant miRNA (Meyers et al., 2008). To classify conserved and novel miRNAs in H. syriacus, the candidates were BLASTN searched against the reference miRNAs in miRBase (release 21.0), allowing two or less mismatches within mature miRNA sequences. Those sequences that were not matched were manually checked with precursor similarity to assess whether they were truly novel miRNAs or belonged to conserved families. Representative and total miRNAs are listed in Table 2 and Supplementary Table S2, respectively. Overall, 93 conserved miRNAs belonging to 33 families and 36 novel miRNAs from 30 families were annotated. The representative structures of miRNA precursors are shown in Fig. 2A. Although most miRNA precursors encode a single mature miRNA, a few encodes two different mature miRNAs such as hsy-miR-n006a-b-5p.2 and hsy-miR-n006a-5p.1, both of which are located 42-nt apart (Fig. 2B). This phenomenon has been reported in several other plants, including A. thaliana, Oryza sativa, Medicago truncatula (Zhang et al., 2010), possibly because of the sequential processing of Dicer-like proteins.

There are many miRNA copies in the H. syriacus genome compared with those in other plant species. For instance, miR166 is a well-conserved miRNA in various plant species (Cuperus et al., 2011). While there are seven copies in the A. thaliana genome, H. syriacus has 50 copies (Supplementary Table S3). This is much greater than that of Glycine max (21 copies), which is known to have many copy numbers of miRNAs (Nozawa et al., 2012). miR171 is another highly conserved miRNA with three copies in A. thaliana, 48 copies in H. syriacus, and 21 copies in G. max genome (Supplementary Table S3). Although the depth for different libraries and pipelines for identification are slightly different, the copy numbers of several miRNAs in H. syriacus are likely greater than those in other species. In contrast, with respect to miR156, there are eight copies in A. thaliana, 27 copies in H. syriacus, but 28 copies in G. max (Supplementary Table S3). A possible reason for this could be whole-genome duplication (WGD) events. After the last speciation, the H. syriacus genome underwent two WGD and diploidization events. Subsequently, the dosage of the individual gene family was differently regulated. Therefore, the H. syriacus species has more gene loci than other Malvaceae plants or Arabidopsis (Kim et al., 2017). Likewise, WGD and diploidization might cause miRNA copy number variations in H. syriacus. Conversely, almost all of the novel miRNAs are present in one or two copies in the genome. The novel miRNA genes might have been generated after speciation, and thus, may experience less WGD events than conserved miRNAs, possibly explaining the reason for the numbers of novel miRNAs being less than those of conserved miRNAs.

miRNA expression patterns of H. syriacus

As normalized to RP10M, expression patterns of conserved and novel miRNAs from four tissues are shown in Fig. 3. First, Figs. 3A and 3B both represent conserved miRNAs with at least two-fold differences between their respective minimum and maximum expressions. Second, the top 10 most highly expressed miRNAs are shown in Fig. 3A, and the others are shown in Fig. 3B. Third, novel miRNAs are shown in Fig. 3C. Finally, some miRNAs that are consistently expressed throughout all four tissues are shown in Fig. 3D. Overall, many conserved miRNAs are expressed higher than novel miRNAs with miR156a-f-3p, miR319a-j-3p, miR390a-o-5p, miR393a-n-5p, and miR396a-i-5p being the most highly expressed (Table 2). The only highly expressed novel miRNA was miR-n001a-b-3p with 92531 reads in leaf tissues. These results indicate that conserved miRNAs are expected to play more significant roles than novel miRNAs in regulatory pathways and have a much wider range of dynamic expression for various regulations.

Many conserved miRNAs show tissue-specific expression patterns. For example, miR156a-g-5p expression was distinctively higher in the leaf and root than in the flower and ovary (Fig. 3A). Another member of miR156, miR156h-p-5p, was not highly expressed in the leaf, but it was highly expressed in the root (Fig. 3A). The differential expression of these may indicate that their biological roles are slightly different, even within the same family members. In contrast, miR172a-h-3p was exclusively expressed in the flower and ovary. Previous studies showed that miR156 is an upstream regulator of miR172 and acts by repressing miR172 through its target, SQUAMOSA PROMOTER-BINDING PROTEIN-LIKE 9 (SPL9) in A. thaliana. This in turn activates miR172 transcription (Wu et al., 2009). In the juvenile phase, miR156 represses SPL, which is necessary for leaf growth (Wu et al., 2009) and lateral root development (Yu et al., 2015). SPL9 consequently activates miR172, which promotes an adult phase transition by repressing its targets, TOE1 and TOE2 (Wu et al., 2009). The expression patterns of these miRNAs (Fig. 3A–3B) are consistent with previous studies (Wu et al., 2009; Yu et al., 2015), indicating that miR156 and miR172 may function similarly in H. syriacus.

miR396a-i-5p expression in the leaf was much higher than that in other tissues (Fig. 3A). miR396 regulates the GROWTH-REGULATING FACTOR (GRF) families to control cell proliferation in the A. thaliana leaf (Rodriguez et al., 2010). In leaf tissues, miR396 expression levels steadily increase during leaf maturation, attenuating cell proliferation (Rodriguez et al., 2010). This is in line with our results that miR396 expression was the highest in the leaf, suggesting that its functional role is conserved in H. syriacus. Moreover, another study showed that when A. thaliana is exposed to UV-B radiation, miR396 is expressed to inhibit leaf growth by GRFs repressions (Casadevall et al., 2013). Because our samples were collected from the open field, miR396 might be highly expressed highly in the leaf.

Similar to miR172, miR319a-j-3p expression was much higher in the flower and ovary than in the leaf and root (Fig. 3A), whereas miR319k-n-3p expression (Fig. 3B) was higher in the root than in any other tissues. Previous studies showed that the loss-of-function miR319 mutants of A. thaliana exhibit defects in petal and stamen development, with narrower and shorter petals and impaired anther formation (Nag et al., 2009; Rubio-Somoza and Weigel, 2013). In addition, each member of the miR319 family is expressed in different tissues. In A. thaliana, miR319a is highly expressed in the petals and stamens during early flower development but not in the leaves, whereas miR319c is expressed in young leaf primordia (Nag et al., 2009). High miR319a-j 3p expression levels in the flower and ovary may indicate a conserved or similar control mechanism that also exists in the flower development of H. syriacus. Furthermore, a recent study showed that transgenic Agrostis stolonifera that overexpressed Osa-miR319a had less root biomass and better salt and drought tolerance than mock controls (Zhou et al., 2013). Another member of the miR319 family, miR319k-n-3p is highly expressed in the root, which might be associated with this phenomenon.

Some miRNAs were consistently expressed in all four tissues (Fig. 3D). For instance, miR166a-ag-3p showed quite high expression patterns in every tissue. Mutation studies demonstrated that miR166 overexpression affected several class III homeodomain-leucine zipper family genes, possibly causing defects in the shoot meristem, leaves, gynoecia, and vascular tissues of A. thaliana (Kim et al., 2005; Williams et al., 2005). In addition, miR166 controls root and nodule development in M. truncatula (Boualem et al., 2008). These results indicated that miR166 is possibly associated with the development of several tissues, which is in accordance with the fact that miR166 is expressed consistently in all tissues.

Novel miRNAs generally showed lower and less dynamic expression than those of conserved miRNAs (Figs. 3A–3C), suggesting their species-specific roles in particular circumstances. However, a few novel miRNAs exhibited dynamic expression patterns; for example, miR-n001a-b-3p and miR-n006 family genes were highly expressed in leaf, whereas miR-n003a-e-3p displayed the highest expression in flower. These results may indicate tissue-specific roles of some novel miRNA families.

miRNA target analysis of H. syriacus

To investigate the potential regulatory roles of miRNAs in H. Syriacus, we predicted the putative target genes (Supplementary Table S4) using psRNAtarget with penalty scores of five or less (Dai and Zhao, 2011). The representative targets of 17 conserved and 28 novel miRNA families are shown (Table 3). Five targets of some conserved miRNAs were validated by 5′ RLM-RACE analysis (Fig. 4) (Park and Shin, 2014).

Noticeable targets of conserved miRNAs include many transcription factors such as SPL, no apical meristem domain-containing proteins, auxin response factors (ARFs), GRFs, and MYBs. Moreover, there are other conserved targets, including auxin signaling F-Box (AFB) and nucleotide-binding leucine-rich repeat (NB-LRR). These results indicate that miRNAs in H. syriacus are likely to have conserved roles in the regulation of growth, development, and disease resistance. For example, SPLs are the major targets of miR156, as an intermediate of the miR156-SPL-miR172 pathway that controls the developmental timing via feedback regulations (Wu et al., 2009). The APETALA2 (AP2)-like are regulated by miR172 via translation inhibition, thereby controlling the flowering time and floral organ identity (Aukerman and Sakai, 2003). ARFs are the targets of miR160 and miR167 and regulate sexual reproduction (Wu et al., 2006) and adventitious rooting (Gutierrez et al., 2009). Other auxin-associated genes, AFBs, under the control of miR393, regulate leaf development via secondary siRNAs (Si-Ammour et al., 2011) and the root system architecture under nitrate treatment (Vidal et al., 2010) in A. thaliana. In another case, TEOSINTE BRANCHED1, CYCLOIDEA, and PCF (TCP) is mainly regulated by miR319, an important player of the morphological integrity of the plant, including shoot lateral organs (Koyama et al., 2007). In addition, miR319 and its target, TCPs, are related to flower and reproductive tissues development (Nag et al., 2009; Rubio-Somoza and Weigel, 2013). NB-LRRs, which are targeted by miR482, recognize pathogen-encoded effectors and initiate an effector-triggered immunity. The NB-LRR domain represents a general character of disease resistance proteins in plants acquired from co-evolutionary arms race between pathogens and their host plants (Park and Shin, 2015).

In this study, we first discovered a novel target of the conserved miRNA, miR477. miR477c-e-5p was predicted to regulate the terpene synthase gene, which was validated by 5′ RLM-RACE analysis (Fig. 4E). Terpene synthase is an enzyme for synthesizing terpene, which is associated with the plant innate immunity (Wang et al., 2012). Thus, we expect these results to provide a foundation for further studies regarding the control of terpene synthesis. The recent genome study of H. syriacus suggested that flowering-related genes were specifically expanded in H. syriacus (Kim et al., 2017). One of the expanded genes is SVP that controls flowering time by repressing FLOWERING LOCUS T expression (Lee et al., 2007). A recent study showed that miR396 triggers mRNA decay of SVP, when the phyllody symptoms 1 effector is treated to Catharanthus roseus, which causes a leafy or abnormal flower phenotype (Yang et al., 2015). These studies suggest the SVP represses flowering and miR396 is an upstream regulator of SVP. We found that expanded SVPs are predicted targets of miR396 in H. syriacus (Table 3), which raises an intriguing possibility that SVPs are associated with a specific flowering pattern of H. syriacus.

As described above, we proposed that WGD might affect the expansion of some miRNA families and protein-coding genes in H. syriacus (Kim et al., 2017). Therefore, we investigated whether isoforms of the expanded miRNA families regulate the same target or other target genes to deepen our comprehension of the impact of WGD on H. syriacus. In the case of miR156 family, miR156a-g-5p showed higher expression in vegetative tissues than that in reproductive tissues and miR156h-p-5p was highly expressed in root, respectively (Fig. 3A). While there are two nucleotide differences between these miRNA isoforms (Supplementary Table S2), their predicted target repertoires are almost overlapped with each other (Supplementary Table S4). Although the repertoires are nearly identical, details are quite different in terms of target affinity. For example, one of the predicted targets (MGA0.9.scaffold2887G00090;1) belonging to SPL family has different penalty scores against each miRNA isoform. The penalty scores are 2.5 against miR156a-g-5p and 1 against miR156h-p-5p, respectively. Likewise, there are many other genes, such as MGA0.9. scaffold4486G00120;1 and MGA0.9.scaffold2429G00200;1 (Supplementary Table S4). Collecting the differences in expression and target affinity of the miRNA isoforms, it indicates the possibility that miRNA target genes are differently regulated by miRNA isoforms in different tissues. Another case supports this possibility. miR319-a-j-3p expression was higher in reproductive tissues (Fig. 3A), whereas miR319k-n-3p expression was higher in root. Similarly, there are two nucleotide differences between these miRNA isoforms (Supplementary Table S2), but their target repertoires are identical and target affinities are different (Supplementary Table S4). MGA0.9.scaffold21723G00010;1 is GAMYB-like gene and it prefer miR-319a-j-3p (penalty score: 1.5) to miR319-k-n-3p (penalty score: 4.5) as a regulator. Consequently, the target gene might be strongly repressed in reproductive tissues and weakly repressed in root. In conclusion, WGD might provide more complicated layer of control of genes in H. syriacus compared with other species.

Novel miRNA targets include pentatricopeptide repeat-containing (PPR) proteins, serine/threonine-protein kinase, receptor-like protein (RLP) and transcription factor MYC2-like. Most novel miRNAs exhibited weaker expressions than conserved miRNAs. The only exception was miR-n001a-b-3p, which is expressed as much as conserved miRNAs (Table 3), indicating that miR-n001a-b-3p may have an important biological role in H. syriacus. A predicted target of miR-n001a-b-3p is the PPR. The PPR family is one of the largest protein families in land plants and plays a role in plant development and growth (Barkan and Small, 2014). The PPR protein typically influences organelle biogenesis and function, including that for chloroplasts or mitochondria (Barkan and Small, 2014). miR-n001a-b-3p is predicted to target many PPRs, unlike other novel miRNAs (Supplementary Table S4). In addition, miR-n001a-b-3p is highly expressed in the leaf compared to that in any other tissues (Fig. 3). Thus, miR-n001a-b-3p might be involved in leaf development or photosynthesis. Furthermore, miR-n007a-5p is expected to target PPRs (Table 3), and its expression pattern was similar to that of miR-n001a-b-3p. Therefore, its function probably resembles the function of miR-n001ab-3p. RLP is another target of novel miRNAs and is a subcategory of receptor-like protein kinase; it recognizes pathogens to trigger plant defense mechanisms or control development (Afzal et al., 2008). In our prediction, miRn006a-5p.1 and miR-n006b-5p.1 target RLPs. Therefore, they might be involved in the immunity of H. syriacus with miR482.

In contrast to conserved miRNAs, we were unable to validate some targets of novel miRNAs via the 5′ RLM-RACE analysis (data not shown). It is in agreement with earlier studies that the targets of novel miRNAs are generally not identifiable (Hwang et al., 2013; Moxon et al., 2008). There are possible reasons as follows: (1) the abundance of most novel miRNAs were not high enough to direct the cleavage of mRNA targets, (2) target mRNAs and miRNAs were not colocalized or coexpressed (Cuperus et al., 2011), (3) the mode of novel miRNA action is repression at the protein level (Aukerman and Sakai, 2003), (4) the gene annotation of H. syriacus was not complete, and (5) most targets of novel miRNAs had higher penalty scores than those of conserved miRNAs (i.e., a relatively poor sequence complementarity) (Hwang et al., 2013).

In this study, we discovered 93 conserved and 36 novel miRNA sequences from a small RNA library of the novel species H. syriacus by implementing and applying an end-to-end pipeline for the analysis. In addition, we identified miRNA targets of H. syriacus in silico. Targets of many conserved miRNAs and several novel miRNAs are consistent with previous studies, showing that H. syriacus shares a lot of biological properties and pathways with model plants. Some miRNAs and their targets showed potential of species-specific regulations, including flower development, terpene synthesis, and disease resistance. Differential expression profiles revealed a tissue specificity of many miRNA sequences, implying the importance of each sequence for orchestrating gene expression at different sites and scales. Future work will entail experimental validation of miRNA sequences and their target cleavage. Finally, functional analysis of the target genes described in this study could provide a deeper interpretation of miRNA-mediated regulation in H. syriacus.