Genome-wide identification of microRNAs associated with the somatic embryogenesis in Eucalyptus

Abstract

Background: MicroRNAs (miRNAs) are a class of small noncoding RNAs with 18-24 nucleotides in length and function in many biological processes in plant. Although Eucalyptus trees are widely planted across the world, our understanding of the miRNA regulation in the somatic embryogenesis of Eucalyptus is still poor. Here we reported for the first time the miRNA profiles of differentiated and dedifferentiated tissues of two Eucalyptus cultivars and identified miRNAs involved in the somatic embryogenesis of Eucalyptus.

Results: Stem and tissue-culture induced callus were obtained from the subculture seedlings of E. camaldulensis and E. grandis x urophylla, and were used as differentiated and dedifferentiated samples, respectively. We generated 346.4 million reads for 12 samples (n=3) and identified 888 miRNA precursors (197 known and 691 novel) which can produce 1,067 mature miRNAs. These miRNAs were mainly distributed in chromosomes Chr03, Chr05 and Chr08, and can produce 46 miRNA clusters. In these samples we detected 998 miRNAs with TPM (transcripts per million reads) > 5 and found that highly expressed miRNAs varied across samples. We identified 327 and 343 differentially expressed miRNAs in the dedifferentiation process of E. camaldulensis and E. grandis x urophylla, respectively. Dysregulated miRNAs shared by the two cultivars might be involved in the development of embryonic callus of Eucalyptus, such as 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. Target prediction and functional analysis showed they might be involved in longevity regulating and plant hormone signal transduction pathways. Then, the expression patterns of these miRNAs were confirmed by qRT-PCR.  

Conclusions: This is the first time to study the miRNAs profiles in the dedifferentiation process of Eucalyptus and it will provide a valuable resource for future studies. More importantly, our findings will improve our understanding of miRNA regulation and molecular mechanisms during the somatic embryogenesis of Eucalyptus, and the output of this study will benefit the Eucalyptus breeding program.

Background

Eucalyptus, a highly diverse genus of the Myrtaceae family, is the most widely planted hardwood due to its increasing importance for fiber and energy in the word. Their remarkable adaptability coupled with fast growth and superior wood properties has driven their rapid adoption for plantation forestry in South America, Africa, Asia, Spain, Portugal, Middle Eastern countries and North America [1]. This genus consists of about 700 species, of which Eucalyptus globulus (E. globulus) and Eucalyptus grandis (E. grandis) are currently the most extensively planted species [2]. Natural regeneration of Eucalyptus mainly relies on seed and the breeding process is generally slow because of the length of the juvenile phase before flowering. Although natural hybridization of Eucalyptus has often been reported to gain some characteristics (e.g., cold resistance), artificial hybridization is still a means of producing a limited number of seeds and the selected hybrids must then be multiplied by vegetative propagation.

Several means of vegetative propagation for Eucalyptus have been reported, such as air layering, grafting, stem cutting and tissue culture [1]. During the in vitro tissue culture-induced organogenesis or somatic embryogenesis, cells of explants are activated to dedifferentiate and re-enter the cell cycle and in some cases shoots or roots are also regenerated after callus formation [3, 4]. The first step of somatic embryogenesis is the callus induction, which includes a number of characteristic events: dedifferentiation of cells, activation of cell division and reprogramming of cell physiology, metabolism and gene expression patterns [5]. The second step of somatic embryogenesis is the plant regeneration, which is related to the ability to regenerate plant (also called as embryogenic potential) [6]. Some studies have been demonstrated to study the molecular mechanisms of somatic embryogenesis in other plants, such as maize [7], rice [8], banana [9] and Catalpa bungei [10]. However, large is not clear about the molecular mechanisms in the somatic embryogenesis of Eucalyptus.

MicroRNAs (miRNAs) are a class of small (~21 nucleotide) noncoding RNAs that function in regulating gene expression at post-transcription level [2, 11]. They have been reported to function in multiple biological processes in plant, such as development, flowering, architecture and response to biotic and abiotic stresses [12]. Some studies have been demonstrated to identify miRNAs in Eucalyptus. For example, Levy found no mutual relationship between alternations in miR156 and miR172 expression in the adventitious root induction [13]; Pappas reported the variable patterns of conservation and diversity with miRNAs across Myrtaceae species [2]; and Lin identified 386 novel miRNAs in E. grandis and predicted their target genes [14]. These studies provide a basis of understanding the miRNAs regulatory network in Eucalyptus. However, little is known about the miRNA regulation in the somatic embryogenesis of Eucalyptus.

Small RNA sequencing has been used to identify known and novel miRNAs in animals and plants. In Eucalyptus, it has been used to characterize the miRNAs and its targets in some studies [2, 14]. In this study, we aimed to identify the miRNA profiles of differentiated (stem) and dedifferentiated (callus) tissues in Eucalyptus camaldulensis Dehnh. (E. camaldulensis, high embryogenic potential) and E. grandis x urophylla (low embryogenic potential). Also, we aimed to characterize differentially expressed miRNAs between differentiated and dedifferentiated tissues, which might be associated with the embryogenic potential of Eucalyptus. Our study will provide a valuable resource for future studies and improve our understanding towards the molecular mechanisms during the somatic embryogenesis and propagation development of Eucalyptus. Our results will benefit the breeding program of Eucalyptus in this field.

Results

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.

Discussion

Although some studies have been demonstrated to study the Eucalyptus miRNAs [2, 13, 14], the miRNA identities and sequences of Eucalyptus are still missing in miRbase (v22). In this study, we conducted the small RNA sequencing and identified a total of 888 miRNA precursors and 1,067 mature miRNAs in the Eucalyptus genome (Additional file 1). Among them, 197 precursors were aligned to other known miRNAs while 691 could be specific to the Eucalyptus. The miRNAs identified in this study can provide a valuable source to study the regulation mechanisms of noncoding RNAs in the somatic embryogenesis of Eucalyptus. Comparison of miRNAs identified this study (Figure 1C) revealed 104 miRNAs specific to the callus tissues of Eucalyptus, including six highly expressed miRNAs (e.g., egd-miR482a-3p, egd-miR482a-5p, egd-miR397a-5p, egd-miR397b-5p, egd-miR159b-3p, egd-miR159a-3p) (Table 1). While 79 miRNAs specifically identified in the callus tissue of DH201-2 (Figure 1C, Additional file 2) might be related to the high potential of embryogenesis.

Among the top 10 highly expressed miRNAs (Table 1), miR482a, miR397a/b, and miR159a/b were shared by the dedifferentiated tissues of both E. camaldulensis and E. grandis x urophylla. In the callus-derived protoplast, both miR482a and miR482b were up-regulated [16], and they were supposed to regulate the disease resistance pathways and the expression of plant protoplast totipotency [17]. miR482 family has been proved to target mRNAs for NBS-LRR disease resistance proteins using transient expression in Nicotiana benthamiana [18]. In this study, no NBS-LRR gene was predicted to be regulated by miR482a, probably due to the limited annotation of Eucalyptus genes. However, we found that miR482a was predicted to target 30 genes encoding disease resistance proteins, such as Eucgr.A00735 (putative disease resistance protein At4g11170), Eucgr.A02363 (RGA1), Eucgr.C02142 (RPS6), Eucgr.H03707 (RML1B), Eucgr.H03723 (TAO1) and 14 genes encoding TMV resistance protein N. Interestingly, MIR397 has been reported to be particularly expressed in undifferentiated tissue and to be overexpressed in undifferentiated rice callus [19, 20], which coincide with our results. MIR397 was reported to target a family of copper-containing oxidase enzymes – laccase, which is proposed to function in the formation of lignin, an integral part of the secondary cell walls of plants [19]. We found a total of 112 laccase genes regulated by MIR397 in Eucalyptus. Also, miR397b-3p was up-regulated in the callus of both E. camaldulensis and E. grandis x urophylla, and miR397a-3p was up-regulated only in the dedifferentiated tissue of E. camaldulensis (Additional file 3). Like MIR397, MIR159 has been studied to regulate the development of callus in multiple plants, such as Dimocarpus longan [21, 22], Moringa oleifera [23] and Larix leptolepis [24]. miR159 has been shown to be induced by ABA and to regulate MYB33 during the seed germination of Arabidopsis [25]. Although in the present study miR159 was not differentially expressed in the dedifferentiated tissue, its high expression indicates a potential role in the dedifferentiation and the callus development.

In addition, we identified 81 up-regulated and 67 down-regulated miRNAs specifically in the dedifferentiated tissue of E. camaldulensis compared to E. grandis x urophylla (Figure 3A). Among them, some miRNAs were detected with relative high expression in the stem or callus, such as miR156a-3p, miR167a-5p, miR482b-3p, miR3488-3p and miR535a-3p (Additional file 2). miR156a and miR482b were down-regulated in the callus of E. camaldulensis (Additional file 3). In Arabidopsis, the repression of miR156a by MYB33, a target of miR159, has the capacity of accelerate the vegetative development through its target SPL9 [26]. In E. grandis and in E. brachyphylla the expression of miR156 was shown to have no relationship with the rooting ability [13]. Thus, the relationship between miR156 and somatic embryogenesis in E. grandis requires further experiments to be explored. The expression of miR167a was found to be up-regulated during the development of rice callus and its target ARF8 was found to be down-regulated [27]. miR535-3p is the passenger miRNA of miR535, which was relevant from their expression (Additional file 2). miR535-5p was down-regulated in both E. camaldulensis and E. grandis x urophylla (Additional file 3). In rice, SPL14, the target of miR535, remains at a low level during the vegetative growth and can strengthen the ability of disease resistance in plant [28]. miR3488 was identified due to its similarity to an animal miRNA. More experiments are required to characterize its presence in Eucalyptus and the association with the somatic embryogenesis.

We also identified a number of novel miRNAs for Eucalyptus (Additional file 1) and some of them were highly expressed in the differentiated and dedifferentiated tissues, such as egd-N-miR1-3p, egd-N-miR2-3p, egd-N-miR4-3p, egd-N-miR5-5p and egd-N-miR7-5p (Table 1, Additional file 2). The abundance of these miRNAs support that they might be specific to the stem or the callus of Eucalyptus. Further, some of them were found to be differentially expressed in this study, such as egd-N-miR1-3p, egd-N-miR2-3p, and egd-N-miR5-5p (Additional file 3). Among them, 99 Eucalyptus novel miRNAs were identified with dysregulation specific in E. camaldulensis, such as egd-N-miR3-3p and egd-N-miR-5p (Additional file 3). Target prediction identified Eucgr.B01315 ((-)-germacrene D synthase) is the target of egd-N-miR3-3p while egd-N-miR5-5p can regulate the expression of Eucgr.C03404 (myosin-17) and Eucgr.E03882 (microtubule-associated protein 70-5). Although it is not clear about their functions in the somatic embryogenesis, our results provide a valuable resource for future studies in Eucalyptus.

Conclusions

In this study, we analyzed the small RNA profiles of differentiated and dedifferentiated tissues of two Eucalyptus cultivars – E. camaldulensis and E. grandis x urophylla. In total, we identified 691 novel and 197 known miRNAs for Eucalyptus, most of which distributed in chromosome Chr03 and Chr05. In the stem and callus we detected 998 miRNAs with TPM > 5, of which 230 shared by all samples. Further, we identified 143 up-regulated and 184 down-regulated miRNAs in the dedifferentiation process of E. camaldulensis, and 128 up-regulated and 215 down-regulated miRNAs in E. grandis x urophylla. Comparison of the differentially expressed miRNAs in dedifferentiated tissues revealed 62 up-regulated and 117 down-regulated miRNAs shared by both E. camaldulensis and E. grandis x urophylla, which might be related to the somatic embryogenesis of Eucalyptus. Further, the up-regulation of 81 miRNAs and down-regulation of 67 miRNAs were specifically identified in E. camaldulensis. These miRNAs were shown to function in the regulation of pathways like “Longevity regulating pathway - multiple species” and “Plant hormone signal transduction”. The expression of miRNAs was confirmed by qRT-PCR and their relationship with somatic embryogenesis was agreed with previous studies in other plants. This is the first time to study the association of miRNAs and somatic embryogenesis in Eucalyptus. The output of this study will provide a valuable resource for future studies in Eucalyptus. More importantly, our findings will benefit the molecular research in this field and the breeding program of Eucalyptus.

Methods

Plants and tissue culture

The plants of E. camaldulensis (high regenerative ability, voucher ID: c0009) and E. grandis x urophylla (low regenerative ability, voucher ID: j0017) were maintained in the experimental fields of Guangxi Forestry Research Institute. The second generation of in vitro tissue-culture induced seedlings of E. camaldulensis and E. grandis x urophylla were maintained on the MS medium supplemented with 20mg/L Ca(NO₃)₂, 0.5 mg/L 6-BA and 0.1 mg/L IAA until 2 to 3 cm height for this project. The second to the third stems from the stem tip of the seedlings were obtained and cut into 0.3 ~ 0.5 cm segments. About 60 segments of each Eucalyptus species were then transferred onto the induction MS medium (supplemented with 20mg/L Ca(NO₃)₂, 1 mg/L KT and 0.5 mg/L 2,4-D) and maintained in dark at 28±2 °C for 10 days. Successful callus samples were used as dedifferentiated tissues. The induction experiment was replicated three times.

RNA isolation and quality control

Total RNA was isolated from the plant tissues (100 mg, in triplicates) using the TRIzol reagent (Invitrogen) according to the manufacturer’s protocol, as previously described [11]. The quantity and quality of total RNA was evaluated by Agilent Bioanalyzer 2100 (Agilent Technologies).

Library preparation of small RNAs

Total RNA (1 µg) of each sample was used to construct the small RNA libraries. Initially, total RNA samples were fractionated on a 15% urea-PAGE gel electrophoresis and a band corresponding to small RNAs (18 – 30 bp) were excised. After small RNAs extracted by centrifugation, they were ligated with the adenylated 3′ adapter. Then, RT primer with barcode was used to anneal the 3’ adenylated adapter, followed by the ligation of 5’ adapter and the reverse transcript reaction. After first strand cDNA synthesis, we amplified the product by 15 cycles and carried out another size selection (103 – 115 bp) from the gel. After gel purification, the PCR product was quantified by Qubit (Invitrogen) and a single strand DNA circle (ssDNA circle), which gave the final miRNA library, was made by pooling multiple samples together.

Small RNA deep sequencing

We first generated the DNA nanoballs (DNBs) with ssDNA circle by rolling circle replication (RCR) to enlarge the fluorescent signals at the sequencing process, as previously described [29]. Then, the DNBs were loaded into the patterned nanoarrays and single-end read of 50 bp were read through on the BGISEQ-500 platform.

Data analysis and miRNA identification

The raw reads were cleaned by SOAPnuke, as previously described [30]. In order to identify miRNAs in the flooded Eucalyptus genome, we merged all the clean reads into one file and filtered reads less than 50 copies. The clean data was aligned to the Eucalyptus genome (http://plantgenie.org) using SOAP2 [31] and potential miRNA precursors were predicted using MIREAP. Then, potential miRNAs were searched against all the plant miRNA sequences in miRBase (v22.1). The family information of Eucalyptus miRNAs was manually identified using their similarities to known miRNAs in other species. miRNA clusters were identified using a 10,000 bp window screening in the genome. miRNAs not aligned to other known miRNAs were considered as novel miRNAs.

Expression profile of small RNAs

Before miRNA expression was profiled, we annotated other kinds of small RNAs in each sample. First, rRNA, tRNA, scRNA, snRNA, snoRNA and srpRNA were annotated by mapping the clean reads to NCBI GenBank and Rfam databases. Then, mRNA fragments were identified by mapping the clean reads to the Eucalyptus mRNA sequences. Next, we aligned the clean reads of each sample to all the miRNA precursors identified before using BLAST (no mismatch). Reads aligned to each mature miRNA were counted. Normalization of miRNA expression was performed using the TPM (transcripts per million reads) method and total clean reads were used as background.

Differential expression analysis

We performed differential expression analysis of miRNAs using edgeR, as previously described [32]. Raw counts of all miRNAs were used as input for edgeR. Some cut-offs were applied to the expression levels and statistical values to select differentially expressed miRNAs, such as TPM > 5, log2 fold change (log2fc) > 1 or log2fc < -1, p-value < 0.05 and false discovery rate (FDR) <0.1.

miRNA target prediction

miRNA target prediction was performed using the psRobot software [33].

Functional analysis

We first annotated the Eucalyptus genes by mapping them to the KEGG and GO databases. Then, enriched GO terms and pathways were identified using two statistical values – p-value (calculated by Fisher’s exact test, < 0.05) and q-value (calculated by the R package ‘qvalue’, < 0.05), as previously described [34].

qRT-PCR

Quantitative real-time PCR (qRT-PCR) was used to confirm the expression levels of miRNAs identified in this study. We randomly selected 7 miRNAs and used the actin gene as internal control.

Declarations

Ethics approval and consent to participate

No specific permits were required for the described field studies. The location is not privately-owned or protected in any way, and the field studies did not involve endangered or protected species.

Consent for publication

Not applicable.

Availability of data and materials

The raw sequencing data can be accessed from the NCBI Sequence Read Archive (SRA) platform (http://trace.ncbi.nlm.nih.gov/Traces/sra/) under the accession number SRA1074080.

Competing interests

The authors declare that there is no conflict of interest.

Funding

This work was financially supported by the Guangxi Science and Technology Program (GuiKeAD18281086), the Guangxi Forestry Science and Technology Project (2016[11], 2014[033]), the National Natural Science Foundation of China (31400522) and the Department of Human Resources and Social Security of Guangxi Zhuang Autonomous Region (GuiCaiSheHan2018112). The authors declare that the funding bodies have no role in the research design, the data collection and analysis, and the manuscript preparation.

Authors' contributions

BC, ZQ and JL conceived and designed the experiments; YZ, YX and XZ performed the experiments; ZQ, JL and HL analyzed the data; ZQ and JL wrote the manuscript; BC revised the manuscript. All the authors have read and approved the final version of manuscript.

Acknowledgements

We acknowledge Mr. Dajie Zhou for his involvement of the small RNA sequencing in BGI-Shenzhen.

Abbreviations

Eucalyptus globulus: E. globulus; Eucalyptus grandis: E. grandis; microRNA: miRNA; TPM: transcripts per million reads.

References

1. Durand-Cresswell R, Boulay M, Franclet A: Vegetative Propagation of Eucalyptus. In: Tissue Culture in Forestry. Edited by Bonga JM, Durzan DJ. Dordrecht: Springer Netherlands; 1982: 150-181.
2. Pappas Mde C, Pappas GJ, Jr., Grattapaglia D: Genome-wide discovery and validation of Eucalyptus small RNAs reveals variable patterns of conservation and diversity across species of Myrtaceae. BMC Genomics 2015, 16:1113.
3. Ikeuchi M, Sugimoto K, Iwase A: Plant callus: mechanisms of induction and repression. Plant Cell 2013, 25(9):3159-3173.
4. Pulianmackal AJ, Kareem AV, Durgaprasad K, Trivedi ZB, Prasad K: Competence and regulatory interactions during regeneration in plants. Front Plant Sci 2014, 5:142.
5. Zimmerman JL: Somatic Embryogenesis: A Model for Early Development in Higher Plants. Plant Cell 1993, 5(10):1411-1423.
6. Dinkova TD, Alejandri-Ramirez ND: MicroRNA Expression and Regulation During Plant Somatic Embryogenesis. In: Epigenetics in Plants of Agronomic Importance: Fundamentals and Applications: Transcriptional Regulation and Chromatin Remodelling in Plants. Cham: Springer International Publishing; 2014: 111-123.
7. Lopez-Ruiz BA, Juarez-Gonzalez VT, Chavez-Hernandez EC, Dinkova TD: MicroRNA Expression and Regulation During Maize Somatic Embryogenesis. Methods Mol Biol 2018, 1815:397-410.
8. Zhang Z, Zhao H, Li W, Wu J, Zhou Z, Zhou F, Chen H, Lin Y: Genome-wide association study of callus induction variation to explore the callus formation mechanism of rice. J Integr Plant Biol 2019, 61(11):1134-1150.
9. Kumaravel M, Uma S, Backiyarani S, Saraswathi MS, Vaganan MM, Muthusamy M, Sajith KP: Differential proteome analysis during early somatic embryogenesis in Musa spp. AAA cv. Grand Naine. Plant Cell Rep 2017, 36(1):163-178.
10. Liu W, Wang C, Shen X, Liang H, Wang Y, He Z, Zhang D, Chen F: Comparative transcriptome analysis highlights the hormone effects on somatic embryogenesis in Catalpa bungei. Plant Reprod 2019, 32(2):141-151.
11. Ji H, Chen M, Greening DW, He W, Rai A, Zhang W, Simpson RJ: Deep Sequencing of RNA from Three Different Extracellular Vesicle (EV) Subtypes Released from the Human LIM1863 Colon Cancer Cell Line Uncovers Distinct Mirna-Enrichment Signatures. PLOS ONE 2014, 9(10):e110314.
12. Axtell MJ, Snyder JA, Bartel DP: Common functions for diverse small RNAs of land plants. Plant Cell 2007, 19(6):1750-1769.
13. Levy A, Szwerdszarf D, Abu-Abied M, Mordehaev I, Yaniv Y, Riov J, Arazi T, Sadot E: Profiling microRNAs in Eucalyptus grandis reveals no mutual relationship between alterations in miR156 and miR172 expression and adventitious root induction during development. BMC Genomics 2014, 15:524.
14. Lin Z, Li Q, Yin Q, Wang J, Zhang B, Gan S, Wu A-M: Identification of novel miRNAs and their target genes in Eucalyptus grandis. Tree Genetics & Genomes 2018, 14(4):60.
15. Griffiths-Jones S, Saini HK, van Dongen S, Enright AJ: miRBase: tools for microRNA genomics. Nucleic Acids Res 2008, 36(Database issue):D154-158.
16. Xu X, Xu X, Zhou Y, Zeng S, Kong W: Identification of protoplast-isolation responsive microRNAs in Citrus reticulata Blanco by high-throughput sequencing. PLoS One 2017, 12(8):e0183524.
17. Wu XM, Kou SJ, Liu YL, Fang YN, Xu Q, Guo WW: Genomewide analysis of small RNAs in nonembryogenic and embryogenic tissues of citrus: microRNA- and siRNA-mediated transcript cleavage involved in somatic embryogenesis. Plant Biotechnol J 2015, 13(3):383-394.
18. Shivaprasad PV, Chen HM, Patel K, Bond DM, Santos BA, Baulcombe DC: A microRNA superfamily regulates nucleotide binding site-leucine-rich repeats and other mRNAs. Plant Cell 2012, 24(3):859-874.
19. Luo YC, Zhou H, Li Y, Chen JY, Yang JH, Chen YQ, Qu LH: Rice embryogenic calli express a unique set of microRNAs, suggesting regulatory roles of microRNAs in plant post-embryogenic development. FEBS Lett 2006, 580(21):5111-5116.
20. Chen CJ, liu Q, Zhang YC, Qu LH, Chen YQ, Gautheret D: Genome-wide discovery and analysis of microRNAs and other small RNAs from rice embryogenic callus. RNA Biol 2011, 8(3):538-547.
21. Lin YL, Lai ZX: Evaluation of suitable reference genes for normalization of microRNA expression by real-time reverse transcription PCR analysis during longan somatic embryogenesis. Plant Physiol Biochem 2013, 66:20-25.
22. Xu X, Chen X, Chen Y, Zhang Q, Su L, Chen X, Chen Y, Zhang Z, Lin Y, Lai Z: Genome-wide identification of miRNAs and their targets during early somatic embryogenesis in Dimocarpus longan Lour. Sci Rep 2020, 10(1):4626.
23. Pirro S, Matic I, Guidi A, Zanella L, Gismondi A, Cicconi R, Bernardini R, Colizzi V, Canini A, Mattei M et al: Identification of microRNAs and relative target genes in Moringa oleifera leaf and callus. Sci Rep 2019, 9(1):15145.
24. Zhang S, Zhou J, Han S, Yang W, Li W, Wei H, Li X, Qi L: Four abiotic stress-induced miRNA families differentially regulated in the embryogenic and non-embryogenic callus tissues of Larix leptolepis. Biochem Biophys Res Commun 2010, 398(3):355-360.
25. Reyes JL, Chua NH: ABA induction of miR159 controls transcript levels of two MYB factors during Arabidopsis seed germination. Plant J 2007, 49(4):592-606.
26. Guo C, Xu Y, Shi M, Lai Y, Wu X, Wang H, Zhu Z, Poethig RS, Wu G: Repression of miR156 by miR159 Regulates the Timing of the Juvenile-to-Adult Transition in Arabidopsis. Plant Cell 2017, 29(6):1293-1304.
27. Jun C, Wang WG, Wang SH, Chen F: The study on Auxin-miR167-ARF8 signal pathway during the growth and development process of rice callus. Journal of Sichuan University 2013.
28. Agirre X, Jimenez-Velasco A, San Jose-Eneriz E, Garate L, Bandres E, Cordeu L, Aparicio O, Saez B, Navarro G, Vilas-Zornoza A et al: Down-regulation of hsa-miR-10a in chronic myeloid leukemia CD34+ cells increases USF2-mediated cell growth. Mol Cancer Res 2008, 6(12):1830-1840.
29. Drmanac R, Sparks AB, Callow MJ, Halpern AL, Burns NL, Kermani BG, Carnevali P, Nazarenko I, Nilsen GB, Yeung G et al: Human genome sequencing using unchained base reads on self-assembling DNA nanoarrays. Science 2010, 327(5961):78-81.
30. Liu W, Chen M, Bai L, Zhuang Z, Fan C, Jiang N, Zhao J, Ma S, Xiang X: Comprehensive transcriptomics and proteomics analyses of pollinated and parthenocarpic litchi (Litchi chinensis Sonn.) fruits during early development. Sci Rep 2017, 7(1):5401.
31. Li R, Yu C, Li Y, Lam TW, Yiu SM, Kristiansen K, Wang J: SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 2009, 25(15):1966-1967.
32. Chen M, Xu R, Rai A, Suwakulsiri W, Izumikawa K, Ishikawa H, Greening DW, Takahashi N, Simpson RJ: Distinct shed microvesicle and exosome microRNA signatures reveal diagnostic markers for colorectal cancer. PLoS One 2019, 14(1):e0210003.
33. Wu HJ, Ma YK, Chen T, Wang M, Wang XJ: PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res 2012, 40(Web Server issue):W22-28.
34. Wei S, Ma X, Pan L, Miao J, Fu J, Bai L, Zhang Z, Guan Y, Mo C, Huang H et al: Transcriptome Analysis of Taxillusi chinensis (DC.) Danser Seeds in Response to Water Loss. PLoS One 2017, 12(1):e0169177.

Tables

Table 1. Top 10 highly expressed miRNAs identified in the differentiated and dedifferentiated tissues of E. camaldulensis and E. grandis x urophylla.

Table 2. Target prediction of differentially expressed miRNAs in E. camaldulensis and E. grandis x urophylla.