Eight libraries were sequenced using the Illumina Hiseq 2000 platform (Table 1, Fig. 1). After filtering the raw data, a total of 535,683,215 clean reads were obtained. Of these reads, 130,705,691 reads with an average GC content of 50.08% were from MA-T1 and MA-T2, 126,431,320 reads with an average GC content of 50.03% were from AE-T1 and AE-T2, 138,068,167 reads with an average GC content of 50.64% were from MA-T3 and MA-T4, and 140,478,037 reads with an average GC content of 49.61% were from AE-T3 and AE-T4 (Table S1). Q30 percentage (percentage of bases with sequencing error rate lower than 1‰) were 93.88% for MA-T1 and MA-T2, 93.94% for AE-T1 and AE-T2, 95.54% for MA-T3 and MA-T4, 92.22% for AE-T3 and AE-T4. These data showed that the RNA-Seq quality was applicable for further analysis.
After the assembly of the clean reads, a total of 57949 unigenes were generated from T1 and T2, which had an average length of 1276 bp and N50 length of 2023 bp (Table S2). A total of 53,699 unigenes were obtained from T3 and T4, which had an average length of 1419 bp and N50 length of 2329 bp (Table S2). In both unigene libraries, unigenes between 300bp and 500bp had the largest percentage (Fig. S1 and Fig. S2). These data indicated the high quality of the assembled unigenes.
A total of 33812 unigenes (58.35%) from T1 and T2, 32469 unigenes (60.46%) from T3 and T4 were annotated to Nr, Swissprot, GO, KEGG, COG and Pfam database. Among the annotated unigenes from T1 and T2, 33243 unigenes (98.32%) were annotated to Nr, 23920 (70.74%) to Swissprot, 19239 (56.90%) to GO, 13752 (40.67%) to KEGG, 10885 (32.19%) to COG and 23397 (69.20%) to Pfam database (Table S3). Among the annotated unigenes from T3 and T4, 32072 unigenes (98.78%) were annotated to Nr, 22773 (70.14%) to Swiss-prot, 18525 (57.05%) to GO, 12913 (39.77%) to KEGG, 10603 (32.66%) to COG and 22845 (70.36%) to Pfam database (Table S3).
DEGs related flower opening were screened from comparisons of MA-T1 vs. MA-T2 and AE-T1 vs. AE-T2 with FDR (False Discovery Rate) <0.1 and FC (Fold Change) ≥2. DEGs related flower closing were filtered from comparisons of MA-T3 vs. MA-T4 and AE-T3 vs. AE-T4 with FDR <0.1 and FC ≥1.5. In the comparison of MA-T1 vs. MA-T2, a total of 5832 unigenes were differentially expressed (2986 up-regulated and 2846 down-regulated), in which 4581 unigenes were annotated against the six public databases (Nr, Swiss-prot, GO, KEGG, COG and Pfam) (Fig. 2; Table S4 and Table S5). In the comparison of AE-T1 vs. AE-T2, a total of 4711 DEGs (3718 annotated) were identified, in which 2910 were up-regulated and 1801 were down-regulated (Fig. 2, Table S4 and Table S5). In the comparison of MA-T3 vs. MA-T4, a total of 975 unigenes (828 annotated) were differentially expressed (493 up-regulated and 482 down-regulated) (Fig. 2, Table S4 and Table S5). In the comparison of AE-T3 vs. AE-T4, a total of 2272 DEGs (1948 annotated) were screened, in which 831 were up-regulated and 1441 were down-regulated (Fig. 2, Table S4 and Table S5).
Table 1
The sampled time and sampled blooming types for transcriptome sequencing
The time points within a day
|
The time of MA starting opening (T1, 9:30)
|
The time of AE starting opening (T2, 14:30)
|
The time of MA starting closing (T3, 18:15)
|
The time of AE starting closing (T4, 20:00)
|
morning-opening and afternoon-closing plant (MA)
|
MA-T1
|
MA-T2
|
MA-T3
|
MA-T4
|
afternoon-opening and evening-closing plant (AE)
|
AE-T 1
|
AE-T2
|
AE-T3
|
AE-T4
|
KEGG pathway and GO enrichment analysis of DEGs
KEGG pathway and GO enrichment analysis of unigenes differentially expressed at different flowering stages provided an overview into the gene function or pathways associated with flower opening and closing. For flower opening process, “amino acid transmembrane transport”, “water transport”, “circadian rhythm” and “fructose and mannose metabolism” were significantly enriched in DEGs of the both comparisons of MA-T1 vs. MA-T2 (SO stage vs. FO stage) and AE-T1 vs. AE-T2 (B stage vs. SO stage) (Fig. 3). “Cell wall biogenesis”, “xyloglucan metabolic process”, “xyloglucan: xyloglucosyl transferase activity”, “plant hormone signal transduction” and “starch and sucrose metabolism” were specifically enriched in up-regulated gene sets of AE-T1 vs. AE-T2 (B stage vs. SO stage) (Fig. 3). In the comparison of MA-T1 vs. MA-T2 (SO stage vs. FO stage), “autophagy” and “acid phosphatase activity” were specifically enriched in up-regulated gene sets (Fig. 3).
As for flower closing process, “cell wall organization”, “circadian rhythm” and “starch and metabolism” were dramatically over-represented in DEGs of both comparisons of MA-T3 vs. MA-T4 (SC stage vs. FC stage) and AE-T3 vs. AE-T4 (FO stage vs. SC stage) (Fig. 4). Of these terms or pathways, cell wall organization term was over-represented in up-regulated gene sets of MA-T3 vs. MA-T4 (SC stage vs. FC stage), and in down-regulated gene sets of AE-T3 vs. AE-T4 (FO stage vs. SC stage) (Fig. 4). Circadian rhythm pathway was enriched in down-regulated gene sets of the both comparisons (MA-T3 vs. MA-T4 and AE-T3 vs. AE-T4). Terms (pathways) specifically enriched in up-regulated gene sets of AE-T3 vs. AE-T4 included carbon-carbon lyase activity, hydrolase activity (acting on glycosyl bonds) and plant hormone signal transduction (Fig. 4). Terms (pathways) specifically enriched in down-regulated gene sets of AE-T3 vs. AE-T4 included glucan metabolic process, pectin catabolic process and pectinesterase activity (Fig. 4). Protein processing in endoplasmic reticulum and ribosome pathway were specifically enriched in down-regulated gene sets and up-regulated gene sets of MA-T3 vs. MA-T4, respectively (Fig. 4).
Transcriptional regulatory mechanism of differential flower opening and closing times between different plants
The purpose of this study was to use RNA-seq on the same time-course transcriptomes of two plants with different within-day flowering course to identify key candidate genes regulating flower opening and closing times. Genes that show different change trends of transcript abundance between MA and AE on the same time-course but different flowering course are valuable for understanding the regulatory mechanisms of flower opening and closing times. The above enrichment analysis and previous reports suggest that flower opening and closing (senescing) would be involved in light input, circadian clock, hormonal regulation, water transport, osmoregulation, cell wall organization and developmental programmed cell death (dPCD). The expression levels and transcript change trends in abundance of these related genes were displayed with heat maps and the detailed data was shown in Table S8.
Light input and circadian clock
In this study, flowers of a single hybrid plant (or a species) start opening and closing at exactly the same time of different days, even in constant darkness (unpublished data). These facts suggested that flower opening and closing times of I. dichotoma, I. domestica and their hybrids were controlled by circadian clock. Furthermore, the flower opening and closing times of I. dichotoma and I. domestica changed with different light-dark cycles (unpublished data), suggesting the involvement of light input in the regulation of flower opening and closing times. The light- and clock- regulation of flower opening and closing times also have been reported in many other species, such as Eustoma grandiflorum and Portulaca umbraticola [8, 9].
In flower opening process (from T1 to T2), two genes associated with circadian clock and light signaling showed different expression profiles between MA and AE. The putative core clock gene EARLY FLOWERING 4-LIKE 4 (ELF4-LIKE 4) was evidently down-regulated in the comparison of MA-T1 vs. MA-T2 but up-regulated in the comparison of AE-T1 vs. AE-T2 (Fig. 5B). Another gene encodes PHYTOCHROME-INTERACTING FACTOR4-like (PIF4-like), whose homologous gene in Arabidopsis thaliana functions as a signaling hub downstream of external signals (light and temperature) and intrinsic signals from the circadian clock [25]. Its transcript abundance significantly increased in MA-T1 vs. MA-T2 (SO stage vs. FO stage) but declined in AE-T1 vs. AE-T2 (B stage vs. SO stage) (Fig. 5A). In flower closing process (from T3 to T4), PIF4-like also presented a different change trend of transcript abundance between MA and AE. A transcript encoding PIF4-like was normally regulated in MA-T3 vs. MA-T4 (SC stage vs. FC stage) but up-regulated in AE-T3 vs. AE-T4 (FO stage vs. SC stage) (Fig. 5C).
Hormones
Plant hormones can orchestrate multiple cellular processes that impact physiology and growth. The regulation of hormones on flower opening has been reported in Iris species [26], as well as flower closing caused by senescence [27].
In flower opening process (from T1 to T2), six gene transcripts associated with auxin synthesis and signaling were detected to present opposite change trends of abundance between MA and AE (Fig. 6A). Two gene transcripts (YUCCA2-like and YUCCA10) involved in auxin biosynthesis [28] were down-regulated in MA-T1 vs. MA-T2 (SO stage vs. FO stage) but up-regulated in AE-T1 vs. AE-T2 (B stage vs. SO stage) (Fig. 6A). A transcript encoding auxin-responsive protein SMUX AUXIN UP RNA 64-like (SAUR64-like), and a transcript encoding auxin-induced protein 22D were also detected to present up-regulation in MA-T1 vs. MA-T2 but down-regulation in AE-T1 vs. AE-T2 (Fig. 6A). Five transcripts encoding auxin efflux carrier component 1 and auxin efflux carrier component 7-like did not show significant differential expression in MA-T1 vs. MA-T2, but were up-regulated in AE-T1 vs. AE-T2 (Fig. 6A). However, the putative auxin efflux carrier component 5 and PIN-LIKES 7-like were down-regulated in the both comparisons of MA-T1 vs. MA-T2 and AE-T1 vs. AE-T2 (Fig. 6A). In Arabidopsis thaliana, the PINFORMED (PIN) efflux carriers function as auxin transporters. Among the PIN family proteins, PIN1 and PIN7 are localized at the plasma membrane and function as auxin efflux carriers [29, 30]. However, PIN5 and PIN-LIKES (PILS) proteins are localized in the endoplasmic reticulum (ER) and facilitate intracellular auxin transport between the cytosol and ER to regulate auxin homeostasis [29, 31]. The putative auxin efflux carrier components and PILS detected in this study belong to PIN family and the function of their proteins may correspond to the PIN proteins of A. thaliana.
Five transcripts involved in gibberellin synthesis and signaling were detected to perform different change trends of abundance between MA and TE in flower opening process (from T1 and T2) (Fig. 6A). Among them, two transcripts encoding ent-kaurenoic acid oxidase (KAO), which are involved in biosynthesis of gibberellin [32] were up-regulated in AE-T1 vs. AE-T2 but normally regulated in MA-T1 vs. MA-T2 (Fig. 6A). The other two transcripts encoding the DELLA protein GAI-like were down-regulated in AE-T1 vs. AE-T2 but did not show significant expressional changes in MA-T1 vs. MA-T2 (Fig. 6A). DELLA proteins function as negative regulators in gibberellin signaling [33].
Only one gene transcript related to cytokinin was detected to have different change trends of abundance between MA and AE in flower opening process (from T1 to T2). The transcript encodes cytokinin dehydrogenase 5-like, a protein associated with degradation of cytokinins [34]. It was up-regulated in MA-T1 vs. MA-T2 but normally regulated in AE-T1 and AE-T2 (Fig. 6A).
Jasmonate (JA)-related gene transcripts almost showed the same change trends of abundance between MA and AE in flower opening processes (from T1 to T2). A gene transcript encoding an enzyme involved in biosynthesis of jasmonates (allene oxide cyclase, AOC) [35], and a transcript encoding ninja-family protein involved in jasmonate signaling [35], were down-regulated in the both comparisons of MA-T1 vs. MA-T2 and AE-T1 vs. AE-T2 (Fig. 6A). Two gene transcripts encoding cytochrome P450 enzyme CYP94B3 were up-regulated in the both comparisons of MA-T1 vs. MA-T2 and AE-T1 vs. AE-T2 (Fig. 6A). The CYP94B3 is reported to cause the hydroxylation of JA-Ile which may switch off JA/JA-Ile signaling [35, 36].
In flower closing (senescing) process (from T3 to T4), only three auxin-related gene transcripts showed different change trends of abundance between MA and AE. Two gene transcripts encoding auxin-responsive proteins (SAUR64-like and SAUR68-like) were normally regulated in MA-T3 vs. MA-T4, but were down-regulated in AE-T3 vs. AE-T4 (Fig. 6B). A gene transcript encoding auxin efflux carrier component 5-like (functioning as intracellular auxin transporter) was down-regulated in MA-T3 vs. MA-T4 but up-regulated in AE-T3 vs. AE-T4 (Fig. 6B). Almost all gene transcripts associated with gibberellin and cytokinin presented the same change trends of abundance between MA and AE in flower senescing process and there were no cytokinin-related DEGs in the two comparisons of MA-T3 vs. MA-T4 and AE-T3 vs. AE-T4 (Fig. 6B). Two gene transcripts encoding CYP94B3 involved in jasmonate metabolism were down-regulated in AE-T3 vs. AE-T4, but did not show significant transcript changes in MA-T3 vs. MA-T4 (Fig. 6B).
Water transport and osmoregulation
Flower opening is generally due to cell expansion driven by intracellular turgor pressure [1, 18]. The increase of intracellular turgor pressure is caused by water flowing into the cell and most water transports across the plasma membrane are mediated by aquaporins [18]. Aquaporins are mostly localized to the plasma membrane or vacuolar membrane, which are respectively called plasma membrane intrinsic proteins (PIPs) and tonoplast intrinsic proteins (TIPs) [37]. In flower opening process (from T1 to T2), two gene transcripts encoding PIP2-5 and PIP2-2 were detected to perform down-regulation in MA-T1 vs. MA-T2 but up-regulation in AE-T1 vs. AE-T2 (Fig. 7A). The gene transcript encoding TIP2-3 showed the same change trend of abundance in the two comparisons (Fig. 7A).
Water flows through aquaporins along the water potential gradient. Hence, water transports can be indirectly regulated by changing osmotic pressure and the change of osmotic pressure can be achieved by changing the cytosolic concentration of osmo-active molecules. The three major categories of osmo-active molecules modulated in osmoregulation are ions, sugars and amino acids [37]. Osmoregulation via ions predominantly occurs through K+ and its counterions [37]. Gene transcripts involved in transports of K+ and H+ did not show different change trends of abundance in the two comparisons of MA-T1 vs. MA-T2 and AE-T1 vs. AE-T2 (Fig. 7A). A gene transcript encoding sugar carrier protein C relating to glucose transmembrane import was down-regulated in MA-T1 vs. MA-T2 but up-regulated in AE-T1 vs. AE-T2 (Fig. 7A). As for amino acid, proline is particularly targeted in osmoregulation [37]. Its concentration is regulated at the levels of biosynthesis, transport and catabolism. Two gene transcripts encoding pyrroline-5-carboxylate synthase (P5CS) were down-regulated in MA-T1 vs. MA-T2 but up-regulated in AE-T1 vs. AE-T2 (Fig. 7A). The rate-limiting step of proline biosynthesis is catalyzed by P5CS [38].
In flower closing process (from T3 to T4), most gene transcripts involved in water transmembrane transport, transport of K+ and H+ and sugar transport were down-regulated in AE-T3 vs. AE-T4, but did not show significant differential expression in MA-T3 vs. MA-T4 (Fig. 7B).
Cell wall-related activity
Cell walls are dynamic structures with rigidity and extensibility. The extensibility of cell walls allows cell expansion which is driven by a strong intracellular turgor pressure [39]. Several gene transcripts whose proteins involve regulating the extensibility of cell walls were detected to show opposite change trends of abundance in the two comparisons of MA-T1 vs. MA-T2 and AE-T1 vs. AE-T2. In flower opening process (from T1 to T2), two gene transcripts encoding fasciclin-like arabinogalactan proteins (AGP2 and AGP12) and a gene transcript encoding xyloglucan endotransglucosylase/hydrolase (XTH) were down-regulated in MA-T1 vs. MA-T2 but up-regulated in AE-T1 vs. AE-T2 (Fig. 8A). AGPs are structural proteins in the cell wall which can regulate wall expansion[40]. XTHs can modulate the interactions between cellulose microfibrils (CMFs) and hemicellulose xyloglucans (XyGs) to strengthen or loosen cell wall [39, 41, 42].
Cell wall degradation is part of petal senescence [19]. In flower closing (senescing) process (from T3 to T4), a large number of gene transcripts relating cell wall activity were detected to present different change trends of abundance between MA and AE. Thirteen gene transcripts were normally regulated in MA-T3 vs. MA-T4 but down-regulated in AE-T3 vs. AE-T4 (Fig. 8B). These gene transcripts included four transcripts encoding AGPs (AGP12 and AGP14), a transcript encoding expansin-A16-like, three transcripts encoding cellulose synthase A (CESA, mediating the biosynthesis of CMFs) [43], a transcript encoding pectin methylesterase (PME, decreasing the crosslinking between pectins and calcium ions to soften of the wall) [44, 45] and four transcripts encoding pectate lyases (PLs, depolymerizing pectic chains so as to loosen cell walls) [46]. Flower senescence was also associated with an increase in the transcript abundance of putative genes encoding β-galactosidase (BGAL), β-glucosidase (BGLU) and polygalacturonase (PG) which are putatively involved in cell wall degradation [19, 47]. In flower senescing process (from T3 to T4), a gene transcript encoding BGAL and two transcripts encoding BGLU were up-regulated in AE-T3 vs. AE-T4 but did not show significant differential expression in MA-T3 and MA-T4. Two gene transcripts encoding PG-like were down-regulated in MA-T3 vs. MA-T4 but up-regulated in AE-T3 vs. AE-T4 (Fig. 8B). In constant, three gene transcripts encoding polygalacturonase inhibitor-like (PGI-like) were up-regulated in MA-T3 vs. MA-T4 but down-regulated in AE-T3 vs. AE-T4 (Fig. 8B).
Execution of senescence
Flower senescence is a highly developmentally regulated process which is induced by programmed cell death (PCD) of petal cells.
Caspases involved in cell death signal transduction of animal cells are cysteine endo-proteases [48]. In flower senescing process (from T3 to T4), several gene transcripts whose proteins are cysteine proteases were detected to perform different change trends of abundance between MA and AE. Two of them encoding senescence-specific cysteine protease SAG39-like were down-regulated in MA-T3 vs. MA-T4 but up-regulated in AE-T3 vs. AE-T4 (Fig. 9). Another transcript encoding metacaspase-9 did not show significant differential expression in MA-T3 vs. MA-T4 but was up-regulated in AE-T3 vs. AE-T4 (Fig. 9). Flower senescence is accompanied by a series of metabolic processes, such as sugar metabolism, protein degradation, remobilization and transport of mineral ions, fatty acid degradation and phospholipid degradation. Several gene transcripts associated with these processes were detected to show up-regulation in AE-T3 vs. AE-T4 but show little transcript changes in MA-T3 vs. MA-T4 (Fig. 9). These gene transcripts included a transcript encoding sucrose synthase 4, three transcripts encoding proteases (cysteine protease and aspartic protease) relating to protein degradation [19], a transcript encoding isocitrate lyase-like (ICL-like) which functions in the fatty acid degradation [19, 49], two transcripts encoding phospholipase D alpha 1 and non-specific phospholipase C2 which function in the phospholipid degradation [19], and a transcript encoding purple acid phosphatase 23 which can remobilize phosphates [19]. Syntaxin are required for the fusion of transport vesicles to target membranes [50]. The two gene transcripts that encode syntaxins were down-regulated in AE-T3 vs. AE-T4 but did not show significant differential expression in MA-T3 vs. MA-T4 (Fig. 9).
RNA-seq validation and expression analysis using qRT-PCR
To verify the reliability of the transcriptome sequencing data, eleven DEGs were selected randomly to perform qRT-PCR on the same sample as those in transcriptome sequencing. The results showed that the expression levels of the selected genes were basically consistent with the RNA-seq results (Fig. 10), which confirmed the validity of the RNA-seq data and the following analysis.