De novo sequencing and comparative transcriptome analysis of adventitious root development induced by exogenous indole-3-butyric acid in cuttings of tetraploid black locust

Background Indole-3-butyric acid (IBA) is applied to the cuttings of various plant species to induce formation of adventitious roots (ARs) in commercial settings. Tetraploid black locust is an attractive ornamental tree that is drought resistant, sand tolerant, can prevent sand erosion and has various commercial uses. To further elucidate the mechanisms of AR formation, we used Illumina sequencing to analyze transcriptome dynamics and differential gene expression at four developmental stages in control (CK) and IBA-treated groups. Results The short reads were assembled into 127,038 unitranscripts and 101,209 unigenes, with average lengths of 986 and 852 bp. In total, 10,181 and 14,924 differentially expressed genes (DEGs) were detected in the CK and IBA-treated groups, respectively. Comparison of the four consecutive developmental stages showed that 282 and 260 DEGs were shared between IBA-treated and CK, suggesting that IBA treatment increased the number of DEGs. We observed 1,721 up-regulated and 849 down-regulated genes in CI vs. II, 849 up-regulated and 836 down-regulated genes in CC vs. IC, 881 up-regulated and 631 down-regulated genes in CRP vs. IRP, and 5,626 up-regulated and 4,932 down-regulated genes in CAR vs. IAR, of which 25 up-regulated DEGs were common to four pairs, and these DEGs were significantly up-regulated at AR. These results suggest that substantial changes in gene expression are associated with adventitious rooting. GO functional category analysis indicated that IBA significantly up- or down-regulated processes associated with regulation of transcription, transcription of DNA dependent, integral to membrane and ATP binding during the development process. KEGG pathway enrichment indicated that glycolysis/gluconeogenesis, cysteine and methionine metabolism, photosynthesis, nucleotide sugar metabolism, and lysosome were the pathways most highly regulated by IBA. We identified a number of differentially regulated unigenes, including 12 methionine-related genes and 12 ethylene-related genes, associated with the KEGG pathway cysteine and methionine metabolism. The GO enrichment, pathway mapping, and gene expression profile analyses revealed molecular traits for root induction and initiation. Conclusion Our study presents a global view of the transcriptomic profiles of tetraploid black locust cuttings in response to IBA treatment and provides new insights into the fundamental mechanisms associated with auxin-induced adventitious rooting. Electronic supplementary material The online version of this article (doi:10.1186/s12864-017-3554-4) contains supplementary material, which is available to authorized users.


Background
Tetraploid black locust (Robinia pseudoacacia L.) is an attractive ornamental tree that has various commercial uses. It exhibits fast growth, drought resistance, salinealkaline tolerance and low soil fertility requirements, and it is the primary tree species used as a wind break and for sand fixation and for soil and water conservation in the Loess Plateau Region [1]. Tetraploid black locust does not root well and thus is difficult to plant in a variety of environments. However, the phytohormone auxin can promote the formation of adventitious roots (ARs) in cuttings [1,2]. The formation of ARs relies on the method used. In agricultural practice, plant loss is usually caused by AR formation and slow-rooting cuttings. Therefore, it is believed that the formation of ARs is necessary for the smooth spreading of cuttings of tetraploid black locust.
AR formation is a highly complex regenerative process that is influenced by numerous internal and external factors, including environmental conditions, phytohormones and nutritional status [3][4][5][6]. Auxin is a crucial phytohormone that promotes AR formation in cuttings [7]. However, the mechanisms underlying the role of auxin in AR formation are only superficially understood, and the lack of details at the molecular level limits improvements to cutting propagation.
The auxin IBA is widely used in woody plant propagation to induce rooting. Although indole-3-acetic acid (IAA) is the primary native auxin in plants, IBA is more effective in promoting ARs [8][9][10][11]. For example, treatment of cuttings with IBA significantly improves rooting rates in Terminalia tomentosa [8], Pinus contorta [12,13], Malus pumila [14,15], and Pinus radiate [16]. Brinker et al. [17] reported that IBA induces the expression of genes involved in cell replication and cell-wall weakening but inhibits genes related to auxin transport, photosynthesis and cell-wall synthesis during P. contorta root initiation. Thus, the processes that occur in cuttings after IBA treatment, and particularly the functions of IBA-regulated genes, should be elucidated. These results also indicate that IBA may directly or indirectly induce the formation and differentiation of root primordia. We previously demonstrated that 5.4 mmol/L IBA significantly increases the rooting rate of tetraploid black locust hardwood cuttings to approximately 80% [1,2], in contrast to approximately 2% AR formation in CK cuttings. To explore the significant impact of IBA on rooting, our recent studies have primarily focused on tetraploid black locust at the anatomical, physiological and biochemical levels [1,2]. Developments in molecular biology and proteomics techniques have allowed us to investigate IBA-induced AR development in tetraploid black locust via homology cloning, quantitative real-time PCR (qPCR), two-dimensional electrophoresis and protein bio-mass spectrometry (MALDI-TOF, Q-TOF), resulting in the isolation and identification of hundreds of IBA-response-related genes and proteins [18][19][20]. The studies cited above represent the first explorations of genes involved in AR formation in tetraploid black locust, but transcriptomic information and identification of many genes related to IBAinduced AR development is scarce. The molecular mechanism of rooting in tetraploid black locust is complex, and the mechanism by which IBA promotes the formation of ARs in cuttings remains unclear due to the lack of transcriptomic and genomic information. Therefore, our study represents a necessary acceleration of the acquisition of transcriptomes related to IBA-induced AR development in tetraploid black locust cuttings. Transcriptomic studies of IBAinduced AR development have been conducted in Camellia sinensis [21], Petunia hybrid [22], Pinus contorta [23], and Morus alba L. [24]. These studies primarily focused on stem cuttings and involved techniques such as Illumina sequencing [21], massively parallel signature sequencing [17], EST analysis [23], microarray analysis, and suppression subtraction hybridization [11]. However, the transcriptome of tetraploid black locust has not been studied yet.
The rapid development of next-generation sequencing (NGS) technology has improved the efficiency and reduced the cost of Illumina/Solexa sequencing technology. The results of Illumina/Solexa sequencing are highly reproducible, both technically and biologically [23,25]. Illumina/Solexa sequencing technology is the most widely used NGS technology for the de novo sequencing and analysis of the transcriptomes of nonmodel organisms.
In this study, we used the Illumina sequencing platform to analyze CK and IBA-treated cuttings of tetraploid black locust in order to identify new genes involved in the IBA-induced formation of ARs in cuttings and to obtain deeper insight into the mechanism of propagation in tetraploid black locust. The application of Illumina next-generation sequencing provides more transcripts to facilitate further genomic studies of tetraploid black locust. This study presents the transcriptome for IBA-treated cuttings and provides a genetic resource for improving woody plant propagation.

IBA induced adventitious root formation in tetraploid locust
We found that approximately 80% of the softwood cuttings form roots after treatment with the optimal concentration of IBA. The process of AR formation in the softwood cuttings involves the following steps ( Fig. 1): First, softwood cuttings after IBA treatment via soaking were inserted into seeding beds (Fig. 1a). During the 7-10 days after IBA treatment, we observed that white calli had formed on the wound surfaces of the soft cuttings (Fig. 1b). During the 15-20 days after the cuttings were treated, tiny AR primordia (RP) formed and subsequently developed into root meristems ( Fig. 1c and d). During the last stage, the ARs on the cuttings were formed and elongated (Fig. 1e).

Illumina sequencing, de novo assembly and sequence analysis
To identify genes involved in IBA-induced AR formation, we subjected cDNA preparations from the basal parts of IBA-treated or CK cuttings to de novo sequencing on the Illumina HiSeq 2000 platform. In total, eight cDNA preparations were sequenced from control cuttings sampled at the first stage (CI), white callus stage (CC), primordia formation stage (CRP), and AR formation stage (CAR), as well as from IBA-treated cuttings sampled at the first stage (II), white callus stage (IC), primordia formation stage (IRP), and AR formation stage (IAR). The total number of raw reads produced for each library ranged from 33.88 million to 46.99 million. The raw data have been submitted to the NCBI repository (https://www.ncbi.nlm.nih.gov/geo/). After filtering, the number of high-quality clean reads per library ranged from 33.74 million to 46.81 million, and the Valid Ratio (Reads) % ranged from 98.96 to 100.00% (Table 1). The short reads were assembled into 127,038 unitranscripts and 101,209 unigenes with average lengths of 986 and 852 bp, total lengths of 125,353,356 and 86,239,985 bp, and N50 lengths of 1,643 and 1,449 bp, respectively (Additional file 1). The number of reads per kilobase of exonic sequence per million of total reads sequenced (RPKM) was used to calculate the transcript abundance in each sample. The average RPKM ranged from 7.41 to 10.44 in CK group and from 7.72 to 9.49 in the IBA treatment group (Table 1). These results indicate that overall transcript abundance was greatly increased in both the CK and IBA-treated groups over the course of AR development. Moreover, the average RPKM of the IBA treatment group was greater than that of the CK at the first stage, white callus stage and primordial formation stage, reflecting a marked increase in gene transcription produced by IBA treatment.
The unigenes were aligned to five public protein databases (Nr, KOGs, KEGG, Pfam, and Swiss-Prot). The majority of the unigenes were annotated using the Nr (74.36%) database (Additional file 2). These results indicate that the Nr database is an informative platform for the functional annotation of tetraploid black locust. Additionally, the Nr database queries revealed that the highest percentage of tetraploid black locust sequences most closely matched sequences from Glycine max (41.1%), followed by Cicer arietinum (22.5%), Medicago truncatula (10.8%), Lotus japonicus (4.2%), Vitis vinifera (2.9%) and Theobroma cacao (1.5%) (Fig. 2). The results indicated that overall homology was highest to model legume plants, implying that the sequences of the tetraploid black locust transcripts obtained were assembled and annotated properly in this study [26].

Transcriptome changes during AR development in tetraploid black locust cuttings
In this study, we used IBA-treated cuttings as a model in which to investigate transcriptome changes during AR c b e a d Fig. 1 Morphological changes in tetraploid black locust cuttings undergoing adventitious root development in a sand bed. a Softwood cuttings before cutting. b White callus appeared 10 days after cutting. c-d Yellow callus appeared and tiny adventitious roots emerged (root primordium) at 15 days after cutting. e Adventitious root formation and elongation at 20 days after cutting. As biological replicates, 10 samples were randomly selected from the groups treated with IBA  35 were not expressed at stage II or stage IC. The cluster analysis also revealed that the abundances of 89.7% of the transcripts detected in the IBA-treated cuttings varied over the course of AR development (Fig. 3). Comparison of the expression patterns of IBA-treated ( Fig. 3) and CK (Additional file 3) cuttings revealed that the 40 groups were common to both the IBA-treated and CK cuttings, and the expression patterns of 97.5% of the genes expressed in CK were similar to those of the genes expressed in the IBAtreated cuttings (the exceptions belonged to groups 17 and 23). Moreover, the expression levels of 282 DEGs (group 23) in the IBA-treated cuttings and of 260 DEGs (group 5) in CK increased continuously. The expression of 61 DEGs increased continuously in both IBA-treated and CK cuttings (Additional file 4). These results suggest that a substantial alteration of gene expression is associated with adventitious rooting. For example, the expression of the beta HLH protein 93 (Trb7468921), early-responsive to dehydration stress protein (ERD4) (Trb4932701), ethylene responsive element binding factor 1 (Trb5562722), peroxidase 2 (Trb4932701), and zinc-finger protein 1 (Trb8690271) DEGs increased continuously in the IBAtreated cuttings. In addition, the expression of the C2H2like zinc finger protein (Trb5430271), cold, circadian rhythm, and RNA binding 1 (Trb6574322), embryospecific protein 3 (Trb4689243), and heat shock protein 60 (Trb7645302) DEGs increased continuously in the CK cuttings. Moreover, the expression of the ACC oxidase 1 (Trb7020711), basic helix-loop-helix (bHLH) DNAbinding protein 5 (Trb8657204), cell wall/vacuolar inhibitor of fructosidase 2 (Trb5894705), LOB domaincontaining protein 41 (Trb6654701), and CAP160 protein (Trb3457711) DEGs increased in both IBA-treated and CK cuttings.

DEGs in IBA effects on adventitious root development
Genes that were differentially expressed in the four developmental stages were identified using IDEG6. Genes were determined to be IBA-regulated if they had fold-change > 2 and P ≤ 0.05 in at least one rooting stage. In total, 14,924 DEGs were observed in the four developmental stages of the IBA-treated cuttings. The genes that exhibited differences in expression between two consecutive rooting stages are shown in Fig. 4. A comparison of tetraploid black locust cuttings in stage II and stage IC revealed that 8,976 genes were differentially expressed, of which 3,948 were down-regulated and 5,028 were up-regulated in stage IC (Fig. 4). Genes showing significant differential expression included 3 zinc finger domain-containing proteins (Trb5503102, Trb5503104, Trb5503101), 2 unknown proteins (Trb7570401, Trb2435501), embryonic abundant protein USP92 (Trb7595601), squamosa promoterbinding-like protein 13-like (Trb4932701), PEBP family protein (Trb7547281), and 4 auxin-induced proteins (Trb6569103, Trb5697601, Trb5697602, Trb5476601). A total of 1,677 DEGs were differentially expressed between stage IC and stage IRP, of which 643 were down-regulated and 1,034 were up-regulated in stage IRP. Genes showing significant differential expression included 2 hypothetical proteins (Trb5933302, Trb5933301), nodulin-26 (Trb14 811401), bidirectional sugar transporter SWEET3-like (Trb6691501), aquaporin TIP1-1 (Trb2089901), a peptide/ nitrate transporter (Trb6088201), a mitochondrial alternative oxidase (Trb6400801), momilactone A synthase-like (Trb5040301), sugar transport protein 13-like, zeatin Oglucosyltransferase-like, and UDP-glycosyltransferase 74B1-like (Trb4128302, Trb3408601, Trb6028401). A total of 4,383 genes were differentially expressed between the IRP stage and the final stage (IAR), of which 1,576 were down-regulated and 2,707 were up-regulated in stage IAR. Genes showing significant differential expression included 2 protease inhibitor-like proteins (Trb5017202 and Trb5017201), 4 disease resistance proteins (Trb2657401, Trb7602201, Trb4619901, Trb4619902), 6 chlorophyll a-b binding protein CP26, chloroplastic-like isoforms (Trb56455601, Trb6596903, Trb6596905, Trb6596902, Trb5674201, Trb6596911), S-adenosylmethionine decarboxylase proenzyme belonging in spermine biosynthesis (Trb4357111, Trb43 57101, Trb4357121), S-adenosylmethionine synthase  Fig. 2 The species distribution of unigene blastx results against the NCBI-Nr protein database and methionine synthase in the S-adenosylmethionine biosynthetic process (Trb48279501, Trb80561501, Trb865153011), 1-aminocyclopropane-1-carboxylate synthase in the ethylene biosynthetic process (Trb34955101), auxin response factor 18-like (Trb20505701), and 6 abscisic acid receptor PYL6-like proteins (Trb6836312, Trb6836315, Trb6653601, Trb6653602, Trb6836313, Trb6833901). Moreover, 6,877, 549, and 2,557 genes showed specific regulation only between stages II and IC, IC and IRP, and IRP and IAR, respectively (Additional file 5). Additionally, 581 DEGs were common between stages II vs. IC and stages IC vs. IRP, 178 DEGs were common between stages IC vs. IRP and stages IRP vs. IAR, and 1,151 DEGs were common to stages II and IC and stages IRP and IAR (Additional file 5). The comparison of the four consecutive developmental stages showed that 368 DEGs were common to all four stages. Among them, ten and six genes were continuously up-regulated and down-regulated, respectively, in all four developmental stages. These genes are mainly involved in methionine metabolism pathways. Several highly induced genes, such as ACC oxidase (Trb7020711), a zinc finger family protein (Trb6191801) and a wound-responsive protein (Trb74240111) were expressed most highly during the stages IRP and IAR (Additional file 5). The functional classification of the differentially expressed unigenes in the four developmental stages is shown in Additional file 6. These results indicate that these genes may promote rooting in tetraploid black locust cuttings. The DEGs included more up-regulated transcripts than down-regulated transcripts, indicating that many genes responded positively to IBA treatment. This result is consistent with previous studies in Arabidopsis and tomato [27][28][29].
DEGs in response to IBA during the four developmental stages of adventitious rooting Analysis using the Z-score method suggested that in each of the four developmental stages, the expression levels of 16,325 genes differed significantly between the transcriptomes of IBA and CK at P < 0.05 and foldchange > 2. Of these, 2570 unigenes, 1721 up-regulated and 849 down-regulated, were differentially expressed in CI vs. II; 1,685 unigenes, 849 up-regulated and 836 down-regulated, were differentially expressed in CC vs. IC; 1,512 unigenes, 881 up-regulated and 631 downregulated, were differentially expressed in CRP vs. IRP; and 10,558 unigenes, 5,626 up-regulated and 4,932 down-regulated, were differentially expressed in CAR vs. IAR (Fig. 4). In addition, 10,588 genes (64.86% of the total) showed differential expression between CK and IBA treatment only during the last stage. The identification of more differentially expressed genes in the final stage might be related to the greater distinction of this developmental stage. Moreover, in all four developmental stages, there were more up-regulated genes than down-regulated genes. These results suggest that IBA treatment increased the number of genes that were upregulated to promote adventitious rooting.

GO enrichment analysis
To discern global patterns of differential transcript abundance over the time course, unigenes with contrasting significance at P ≤0.05 were further filtered to include only those with values greater than fold-change > 2 in a comparison of unigene abundance between the samples [32]. For the groups of up-regulated and downregulated genes, we applied WEGO to compare the GO classifications of these genes [33]. Comparing the IBAtreated cuttings to the CK cuttings at each of the four developmental stages, the results showed significantly more up-regulated GO classifications than downregulated GO classifications at all developmental stages (Fig. 5). Further, when comparing pairs of consecutive The GO analysis of the DEGs in the four developmental stages between IBA-treated and CK revealed that most of the encoded products were associated with the following GO categories: regulation of transcription, transcription of DNA dependent, integral to membrane and ATP binding. The most common categories associated with the AR stage were integral to membrane and ATP binding activity (Additional file 8). Several genes showed highly significant differences, including cytokinin dehydrogenase (Trb5994202), 1-aminocyclopropane-1-carboxylate oxidase (Trb7341803), 4-coumarate: CoA ligase-like 9 (Trb6655606), ubiquitin protein ligase RGLG1-like (Trb5586102), S-adenosylmethionine synthase (Trb5508801 and Trb9018001), spermidine synthase, S-adenosylmethionine decarboxylase (Trb6300002 and Trb4357121) and zeaxanthin epoxidase (Trb72 43506). These auxin-responsive genes have important functions in AR formation (Additional file 9). In addition, highly significant differences in gene expression were detected, most of which were distributed among the following GO terms: cytokines in metabolic process, ethylene biosynthetic process, jasmonic acid biosynthetic process, auxin metabolic process, S-adenosylmethionine biosynthetic process, spermidine biosynthetic process and abscisic acid biosynthetic process.
The top 50 most significantly up-and down-regulated GO categories during the IBA-treated development process are listed in Additional file 10. The results revealed that the GO classifications associated with organism development, such as AR development, xylem development, phloem development, post-embryonic root development, and organ development, were significantly up-regulated in CI compared with II but significantly down-regulated at stage IRP compared with IAR. Hormone-related pathways, such as ethylene-mediated signaling pathway and auxin metabolic process, were significantly up-regulated in CI compared with II, whereas genes with the GO classifications abscisic acid biosynthetic process, jasmonic acid biosynthetic process and cytokinin metabolic process were significantly upregulated at stage IRP compared with IAR.

KEGG pathway enrichment analysis
To further determine which biological pathways were significantly (P ≤ 0.05) modulated during AR formation, KEGG pathway enrichment was performed using the KEGG Automatic Annotation Server (KAAS) [34] to reveal KEGG pathway enrichment in the transcriptomes of IBA-treated and CK cuttings at four developmental stages. Between stages CI and II, thirteen KOs were significantly down-regulated or up-regulated, including glycolysis/gluconeogenesis, cysteine and methionine metabolism, photosynthesis, amino sugar and nucleotide sugar metabolism, lysosome, and starch and sucrose metabolism. Between stages CC and IC, five KOs were significantly down-regulated or up-regulated, including alanine, aspartate and glutamate metabolism, bacterial secretion system, butanoate metabolism, and starch and sucrose metabolism. Further, when comparing stages CRP and IRP, seven KOs were significantly downregulated or up-regulated, including Alanine, aspartate and glutamate metabolism, beta-Alanine metabolism, and Type I diabetes mellitus. Meanwhile, comparing stages CAR and IAR, twenty-five KOs were significantly down-regulated or up-regulated, including glycolysis/ gluconeogenesis, carbon fixation in photosynthetic organisms, alzheimer's disease, starch and sucrose metabolism, cysteine and methionine metabolism, arginine and proline metabolism, MAPK signaling pathway, ABC transporters, valine, leucine and isoleucine degradation, antigen processing and presentation, nitrogen metabolism, tryptophan metabolism and selenoamino acid metabolism (Additional file 11). These results suggest that significant metabolic changes occur during the period of AR formation.
To gain insight into the differential KOs specifically induced in cuttings by IBA at the four developmental stages, the pairs of consecutive stages in the IBA-treated group were compared: II vs. IC, IC vs. IRP, and IRP vs. IAR. Forty KOs were significantly down-regulated or upregulated between stage II and stage IC, including glycolysis/gluconeogenesis, carbon fixation in photosynthetic organisms, pentose phosphate pathway, pyruvate metabolism, amino sugar and nucleotide sugar metabolism, cysteine and methionine metabolism, and alanine aspartate and glutamate metabolism. These results indicate that glycolysis/gluconeogenesis and carbon fixation in photosynthetic organisms were significantly upregulated by IBA treatment in stage IC relative to stage II. In stage IRP compared with stage IC, the KO Phenylalanine metabolism was significantly down-regulated, and the group Starch and sucrose metabolism was significantly up-regulated. Phenylpropanoids contribute to plant defenses as inducible chemical barriers or as signaling molecules [35][36][37]. From stage IRP to stage IAR, fifteen KOs were significantly down-regulated or upregulated, including alanine, aspartate and glutamate metabolism, carbon fixation in photosynthetic organisms, glycolysis/gluconeogenesis, arginine and proline metabolism, cysteine and methionine metabolism, pentose phosphate pathway, oxidative phosphorylation, and photosynthesis. cysteine and methionine metabolism, as well as ethylene pathway and associated polyamines, might play important roles in IBA-induced adventitious rooting [38]. Furthermore, the results further suggest that IBA increased the number of KOs and genes exhibiting changes in expression during the early stages of rooting.

Verification of DEGs during the four developmental stages of adventitious rooting
RNA was extracted from IBA-treated and CK cuttings at the four selected stages of AR development and used as the template for q-PCR-based validation of the sequence-based transcription profiles of 21 differentially expressed candidate unigenes. Detailed information on these genes is presented in Additional file 12. The selected genes were associated with methionine metabolism, plant hormone signal transduction, phenylalanine metabolism and other enzymatic processes (Fig. 6). Linear regression [(RNA-seq value) = a(q-PCR value) + b] analysis revealed an overall correlation coefficient of 0.71 according to Villacorta-Martín et al. and Yu et al. [39,40]. The q-PCR analysis confirmed that the RNAseq approach provided reliable data regarding differential gene expression during the AR developmental stages of tetraploid black locust cuttings. The genes threonine kinase (Trb6720901), auxinrepressed protein (Trb2435601), ethylene responsive transcription factor (Trb6253801), auxin-induced in root cultures protein (Trb6310802), auxin-responsive protein IAA (Trb8029901), and SAMDC (Trb7651101) beta-Dxylosidase (Trb6237301) remained highly expressed during AR development for all cutting stages, and the expression of these genes was higher in IBA-treated cuttings than in CK cuttings in all four stages. We observed the highest expression of flowering promoting factor protein (Trb5566601) in the callus induction phase of IBA-treated cuttings, as well as high expression in the AR formation phase, and the expression of this gene in the AR formation phase was higher in IBA-treated cuttings than in CK cuttings.
We observed the highest expression levels of heat shock cognate protein (Trb5873209), zinc finger protein (Trb6191801), ACO (Trb7020711), heat shock cognate protein (Trb5873209), SAMS (Trb6301201) woundinduced protein (Trb7602201), and aspartic protease in guard cell (Trb7014201) in the AR formation phase of IBA-treated cuttings, and the expression of these genes in the AR formation phase was higher in IBA-treated cuttings than in CK cuttings. The most highly expressed protein detected was phloem protein (Trb6463501) in the root primordia formation phase of IBA-treated cuttings, as well as high expression in the AR formation phase, and expression of this gene in the AR formation phase was higher in IBA-treated cuttings than in CK cuttings. We observed the highest expression levels of Phloem WRKY transcription factor (Trb6276401), 2,4-D inducible glutathione S-transferase (Trb5981501), and peroxidase (Trb6791601) during the initiation formation phase of IBA-treated cuttings, as well as high expression in the AR formation phase, and expression of these genes in the AR formation phase was higher in IBAtreated cuttings than in CK cuttings. These results demonstrate that IBA might directly or indirectly regulate the expression of the above genes during AR development in tetraploid black locust. Several members of the threonine kinase, auxin-repressed protein, ethylene responsive transcription factor, auxin induced in root cultures protein, auxin-responsive protein IAA, ACO and SAMDC beta-D-xylosidase families have been identified and shown to mediate adventitious rooting [22,41]. However, during the AR formation phase, the expression levels of the genes cytochrome p450 (Trb5844701), Delta-l-pyrroline-5-carboxylate synthetase (Trb6734101), and transcription factor bHLH (Trb6501202) were higher in CK cuttings than in IBA-treated cuttings. These results show that these genes were less affected by IBA treatment during the stages of AR development. Instead, cytochrome p450 and transcription factor bHLH genes are involved in responses to stress, such as drought and high salinity, thus leading to efficient adventitious rooting [11].

Unigene determination of transcriptome sequencing during AR development in tetraploid black locust cuttings
Illumina RNA-seq technology has been extensively used for model plant transcriptome sequencing [42] with reference genome data or for non-model plants [43] without reference genomic information. In this study, the Illumina HiSeq 2000 platform was used to perform a de novo transcriptome sequencing analysis of the tetraploid black locust cuttings to better understand the gene expression changes during adventitious rooting. Pooled RNA samples from IBA-treated and CK cuttings sampled at four time points after AR excision were used to construct cDNA libraries for deep sequencing. In this sequencing, approximately 33.74 million to 46.81 million paired-end clean reads were obtained from the IBAtreated and CK cuttings at the four time points. After de novo assembly, we obtained 101,209 unigenes with a mean length of 986 bp, which is longer than has been reported previously in studies using the same technology [26,28,43,44].
We identified a total of 10,181 and 14,924 DEGs (fold change > 2) were detected in at least one of the four stages in the CK and IBA treated cuttings, respectively. The genes that exhibited differences in expression between two consecutive rooting stages are shown in Fig. 4. A comparison of tetraploid black locust cuttings in stage II and stage IC revealed that 8,976 genes were differentially expressed, of which 3,948 were down-regulated and 5,028 were up-regulated in stage IC. Using a DNA microarray method, Rigal et al. [30] studied gene expression changes during adventitious rooting in the model Fig. 6 q-PCR validation of differential expression. Transcript levels of 21 genes in CK (black column) and IBA (gray column). The y-axis shows the relative gene expression levels as analyzed by q-PCR. Bars represent the standard error (n = 3). A Comparison of the gene expression ratios obtained from RNA-seq data and q-PCR tree Populus trichocarpa. Their results indicated that 5,781 genes were differentially expressed in the organization of the AR primordium; 6,538 genes were differentially expressed during primordium differentiation; and 1,146 genes were differentially expressed between these two stages [9]. In another similar study using cDNA microarrays, Brinker et al. [17] identified 220 genes whose expression changed significantly during root development in hypocotyl cuttings of Pinus contorta [6]. The results obtained suggest that RNA-Seq is a sensitive, low-cost, and accurate method for deepsequencing the transcriptomes of plant without available genomic information, and this method was able to identify more DEGs during the early stages of adventitious rooting relative to the results of DNA microarrays. This technology also enables the precise elucidation of transcripts in the samples.
The associations between the GO terms and the lists of DEGs were investigated using GO functional enrichment analysis. Significantly enriched GO terms for the DEGs in the II vs. CI, IC vs. CC, IRP vs. CRP and IAR vs. CAR comparisons included cellular component, cell part, membrane and membrane-bound organelle (Fig. 5). Among biological processes, the highest number of unigenes belonged to cellular metabolic process, and the top three classes of genes were primary metabolic process, macromolecular metabolic process and response to stimulus. Additionally, the molecular function and cellular metabolic process terms contained the highest numbers of unigenes, followed by metabolic process, macromolecular metabolic process and response to stimulus. The GO functional enrichment analysis and unigene expression abundance data will provide useful information for the identification of genes involved in AR development in tetraploid black locust.
Genes involved in cysteine and methionine metabolism were significantly regulated by IBA during AR formation We further examined the genes encoding proteins involved in Cysteine and methionine metabolism during the process of adventitious rooting. We identified a number of unigenes (fold change > 2), including 12 methionine-related genes, and 12 ethylene-related genes, associated with the KEGG pathway Cysteine and methionine metabolism (Additional file 11). Among those 62 genes, 60 were identified as auxin-related. The genes down-regulated at stage II included a threonine synthase gene and four S-adenosylmethionine decarboxylase genes; between stages II and IC, the down-regulated genes included six S-adenosylmethionine synthetase-like genes and an aspartate kinase-like gene. Compared with stages II and IC, a total of 21 genes, including 11 upregulated and 10 down-regulated, were identified as differentially regulated at stages IRP and IAR. The up-regulated genes were primarily members of the Sadenosylmethionine synthetase family, while the down-regulated genes were mostly members of the adenosylhomocysteinase family. S-adenosylmethionine decarboxylase has been known to function as an auxin carrier complex in cellular auxin efflux and influx [44]. These results indicate that most of the genes related to methionine metabolism were significantly up-regulated by IBA treatment. The upregulation of aspartate aminotransferase and the DNA methyltransferase during stages II and IC suggests that auxin transport occurs in these stages. In other studies, the expression levels of spermidine synthase and S-adenosylmethionine decarboxylase, which are essential for AR formation [45], were up-regulated by IBA treatment. Among the ethylene-related genes, 17 DEGs were identified at stage IRP, with 15 upregulated and two down-regulated, and 11 DEGs were identified at stage IAR, with eight up-regulated and three down-regulated. Compared with stages II and IC, stages IRP and IAR showed a total of 20 DEGs, with 9 up-regulated and 11 down-regulated.

The DEGs in each developmental stage in the CK and IBA treatments
To further understand the roles of the DEGs in each developmental stage in the CK and IBA-treated cuttings, we further analyzed the expression levels of a total of 25 genes that showed differential expression in the samples. Among these, several genes were differentially expressed in the AR formation phase, with higher expression in IBA-treated cuttings than in CK, including ACC oxidase 1 (ACO1). Interestingly, the expression of these genes is induced by exogenous auxin in the rooting-competent cuttings of two distantly related forest species, and there are some similar reports that ACO genes lead to AR formation in mung bean and tomato [46,47]. It appears that IBA-induced ethylene production may contribute to the stimulation of adventitious rooting [47]. Arabinogalactan protein 22 (AGP22) and other AGPs are extracellular proteoglycans that are implicated in many plant growth and developmental processes; for example, AtAGP30 is a non-classical AGP core protein from Arabidopsis that is expressed only in roots [48]. Laskowski et al. report that FQR1 is a novel primary auxin-response gene that encodes a flavin nucleotidebinding flavodoxin-like quinone reductase, and accumulation of FQR1 mRNA can be detected in roots and root cultures [49]. AtMRP5, as a newly identified member of the ABC transporter superfamily, controls root development in Arabidopsis thaliana, according to Gaedeke et al. [50]. Further, Zhang et al. and To et al. have also reported that the genes ARR8/ARR9 are expressed in the root and that they act partially redundantly to negatively regulate the response of roots to exogenous cytokinin during root meristem initiation [51,52]. Plants possess a superfamily of arm-repeat proteins, which were reported to be involved in abscisic acid signaling by Kim et al. [53] and in lateral root development by Coates [54].
The DEGs detected encompassed ten gene families. Among the auxin-responsive factor (ARF) family proteins, ARF2 and ARF5 were up-regulated both in CI compared with II and in CAR compared with IAR. Two auxin-responsive factor genes have been found to be involved in auxin signaling and to regulate adventitious rooting [55]. Zinc finger family proteins, specifically of the basic leucine zipper and C2H2 types, were up-regulated in CI compared with II and in CRP compared with IRP. Similar reports have indicated that basic leucine zipper genes promote the modulation of meristems and primordial development. The two wound-responsive family genes and two WRKY transcription factor family genes were significantly down-regulated in CAR compared with IAR and were significantly regulated by IBA. However, a conflicting report by Li et al. indicated that zinc finger proteins and WRKY were down-regulated by IBA treatment during the early stages of adventitious rooting [47]. These results indicate that these differentially expressed genes were directly or indirectly regulated by IBA.

Conclusion
Transcriptome sequence data for CK and IBA-treated tetraploid black locust cuttings at four developmental stages were obtained using the Illumina sequencing method, with subsequent de novo assembly. Despite the economic importance of this tree, the tetraploid black locust genome is not publicly available, and sequence data for tetraploid black locust are limited. Our study generated the first large-scale transcriptome dataset of tetraploid black locust for CK and IBA-treated cuttings. Additionally, the types and quantities of the genes expressed in CK and IBA-treated tetraploid black locust cuttings, along with their functions, classifications, and metabolic pathways, were revealed for the first time. In total, 10,181 and 14,924 DEGs were detected in CK and IBA-treated cuttings, respectively. Our study presents a global view of transcriptome dynamics and differential gene expression analysis for IBA-treated tetraploid black locust cuttings over four developmental stages and provides new insights into the fundamental mechanisms associated with auxin-induced adventitious rooting. Our data constitute a valuable resource for genomic investigations of AR formation in tetraploid black locust cuttings for improving rooting of difficultto-root varieties.

Plant materials and IBA-treated adventitious root development
Cuttings of tetraploid black locust (Robinia pseudoacacia L.) were collected from a 1-year-old seed bearer at the nursery of Northwest Agriculture and Forestry University, Yangling, China. The preparation, experimental plan and IBA-treated cuttings were performed according to the methods described by Quan et al. [18,19]. Subterminal parts of stems 15 cm in length and 10-12 mm in diameter were collected from the same cultivar. The basal 2.5 cm of each cutting was then placed for 4 h in water for the CK or in 5.4 mM IBA as the auxin treatment. The cuttings of tetraploid black locust were subsequently placed on a bench in a glasshouse, and 5 cm portions of the basal parts of cuttings were buried in sand. Cuttings were managed and randomly selected according to the method described by Quan et al. [20]. The 2-cm-long basal stem regions of each cutting, where ARs develop, were cut and collected separately at 0, 15, 20, and 25 days after planting, and the harvested samples were designated control cuttings sampled at stages CI, CC, CRP and CAR for the CK and at stages II, IC, IRP and IAR for the IBA-treated cuttings. Ten cuttings randomly selected from each group were pooled. For each treatment and time point, three biological replicates, each consisting of cuttings from ten base stems, were collected, frozen immediately in liquid nitrogen, and stored at −80°C prior to RNA extraction.

RNA separation, cDNA library preparation and Illumina sequencing
Using the manufacturer's protocol (TIANGEN Biotech, Beijing), the RNA prep Pure Plant Kit was applied to isolate total RNA from each sample. To prevent contamination with genomic DNA, RNA samples were processed with RNase-free DNase I (Takara, Japan). We used an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) to obtain and quantify RNAs. Thereafter, we used ethidium bromide staining and denatured agarose gel electrophoresis, to inspect the completeness of RNAs. Subsequently, the analysis only included RNA integrity numbers ≥7.5, RNA 28S:18S ratios higher than one and RNA samples with A260/A280 ratios between 1.9 and 2.10.
Based on the manufacturer's guidelines, Dynal oligo (dT) 25 beads were used to separate poly(A) + RNA from total RNA of tetraploid black locust in order to conduct Illumina sequencing. We purified the poly (A) + RNA and then cleaved the mRNA into short pieces in a cleavage buffer. Using the random hexamer primer N6, SuperScript III reverse transcriptase and short pieces as templates, we synthesized the first-strand cDNAs. Next, RNase H, buffer, DNA polymerase I and dNTPs were used for second-strand cDNA synthesis. End repair of the double-stranded cDNAs was performed using T4 polynucleotide kinase, T4 DNA polymerase and DNA polymerase I Klenow fragment. Then, T4 DNA ligase was used to ligate the fragments to adapters. The Qia-Quick PCR extraction kit was used to isolate the ligated fragments (200 ± 25 bp), which were then eluted in EB buffer and analyzed using agarose gel electrophoresis. Then, we chose appropriate segments as samples for PCR amplification. Eventually, we used an Illumina HiSeq™ 2000 to sequence the library at LC Sciences in Hangzhou, China. All technical steps were performed twice.

De novo assembly and functional annotation
We discarded empty reads, adapter sequences, and poor sequences from the clean raw reads (with unknown sequences (N) or sequences less than 25 bp). Next, short reads were processed in triplicate for de novo assembly of transcriptomes. Then, the triplicate samples were used to collect clean reads into non-repeating transcripts. To enhance the quality of the reads, we also deleted short sequences (<200 bp). Using the Kyoto Encyclopedia of Genes and Genomes (KEGGs) database with an E-value cut-off of 10 −5 and the Clusters of Orthologous Groups of proteins (COGs), Swiss-Prot, Pfam and NCBI non-redundant protein (Nr) databases, we identified the most important sequences through BLAST searches and annotation. We used Blast2go software to perform functional annotation of Gene Ontology (GO) (www.geneontology.org).

Identifying differentially expressed genes
We applied RPKM (reads per kilobase per million reads), as the standard value of gene expression [56], to analyze differential gene expression. Using the IDEG6 web tool (http://telethon.bio.unipd.it/bioinfo/IDEG6_form/) [57], we performed statistical comparisons between the CK samples and the IBA-treated samples in terms of RPKM values. In many experiments, the importance of changes in gene expression was assessed based on the P-value, limiting the false discovery rate (FDR) at 0.01. In comparisons between libraries, genes showing fold-change >2 with P ≤0.05 were defined as differentially expressed. For up-and down-regulated unigenes, we used WEGO to compare their GO classifications [33].

Real-time quantitative PCR (RT-qPCR)
For RT-qPCR, we chose 21 candidate genes (Additional file 12) that showed differential expression to verify the Illumina sequencing data. Reverse transcription was performed using the PrimeScript RT reagent kit (Takara, Dalian, China) and 1 μg of RNA. The products of reverse transcription were diluted to 150 ng/μl in double-distilled water. Primer Premier 6.0 was used to generate antisense and sense primers. Using TaKaRa SYBR Premix Ex Taq II (Perfect Real Time), we performed RT-qPCR in a 20 μl reaction volume with a RC96 Real-Time PCR detection system (Rocipe Laboratories, Hercules, CA, USA). The 18S rRNA gene was used as the endogenous reference gene. After 2 min of predenaturation at 94°C to measure fluorescence, RT-qPCR was performed with forty-five cycles of denaturation at 94°C for 10, annealing at 60°C for 15 s and elongation at 72°C for 30 s. To ensure the absence of primer dimers, the melting curve was analyzed in the first cycle. The 2 -△△Ct method [14] was used to confirm the relative expression of target genes. Each sample was analyzed in triplicate. We used Origin (Version 8.0) to analyze the data and then log-transformed all data obtained from the q-PCR analysis.