Transcriptomic profiling provides molecular insights into hydrogen peroxide-induced adventitious rooting in mung bean seedlings

Background Hydrogen peroxide (H2O2) has been known to function as a signalling molecule involved in the modulation of various physiological processes in plants. H2O2 has been shown to act as a promoter during adventitious root formation in hypocotyl cuttings. In this study, RNA-Seq was performed to reveal the molecular mechanisms underlying H2O2-induced adventitious rooting. Results RNA-Seq data revealed that H2O2 treatment greatly increased the numbers of clean reads and expressed genes and abundance of gene expression relative to the water treatment. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses indicated that a profound change in gene function occurred in the 6-h H2O2 treatment and that H2O2 mainly enhanced gene expression levels at the 6-h time point but reduced gene expression levels at the 24-h time point compared with the water treatment. In total, 4579 differentially expressed (2-fold change > 2) unigenes (DEGs), of which 78.3% were up-regulated and 21.7% were down-regulated; 3525 DEGs, of which 64.0% were up-regulated and 36.0% were down-regulated; and 7383 DEGs, of which 40.8% were up-regulated and 59.2% were down-regulated were selected in the 6-h, 24-h, and from 6- to 24-h treatments, respectively. The number of DEGs in the 6-h treatment was 29.9% higher than that in the 24-h treatment. The functions of the most highly regulated genes were associated with stress response, cell redox homeostasis and oxidative stress response, cell wall loosening and modification, metabolic processes, and transcription factors (TFs), as well as plant hormone signalling, including auxin, ethylene, cytokinin, gibberellin, and abscisic acid pathways. Notably, a large number of genes encoding for heat shock proteins (HSPs) and heat shock transcription factors (HSFs) were significantly up-regulated during H2O2 treatments. Furthermore, real-time quantitative PCR (qRT-PCR) results showed that, during H2O2 treatments, the expression levels of ARFs, IAAs, AUXs, NACs, RD22, AHKs, MYBs, PIN1, AUX15A, LBD29, LBD41, ADH1b, and QORL were significantly up-regulated at the 6- and/or 24-h time points. In contrast, PER1 and PER2 were significantly down-regulated by H2O2 treatment. These qRT-PCR results strongly correlated with the RNA-Seq data. Conclusions Using RNA-Seq and qRT-PCR techniques, we analysed the global changes in gene expression and functional profiling during H2O2-induced adventitious rooting in mung bean seedlings. These results strengthen the current understanding of H2O2-induced adventitious rooting and the molecular traits of H2O2 priming in plants. Electronic supplementary material The online version of this article (doi:10.1186/s12864-017-3576-y) contains supplementary material, which is available to authorized users.


Background
The plant root system is composed of primary roots, lateral roots, and adventitious roots. Adventitious roots refer to roots that form from any non-root tissue. These roots can generates from the tissues of old roots, stems or leaves during normal development to replenish and strengthen the roles of primary roots and in response to stresses [1]. The formation of adventitious roots is an important technique for commercial propagation of homogeneous horticultural plants and forest trees [2]. In addition, the process of adventitious root formation provides an ideal experimental system with which to study the important physiological and molecular events that occur during the tissue dedifferentiation and morphogenesis [3]. Many exogenous and endogenous factors that play important roles in adventitious rooting have been extensively studied using physiological and molecular methods. Auxin is known to play a critical role in inducing cell dedifferentiation and root primordia formation in cuttings. Other molecules, such as hydrogen peroxide (H 2 O 2 ) and nitric oxide (NO) that produced by cells in response to various stresses, also act as signalling molecules that are involved in modulating of the induction and initiation of adventitious root formation [4]. H 2 O 2 and NO are considered as key components of the molecular control of adventitious root formation in cuttings responding to wounding and auxin [5]. Li et al. reported that exogenous application of H 2 O 2 promoted adventitious rooting in mung bean and cucumber hypocotyl cuttings and that H 2 O 2 acts as a downstream signalling molecule in auxin-induced adventitious rooting [6,7]. Yang demonstrated that H 2 O 2 acts as a downstream modulator in SA-triggered adventitious root formation in mung bean seedlings [8]. However, the molecular mechanisms by which H 2 O 2 promotes adventitious root development remain elusive.
In our previous studies, we reported that exogenous H 2 O 2 strongly promoted adventitious rooting in mung bean hypocotyl cuttings [6]. Recently, we characterized the mung bean transcriptome and gene expression profiles during the induction and initiation stages of adventitious rooting in water-treated [29] or IBA-treated [30] mung bean hypocotyl cuttings using RNA-Seq technology. In the present study, RNA-Seq technology and qRT-PCR were exploited to reveal the transcriptome of H 2 O 2 promotion adventitious rooting in mung bean hypocotyl cuttings. Our main objective of the present study was to analyse the gene profiles and metabolic pathways that specifically responded to exogenous H 2 O 2 treatment and further to uncover the molecular basis of H 2 O 2 priming and promotion of adventitious rooting in plants.

Results and discussion
RNA sequencing, reads mapping, de novo assembly, unigene determination, and gene expression abundance analysis In this experiment, the hypocotyl tissues that were treated with 10 mM H 2 O 2 for 6 h or 24 h, designated as HO6 and HO24, which represent the induction stage and initiation stage of adventitious root development in mung bean seedlings, respectively [29,30], were separately sampled for total RNA extraction, cDNA library construction, and Illumina Solexa RNA paired-end sequencing. The sequencing data and mapping results are listed in Table 1. This sequencing generated approximately 93% of the high-quality reads from the samples, and approximately 93% of the quality reads mapped uniquely to the reference sequences that were constructed in our previous study [29,30]. Using the Chrysalis cluster module of TRINITY, we identified a total of 78,697 unigenes from all the samples [29,30]. The sequencing data were used for bioinformatics analysis with the data from the same batch samples treated with water for 6 h (Wat6) or 24 h (Wat24) that were reported in our previous study [29]. Data analysis showed that the number of clean reads and expressed genes were increased by 33.47% in HO6 and 2.98% in HO24 and by 9.21% (6142 genes) in HO6 and 7.12% (4602 genes) in HO24 compared with Wat6 and Wat24, respectively [29]; and the expressed gene number were increased by 4.71% and 2.0% compared with 6-and 24-h IBA treatments, respectively [30]. This result indicates that H 2 O 2 treatment greatly increased the number of expressed genes during two time points in comparison to the water and IBA treatments, particularly in the sample treated for 6 h. Cluster analysis of the samples showed that the HO24 and HO6 treatments separately grouped with Wat24 and Wat6, suggesting that the changes in gene expression abundance mainly depended on the time duration rather than the treatments during adventitious rooting in mung bean seedlings (Fig. 1a). Furthermore, gene expression abundance was increased in HO6 and HO24 relative to Wat6 and Wat24 (Fig. 1b).

Functional enrichment analysis of the genes expressed during H 2 O 2 -induced adventitious rooting
The functional enrichment results of the unigenes from the Clusters of Orthologous Groups for Eukaryotic Complete Genomes (KOG) database were clustered using Blast2GO software [31] and WEGO software [32] and designated as GO terms. The differentially expressed unigenes with 2fold changes (log2 > 1) between the samples under each GO category were further filtered. The statistical analysis results are listed in Additional file 1 and Table 2. Table 2 shows that the total number of GOs enriched in HO6 is 39.6% higher than that in HO24 and that there are 78.5% and 10.9% more up-and down-regulated GOs in HO6 than in HO24, respectively, indicating that a profound change in gene function occurred during the 6-h H 2 O 2 treatment (Table 2). There are 21.7% more up-regulated GOs than down-regulated GOs in HO6 and 32.3% more downregulated GOs than up-regulated GOs in HO24, suggesting that H 2 O 2 mainly enhanced gene expression levels at the 6-h time point but reduced gene levels at the 24-h time point.
We further analysed the ratio of differentially expressed genes (DEGs) in each up-and down-regulated GO category (Fig. 2). In the up-regulated GO group, nearly all of the DEG ratios in HO6 were greater than those in HO24, while the opposite result was observed in the down-regulated GO group, suggesting that most of the DEGs were significantly up-regulated by the 6-h H 2 O 2 treatment. The DEG ratios under the GO terms extracellular region and extracellular region part in cellular component, as well as electron carrier activity, antioxidant activity, enzyme regulator activity, and nutrient reservoir activity in molecular function were strongly increased from HO6 to HO24 in the up-regulated GO group, suggesting the DEGs in these GO terms were significantly up-regulated from 6 h to 24 h of treatment. In the downregulated GO group, DEG ratios under the GO terms cell killing in biological process and extracellular region part in cellular component were significantly increased in HO6, while DEG ratios under the GO terms antioxidant activity, structural molecule activity, and nutrient reservoir activity in molecular function were significantly increased in HO24, suggesting that the DEGs in these GO terms were down-regulated by H 2 O 2 treatment for 6 h and 24 h, respectively.
The most significantly up-and down-regulated GOs that specifically responded to H 2 O 2 treatment are listed in Additional file 1 and Table 3. In Summary, H 2 O 2 treatment for 6 h significantly up-regulated the GOs associated with the functions cellular component movement, transcription, DNA synthesis, cell cycle, hydrolase activity, and response to stress and down-regulated the GOs associated with oxidoreductase, cellular respiration and oxidative phosphorylation, and lipid transport. H 2 O 2 treatment for 24 h significantly up-regulated the GOs associated with oxidoreductase activity, transcription, isoflavonoid biosynthetic process, and ethylene biosynthetic process and down-regulated the GOs related to protein synthesis, cellular biosynthetic process, and cell wall organization. Interestingly, 6-and 24-h IBA treatment mainly up-regulated the GOs associated with auxin signalling, ribosome assembly and protein synthesis, and

KEGG pathway enrichment analysis
The KEGG pathways that were significantly regulated by H 2 O 2 were enriched using the KEGG Automatic Annotation Server (KAAS) [33]. The results showed that there were 58.1% more up-regulated KOs (204) than down-regulated KOs (129) in HO6, suggesting that H 2 O 2 treatment for 6 h mainly up-regulated certain metabolic pathways (Table 2). When compared with the water treatment, pathways such as glutathione metabolism, biosynthesis of amino acids, protein processing, DNA replication, purine metabolism, cutin, suberine and wax biosynthesis, and starch and sucrose metabolism were upregulated, while pathways including oxidative phosphorylation, phenylalanine metabolism, photosynthesis, flavonoid biosynthesis, ribosome, plant hormone signal transduction, and phenylpropanoid biosynthesis were down-regulated by H 2 O 2 treatment for 6 h or 24 h ( Table 4). The enriched DEGs under the glutathione metabolism pathway are mainly genes coding for glutathione S-transferase, which is involved in cell redox homeostasis and plays a vital role in the responses to various stresses [34]. Phenylpropanoids and flavonoids are major plant secondary metabolites that have a wide variety of functions as structural, signalling, and stress response molecules [35]. The derivatives from the phenylpropanoid pathway are also involved in the regulation of cell division and differentiation [36] and stimulate in vitro rooting [37]. Limonene and pinene are involved in terpenoid biosynthesis. These results indicate that H 2 O 2   (Fig. 1d, e, and f; Fig. 3a (Fig. 3c).

Gene expression patterns during H 2 O 2 -induced adventitious rooting
To gain insight into the DEGs patterns that specifically responded to H 2 O 2 treatment during two stages of adventitious rooting, the unigenes defined as unknown, hypothetical protein, uncharacterized protein, predicted protein, and unnamed protein in the database were filtered out; thus only the genes that had been assigned to proteins in the database remained. The remaining DEGs were further filtered, retaining those with RPKM ≥ 10 and fold change ≥ 2.

Top significantly regulated DEGs during root induction stage with H 2 O 2 treatment
Sixteen DEGs were up-regulated 8-to 55-fold in HO6 vs. Wat6, including 11 small heat shock protein (HSP)like genes, a BAG family molecular chaperone regulator 6-like, a germin-like protein, a multiprotein-bridging factor 1c-like, and an extensin class 1 protein (Additional file 2). These proteins are known to be involved in the response to various stresses, including H 2 O 2 [35]. Small heat shock protein and chaperone regulator, which is a modulator of chaperone activity, are involved in protein assembly and folding in response to stress [35]. Multiproteinbridging factor is a transcriptional coactivator that functions in the tolerance to heat and osmotic stress by partially activating the ethylene-response signal transduction pathway [35]. Extensin serves as structural constituent of the cell wall and is expressed in response to wounding [39]. A total of 10 DEGs were down-regulated 8-to 13-fold in HO6 vs. Wat6, including two peroxidase-like isoform genes, a cysteine-rich repeat secretory protein, a yieldin precursor, an auxin-responsive protein IAA12-like, a photosystem II protein D1, a 7-ethoxycoumarin O-deethylaselike (flavonoid 3′-monooxygenase), an E3 ubiquitin-protein ligase PUB23-like, and a polygalacturonase PG1 precursor gene (Additional file 3). The genes coding for cysteine-rich repeat secretory protein, peroxidase, and E3 ubiquitinprotein ligase are all involved in the stress response [35].
The isoform of peroxidase is also involved in the synthesis and degradation of lignin, suberization, and auxin catabolism [35]. Yieldin and polygalacturonase PG1 are glycoside hydrolases (GH) that are involved in cell wall loosening [35]. Photosystem II protein D1 functions in pathway of photosynthetic electron transport in photosystem II [35]. The auxin-responsive protein IAA12, a member of AUX/IAA family, acts as a repressor of auxin signalling and plays an important role during the different phases of adventitious root formation [5]. These results indicate that, during induction of adventitious   [35]. N-hydroxycinnamoyl/benzoyltransferase catalyzes the formation of phytoalexin, which are essential for the expression of disease resistance [35].
Metalloendoproteinase may play a role in the degradation and remodeling of the extracellular matrix during development or in response to stresses [35]. Zinc finger protein STOP1 acts as a transcription factor that was showed to be involved in aluminum tolerance [35]. Nitronate monooxygenase acts as a ferredoxin-dependent glutamate synthases, which have been implicated in a number of functions [35]. GST, isoflavone synthase, galactose oxidase, isoflavone reductase, pterocarpan reductase, mannitol dehydrogenase, formate dehydrogenase, polyphenol oxidase, and cationic peroxidase play vital roles in the responses to various stresses and cell redox homeostasis and were highly induced by exogenously applied H 2 O 2 [19] and IBA during IBA-induced adventitious rooting [30]. GSTs are widely studied versatile proteins that are involved in hydroxyperoxide detoxification [40], cellular redox homeostasis [40], and stress signalling proteins [41]. GSTs also function as non-catalytic auxin carriers, thus contributing to hormone homeostasis [40]. Isoflavone synthase belongs to the cytochrome P450 family enzymes, which are important for the biosynthesis of several compounds such as hormones, defensive compounds and fatty acids [35]. Pterocarpan reductase (also isoflavone reductase), isoflavone reductase homologue A622, and isoflavone reductase are involved in the pathway of pterocarpan phytoalexin biosynthesis and cell redox homeostasis [35]. The auxin-induced protein, an Aux/IAA protein family member, and F-box protein PP2-B15 are involved in auxin signaling pathways [35]. Glutamine amidotransferase catalyses purine biosynthesis from glutamine and is required for chloroplast biogenesis and cell division [35]. Limonoid UDP-glucosyltransferase, hydroquinone glucosyltransferase, and UDP-glycosyltrans ferase 73B3 belong to a large family of membrane-bound microsomal enzymes which catalyse the transfer of glucuronic acid to a wide variety of exogenous and endogenous lipophilic substrates. These enzymes are of major importance in the detoxification and subsequent elimination of xenobiotics such as drugs and carcinogens [35]. Anthranilate N-methyltransferase, isoliquiritigenin 2′-O-methyltransferase, embryo-abundant protein EMB, and theobromine synthase are methyltransferase family enzymes, mostly S-adenosylmethionine (SAM)-dependent methyltransferase, which utilise the ubiquitous methyl donor SAM as a cofactor to methylate small molecules, lipids, proteins, and nucleic acids, and therefore involved in many essential cellular processes including biosynthesis, signal transduction, protein repair, chromatin regulation and gene silencing [35]. ERF098 may be involved in the regulation of gene expression by stress factors and by components of stress signal transduction pathways [35]. Quinate hydroxycinnamoyltransferase is an acyltransferase involved in the biosynthesis of lignin [35]. The other proteins such as sugar transport protein 13, basic 7S globulin (an aspartic-type endopeptidase), lysine histidine transporter 1, syntaxin-24 (a vesicle trafficking protein), glutamate decarboxylase 1, and nonsymbiotic haemoglobin (an electron transfer) are involved in nutrients and energy metabolic processes. Interestingly, the gene coding for MATE efflux family protein FRD3-like was up-regulated 1789-fold after 24 h of treatment. This protein belongs to a secondary transporter family [42,43], which plays important roles in a wide range of biological processes in plants [44] and which has been shown to mediate the citrate efflux and iron transport to confer plant tolerance to Al toxicity [45]. Seven genes coding for small heat shock protein-like genes and a BAG family molecular chaperone regulator 6-like protein gene were highly up-regulated after 24 h of treatment. These proteins may act as protein chaperones that can protect other proteins against heat-induced denaturation and aggregation and are involved in plant defence processes [35]. Briefly, the above results indicate that, during initiation stage of adventitious rooting, H 2 O 2 treatment greatly up-regulated the expression of genes coding for proteins that function in various stress responses, thus improving cell tolerance to stress for further development.
The auxin transporter-like protein 1 acts as a carrier protein that is involved in proton-driven auxin influx and auxin gradient formation [35]. The auxin-responsive GH3 product belongs to IAA-amido synthetases that negatively regulate the levels of free IAA by conjugating IAA to amino acids [46]. The terminal flower 1a protein was shown to prevents the expression of ' APETALA1' and 'LEAFY' which are transcription factors that promotes early floral meristem identity is involved in cell differentiation [35]. The protein suppressor of PHYA-105 1-like was involved in suppression of photomorphogenesis in dark-grown seedlings and is required for normal elongation growth of adult plants [35]. The gibberellic acid-stimulated protein 1 has some role in plant development. The 21 kDa protein, which has pectinesterase inhibitor domain, and expansin-like B1 are implicated in the regulation of cell wall extension [35]. The above results indicate that, in comparison with 6-h treatment, a large number of genes were greatly up-regulated by 24-h H 2 O 2 treatment. These genes are mainly associated with stress response, redox homeostasis, cell wall loosening and remodeling, auxin signalling pathway, ethylene and gibberellin signalling pathways, protein and lipid degradation, and cell differentiation and elongation. Notably, the genes coding for endoglucanase 17-like, expansin-like B1like, MATE efflux family protein FRD3-like, and pectinesterase 2-like were highly up-regulated 2193-, 544-, 526-, and 182-fold, respectively, in HO24 vs. HO6, suggesting that cell wall loosening and remodeling and stress response were the profound molecular processes modulated by H 2 O 2 from 6-to 24-h. Six genes were down-regulated 8-to 134-fold in both HO24 vs. Wat24 and HO24 vs. HO6, including a cationic peroxidase 1-like, a germin-like protein subfamily 1 member 7-like, a CPRD49 (GDSL esterase/lipase), a 2-oxoglutarate/Fe(II)-dependent dioxygenaselike, a casparian strip membrane protein 3, a protease inhibitor, and a lipoxygenase gene (Additional file 7). The 2-oxoglutarate/Fe(II)-dependent dioxygenase belong to dioxygenase domain enzymes that is involved in the formation of plant hormones, such as ethylene, gibberellins, anthocyanidins and pigments such as flavones [35]. Lipoxygenases (LOXs) were involved in growth, development and oxidative stress [47] and were upregulated in tissue cultures for adventitious rooting [48]. LOXs oxidize linolenic acid present in membranes, resulting in the lipoperoxidation of membranes [49,50]. The down-regulation of this gene may alleviate the lipoperoxidation of membranes. Casparian strip membrane proteins (CASPLs) have recently been shown to form membrane scaffolds and direct the local modification of the cell wall in Arabidopsis thaliana [51].

Genes involved in plant hormone signalling were significantly regulated during H 2 O 2 -induced adventitious rooting
Plant hormones are known to play vital roles in adventitious rooting; however, we do not known whether H 2 O 2 promotes adventitious rooting by regulating the expression of plant hormone signalling-related genes. Thus, these genes were searched among all DEGs, and a total of 137 genes with both fold change > 2 and RPKM ≥ 5 were identified, including 21 ABC transport family members, 58 auxin pathway genes, 39 ethylene pathway genes, six cytokinin pathway genes, ten gibberellin-related genes, and three abscisic acid-related genes (Additional file 8 and Table 5).
Eight and three genes of ATP-binding cassette (ABC) transporter family members were up-regulated at 6 h and 24 h time points, respectively, but eight genes were down-regulated from 6 h to 24 h, suggesting that H 2 O 2 increased the expression of ABC transport family genes. ABC transporters have been known to participate in the export or import of a wide variety of molecules for growth and development under abiotic stress [52] and to act as an auxin carrier complex in cellular auxin efflux and influx [53].
In the group of auxin-related genes, a total of 58 genes were identified, including 26 auxin-induced proteins, ten auxin response proteins (AUX/IAA), nine indole-3-acetic acid (IAA)-amido synthetase GH3s, seven auxin response factor (ARF), four auxin transporter-like protein (LAX), an auxin efflux carrier (PIN1), and an auxinbinding protein ABP19a (ABP) gene (Additional file 7). Of the auxin-induced protein genes, 14 auxin-induced protein 5NG4-like genes were up-regulated at 6 h but down-regulated at 24 h and from 6 h to 24 h. Among these genes, two auxin-induced protein PCNT115-like genes were highly up-regulated. The auxin-induced protein 5NG4-like proteins belong to the plant drug/metabolite exporter (DME) family that consists of a group of plant WAT1 (walls are thin1)-related proteins, which has been shown to be involved in secondary cell wall deposition [54]. The AUX/IAAs were down-regulated at 6 h and from 6 h to 24 h. The ARFs were up-regulated at 6 h and 24 h but down-regulated from 6 h to 24 h. The IAA-amido synthetase GH3 genes were downregulated at 6 h but up-regulated at 24 h and from 6 h to 24 h. PIN1 and LAX1 were up-regulated 6-to 20-fold from 6 h to 24 h. ABP was up-regulated 136-fold from 6 h to 24 h. The auxin-inducible Aux/IAA family proteins function as repressors of early auxin response genes at low auxin concentrations [55]. GH3s catalyse the synthesis of IAA-amino acid conjugates, providing a mechanism for the plant to cope with the presence of excess auxin [35]. PIN1 and LAX1 are involved in  proton-driven auxin influx and mediate auxin gradient formation [56]. The ABP is probably an additional auxin receptor that mediates rapid cellular auxin effects through a non-transcriptional auxin response pathway [35]. These results suggest that H 2 O 2 increased the expression of auxin transporter and receptor genes and down-regulated ARFs and AUX/IAAs from 6 h to 24 h after the treatment, resulting in the initiation of auxin signalling during H 2 O 2 -promoted adventitious rooting. This modulation mechanism of auxin signalling is similar to the IBA-induced process [30,57] and was reported during adventitious rooting in petunia cuttings [5].
In the group of ethylene-related genes, 24 AP2-like ethylene-responsive transcription factor genes (AP2/ERFs), six 1-aminocyclopropane-1-carboxylate oxidase (ACO) genes, two 1-aminocyclopropane-1-carboxylate synthase (ACS) genes, and five ethylene-insensitive protein-like genes were identified. H 2 O 2 treatment up-regulated the expression of ACO and ACS genes at 24 h but downregulated their expression from 6 h to 24 h, indicating an increase in ethylene biosynthesis at 6 h compared with water treatment. Ethylene-insensitive protein 2 (EIN2), a central factor in signalling pathways regulated by ethylene, is involved in various processes including development, plant defences, senescence, nucleotide sugar flux, and tropisms [35]. The gene coding for ethylene-overproduction protein 1-like (ETO1), which functions as an essential regulator of the ethylene pathway by regulating the stability of ACS enzymes, was down-regulated from 6 h to 24 h [35]. Most AP2/ERF genes were down-regulated at 6 h but upregulated from 6 h to 24 h. These genes may be involved in the regulation of gene expression by stress factors and by components of stress signal transduction pathways [35] and have also been shown during adventitious rooting in petunia [58]. Ethylene and its crosstalk with auxin have been showed to be required for adventitious root formation [58,59]. These results suggest that the H 2 O 2 -induced process leading to adventitious rooting may modulate ethylene gene expression, and this mechanism has been revealed in IBA-induced adventitious rooting [30].
Six cytokinin-related DEGs were detected. Of these DEGs, a cytokinin-induced message gene that code for expansin-like B1 [35], was down-regulated 11-fold at 6 h but up-regulated 62-fold from 6 h to 24 h. A gene for cytokinin dehydrogenase 3 that degrades the phytohormone cytokinin [60], two genes for cytokinin hydroxylase-like that catalyses the biosynthesis of trans-zeatin in plants [35], and a cytokinin riboside 5′-monophosphate phosphoribohydrolase LOG1 that converts inactive cytokinin nucleotides to the biologically active free-base forms [35] were downregulated from 6 h to 24 h. This result indicates that H 2 O 2 treatment inhibited cytokinin activity via the regulation of cytokinin-related gene expression. The same result was shown in IBA-induced adventitious rooting [30].
Nine DEGs, including a gibberellic acid-stimulated protein 1, seven gibberellin oxidase genes, and a GIR1 gene, were down-regulated from 6 h to 24 h, while a gibberellin receptor GID1B-like was up-regulated during this period. Both gibberellic acid-stimulated protein 1 and GIR1 are the GASA gibberellins-regulated cysteinerich proteins that are involved in cell division and/or elongation and that respond to biotic and abiotic stresses [35,61,62]. The gibberellin oxidases are key enzymes in the biosynthesis of gibberellins and are involved in the production of bioactive GA for vegetative growth and development [35]. Similar to IBA treatment, H 2 O 2 decreased the expression of the GA synthesis genes, leading to inhibition of GA synthesis [30].
In summary, the above results indicate that, during the process of H 2 O 2 -induced adventitious rooting, H 2 O 2 initiated auxin and ethylene signalling and repressed cytokinin and gibberellins signalling for improving rooting.
Notably, Aux/IAA proteins are short-lived transcriptional factors that function as repressors of early auxin response genes at low auxin concentrations [35]. H 2 O 2 treatment down-regulated the AUX/IAA genes but upregulated the ARF genes, suggesting the involvement of H 2 O 2 in the auxin signalling, which has been shown to be a key regulator of adventitious rooting [1]. Furthermore, LBD family genes are up-regulated by auxin [30,63,64] and have been found to act early in auxin signalling and to regulate adventitious rooting [5]. Of those genes, LOB16 and LOB29 were upregulated 16-to 32-fold by IBA treatment [30]. In this study, five LBD family genes were up-regulated 4-to 31-fold by H 2 O 2 treatment from 6 h to 24 h, of which LOB16 was up-regulated 31-fold. Most of the AP2/ERF family genes were down-regulated at both 6 h and 24 h but up-regulated from 6 h to 24 h. For example, genes for ERF LEP and ERF086-like were up-regulated by 202-fold and 60-fold, respectively. These two genes have been shown to act as cell divisionpromoting factors involved in various processes of plant developments such as leaf blade differentiation and inflorescence branching [65]. ERF086 is also involved in the control of cell division patterns during early lateral root primordium development [65]. Although most of the AP2/ERF family genes were down-regulated at both 6 h and 24 h, the ERF At1g16060-like and ERF098-like were up-regulated by 7-and 238-fold at 6 h and 24 h, respectively. These two genes have been known to be involved in the regulation of gene expression by stress factors and by components of stress signal transduction pathways [35]. These results imply that H 2 O 2 treatment initiated the ethylene signalling pathway, leading to adventitious rooting via up-regulating both the expression of stress responsive ERF genes at 6 h and 24 h and cell developmental ERF genes from 6 h to 24 h. Recently, a number of studies have reported that the AP2/ERF family genes were the most highly regulated TFs [66] and were key endogenous regulators of adventitious rooting in petunia [58] and poplar [67] and IBAinduced adventitious rooting in mung bean [30]. In addition, the WRKY and NAC family genes were highly up-regulated at 6 h and 24 h. These two families of genes have been known to be involved in plant defence responses. NAC proteins are also involved in developmental processes, including the formation of the shoot apical meristem, floral organs and lateral shoots, as well as in plant hormonal control [35]. Some of NAC proteins, such as transcription factor JUNGBRUNNEN 1-  like, are involved in modulating cellular H 2 O 2 levels and enhancing tolerance to various abiotic stresses through the regulation of DREB2A [68].
Heat shock proteins (HSPs)-coding genes were significantly regulated by H 2 O 2 treatment Interestingly, we found that a large number of genes coding for heat shock proteins (HSPs) and heat stress transcription factors (HSFs) were significantly upregulated during H 2 O 2 treatment. Therefore, we screened the genes with both fold change > 2 and RPKM ≥ 5, and 29 HSPs and 5 HSFs were identified (Table 7). Amongst these genes, 25 HSP genes were up-regulated 2-to 98-fold at 6 h, 20 HSP genes were up-regulated 2-to 94-fold at 24 h, and 29 HSP genes were up-regulated 2-to 4-fold from 6 h to 24 h. Heat shock proteins are known to act as stress proteins, and their up-regulation is described more generally as part of the stress response [69,70]. Many researchers have reported that this family of proteins was produced by cells in response to exposure to stressful conditions, including heat shock, cold, UV light, nitrogen deficiency, water deprivation, and wounding [69,70], as well as exogenously applied H 2 O 2 [14]. Many HSP members function as intra-cellular chaperones that can protect other proteins against stress-induced denaturation and aggregation by stabilizing new proteins to ensure correct folding or by helping to refold proteins that were damaged by stress [69,70]. Clearly, our result indicates that the upregulation of HSP genes and HSF genes is an important mechanism under H 2 O 2 -induced adventitious rooting in mung bean seedling.

Validation of gene expression by qRT-PCR
The RPKM values of the DEGs obtained from the transcriptome were validated using real-time quantitative PCR (qRT-PCR). For qRT-PCR, we selected a total of 36 interesting DEGs; the results are shown in Fig. 5. Detailed information regarding these genes is presented in Additional file 10. We selected the housekeeping gene CPY20 as an internal reference gene for qRT-PCR measurement. The expression levels of most of the genes, as measured by qRT-PCR, showed a strong correlation to the RNA-Seq data (more than 83% of the dataset had a correlation coefficient r > 0.9). The qRT-PCR results showed that the expression levels of 26 genes, including ARF1 , ARF2, ARF3, ARF6, ARF8, ARF18,  ARF19, IAA8, IAA9, IAA14, IAA26, LAX4, AUX22E,  AUX22B, LBD24, NHX2, NAC, NAC21, NAC25,  The genes involved in the auxin signalling pathway are required for the formation of adventitious roots. The auxin efflux carrier PIN family members and auxin influx carrier AUX family members mediate auxin polar transport and control auxin distribution to establish and maintain auxin concentration gradients for various developmental processes [71] including adventitious rooting [57]. The ARFs [72], GH3-like proteins [73], and LBDs [48,74,75] have also been shown to be essential for adventitious root formation. NAC proteins are involved in plant developmental processes [76,77] and defence responses [78]. For example, NAC21 mediates auxin signalling to promote lateral root development [79]. AHK1 functions as an osmosensor histidine kinase and is required for the regulation of desiccation processes [80]. AHK2 acts as a cytokinin receptor and, together with AHK3, functions as a redundant negative regulator of drought and salt stress responses and abscisic acid signalling [81]. These proteins were shown to be involved in the water stress response during the early vegetative stages of plant growth and regulation of meristem development [82]. Arabidopsis mutants lacking AHK2, AHK3, and AHK4 exhibited enhanced adventitious root growth [83]. ADH1 (Alcohol dehydrogenase 1) and NHX2 (Na + /H + exchanger protein) are involved in the osmotic stress response [84]. Dehydration-responsive protein RD22 was up-regulated by dehydration, salt stress and abscisic acid and can enhance drought tolerance [85]. The Arabidopsis homologue of quinone oxidoreductaselike protein (QORL) belongs to chloroplastic NADPHdependent alkenal/one oxidoreductase (AOR) which contributes to detoxifying lipid peroxide-derived reactive carbonyls (RCs) produced under oxidative stress and to protect respiration and growth [86]. PERs are involved in the removal of H 2 O 2 , oxidation of toxic reductants, biosynthesis and degradation of lignin, suberization, auxin catabolism, and responses to environmental stresses, such as wounding, pathogen attack and oxidative stress. Briefly, these q-PCR results confirmed that H 2 O 2 -priming up-regulated the expression of several auxin signallingrelated genes, such as ARFs, AUX/IAAs, and PINs, and several genes from the LOB and MYB transcription factor families, particularly promoting the expression of various stress response-and oxidative-related genes, such as PERs and QORL, to further increase the stress tolerance of the cells.
In summary, all the above results indicate that a majority of most highly up-or down-regulated genes were similar during both H 2 O 2 -and IBA-induced adventitious rooting, suggesting the similar regulation pattern of gene expression in the processes [30]. This supports the concept that H 2 O 2 is an important component of auxin-mediated adventitious root formation in cuttings [5]. For example, both H 2 O 2 and IBA treatments up-regulated the genes coding for auxin signalling pathway proteins GH3s, ABP, and auxin-induced proteins, whereas down-regulated Aux/IAA and ARFs genes; up-regulated genes coding for ethylene signalling pathway proteins AP2/ERFs, ACS, and ACO, while downregulated some members of AP2/ERFs; up-regulated cytokinin-induced message genes but down-regulated cytokinin synthesis genes; up-regulated the cell wall related expansins and polygalacturonse genes while downregulated CASPL genes; up-regulated the most of LOBs and WRKYs transcription factor genes but down-regulated the most of bHLHs and MYBs genes; up-regulated GST, PODs, and polyphenol oxidase genes that are related to cell redox homeostasis but down-regulated some members of POD genes; up-regulated GDSL esterases and sugar transporter genes involved in metabolism but down-regulated LOX genes. However, some significant differences on gene expression levels were also observed between H 2 O 2 and IBA treatments. For example, major of the most highly regulated genes by H 2 O 2 were those involved in stress response, for IBA, were those involved in auxin signalling, and ethylene signaling pathways, as well as LOB family and meristems and primordial development-related. In addition, a lot of genes involved in secondary metabolism and flavonoid biosynthesis were down-regulated by IBA. Many of small HSPs and HSFs genes were up-regulated by H 2 O 2 . In short, IBA application primarily initiated the auxin signalling processes, whereas, H 2 O 2 -priming initially started the stress response processes, further leading to adventitious rooting in mung bean seedlings [30].

Conclusions
It is well known that H 2 O 2 functions as an important regulator molecule that is involved in many signalling Validation for the RNA-Seq data using qRT-PCR. The qRT-PCR data were compared with RNA-Seq data and the Con, Wat 6, and Wat 24 data were from our previous study [29]. Bars represent the mean (±SE) of three replicates. Different letters (a, b, and c) represent statistically significant differences (P <0.01) among the qRT-PCR data, which were analysed using Student's t-test By integrating the results of GO and KO enrichment and analyses of most highly differentially expressed genes in response to H 2 O 2 , we proposed a model including the most highly regulated pathways by H 2 O 2 -priming that might be involved in the induction and initiation of adventitious roots in mung bean seedlings (Fig. 6).

Plant material and culture conditions
After mung bean [Vigna radiata (L.) R. Wilczek] seeds were surface-sterilized in a 6% NaClO solution for 15 min and rinsed three times in sterile distilled water, they were covered with a thin layer of sterilized perlite and germinated in Petri dishes in a growth chamber at 25 ± 1°C for 36 h in the dark and then at 25 ± 1°C with a 14-h light/ 10-h dark photoperiod under white fluorescent lamps (PAR of 100 μM m −2 s −1 ). The five-day seedlings that were 5 cm in height were used as experimental material. After the primary roots were removed from the bases of the hypocotyls, the resulting seedlings were incubated in 50-mL beakers (10 per beaker), each containing 40 mL of sterilized distilled water or 10 mM H 2 O 2 under the same aseptic conditions as the seedling culture. To evaluate the gene expression profile during adventitious rooting, the basal part (5 mm in length) of each hypocotyl, where the adventitious roots would form, was harvested separately at the 0, 6, and 24 h incubation time points and used as samples for RNA extraction. Each sample contained ten stem cutting bases in three biological replicates in parallel. The samples were separately marked as Con, Wat6, HO6, Wat24, and HO24, immediately frozen in liquid nitrogen and stored at −80°C until further analysis.
Total RNA extraction, cDNA library construction and Illumina RNA-Seq For total RNA extraction, 50 mg of the samples was fully ground in liquid nitrogen, and total RNA extraction was performed using a plant total RNA Kit according to the manufacturer's protocol (kit SK8631, Sangon, Shanghai, China). The resulting total RNA was dissolved in 50 μL of RNase-free water, and stored at −80°C. The RNA integrity was quantified with RNA integrity number (RIN) values of 8.1-9.9 using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), and RNA concentration was determined using a NanoDrop ND-1000 Spectrophotometer.
Equal amounts of total RNA from each sample were used to purify the Poly (A) mRNAs using Oligo(dT) 25 beads (Invitrogen) according to the manufacturer's instructions. The Fragment Mix reactive system was used to fragment the purified mRNA at 94°C for 4 min. The firststrand cDNA was synthesized from fragmented mRNA templates using Superscript II reverse transcriptase (Invitrogen), random hexamer primers, and First Strand Master Mix at 25°C for 10 min, 42°C for 50 min, and 70°C for 15 min, with a final hold at 4°C. The RNA template was then removed and second-strand cDNA was synthesized using Second Strand Master Mix (Invitrogen). After the double-strand cDNA was synthesized, it was purified with Agencourt AMPure XP Beads (Agencourt), and the 3′ ends were repaired using End Repair Control, followed by purification with AMPure XP beads. Subsequently, an A-base was added to the blunt end fragments using Klenow exo (M0212L, NEB) to perform the adenylation of the 3′ ends of the cDNA fragments. Thereafter, Illumina indexing adapters were ligated to the cDNA fragments using T4 Ligase (Fermentas) according to the standard protocol. After the ligated cDNAs were purified twice using AMPure XP Beads, they were enriched and amplified using selective PCR according to the following procedure: 98°C for 30 s; 15 cycles of 98°C for 10 s, 60°C for 30 s, 72°C for 30 s, and 72°C for 5 min; holding at 4°C , followed by purification with AMPure XP beads. The quality and concentration of the cDNA library were confirmed using an Agilent 2100 Bioanalyzer and Qubit 2.0 (Life Technologies). The resulting cDNA libraries were then paired-end sequencing on an Illumina HiSeq 2000 system and performed at Sangon Biotech Co., Ltd. (Shanghai, China). The sequencing quality was checked using FastQC V0.10.1 with an ASCII Q-score offset of 33.

De novo assembly and sequence clustering
For de novo assembly, the adaptor sequences, low quality reads, and reads with unknown sequences > 5% were removed from the raw sequence, and then the resulting high quality sequence reads were de novo assembled using the Trinity paired-end assembly method [31]. The assembled sequences were clustered with Chrysalis. The longest sequences that could not be extended on either end within each clustered loci were obtained and defined as unigenes [29]. Finally, similarity alignment of the resulting unigenes were performed against the public protein and nucleotide sequence databases NR, SWISS-PROT, TrEMBL, Pfam, and CDD with similarity set at >30% and an E-value ≤ 1e-5 using BLASTx locally installed BLAST+ v2.2.27 software [87] and MEGA-BLAST, respectively. BLAST annotations were filtered using either subject or query coverage (>30%) and sequence identity (>50% for MEGABLAST and >30% for BLASTx). The results that presented the best alignment were used to identify the sequence direction and to predict the coding regions using BLASTx searches against protein databases, with the priority order of NR, SWISS-PROT, KEGG and KOG if conflicting results were obtained. The ESTScan software [88] was used to analyse the unigenes that did not align to any of the above databases. The assembled unigenes were deposited in the Transcriptome Shotgun Assembly Sequence Database (http://www.ncbi.nih.gov/genbank/tsa.html) at DDBJ/ EMBL/GenBank under the sequence read archive SRR 1653637 and the accession numbers GBXO01 000001-GBXO01078617.

GO and KEGG pathway enrichment of the unigenes
To understand which functional categories of DEGs were modulated during H 2 O 2 -induced adventitious rooting in mung bean seedlings, the unigenes were further aligned against the KOG (Clusters of Orthologous Groups for eukaryotic complete genomes) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway annotations using BLASTALL and KEGG Automatic Annotation Server (KAAS) software with an E-value ≤ 1e-5. The resulting blast hits were processed using Blast2GO software (version 2.3.5) [31] with an E-value threshold of 1e-5 to retrieve associated Gene Ontology (GO) terms, and WEGO software was used for achievement of GO classification [31]. To further enrich the pathway annotations, unigenes were submitted to the KAAS [33], and the single-directional best-hit information method was selected. To identify the enriched pathways, the phyper test was used to measure the relative coverage of the annotated KEGG orthologous groups of a pathway against the transcriptome background, and the pathways with a p-value ≤ 0.05 were classified as enriched. The difference analysis between a pair samples was performed between HO6 and Wat6, HO24 and Wat24, and HO24 and HO6 and was represented as HO6 vs. Wat6, HO24 vs. Wat24, and HO24 vs. HO6, respectively.

Unigene expression analysis and DEG confirmation
To calculate unigene expression levels, the unigene sequences were mapped to the reference sequences [29,30,89]. The clean reads were mapped back to the assembled unigenes using Bowtie 2 in the end-to-end alignment mode, and the number of clean reads that were mapped to each unigene was calculated and then normalized to RPKM (reads per kb per million reads) using ERANGE3.1 software [90]. The DEGseq R package was used to analyse the unigene expression levels [38] with the MARS (MA-plot-based method with Random Sampling) model. To further know the DEGs that were specifically regulated during H 2 O 2 -induced adventitious rooting, we screened the DEGs and their expression levels between HO6 and Wat6, HO24 and Wat24, and HO24 and HO6. The DEGs between each pair of samples were screened using the Audic-Claverie algorithm [91] with an FDR threshold of ≤0.001 and an absolute value of log2 ≥ 1, and the Benjamini-Hochberg correction was used to perform multiple test corrections for the p-value and FDR [92].

Quantitative reverse transcription PCR (qRT-PCR) validation
To validate the results from RNA-Seq, qRT-PCR was performed. For q-PCR, each of three biological replicates (n = 10) same as those subjected to Illumina RNA-Seq were used for total RNA extraction. Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. RNA integrity was examined with an Agilent Bioanalyzer 2100 (Agilent Technologies). In brief, cDNA samples preparation was performed as described in our previous study [30]. In this experiment, 36 genes that had been shown to be involved in auxin-related development and the stress response were selected for q-PCR analysis. The genes CPY20, eIF5A, and ACTIN (Actin-related protein 4) were selected as internal reference genes according to the RNA-Seq results and a previous report [93]. The genespecific primer pairs were designed using Primer Premier 5.0 software (Applied Biosystems, Foster City, CA, USA) according to the sequences from RNA-Seq (Additional file 10). Real-time PCR was run using a LightCycler 480 II (Roche Applied Science) and ABI SYBR Green PCR Master Mix (ABI, Foster, USA). The PCR cycling conditions were as follows: 95°C for 3 min, 40 cycles at 95°C for 15 s, and 60°C for 40 s. For negative controls, 'No cDNA' samples (water) and 'no RT' samples were included. All reactions were performed in triplicate. To evaluate the primer specificity for each primer set, melting curve analysis was conducted to verify the presence of a single melting peak. All q-PCR reactions were carried out with three independent biological replicates. The relative expression levels of the selected genes were calculated in relation to the reference gene using the comparative threshold cycle method with the delta-delta Ct method [94]. Statistical analyses were conducted with Student's t-tests at the P < 0.05 level of significance.