Global miRNA identification in Eucalyptus
Due to the unavailable access of Eucalyptus miRNAs, we first identified miRNAs in the E. grandis genome using the small RNA sequencing data of differentiated (stem) and dedifferentiated (tissue-culture induced callus) of two Eucalyptus cultivars (E. camaldulensis stem: A1; E. camaldulensis callus: A2; E. grandis x urophylla stem: B1; E. grandis x urophylla callus: B2). Initially, we obtained a total of 346.4 million raw reads (average 28.9 million reads) and 304.2 million clean reads (average 25.4 million reads) after data cleaning. Using the clean data, we identified 888 miRNA precursors (producing 1,067 mature miRNAs, Additional file 1), including 691 novel and 197 known miRNAs. Length distribution (Figure 1A) showed that 786 Eucalyptus mature miRNAs were in the length of 21 nt, which took 73.7% of the total miRNAs. Figure 1A also showed that the first nucleic acid of the miRNAs was enriched with uracil (U), which agrees with the base bias of plant miRNAs in miRBase [15]. Next, we analyzed the chromosome distribution of the Eucalyptus miRNAs. Figure 1B showed that chromosome Chr03 can produce the most miRNAs (142 miRNA precursors), followed by chromosomes Chr05 (129 miRNA precursors) and Chr08 (92 miRNA precursors). Based on the miRNA locations in the Eucalyptus genome, we identified 46 miRNA clusters whose members are located within 10 kb (Additional file 1) and have the potential of being transcribed together. We also found that 55, 52 and 782 miRNA precursors were located in the exon, intron and intergenic regions, respectively. The redundant one was egd-MIR169aj, which can be produced by two chromosomes – Chr03 and Chr09 (Additional file 1). These miRNAs were used as the reference for downstream analyses, including the miRNA expression profiling and target prediction.
miRNA expression profiles
We next profiled the miRNA expression in differentiated and dedifferentiated tissues of E. camaldulensis and E. grandis x urophylla. After lowly expressed miRNAs (average TPM < 5) were filtered, we identified 998 miRNAs in all samples (Figure 1C, Additional file 2), of which 500, 547, 505 and 620 distributed in A1, A2, B1 and B2, respectively. Venn diagram (Figure 1C) showed that 230 miRNAs were detected in all samples while 81, 79, 60 and 163 miRNAs were specifically identified in A1, A2, B1 and B2, respectively. The distribution of miRNAs in different expression levels revealed that more than 60% of the miRNAs were less than 5 TPM (Figure 1D) while 31~38 miRNAs were expressed more than 1,000 TPM in these samples. We showed in Table 1 the top 10 highly expressed miRNAs in these samples. It is interesting that the only the top 3 highly expressed miRNAs were shared by all four samples and that MIR166 family members were highly expressed in the dedifferentiated tissues of both E. camaldulensis and E. grandis x urophylla.
Differentially expressed miRNAs in E. camaldulensis
We next compared the miRNA profiles of differentiated and dedifferentiated tissues in E. camaldulensis. Compared to the differentiated tissue we first identified 143 up-regulated and 184 down-regulated miRNAs in the dedifferentiated tissue of E. camaldulensis (Figure 2A,Additional file 3). Among them, we found that 3 MIR167, 26 MIR169, 7 MIR3627 and 3 MIR159 members were up-regulated and that 11 MIR156, 4 MIR164, 9 MIR166, 8 MIR171, 7 MIR395 and 43 MIR399 were down-regulated. Target prediction showed that 310 of the differentially expressed miRNAs have the potential of targeting 1,401 genes in E. camaldulensis (Table 2). For example, 135 up-regulated miRNAs can target 819 genes including Eucgr.A00998 (transcription factor PIL1), Eucgr.B03766/Eucgr.G03183 (transcription factor GAMYB), Eucgr.C00823 (growth-regulating factor 1) and Eucgr.B03050 (AP2/ERF and B3 domain-containing transcription factor RAV1); and 175 down-regulated miRNAs in dedifferentiated tissue can target 608 genes including auxin response factors (e.g., Eucgr.D00264, Eucgr.F04380, Eucgr.G02838, Eucgr.J00923, Eucgr.K01240), Eucgr.K02897 (transcription factor MYB1R1), Eucgr.F02691 (ethylene-responsive transcription factor ABR1) and Eucgr.D02118 (translation initiation factor IF-2).
Next, we analyzed the potential pathways and biological processes involved by the differentially expressed miRNAs in the dedifferentiated tissue of E. camaldulensis. KEGG pathway enrichment (Figure 2B) showed that 167, 62, 61, 37, 10, 9 and 3 target genes were involved by “Ras signaling pathway” (ko04014), “Plant hormone signal transduction” (ko04075), “Longevity regulation pathway – multiple species” (ko04213), “Flavonoid biosynthesis” (ko00941), “Zeatin biosynthesis” (ko00908) and “Aflatoxin biosynthesis” (ko00254), respectively. Gene ontology annotation (Additional file 4) showed that the target genes of differentially expressed miRNAs in dedifferentiated tissue have the potential of regulation in the biological processes such as “GO:0009809~lignin biosynthetic process”, “GO:0009744~response to sucrose’, “GO:0009909~regulation of flower development” and “GO:0010075~regulation of meristem growth”.
Differentially expressed miRNAs in E. grandis x urophylla
In E. grandis x urophylla we identified 128 up-regulated and 215 down-regulated miRNAs in dedifferentiation process (Figure 2C, Additional file 3). The callus up-regulated miRNAs include 13 MIR169, 4 MIR3627, 3 MIR390 and 4 MIR482 while the callus down-regulated miRNAs include 21 MIR156, 3 MIR160, 11 MIR166, 44 MIR169, 12 MIR171 and 43 MIR399. We found that 314 differentially expressed miRNAs in the dedifferentiated tissue of E. grandis x urophylla were predicted to target 1,550 genes (Table 2). In details, 105 up-regulated miRNAs target 783 genes, including Eucgr.A00998 (transcription factor PIL1), Eucgr.C00360 (transcription factor bHLH69), Eucgr.C00823 (growth-regulating factor 1), Eucgr.C01912 (probable WRKY transcription factor 19), Eucgr.D00264 (auxin response factor 6) and Eucgr.J02212 (transcription factor MYB23), and 209 down-regulated miRNAs target 810 genes, such as auxin response factors (e.g., Eucgr.F04380, Eucgr.G02838, Eucgr.K01240, Eucgr.J00923, Eucgr.K03433, Eucgr.D00264), Eucgr.F03952 (factor of DNA methylation 1), Eucgr.H01993 (myb family transcription factor EFM), Eucgr.G01493 (transcription factor bHLH95), transcription factor GAMYB (Eucgr.B03766, Eucgr.G03183), Eucgr.K03562 (transcription factor MYB108), Eucgr.C00212 (transcription factor IBH1) and Eucgr.I02059 (truncated transcription factor CAULIFLOWER A). Functional analysis (Figure 2D) showed that the differentially expressed miRNAs in the callus tissue of E. grandis x urophylla were involved by pathways such as “Aflatoxin biosynthesis” (ko00254), “Brassinosteroid biosynthesis” (ko00905), “ABC transporters” (ko02010), “Spliceosome” (ko03040) and “Estrogen signaling pathway” (ko04915). GO analysis (Additional file 4) of the target genes revealed that “GO:0009809~lignin biosynthetic process”, “GO:0006952~defense response”, “GO:0006355~regulation of transcription, DNA-templated”, “GO:0009658~chloroplast organization” and “GO:0009867~jasmonic acid mediated signaling pathway” were significantly enriched by the differentially expressed miRNAs in dedifferentiation process of E. grandis x urophylla. Based on the functional analysis, we found that the biological processes and pathways regulated by the differentially expressed miRNAs in E. camaldulensis and E. grandis x urophylla were different, which might indicate the different embryogenetic potentials of these two Eucalyptus cultivars.
miRNAs related to somatic embryogenesis
To investigate the miRNAs related to the embryogenetic potential of Eucalyptus, we next compared the differentially expressed miRNAs in E. camaldulensis and E. grandis x urophylla. It showed that they shared 62 up-regulated and 117 down-regulated miRNAs (Figure 3A, Additional file 3), including miRNAs from families MIR156, MIR159, MIR160, MIR164, MIR166, MIR169, MIR171, MIR399 and MIR482. We also identified 81 up-regulated (e.g., miR159c-3p, miR167a-5p, miR397a-3p, miR397c-5p, miR397d-3p, miR397d-5p, N-miR1-5p and N-miR5-5p) and 67 down-regulated (e.g., miR482b-3p, N-miR3-3p, miR156a-3p, N-miR40-3p and N-miR18-5p) miRNAs specific to E. camaldulensis (Figure 3A). The heat map (Figure 3B) confirmed the dysregulation of these miRNAs specific to E. camaldulensis and showed the expression patterns across the differentiated and dedifferentiated tissues of E. camaldulensis and E. grandis x urophylla. In Figure 3C, we listed the expression levels of 6 miRNAs (e.g, miR399a-5p, N-miR83-5p, N-miR125-3p, miR395a-3p, N-miR40-3p and miR156a-3p) specifically dysregulated in E. camaldulensis. Functional analysis of the targets of E. camaldulensis specifically dysregulated miRNAs revealed that they were involved in the pathways including “MAPK signaling pathway” (ko04010), “Longevity regulating pathway - multiple species” (ko04213), “Estrogen signaling pathway” (ko04915), “Plant hormone signal transduction” (ko04075) and “Brassinosteroid biosynthesis” (ko00905).
qRT-PCR
Next, we used tailed miRNA qRT-PCR strategy to validate the differential expression of seven randomly selected miRNAs in the differentiated and dedifferentiated tissues of E. camaldulensis and E. grandis x urophylla. The primer sequences for these miRNAs can be found in Additional file 5. Comparison of small RNA sequencing and qRT-PCR results is shown in Figure 4. Overall, 12 (85.7%) out of 14 events were agreed by both small RNA sequencing and qRT-PCR. Using deep sequencing, egd-miR482b-3p, egd-miR156a-3p and egd-N-miR1-5p were found with down-regulation only in DH201-2 while egd-miR390a-3p, egd-miR396d-3p and egd-N-miR6-3p were found with dysregulation only in GL9. The specific regulations of these miRNAs were confirmed qRT-PCR. High agreement of miRNA expression patterns in small RNA sequencing and qRT-PCR indicate that the miRNAs identified in this study might be associated with the somatic embryogenesis and embryogenetic potential of Eucalyptus, which requires more experiments.