Research article | Open | Published:
De novo transcriptome analysis of petal senescence in Gardenia jasminoides Ellis
BMC Genomicsvolume 15, Article number: 554 (2014)
The petal senescence of ethylene insensitive species has not been investigated thoroughly while little is known about the temporal and tissue specific expression patterns of transcription factors (TFs) in this developmental process. Even less is known on flower senescence of the ornamental pot plant Gardenia jasminoides, a non climacteric flower with significant commercial value.
We initiated a de novo transcriptome study to investigate the petal senescence in four developmental stages of cut gardenia flowers considering that the visible symptoms of senescence appear within 4 days of flower opening. De novo assembly of transcriptome sequencing resulted in 102,263 contigs with mean length of 360 nucleotides that generated 57,503 unigenes. These were further clustered into 20,970 clusters and 36,533 singletons. The comparison of the consecutive developmental stages resulted in 180 common, differentially expressed unigenes. A large number of Simple Sequence Repeats were also identified comprising a large number of dinucleotides and trinucleotides. The prevailing families of differentially expressed TFs comprise the AP2/EREBP, WRKY and the bHLH. There are 81 differentially expressed TFs when the symptoms of flower senescence become visible with the most prevailing being the WRKY family with 19 unigenes. No other WRKY TFs had been identified up to now in petal senescence of ethylene insensitive species. A large number of differentially expressed genes were identified at the initiation of visible symptoms of senescence compared to the open flower stage indicating a significant shift in the expression profiles which might be coordinated by up-regulated and/or down-regulated TFs. The expression of 16 genes that belong to the TF families of WRKY, bHLH and the ethylene sensing pathway was validated using qRT – PCR.
This de novo transcriptome analysis resulted in the identification of TFs with specific temporal expression patterns such as two WRKYs and one bHLH, which might play the role of senescence progression regulators. Further research is required to investigate their role in gardenia flowers in order to develop tools to delay petal senescence.
One of the main issues that floriculture industry has to confront is flower senescence, a term that signifies all those processes that follow physiological maturity and lead to the death of a whole plant, tissue or cell and represents the last stage of flower development. A series of events take place during flower senescence resulting in highly regulated, genetically programmed and developmentally controlled morphological, physiological and biochemical changes. By the end of this process, the flower petals wilt, lose colour and abscise .
Plants are classified as climacteric or non-climacteric, depending on their ethylene production rate. Climacteric flowers such as carnation and petunia are characterized by the climacteric ethylene and respiration rate production that promotes senescence, while treatment with exogenous ethylene results in acceleration of senescence . In addition, an inhibitor of prolyl 4 hydroxylase activity, pyridine 2,4 dicarboxylate, suppressed the climacteric ethylene production in cut carnation flowers . In contrast to carnation and petunia, the flowers of lilies (Lilium spp.), gladiolus, tulip (Tulipa spp.), iris and morning glory (Ipomoea nil) are classified as non-climacteric since they are not responsive to ethylene and exogenous application has little or no effect on petal senescence . Ethylene acts as a trigger to flower senescence in long-lived flowers, possibly as a mechanism to terminate flower life after pollination. In short-lived flowers on the contrary, such a mechanism would not be beneficial since the life of individual flowers is very short .
Genes that are associated with ethylene sensing and sensitivity have also been found in ethylene insensitive plants such as the two homologues of Arabidopsis ethylene receptor ERS1 which were found in gladiolus . Over the last years several signalling components have been found to be associated with petal senescence in climacteric species, such as the F-box proteins EBF1/EBF2 that recognize and degrade EIN3 through the ubiquitin/26S proteasome pathway [6, 7]. They also interact with EIN2 to function in the modulation of ethylene signalling  as well as with EIN3, EIL1 and ASK1 which is a component of the SCF ubiquitin ligase complex.
In Arabidopsis petals, 316 TFs showed differential expression patterns of which 130 changed expressions only in petal but not in leaf or silique senescence . These TFs are members of 47 gene families while the three most represented families with up-regulated patterns of expression were the ERF, NAC and WRKY .
The EIN/EIL transcription factors (TFs) are targeting genes with an Ethylene Responsive Element in their promoter. They belong to the AP2/EREBP-type TFs and are correlated with stress response [10–12], ripening and senescence [13, 14]. They are also correlated with the effect of sucrose in delaying senescence .
In addition to ERFs there are other families of TFs associated with senescence such as the WRKYs and basic helix loop helix (bHLH). The WRKY TF family constitutes a large family with at least one WRKY domain of approximately 60 amino acids. They are unique to the plant kingdom and some of its members have a very important role in the regulation of leaf senescence, response of the plant to bacterial infection, signaling pathways and many biotic and abiotic stresses [16–18]. WRKY target genes with a W-box (TTGAC or TTGACC/T) in their promoter [19, 20] regulating, among others, the expression of senescence associated with related genes [16, 21]. Members of the WRKY family might constitute a hub transcription factor during senescence via mediation of jasmonic acid and salicylic acid signaling [22, 23].
The bHLH proteins can be found both in mammals and plants  and are known to interact with other proteins and transcription factors such as MYB [25, 26]. They were also associated with a series of developmental phenomena in the plant, such as the development of root hair and leaf trichomes [27, 28], cell proliferation and cell differentiation [24, 27].
However, little is known about the expression of TFs in petal senescence of ethylene insensitive species. In alstroemeria 21 TFs were identified among 2000 ESTs , while other TFs were identified in daffodil  and iris , while in morning glory (Ipomoea nil) suppression of a leucine-rich repeat transmembrane protein kinase displayed accelerated petal senescence .
Recently, there is an increased interest for the ornamental pot plant Gardenia jasminoides Ellis, since apart from its wide use as a floral plant, it is of high pharmaceutical value considering that contains the substance geniposide that can be transformed into the anti-inflammatory and anti-angiogenesis agent genipin [32–34]. The physical map of gardenia is not available yet, while molecular studies are restricted to discriminate gardenia species for systematic reasons or to study the phytogeography of the wild and commercial populations [35–40]. Gardenia is considered a non-climacteric species with flowers that do not produce detectable levels of ethylene therefore is not included in the list of climacteric flower crops . We initiated a study to investigate the petal senescence of cut gardenia flowers at the transcriptome level considering that the visible symptoms of senescence appear within four days. Towards this direction we used Next Generation Sequencing Technologies, in particular Illumina HiSeqTM 2000 for de novo sequencing and characterization of gardenia petals transcriptome at four stages of senescence progression.
The senescence progress of gardenia flowers is completed within four days, therefore it was divided in four developmental stages. The first stage (A) comprises closed buds ready to open and the second stage (B) open flowers with the outer whorl of petals at a 90° angle to the flower stalk. At the third developmental stage (C) the flowers are fully open, while at the fourth stage (D) the yellow discoloration of petals was initiated as well as the appearance of brown spots (Figure 1). Total RNA was extracted from each development stage and used for de novo transcriptome sequencing with Illumina HiSeqTM 2000.
cDNA sequence generation, de novo assembly and quantification of gene expression
Sequencing of the gardenia transcriptome resulted in a total of 50,335,672 reads that were obtained after cleaning the raw data (Additional file 1). De novo assembly of these reads resulted in 102,263 contigs with mean length of 360 nt that generated 57,503 unigenes with a mean length of 796 nt. These were further clustered into 20,970 clusters and 36,533 singletons. Quantification of gene expression for the four developmental stages of gardenia flower senescence, the uniquely mapped reads for each stage and the corresponding unigenes are shown in Additional file 2. The E-value distribution of the top hits in NCBI non-redundant (Nr) database, the similarity distribution as well as the species distribution are included in Additional file 3.
Gene Ontology annotation of gardenia transcriptome
The gardenia genes were classified according to the Gene Ontology annotation in three categories, namely Biological Process, Cellular Component and Molecular Function (Figure 2).
In the category “Biological Process” most genes were associated with “biological regulation”, “cellular process”, “metabolic process”, “regulation of biological process”, “response to stimulus”, “single-organism process” and “signaling”. For the category “Cellular Component”, the dominant subcategories were those genes associated with “cell” and “cell part”, “membrane”, “organelle” and “organelle part”. In the category “Molecular Function”, “binding” and “catalytic activity” were the prevailing terms. Overall, in the present study, genes that belong to the categories “cellular process”, “metabolic process”, “response to stimulus” (biological process), “cell”, “cell part”, “organelle” (cellular component) and “binding”, “catalytic activity”, “metabolic activity” (molecular function) were the most highly represented. The above mentioned gene categories were the most highly represented in other transcriptome studies as well, such as in safflower and in chrysanthemum [41, 42].
Clusters of Orthologus Groups (COG) annotation
Search against the COG database resulted in the classification of 13,462 unigenes in 24 COG categories. As it can be seen in Figure 3, the “General Function prediction only” represents the largest group with 4,333 unigenes followed by “Transcription” (2,227 unigenes), “Replication, recombination and repair” (2,089 unigenes), “Posttranslational modification, protein turnover, chaperones” (1,984 unigenes) and “Signal transduction mechanisms” (1,782 unigenes).
Differentially expressed unigenes
The unigenes that are differentially expressed in pairs of two consecutive stages of senescence progress are shown in Table 1. As the gardenia flower develops from closed bud of stage A to early open flower of stage B, 4,360 unigenes are differentially expressed, of which 3,402 are down-regulated and 958 are up-regulated (Table 1). Among them, there are 3,162 unique unigenes which are differentially expressed only in these two developmental stages but not in the other consecutive stages (Figure 4). The progression of the partially open into a fully open flower at stage C is characterized by 2,464 differentially expressed unigenes of which 1,654 are down-regulated and 810 are up-regulated (Table 1). Among them, there are 1,129 unique, differentially expressed transcripts (Figure 4). The developmental stage D with visible symptoms of senescence in comparison to stage C is characterized by 2,643 differentially expressed unigenes of which 1,061 are down-regulated and 1,582 are up-regulated (Table 1). Among them, there are 1,654 unique, differentially expressed transcripts (Figure 4).Moreover, 682 differentially expressed transcripts are common between the stages A and B and those of B and C (Figure 4). Furthermore, 473 differentially expressed transcripts are found in stages B and C, and C and D, while 336 differentially expressed transcripts are common in the stages A and B, and C and D. The comparison of the consecutive developmental stages resulted in 180 common, differentially expressed unigenes. Among them, 21 unigenes were up-regulated, while 66 were down-regulated in all four developmental stages.
As petal senescence progresses, the number of down-regulated unigenes decreases, while that of up-regulated increases only at the stage D of visible symptoms of senescence (Table 1).
The comparison between stage A and the other developmental stages B, C and D indicates that as senescence progresses, the number of differentially expressed unigenes increases (Table 2). Specifically, 4,360, 7,658 and 10,401 unigenes are differentially expressed between stages A and B, C and D, respectively, indicating significant variation in gene expression profiles among developmental stages.The functional classification of the differentially expressed unigenes in the four developmental stages was performed as shown in Figures 5, 6 and 7.
Identification of simple sequence repeats
A total of 9,549 SSRs were indentified in 7,641 (13.3%) transcripts of gardenia (Additional file 4). Among them, 1,398 (18.3%) transcripts contained more than one SSR, while the largest category of the SSRs was mononucleotides (3129, 32.7%), followed by trinucleotides (2,937, 30.7%) and dinucleotides (2,853, 29.8%). The A/T (30.9%) accounted for the 94.5% of the mononucleotide repeats, the AG/CT (19.2%) accounted for the 69.5% of dinucleotides, while the AAG/CTT (8%), AGG/CCT (7.7%) and AAC/GGT (3.5%) accounted for the 62.4% of the trinucleotides. A small number of tetranucleotides (230, 2.4%), pentanucleotides (209, 2.2%) and hexanucleotide (190, 2.0%) were also identified.
Differentially expressed transcription factors during petal senescence
The 202 differentially expressed unigenes that encode TFs within the four developmental stages are shown in Figure 8. Among them, there are 107 TFs with differential expression patterns during progression to stage B. The prevailing family of transcription factors is the “AP2/EREBP” family with 11 unigenes, followed by the “WRKY” transcription factor family with 7 unigenes and the “bHLH” TFs with 6 unigenes, while there are 44 unknown/uncharacterized TFs (Figure 9A-B). Additional differentially expressed TFs include the “GATA” and “GRAS” with 5 unigenes, the “MYB” family with 4 unigenes and the “heat shock / stress”, “MADS”, “homeobox-leucine zipper”, “lycine-specific dimethylase” and “YABBI” families with 3 unigenes.As senescence progresses from partial to fully open flower, 58 TFs are differentially expressed comprising the “WRKY” family with 6 unigenes, followed by the “MYB” and the “AP2/EREBP” family with 4, the “bHLH” TFs with 3, and the “auxin response”, “GATA” and “kinesine” with 2 unigenes (Figure 9B-C).There are 81 differentially expressed TFs when the symptoms of flower senescence become visible at stage D. The most prevailing TF families include the “WRKY” family with 19 unigenes, the “bHLH” with 10, the “AP2/EREBP” with 12 while “auxin response” and “DNA binding” comprise 4 and 3 unigenes, respectively (Figure 9C-D).The clusters a, b, e and h mainly comprise WRKY TFs and their expression levels are higher in stage D while the cluster c mainly comprise bHLH TFs and their transcript abundance decreases after stage B (Figure 8). The GRAS and GATA TFs are grouped in cluster d and their expression decreases during senescence progression. AP2/EREBP TFs are mainly present in five clusters, b, e, f, g and i showing different patterns of expression. In clusters b and e the expression of the unigenes is higher in stage D, in clusters f and g is higher in the intermediate stages B or C and in cluster i the expression of the unigenes is high at the first stage and then decreases throughout senescence (Figure 8).
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway mapping
The 21,614 unigenes were assigned to 128 KEGG pathways. The most abundant pathways were ‘Metabolic pathways’ (5,001), ‘Biosynthesis of secondary metabolites’ (2,497), ‘Plant-pathogen interaction’ (1,146), ‘Plant hormone signal transduction’ (990), ‘Spliceosome’ (750), ‘RNA transport’ (676), ‘Endocytosis’ (603) and ‘Starch and sucrose metabolism’ (578). The prevailing categories, based on transcript abundance, between stages A and B were ‘Metabolic pathways’, ‘Biosynthesis of secondary metabolites’, ‘Ribosome’ and ‘Plant-pathogen interaction’, while between stages B and C were ‘Metabolic pathways’, ‘Biosynthesis of secondary metabolites’, ‘Starch and sucrose metabolism’, ‘Plant hormone signal transduction’ and ‘Plant-pathogen interaction’. The prevailing categories between stages C and D included ‘Metabolic pathways’, ‘Biosynthesis of secondary metabolites’, ‘Plant-pathogen interaction’, ‘Plant hormone signal transduction’ and ‘Starch and sucrose metabolism’, while between the stages A and D included ‘Metabolic pathways’, ‘Biosynthesis of secondary metabolites’, ‘Plant-pathogen interaction’, ‘Plant hormone signal transduction’ and ‘Ribosome’.There are 19 WRKY unigenes which are differentially expressed when the visible symptoms of senescence appeared, while 13 among them are unique to this developmental stage indicating specific involvement on the initiation of senescence symptoms. Only two WRKY and one bHLH unigenes showed a progressive increase and decrease, respectively in expression throughout petal senescence (Figure 10).
Ethylene sensing in the non-climacteric flower of gardenia
Differentially expressed unigenes were assigned to KEGG pathways associated with ethylene sensing (Figure 11). The CTR1 and the MPK6 unigenes are up-regulated, while the EIN3 and the ERF1/ERF2 are down-regulated during the transition from the stage A of closed buds ready to open to stage B of open flowers with the outer whorl of petals at a 90° angle to the flower stalk (Figure 11). As senescence progresses to fully open flowers, two CTR1 unigenes are still up-regulated, while one is-down regulated. In addition, two EBF1/EBF2 unigenes are up-regulated (Figure 11). When visible symptoms of senescence appear in the petals, four ERF1/ERF2 unigenes were up-regulated (Figure 11). These results indicate progressive up-regulation of ethylene sensing components during senescence development.
Validation of expression of selected unigenes using real time PCR
The patterns of expression of selected unigenes during the four stages of petal senescence were determined by real-time PCR using two reference genes, actin (CL715.Contig1) and PP2A (Protein Phosphate 2A) (Unigene23262). Two WRKY and one bHLH unigenes were up- and down-regulated, respectively throughout senescence showing similar to RNAseq patterns of expression (Figure 10). In addition, the transcript abundance of five differentially expressed WRKYs (Figure 12) and three differentially expressed bHLHs (Figure 13) were also determined during senescence. Their patterns of expression are similar to those determined using RNAseq indicating that quantification of expression profiles with this technological platform might be considered reliable. Moreover, five ethylene sensing components such as CTR1, MPK6, EIN3 and two ERF1/2 were further characterized by quantifying their mRNA levels during petal senescence with real-time PCR (Figure 14). Each one of the five ethylene signal transduction pathway components showed a distinct pattern of expression confirming the RNAseq analysis (Figure 14).
A large scale transcriptome analysis of petal senescence in a non-climacteric species, Gardenia jasminoides Ellis, was performed. Four distinct developmental stages of senescence were analyzed leading to the identification of 57,503 unigenes which were annotated against known sequences in various databases, whereas differentially expressed unigenes were also identified. The relative transcript levels of 16 differentially expressed unigenes were further validated by using real time PCR.
Annotated unigenes in all databases (NR, NT, Swiss-Prot, KEGG, COG, and GO) comprised 68.6% of the total number of unigenes, which should be considered adequate taking into account the lack of sequencing data on gardenia species. In other similar studies in other species, such as bamboo and cucumber, the annotated unigenes represented 55% and 72% of the total number of unigenes, respectively [43, 44]. The E-value distribution of the top hits in NR database also showed high homology, even in the absence of a reference genome.
The functional classification of gardenia genes according to GO annotation identified genes associated to ‘Biological Process’, ‘Cellular Component’ and ‘Molecular Function’ categories. Genes annotated to ‘Biological process’ were identified during flower development of Eustoma grandiflorum  and in three genotypes of Eichornia paniculata, while those involved in metabolic process were also identified in carnation, snapdragon, alstroemeria, tobacco and Arabidopsis [9, 29, 47–49] indicating significant metabolic activity during petal senescence progression.
Gene categories associated with ‘Cellular Component’ were also identified during flower senescence of bamboo  implicating the “membrane” associated genes to alterations in membrane permeability due to changes in composition and structure of the lipid bilayers . Gene categories associated with ‘Molecular Function’ were also found to be highly represented during flower formation and petal senescence in safflower, carnation, orchids and wallflower [41, 51–53]. Search against the COG database that briefly shows the proteomic profile indicated similarities between gardenia petal senescence and Dendrocalamus latiflorus flower development .
Microsatellites or Simple Sequence Repeats (SSRs) are considered reliable and informative markers in plant research. In this study a total number of 9,549 SSRs were identified in 7,641 sequences. The majority of SSRs were mononucleotides followed by trinucleotides and dinucleotide repeats, as is the case in chickpea . These SSRs can be used for fingerprinting and breeding efforts within gardenia species.
A large number of genes are differentially expressed among the four different developmental stages of senescence. Among them, a plethora of transcription factors such as the WRKY, AP2/EREBP and bHLH were the most highly represented.
The WRKY family comprises the highest number (19) of differentially expressed unigenes in the final stage of gardenia petal senescence indicating significant involvement in this process compared to the other TF families. No other WRKY TFs had been identified up to now in petal senescence of ethylene insensitive species . It is of particular interest the fact that there are 1,582 differentially expressed genes at the initiation of visible symptoms of senescence compared to the open flower stage indicating a significant shift in the expression profiles, which might be coordinated by up-regulated TFs in these two developmental stages. Such candidate TFs might be the two WRKYs (‘CL7516.Contig2’ and ‘Unigene25021’) exhibiting expression profiles of up-regulation throughout petal senescence. These patterns were also confirmed using real time PCR (Figure 10). Moreover, another five WRKY TFs were shown to be up-regulated mainly at the last stage of senescence progression according to real time PCR analysis. The WRKY TFs are well known to be up-regulated during senescence in various organs in Arabidopsis, while petal and silique senescence was associated with the expression of more members of this gene family . In addition, WRKY TFs were also found to be highly abundant during bamboo flower development .
In gardenia petals, the bHLH unigenes are expressed in a stage specific manner and the expression pattern is distinct related to the other TFs while only one unigene (‘CL1446.Contig1’) is differentially expressed throughout senescence. The expression of this unigene is progressively down-regulated according to RNAseq and real-time PCR analysis (Figure 10). Moreover, there are ten differentially expressed unigenes at the stage of visible symptoms of senescence appearance; two were up- and eight were down-regulated (Table 3). Real time PCR analysis confirmed the expression profiles of three of them; two showed an increase in transcript abundance and one a decrease (Figure 13). The bHLH TFs are known to be associated with flower organogenesis, floral development and promotion of the proliferation inside the flower on a tissue-specific way [55–57]. They are also known to interact with other proteins and TFs such as the MYB families [25, 26] which are also differentially expressed throughout gardenia flower senescence (Figure 8).
Although gardenia is a non-climacteric flower, 12 AP2/EREBP TF unigenes were found to be differentially expressed with 11 of them up- and one down-regulated when the visible symptoms of senescence appeared. In addition to these positive TF regulators of ethylene response, other ethylene signaling components seem to be up-regulated in earlier stages of flower senescence such as the CTR1 and EBF1/EBF2 unigenes, which serve as negative regulators of ethylene response. These results suggest a pattern of expression for ethylene sensing and response genes in which the negative regulators are expressed early in flower senescence, while the positive regulators at later stages in order to possibly accelerate this developmental program.
The expression patterns of these five ethylene sensing components were validated using real time PCR suggesting involvement of ethylene in gardenia petal senescence.
Overall, the mRNA levels of several genes of interest by RNAseq analysis were further validated by real time PCR analysis indicating the reliability of RNAseq data.
The transcriptome analysis of gardenia petal senescence resulted in the identification of several TFs which might serve as targets for manipulation of their expression using genetic engineering approaches.
The de novo transcriptome analysis of Gardenia jasminoides Ellis resulted in the identification and functional classification of differentially expressed unigenes during petal senescence as well as in a large number of SSR markers. The assignment of these unigenes in KEGG pathways revealed potential involvement of ethylene sensing components in this developmental program. Moreover, differentially expressed transcription factors such as two WRKYs and one bHLH were identified showing specific temporal expression patterns which might play the role of senescence progression regulators. Their expression patterns as well as of other members of these gene families were further validated using real time PCR approaches. However, further research is required to investigate the specific role of transcription factors in gardenia flowers in order to develop molecular tools for delaying petal senescence.
Thirty five pot-plants of Gardenia jasminoides Ellis were cultivated in the greenhouse of the Laboratory of Floriculture at the Aristotle University of Thessaloniki, Greece. All plants were propagated by terminal shoot cuttings (7-8 cm in length), from a single mother plant that was provided by the cooperative “Gardenia Growers Group” of Magnesia, Greece. The basal portion (1 cm) of each cutting was dipped into 500 ppm of K-IBA solution and then planted in 2:1 peat:perlite substrate, in a fog propagation system (RH 95-98%). Two months later, the rooted cuttings were transplanted in 1,3 L pots filled with peat and grown identically in the greenhouse for one year, according to the commercial cultivation techniques.
Floral stages, RNA extraction and Illumina HiSeqTM 2000
Twenty five buds or flowers per stage were randomly selected from the above plants and harvested. Since the senescence progression of gardenia flowers is completed within four days, we divided the phenomenon in four developmental stages. The first stage (A) comprises closed buds ready to open and the second stage (B) open flowers with the outer whorl of petals at a 90° angle to the flower stalk. At the third developmental stage (C) the flowers are fully open, while at the fourth stage (D) the yellow discoloration of petals was initiated as well as the appearance of brown spots (Figure 1). Immediately after harvest, the petals were removed with a single-use blade and ground to a fine powder with liquid nitrogen. Total RNA was extracted according to the method of Bachem et al. 1998  and quantified using the IMPLEN GmbH Nanophotometer/NanoPhotometerTM Pearl. RNA samples were stored at −80°C.
The isolated RNA was pooled for each floral stage and 400 μg of the pooled RNA (per stage) were diluted in 2.5 volumes of absolute ethanol and 0.1 volume of sodium acetate (pH 5.2) in DEPC-treated ddH2O. We performed de novo transcriptome sequencing in the four samples, using Illumina HiSeq™ 2000 (Beijing Genomics Institute - BGI, Hong Kong).
Library preparation and sequencing
The RNA samples were treated with DNase I to ensure that are DNA-free. Subsequently the mRNAs were enriched by using the oligo (dT) magnetic beads and then fragmented into short pieces (about 200 bp). These fragments were used to synthesize the dsDNA using random hexamer primer for the first strand and DNA polymerase I for the second strand. Double stranded cDNA was purified with magnetic beads and subjected to end reparation and 3′ single adenylation. Sequencing adaptors were then ligated to the adenylated fragments, which were subsequently enriched by PCR amplification. The Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System was used to qualify and quantify the sample library during the QC step. The library products were sequenced using Illumina HiSeqTM 2000.
Data processing and de novo assembly
Raw reads that are obtained from the sequencer often contain low quality reads and/or adaptor sequences. Therefore, a preprocessing step of data cleaning was required to obtain the “clean reads” that would be used in the next steps of data analysis. This step included the adaptor removal as well as the application of a stringent filtering criterion to remove reads with >10% unknown bases and remove low quality reads, i.e. reads that comprise at least 50% low-quality bases. The reads that were obtained after the above mentioned filtering criteria were considered clean reads. These reads were used for the de novo assembly using Trinity  software that comprises three independent software modules, namely Inchworm, Chrysalis and Butterfly. Inchworm assembles the reads into longer transcripts in order to form contigs, Chrysalis groups the contigs into clusters that represent the transcriptional complexity of the gene and Butterfly reports the full-length transcripts for alternatively spliced isoforms. The output of Trinity are sequences called unigenes that can either form clusters in which the similarity among sequences is >70%, or unigenes that are unique and form the singletons. Lastly, all unigenes are aligned to protein databases including NR, Swiss-Prot, KEGG and COG (blastx, E-value <10-5) and the optimal result is used the sequence direction of unigenes. In the case of conflict among the results of the databases, the priority order of NR, Swiss-Prot, KEGG and COG is used. When the unigene is not aligned in any of the databases, ESTScan  was used to detect possible coding regions of the gene as well as identify the sequencing direction.
Quantification of gene expression, differentially expressed genes and gene annotation
The gene expression levels are determined by calculating the RPKM (Reads Per kb per Million reads) values using the formula RPKMA = 106 * CA (N * LA/103)-1, where CA is the number of reads that align uniquely to gene A, N is the total number of reads that align uniquely to all genes, and LA is the number of bases of this gene. Fold changes between (pairs of) stages were computed as the ratio of RPKM values.
Identification of differentially expressed genes in a pair of conditions 1 and 2 is based on the p-value:
where N1 and N2 represent the total number of reads and x and y represent the number of transcripts for a particular gene A, in conditions 1 and 2, respectively. False Discovery Rate (FDR) is the ratio of false positive over true differentially expressed genes and is set to a threshold of FDR ≤ 0.001, while the absolute value of log2Ratio was set to 1 or larger. Differentially expressed genes are then carried out into GO functional analysis and KEGG Pathway analysis.
Quantitative real-time PCR analysis
Expression analysis of the selected genes was assessed using quantitative RT- PCR. 2 μg of total RNA of each sample was converted into single stranded cDNA using TaKaRa PrimeScript Reverse Transcriptase (Cat. #2680A). The cDNA products were diluted 100-fold in DEPC ddH2O before used as templates. The reaction was performed on a Step One Plus Real-Time PCR System (Applied Biosystems, USA), using the Thermo Scientific Maxima SYBR Green/ROX qPCR Master Mix (2x) (Cat. #K0222). 20 μL of the reaction system contained: 10 μL of SYBR Green qPCR Master Mix (2×), 2 μL of each of the forward and reverse primers, 1 μL of water and 5 μL of cDNA template (1 ng/μL). The reaction conditions were performed as follows: 95°C for 10 min, followed by 40 cycles of 95°C for 15 sec, 60°C for 30 sec and 72°C for 30 sec. At least two biologically independent replicates were used for each sample. Average Ct values of actin and PP2A were used as internal reference genes for the normalization of the expression levels of the selected transcripts. Relative gene expression levels were calculated using the 2-ΔΔCt. The primer sequences used for qRT – PCR are listed in Additional file 5.
Availability of supporting data
This Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GAQP00000000.
Arora A: Chapter 4. Biochemistry of Flower Senescence. Postharvest Biol Technol Fruits, Veg Flowers. Edited by: Paliyath G, Murr DP, Handa AK, Lurie S. 2008, Iowa, USA: Wiley-Blackwell Publishing, 51-85.
Rogers HJ: From models to ornamentals: how is flower senescence regulated?. Plant Mol Biol. 2013, 82 (6): 563-574.
Vlad F, Tiainen P, Owen C, Spano T, Daher FB, Oualid F, Senol NO, Vlad D, Myllyharju J, Kalaitzis P: Characterization of two carnation petal prolyl 4 hydroxylases. Physiol Plant. 2010, 140: 199-207.
Shibuya K: Molecular Mechanisms of Petal Senescence in Ornamental Plants Kenichi Shibuya. J Japanese Soc Hortic Sci. 2012, 81: 140-149.
Arora A, Watanabe S, Ma B, Takada K, Ezura H: A novel ethylene receptor homolog gene isolated from ethylene-insensitive flowers of gladiolus (Gladiolus grandiflora hort.). Biochem Biophys Res Commun. 2006, 351: 739-744.
Guo H, Ecker JR, Bleecker ACC: Plant responses to ethylene gas are mediated by of EIN3 transcription factor. Cell. 2003, 115: 667-677.
Binder BM, Walker JM, Gagne JM, Emborg TJ, Hemmann G, Bleecker AB, Vierstra RD: The arabidopsis EIN3 binding F-Box Proteins EBF1 and EBF2 have distinct but overlapping roles in ethylene signaling. Plant Cell. 2007, 19 (2): 509-523.
Qiao H, Chang KN, Yazaki J, Ecker JR: Interplay between ethylene, ETP1 / ETP2 F-box proteins, and degradation of EIN2 triggers ethylene responses in Arabidopsis. Genes Dev. 2009, 23 (4): 512-521.
Wagstaff C, Yang TJW, Stead AD, Buchanan-Wollaston V, Roberts JA: A molecular and structural characterization of senescing Arabidopsis siliques and comparison of transcriptional profiles with senescing petals and leaves. Plant J. 2009, 57: 690-705.
Alonso M, Ecker JR, Salinas J, Novillo F: CBF2/DREB1C is a negative regulator of CBF1/DREB1B and CBF3/DREB1A expression and plays a central role in stress tolerance in Arabidopsis. Proc Natl Acad Sci U S A. 2004, 101 (11): 3985-3990.
Sakuma Y, Maruyama K, Osakabe Y, Qin F, Seki M, Shinozaki K: Functional Analysis of an Arabidopsis Transcription Factor, DREB2A, Involved in Drought-Responsive Gene Expression. Plant Cell. 2006, 18 (May): 1292-1309.
Sakuma Y, Maruyama K, Qin F, Osakabe Y, Shinozaki K: Dual function of an Arabidopsis transcription factor DREB2A in water-stress-responsive and heat-stress-responsive gene expression. 2006
Sherif S, Mila I, Bouzayen M, Jayasankar S: Molecular characterization of seven genes encoding ethylene-responsive transcriptional factors during plum fruit development and ripening. J Exp Bot. 2009, 60: 907-922.
Yin X-R, Allan AC, Chen K, Ferguson IB: Kiwifruit EIL and ERF genes involved in regulating fruit ripening. Plant Physiol. 2010, 153: 1280-1292.
Liu J, Li J, Wang H, Fu Z, Liu J, Yu Y: Identification and expression analysis of ERF transcription factor genes in petunia during flower senescence and in response to hormone treatments. J Exp Bot. 2011, 62: 825-840.
Duan M, Nan J, Liang Y, Mao P, Lu L, Li L, Wei C, Lai L, Li Y, Su X: DNA binding mechanism revealed by high resolution crystal structure of Arabidopsis thaliana WRKY1 protein. Nucleic Acids Res. 2007, 35: 1145-1154.
Ulker B, Shahid Mukhtar M, Somssich IE: The WRKY70 transcription factor of Arabidopsis influences both the plant senescence and defense signaling pathways. Planta. 2007, 226: 125-137.
Zhu X, Liu S, Meng C, Qin L: WRKY Transcription Factors in Wheat and Their Induction by Biotic and Abiotic Stress. 2013
Maleck K, Levine A, Eulgem T, Morgan A, Schmid J, Lawton KA, Dangl JL, Dietrich RA: The transcriptome of Arabidopsis thaliana during systemic acquired resistance. Nat Genet. 2000, 26: 403-410.
Dong J, Chen C, Chen Z: Expression profiles of the Arabidopsis WRKY gene superfamily during plant defense response. Plant Mol Biol. 2003, 51: 21-37.
Gregersen PL, Culetic A, Boschian L, Krupinska K: Plant senescence and crop productivity. 2013
Devoto A, Balbi V: Jasmonate signalling network in Arabidopsis thaliana: crucial regulatory nodes and new physiological scenarios. New Phytol. 2008, 177 (2): 301-318.
Miao Y, Zentgraf U: The Antagonist Function of Arabidopsis WRKY53 and ESR/ESP in Leaf Senescence Is Modulated by the Jasmonic and Salicylic Acid Equilibrium. Plant Cell. 2007, 19 (March): 819-830.
Atchley WR, Terhalle W, Dress A: Positional dependence, cliques, and predictive motifs in the bHLH protein domain. J Mol Evol. 1999, 48: 501-516.
Husbands A, Bell EM, Shuai B, Smith HMS, Springer PS: LATERAL ORGAN BOUNDARIES defines a new family of DNA-binding transcription factors and can interact with specific bHLH proteins. Nucleic Acids Res. 2007, 35: 6663-6671.
Zimmermann IM, Heim MA, Weisshaar B, Uhrig JF: Comprehensive identification of Arabidopsis thaliana MYB transcription factors interacting with R / B-like BHLH proteins. Plant J. 2004, 40: 22-34.
Morohashi K, Zhao M, Yang M, Read B, Lloyd A, Lamb R, Grotewold E: Participation of the Arabidopsis bHLH Factor GL3 in trichome initiation regulatory events. Plant Physiol. 2007, 145 (3): 736-746.
Zhao M, Morohashi K, Hatlestad G, Grotewold E, Lloyd A: The TTG1-bHLH-MYB complex controls trichome cell fate and patterning through direct targeting of regulatory loci. Development. 2008, 135 (11): 1991-1999.
Wagstaff C, Bramke I, Breeze E, Thornber S, Harrison E, Thomas B, Buchanan-wollaston V, Stead T, Rogers H: A specific group of genes respond to cold dehydration stress in cut Alstroemeria flowers whereas ambient dehydration stress accelerates developmental senescence expression patterns. J Exp Bot. 2010, 61: 2905-2921.
Master LD, Hunter DA, Steele BC, Reid MS: Identification of genes associated with perianth senescence in Daffodil (Narcissus pseudonarcissus L. “Dutch Master”). Plant Sci. 2002, 163: 13-21.
Doorn WG, Balk PA, Houwelingen AM, Hoeberichts FA, Hall RD, Vorst O, Schoot C, Wordragen MF: Gene expression during anthesis and senescence in Iris flowers. Plant Mol Biol. 2003, 53: 845-863.
Koo H, Lim K, Jung H, Park E: Anti-inflammatory evaluation of gardenia extract, geniposide and genipin. J Ethnopharmacol. 2006, 103: 496-500.
Nyon K, Choi Y, Jung H, Hyuk G, Park J, Moon S, Cho K, Kang C, Kang I, Sook M, Lee EH: International Immunopharmacology Genipin inhibits the inflammatory response of rat brain microglial cells. Int Immunopharmacol. 2010, 10: 493-499.
Qi PX, Nun A, Wickham ED: Reactions between β-Lactoglobulin and Genipin: Kinetics and Characterization of the Products. 2012
Han J, Chen S, Zhang W, Wang Y: Molecular ecology of Gardenia jasminoides authenticity. Chinese J Appl Ecol. 2006, 17: 2385-2388.
Han J, Zhang W, Cao H, Chen S, Wang Y: Genetic diversity and biogeography of the traditional Chinese medicine, Gardenia jasminoides, based on AFLP markers. Biochem Syst Ecol. 2007, 35: 138-145.
Han J-P, Chen S-L, Zhang W-S, Li X-Y: Study on genetic diversity and differentiation of Gardenia jasminoides Ellis using RAPD markers. Chinese Pharm J. 2007, 42: 1774-1778.
Criley RA, Roh MS, Kikuchi M, Manshardt RM: Acta Hortic 27th Int Hortic Congr Symp 5 - Int Symp Ornamentals 2006. A comparison of Gardenia augusta cultivars using isozymes and RAPD markers. 2008, 766: 461-468.
Lei L, Wang Y, Zhao A-M, Zhu P-L, Zhou SL: Genetic relationship of Gardenia jasminoides among plantations revealed by ISSR. Chinese Tradit Herb Drugs. 2009, 40: 117-120.
Tsanakas GF, Polidoros AN, Economou AS: Scientia Horticulturae Genetic variation in gardenia grown as pot plant in Greece. Sci Hortic (Amsterdam). 2013, 162: 213-217.
Lulin H, Xiao Y, Pei S, Wen T, Shangqin H: The first illumina-based de novo transcriptome sequencing and analysis of safflower flowers. J Climate. 2012, 7: 1-11.
Wang H, Jiang J, Chen S, Qi X, Peng H, Li P, Song A, Guan Z, Fang W, Liao Y, Chen F: Next-generation sequencing of the Chrysanthemum nankingense (Asteraceae) transcriptome permits large-scale unigene assembly and SSR marker discovery. PLoS One. 2013, 8: e62293-
Zhang X, Zhao L, Larson-rabin Z, Li D, Guo Z: De Novo Sequencing and Characterization of the Floral Transcriptome of Dendrocalamus latiflorus (Poaceae: Bambusoideae). PLoS One. 2012, 7 (8): e42082-
Guo S, Zheng Y, Joung J, Liu S, Zhang Z: Transcriptome sequencing and comparative analysis of cucumber flowers with different sex types. 2010
Kawabata S, Li Y, Miyamoto K: Scientia Horticulturae EST sequencing and microarray analysis of the floral transcriptome of Eustoma grandiflorum. Sci Hortic (Amsterdam). 2012, 144: 230-235.
Ness RW, Siol M, Barrett SCH: De novo sequence assembly and characterization of the floral transcriptome in cross- and self-fertilizing plants. BMC Genomics. 2011, 12: 298-
Müller GL, Drincovich MF, Andreo CS, Lara MV: Role of photosynthesis and analysis of key enzymes involved in primary metabolism throughout the lifespan of the tobacco flower. J Exp Bot. 2010, 61: 3675-3688.
Hoeberichts FA, Van Doorn WG, Vorst O, Hall RD: Sucrose prevents up-regulation of senescence-associated genes in carnation petals. J Climate. 2007, 58: 2873-2885.
Bey M, Stu K, Fellenberg K, Schwarz-sommer Z, Sommer H: Characterization of Antirrhinum Petal Development and Identification of Target Genes of the Class B MADS. J Climate. 2004, 16: 3197-3215.
Thompson J: The Molecular Basis for Membrane Deterioration during Senescence. 1988, San Diego: Academic, 51-83.
Tanase K, Nishitani C, Hirakawa H, Isobe S, Tabata S, Ohmiya A, Onozaki T: Transcriptome analysis of carnation (Dianthus caryophyllus L.) based on next-generation sequencing technology. BMC Genomics. 2012, 13: 292-
Chang Y-Y, Chu Y-W, Chen C-W, Leu W-M, Hsu H-F, Yang C-H: Characterization of Oncidium “Gower Ramsey” transcriptomes using 454 GS-FLX pyrosequencing and their application to the identification of genes associated with flowering time. Plant Cell Physiol. 2011, 52: 1532-1545.
Price AM, Aros Orellana DF, Salleh FM, Stevens R, Acock R, Buchanan-Wollaston V, Stead AD, Rogers HJ: A comparison of leaf and petal senescence in wallflower reveals common and distinct patterns of gene expression and physiology. Plant Physiol. 2008, 147: 1898-1912.
Garg R, Patel RK, Tyagi AK, Jain M: De novo assembly of chickpea transcriptome using short reads for gene discovery and marker identification. DNA Res. 2011, 18: 53-63.
Zahn LM, Ma X, Altman NS, Zhang Q, Wall PK, Tian D, Gibas CJ, Gharaibeh R, Leebens-Mack JH, Depamphilis CW, Ma H: Comparative transcriptomics among floral organs of the basal eudicot Eschscholzia californica as reference for floral evolutionary developmental studies. Genome Biol. 2010, 11: R101-
Zhang W, Sun Y, Timofejeva L, Chen C, Grossniklaus U, Ma H: Regulation of Arabidopsis tapetum development and function by DYSFUNCTIONAL TAPETUM1 (DYT1) encoding a putative bHLH transcription factor. Development. 2006, 133: 3085-3095.
Heisler MG, Atkinson A, Bylstra YH, Walsh R, Smyth DR: SPATULA, a gene that controls development of carpel margin tissues in Arabidopsis, encodes a bHLH protein. Development. 2001, 128: 1089-1098.
Bachem CW, Oomen RJ, Visser RG: Transcript Imaging with cDNA-AFLP: A Step-by-Step Protocol. Plant Mol Biol Report. 1998, 16: 157-173.
Grabherr M, Haas B, Yassour M, Levin J, Thompson D, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, Palma F, Birren B, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A: Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotechnol. 2011, 29: 644-652.
Iseli C, Jongeneel C, Bucher P: ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences. Proc Int Conf Intell Syst Mol Biol. 1999, Palo Alto, CA: AAAI, 138-148.
Livak KJ, Schmittgen TD: Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001, 25: 402-408.
The authors would like to thank the cooperative “Gardenia Growers Group” of Greece for providing the plant material and the Assoc. Prof. Dr. Danae Venieri of the School of Environmental Engineering (Technical University of Crete) for providing the real - time equipment for the analysis. We would also like to thank Thodhoraq Spano for excellent technical support. GFT also acknowledges the State Scholarships Foundation of Greece (IKY) for granting a doctoral scholarship and MAICh as the hosting Institution.
The authors declare that they have no competing interests.
Conceived and designed the experiments: PK and ASE. Performed the experiments: GFT. Analyzed the data: MEM. Wrote the paper: GFT, MEM, PK, ASE. All authors read and approved the final manuscript.