The rooting phenotype of blueberry green cuttings varied significantly between control and IBA treatment (Fig. 1a-b). The rooting percentage and average root number per cutting were both significantly promoted by exogenously-applied IBA (Fig. 1c-d). There was no significant difference in average root length between these two treatments (Fig. 1e), suggesting that exogenous auxin IBA could enhance the adventitious root formation without influencing their average lengths.
Microstructure observation of adventitious root formation
To clear and definite the AR developmental process against the stem of blueberry green cuttings, the microstructure was observed. The results showed that callus tissues started to form at 7 days after cutting (DAC) in IBA treatment (Fig. 2a), then adventitious root primordium initiation (rpi) was produced at 14 DAC (Fig. 2b), and developed into adventitious root primordium (rp) at 21 DAC (Fig. 2c), then rp was developed to adventitious root (AR) (Fig. 2d) and finally outgrowth from the stem at 35 DAC (Fig. 2e), suggesting the blueberry AR formation was initiated from non-root pericycle cells as a type of auxin-induction.
Analysis of plant hormones in blueberry green cuttings during ARs formation
Changes of plant hormone during blueberry AR formation were analyzed. The level of IAA exhibited a significantly increase at 21 day after cutting (DAC) under IBA treatment (Fig. 3a-b), which might indicating that the higher IAA might contribute to initiate adventitious root primordium (rpi). No significant difference in GA3 content was observed in green cuttings between control and IBA treatment (Fig. 3c). The content of ABA and zeatin under IBA treatment was kept at a lower level during the whole trail, and a decreasing tendency of ABA and zeatin was recorded in control treatment and got a similar level as IBA treatment after 21 DAC (Fig. 3d), suggesting that a lower level of ABA and zeatin may facilitate the ARs formation.
De novo sequencing, Assembly, and Gene Annotation in RNA-seq of blueberry stem
The basal stem of cuttings treated with or without IBA were sampled every 7 days after insertion, and ten samples (CK7, CK14, CK21, CK28, CK35, T7, T14, T21, T28 and T35) were subjected to total RNA extraction and RNA-seq analysis. High-throughput sequencing generated 41.45 - 44.46 million (M) pairs of 150 bp raw reads from each library. After a stringent quality filtering process, 212 million clean reads (99.54% of the raw data) remained, with a Q30 percentage (an error probability lower than 0.1%) ranging from 92.53% to 94.00%. The number of clean reads per library ranged from 41.27 to 44.25 M (Table 1). The total clean reads were de novo assembled into transcripts by trinity, and 672606 transcripts and 308719 unigenes were assembled with an average length of 735.87 bp (N50=1082) and 617.30 bp (N50=844), respectively (Table 2). All unique sequences were further annotated based on Blastx searches with a cut-off E-value of 10-5 against five databases, total unigenes of 107111 (34.70%), 53859 (17.45%), 6696 (2.17%), 101729 (32.95%) and 89910 (29.12%) were respectively annotated by NCBI non-redundant (NR) database, gene ontology (GO) database, kyoto encyclopedia of genes and genome (KEGG) database, evolutionart genealogy of genes: Non-supervised orthologous groups (eggNOG) database and Swiss-Prot database (Table 3).
Based on the NR annotations, 25.33% of the annotated sequences had strong homology (E-value < 10-60), 21.80% of the annotated sequences showed strong homology (10-60< E-value < 10-30), and an additional 82.87% of the annotated sequences showed homology (10-30< E-value < 10-5) to available plant sequences (Additional file1: Fig. S1a). The similarity distribution was comparable, with 26.41% of the sequences having similarities higher than 80%, while 73.5% of the hits had similarities of 0-80% (Additional file1: Fig. S1b). With respect to species, 72.99% of the unique genes showed high matches with sequences form other species, followed by Vitis vinifera (9.3%), Oryza sativa Japonica Group (5.74%) and Coffea canephora (2.83%) (Additional file1: Fig. S1c).
GO analysis was also performed with the Blast2GO software. Out of 308719 unigenes, 263143 were classified into the “biological process”, “cellular component” and “molecular function” categories (Additional file2: Fig. S2). This classification provided some information on the percentage of blueberry unigenes in different signal transductions, catabolic and anabolic process. For the biological process category, the majority of unigenes was grouped into “metabolic process” (27.17%), “cellular process” (27.63%) and “single-organism process” (18.03%), which accounted for about 72.83%. In the cellular component category, the unigenes were mainly distributed into “cell” (21.20%), “cell part” (21.06%), “membrane” (15.31%), “membrane part” (10.95%) and “organelle” (15.20%), accounting for about 83.72%. For the molecular function category, a large number of unigenes was distributed into “catalytic activity” (10.43%) and “binding” (10.09%), which accounted for about 86.26% (Additional file2: Fig. S2).
Ten DEGs libraries were further analyzed based on RSEM quantitative software, and the FPKM (fragments per kb per million fragments, FPKM) of all unigenes were calculated at each sampled period in blueberry between control (CK) and IBA treatment (T). Differences in gene expression was examined using the threshold of |log2FoldChange| > 1 at P-value ≤ 0.05. The DEGs were identified by pairwise comparisons of ten libraries, i.e. CK7 vs. T7, CK14 vs. T14, CK21 vs. T21, CK28 vs. T28, and CK35 vs. T35 (Additional file3: Fig. S3a). A total of 14970 DEGs were detected between CK and T libraries, of which there were 7467 up-regulated and 7503 down-regulated genes (Additional file3: Fig. S3b). For each sampled period, 3252 DEGs was detected between CK7 and T7 libraries, these unigenes were directly affected by IBA treatment and might be associated with callus tissue formation (Fig. 2a, Additional file3: Fig. S3b). In CK14 vs. T14, there are 3999 DEGs were indentified, which might be related with rpi formation (Fig. 2b, Additional file3: Fig. S3b). In CK21 vs. T21, total 1488 DEGs was determined, which were contributed to rp formation (Fig. 2c, Additional file3: Fig. S3b). In CK28 vs. T28, 2648 DEGs were indentified, which were contributed to AR formation (Fig. 2d, Additional file3: Fig. S3b). And, in CK35 vs. T35, total 3583 DEGs were detected, which might be contributed to AR outgrowth and development (Fig. 2e, Additional file3: Fig. S3b). However, there were no commonly up-regulated or down-regulated DEGs at all sampled periods, as illustrated in the Venn diagram (Additional file3: Fig. S3c-d), suggesting DEGs might play special roles during AR formation, outgrowth and development progress.
DEGs enriched in the auxin signaling pathway
Auxin was proved to play key roles in promoting AR formation, therefore the DEGs in auxin-signaling pathway were annotated and further analyzed. A total of 29 auxin-related DEGS were mapped (Fig. 4). Of these DEGs, there were ten unigenes belongs to auxin-responsive protein including four ARFs and six SAURs family. Six auxin transporter-like genes were also identified including three influx carriers AUX22, AUX-LIKE 3 (LAX3), LAX5 and three exflux carriers PIN-LIKE 6 (PIL6s). All these DEGs were up-regulated during AR development in some extent (Fig. 4). To verify the result of comparative transcription analysis, two auxin responsive factor ARFs family (i.e. ARF7and ARF9) and six auxin transporter genes (i.e. AUX22, LAX3, LAX5 and PIL6a-6c) were selected to identify the DEGs expression profiles by qRT-PCR. RNA_seq showed that ARF7 and ARF9 exhibited a significant transcript peak at 28 DAC in IBA treatment, while both these two ARFs showed a decrease trend in control (Fig. 5a, Additional file4: Fig. S4). RNA_seq data revealed that AUX22 was kept at a lower expression in control, but it showed a significant up-regulation at 7 DAC – 21 DAC in IBA treatment (Fig. 5b). Its homologous genes LAX3 and LAX5 showed a decrease tendency in control, it was significantly un-regulated at 14 DAC and 28 DAC in IBA treatment (Additional file5: Fig. S5), suggesting these genes might play an important role in auxin transporting and contributed to auxin asymmetrical distribution which was facilitated to induce the adventitious root primordium initiation at early stage of AR formation. As for three auxin efflux carriers PIL6a-6c, in control treatment these three genes was low expression at 7 DAC – 21 DAC, and they were up-regulated until 28 DAC, whereas when applying IBA, PIL6a was significantly up-regulated at 14 DAC (Fig. 5c), PIL6b and PIL6c showed a remarkable increase from 7 DAC and peaked at 35 DAC (Additional file6: Fig. S6), suggesting auxin enhance expression of PIL6s and thus promote the AR formation. qRT-PCR data indicated a similarity with transcript information.
Transcription factor LATERAL ORGAN BOUNDARIES DOMAIN/LOB domain-containing proteins (LBDs) were suggested to regulate AR formation as downstream target genes of ARFs family. In present work, thirteen homolog genes of LBDs were identified based on NR annotation (Fig. 4), RNA-seq data revealed that all these LBDs were up- or down-regulated at different stage in a certain some extent, especially at 28 DAC. And, 4 out of these 13 LBDs (i.e. LBD16, 23, 29 and 37) were representatively selected to confirm their expression profiles by qRT-PCR, it indicated that LBD16 and LBD37 showed a good reproducibility with RNA-seq data, while qRT-PCR indicated LBD23 exhibited another up-regulated at 14 DAC in IBA treatment, LBD29 showed a continuous up-regulation during 21 DAC - 35 DAC in IBA treatment and were also up-regulated at 14 DAC and 35 DAC in control (Fig. 5d, Additional file7: Fig. S7).
DEGs involved in root primordium formation
Based on the NR annotation, six rooting-related DEGs were obtained including four homolog genes of lateral root primprdium 1 (LRP1), one putative root meristem growth factor 9 (RGF9), and one dormancy-associated protein homolog 3 (DRMH3) (Fig. 4). RNA_seq data showed that LRP1 was significantly up-regulated at 28 DAC (Fig. 5e), its homologous genes LRP1-like, LRP-type1, LRP-like2 exhibited a continuous up-regulation after 14 DAC and peaked at 28 DAC in IBA treatment (Additional file8: Fig. S8). An obvious up-regulation of RGF9 was observed at 7 DAC, 14 DAC and 28 DAC (Fig. 5f), DRMH3 showed a significantly up-regulation at 28 DAC in IBA treatment (Fig. 5g). The qRT-PCR analysis of RGF9 indicated a good consistency with RNA_seq data, the expression of DRMH3 in IBA treatment showed a good consistency with RNA_seq data, but it was continuously expressed at a high level in control treatment (Fig. 5g).
Putative gene regulatory networks that control blueberry AR formation
According to the known regulative networks reported previously in root formation in Arabidopsis and other plant species, the regulative pathway that controls blueberry AR formation was derived. It was speculated that IBA would induce the expression of auxin responsive factors ARF7/9 to percept auxin signaling, whereafter ARF7/9 directly or indirectly affected their downstream target LBDs genes to establish AR founder cells with nuclei migration, and then auxin polar carriers including influx carriers AUX22 or LAX3/5 and efflux carriers PIL6s would be up-regulated to facilitate the establishment of auxin asymmetric distribution and thus include the AR primordium formation, finally AR primordium transform to AR apical meristem and outgrowth form the cuttings under the effect of LRP1, RGF9, DRHM3 and other genes (Fig. 6).