Transcriptome sequencing and expression profiling of genes involved in daylily (Hemerocallis citrina Borani) flower development

Abstract

Daylily (Hemerocallis citrina Borani), a widespread ornamental plant around the world, has been considered as an ideal model to study flower development recently. To initiate the research to comprehensively investigate the molecular mechanism of the flower development of daylily, we performed a transcriptome sequencing study to establish the transcriptional profiles of this plant at different flower developmental stages. Briefly, we divided the flower development of daylily into six stages, and performed transcriptome sequencing for the samples collected from each stage. A total of 56,598 unique transcripts were obtained, of which 35,424 ones (about 62.6%) encode proteins with known functions based on homology searches against well-studied model plants. The expression profiles of these genes showed that 7,395 genes were differentially expressed during flower development. Further functional enrichment analysis for these differentially expressed genes revealed that they are enriched in the processes of translation, oxidation-reduction, and cell wall metabolism, suggesting the critical roles of these processes in daylily flower development.

Keywords:

• Hemerocallis citrina
• RNA-Seq
• flower development
• ROS
• cell wall

Introduction

Daylily (Hemerocallis citrina Baroni), a popular herbaceous perennial plant, is widely cultivated around the world as an ornamental plant. Daylily has distichously arranged long, linear leaves, and showy flowers, which are borne on short pedicels [1]. The flower buds are commonly used as a vegetable in Asia [2]. Recent studies showed that the flowers of daylily have more antioxidant compounds than other plants, making it an excellent nutraceutical food and sparking a potential for medical use [3,4].

In the past few decades, cumulative studies have been performed to explore the mechanism of daylily flower development, from morphology to biochemistry, physiology and genetics. More and more evidence shows that daylily is an ideal model for plant flower development [5]. Bieleski and Reid [6] reported the increase of sugar and ion efflux in daylily flower senescence and they also confirmed the significant sugar content changes in daylily flower development. Panavas et al. [7] reported that the ROS increase due to reduced effectiveness of protective enzymes, leads to programmed cell death (PCD), which plays important roles in daylily flower development. They also showed that up-regulation of activities of wall-based enzymes with similar functions, such as polygalacturonase and β-galactosidase, may contribute to flower senescence [8]. In addition, Nitta et al. [9] reported that a daylily gene, similar to TOC1 in Arabidopsis thaliana, controls flower opening; and a major locus controls the timing of flower closing. The molecular mechanism by which these two genetic loci regulate the flower development [10], however, is still unclear. Here, to explore the molecular mechanism of daylily flower development from the perspective of genomics, we studied the transcriptional profiles of daylily flower at different developmental stages.

In this study, with the goal to identify essential genes involved in the regulation of daylily flower development, transcriptome sequencing was performed to compare the transcriptional profiles of six flower developmental stages. It showed that the majority of the up-regulated genes are involved in the processes of translation, oxidation-reduction and cell wall metabolism, suggesting the critical roles of these processes in flower development.

Materials and methods

Plant material and growth

Daylily plants (kindly supplied by Prof. Chen Wang at College of Life Science and Technology, Harbin Normal University), were cultivated in the experiment station of Harbin Normal University. All samples were collected on June 26th, 2019, from one daylily clone. The samples were classified into six stages based their size: stage A (2.1-2.5 cm), stage B (4.1-4.5 cm), stage C (5.1-5.5 cm), stage D (7.1-8.0 cm), stage E (full opening) and stage F (flower senescence). The details are shown in Supplemental Figure S1. Two samples were collected for each stage. All samples were flash frozen in liquid nitrogen and then stored at −80 °C.

Construction and sequencing of RNA-seq library

Total RNA was extracted using the RNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA) following the manufacturer’s instructions as our previous study described [11]. RNA samples were sent to BGI-Shenzhen Co. Ltd (Shenzhen, China), and RNA-seq was performed on the BGI-Seq500 platform according to the manufacturer’s instructions to generate 150-bp paired-end reads as described in a previous study [12].

Sequences assembly and annotation of RNA-seq data

Adapter sequences and low-quality reads were removed from raw sequences to produce clean data, which was performed by BGI Co. Ltd. These clean data were submitted to NCBI SRA database (Accession no.: SRP243199). The clean reads were then assembled de novo into contigs using Trinity with default parameters [13]. The contigs were further clustered into unigenes using CORSET with default parameters as described [14]. All unigenes were BLAST-searched against combined databases of Arabidopsis thaliana, rice (Oryza sativa), maize (Zea mays), orchid (Phalaenopsis equestris) and asparagus (Asparagus officinalis) protein sequences using BLASTX program for functional annotation (E-value cutoff: 1e-5) [15]. Next, daylily unigenes were assigned with functional annotation based on their corresponding homologs in the combined database, including Gene Ontology (GO) and KOG (EuKaryotic Orthologous Groups) annotations [16,17]. In addition, daylily unigenes were scanned using the iTAK pipeline to identify transcription factors, as Jin et al. [18] described.

Expression quantification and differential expression analysis

All RNA-seq reads were mapped to daylily unigenes using TopHat software [19], and the expression levels of each unigene were estimated as FPKM values (fragments per kilobase of exon per million fragments mapped) using Cufflinks software [20]. Differentially expressed unigenes were identified using Cufflinks software with the time-course model as the protocol described. The GO annotation results of these differentially expressed unigenes involved in flower development were retrieved, and GO functional enrichment analyses were performed using the topGO package in R [21].

Results and discussion

Sequencing and de novo assembly of the daylily transcriptome

In order to establish the transcriptional profiles of daylily at different flower developmental stages, we collected two samples from each of the classified six developmental stages. RNA was extracted from each sample, and RNA-seq was performed by using next-generation sequencing. In total, around 291 million reads were generated from these twelve samples. We have submitted the data to NCBI SRA database (Accession numbers: SRP243199). After these data were de novo assembled by using Trinity software, and the redundant transcripts were eliminated using CORSET, we obtained 56,598 transcripts (genes) from these reads. The mean length of these transcripts was 1,246 bp, and the N50 value was 1,735 bp (Table 1). To determine their identities, we performed a BLAST search to compare these transcripts with the database established by combining the protein sequences of the well-studies plant organisms: Arabidopsis, rice, maize, orchid and asparagus. It showed that, of these transcripts, 35,424 ones (about 62.6%) were identified as homologs with significant hits to the well-studied plant organisms. As expected, it also showed that the longer transcripts had more hits than the short ones (Figure 1), implying that the longer transcripts were more probably annotated. In other words, the identities of the longer transcripts have more confidence, which is consistent with previous reports in other plants [11,12].

Figure 1. Length distributions of daylily unique transcripts.

Note: The X-axis lists the length of transcripts, and Y-axis lists the numbers of transcripts. Red bars represent the total number of transcripts at specific length ranges, and the adjacent blue bars represent the ones annotated.

Table 1. Summary of de novo assembled the daylily transcriptome.

Data type	Number
Total sequence	56598
Number of sequences in 201-500 bp	12106
Number of sequences in 500-1000 bp	19220
Number of sequences more than 1000 bp	25272
Minimal length (bp)	201
Maximal length (bp)	15161
N50 (bp)	1735
Average length (bp)	1246

Functional annotation analysis of daylily unigenes

To identify the biological processes that the transcripts are involved in, we further used the GO terms and KOG annotations of their homologs to assign functions to the transcripts of daylily. It showed that 23,019 genes (about 40.67%) can be properly grouped into the following GO processes (Figure 2 and Supplemental Figure S2): DNA metabolic process (GO:0006259, 699 unigenes), lipid metabolic process (GO:0006629, 689 unigenes), transport (GO:0006810, 93 unigenes), cell–cell signaling (GO:0007267, 221 unigenes), multicellular organism development (GO:0007275, 170 unigenes), circadian rhythm (GO:0007623, 54 unigenes), cell death (GO:0008219, 25 unigenes), and flower development (GO:0009908, 42 unigenes). These biological processes provide valuable information for understanding daylily flower development. For example, DNA and lipid metabolic processes imply DNA and lipids undergo tremendous changes in daylily flower development, which is consistent with previous reports [4,8]. Similarly, transport, cell–cell signaling, circadian rhythm and cell death are all considered critical processes in daylily flower development, which have also been reported in previous studies but with insufficient evidence. The grouping of a significant number of transcripts into the GO term “flower development” shows the significance of these genes in daylily flower development. Meanwhile, KOG annotations were also assigned to all the transcripts, and it showed that the three top KOG terms were signal transduction mechanisms (2,930 genes, 13.98%), general function prediction only (2,887 genes, 13.77%), and posttranslational modification, protein turnover and chaperones (2,171 genes, 10.35%), which are consistent with the GO annotations. In addition, plant transcript factors (TFs) were also identified among these transcripts; in total 2281 TF genes were identified and classified into 74 TF gene families according to the classification in PlantTFDB. The most abundant TF family was the MYB (166 genes), followed by the AP2/ERF (111 genes), NAC (110 genes), C2H2 (105 genes) and bHLH (104 genes) (Supplemental Figure S3). These TF families have been shown to be key in the regulation of flower development in well-studied model plants [22–24], suggesting that the daylily has similar regulatory mechanisms.

Figure 2. GO annotation results of daylily unique transcripts.

Note: GO terms are listed along the X-axis, and the numbers of transcripts grouped in the GO terms are shown on the Y-axis. Red bars represent Biological processes; Green bars represent Molecular functions; and blue bars represent cellular components.

Expression analysis of genes involved in daylily flower development

To determine which genes play a role in the daylily flower development process, five comparisons between adjacent flower development stages were performed using Cufflinks software. In total, 7,395 differentially expressed genes (DEGs) were identified with a false discovery rate of 0.01. Specifically, there were 129 genes involved in the transition from stage A to B, 348 genes involved in stage B to C, 1551 genes involved in stage C to D, 3964 genes involved in stage D to E, and 2673 genes involved in stage E to F. These genes were clustered into four groups based on their expression levels (Figure 3). Group I was highly expressed at stages A, B, C and D, but less expressed at stages E and F; Group II showed the opposite trend to Group I; Group III had low expression at all stages except stage F, suggesting their special functions at stage F; whereas Group IV was highly expressed at stage D, E and F. These patterns are consistent with previous studies [9], which considered flower opening and flower senescence as the most import processes during daylily flower development with the identification of the largest number of genes related to the regulation of the two processes.

Figure 3. Heatmap showing expression profiles of differentially expressed genes during daylily flower development.

Note: The heatmap was generated using R gplots package; the FPKM values of each gene were normalized and used as input. The color transition from dark blue to dark red represents the expression level of genes from the lowest to the highest.

To determine the biological processes of these DEGs, GO annotation enrichment analysis was performed using the software package topGO. The analysis had 128 GO terms significantly enriched (Supplemental Table S1). It has been reported that programmed cell death (PCD) plays important role in daylily flower development [25], which involves many complex biochemical pathways, including protein synthesis (translation) and proteolysis (protein disassembly). PCD is involved in cell death in flower tissues and is considered to be the major role for regulating flower development states in daylily [7,8,25,26]. Two biological processes related to PCD were identified in this study. They are cell wall metabolism and reactive oxygen species (ROS) metabolism. Our finding suggests that daylily autonomously regulates enzymes that modify cell walls, such as polygalacturonases and b-galactosidases. Also, it regulates ROS related enzymes, such as oxidoreductase. In addition, it promotes the PCD process, such as inhibition of translation, protein disassembly [27], and DNA digestion [28]. In the “biological processes” category, the main terms, including translation (GO:0006412), oxidation-reduction process (GO:0055114), cell wall organization or biogenesis (GO:0071554), response to auxin (GO:0009733), protein complex disassembly (GO:0043241) and response to auxin (GO:0009733), suggest that the genes involved in these processes potentially contributed to daylily flower development. These findings have well confirmed previous suggestions; for example, genes involved in cell-wall metabolism (GO:0071554) have obvious changes of expressional profiles during flower development (Figure 4), implying their regulatory function in flower opening and senescence. In addition, GO enrichment analysis were enriched in terms including hydrolase activity (GO:0016798), oxidoreductase activity (GO:0016491), transporter activity (GO:0005215), ribosome (GO:0005840) and cell wall (GO:0005618) from the categories of “molecular function” and “cellular component”, which also confirmed their function in the PCD process in flower development. However, the major locus that regulates daylily flower development, as Nitta et al. [9] described, is still unknown. It is reasonable to speculate that some transcription factors related to flower development regulation might be encoded by this locus.

Figure 4. Heatmap showing expression profiles of differentially expressed genes involved in cell wall metabolism during flower development.

Note: The genes that were grouped into the GO term (GO:0071554) were identified as cell–wall metabolism genes. Their expressional levels (FPKM values) were normalized, and used as input to generate the cluster plot. Color annotations are the same as those in Figure 3.

Conclusions

In the present study, we performed RNA-seq to establish the transcriptional profiles of different flower developmental stages of daylily. The transcriptome sequencing data were assembled into 56,598 transcripts, and 23,019 genes were properly annotated. Importantly, 7,395 genes were identified as differentially expressed genes during daylily flower development. These genes were enriched in translation, oxidation reduction process, cell wall metabolism, and ROS metabolism, implying their critical roles in daylily flower development.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Natural and Science Foundation of China (grant number 31302019) and China Postdoctoral Science Foundation (grant number 2016T90307 and 2015M571430).