6-BA treatment effects bud size and flowering intensity
Two-years of observations indicated that a notably higher percentage of flowering was evident on short shoots than medium or long shoots (Fig. 1a); and that a higher number of flowers were present in trees treated with 6-BA than in untreated, control trees (Fig. 1a, b). The length, width, and fresh weight of terminal buds were also significantly greater in the 6-BA-treated trees, relative to the control trees (Fig. 2a, b). Scanning electron microscopy of buds at 100 DAFB indicated that the 6-BA treatment did not affect the structure of the floral primordia (Fig. 2c, d). Collectively, the data indicated that the 6-BA treatment affected bud size and flowering intensity.
RNA-Seq analysis of 6-BA treated buds and control buds
The total number of raw reads obtained from each sample ranged from 44.71 million to 69.26 million. After filtering, the number of clean reads per library ranged from 43.5 million to 65.6 million. Greater than 86% of the clean reads of most samples could be uniquely mapped to the apple genome (Additional file 1: Table S1). Approximately 28,000 genes with an average FPKM ≥1 in at least one sample were considered to be expressed in each sample (Fig. 3a). A Pearson correlation analysis indicated that gene expression levels between the biological replicates were highly related, having an R 2 > 0.92 (Additional file 2: Figure S1). DEGs were identified from a comparison between 6-BA treated and control buds at each time point (Fig. 3b; Additional file 3: Table S2). A total of 8015 non-redundant DEGs were obtained by pairwise comparisons of B30 vs. C30, B50 vs. C50, and B70 vs. C70. In the B30 vs C30 comparison, 3803 were up-regulated and 1989 were down-regulated. In the B50 vs. C50 comparison, 1215 were up-regulated and 2035 were down-regulated; and in the B70 vs. C70 comparison, 1070 were up-regulated and 946 were down-regulated (Fig. 3b, c).
Twelve DEGs representing 6-BA responsive genes were selected for RT-qPCR analysis. The results indicated that the RNA-seq and RT-qPCR data were similar and produced similar patterns of expression, thus confirming the expression levels determined by RNA-seq (Additional file 4: Figure S2).
Functional classification of 6-BA responsive genes in apple buds during floral transition
The DEGs obtained from the B30 vs C30, B50 vs C50, and B70 vs C70 comparisons were subjected to KEGG pathway enrichment analysis (Additional file 5: Table S3). The most significantly enriched KEGG in the B30 vs. C30 DEGs was plant hormone signal transduction, following by fructose and mannose metabolism, galactose metabolism, metabolic pathways, biosynthesis of secondary metabolites, and starch and sucrose metabolism. Five KEGG pathways were identified to be significant in the DEGs obtained from the B50 vs. C50 comparison, including plant hormone signal transduction, photosynthesis, photosynthesis - antenna proteins, DNA replication, and plant-pathogen interaction. Lastly, three KEGG pathways were determined to be significantly enriched in the DEGs from the B70 vs. C70 comparison, including protein processing in the endoplasmic reticulum, carotenoid biosynthesis and plant hormone signal transduction.
Cluster analysis of 6-BA responsive genes during floral transition
Hierarchical cluster (H-cluster) analysis was performed to provide greater insight into the time-course profiles of the 8015 non-redundant DEGs. As a result, all of the DEGs were classified into 6 subclusters in both the 6-BA treatment and control groups (Fig. 4a, b; Additional file 6: Table S4). Genes with similar expression patterns are most likely regulated in a similar manner [31]. Therefore, the expression pattern of SOC1 homologs (Fig. 4c), which are well known key flowering pathway integrators [32], were further examined. Two SOC1 homologs were classified in subcluster 6 of the control group (Fig. 4a). KEGG enrichment indicated that photosynthesis (28.4%), photosynthesis - antenna proteins (38.6%), plant hormone signal transduction (10.1%), flavonoid biosynthesis (17.3%), and carbon fixation in photosynthetic organisms (13.2%) were over-represented in subcluster 6 of the control group (Fig. 4d). In the 6-BA treated group, two SOC1 transcripts exhibited the expression pattern represented by subcluster 3 (Fig. 4b, d). The top five enriched pathways in this subcluster included hormone signal transduction (8.8%), starch and sucrose metabolism (8.6%), mismatch repair (11.1%), galactose metabolism (11.0%), and arginine and proline metabolism (9.1%) (Fig. 4d). It is plausible that genes which were classified in the same subcluster as SOC1 may play an important role in floral transition.
Weighted gene co-expression network identifies biological processes and candidate genes associated with floral transition
Eight modules were identified by a weighted gene co-expression network analysis (WGCNA) with 8015 non-redundant DEGs (Fig. 5a). Module-trait relationship analysis revealed that the expression patterns of two of the flowering genes (SOC1 and SPL5), and two of the signaling components (ARR2 and TPS1), were highly related to the ‘turquoise’ modules. The ‘turquoise’ module genes were significantly enriched in hormone signal transduction pathway genes (Fig. 5c), containing a total of 84 genes. In addition, photosynthesis, fructose, and mannose metabolism and metabolic pathways were also over represented. Among the flowering-related genes, only SOC1 contained a lot of edges in the ‘turquoise’ module (Additional file 7: Table S5). The expression patterns of 80 of the 84 genes were highly connected with SOC1 (weight ≥ 0.1) (Fig. 5d). Key genes in sugar metabolism and signaling pathways, such as INV1, INVB, INVE, SPS2, SPS4, IDD4, IDD5, IDD7, IDD12, IDD14, TPS1, and TPS5 were also highly associated with SOC1 (Additional file 7: Table S5). The ‘blue’ module, however, was strongly positively correlated with SPL9 but negatively correlated with GAI; while the ‘black’ module was strongly negatively correlated with SPL9 but positively correlated with GAI (Fig. 5b). Sugar-related metabolic pathways were significantly enriched in the ‘blue’ module; as was the plant hormone signal transduction pathway (Fig. 5c). The Cytoscape representation of the 44 genes in the enriched sugar-related metabolic pathways and hormone signal transduction pathways in the ‘blue’ module had edges with SPL9 (weight ≥ 0.1) (Fig. 5e).
Hormone- and sugar-related genes during floral transtion
Based on their functional classification and the WGCNA analysis, genes whose expression was affected by the 6-BA treatment were mainly associated with hormone signal transduction pathways and sugar metabolism pathways. Particularly details of the changes in the expression of CK and GA signaling componets were presented in Fig. 6 and Additional file 8: Table S6. Three transcripts encoding homologs of DELLA proteins (GAI), one gene annotated as a GA receptor, (GID1C) and two GA-regulated transcription factors (PIF3) exhibited higher levels of expression in 6-BA treated buds at 30 and 70 DAFB, relative to untreated control buds; but lower expression at 50 DAFB (Fig. 6a). In contrast, one GA regulated transcription factor, GAMYB (GAM1), and a GID1 transcript (GID1B) exhibited the opposite response to the 6-BA treatment, relative to the response observed for the GA signaling components described above (Additional file 3: Table S2). The putative CK signaling components, cytokinin receptors (AHK4 and AHK2), three phosphotransmitters (AHP1), and four response regulators (ARR2, ARR9, and ARR6) were significantly upregulated in 6-BA treated buds at 30 and 70 DAFB but downregulated at 50 DAFB (Fig. 6c), relative to untreated controls. Two brassinosteroid signal transduction components (BSK and BZR1) and three cell cycle regulators (two CYCD3–1 and one CYCD3–2) were also present in the ‘turquoise’ module and exhibited a similar pattern of expression pattern as the cytokinin signaling components (Additional file 8: Table S6).
As presented in Fig. 6b, CK levels were not significantly different between the 6-BA treated and control buds at 30 DAFB but were significantly greater in 6-BA treated buds than in control buds at 50 DAFB. Not noly CK biosynthesis genes (IPT1) but also CK degradation gene (CKX1 and CKX7) exhibited higher expression levels in 6-BA treatment at 30 and 70 DAFB and lower transcription levels at 50 DAFB (Fig. 6c). The level of GAs was lower in 6-BA treated buds at 30 and 50 DAFB compared to control buds. At the same time, however, higher expression level of GA degradation genes (GA2OX2) at 30 and 50 DAFB and lower expression level of GA biosynthesis genes (GA20OX2, KAO, KO, KS, and CPS) at 50 DAFB were observed (Fig. 6b, Additional file 9:Table S7).
The levels of different sugars were also measured (Fig. 7b). Sucrose levels were higher in the 6-BA-treated samples than in the controls at 30 and 70 DAFB. The level of glucose in 6-BA treated buds was higher than in the control buds at 30 and 50 DAFB, while the concentration of fructose was significantly higher at 30 DAFB. Starch content was significantly lower in 6-BA treated buds than in the control buds at 30 DAFB, but similar in buds of both groups at all of the other sampling time points. Genes associated with the synthesis of fructose and glucose, including INV1, INV4, SPS2, SPS4, Glyco, and FKB, were significantly upregulated in the 6-BA treated buds at 30 and/or 50 DAFB (Fig. 7a). SSY2 genes, however, were downregulated in 6-BA treated buds at 30 and 70 DAFB (Fig. 7a). These genes may be partially responsible for the higher level of fructose and glucose, and lower level of starch, observed in the 6-BA treated buds. Two putative TPS transcripts (TPS1and TPS5), exhibited a higher level of expression in 6-BA treated buds than control buds at 30 DAFB. Another three TPP transcripts (TPPA, TPP4 and TPPD) that convert T6P into trehalose and Pi [33], were noticeably upregulated by the 6-BA treatment at 30 DAFB but downregulated at 50 DAFB (Fig. 7a). These transcript patterns of expression indicated that T6P metabolism and sucrose availability was more active in 6-BA treated buds during floral transition.
Homologs of flowering genes involved in six different flowering pathways were identified in the transcriptome data. All of them were found in the ‘turquoise’ module (Fig. 7c), except for COL3, SPL9, SPL12, and AGL9. All of SPLs genes (SBP2, SPL5, 6, 8, 14, and 15) were upregulated in 6-BA treated buds at 30 DAFB, except for SPL9 and SPL12, which were downregulated. COL, FD, and GI genes also exhibited different levels of expression in 6-BA treated buds vs. control buds. The expression of FT gene, however, which functions as an integrator in the photoperiod-regulated flowering pathway, did not exhibit differences in expression between the 6-BA and control buds (Additional file 9: Table S7). Genes that function as important mediators in the thermosensory, autonomous, and vernalization regulated flowering pathways, such as SVP, FCA, and VRN, were also detected in the transcriptome. Only minor changes were observed, however, in response to the 6-BA treatment and differences between treated and control buds did not reach a level of significance (Additional file 9: Table S7). Some flowering pathway integrator genes, such as SOC1 and AP1; as well as other members in MADS box family, were significantly upregulated in 6-BA treated buds at 30 and 70 DAFB.
Transcription factors that responded to the 6-BA treatment during floral transition
The expression of 573 non-redundant transcription factors was induced by the 6-BA treatment over the time course of floral transition (Fig. 8a). The identified transcription factors were classified into 56 distinct transcription factor families (Additional file 10: Table S8). Members of the SBP, bZIP, CO-like and MADS-box family of transcription factors have been reported to control flowering time (Yu et al., 2012, Shim et al., 2017). Most of these TFs exhibited a > 2-fold higher expression level in 6-BA treated buds than what was observed in control buds at 30 DAFB (Fig. 8c). Interestingly, members of the bHLH, NAC, HB and GATA transcription factor families were also found to be sharply upregulated in 6-BA treated buds at 30 DAFB (Fig. 8b). A total of 4, 6, 3, and 2 members of the bHLH, NAC, HB, and GATA transcription factor families, respectively, exhibited > 8-fold higher transcription levels in treated buds, relative to control buds, at 30 DAFB (Fig. 8c). It is plausible that these transcription factors may play an important role in directing the apple shoot apical meristem to flower development.