Morphology of CH and CL flowers in V. prionantha.
V. prionantha is a perennial plant that can develop both CH and CL flowers at different times during the flowering season. This is quite similar to observations in V. philippica [9]. CH flowers of V. prionantha had also five large and purple petals, with the lowest petal protruding slightly at the base into a spur, and five stamens forming a cone surrounding the pistil (Fig. 1A and B). Each stamen of V. prionantha had four pollen sacs, and the lowest two stamens were attached to noticeable nectar glands (Fig. 1C). The style of the pistil was erect and exceeded the stamens (Fig. 1D). In contrast, CL flowers of V. prionantha had two stamens without nectar glands, each stamen had two pollen sacs, and the five petals were all undeveloped as primordial structures (Fig. 1E–G). The style of the pistil was curved and bent to the two big stamens (Fig. 1H). Under certain conditions, intermediate CL (inCL) flowers developed in V. prionantha, and they displayed variable characteristics. Typical ones included one to three poorly developed petals, two to five developed stamens, and a curving style of the pistil (Fig. 1I–K). Therefore, the number and size of petals and stamens in CL and inCL flowers were smaller than those of normally developed CH flowers in V. prionantha. In addition, the filaments of CL and inCL stamens were distinct compared to that of CH flowers in V. prionantha, in which they were undetectable (Fig. 1C,G and J).
Photoperiod affects CH–CL flower development
Under natural conditions, V. prionantha developed complete CH flowers in the early spring and a mixture of inCL and CH flowers in the late spring, whereas complete CL flowers was produced in the summer and early autumn, suggesting that dimorphic flower development in this species might be regulated by the photoperiod. We thus used three photoperiods (10-, 12-, and 16-h daylight) to test this assumption in V. prionantha. The results demonstrated that complete CL flowers were formed under 16-h daylight, and both inCL and CH flowers were simultaneously developed at 10–12-h daylight in V. prionantha. Under a photoperiod of 10-h daylight, > 80% of the floral buds were developed to be CH flowers (Fig. 1L). The CH-CL formation pattern of V. prionantha in response to photoperiods is also similar to that in V. philippica [9]. In addition, flowering was induced as the photoperiod was extended in V. prionantha (Fig. 1M).
Upon naked-eye visualization of the floral buds, roughly five floral developmental stages (S1–S5) were defined in V. prionantha for both CH flowers under 10-h daylight (Fig. 1N–R) and CL flowers under 16-h daylight conditions (Fig. 1S–W). At S1, CH flowers had five obvious petals and stamens (Fig. 1N), whereas CL flowers had two stamens, and other stamens and all petals were retained as organ primordial structures (Fig. 1S). At S2, four whorl floral organs of CH flowers continued to develop, and the style of the pistil was higher than that of the stamens. The nectar glands at the base of the two stamens and spur of lower petals began to appear, while the style of the pistil in CL flowers began to bend to the two developed stamens (Fig. 1O and T). At S3, the nectar glands at the base of the two stamens were more obvious in the CH flowers, while the five petals and the other three stamens remained in the organ primordial state, and the filaments of stamen were more distinct in the CL flowers (Fig. 1P and U). At S4, the petals turned pink in the CH flowers, the filaments of the stamens became longer and finer, and the stamens were closer to the pistil in the CL flowers (Fig. 1Q and V). At S5, the petals turned purple, and the stamen cap became yellow in the CH flowers, while in the CL flowers, the stamen cap also became yellow, and the development of the three stamens and all petals was completely arrested as primordial structures (Fig. 1R and W).
To further reveal the effect of the photoperiod on CH and CL flower development, the photoperiod of V. prionantha plants was interconverted between 16- and 10-h daylight. In both cases of the “10-h to 16-h” daylight shift (Figure S1A–E) and “16-h to 10-h” daylight shift (Figure S1F–J), the CL or CH flower buds that had been induced in their original photoperiod continued to develop as their pre-programmed developmental trajectory in the new environment. Statistical analyses demonstrated that, within 8 days after the environment shift, the newly induced flower buds showed a series of inCL flower types, whereas buds initiated 8 days after the environment shift were directly transformed into CL flowers or CH flowers corresponding to the ultimate conditions of 16- or 10-h daylight (Figure S1K–P). These results indicate that the photoperiod regulates CH–CL flower development and transition in V. prionantha, and the long-day (LD) conditions can inhibit the development of stamens and petals. To understand the molecular basis for the dimorphic flower development or the CH–CL transition, we performed comparative transcriptomics between CL and CH flowers in V. prionantha.
Transcriptome sequencing and functional annotation
We first sequenced the full-length transcriptome of V. prionantha of whole plants, including the roots, stems, leaves, and involved floral buds, using the Pacific Biosciences Sequel platform, and generated 1,021,145 reads in total (Fig. 2A, Table S1). There were 155,918 full-length non-chimeric and clean reads, and 42,765 consensus isoforms were identified (Fig. 2A). The mean and N50 of the isoform length were 2346 bp and 3669 bp, respectively (Fig. 2A, Figure S2A). Functional annotations of the longest open reading frame (ORF), as the coding sequence (CDS), were performed using BLAST, Blast2GO, and InterProScan5, and 98.15% of transcripts were annotated (Fig. 2A, Tables S2 and S3). The NR database had the highest annotation rate, accounting for 96.1% of all transcripts, and the GO database had the lowest annotation rate, accounting for 35.39% of all transcripts (Table S3). The total number of transcripts annotated by all seven databases was 17,102, accounting for 40.74% of all transcripts (Figure S2B). With NR annotation, approximately 74% of all transcripts of V. prionantha were similar to Jatropha curcas, Populus trichocarpa, P. euphratica, and Ricinus communis (Figure S2C). To reveal the evolutionary relationship of V. prionantha, we first used Orthofinder to identify orthologous genes in the amino acid sequences of 17 other plant species, including P. deltoids, P. trichocarpa, and Arabidopsis thaliana, and a total of 24,585 gene orthogroups were constructed (Table S4). Phylogenetic analysis with these orthologous sequences using H. vulgare and O. sativa as an outgroup supported the notion that V. prionantha was the closest relative to P. deltoides and P. trichocarpa (Figure S3). We then evaluated the expansion and contraction of certain gene families and found that a total of 2555 families were expanded (73 significant, P ≤ 0.05), while 4586 gene families were contracted (120 significant, P ≤ 0.05) in V. prionantha (Figure S3). GO and KEGG enrichment analyses of the 1405 genes in the significantly expanded gene families (Table S5) revealed GO (Figure S4A, Table S6) and KEGG (Figure S4B, Table S6) categories/pathways, including ethylene-responsive proteins, lipoxygenase, and jasmonate ZIM domain protein, that might be involved in stress response and flower development.
RNA-seq between CH and CL flowers
To observe general variation patterns at the transcriptomic level between CH and CL flowers, three developmental stages (CH1-3 or CL1-3) of each flower type, corresponding to S1, S3, and S5, were subjected to RNA-seq analysis. We sequenced the 18 samples using the BGISEQ-500 platform and generated approximately 23.98 Mb RNA clean reads per sample on average (Table S7). By mapping them to the identified unique isoforms, 67.62% of the RNA clean reads were reasonably aligned (Table S8). Moreover, most transcripts were completely covered with uniform read coverage alone in the body using the RNA-seq data (Figure S5A and B). Gene expression correlation with samples showed that the samples in the three biological repeats at the same developmental stage had a high correlation, whereas the samples in different developmental stages had a relatively low correlation (Fig. 2B). Principal component analysis (PCA) also showed that three repeated samples at the same developmental stage were clustered together, whereas the samples at different stages were far away from each other (Fig. 2C). Thus, all generated RNA-seq data could be used for further analysis.
GO enrichment of stage-specific DEGs between CH and CL flowers
To identify genes related to dimorphic floral development in V. prionantha, we performed a DEG analysis by comparing the CH and CL flowers for the developmental stages (S1, S3, and S5) with DESeq2. A total of 162, 673, and 1716 genes were respectively found to be differentially expressed in S1, S3, and S5 (Table S9), and a Venn diagram revealed the distribution and relationship of these DEGs among these three developmental stages (Figure S6A). To further reveal the contribution of the DEGs’ variation to the CH-CL floral transition, we investigated the potential biological functions of the stage-specific DEGs through GO enrichment analysis.
In S1, eight enriched GO terms were constructed with 162 DEGs between the CH and CL flowers, and the most dominant was transcription regulation (Fig. 3A, Figure S6A, Table S10). Of these, 5 genes were downregulated and 14 genes were upregulated in CL flowers relative to CH flowers. The genes downregulated in CL flowers included four B-class MADS-box genes PISTILLATA (PI, Vp42765), APETALA3 (AP3, Vp31710, Vp32333, Vp33366), and one F-box transcription factor (TF) gene UNUSUAL FLORAL ORGANS (UFO, Vp21783) (Figure S6B), which might be critical for promoting the initiation and development of petals and stamens (Table S11). While the 14 upregulated TF genes in CL flowers included SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1, Vp24825), FLAVIN-BINDING KELCH REPEATF-BOX1/ ADAGIO3 (FKF1/ADO3) (Vp13028), ETHYLENE RESPONSE FACTOR53 (ERF053) (Vp15804, Vp19048), and ERF4 (Vp23464) (Figure S6B), which might be related to the circadian clock, carbohydrate synthesis, and defense response (Table S11).
A total of 23 enriched GO terms were identified for the 571 DEGs specific in S3 relative to S1, which mainly included defense response, carbohydrate metabolic process, nitrate assimilation, sporopollenin biosynthetic process, calcium ion transport, fatty-acyl-CoA biosynthetic process, glucose catabolic process, and jasmonic acid-mediated signaling pathways (Fig. 3B, Figure S6A, Table S10). Notably, 17 DEGs were involved in the defense response, and they were all upregulated in CL flowers compared with CH flowers (Fig. 3B). Of these, 16 DEGs involved in the carbohydrate metabolic process included 7 downregulated genes and 9 upregulated genes in CL flowers (Fig. 3B). In addition, genes in sporopollenin biosynthetic, fatty-acyl-CoA biosynthetic, and glucose catabolic processes were downregulated in CL flowers compared with CH flowers, and genes in nitrate assimilation and calcium ion transport were upregulated in CL flowers (Fig. 3B).
A total of 33 enriched GO terms were identified for the 1331 specific DEGs in S5 and included genes in the glycolytic process, carbohydrate metabolic process, and fatty acid biosynthetic process (Fig. 3C, Figure S6A, Table S10). Most of these genes were downregulated in CL flowers compared with CH flowers, whereas the genes involved in temperature compensation of the circadian clock were upregulated in CL flowers compared with CH flowers (Fig. 3C).
Taken together, the number of DEGs increased as the flower developed. The differentially expressed TF genes occurred at the earlier stages, and then the genes involved in metabolic processes related to sugar and fatty acids became differentially expressed between CH and CL flowers at late developmental stages, implying the essential role of these DEGs and their programmed regulation in the divergence of CH and CL flowers in V. prionantha.
Expression patterns of MADS-box genes in V. prionantha.
Downregulation of B-class MADS-box genes induces CL flowers in V. philippica [9]. The downregulated genes enriched in the regulation of transcription were mainly B-class MADS-box genes. We therefore characterized this gene family and obtained 23 MIKC-type MADS-box genes from the whole transcript data of V. prionantha. Phylogenetic analysis showed that the MADS-box genes of V. prionantha always clustered with their homologous genes in A. thaliana (Figure S7), suggesting their orthology. An expression profile of the identified MADS-box genes in three developmental stages of CH and CL flowers in V. prionantha showed that a high level of expression of the four B-class MADS-box genes VprAP3 (Vp31710, Vp32333, Vp33366) and VprPI (Vp42765) was observed in CH flowers relative to CL flowers (Fig. 4A, Table S12). In addition, the level of expression of the C-class MADS-box gene AGAMOUS (AG, Vp28027) was also higher in the CH flowers. However, the photoperiod-induced genes, such as SOC1 (Vp24825) and FRUITFULL (FUL, Vp29424), D-class MADS-box gene SEEDSTICK (STK, Vp30579), and AGAMOUS-LIKE6 (AGL6, Vp28883, Vp37510), showed lower expression levels in CH flowers than in CL flowers (Fig. 4A, Table S12).
We focused particularly on genes co-expressed with VprAP3/VprPI. A total of 37 genes showed the highest correlation (PCC > 0.9) with VprAP3/VprPI expression (Fig. 4B, Table S13). GO enrichment analysis showed that these genes were mainly involved in the fatty acid biosynthetic process, glycolytic process, carbohydrate metabolic process, glucose metabolic and catabolic process, positive regulation of flower development, and pollen development (Fig. 4C, Tables S14 and S15). Moreover, the expression of most genes co-expressed with VprAP3 and VprPI was higher in CH flowers of the three developmental stages than in CL flowers (Figure S8). The results suggested that glycolytic, fatty acid biosynthetic, and carbohydrate metabolic processes were putatively controlled by VprAP3/VprPI for dimorphic flower development. In addition, VprAP3/VprPI were also co-expressed with GAIP-B-like (Vp14987) and YABBY5 (YAB5, Vp31665) genes that regulate the development of petals, stamens, and nectaries (Fig. 4C, Table S15).
Expression profile of genes related to sugar–fatty acid processes in V. prionantha.
Genes involved in sucrose metabolic, glycolytic, and fatty acid biosynthetic processes were highly expressed in CH flowers compared to CL flowers, especially at S5 (Fig. 3C, Table S10). A higher expression level of sucrose metabolic process genes, i.e., SUCROSE-PHOSPHATE SYNTHASE 1 (SPS1, Vp9177, Vp9309) and SPS3 (Vp10032), which encode the key enzymes for sucrose synthesis, was observed in CH flowers than in CL flowers (Fig. 5A, Table S16). In the GO term of the glycolytic process, 16 genes were significantly differentially expressed, and most of them were upregulated in CH flowers at three developmental stages compared to CL flowers (Fig. 5B, Table S16).
Twenty-eight genes in the fatty acid biosynthetic process were significantly differentially expressed between the CH and CL flowers (Fig. 5C). The expression level of genes encoding biotin carboxyl carrier protein 2 (BCCP2, Vp25917), ketoacyl-CoA synthase (KCS, Vp15066, Vp19527, Vp30421, Vp19022, and Vp13523), and long-chain acyl-CoA synthetase (LACS, Vp13395, Vp14397, Vp15333, and Vp17649) participating in fatty acid biosynthesis was stably higher in CH flowers than in CL flowers during development, while LINOLEATE 9S-LIPOXYGENASE 5 (LOX1.5, Vp10416) was less expressed in CH flowers than in CL flowers (Fig. 5C, Table S16). Moreover, the expression of genes involved in lignin and cutin synthesis, such as COUMARATE—COA LIGASE (4CLL1, Vp19273, and Vp19820), was also constitutively higher in CH flowers than in CL flowers during development (Fig. 5C, Table S16). These results suggest that CH floral development might require more fatty acids, which could promote stamen fertility, increase the osmotic pressure of stamens and petals, and absorb and transport metal ions in floral organs (Table S16).
Additionally, there were 45 DEGs in the carbohydrate metabolic process, of which 14 genes were downregulated and 31 genes were upregulated in CH flowers compared to CL flowers (Fig. 5D). Most gene-related glycosidase families, including xylosidase (Xyl, Vp13859, Vp14913, Vp15375, Vp15991, and Vp16204), Glucan endo-1,3-beta-glucosidase (GLC1, Vp25276), and cellulose1 (CEL1, Vp40938), were stably downregulated in CH flowers compared to CL flowers, whereas polygalacturonase (PG) genes, including Vp34857, Vp39992, Vp24046, and Vp14690, were upregulated in CH flowers (Fig. 5D, Table S16).
Gene co-expression networks during CH–CL flower development
To facilitate the understanding of the regulatory network of flowering interconversion, we built weighted gene co-expression networks with all DEGs (2064 genes) among the three different developmental stages of V. prionantha. We chose a power of β = 10 based on the scale-free topology criterion to generate a hierarchical clustering tree (Figure S9A and B). A total of five co-expression modules (mergeCutHeight = 0.2) M1–M5 marked with different colors were ultimately identified, and they had a gene number in a module ranging from 71 to 969 (Fig. 6A). Correlations between the modules’ eigengenes (MEs) and genotypes/stages showed that all modules were obviously correlated with flower phenotypes, and M1 (383 genes) and M2 (478 genes) modules exhibited the most significant correlation (Fig. 6B). Significant correlations were also detected between stage-specific gene expression and module membership (MM) in each module (Figure S10), particularly in M1 and M2 modules (Figure S10A and B). Therefore, genes in M1 and M2 might play a role in CH–CL flower development. GO annotations revealed that most genes in M1 were involved in defense response, cellulose biosynthetic process, response to oxygen-containing compounds, response to karrikin, regulation of transcription, lipid metabolic process, and carbohydrate metabolic process (Fig. 6C, Table S17). The top 10 enriched GO terms in M2 were related to the response to jasmonic acid, response to fungus, response to blue light, regulation of transcription, protein dephosphorylation, positive regulation of LD photoperiodism, and flowering (Fig. 6D, Table S18).
To further identify the key genes essential for CH–CL flower development, we extracted hub genes in the M1 and M2 modules. The top 20 genes with the highest maximum clique centrality (MCC) values were considered hub genes (Fig. 6E and F, Tables S19 and S20). For instance, genes in M1, including HOMEOBOX-LEUCINE ZIPPER PROTEIN 12 (ATHB-12), PHYTOCHROME A SIGNAL TRANSDUCTION 1 (PAT1), MYB44, SENESCENCE-ASSOCIATED GENE 101 (SAG101), and GERMIN-LIKE PROTEIN 6 (GLP6), were responsible for abiotic and biotic stress in other plant species (Fig. 6E, Table S21), while genes in M2, including ERF4 and XYLOGLUCAN ENDOTRANSGLUCOSYLASE/HYDROLASE 23 (XTH23), were also related to stress responses (Fig. 6F, Table S21). Notably, a hub gene Vp8612 encoding V. prionantha GIGANTEA (VprGI) was found in M2, and this gene belonged to the GO terms of positive regulation of LD photoperiodism flowering (Fig. 6F, Table S21). These genes were expressed at a low level in CH flowers compared with CL flowers (Figure S11A and B). We then focused particularly on the VprGI co-expression network in V. prionantha. The expression of 17 genes showed the highest correlation (PCC > 0.8, MR ≤ 10) with VprGI expression (Fig. 6G, Table S22), and they included V. prionantha CONSTANS (VprCO), PSEUDO-TYPE RESPONSE REGULATORS1 (APRR1), LATE ELONGATED HYPOCOTYL (LHY), EXOCYST COMPLEX COMPONENT 70B1 (EXO70B1), ERF4, PBS1-like kinases 8 (PBL8), and PXY-correlated (PXC3), which are involved in the circadian clock and defense response progress (Table S23). The expression of these genes, except for LHY, was at low levels in CH flowers relative to CL flowers (Fig. 6H). In addition, PHOSPHOMETHYLPYRIMIDINE SYNTHASE (THIC), which plays an important role in linking the circadian clock and metabolic rhythm, was also at low levels in CH flowers relative to CL flowers in V. prionantha (Table S23).
The results suggest that LD light could induce the co-expression of floral inducers, such as VprGI and VprCO, and coincidently stress-response genes during the induction of CL flowers in V. prionantha.
Expression of key TF genes correlated with CH–CL development
Six key candidate regulatory genes (GI:Vp8612, CO:Vp14425, SOC1:Vp24825, UFO:Vp21783, AP3:Vp33366, and PI:Vp42765) in V. prionantha that might be essential for CH–CL floral transition were targeted and designated VprGI, VprCO, VprSOC1, VprUFO, VprAP3, and VprPI, respectively. Their expression variations during CH and CL development were validated using qRT-PCR. The expression of VprGI, VprCO, and VprSOC1 was at high levels in CL flowers under 16-h daylight relative to those subjected to 12- and 10-h daylight, whereas the expression of VprUFO, VprAP3, and VprPI was at high levels in CH flowers under 10-h daylight relative to inCL and CL flowers under 12- and 16-h daylight (Fig. 7A–C, Table S24). These results were consistent overall with the transcriptome data. We also investigated the expression of these genes in V. prionantha plants under 10-h daylight after they were switched to 16-h daylight conditions. The expression of VprGI, VprCO, and VprSOC1 in the first and second stages of floral buds were elevated with continuous increasing illumination and became significantly upregulated after a one-day shift from 10- to 16-h daylight, while expression variation of VprUFO, VprAP3, and VprPI showed the opposite trend (Fig. 7D, Table S24).