Induction of somatic embryogenesis
The procedure of somatic embryogenesis and regeneration in H. brasiliensis was established (Fig. 1) as described previously . The immature anthers were cultured in solid Murashige and Skoog (MS) medium supplemented with 2, 4-dichlorophenoxyacetic acid (2, 4 -D), kinetin (KT) and naphthylacetic acid (NAA) for 50 days. At the end of the period, the embryogenic calluses (ECs) were obtained. ECs were placed in the MS medium containing indole-3-acetic acid (IAA) and gibberellic acid (GA3) for embryo induction. After 40 days, primary embryos (PEs) were collected. The PEs were transferred to MS medium containing 6-benzyl aminopurineand AgNO3 for growing. After 40 days, two different embryos based on their phenotype (cotyledonary embryo (CE), abnormal embryo (AE)) were observed in the culture medium. . We observed a significant difference between CEs and AEs in phenotype. The CEs and AEs were placed on half-strength MS medium containing IAA and BA. The CEs turned stronger into the mature cotyledonary embryo (MCE) 20 days later, whereas the AEs turned brown and grown up into withered abnormal embryo (WAE). After 30 days, the MCEs grew into complete seedlings, whereas the WAEs turned black and died. Based on the above phenotypic observation, six different samples during SE were selected for further study.
Transcriptome analysis of rubber tree SE
High-throughput sequencing generated 915,535,874 raw reads in EC, PE, CE, AE, MCE and WAE samples. A total of 887,852,416 clean reads were retained by filtering the reads with adaptor sequences and ambiguous “N” base. The percentage of quality score above 30 (Q30) was 97.92% and the GC percentage was 43 % (Table 1). On average, 85.92% of the clean reads were mapped to H. brasiliensis genome.
All unigenes were annotated by the blast search against the public databases using BLASTx (E-value–5 ≤ 10). All 36937 unigenes were annotated in 4 databases involved in the Clusters of Orthologous Groups of proteins (COG) database, the Gene Ontology (GO) database, the clusters of euKaryotic Orthologous Groups (KOG) database and the Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups (eggNOG) database (Table 2). According to the COG functional classification, the 13421 unigenes were categorized into 25 COG categories. The four most highly represented COG categories were “general function prediction only” (20.57%), “transcription” (11.75%), “replication, recombination and repair” (11.53%) and “signal transduction mechanisms” (10.51%)（Fig. 2）. In addition, 19619, 20954 and 36362 unigenes were successfully annotated in GO, KOG, eggNOG, respectively (Fig. S1, S2, S3).
Global analysis of gene expression during rubber tree
A Venn diagram was created to find the overlapped genes in the four different developmental stages of H. brasiliensis SE (Fig.3a). A total of 25841 genes overlapped in the four stages. Among them, 155 genes overlapped between EC and PE; 290 genes overlapped between PE and CE; 193 genes overlapped between CE and MCE. A total of 388, 297, 152 and 582 genes were uniquely expressed in EC, PE, CE and MCE respectively. Another Venn diagram was also created to find the overlapped genes in the comparisons of PE, AE and CE of H. brasiliensis SE (Fig.3b). As shown in Fig.3b, 662 genes were exclusive to PE vs. AE. 1369 genes were exclusive to PE vs. CE. Moreover, 365 genes were found in AE vs. CE. To evaluate the differences of molecular response among all samples, the expression level of the unigenes was calculated by the expected number of Fragments Per Kilobase of transcript sequence per Million base pairs sequenced (FPKM). The top 20 expressed genes from EC, PE, CE and MCE libraries were shown in Table 3. Some of them including glutathione S-transferase（GST）, lipid-transfer protein（LTP）, peroxidase（POD）, indole-3-acetic acid-amido synthetase GH3.1, ADP-ribosylation factor, catalase isozyme, and polyubiquitin, were highly expressed in four stages.
In order to reveal the potential key factors and deeply understand the regulatory network of SE, the unigenes of each library of H. brasiliensis SE were compared under the condition of −1.0 ≥ Log2 [Fold Change (FC)] ≥ 1.0 and False Discovery Rate (FDR) < 0.01. A total of 9415 DEGs were obtained in EC vs. PE, PE had 5260 up-regulated and 4155 down-regulated genes. In PE vs. CE, CE had 1483 genes up-regulated and 2366 genes down-regulated. In CE vs. MCE, 6449 DEGs were obtained, of which 4016 DEGs were up-regulated, whereas 2433 DEGs were down-regulated. The 2820 DEGs were found in PE vs. AE with 1300 up-regulated and 1520 down-regulated DEGs. In AE vs. WAE, 5590 DEGs were obtained, of which 3318 DEGs were up-regulated, whereas 2272 DEGs were down-regulated. In AE vs. CE, 1536 DEGs were found with 556 up-regulated and 980 down-regulated DEGs. The 3307 DEGs were found between WAE vs. MCE with 1938 up-regulated and 1369 down-regulated DEGs (Fig. 4).
GO analysis of DEGs
To further demonstrate the unigenes functions, GO assignments were carried out using the Blast2GO program. In AE vs. CE, 843 DEGs were classified into three major categories: biological processes (BP), cellular components (CC) and molecular function (MF). A total of 41 GO subcategories were enriched over three major functional categories. The main subcategories are shown in Fig. 5a. The six major subcategories of BP were metabolic process, cellular process, single-organism process, biological regulation, localization and response to stimulus. The five major subcategories of CC were membrane, cell, cell part, organelle and membrane part. The four major subcategories of MF were binding, catalytic activity, transporter activity and nucleic acid binding transcription factor activity. In WAE vs. MCE, 1927 DEGs were also classified into BP, CC and MF and subcategorized in 41 GO (Fig 5b). The major subcategories of three categories were consistent with the result in AE vs. CE.
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway of DEGs
There were 376 DEGs in AE vs. CE, which were assigned to 46 KEGG pathways (Fig. 6a). The most representative pathways were phenylpropanoid biosynthesis (25 unigenes, Fig. S4A), plant hormone signal transduction (21 unigenes, Fig. S4B), starch and sucrose metabolism (20 unigenes, Fig. S4C), phenylalanine metabolism (19 unigenes), carbon metabolism (15 unigenes), biosynthesis of amino acid (14 unigenes) and glutathione metabolism (14 unigenes).
In WAE vs. MCE, the 771 DEGs were assigned to 57 KEGG pathways (Fig. 6b). The most represented pathways were phenylpropanoid biosynthesis (63 unigenes, Fig. S5A), starch and sucrose metabolism (49 unigenes, Fig. S5B), plant hormone signal transduction (46 unigenes, Fig. S5C), carbon metabolism (31 unigenes), photosynthesis (30 unigenes), phenylalanine metabolism (29 unigenes) and cyanoamino acid metabolism (29 unigenes). The results indicated that phenylpropanoid biosynthesis, phytohormones signaling pathway, and sucrose and starch metabolism played importance roles during H. brasiliensis late SE.
Differential expression of hormone signal transduction related genes between CE and AE
Various phytohormones induced SE and regeneration in several plants have already been reported. For instance, auxin was used alone or in combination with other plant growth regulators on plant SE induction [43, 44]. To further understand hormone regulation, FPKMs of hormonal signal transduction related genes were analyzed (Fig. 7a and Table S1). Among auxin signal transduction related genes, AUX-like5, IAA9-like, IAA28-like and GH3.1 were up-regulated in CE. SAUR71-like were highly expressed in AE than in CE. AUX22D-like, AUX28-like, AUX-like1, AUX-like2, SAUR32-like, IAA14-like and IAA27-like were highly expressed in MCE. ARF5-like was lowly expressed in CE but highly expressed in MCE. These genes participated in the auxin signaling pathway, which was important for cell enlargement and plant growth (Fig 7b).
Among abscisic acid (ABA) signal transduction related genes, PYL2-like was down-regulated in CE. PYL4-like was down-regulated in AE. Among jasmonic acid (JA) signal transduction related genes, JAZ7 was highly expressed in CE than in AE. JAZ5 was up-regulated in AE. Among ethylene (ET) signal transduction related genes, RAP2-3 was up-regulated in CE and in AE. RAP2-12-like and WRI1-like were highly expressed in CE. ERF4-like was up-regulated in MCE. ERF018-like was only up-regulated in AE. All the genes involved in the hormones signaling transduction pathways, including auxin, ABA, JA, ET, suggested that these hormones had an indispensable role in their complicated crosstalk process during H. brasiliensis late SE. In vitro studies have suggested the role of various regulatory genes in embryogenic transition that are triggered by plant hormones . The dynamic changes of these genes expression were critical for development of SEs.
Differential expression of TFs and SE-related genes between CE and AE
Transcription factors (TFs) play important roles in hormone signaling and stress responses as multifunctional regulators in both zygotic embryo and SE. Some of these TFs have been used as markers of totipotency in plant species . In the present study, we show that several TFs might play an important role during late SE of H. Brasiliensis. In this regard, 219 TFs were identified. The following TFs families were overrepresented: WRKY, MYB, MADS-box, AP2/ERF, bHLH. The expression profiles of 19 TFs in CE, AE, MCE and WAE are shown in Fig.8a and Table S2. WRKY40 and WRKY70 were up-regulated in CE and down-regulated in AE. WRKY23 were highly expressed in AE than in CE. MYB26-like and MYB98-like were up-regulated in AE. MYBS3-like and MYB1R1-like were up-regulated in MCE. AGL11 and AGL15 were up-regulated in AE. BBM2 was highly expressed in AE. AIL6 was highly expressed in CE than in AE. bHLH93-like was highly expressed in CE. The expression of bHLH94-like was up-regulated in AE. The results implied these TFs may play a key role in H. brasiliensis late SE.
Some SE-related genes, such as CAM , SERK [47, 48], LEA [49, 50], have been identified to play a vital role during plant embryogenesis. CML13 and CML36 were up-regulated in CE but down-regulated in AE. CAM-5-like and CAM (LOC110641724) were up-regulated in AE but had not changed in CE. CAM-7 was up-regulated in CE but down-regulated in AE. SERK1 was up-regulated in CE. LEAD-34-like and SERK2-like showed higher expression in AE than in CE. LEAD-29-like was up-regulated in MCE. The dynamic variation of the FPKM of these somatic embryogenesis-related genes suggested that they were critical for H. brasiliensis late SE.
Differential expression of histone modifications related genes between CE and AE
The plant growth regulators and abiotic stress can contribute SE. In the meantime, these factors may induce epigenetic modifications . Histone modification is one of the most important epigenetic modifications and plays a key role in the regulation of gene expression . Therefore, the expression levels of histone modifiers were analyzed and shown in Fig. 8b and Table S3. CURLY LEAF (CLF), encoding one of polycomb repressive complex 2 (PRC2) catalytic subunit that repress gene expression through trimethylating histone H3 at lysine 27 (H3K27me3), was higher expression in AE than in CE. The histone H3 lysine 9 methyltransferase genes (SUVH1-like, SUVH3-like, SUVH4-like and SUVH9), SUVR3-like, EZA1-like and ASHH3-like were expressed at a higher level in CE. In addition, histone demethylation related genes, LSD1-homolog 1-like was highly expressed in CE. LSD1-homolog 2 was up-regulated in MCE. The increased expression of genes in CE or MCE suggested that it is likely to have a function during late SE.
The acetylation of histones is believed to promote open chromatin state and activate gene transcription. Ten of the eleven genes related to histone acetylation showed significant differential expression in CE vs. AE. HAG6, HAC12-like, MCC1 and GCN5-like were up-regulated in CE. HAG11, HAG16, HAG18 and HATB-like were up-regulated in AE. 7 of the 13 genes related to histone deacetylation showed an obvious difference of expression in CE vs. AE. HDAC15-like and HDAC19 were highly expressed in CE. HDAC6-like, HDAC9 and SAP18-like were up-regulated in AE.
The histone phosphorylation related genes (Aurora-1, Aurora-2 like, Aurora-3 and Aurora-3 like) were highly expressed in AE than in CE. Plant Auroras can be divided into two categories according to the functions of Auroras. The alpha Auroras (Aurora 1 and Aurora 2) involve in controlling formative divisions throughout plant development. The beta Aurora (Aurora 3) associate with chromosome separation . These genes highly expressed in AE can be used as candidate genes for in-depth study in vitro embryogenesis.
qPCR verification of selected DEGs
To verify the reliability of transcriptome data, twenty genes related to SE were selected to carry out expression level analysis using qRT-PCR across 6 different tissues of H. brasiliensis (Fig. 9). Based on the transcriptome data analysis of H. brasiliensis SE, ARF4-like, GST, I2’H-like, PRX5-like, RBX1a-like, WRKY40 and WRKY70 were highly expressed in CE than in AE. E2 20-like, two EP3-likes, ERF9-like, FLC-like, five H3.2 genes, H3.2-like, MYB98-like and U17-like were lowly expressed in CE than in AE. The qPCR results validated the expression levels of 19 genes which were highly consistent with transcriptome data.