Morphological characterization of somatic callus
To study the somatic embryogenesis of Eucalyptus, the stem (differentiated) segments of both E. camaldulensis (A1) and E. grandis x urophylla (A2) were tissue-culture induced to callus (B1 for E. camaldulensis and B2 for E. grandis x urophylla). In these two Eucalyptus species, the morphologies (Figure 1A) of the tissue culture induced callus were similar and their growth curves (Figure 1B) showed S-shaped during dedifferentiation. The weight and appearance of the stem had no significant change during the first two days of induction. After 3 to 6 days induction, the two ends of the stem started to be induced to callus, which showed that both ends turned pale yellow, expanded and thickened like dumbbell (Figure 1A). Then, it fell into the rapid growth period of 7 to 12 days induction. The dedifferentiated callus gradually extended from both ends to the middle and the whole stem was completely dedifferentiated into pale yellow and moist callus. After 13 days of induction, it entered into the slow growth stage. The callus which had been cultured in the induction medium for 10 days were used for transcriptome sequencing.
Transcriptome sequencing and gene expression profiles
Next, we analyzed the transcriptome profiles of differentiated and dedifferentiated tissues of E. camaldulensis and E. grandis. Table 1 showed an overview of the transcriptome sequencing and the numbers of genes identified in each sample. Initially, we obtained 528 million reads for all 12 samples (in triplicates, n=3) with an average of 44 million reads. Then, the reads were aligned to the E. grandis genome and it showed that 64.66% to 74.84% of the total reads were mapped. Out of the total 36,349 E. grandis genes, we identified 18,777 to 20,240 genes with TPM > 1 in the differentiated and dedifferentiated tissues of Eucalyptus (Additional file 1). When we filtered genes using TPM < 5, only 12,650 to 14,464 genes were identified. Figure 1C showed that more than 60% of the genes were lowly expressed (TPM < 5) and only 0.165% to 0.272% of the Eucalyptus genes were expressed higher than 1,000 TPM. Interestingly, highly expressed genes varied in differentiated and dedifferentiated tissues of Eucalyptus. Figure 1D showed the comparison of top 10 highly expressed genes in all samples according to the average TPM. Dedifferentiated tissues of both E. camaldulensis and E. grandis x urophylla shared 8 highly expressed genes including Eucgr.H01085 (ethylene-responsive transcription factor ERF071), Eucgr.G03106 (wound-induced protein 1) and Eucgr.F00114 (zinc finger protein ZAT10). We next compared the genes identified in all the samples (TPM > 1) and found 10,374 genes shared (Figure 1E). In addition, 628, 459, 671 and 580 genes were specifically identified in A1, A2, B1 and B2, respectively. Based on the gene expression we analyzed the correlation between samples. It was revealed that replicates were performed well and that differentiated tissues were distinguishable from dedifferentiated tissues in E. camaldulensis and E. grandis x urophylla (Figure 1F).
DEGs in E. camaldulensis
To study the possible functional genes/pathways in the SE of Eucalyptus, we first identified DEGs in the dedifferentiated tissue of E. camaldulensis compared to the differentiated tissue. Using edgeR we identified 4,690 up-regulated and 4,539 down-regulated genes (Figure 2A, Additional file 2). It showed that Eucgr.H02264 (probable indole-3-acetic acid-amido synthetase GH3.1), Eucgr.D02625 (phosphoribulokinase, chloroplastic), Eucgr.J03055 (hypothetical protein), Eucgr.H02600 (protein SRG1) and Eucgr.B03016 (LOB domain-containing protein 40) were the top 5 up-regulated genes while Eucgr.J00794 (DNA-damage-repair/toleration protein DRT100), Eucgr.I02271 (endochitinase), Eucgr.A01080 (glycine-rich RNA-binding protein), Eucgr.D00703 (beta-galactosidase 8) and Eucgr.A01881 (trans-resveratrol di-O-methyltransferase) were the top 5 down-regulated genes, according to the FDR (false discovery rate) values (Additional file 2). We next analyzed the KEGG pathways enriched by the DEGs and found that “plant hormone signal transduction” (ko0475), “RNA transport” (ko03013) and “carbon metabolism” (ko01200) were the top 3 significant pathways. GO enrichment analysis (Figure 2B) showed that together with the GO terms enriched by more than 10% of the DEGs, “reproduction” and “reproduction process” were found with interest.
Then, we analyzed some gene groups might be involved in the dedifferentiation process of E. camaldulensis. Initially, 27 up-regulated and 22 down-regulated genes were found to associate with ethylene (Table 2), including 7 AP2-like ethylene-responsive TFs, 41 ethylene-responsive TFs and 1 ethylene response sensor (Additional file 2). In addition, 41 DEGs encoding auxin related products were found in the callus of E. camaldulensis compared to the stem (Table 2), including 7 auxin response factors, 3 auxin-binding proteins, 2 auxin-induced in root clusters proteins and 16 auxin-responsive proteins (Additional file 2). Furthermore, we found 111, 79, 36, 39 and 272 DEGs (Table 2) encoding RP, ZFP, HSP, histone and TF, respectively, might be involved in the dedifferentiation process of E. camaldulensis. Notably, 8 up-regulated and 3 down-regulated embryogenesis related genes were identified in the callus of E. camaldulensis compared to the stem, such as late embryogenesis abundant protein, somatic embryogenesis receptor kinas 1 and embryogenesis-associated protein EMB8 (Additional file 2). Further, based on the GO annotation we identified 274 DEGs (154 up-regulated and 120 down-regulated) related to cell wall (Table 2). Details of these cell wall related DEGs in different categories, such as “GO:0009505~plant-type cell wall” and “GO:0009834~plant-type secondary cell wall biogenesis”, can be accessed in the additional file 3. Among them, 9 DEGs encoding expansin were up-regulated and 1 gene encoding PME53 (probable pectinesterase 53) were down-regulated.
DEGs in E. grandis x urophylla
We next identified 4,200 up-regulated and 4,708 down-regulated genes in the dedifferentiated tissue of E. grandis x urophylla compared to the differentiated tissue (Figure 2C). According to the FDR, top 5 DEGs include Eucgr.B03016 (LOB domain-containing protein 40), Eucgr.H02264 (probable indole-3-acetic acid-amido synthetase GH3.1), Eucgr.I02271 (endochitinase), Eucgr.L02894 and Eucgr.L02534 (Additional file 2). Notably, pathway analysis identified “Plant hormone signal transduction” (ko04075) significantly enriched by 260 DEGs (q-value: 8.36E-14), which indicates that these DEGs may play key roles in the SE of E. grandis x urophylla. GO analysis (Figure 2D) showed that most of the GO terms involved by the DEGs were shared by E. grandis x urophylla and E. camaldulensis but some might be specific to E. grandis x urophylla, such as biological adhesion, virion and nucleoid.
Table 2 showed the numbers of DEGs from the seven gene groups identified in the dedifferentiation process of E. grandis x urophylla. Interestingly, we found more DEGs (41 up-regulated and 8 down-regulated) encoding ribosomal proteins in E. grandis x urophylla than E. camaldulensis, especially some RP genes which were up-regulated in E. grandis x urophylla but down-regulated in E. camaldulensis. In addition, compared to E. camaldulensis more down-regulated ethylene related genes and more up-regulated HSP and histone genes were found in E. grandis x urophylla (Table 2). Figure 2E showed an overview of all the DEGs identified in the dedifferentiated tissues compared to the differentiated tissues in both E. camaldulensis and E. grandis x urophylla. It showed that some DEGs were specifically identified in E. camaldulensis or E. grandis x urophylla, which might relate to the regenerative ability of Eucalyptus.
Regenerative ability associated genes
We next compared the DEGs identified in the dedifferentiation process of E. camaldulensis and E. grandis x urophylla. It showed in the upper panel of Figure 3A that they shared 2,687 up-regulated genes in the dedifferentiated tissue compared to differentiated tissue, including Eucgr.H01085 (ethylene-responsive transcription factor ERF071), Eucgr.A01538 (fructose-bisphosphate aldolase 6, cytosolic), Eucgr.G03106 (wound-induced protein 1), Eucgr.K02614 (NDR1/HIN1-Like protein 3), Eucgr.F00114 (zinc finger protein ZAT10) and Eucgr.H03082 (early nodulin-75) (Additional file 2). There were 2,003 and 1,513 up-regulated genes specifically identified in the callus of E. camaldulensis and E. grandis x urophylla, respectively (upper panel of Figure 3A). Then, we compared to down-regulated genes and found that they shared 2,581 genes (lower panel of Figure 3A), including Eucgr.J00025 (heat shock cognate 70 kDa protein 2), Eucgr.B01596 (probable xyloglucan endotransglucosylase/hydrolase protein 23), Eucgr.A01080 (glycine-rich RNA-binding protein), Eucgr.F00590 (snakin-2) and Eucgr.H03983 (major allergen Pru ar 1) (Additional file 2). A total of 1,958 and 2,127 down-regulated genes were specifically identified in E. camaldulensis and E. grandis x urophylla, respectively (lower panel of Figure 3A).
Interestingly, we identified 6, 17, 12, 83, 28, 10, 9, 122 and 98 DEGs related to the embryogenesis, ethylene, auxin, RP, ZFP, HSP, histone, cell wall and TF, respectively, only in the dedifferentiation process of E. camaldulensis (Table 2, Additional file 2). The 6 embryogenesis related genes include 5 up-regulated genes (2 SE receptor kinase, 3 LEA) and 1 down-regulated gene encoding embryogenesis-associated protein EMB8. Among the 17 DEGs encoding auxin related products (Figure 3B) specifically identified in E. camaldulensis, Eucgr.K02992 (auxin transporter-like protein 4) and Eucgr.C02984 (auxin-responsive protein IAA26) were down-regulated and Eucgr.H03171 (auxin-induced protein 22D) was up-regulated in E. camaldulensis but down-regulated in E. grandis x urophylla (Additional file 2). All the DEGs related to ethylene and specifically identified in the callus of E. camaldulensis were found to encode ER TFs, and 7 had reverse regulation in E. camaldulensis and E. grandis x urophylla (Figure 3C). Interestingly, heat maps (Figure 3D and Figure 3E) showed that most of the E. camaldulensis specific DEGs encoding HSP and RP were down-regulated in the dedifferentiated tissue compared to the differentiated tissue. Although some of the DEGs encoding histone were more or less changed in E. grandis x urophylla, the difference of their expression was not significant as that in E. camaldulensis (Figure 3F). Importantly, we found another 4 TF subfamilies including ASIL2, bHLH, MYB and WRKY were dysregulated only in E. camaldulensis (Figure 3G). Further, the expression of MYB and WRKY TF genes were elevated during the dedifferentiation process of E. camaldulensis. The cell wall related DEGs specific to E. camaldulensis were involved in multiple process, such as “GO:0005199~structural constituent of cell wall”, “GO:0009664~plant-type cell wall organization” and “GO:0009832~plant-type cell wall biogenesis” (Additional file 3). It is interesting that Eucgr.F01000 (formin-like protein 5, log2FC=1.65, p=5.67E-05) were the only DEG involved in the “GO:0005199~structural constituent of cell wall”. In addition, some other gene families were identified to be specifically differentially expressed in the SE of E. camaldulensis (Table 2, Additional file 3), such as 3 ABA related genes, 3 arabinogalactan protein genes, 4 ABC transporter genes and 21 abscisic stress-ripening protein genes.
qRT-PCR validation
We used qRT-PCR to confirm the expression patterns of DEGs in the dedifferentiation process of E. camaldulensis and E. grandis x urophylla. We randomly selected 9 genes (Eucgr.A00971, Eucgr.A01091, Eucgr.B03715, Eucgr.C03048, Eucgr.D01811, Eucgr.F00490, Eucgr.F01164, Eucgr.H03077 and Eucgr.K01605) for the qRT-PCR experiment and the H2B gene was used as the internal control. The primer sequences for all these genes can be accessed in the Additional file 4. Each gene was replicated three time in every sample, so we performed 9 reactions in total for the differentiated and dedifferentiated tissues. Log2 fold change (log2FC) and log2 relative normalized expression (log2RNE) were used to present the gene changes detected by transcriptome sequencing and qRT-PCR, respectively. Overall, 14 (77.8%) out of 18 events were agreed by both qRT-PCR and deep sequencing (Figure 4). Interestingly, the up-regulation of SE marker gene SERK1 was detected by qRT-PCR in both Eucalyptus species while transcriptome sequencing only detected its up-regulation in E. camaldulensis. However, other DEGs encoding SERK1 were identified with up-regulation in E. grandis x urophylla (Additional file 2). It is notable that the increase of WRKY TF (Eucgr.D01811) and the decrease of RP gene (Eucgr.A00971) in E. camaldulensis were confirmed by qRT-PCR. The dysregulation of Eucgr.B03715, Eucgr.C03048 and Eucgr.F00490 in E. grandis x urophylla were also confirmed. High agreement of gene expression patterns in transcriptome sequencing and qRT-PCR indicate that the genes identified in this study might be associated with the regenerative ability and the SE of Eucalyptus, which requires future experiments to be explored.
SNPs and indels
We next identified gene variants (e.g., SNPs and small indels) in E. camaldulensis and E. grandis x urophylla using the transcriptome sequencing data. Initially, we obtained 97,504 and 75,582 variants in the differentiated and dedifferentiated tissues of E. camaldulensis, respectively (Table 3). After the variants supported by <100 reads were filtered, we identified 97,974 variants for E. camaldulensis. Likewise, 72,208 and 66,311 variants were found in the differentiated and dedifferentiated tissues of E. grandis x urophylla, respectively, and they produced 78,977 variants after filtering variants with low supportive reads (Table 3). Comparison showed 49,527 variants shared by E. camaldulensis and E. grandis x urophylla, and 48,447 variants were specifically identified in E. camaldulensis. Then, we annotated the E. camaldulensis specific variants using ANNOVAR and found that 13,434 variants were functional, such as nonsynonymous, frameshift insertion, frameshift deletion, stop gain and stop loss variants (Table 3, Additional file 5). Interestingly, these 13,434 variants were derived from 4,723 Eucalyptus genes, such as annexin (Eucgr.F02423, Eucgr.H00564), ARF guanine-nucleotide exchange factor GNOM (Eucgr.B03196), AP2-like ER TF (Eucgr.I00278, Eucgr.J02113), auxin response factors (Eucgr.G00076, Eucgr.C02178, Eucgr.C03293, Eucgr.J00923, Eucgr.F02090, Eucgr.D00264, Eucgr.E00888) and wall-associated receptor kinas-like (Eucgr.I01022). KEGG pathway analysis showed that 85 pathways including “longevity regulating pathway” (ko04211) and “Plant hormone signal transduction” (ko04075) were enriched by the of the E. camaldulensis specific mutated genes (Additional file 6).