De novo transcriptome analysis of petal senescence in Gardenia jasminoides Ellis
© Tsanakas et al.; licensee BioMed Central Ltd. 2014
Received: 19 February 2014
Accepted: 11 June 2014
Published: 4 July 2014
Skip to main content
© Tsanakas et al.; licensee BioMed Central Ltd. 2014
Received: 19 February 2014
Accepted: 11 June 2014
Published: 4 July 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.
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.
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].
Up- and down-regulated transcripts during the progression of senescence
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.
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.
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).
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’.
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 .
Common and unique WRKY and bHLH unigenes across the stage comparisons
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.
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).
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.
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.
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.
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.
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.
This Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GAQP00000000.
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.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.