Adventitious root formation during in vitro and ex vitro rooting
During in vitro rooting, compared with the control group, NAA did not induce ARs whereas IBA or IAA could. High concentrations of IBA and IAA inhibited rooting percentage (Fig. 1a), root number (Fig. 1b) and root length (Fig. 1c). Highest rooting percentage (72.50%) was obtained by 0.5 μM IAA 60 d after treatment. During ex vitro rooting, NAA, IBA and IAA significantly increased rooting percentage (Fig. 1d), root number (Fig. 1e) and root length (Fig. 1f) compared with the control group. A low concentration of IBA (1 mM) most effectively induced rooting, resulting in the highest rooting percentage (98.33%) and root length (2.72 cm) 25 d after treatment.
Ex vitro rooting induced the highest rooting percentage (98.33%) at 25 d, while highest rooting percentage during in vitro rooting was 72.50% at 60 d. Thus, ex vitro rooting induced AR from E. excelsum plantlets earlier (faster) than in vitro rooting. The samples collected from ex vitro rooting were used for the next analysis. Eight days after the 1 mM IBA ex vitro rooting treatment, AR primordia were evident, and ARs emerged from the epidermis after 10 d (Fig. 2). ARs elongated, rooting percentage was almost 100% by 25 d, and plantlet survival reached 100%.
IAA and H2O2 content analysis
In the 1 mM IBA treatment during ex vitro rooting, IAA content in stem bases increased gradually from 0 to 8 d, then dropped at 10 and 12 d. A sharp increase in IAA content was observed at 8 d (Fig. 3a). The trend of H2O2 content was different from that of IAA content (Fig. 3b). H2O2 accumulated rapidly after treatment, peaked at 2 d, then sharply decreased at 8 d. The highest content of IAA and lowest content of H2O2 at 8 d corresponded to the timing of AR primordia formation.
De novo assembly and sequence analysis
To identify genes involved in AR induction of E. excelsum plantlets during ex vitro rooting, 21 cDNA libraries were prepared from three repeat mRNA samples collected from 0 (ER0), 2 (ER2), 4 (ER4), 6 (ER6), 8 (ER8), 10 (ER) and 12 (ER12) d after 1 mM IBA treatment (Table 1). The total number of raw reads produced for each library ranged from 42,845,192 to 52,575,518 with Q20 > 97.48%. The sequenced raw data was submitted to the SRA at the NCBI database with the following accession numbers: SRR14278060, SRR14278059, SRR14278048, SRR14278046, SRR14278045, SRR14278044, SRR14278043, SRR14278042, SRR14278041, SRR14278040, SRR14278058, SRR14278057, SRR14278056, SRR14278055, SRR14278054, SRR14278053, SRR14278052, SRR14278051, SRR14278050, SRR14278049, and SRR14278047. After filtering, the clean reads per library ranged from 39,721,158 to 49,142,018 with the percentage of clean reads > 91.07% (Table 1). Trinity software was used to assemble clean reads and obtain transcripts and unigenes for subsequent analysis. The quality of transcripts and length distribution of unigenes are shown in Fig. S1.
The unigenes were processed in six databases to perform best hits by Blast with E values < 10-5, and inferred putative functions of the sequences were assigned. A total of 52,188 (40.15%) unigenes were matched to known genes in the NR database, 23,159 (17.82%) sequences to Pfam and 37,827 (29.10%) sequences in the Swiss-Prot database (Table 2). The NR database queries revealed that the annotated unigenes were assigned with a best score to sequences from the top seven species (Fig. 4): Vitis vinifera (21.76%), Theobroma cacao (4.34%), Coffea canephora (4.01%), Nelumbo nucifera (3.88%), Sesamum indicum (3.13%), Ziziphus jujuba (2.67%) and Manihot esculenta (2.25%).
The annotation of GO terms revealed that 24,939 unigenes (19.19%) were assigned to biological processes, molecular functions, and cellular components (Fig. S2). Most annotated unigenes in biological processes were involved in “cellular process”, “metabolic process”, and “single-organism process”. In the cellular component category, most annotated unigenes were annotated as “cell”, “cell part” and “membrane”. In the molecular functions, most annotated unigenes were categorized as “binding”, “catalytic activity” and “transporter activity”.
A total of 22,160 unigenes (17.05%) and 33 pathways were assigned based on metabolism, genetic information processing, environmental information processing, cellular processes and organismal systems pathway (Fig. S3). On the basis of KEGG analysis, most unigenes were annotated into “carbohydrate metabolism” of metabolism, “translation” of genetic information processing, “signal transduction” of environmental information processing, “transport and catabolism” of cellular processes, and “endocrine system” of organismal system.
The possible functions of unigenes were predicted and classified by alignment to the eggNOG database. A total of 50,632 unigenes (38.95%) were distributed into 25 categories (Fig. S4). Among them, the NOG category “general function prediction only” represented the largest group, followed by “function unknown”, “signal transduction mechanisms”, and “posttranslational modification, protein turnover, chaperones”.
DEGs in response to IBA-induced ex vitro rooting
Hierarchical clustering was used to analyze the expression patterns of DEGs in ER0, ER2, ER4, ER6, ER8, ER10 and ER12 libraries with three biological replicates. These DEGs were divided into nine main clusters (Fig. S5). DEGs in clusters 1, 2, 3 and 4 always showed high expression in the ER0 library with different trends in the other six libraries. The remaining five clusters represented DEGs with high expression levels induced by IBA treatment. The highest number of up-regulated genes was observed in the ER2 library (7,364) and fewest in the ER8 library (5,649) (Fig. 5a). Upset plot diagram analysis showed that 4,635 unigenes maintained differential expression after IBA-induced treatment from 2 to 12 d (Fig. 5b).
GO enrichment analysis
According to GO enrichment analysis, the degree of enrichment was measured based on the rich factor (higher rich factor represents greater enrichment), the FDR value (range from 0 to 1; a score close to 0 indicates more significant enrichment) and the number of genes enriched to a GO term. The significant enrichment GO terms of DEGs showed a few differences in the six libraries (Fig. S6). In the ER2 library, the significantly enriched terms were “monooxygenase activity”, “response to auxin” and “oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen”. In the ER4 library, the significantly enriched terms were the same as the ER2 library, but the term “response to auxin” was replaced by “photosynthesis, light reaction”. In the ER6 and ER8 libraries, “monooxygenase activity” and “photosystem” were listed as significantly enriched terms. In the ER10 library, “photosynthesis, light harvesting”, “chlorophyll binding” and “photosystem I” represented the main significantly enriched terms. In the ER12 library, the significantly enriched terms were “hydrogen peroxide metabolic process”, “hydrogen peroxide catabolic process” and “phenylpropanoid metabolic process”.
Besides “hydrogen peroxide metabolic process” and “hydrogen peroxide catabolic process”, several terms related to adversity stress were also identified in the GO enrichment analysis (Fig. 6). The terms “hydrogen peroxide metabolic process” and “hydrogen peroxide catabolic process” shared the same number and type of DEGs, most of which, including DEGs for the “response to oxidative stress” term, were up-regulated at 8 d after ex vitro treatment with IBA (Table S2). Most DEGs were associated with the term “response to abiotic stimulus”, followed by “response to oxidative stress” and “response to biotic stimulus”, while the fewest DEGs were associated with the term “response to hydrogen peroxide” (Fig. 6a). Most of the DEGs in the terms “response to abiotic stimulus” and “response to biotic stimulus” were up-regulated throughout the entire process of AR formation (Table S2). The terms “hydrogen peroxide metabolic process”, “hydrogen peroxide catabolic process” and “response to oxidative stress” encompassed 43 DEGs simultaneously (Fig. 6b), and these were mainly identified as genes related to cationic peroxidase and peroxidase (Table S2).
KEGG enrichment analysis
KEGG pathway enrichment analysis was performed in addition to GO enrichment analysis, and the pathways significantly enriched in each stage were similar (Fig. S7). DEGs were extremely enriched in “photosynthesis − antenna proteins”, “diterpenoid biosynthesis”, “brassinosteroid biosynthesis”, “flavone and flavanol biosynthesis” and “flavonoid biosynthesis” pathways in the ER2, ER4, ER6, ER8, ER10 and ER12 libraries. Most DEGs were enriched in “plant hormone signal transduction” and “phenylpropanoid biosynthesis” pathways.
DEGs enriched in plant hormone signal transduction pathway
KEGG pathway enrichment analysis showed that many DEGs were enriched in the “plant hormone signal transduction” pathway, and 132 up-regulated DEGs involved in auxin, cytokinine, gibberellin, abscisic acid, ethylene, brassinosteroids, jasmonic acid and salicylic acid signal transduction were identified, including genes encoding auxin-induced protein (AUX), auxin-responsive proteins (SMALL AUXIN UP RNA, SAUR; Auxin/Indole-3-Acetic Acid, AUX/IAA; Indole-3-acetic acid-amido synthetase Gretchen Hagen3, GH3), auxin transporter-like protein (AUX1/LAX), ethylene-responsive transcription factor (ERF), and others (Table S3). The up-regulated expression of 24 DEGs was maintained after the IBA-induced ex vitro treatment, mainly SAUR50, SAUR32, SAUR36, IAA11, IAA13, IAA16, GH3.1, GH3.5, GH3.6, and ERF.
A total of 42 up-regulated genes related to auxin-responsive proteins, including 13 IAA, 19 SAUR and 10 GH3 genes, were identified in the significantly enriched pathway “plant hormone signal transduction” (Fig. 7). Most of those genes were differentially expressed in early response to IBA treatment, at 2 and 4 d. Four IAA and three SAUR genes sustained up-regulated expression after IBA-induced treatment. IAA11 (TRINITY_DN5977_c2_g1) showed high expression with log2(FC) > 4 at ER2, ER4, ER6, ER8, and ER10 stages, and even with log2(FC) > 8 at the ER2 stage. IAA30 (TRINITY_DN1790_c2_g1) showed high expression with log2(FC) > 4 at all stages, and even with log2(FC) > 8 at the ER2 and ER4 stages. Some DEGs were specifically differentially expressed at certain stages, such as TRINITY_DN14445_c0_g2 at ER2, TRINITY_DN58019_c0_g2 at ER10, and TRINITY_DN1958_c1_g2, TRINITY_DN25899_c0_g1, TRINITY_DN58019_c0_g1 and TRINITY_DN9289_c0_g1 at ER12. Five GH3 genes were up-regulated at all stages, demonstrating log2(FC) > 4 for TRINITY_DN2748_c0_g1 and TRINITY_DN2800_c0_g1, and log2(FC) > 8 for TRINITY_DN5152_c0_g1.
The 12 up-regulated genes related to auxin-induced proteins, identified as AX6B, AX10A, A10A5, AX15A, AUX22, AUX22D, and AUX28, were sharply enriched in the “plant hormone signal transduction” pathway (Fig. 8a). Most of those genes were also differentially expressed at an early stage (2, 4 d) after IBA treatment. Among these genes, TRINITY_DN39834_c0_g1 and TRINITY_DN31248_c0_g1 were extremely highly differentially expressed with log2(FC) > 8 in the ER0 vs ER2 and ER0 vs ER4 comparisons. Four DEGs were up-regulated with log2(FC) > 1 by IBA treatment under ex vitro rooting in the ER0 vs ER2, ER0 vs ER4, ER0 vs ER6 and ER0 vs ER8 comparisons while TRINITY_DN5677_c0_g1 maintained up-regulated expression at all stages.
In addition, five up-regulated LAX genes (Fig. 8b) were significantly enriched in the “plant hormone signal transduction” pathway. Only TRINITY_DN9654_c0_g1 was up-regulated at all stages while the other four LAX genes were up-regulated at ER10 and ER12, except for TRINITY_DN380_c4_g1, which was up-regulated at ER6.
QRT-PCR analysis of gene expression
To further validate the results from the RNA-seq data, 10 candidate DEGs were selected for qRT-PCR analysis of E. excelsum samples that were collected 0, 2, 4, 6, 8, 10 and 12 d after 1 mM IBA treatment during ex vitro rooting. In the seven time points, the expression trend of the unigenes from qRT-PCR and RNA-seq analysis were largely consistent (Fig. 9). These results demonstrate that the transcriptome data accurately reflects the ex vitro response of IBA-induced AR formation of E. excelsum plantlets.