An overview of high-throughput sequencing data sets
To identify the miRNA in P. heterophylla, a small RNA library from four organs (leaves, stem, tuberous root cortex and xylem) was constructed and sequenced independently. A total of 180 million raw reads were obtained from these 12 libraries. These raw data were filtered through several criteria to identify conserved and specific miRNAs.
The length distribution of the small RNAs was shown in Fig 1. The majority of reads were 21 to 24 nt in length. The 24-nt RNA is the most abundant class in the sRNA library and comprised 35.7% of root bark, 32.7% of root xylem, 31.5% of leaf and 31.8 % stem. Similar results were reported in potato, Arabidopsis thaliana, tomato, pepper, cucumber [29], but in contrast to grapevine, wheat [30].
The transcriptome of P. heterophylla
We obtain 603, 398, 664 raw reads and 598, 156, 596clean reads from four organs with the base average error rate below 0.02%. Trinity program was used for de novo assembling of clean data. 138, 888 unigenes and 207, 390 transcripts were obtained. The length of unigenes ranged from 201 to 14, 793 bp and the N50 length was 1, 862. 71, 215 unigenes (51.28%) and 121, 508 Transcripts (58.59%) were annotated based on the public databases. These sequences were used for identifying the miRNAs.
Identification of known miRNAs in P. heterophylla
In order to identify the conserved miRNAs in P. heterophylla, the small RNA were submitted to the miRBase 21 (http://www.mirbase.org/). Sequence with high identify to other species were considered as known miRNA. A total of 163 miRNAs belonging to 66 families were identified (Additional file 1). Among these families, the miRNA169 has 13 members distinguished based on nucleotide differences and followed by miRNA166 (12 members each); miRNA156 (11 members each); miRNA395 (9 members each); miRNA160 (8 members each) and miRNA399 (8 members each). The remaining miRNA families comprised 1 to 6 members each (Fig 2).
Of the 66 conserved miRNA families, miRNA396, miRNA159, miRNA319 and miRNA166 showed high abundance with total TPM >100 000 (Additional file 1). The highest read abundance (7 202 956 normalized reads) was detected for miR396. However, 37 miRNA families with total TPM less than 10 were found, such as miR1871, miR829, miR828 and miR2873. These results suggested variations in expression levels of miRNAs and their highly conserved function.
Identification of novel candidate miRNAs in P. heterophylla
The non-conserved small RNA reads were mapped to unigenes in transcriptome of P. heterophylla and a total of 303 novel miRNA candidates were identified (Additional file 2). A standard stem-loop hairpin secondary structure (SS) was found in all miRNA precursors. The secondary structure of precursors of novel miRNAs were predicted by using Mfold web server [31] and a few were showed in Figure 3. These miRNA precursors had folding free energies ranging from -48 to -18 kcal/mol (average, -27.28 kcal/mol). The length of predicted precursor of these novel miRNAs were 58-101 nt. The MFEI values of all miRNAs were above 0.85 and these sequences were considered as miRNA[32]. Of the 303 novel miRNA, 4 novel miRNAs had high expression level with TPM >1 000. The highest read abundance (2 662 normalized reads) was detected for Phn-m0190. 73 novel miRNA with total TPM less than 10 were found, such as Phn-m0032, Phn-m0045, Phn-m0050, Phn-m0058 and Phn-m0076.
Prediction of miRNA target genes
To explore the biological functions of the identified miRNAs, Target Finder program were used for searching the complimentary mRNA sequences from transcriptome sequence data of P. heterophylla above. Among the conserved targets, numbers of transcripts coding for transcription factors were identified (Additional file 3), such as GAMYB-like (targeted by miR159), APETALA2 and bHLH (target of miR172), PHAVOLUTA-like HD-ZIPIII protein (target of miR166), nuclear transcription factors – YA3, YA5, Zinc finger CCCH domain-containing protein and WRKY (target of miR169), NAC domain containing protein (targeted by miR164), Squamosa promoter-binding protein and MADS (regulated by miR156), TCP 4(targeted by miR319), Auxin response factor (regulated by miR160). This result suggested that they may participate in post-transcriptional regulation and transcriptional networks.
For novel miRNAs, 305 potential target genes involving with protein kinase were identified (Additional file 4), such as serine/threonine-protein kinase-like protein CCR4 (TRINITY_DN90887_c1_g2), L-type lectin-domain containing receptor kinase S.7 (TRINITY_DN87850_c4_g1), LRR receptor-like serine/threonine-protein kinase (TRINITY_DN93055_c2_g2 and TRINITY_DN93285_c2_g1) and plastidial pyruvate kinase 4 (TRINITY_DN88253_c0_g1). 162 potential target unigenes coding transcription factors were also observed, including MYB (44 unigenes), bHLH (34 unigenes), WRKY (15 unigenes) and MADS (12 unigenes). In addition, 34 potential target genes were founded to be involving with the synthesis and biotransformation of NADPH. These results indicated that novel miRNAs in P. heterophylla may participate in several diverse physiological and metabolic processes, including the regulation of plant metabolism, transport, cell growth and maintenance, and stress responses.
In addition, multiple potential regulation patterns were discovered in the same unigenes (Additional file 5). 8 binding site of Phn-m051-5p and Phn-m052-5p were identified in the mRNA sequence of ubiquitin-conjugating enzyme E2 24 (TRINITY_DN90017_c1_g1). transmembrane protein 97 (TRINITY_DN86988_c0_g5) have 11 potential regulation patterns target by miRNA399 families, Phn-m0052-3p, Phn-m0051-5p, Phn-m0052-5p and Phn-m0220-5p. Three HD-ZIPIII protein, ATHB8 (TRINITY_DN91558_c0_g5), ATHB15 (TRINITY_DN91558_c0_g3) and REVOLUTA (TRINITY_DN91558_c0_g7), have 9 potential regulation patterns target by miRNA166 families, Phn-m0140-3p and Phn-m0141-3p. Furthermore, several unigenes were only targeted by novel miRNA, such as multicopper oxidase LPR2 (TRINITY_DN90101_c3_g6), L-ascorbate oxidase homolog (TRINITY_DN92548_c7_g1), ribosomal RNA-processing protein 7 homolog A (TRINITY_DN90150_c1_g2) and pyrrolidone-carboxylate peptidase (TRINITY_DN87378_c4_g1). To summarize, the novel miRNA in P. heterophylla show an extreme important role in regulation of growth and development.
GO and KEGG analysis
GO analysis was performed to explore the potential regulation of predicted target transcripts of miRNAs in P. heterophylla. Based on their functional annotations, the target genes were into 122 categories (Additional file 6). In particular, 6 categories have higher enrichment, such as protein kinase activity, protein serine/threonine kinase activity, phosphotransferase activity, alcohol group as acceptor, calcium ion binding, kinase activity and response to acid chemical (Fig 4 ). KEGG pathway enrichment analysis showed that the target genes of the miRNAs were mainly involved in 8 pathways (P≤0.05): Plant-pathogen interaction, MAPK signaling pathway-plant, plant hormone signal transduction, ubiquitin mediated proteolysis, phosphatidylinositol signaling system, phenylpropanoid biosynthesis, inositol phosphate metabolism and diterpenoid biosynthesis (Additional file 7 and Fig 5). These results suggested that the miRNAs in P. heterophylla may play an important role in regulating cell signal transduction.
Role of miRNAs in ground and underground parts
For a better understanding of miRNAs in morphogenesis of organs, some tissue-special miRNAs were analyzed by combining with the transcriptome (Additional file 8). 6 conserved miRNAs was found to be most abundant in ground parts and could rarely be detected in underground parts, including 2 member of miR172 family (miR172e-5p and miR172h-5p), miR156a−3p, miR167g-3p, miR397b and miR398a-3p. 9 novel miRNA showed special expression in ground parts, including Phn−m021−5p, Phn−m0028−3p, Phn−m0084−3p, Phn−m0121−5p, Phn−m0147−3p, Phn−m0170−3p, Phn−m0173−5p, Phn−m0292−5p and Phn−m0246−5p (Fig 6). 115 potential target genes were identified for 18 miRNAs (Additional file 9), including 33 coding unigenes and 82 non-coding unigenes. 16 unigenes target by Phn−m0170−3p and 6 genes target by Phn−m0173−5p were identified. CRK7 (TRINITY_DN92623_c1_g1), PUB6 (TRINITY_DN88470_c4_g3) and BBX19 (TRINITY_DN86966_c1_g1) showed significantly negative correlation with the expression of Phn−m0173−5p. SODC (TRINITY_DN88258_c4_g4) was negatively regulated by miR398a-3p.
12 conserved miRNAs showed abundant expression in underground parts (root bark and xylem), while be restricted to ground parts (stem and leaf), including 3 member of miR164 family (miR164b, miR164c and miR164e), 4 member of miR169 family (miR169a, miR169a−3p, miR169j and miR169y), 2 member of miR1863 family (miR1863a and miR1863b), miR156b-3p, miR172c and miR6470. 13 novel miRNAs also high expression in underground parts (root bark and xylem), including Phn−m0008−3p, Phn−m0274−3p, Phn−m0103−5p, Phn−m0080−5p, Phn−m0114−5p, Phn−m0046−5p, Phn−m0238−3p, Phn−m0151−5p, Phn−m0167−5p, Phn−m0286−5p, Phn−m0035−5p, Phn−m0046−5p and Phn−m0113−3p. 107 potential target genes were identified for 18 miRNAs, including 37 coding mRNA and non-coding mRNA. It is founded that the expression profiles of two conserved miRNA, miRNA169a and miRNA169a, and Phn-m0046-5p were negatively correlated with NF-YA5 (TRINITY_DN85133_c0_g1) and PP2C08 (TRINITY_DN86665_c2_g1) respectively, which is in accordance with the gene-silencing function. Additionally, the expression of several non-coding unigenes showed significantly positive and negative correlation with these miRNAs.
Role of miRNAs in development of leaf
Interestingly, lot of miRNAs showed special expression profile in a certain organs. 9 conserved miRNAs and 2 novel miRNAs and could rarely be detected in leaf (Fig 7 and Additional file 10), including 5 member of miR319 family (miR319a, miR319a-3p.2-3p, miR319b-5p, miR319d-5p and miR319i), miR164f, miR171d, miR390b, miR394a-5p, Phn−m0151−5p and Phn−m0113−3p. 36 potential target genes were identified for 11 novel miRNAs, including 13 coding mRNA and 23 non-coding mRNA. In particular, the TCP (TRINITY_DN90304_c1_g3, TRINITY_DN89893_c1_g2 and TRINITY_DN89893_c1_g1) show significantly negative correlation with 3 miRNA319 families, miR319a, miR319a-3p.2-3p and miR319i. Interestingly, 3 non-coding mRNA have extremely positive correlation with the miRNA319 family. These results suggested that miR156a, miR159 family and miR319 family play an important role in development of leaf.
However, 19 conserved miRNAs and 2 novel miRNAs show abundant expression in leaf, including 3 members of miR156a (miR156a, miR156a−3p and miR156a−5p), 3 members of miR159 (miR159a, miR159d and miR159c), 2 members of miR167 (miR167h−5p and miR167d), 2 members of miR168(miR168b−5p and miR168a), 4 members of miR395 (miR395a, miR395c, miR395h and miR395k), 2 members of miR396 (miR396a−3p and miR396b−3p), miR171c−5p, miR172b, miR2111f, Phn−m0026−3p and Phn−m0222−5p. 88 potential target genes were identified for 21 novel miRNAs, including 29 coding mRNA and 59 non-coding mRNA. miR156a targeted 10 unigenes and GAMYB (TRINITY_DN91612_c2_g3 and TRINITY_DN89557_c1_g5) may regulated by miR159 families. Though none of miRNA showed negative correlation with the target unigenes, some unigenes have positive correlation with miRNAs, such as GAMYB (target by miR159a, miR159c and miR159d), EAF6 (target by miRNA156a) and receptor-like protein kinase At2g23200 (target by miR172b).
Role of miRNAs in development of stem
3 conserved miRNAs and 12 novel miRNAs have abundance transcripts in stem (Fig 8 and Additional file 11), including miR828b−5p, miR169d, miR156e, Phn−m0140−3p, Phn−m0208−5p, Phn−m0098−3p, Phn−m0282−5p, Phn−m0082−3p, Phn−m0291−5p, Phn−m0077−5p, Phn−m0252−3p, Phn−m0038−5p, Phn−m0295−3p, Phn−m0044−5p, Phn−m0024−3p. 239 potential target genes were identified for 15 miRNAs, including 75 coding mRNA and 164 non-coding mRNA. 36 coding mRNA target by Phn-m0082-3p, 10 coding mRNA target by Phn-m0077-5p and 7 coding mRNA target by Phn-m0024-3p. The expression of snakin-1 (TRINITY_DN90759_c4_g6) and PUB3 (TRINITY_DN90343_c1_g6) has negative correlation with Phn-m0082-3p and Phn-m0024-3p, respectively. However, miR156l-5p and Phn−m0224−3p showed low expression in stem whereas relatively high expression was observed in underground parts and leaf. The expression abundance of miR156l-5p showed negative correlation with SPL5 (TRINITY_DN90116_c0_g1).
Role of miRNAs in development of tuberous root
To explore the function of miRNAs in development of tuberous root, tissue-special miRNAs in xylem and bark were identified. 12 novel miRNAs were especially up-regulated in root xylem (Fig 9 and Additional file 12), including Phn−m0262−3p, Phn−m0217−3p, Phn−m0249−3p, Phn−m0261−5p, Phn−m0101−5p, Phn−m0198−3p, Phn−m0279−3p, Phn−m0274−5p, Phn−m0100−5p, Phn−m0184−3p, Phn−m0092−5p, Phn−m0214−5p. 96 potential target genes were identified for 12 novel miRNAs, including 36 coding mRNA and 60 non-coding mRNA. Though 24 coding unigenes was targeted by Phn−m0101−5p and 6 coding unigenes was targeted by Phn−m0279−3p, none of miRNA showed negative correlation with the target unigenes.
However, 16 conserved miRNAs and 5 novel miRNAs specially showed low expression in root xylem, including 7 member of miR166 family (miR166a−3p, miR166b-5p, miR166g-5p, miR166i-3p, miR166k, miR166k−3p and miR166n), 2 member of miR160 family (miR160b-3p and miR160d-5p), 2 member of miR171 family (miR171c-3p and miR171e-3p), 2 member of miR399 family (miR399a and miR399h), miR159b−5p, miR165a-3p, miR393a-5p, Phn−m0225−3p, Phn−m0162−5p, Phn−m0086−5p, Phn−m0042−5p and Phn−m0022−5p. 77 potential target genes were identified for 21 miRNAs, including 38 coding mRNA and 39 non-coding mRNA. Three HD-ZIP III proteins were extreme negatively regulated by miR166 family. It is seem that SCL6, BIR2, ARFR and WDL5, were negatively regulated by miR171 family, Phn-m0162-5p, miR160d-5p and miR166g-5p in xylem of tuberous root, respectively. These results suggested that these miRNAs play an important role in in development of tuberous root.
In addition, 7 conserved miRNAs were especially up-regulated in root bark (Fig 10 and Additional file 13), including 3 member of miR160 family (miR160f−5p, miR160g and miR160h), 2 member of miR390 family (miR390-3p and miR390d-3p), miR169w, miR171a. However, none of miRNA showed negative correlation with its target gene.
Role of miRNAs in hormone signal transduction
Hormone plays a curial role in the morphogenesis of plant organs. In this study, we analyzed regulation function of miRNA on hormone signal transduction by searching the target genes of miRNAs. The results showed that a large number of genes in most of the hormone signal transduction pathways were targeted by identified miRNAs in P. heterophylla, mainly involving with signal transduction of Auxin, Zeatin, Abscisic acid and Jasmonic acid (Fig 11 andAdditional file 14). Five genes in Auxin signal transduction pathway, TIR1, Aux/IAA and ARF, may heavily be regulated by miRNAs in P. heterophylla. SAUR is a multiple family and targeted by 10 miRNAs, including 2 known miRNAs and 8 novel miRNAs. In CTK signaling transduction, there were several binding sites of miRNAs in the sequences of AHP, B-ARR and A-ARR. Interestingly, PYL, PP2C and SnRK2 in ABA signal transduction pathway may strongly be regulated by miRNAs, mainly composed of novel miRNAs. We also founded that few novel miRNAs may regulate the expression of COI1 and JAZ in jasmonic acid signaling transduction pathway. In addition, gibberellin receptor GID2 in gibberellin signal transduction was targeted by Phn-m0174-5p. BSK in Brassinosteroid signaling transduction pathway was targeted by 3 known miRNAs and 5 novel miRNAs. A few miRNAs were identified involving with the ethylene signal transduction, such as MPK6 (targeted by Phn-m0087-5p, Phn-m0154-3p and Phn-m0281-5p), EBF1/2 (targeted by Phn-m0033-5p), and EIN3 (targeted by mir166e-5p and Phn-m0122-3p).
Role of miRNAs in interaction between plant and entophytes
For a better understanding of the role of miRNAs in interaction between plant and endophytes, mRNA sequences of 4 microorganisms (Fusarium oxysporum, Leptosphaeria maculans, Microbotryum violaceum, Rhodosporidium toruloides) were downloaded from Genbank and used for searching the target genes of miRNAs of P. heterophylla (Fig 12 andAdditional file 15). 39 genes of from F. oxysporum were targeted by 30 miRNAs of P. heterophylla, including 15 known miRNAs and 15 novel miRNAs. 6 genes were targeted by miR156h, followed by miR414(4 genes) and miR171a(3 genes)。Furthermore, some important genes of F. oxysporum may be regulated by miRNAs of P. heterophylla, such as pectin lyase F (target by Phn-m0104-3p), 30S ribosomal protein S17-likeprotein (miR171e-3p) and myo-inositol transport protein ITR1 (target by miR171a). Similar results were observed in L. maculans,M. violaceum,R. toruloides。Interestingly, there were 37, 43, 27 genes in L. maculans,M. violaceum,R. toruloides targeted by miR414, respectively. This result suggested that miR414 may play an important role in interaction between P. heterophylla and endophytes.
miRNA-induced cleavage of predicted targets
Cleavage site of predicted targets of a conserved miRNA and a novel miRNA in
- heterophylla were further verified by 5′ RLM-RACE (Fig 13). The predicted targets of miR166k and Phn-m0121-5p were HD-ZIP REVOLUTA and delta(24)-sterol reductase, respectively. All the fragment of 5’RLM-RACE were analysed on agarose gel, purified, cloned and sequenced. The cleavage site was found to be between 2–3 bases, from the 5′ end of the miRNA, for HD-ZIP REVOLUTA (TRINITY_DN91558_c0_g3). However, Delta (24)-sterol reductase (TRINITY_DN86257_c0_g3) was found to be cleaved between 5–6 from the 5′end of the miRNA binding site.