Phenotypic analysis
The rooting phenotype of blueberry green cuttings varied significantly between control and IBA treatment (Figure 1a, 1b). The rooting percentage and average root number per cutting were both significantly promoted by exogenously-applied IBA (Figure 1c, 1d). There was no significant difference in average root length between these two treatments (Figure 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 (c) tissues started to form at 7 days after cutting (DAC) in IBA treatment, then adventitious root primordium initiation (rpi) was produced at 14 DAC, and developed into adventitious root primordium (rp) at 21 DAC, the rp was finally developed to adventitious root (r) and outgrowth from the stem (Figure 2), 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 (Figure 3a, 3b), 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 (Figure 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 (Figure 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 stem from base of cuttings 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 bp (N50=1082) and 617 bp (N50=844), respectively. 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. 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 (Figure 4a). 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% (Figure 4b). 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%) (Figure 4c).
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 (Figure 5a). A total of 14970 DEGs were detected between CK and T libraries (Figure 5b), of which there were 7467 up-regulated and 7503 down-regulated genes. 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. In CK14 vs. T14, there are 3999 DEGs were indentified, which might be related with rpi formation. In CK21 vs. T21, total 1488 DEGs was determined, which were contributed to rp formation. In CK28 vs. T28, 2648 DEGs were indentified, which were contributed to AR outgrowth. And, in CK35 vs. T35, total 3583 DEGs were detected, which might be contributed to AR development (Figure 5b, Figure 2). However, there were no commonly up-regulated or down-regulated DEGs at all sampled periods, as illustrated in the Venn diagram (Figure 5c), 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 (Figure 6). 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 (Figure 6). 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 (Figure 7, 8, 9). RNA_seq showed that ARF7 and ARF9 exhibited a significant transcript peak at 28 d in IBA treatment, while both these two ARFs showed a decrease trend in control (Figure 7). RNA_seq data revealed that AUX22 was kept at a lower expression in control, but it showed a significant up-regulation at 7d – 21d in IBA treatment (Figure 8). Its homologous genes LAX3 and LAX5 showed a decrease tendency in control, it was significantly un-regulated at 14 d and 28 d in IBA treatment (Figure 8), suggesting these genes might play an important role in auxin transporting and contributed to auxin asymmetrical distribution which were 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 – 21 DAC, and they was up-regulated until 28 DAC, whereas when applying IBA, PIL6a was significantly up-regulated at 14 DAC, PIL6b and PIL6c showed a remarkable increase from 7 DAC and peaked at 35d (Figure 9), suggesting auxin enhance expression of PIL6s and thus promote the AR formation. qRT-PCR data indicated a similarity with transcript information (Figure 7, 8, 9).
Transcription factor LATERAL ORGAN BOUNDARIES DOMAIN/LOB domain-containing proteins (LBDs) were suggested to regulate AR formation as downstream target genes of ARFs family [19, 20]. In present work, thirteen homolog genes of LBDs were identified based on NR annotation (Figure 6), 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 21d -35d in IBA treatment and were also up-regulated at 14d and 35d in control (Figure 10).
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). RNA_seq data showed that LRP1 was significantly up-regulated at 28 d, LRP1-like, LRP-type1, LRP-like2 exhibited a continuous up-regulation after 14d and peaked at 28 d in IBA treatment (Figure 11). An obvious up-regulation of RGF9 was observed at 7d, 14d and 28d (Figure 12), while DRMH3 showed a significantly up-regulation at 28 d (Figure 12). The qRT-PCR analysis of RGF9 and DRMH3 indicated a good consistency with RNA_seq data (Figure 11, 12).