Morphological differentiation of shoot apexes during floral transition
The morphological differentiation of L. gratissima shoot apexes was observed through paraffin sections. The results showed that 0 d to 7 d under the SD treatment (SD0 to SD7) was the vegetative growth stage (undifferentiated stage), in which the tip of the growth cone in the bud was narrow and pointed and surrounded by leaf primordia (Fig. 2a-c). At 10 d after the initiation of the SD treatment (SD10), the bract primordial differentiation stage began in which the growth cone of the bud appeared hemispherical, which was an important marker of the transition from vegetative growth to reproductive growth. Subsequently, the two sides of the hemispherical growth cone elongated to form two protrusions, i.e., bract primordia (Fig. 2d-f). At 13 d after the initiation of the SD treatment (SD13), the inflorescence primordial differentiation stage began. At this stage, the growth cone in the bract primordia elongated to form three hemispherical protrusions, i.e., inflorescence primordia. Simultaneously, the lateral base of the bract primordia differentiated into lateral inflorescence primordia. Next, bilateral protrusions at each hemispherical inflorescence primordium differentiated into bract inflorescences (Fig. 2g-i). At 19 d after the initiation of the SD treatment (SD19), the floret primordial differentiation stage began and a single inflorescence primordium in the bract primordia gradually widened to become floret primordia at the tip of the bud (Fig. 2j-l). These results showed that the floral transition period began 10 d after the initiation of the SD treatment, and the selection of time points before and after this period could facilitate the physiological study of floral transition. Therefore, the LD treatment was used as a control in this study, and 7 d (SD7), 10 d (SD10), 13 d (SD13), and 19 d (SD19) in the SD treatment were selected to study the physiological and molecular regulation patterns of floral transition.
Dynamic changes in endogenous substance content during floral transition
Contents of soluble sugars and endogenous hormones (gibberellin [GA3], indoleacetic acid [IAA], abscisic acid [ABA], and zeatin [ZT]) were measured at 0 d [SD0], 7 d [SD7], 10 d [SD10], 13 d [SD13], and 19 d [SD19] after the initiation of the SD treatment. The Kruskal-Wallis H test results showed that except for GA3, which could not be detected because it was below the limit of quantitation (0.1 ng/ml), there were significant differences in the contents of the other substances among the five stages (adjusted p < 0.05; Fig. 3). Soluble sugar, ZT, and IAA reached their peaks at SD0, which were 28.86 ± 0.67 mg g-1 FW, 2.15 ± 0.30 ng g-1 FW, and 0.69 ± 0.04 ng g-1 FW, respectively. Additionally, soluble sugar and ZT decreased from SD0 to SD19. Interestingly, IAA showed an increase in SD13 before decreasing. Similarly, ABA initially increased from SD0 to SD13 and subsequently declined.
Additionally, post-hoc results showed that there were extremely significant differences in the pairwise comparisons between the five timepoints for ABA (p < 0.001). IAA only showed no significant differences between SD7 and SD13 (p > 0.1). Soluble sugar did not show any significant differences between SD0 and SD7 (p > 0.1). ZT did not show any significant differences between SD0 and SD10, between SD7 and SD10, or between SD13 and SD19 (p > 0.1). From these results, it can be seen that ABA levels changed rapidly and dynamically over these five stages, whereas ZT levels exhibited little change over this same time period. Changes in soluble sugar content mainly occurred in later periods (SD13 to SD19). In contrast to these substances, IAA changes were relatively constant (Fig. 3).
RNA-seq and qRT-PCR identification of differentially expressed genes (DEGs)
Transcriptomes were generated for three biological replicates from the SD and LD treatments at each of the time points corresponding to the four stages of bud differentiation in SD treatment plants, i.e., at 7 d, 10 d, 13 d, and 19 d (Fig. 2), yielding a total of 24 transcriptomes. A total of 1,236,426,670 raw sequencing reads were generated from 24 samples, 1.2 × 109 high-quality clean reads (181 Gb) were obtained after filtering, with mean Q20, Q30, and GC contents of 99.11%, 97.18%, and 43.53%, respectively (Additional file 1: Table S1). A total of 79,870 unigenes (≥ 200 b) were generated from de novo assembly, and the N50 length was 2,118 bp (Additional file 2: Table S2). Among these unigenes, 35,725 unigenes (44.73%) were successfully annotated to at least one database (Additional file 3: Fig. S1).
With RNA from the same 24 samples used for transcriptome generation, qRT-PCR was conducted for nine flowering-related unigenes identified in through RNA-seq, including COP1 (Unigene0031506), ZTL (Unigene0041339), FKF1 (Unigene0038380), GI (Unigene0051409), ELF3 (Unigene0051761), PRR1 (Unigene0045946), PRR7 (Unigene0003564), PRR5 (Unigene0047475), and LHY (Unigene0035686). The results of qRT-PCR showed that except for PRR5 (Unigene0047475), the expression patterns of the other eight genes were generally consistent with the RNA-seq data (Additional file 4: Fig. S2), indicating that the transcriptome data generated in this study were reliable and valid.
A total of 113 (LD7-vs.-SD7), 420 (LD10-vs.-SD10), 483 (LD13-vs.-SD13), and 464 (LD19-vs.-SD19) differentially expressed genes (DEGs) were obtained by comparing the LD and SD treatments (Additional file 5: Fig. S3; Additional file 6: Table S3). A total of 1,226 DEGs were identified from these four comparisons, of which five DEGs were shared by four comparisons, and 250 DEGs were present in at least one comparison. There were 110, 302, 288, and 276 stage-specific DEGs in LD7-vs.-SD7, LD10-vs.-SD10, LD13-vs.-SD13, and LD19-vs.-SD19, respectively (Additional file 5: Fig. S3).
Functional classifications of DEGs
MapMan is an effective tool for systematic analysis of plant transcriptome metabolic pathways and other biological processes [12]. We employed MapMan to overview transcriptional changes in regulatory, metabolic, and cellular response-related genes (Additional file 7: Table S4). In “regulation overview”, more DEGs were detected in the other three comparisons contrasted with LD7-vs.-SD7, showing that the physiological and molecular characteristics after flower bud differentiation (SD10, SD13, and SD19) were significantly different from that before flower bud differentiation (SD7). In the IAA metabolic subclass, more DEGs were upregulated in LD19-vs.-SD19 compared with LD7-vs.-SD7, LD10-vs.-SD10, and LD13-vs.-SD13 (Additional file 8: Fig. S4). However, IAA content was the highest at SD0 and lowest on SD10 and continuously decreased from SD13 to SD19 (Fig. 3). Therefore, IAA was not a key factor mediating floral transition in L. gratissima. ABA metabolism-related DEGs were significantly upregulated in all four comparisons (Additional file 8: Fig. S4), and ABA levels were overall increasing in the process of floral transition (Fig. 3), demonstrating that ABA could promote floral transition in L. gratissima. In “minor CHO metabolism”, trehalose biosynthesis-related DEGs were only upregulated in LD7-vs.-SD7 (Additional file 9: Fig. S5). “Cellular response overview” showed that more development-related DEGs were upregulated in LD10-vs.-SD10 compared with the other three comparisons (Additional file 10: Fig. S6), indicating that these DEGs promoted floral transition in L. gratissima.
Co-expression module analysis for DEGs
Weighted gene co-expression network analysis (WGCNA) is a systems biology method for analyzing the correlation relationships between genes in multiple samples [13]. In this study, the results of WGCNA showed that 1,226 DEGs in eight samples were clustered in 11 different co-expression modules (labelled with different colors) (Fig. 4a). It is noteworthy that four out of 11 co-expression modules significantly correlated with a single sample (r > 0.9, p < 0.05; Fig. 4b; Additional file 11: Table S5). For example, the largest module (black module) included 247 (20.15%) SD19-specific DEGs (Fig. 4b; Additional file 11: Table S5a).
We further conducted GO enrichment analysis on 11 co-expression modules, and only the greenyellow module was not significantly enriched for any GO terms (Additional file 12: Table S6). Some GO terms were specifically identified in only a single module. For example, 120 specific GO terms were identified in the black module, which mainly involved signal transduction and negative regulation of metabolic processes, and 34 module-specific GO terms were identified in the brown module, which was mainly associated with growth and development (Additional file 12: Table S6). However, several GO terms, including “response to organic substance” and “response to a stimulus”, appeared in multiple modules (Additional file 12: Table S6), indicating possible module-gene interactions. Overall, the extensively enriched GO terms showed that floral transition in L. gratissima activated functional activities at the systems level.
The 11 modules were divided into seven categories based on the correlations between modules (Fig. 4c). The heatmap showed that there was a high correlation between the blue, magenta, pink, and tan modules, in which the genes were highly expressed in SD7 and SD10 (Fig. 4b, c), and were significantly enriched in multiple GO terms involving secondary metabolite biosynthesis, signal transduction, and regulation of developmental processes (Additional file 12: Table S6).
Identification of DEG expression patterns associated withfloral transition in L. gratissima
According to the above functional classifications and WGCNA of these DEGs, and flowering-related genes previously reported in model plants (such as A. thaliana) [5, 7], a total of 146 unigenes were identified as homologous genes related to floral transition in L. gratissima, involving several flowering pathways: sugar metabolism, hormone metabolism and signal transduction, photoperiod, ambient temperature, aging pathways, as well as floral integrator and floral meristem identity genes. Among these floral transition-related homologous genes, stage-specific DEGs, and common DEGs in LD7-vs.-SD7, LD10-vs.-SD10, LD13-vs.-SD13, and LD19-vs.-SD19 are listed in Additional file 13: Table S7.
The expression pattern of sugar signal-related homologs
The sugar signal pathway, which responds to the sugar budget in plants, is one of the important pathways mediating the transition from vegetative to floral meristems [5]. A total of 29 (19.86%) DEGs associated with sugar signal-related genes were identified, involving 23 sugar signal-related homologs. These genes expressed differently in different development stages of L. gratissima. For example, HEXOKINASE (HK) homologs (Unigene0044869 and Unigene0044870) were significantly upregulated in LD7-vs.-SD7 and LD13-vs.-SD13, and a BETA-GLUCOSIDASE 24 (BGLU24) homolog (Unigene0013088) was significantly upregulated in LD10-vs.-SD10. Meanwhile, Unigene0009721 and Unigene0041893, homologs of GALACTINOL SYNTHASE 2 (GOLS2) and RAFFINOSE SYNTHASE (RFS) participating in raffinose synthesis, were upregulated in LD7-vs.-SD7. In addition, TREHALOSE-6-PHOSPHATE SYNTHSE (TPS) homologs (Unigene0019787, Unigene0024389, Unigene0013555, Unigene0054604, Unigene0004913, and Unigene0062998) were upregulated at various stages, and SWEET16 homolog (Unigene0012661) was significantly upregulated in LD7-vs.-SD7 and LD10-vs.-SD10 (Fig. 5e; Additional file 14: Table S8). Hence, these genes may directly or indirectly participate in floral transition in L. gratissima.
The expression patterns of phytohormone metabolism and signal transduction homologs
Many studies have demonstrated that various phytohormones participate in the regulation of floral transition [7, 14-16]. A total of 20 (13.70%) DEGs associated with phytohormone metabolism were identified, and these involved 16 phytohormone metabolism homologous genes and were related to nine phytohormone metabolism pathways. Among these genes, GIBBERELLIN 2-BETA-DIOXYGENASE 1 (GA2OX1) homologs (Unigene0030732) and GA2OX8 homologs (Unigene0073113), which are involved in GA metabolism, were significantly upregulated in LD10-vs.-SD10 and/or LD19-vs.-SD19. Meanwhile, Unigene0034382 (CYP707A1 homolog) and Unigene0042754 and Unigene0042755 (NCED1 homologs), respectively, encoding abscisic acid (ABA) 8'-hydroxylase 1 and nine-cis-epoxycarotenoid dioxygenase, were significantly upregulated in LD10-vs.-SD10. In addition, a homolog (Unigene0035296) of YUC4, encoding indole-3-pyruvate monooxygenase, which mediates auxin biosynthesis, was significantly upregulated in LD19-vs.-SD19. Additionally, genes encoding cytokinin (CK) dehydrogenase 7 (CKX7; Unigene0036599) and cytokinin dehydrogenase (CYP735A1; Unigene0029738) were significantly downregulated in LD19-vs.-SD19. CYTOCHROME P450 734A1 (CYP734A1) homolog (Unigene0036368), which participates in brassinolide (BR) biosynthesis, was upregulated in LD10-vs.-SD10 and LD19-vs.-SD19; the jasmonate (JA) metabolism-related JASMONATE O-METHYLTRANSFERASE (JMT) homolog (Unigene0020912) and the salicylic acid (SA) metabolism-related UDP-GLYCOSYLTRANSFERASE 74F1 (UGT74F1) homolog (Unigene0004033) were downregulated in LD19-vs.-SD19. A homolog of CAROTENOID CLEAVAGE DIOXYGENASE 7 (CCD7, Unigene0069349) involving in strigolactone (SL) biosynthesis, was also identified and showed significant downregulation in LD10-vs.-SD10 (Fig. 5c; Additional file 14: Table S8).
A total of 39 (25.85%) DEGs associated with phytohormone signal transduction were identified and involved 30 phytohormone signal transduction homologs that were associated with signal transduction for nine hormones. Among these DEGs, GID1B homologs (Unigene0032780, Unigene0032781, and Unigene0063035), encoding a gibberellin receptor, were upregulated in LD10-vs.-SD10, whereas an RGL3 homolog (Unigene0071862), encoding a DELLA protein, was significantly downregulated in LD19-vs.-SD19. The ABA signal transduction-related EID1-LIKE F-BOX PROTEIN 3 (EDL3) homolog (Unigene0018152) was upregulated in LD10-vs.-SD10, and SAUR71 homologs (Unigene0021953 and Unigene0025106), encoding auxin-responsive protein, were upregulated in LD13-vs.-SD13 and LD19-vs.-SD19. Moreover, in the CK signaling pathway, homologs of AHPs (Unigene0034629, Unigene0004315, and Unigene0034630), encoding histidine-containing phosphotransfer protein, were highly expressed in SD10, SD13, and SD19, and an ARR6 homolog (Unigene0049441), encoding a two-component response regulator, was upregulated in LD19-vs.-SD19. In addition, a BRI1 homolog (Unigene0024976) in the BR signaling pathway was significantly upregulated in LD7-vs.-SD7; homologs of MYC4 (Unigene0009399) and TIFYs (Unigene0022959 and Unigene0019294) in the JA signaling pathway were upregulated in LD10-vs.-SD10; a DWARF14 (D14) homolog (Unigene0028658), participating in SL signal transduction, was upregulated in LD7-vs.-SD7 but downregulated in LD19-vs.-SD19 (Fig. 5d; Additional file 14: Table S8).
Expression patterns of genes associated with photoperiod pathways
The photoperiod flowering pathways in plants include the photosensory pathway, the circadian clock, and the systemic effector [17]. A total of 10 (6.84%) photoperiod-related homologs were identified. Among these homologs, CHLOROPHYLL A-B BINDING PROTEIN (CAB40, Unigene0075619) was downregulated in LD19-vs.-SD19, whereas CONSTANS-LIKE 12 (COL12, Unigene0039617) and FD (Unigene0027311) were upregulated in LD13-vs.-SD13. Meanwhile, homologs of the FLAVIN-BINDING KELCH REPEAT F-BOX PROTEIN 1 (FKF1, Unigene0038380) and the PSEUDO-RESPONSE REGULATOR 7 (PRR7, Unigene0003564) were both downregulated in LD13-vs.-SD13, and a LUX homolog (Unigene0011585) was downregulated in LD10-vs.-SD10, whereas homologs of the NF-Ys (Unigene0025001, Unigene0002375, and Unigene0033157) were upregulated at one or more stages (Fig. 5a; Additional file 14: Table S8).
Expression patterns of genes associated with the ambient temperature pathway
Plant responses to photoperiod and temperature are coupled [18, 19]. The photoperiod-induced floral transition could also affect the expression of a series of ambient temperature-related genes in plants. We identified 28 (19.18%) ambient temperature-related DEGs involving 18 homologs, primarily including the HEAT STRESS TRANSCRIPTION FACTORS (HSFs), HEAT SHOCK PROTEIN/ COGNATE (HSPs), and pEARLI1, most of which were highly expressed at several stages under LD (Fig. 5b; Additional file 14: Table S8).
Expression patterns of aging pathway-related, floral integrator, and floral meristem identity genes
The aging pathway is an endogenous flowering pathway in plants [20]. SQUAMOSA PROMOTER-BINDING-LIKE PROTEIN 4 (SPL4) homologs (Unigene0024429 and Unigene0024430) in the aging pathway were upregulated in LD10-vs.-SD10, LD13-vs.-SD13, and LD19-vs.-SD19 (Fig. 5f; Additional file 14: Table S8).
Floral integrators combine environmental and endogenous signals to mediate flowering in plants [5]. The floral integrator gene SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1) homologs (Unigene0039572 and Unigene0039575) were upregulated in LD10-vs.-SD10, LD13-vs.-SD13, and LD19-vs.-SD19, whereas the AGAMOUSLIKE24 (AGL24) homolog (Unigene0049016) was downregulated in LD19-vs.-SD19 (Fig. 5g; Additional file 14: Table S8).
Genetic networks regulating floral transition in plants ultimately activated floral meristem identity genes, thereby causing the transformation from vegetative to floral meristems [21]. A total of 15 (10.27%) related DEGs were identified, involving nine floral meristem identity genes (Additional file 14: Table S8). Among these genes, homologs of AGL8/FUL (AGL8, also known as FUL; Unigene0019277, Unigene0004737, Unigene0042052, Unigene0042053, and Unigene0042058), APETALA 1 (AP1; Unigene0019278, Unigene0019279, and Unigene0031106), LEAFY (LFY; Unigene0030979 and Unigene0030980), and SEPALLATAs (SEPs; Unigene0000607, Unigene0034045, and Unigene0025130) were upregulated in one or more developmental stages, whereas homologs of SHORT VEGETATIVE PHASE (SVP, Unigene0049018) and TERMINAL FLOWER 1 (TFL1, Unigene0026727) were downregulated in LD19-vs.-SD19, LD10-vs.-SD10, and LD13-vs.-SD13 (Fig. 5h; Additional file 14: Table S8).
Co-expression network of floral transition-related genes
A co-expression network constructed using 126 floral transition-related DEGs with edge weights > 0.1 showed ten hub genes with great connectivity, including homologs of GLYCERALDEHYDE-3-PHOSPHATEDEHYDROGENASE (GAPDH, Unigene0005846), AKR1B1 (Unigene0076531), PKM (Unigene0073914), ENOLASE1 (ENO1; Unigene0011083 and Unigene0011084), MED37E (Unigene0051600), L-LACTATE DEHYDROGENASE A CHAIN (LDHA, Unigene0009368), HSP83A (Unigene0031524), FUL (Unigene0042052), and SEP4 (Unigene0025130) (Additional file 15: Fig. S7). The genes with the highest network degree were GAPDH (Unigene0005846), AKR1B1 (Unigene0076531), and PKM (Unigene0073914), which participated in sucrose and starch catabolism (Additional file 14: Table S8).