RNA-sequencing reveals early, dynamic transcriptome changes in the corollas of pollinated petunias
BMC Plant Biology volume 14, Article number: 307 (2014)
Pollination reduces flower longevity in many angiosperms by accelerating corolla senescence. This response requires hormone signaling between the floral organs and results in the degradation of macromolecules and organelles within the petals to allow for nutrient remobilization to developing seeds. To investigate early pollination-induced changes in petal gene expression, we utilized high-throughput sequencing to identify transcripts that were differentially expressed between corollas of pollinated Petunia × hybrida flowers and their unpollinated controls at 12, 18, and 24 hours after opening.
In total, close to 0.5 billion Illumina 101 bp reads were generated, de novo assembled, and annotated, resulting in an EST library of approximately 33 K genes. Over 4,700 unique, differentially expressed genes were identified using comparisons between the pollinated and unpollinated libraries followed by pairwise comparisons of pollinated libraries to unpollinated libraries from the same time point (i.e. 12-P/U, 18-P/U, and 24-P/U) in the Bioconductor R package DESeq2. Over 500 gene ontology terms were enriched. The response to auxin stimulus and response to 1-aminocyclopropane-1-carboxylic acid terms were enriched by 12 hours after pollination (hap). Using weighted gene correlation network analysis (WGCNA), three pollination-specific modules were identified. Module I had increased expression across pollinated corollas at 12, 18, and 24 h, and modules II and III had a peak of expression in pollinated corollas at 18 h. A total of 15 enriched KEGG pathways were identified. Many of the genes from these pathways were involved in metabolic processes or signaling. More than 300 differentially expressed transcription factors were identified.
Gene expression changes in corollas were detected within 12 hap, well before fertilization and corolla wilting or ethylene evolution. Significant changes in gene expression occurred at 18 hap, including the up-regulation of autophagy and down-regulation of ribosomal genes and genes involved in carbon fixation. This transcriptomic database will greatly expand the genetic resources available in petunia. Additionally, it will guide future research aimed at identifying the best targets for increasing flower longevity by delaying corolla senescence.
The longevity of individual flowers is genetically programmed to allow for efficient reproduction while limiting energy costs associated with maintaining the petals ,. In many angiosperms, pollination reduces flower longevity and initiates global gene expression changes that lead to flower senescence ,. Pollination-induced senescence of the corolla allows for nutrients to be recycled from the petals to the developing ovary ,. In petunias, ethylene biosynthesis is induced by pollination, and the application of exogenous ethylene accelerates senescence . Ethylene in wild type petunias can be measured from pollinated styles within an hour after pollination. This initial ethylene production is not sufficient to induce corolla senescence, but is followed by ethylene biosynthesis in the corolla, which then induces petal wilting ,,. In an effort to extend flower longevity, transgenic approaches have been utilized to alter ethylene perception in petunia. These experiments have created ethylene insensitive petunia flowers that last approximately twice as long as wild type flowers and do not undergo accelerated senescence after pollination ,,,.
Pollen is thought to contain a signaling factor(s) that triggers petal senescence in ethylene-sensitive species . Relatively large amounts of 1-aminocyclopropane-1-carboxylic acid (ACC) and auxin are found in petunia pollen, but experimental evidence has shown that only excessive amounts of these substances are able to increase ethylene production and accelerate flower senescence ,. Other factors such as short-chain fatty acids and changes in electrical potential may play a larger role in pollination-induced petal senescence, either by acting as a signaling factor or by increasing ethylene sensitivity ,. While pollination induces ethylene production and leads to senescence in ethylene-sensitive flowers, it remains unclear how pollination is linked to ethylene biosynthesis. Rather than blocking downstream ethylene-induced responses to delay flower senescence, inhibiting pollination signals that lead to ethylene biosynthesis may provide an alternative means of extending flower longevity.
Transcriptomic approaches, including microarrays and RNA-sequencing (RNA-seq), have been used to profile gene expression changes during flower petal development and senescence in multiple species -. A large percentage of the genes that are up-regulated during senescence encode enzymes involved in degradation and transport. The systematic degradation of proteins, nucleic acids, lipids, and cell wall components allows for the remobilization of sugars and other nutrients before the death of the petal cells . A suppressive subtractive hybridization experiment in Alstroemeria flowers showed that genes involved in cell wall synthesis, protein synthesis, metabolism, and signaling were most abundant in the petals of younger flowers, while those involved in macromolecule breakdown were highest at the later stages . Pollination-induced senescence involves similar processes and can reduce flower longevity of Ophrys (orchid) to five or six days. In orchid labella, genes involved in macromolecular breakdown, stress and defense, and nutrient remobilization are differentially expressed after pollination. Floral scent and pigment genes are down-regulated by two days after pollination .
While microarrays have been utilized to study gene expression changes in petunia ,, to our knowledge, genome-wide expression profiling using RNA-sequencing (RNA-seq) has not been performed in petunia flowers. Microarrays are able to measure gene expression changes, but are limited by the availability of Expressed Sequence Tags (ESTs). Additionally, highly expressed genes can saturate the microarrays and reduce the accuracy of gene expression data, especially for lower expressed genes. RNA-seq experiments can provide a global overview of gene expression during corolla senescence without any a priori genetic data. The recent reductions in sequencing costs have made this technology more readily accessible to researchers. RNA-seq is particularly useful for identifying genes and their isoforms, and it can measure gene expression levels that have more than an 8,000-fold difference ,.
This experiment was designed to profile early gene expression changes in petunia corollas following pollination, with the goal of identifying the signaling pathways that are involved in initiating corolla senescence. Another objective was to generate an assembled and annotated RNA-seq transcriptome for petunia corollas. Data from this experiment will provide a valuable addition to the molecular resources available for petunia. This research will guide the future selection of promising candidate genes for extending flower longevity by delaying corolla senescence.
Results and discussion
Pollen tube growth and ethylene biosynthesis of post-anthesis petunia flowers
Pollination accelerates the senescence of petunia flowers. Inducing flower senescence by pollination synchronizes the senescence program and allows for the collection of corollas that are at a very similar stage of senescence . A characterization of pollen tube growth, ethylene production, and visual senescence symptoms in Petunia × hybrida `Mitchell Diploid’ was conducted to identify the best time points for RNA-seq library construction. The goal was to identify genes and pathways involved in early senescence signaling within the corolla, so time points before fertilization, climacteric ethylene production from the petals, and visual corolla wilting were desired.
Pollinated corollas were slightly less turgid (i.e. soft to the touch) at 36 hap and were visibly wilted by 48 hap. Corollas of pollinated flowers from 0 - 24 hours after pollination (hap) were morphologically indistinguishable from each other and from unpollinated flowers of the same age. Previous studies have shown that unpollinated flowers are not wilted until around 192 h . Pollen tube growth was measured at various times after pollination. Pollen tubes maintained a relatively steady, linear growth rate and reached the end of the style after 24 hap, but before 36 hap (Figure 1A). Ethylene biosynthesis from styles and corollas was measured separately at specific times after pollination. In the initial measurements, ethylene production could be detected from pollinated styles, and ethylene peaked at 12 and 24 hap, with a slight decline at 18 hap. Ethylene production sharply declined at 36 and 48 hap (Figure 1B). In pollinated corollas, ethylene was first detectable at 18 hap, though at very low levels (2.3 nl g-1 h-1). Ethylene production peaked at 36 hap, followed by a sharp decline at 48 hap (Figure 1C). Previous studies have demonstrated that ethylene, ACC synthase, and ACC increase within the first hap, predominantly in the stigma ,; however, this initial ethylene production (within the first seven hours) is not sufficient to induce petal wilting. Pollination, therefore, requires additional factors to induce ethylene production in the corollas that leads to petal senescence .
Petunia corolla EST library construction and evaluation
Strand-specific RNA-sequencing libraries were constructed from corolla mRNA of unpollinated and pollinated flowers at 12, 18, and 24 hours after flower opening. Using the Illumina HiSeq platform, we generated a total of 488,762,314 paired-end reads that were 101 bp in length from 18 libraries. Reads per library ranged from 11,502,467 to 47,030,266, with a mean of 27,153,462 (Table 1). After preprocessing and quality trimming, the remaining 471,116,383 paired-end reads were used for de novo transcriptome assembly. We chose Trinity for de novo assembly because it has been shown to be more accurate than other programs, including Trans-ABySS and SOAPdenovo-trans ,. A total of 161,974 contigs were generated using Trinity , and they had an N50 of 2,181 bp (Figure 2A).
To evaluate the accuracy of the assembly, the contigs were compared to 404 complete Petunia × hybrida coding sequences (CDS) available in GenBank (www.ncbi.nlm.nih.gov). From the GenBank-obtained sequences, 164 (41%) were 90-100% identical to the de novo assembled contigs (Figure 2B). The ortholog-hit ratio (OHR)  was calculated using the Solanum lycopersicum ITAG2.3 protein database, and 44% of the contigs had an OHR between 0.8 and 1.2 (Figure 2C). Together, these comparisons indicate that the de novo assembly was robust and accurate.
To generate an EST library, the 162 K contigs were screened for ORFs using TransDecoder, and 37,939 contigs contained putative ORFs larger than 100 amino acids. Additionally, we added 619 contigs that had an OHR greater than 0.8 and did not share the same component identification number that was assigned by Trinity. This was done to prevent removal of contigs that had a putative S. lycopersicum ortholog. Finally, contigs of high similarity to each other (threshold of 90%) were removed using CD-HIT-EST. This threshold was selected to increase the number of uniquely mapped reads during expression analysis, and resulted in an expressed sequenced tagged (EST) library of 33,292. A total of 26,006 genes met specific annotation thresholds and were successfully annotated using Blast2GO. Our data represents the first RNA-seq generated transcriptome from petunia corollas.
Differential gene expression identifies many pollination-associated gene changes
Expression data was generated by aligning the preprocessed, quality-trimmed reads to the EST library. Approximately 84% of the reads from all libraries mapped to the EST library. We used the principle component analysis (PCA) function within the R package DESeq2  and the average linkage cluster tree analysis within the weighted gene network correlation analysis (WGCNA) R package , to screen for outlying libraries (Figure 3). PCA revealed that the libraries were segregated horizontally (PC1) based on the time of sample collection. Vertical segregation (PC2) occurred between pollinated and unpollinated samples at 18 and 24 hap. The linkage cluster tree revealed that libraries P18 r2 and P24 r3 did not group with their corresponding biological replicates. The correlation between the biological replicates of the libraries was calculated and visualized using scatterplots. All biological replicates had a strong correlation (R 2 value above 0.9) except for libraries P18 r2 and P24 r3 (Additional file 1). Based on these results, outlying libraries P18 r2 and P24 r3 were removed from further analysis. Library P18 r2 had 11.5 M reads, which is 58% lower than the average library reads. Reduced sequencing depths in RNA-seq experiments result in less reliable gene expression data, especially for low-expressed genes . The other outlying library (P24 r3) had good sequencing coverage, but did not group with the other pollinated 24 hour replicates. This may have resulted from differences in pollen load, pollen viability, or stigma damage during emasculation ,.
DESeq2 was used to identify significant pollination-associated gene changes in petunia corollas. Using normalized count data, 2,878 significant (FDR <0.05) differentially expressed genes were identified after comparing the pollinated and unpollinated treatments (P/U). Additionally, pairwise comparisons between libraries of the same time points were made (e.g. 12 hour pollinated versus 12 hour unpollinated; 12-P/U). The 12-P/U list contained 618 differentially expressed genes, 18-P/U had 2,644, and 24-P/U had 248 (Additional file 2). A total of 4,746 non-redundant (i.e. genes that were differentially expressed in more than one pairwise comparison were only counted once), pollination-associated genes were identified from these pairwise comparisons (Figure 4). These data showed that thousands of gene changes occurred well before the corollas displayed any visual symptoms of senescence (i.e. wilting) and before the pollen tubes have reached the base of the style (Figure 1A).
The total number of gene changes demonstrates the complex, dynamic, and orchestrated processes of initiating petal senescence in petunia. These findings are in line with other flower development studies. For example, RNA-seq data from developing Chimonanthus praecox (wintersweet) flowers had 2,706 differentially expressed genes between bud and open flowers and 1,466 between open and senescent flowers . More than 5,400 differentially expressed genes were identified in Rosa chinensis between open and senesced flowers . A microarray experiment in orchid (Ophrys fusca) compared pollinated and unpollinated labella and found that 8.2% of the printed cDNA clones were differentially expressed within 48 hours after pollination. These gene changes occurred long before visual cues of senescence were visualized at 5 to 6 days after pollination . Together these data demonstrate the highly dynamic nature of transcriptomic data in senescing flowers. Similarly, transcriptomic studies in leaves have identified thousands of genes that show either increased or decreased expression during leaf senescence ,.
Weighted gene correlation network analysis identified three pollination-specific modules
The differential gene expression analyses identified significant changes in thousands of genes after pollination. We hypothesized that many pollination-associated genes may be acting together in a network to regulate senescence in the corolla. Genes that form protein complexes often share similar expression patterns . To test this hypothesis, WGCNA was used to identify gene clusters (modules) that have highly correlative expression patterns. With a stringency threshold of 0.75, a total of 10 modules were identified from petunia corollas using WGCNA (Additional file 3). Three of these modules had expression patterns that were associated with pollination (i.e. changes in expression profiles appeared in only one treatment for at least one time point), and these included red (Module I), cyan (Module II), and grey60 (Module III). Heatmaps of the modules were generated to visualize the gene expression patterns over time (Figure 5). Module I had increased gene expression across all times points (12, 18 and 24 h) in corollas of pollinated flowers. This module had 1,303 genes, 75% of which also belong to the DESeq2 P/U differentially expressed gene list (Figure 6). Module II consisted of 780 contigs and was the smallest. This module's expression (i.e. eigengene) was similar at 12 and 24 h, but expression was up-regulated at 18 h in pollinated corollas. It had 348 genes (45%) in common with the 18-P/U differentially expressed genes. The largest of the modules was Module III, containing 1,359 genes. It was similar to Module II (Figure 5B, C) in expression patterns and contained 803 (59%) genes in common with the 18-P/U differentially expressed genes (Figure 6).
The WGCNA and DESeq2 analyses both identified two main expression patterns (i.e. genes that were differentially expressed in pollinated corollas and genes that were differentially expressed at 18 hap) when comparing corollas from pollinated and unpollinated flowers at the same developmental age. Pollination induced changes in gene expression that occurred prior to fertilization and ethylene biosynthesis in the corollas. After pollination, it took more than 24 h for pollen tubes to reach the bottom of the style (Figure 1A). Therefore, a signal(s) must precede fertilization to elicit the expression changes in the corolla that lead to accelerated petal senescence. Pollination signaling may involve ACC, auxin, ethylene, short-chain fatty acids, or electrical pulses ,,. Although ethylene production did not peak until 36 hap in corollas, the styles produced ethylene within the first hour after pollination and continued for 48 h. Inhibiting ethylene production or perception in the style with aminoethyoxyvinylglycine (AVG) or diazocyclopentadiene (DACP), respectively, prevents pollination-induced corolla senescence ,. These results suggest that ethylene signaling within the gynoecium is required for the corollas to respond to pollination. However, the ethylene from pollinated styles that are immediately severed from the flower, but left in the corolla, is not sufficient to accelerate senescence , suggesting that additional factors must be transmitted to the corolla to induce senescence. Wounding also results in elevated ethylene production from petunia stigmas, and at 17 hours after the stigma wounding, petal wilting can no longer be delayed by removing the damaged stigmas . This suggests that the necessary signals for stigma-induced, flower senescence are in place within the first 17 hours after stigma wounding. Short-chain fatty acids that are produced in the gynoecium and transported to the corolla within 12 h of pollination have been shown to increase ethylene sensitivity in corollas, and this may be a component of the pollination signaling ,.
Validation of RNA-seq data by quantitative PCR
To confirm the gene expression patterns identified by the RNA-seq data, the transcript levels of five genes were examined by quantitative PCR (Figure 7). Three of the genes (comp31514_c0_seq2, comp39985_c0_seq4, and comp18014_c0_seq1) were from Module III (grey60), which was characterized by higher expression at P18 compared to U18. Quantitative PCR analysis confirmed large differences in transcript abundance between P18 and U18, with much smaller changes between P12 and U12 and P24 and U24. Two additional genes that were identified as differentially regulated between pollinated and unpollinated libraries (P/U) by DESeq2 analysis (comp40361_c0_seq2 and comp47181_c0_seq6) also showed very similar patterns of transcript abundance as determined by RNA-seq and qPCR. All of the gene expression patterns were confirmed to be consistent with the RNA-seq data.
Enriched GO terms suggest involvement of plant hormones within 12 hap
To identify the biological relevance of the pollination-associated gene changes, gene ontology (GO) was used to determine the biological processes, cellular components, and molecular functions of the differentially expressed genes  (Additional file 4). At 12 hap, 35 enriched GO terms were identified (FDR <0.05). Many of these terms involve plant hormones like abscisic acid (ABA), auxin, jasmonic acid (JA), and salicylic acid (SA). Of note are the response to auxin stimulus and response to 1-aminocyclopropane-1-carboxylic acid (ACC) GO terms. Both auxin and ACC are found in relatively high concentrations in pollen , and the corolla may be responding to hormonal signals that are transmitted through the gynoecium. At 18 hours after pollination, 154 enriched GO terms were identified including the ethylene signaling pathway. This coincided with the initiation of ethylene production from the corollas. Three of the molecular function GO terms involve autophagy. Autophagy is a catabolic process that involves transporting cellular components to the vacuole for further degradation and nutrient recycling . No enriched terms were identified at 24 hap, but 368 enriched terms were identified when comparing pollinated to unpollinated (P/U) corollas at any time (12, 18, and 24 h). Enriched terms consisted of sucrose metabolic process, response to chitin, and response to wounding. The number of GO terms (557 in total of which 508 were unique) reflects the breadth of changes that occur between 12 and 24 hap in corolla tissue (Additional file 4).
KEGG enrichment identifies pollination responsive pathways in the corolla
To identify the molecular pathways associated with pollination-induced corolla senescence, the significant DESeq2 and WGCNA genes were searched against A. thaliana proteins using BLASTx . Top BLASTx hits were considered as the putative A. thaliana orthologs (Additional file 2). These hits were mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. A total of 15 unique, enriched KEGGs were identified (Table 2 and Additional file 5). The KEGG pathways provided insight into potential biological pathways that function in the corollas of pollinated flowers. For example, eight of the KEGGs were involved in metabolism, including carbohydrate, energy, and lipid metabolism as well as the metabolism of terpenoids and polyketides. Other KEGGs were categorized under transport and catabolism, translation, signal transduction, and environmental adaptation.
Four enriched KEGG pathways were identified in pollinated corollas
Four unique, enriched KEGG pathways were identified from the P/U genetic changes identified in DESeq2 and the WGCNA Module I (red). They included Plant-pathogen interactions, Starch and sucrose metabolism, Pentose and glucuronate interconversions, and Plant hormone signal transduction (Table 2). The genes within these KEGG pathways are associated with pollination and may contain key signaling components and molecular events that lead to flower senescence.
The Plant-pathogen interaction pathway was enriched in both P/U and in Module I. Genes encoding enzymes that are putatively involved in defense have been reported to be up-regulated during the senescence of many different flowers ,. In our analysis, 35 P/U genes mapped to this pathway, and 20 genes mapped from Module I (Figure 8 and Additional file 5). The majority (80%) of these genes are predicted to interact with Ca2+ in this pathway, and some examples include calcium-dependent protein kinases, putative calcium binding proteins, and calmodulin-like proteins. While the role of calcium in the corollas of pollinated petunias has not been delineated, many studies have been performed to understand the role of Ca2+ signaling within the pollinated gynoecium. Changes in electrical potential have been observed , and calcium signaling is integral to pollen germination, pollen tube growth, and fertilization -. Our data suggest that Ca2+ signaling is continued from the style to the corolla and may be important for relaying pollination signals to the petals to initiate corolla senescence. A CBL-interacting kinase (CIPK) is up-regulated in an ethylene-dependent manner early in the senescence of carnation flowers . CIPK regulates phosphorylation cascades that transmit Ca2+ signals, and it was hypothesized from these studies that calcium signaling was involved in carnation petal senescence.
Pollination and fungal infection share striking similarities in terms of biological responses, and both processes result in cell death ,. X-ray microanalysis revealed that both pollen tubes and fungal hyphae penetration result in the accumulation of Ca2+ at the interaction sites . Two well-known microbe-associated molecular pattern (MAMP) LRR receptor-like serine-threonine protein kinases, flagellin insensitive 2 (FLS2) and EF-Tu receptor (EFR), were both up-regulated following pollination (Figure 8). Activation of these receptors results in changes in ion flux, reactive oxygen species formation, MAP kinase activation, and ethylene production . It has been hypothesized that pathogen-related proteins are up-regulated during senescence to protect the senescing tissue from pathogenic attack , but petunia pollen tubes may contain an elicitor-like motif that activates FLS2 and EFR. Altering or eliminating these elicitors from pollen may prevent or delay pollination-induced senescence. Alternatively, increased expression of these genes may be a result of elevated ethylene levels. EIN3 and EIL have been shown to activate transcription of FLS2 in Arabidopsis .
The Starch and sucrose metabolism pathway involves the catabolism of carbohydrates. The P/U list had 30 genes map to this KEGG pathway and Module I had 19 (Table 2 and Additional file 5). Many of these genes are involved in the conversion of UDP-D-galacturonate to D-galacturonate, which interacts with ascorbate metabolism. There are also many pectinesterase genes involved in the catabolism of pectin (Additional file 5). Soluble carbohydrates move from senescing to non-senescing flowers in gladiolus . Sugars, particularly sucrose, increase in the phloem of Ipomoea and Hemerocallis (daylily) petals as the flowers open, mature, and senesce ,. Labeling studies in carnations demonstrate that sucrose moves in the phloem from the petals to the gynoecium during senescence, and that this remobilization is accelerated by ethylene treatment . The enrichment of this KEGG pathway suggests that a similar process involving the movement of carbohydrates to sinks, like the developing ovules, may also occur following pollination in petunia. Sucrose has profound effects on extending flower longevity, and has been implicated in the stability of EIN3 in Arabidopsis . The application of sucrose to cut carnation flowers delays petal senescence and the up-regulation of genes involved in ethylene signaling . The competition for carbohydrates also regulates the timing of senescence in ethylene-insensitive flowers like lilies (Lilium), where flower senescence is observed once the carbohydrate content of the tepals is reduced by ~50% .
The Pentose and glucuronate interconversions KEGG pathway was enriched in Module I, but not in the P/U gene list. A total of 12 mapped genes were found within Module I, nine of which involve pectin degradation (Table 2). This pathway contained five pectinesterase proteins, a polygalacturonase, and a UDP-glucose 6-dehydrogenase that overlap with the Starch and sucrose metabolism pathways. The other four Pentose and glucuronate interconversion-associated proteins are putative pectin lyase proteins. These enzymes are involved in cell wall loosening and have been shown to increase free Ca2+ levels as the calcium-cross-linked bridges are lysed (Additional file 5) ,,. Galactose loss is the main feature of cell wall changes during the senescence of petunia, Sandersonia and carnation flowers -.
Plant hormones are an integral part of flower senescence , and the Plant hormone signal transduction KEGG pathway was enriched in the P/U gene list. Ethylene-sensitive crops, like petunia, produce ethylene after pollination , and the application of ethylene synthesis and perception inhibitors delays flower senescence when applied to the whole flower or to the pollinated gynoecium . Transgenic petunias containing the mutant allele etr1-1  or a down-regulated EIN2  gene do not exhibit accelerated senescence after pollination, proving that ethylene is required for post-pollination signaling between the gynoecium and the corolla . Nearly all of the 37 genes from the P/U gene list that mapped to this pathway were up-regulated, and they involved members of every major plant hormone pathway (Figure 9). Abscisic acid (ABA), ethylene, and jasmonic acid (JA) lead to genetic changes that promote senescence, while exogenous applications of cytokinin and gibberellin slow senescence ,-. The complex interplay of these plant hormones in petunia senescence is not fully understood, and protein-protein networks can provide preliminary information about potential targets for further analysis. We utilized STRING (string-db.org) to view the plant hormone protein interactions within a network using TAIR (The Arabidopsis Information Recourse) codes obtained via BLASTx restricted to A. thaliana. This network suggests that there are direct or indirect interactions between the plant hormones ethylene and JA as well as ABA and auxin (Figure 10). Of interest is the transcription factor APETALA 2 (AP2), which shares edges with ABA, auxin, ethylene, and salicylic acid (SA). This transcription factor belongs to the AP2/ERF family and is important for flower and seed development . In Brassica napus, BnAP2 RNAi lines have delayed flower abscission and senescence . In tomato, AP2 RNAi lines have altered levels of ethylene biosynthesis and differences in the timing of fruit ripening .
Eleven KEGGs are enriched at 18 hap
Large gene changes were observed specifically at 18 hap, and 11 enriched KEGG pathways were identified (see KEGGs designated with Module I and 18-P/U in Table 2). Most of the genes within the alpha-Linolenic acid metabolism, Endocytosis, Limonene and pinene degradation, Peroxisome, and Regulation of autophagy KEGG pathways were up-regulated following pollination, while the Carbon fixation in photosynthetic organisms, Ribosome, and Ribosome biosynthesis in eukaryotes KEGG pathways were down-regulated (Table 2, Additional file 5, and Additional file 6).
The Regulation of autophagy was one of the most significantly (FDR = 5.05 × 10-6) enriched KEGG pathways, with 10 of 15 genes mapped. Autophagy is an intracellular degradation process where cytoplasmic constituents are transported to the vacuole for degradation so that nutrients can be remobilized ,. Genes mapped throughout this pathway, including processes involving autophagy induction, vesicle nucleation, and vesicle expansion and completion (Figure 11). A previous high-resolution temporal profiling of Arabidopsis leaf senescence also identified 15 autophagy genes that were up-regulated during senescence . Shibuya et al.  reported that the transcript abundance of PhATG8 in petunia increases with ethylene exposure, and this coincides with increased nitrogen content within the ovary. The nitrogen content of `Mitchell Diploid' petunia corollas decreases by greater than 60% during pollination-induced senescence . Recently, pulse/chase experiments with 15 N have shown that nitrogen remobilization is reduced in atg mutants, and that this decreases biomass and yield ,. Autophagy also clearly has a role in longevity, because atg mutants consistently display early leaf senescence. Evidence suggests that autophagy has both pro-survival and pro-death roles during plant development, but it is unclear how this dual function is regulated .
The petunia autophagy genes APG5, APG7, APG8H, ATG1C, ATG13, ATG6, ATG8C, and ATG8F are up-regulated at 18 h in corollas of pollinated flowers (Additional file 5). The 18-hour time point was collected after approximately 5 hours of darkness. This suggests that many pollination-induced autophagy genes may be regulated by darkness or may be functioning during the night. Rubisco degradation via autophagy occurs in early stages of dark-induced senescence . During this process, Rubisco is remobilized to the vacuole, and a decrease in chlorophyll can be measured after just one day of darkness. However, in Arabidopsis atg5 mutants, Rubisco is not remobilized in darkened leaves ,. Changes in the expression of autophagy genes have also been reported during starch degradation in darkened Arabidopsis leaves . Similarly, the remobilization of Rubisco from chloroplasts in the petunia corolla, which are primarily located in the tube, may be occurring via autophagy. Monodansylcadaverine (MDC) staining has been used to visualize the accumulation of acidic bodies during the pollination-induced senescence of petunia petals, but MDC staining is not specific to autophagosomes . While all of these studies have provided compelling evidence for the involvement of autophagy in corolla senescence, additional morphological studies are needed to confirm the accumulation of autophagosomes in senescing petunia corollas .
Five transcription factor families have more than 20 members in pollinated corollas
Among the 4,746 differentially expressed genes we identified 301 putative transcription factors from 42 different families. More than 20 members from each of the following transcription factor families were identified: ERFs, NAC, bZIP, HD-Zip and WRKY (Figure 12). These transcription factor families have been implicated in senescence, abiotic stress responses, and plant hormone studies -. For example, in the bZIP transcription factor superfamily there is a subset of Arabidopsis transcription factors, termed S-type, that are master transcriptional and translational regulators of enzymes involved in amino acid catabolism under sucrose starvation and senescence-induced nutrient translocation ,. In soybean (Glycine max) leaves, GmbZIP53A is up-regulated during sucrose starvation and may indirectly promote autophagy under those conditions . As stated previously, sucrose has been shown to extend flower longevity and prevent ethylene production. In carnation, sucrose down regulates a key ethylene transcription factor DcEIL3, which is necessary for the initiation of the ethylene response . Transcriptome changes in corollas from transgenic petunias with an inducible etr1-1 were identified using microarrays . This etr1-1 transgene allows for ethylene insensitivity to be controlled by applying dexamethasone (DEX). A comparison between etr1-1-induced corollas and non-induced corollas revealed that B-box zinc finger, bHLH DNA-binding, homeodomain-like (HD), MADS-box, MYB domain, and NAC domain proteins are down regulated in ethylene insensitive corollas. These findings reflect our results, in that many of these proteins were found to be transcriptionally up-regulated after pollination (Figure 12). Altering the expression of transcription factors can have profound effects on flower longevity. For example, the down regulation of PhHD-Zip delays corolla senescence in both pollinated and unpollinated petunia flowers . More research is needed to identify the functional role and interactions of these transcription factors for further manipulation of flower longevity.
Uncharacterized genes may play an integral role in pollination and future research
The KEGG enrichment analysis provided a wealth of biological relevance through the identification of 15 uniquely enriched pathways. This provided a meaningful framework for the specific biological activities that are involved in pollination-induced corolla senescence in petunias. These enriched KEGG pathways still only represent a minority of the differentially expressed genes. The remaining differentially expressed genes likely also hold important biological relevance, particularly for those genes and pathways that might be unique to petunia. This analysis demonstrates the power of next generation sequencing to capture a global overview of thousands of gene expression changes in a single experiment. Understanding the relevance of these genes is currently the rate-limiting step. This technology provides fundamental data upon which more hypothesis-driven experiments can be organized and conducted.
Pollination induces many hormonal, physiological, and molecular changes in petunia corollas that lead to senescence. Gene expression in the corollas was already altered by 12 hap, and 618 differentially expressed genes were identified. These changes occurred well before fertilization, ethylene biosynthesis from the corolla, and petal wilting. At 18 hap, large changes in gene expression were measured and an additional 2,137 genes were identified as being differentially expressed. The enriched GO term analysis suggested that at 12 hap, the corollas were responding to auxin and ACC, which are found in high abundance in pollen. KEGG enrichment identified 15 pathways, 11 of which were involved in metabolic pathways and autophagy regulation. The sequence data from this experiment will make a valuable contribution to the genomic resources available in petunia and will enable researchers to identify the genes involved in regulating flower senescence. While senescence studies have demonstrated that the initiation and timely progression of senescence requires transcription, senescence is also controlled post-transcriptionally. Previous studies in petunia have shown that genes expression patterns do not always correlate with protein expression . Combined genomic, proteomic, and metabolomic approaches will be required to gain a comprehensive understanding of petal senescence.
The plants used in this study were Petunia × hybrida (Hook.) Vilm. `Mitchell Diploid' , a doubled haploid derived from a P. axillaris/P.hybrida `Rose of Heaven' hybrid . Seeds were originally obtained from Dr. David Clark (University of Florida). Petunia seeds were sown in plug trays on soil-less media (Pro-mix BX, Premier Horticulture, Quebec, Canada) and grown under fluorescent, full-spectrum lights. After four weeks, plants were transplanted into 16-cm pots and moved to a greenhouse with temperatures set at 24/16°C (day/night) and a 13 h photoperiod. Supplemental lighting was supplied by high pressure sodium and metal halide lights (GLX/GLS e-systems GROW lights, PARSource, Petaluma, CA, USA) to maintain light levels above 300 μmol M-2 s-1.
Pollen tube growth measurements
To prevent self-pollination, anthers were removed 1 d before flower opening. On the following day, emasculated flowers were pollinated between 8:00 and 8:30 AM. Four styles were collected at 0, 12, 18, 24, and 36 hours after pollination (hap) and submerged in a 3:1 ratio of ethanol and acetic acid to fix the tissue overnight. They were then rinsed with 1 M potassium phosphate buffer (pH 7.0) followed by submersion in 1 N sodium hydroxide for 24 h. Finally, the styles were triple rinsed in sterile dH2O and stained with 0.1% aniline blue overnight. Styles were fixed on glass slides and visualized under an inverted epifluorescence Leica DM IRB microscope (Wetzlar, Germany) equipped with a Q Imaging Retiga 2000 cooled digital camera (Burnaby, BC, Canada). The lengths of the pollen tubes were measured using ImageJ .
Three biological replicates of two pollinated and unpollinated flowers were collected and photographed at 0, 12, 18, 24, 36 and 48 h. The corollas and styles from those flowers were collected and sealed in 22 ML and 7 ML glass vials, respectively. After a 30 Minute incubation period, 1 ML of the headspace was withdrawn from each vial through a rubber septum in the lid. The samples were injected into a gas chromatograph equipped with a flame ionization detector and separated on a stainless steel column packed with HayeSep R (Varian 3800, Agilent, Santa Clara, CA, USA). The flow rate of the carrier gas (He) was 30 ML Min-1.
RNA extraction and library preparation
Flowers were emasculated as described previously. Four unpollinated and four pollinated flowers were harvested at 12, 18, and 24 h after flower opening. Three biological replicates were collected for each treatment-time combination. Corollas were detached from the receptacle (which removed the ovary and style), filaments were removed, and corollas were rinsed with sterile dH2O to remove any pollen contamination. Total RNA was extracted from the corollas using Trizol reagent (Invitrogen, Carlsbad, CA, USA) followed by an additional purification step using mini spin columns (Qiagen, Valencia, CA, USA). The quality of the RNA was determined using an Agilent 2100 Bioanalyzer RNA 6000 Nano kit (Agilent, Santa Clara, CA, USA) and it was quantified using a Qubit 2.0 fluorometer RNA Assay Kit (Invitrogen Inc. Carlsbad, CA, USA). A total of 5 μg of RNA was used to create each strand-specific RNA-seq library. Eighteen libraries (3 time points × 2 treatments × 3 biological replications) with six unique barcodes were prepared following the protocol of Zhong et al. , including the modification using the universal adaptor system. The libraries were sequenced at the Genomics Resources Core Facility at Weill Cornell Medical College (New York, NY, USA). Paired-end sequences (101 bp) were generated using three lanes of an Illumina HiSeq2000 flow cell (Ilumina Inc. San Diego, CA, USA). Individual biological replicates for each library were run in separate lanes on the flow cell.
Sequence quality assessment and de novo assembly
Sequence qualities were assessed before and after trimming using FastQC version 0.10.1 (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc). Reads with a Phred quality score less than 20, and sequences shorter than 40 bp, were removed using trim-fastq.pl version 1.2.2 . This resulted in two files that contained proper paired-end sequences and one file that contained sequences that lost the mate due to the preprocessing. Adaptors, barcodes, polyA, and polyT ends were trimmed using cutadapt version 1.2.1 . After trimming, paired-end sequences were normalized to a maximum depth of 1,500 and assembled using Trinity r2012-10-05 . To create an EST database for further analysis, the resulting contigs were screened for putative open reading frames (ORF) using the TransDecoder utility from Trinity. Additionally, contigs that had both an ortholog hit ratio  of more than 80% to the Solanum lycopersicum ITAG2.3 protein database (using BLASTx) and a unique component (comp#) and subcomponent (c#) were added to the EST library. Finally, CD-HIT-EST  was used to remove contigs that had 90% or greater sequence identity to each other.
EST library annotation
The EST library was annotated using Blast2GO version 2.7.0 . The translated Basic Local Alignment Search Tool (BLASTx)  was used to obtain top hits from the SwissProt database  for each contig using a minimum E-value threshold of 1.0 × 10-3. The remaining contigs with no BLASTx hits were aligned against the non-redundant (NR) database from National Center for Biotechnology Information (NCBI). Following the BLASTx step, finalized annotations for each gene were filtered with an E-value of 1.0 × 10-6, gene ontology (GO) terms were added, and conserved domains were identified using the InterPro Scan tool  for each contig. To obtain A. thaliana-specific annotations, the Arabidopsis proteome database was downloaded from UniProt and BLASTx was performed locally.
Expression analysis from read mapping
Burrows-Wheeler Aligner (BWA) version 0.7.5a-r405  was used to align the unprocessed Illumina reads to the EST library using the default alignment stringency. Paired-end and single reads that resulted from the pre-processing step as mentioned above were used to calculate the expression profile of each contig within a library. Sam2counts.py (https://github.com/vsbuffalo/sam2counts/blob/master/sam2counts.py) generated count tables of the reads that aligned to the EST library. Only uniquely-mapped reads were used for differential gene expression analysis. The R package DESeq2 version 1.4.1  was used to determine the significant differentially expressed genes. In this package, principle component analysis (PCA) was used to screen for outliers among the libraries. A base mean threshold of ten was set to eliminate contigs with few counts, since contigs with very low reads typically have inaccurate expression patterns due to sampling error. Comparisons of all pollinated and unpollinated corollas (P/U), 12-h pollinated with 12-h unpollinated (12-P/U), 18-h pollinated and unpollinated (18-P/U), and 24-h pollinated and unpollinated (24-P/U) were made. An adjusted p-value (using the Benjamini & Hockberg adjustment) of 0.05 was used as the statistical cutoff for differentially expressed genes. A Venn diagram was used to visualize overlapping genes between comparisons .
Quantitative PCR validation of gene expression patterns
The expression patterns of five genes (comp31514_c0_seq2, comp39985_c0_seq4, comp18014_c0_seq1, comp40361_c0_seq2, and comp47181_c0_seq6) were analyzed using quantitative real-time PCR (qPCR). RNA from four biological replicates of each treatment and time point were included in the qPCR analysis. cDNA was synthesized from 2 μg total RNA using Omniscript RT Kit (Qiagen, Valencia, CA). Primers were designed to amplify the specific transcripts using IDT Primer Quest (Additional file 7). Quantitative PCR was performed in a 20 μl reaction volume using iQ SYBR Green Supermix (Bio-Rad, Hercules, CA) with 1 μl cDNA as template as described previously . Each reaction was performed in triplicate. Amplification specificity was determined by melt curve analysis. Amplification efficiencies of the target genes and reference genes were similar. Fold change for each target gene, normalized to PhACTIN, was calculated relative to expression in the U12 sample using the 2-ΔΔ Cq method.
Weighted gene correlation network analysis (WGCNA)
The R package WGCNA version 1.36 , was used to identify modules within the data set and to create dendrograms and heatmaps. A soft threshold value, power of 16, was used to transform the adjacency matrix to meet the scale-free topology criteria for optimal clustering. The libraries were clustered to identify outlier libraries using an average linkage hierarchical cluster tree based on Euclidean distance. Modules were grouped using a stringency threshold of 0.75. The code for the WGCNA analysis is available at the GetHub repository (https://github.com/wijerasa/WGCNA_Analysis). The pollination-specific modules were identified as Module I (red), Module II (cyan), and Module III (grey60) (Figure 5 and Additional file 3).
GO and KEGG enrichment
An enriched GO term analysis was conducted using a Fisher's Exact Test on the differentially expressed genes in Blast2GO. This test includes a correction for multiple testing  to reduce the false discovery rate (FDR). GO terms with a Term Filter Value of above 0.05 were excluded. TAIR codes for the EST library were obtained from a BLASTx that was restricted to A. thaliana. DESeq2 and WGCNA module contigs were mapped directly in KEGG mapper (http://www.genome.jp/kegg/mapper.html). The hypergeometric function in R was used to test for enriched pathways, and P-values were adjusted using an FDR correction.
Protein-protein interactions visualized in STRING
To visualize the plant hormone protein interactions, TAIR codes for the P/U genes were uploaded into STRING to identify the genes belonging to the Plant hormone signal transduction KEGG pathway and their nearest interacting partner (stringency of interactions set at high confidence level of 0.70) . The genes and their interacting partners were then uploaded into STRING and a network image was generated. The confidence view, which displays edges as blue lines, was selected and the image was exported. Colors of the circles were altered in Photoshop CS6 (Adobe, San Jose, CA, USA).
Transcription factor analysis
Transcription factors within the 4,746 unique DESeq2 differentially expressed genes were identified by matching TAIR codes to the A. thaliana transcription factor database (Plant TFDB v3.0 ).
Availability of supporting data
RNA-seq data were deposited with the Sequence Read Archive (SRA) database at NCBI (BioProject ID: PRJNA259884). Code for the WGCNA analysis can be accessed at the GitHub repository (https://github.com/wijerasa/WGCNA_Analysis).
SB participated in experimental design, collected tissue samples, participated in library construction and data analysis, and drafted the manuscript. SW participated in all bioinformatics analyses. AW participated in experimental design and bioinformatics analyses. LC conducted the qPCR validation experiments and edited the manuscript. TM participated in experimental design and data analyses. MJ participated in experimental design, data interpretation, manuscript writing and editing. All authors read and approved the final manuscript.
Ashman TL, Schoen DJ: How long should flowers live?. Nature. 1994, 371 (6500): 788-791.
Chapin L, Jones M: Nutrient remobilization during pollination-induced corolla senescence in petunia. Acta Hortic. 2007, 755: 181-190.
Stead A: Pollination-induced flower senescence: a review. Plant Growth Regul. 1992, 11 (1): 13-20.
Jones ML: Ethylene signaling is required for pollination-accelerated corolla senescence in petunias. Plant Sci. 2008, 175 (1-2): 190-196.
Jones ML: Mineral nutrient remobilization during corolla senescence in ethylene-sensitive and -insensitive flowers. AoB Plants. 2013, 5: plt023-
Wilkinson JQ, Lanahan MB, Clark DG, Bleecker AB, Chang C, Meyerowitz EM, Klee HJ: A dominant mutant receptor from Arabidopsis confers ethylene insensitivity in heterologous plants. Nat Biotechnol. 1997, 15 (5): 444-447.
Jones M, Stead A, Clark D: Petunia flower senescence. Petunia. Edited by: Gerats T, Strommer J. Springer, New York; 2009:301-324.
Hoekstra FA, Weges R: Lack of control by early pistillate ethylene of the accelerated wilting of Petunia hybrida flowers. Plant Physiol. 1986, 80 (2): 403-408.
Shibuya K, Barry KG, Ciardi JA, Loucas HM, Underwood BA, Nourizadeh S, Ecker JR, Klee HJ, Clark DG: The central role of PhEIN2 in ethylene responses throughout plant development in petunia. Plant Physiol. 2004, 136 (2): 2900-2912.
Langston BJ, Bai S, Jones ML: Increases in DNA fragmentation and induction of a senescence-specific nuclease are delayed during corolla senescence in ethylene-insensitive (etr1-1) transgenic petunias. J Exp Bot. 2005, 56 (409): 15-23.
Gilissen LJW, Hoekstra FA: Pollination-induced corolla wilting in Petunia hybrida rapid transfer through the style of a wilting-inducing substance. Plant Physiol. 1984, 75 (2): 496-498.
Singh A, Evensen KB, Kao TH: Ethylene synthesis and floral senescence following compatible and incompatible pollinations in Petunia inflata . Plant Physiol. 1992, 99 (1): 38-45.
Whitehead CS, Halevy AH: Ethylene sensitivity: the role of short-chain saturated fatty-acids in pollination-induced senescence of Petunia hybrida flowers. Plant Growth Regul. 1989, 8 (1): 41-54.
Liu D, Sui S, Ma J, Li Z, Guo Y, Luo D, Yang J, Li M: Transcriptomic analysis of flower development in wintersweet (Chimonanthus praecox). PLoS One. 2014, 9 (1): e86976-
Hoeberichts FA, van Doorn WG, Vorst O, Hall RD, van Wordragen MF: Sucrose prevents up-regulation of senescence-associated genes in carnation petals. J Exp Bot. 2007, 58 (11): 2873-2885.
Wang Y, Huang H, Ma Y, Fu J, Wang L, Dai S: Construction and de novo characterization of a transcriptome of Chrysanthemum lavandulifolium: analysis of gene expression patterns in floral bud emergence. Plant Cell Tiss Org. 2014, 116 (3): 297-309.
Wang H, Stier G, Lin J, Liu G, Zhang Z, Chang Y, Reid MS, Jiang C: Transcriptome changes associated with delayed flower senescence on transgenic petunia by inducing expression of etr1-1, a mutant ethylene receptor. PLoS One. 2013, 8 (7): e65800-
Jones ML: Changes in gene expression during senescence. Plant Cell Death Processes. Edited by: Nooden L. Elsevier Science, San Diego, California; 2004:51-72.
Monteiro F, Sebastiana M, Figueiredo A, Sousa L, Cotrim HC, Pais MS: Labellum transcriptome reveals alkene biosynthetic genes involved in orchid sexual deception and pollination-induced senescence. Funct Integr Genomics. 2012, 12 (4): 693-703.
Breeze E, Wagstaff C, Harrison E, Bramke I, Rogers H, Stead A, Thomas B, Buchanan-Wollaston V: Gene expression patterns to define stages of post-harvest senescence in Alstroemeria petals. Plant Biotechnol J. 2004, 2 (2): 155-168.
Price AM, Orellana DFA, 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 (4): 1898-1912.
van Doorn WG, Balk PA, van Houwelingen AM, Hoeberichts FA, Hall RD, Vorst O, van der Schoot C, van Wordragen MF: Gene expression during anthesis and senescence in Iris flowers. Plant Mol Biol. 2003, 53 (6): 845-863.
van Doorn WG, Woltering EJ: Physiology and molecular biology of petal senescence. J Exp Bot. 2008, 59 (3): 453-480.
Wang Z, Gerstein M, Snyder M: RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009, 10 (1): 57-63.
Wang Y, Ghaffari N, Johnson CD, Braga-Neto UM, Wang H, Chen R, Zhou H: Evaluation of the coverage and depth of transcriptome by RNA-seq in chickens. BMC Bioinformatics. 2011, 12: S5-
Bai S, Willard B, Chapin LJ, Kinter MT, Francis DM, Stead AD, Jones ML: Proteomic analysis of pollination-induced corolla senescence in petunia. J Exp Bot. 2010, 61 (4): 1089-1109.
Jones ML, Chaffin GS, Eason JR, Clark DG: Ethylene-sensitivity regulates proteolytic activity and cysteine protease gene expression in petunia corollas. J Exp Bot. 2005, 56 (420): 2733-2744.
Kovaleva LV, Timofeeva GV, Zakharova EV, Voronkov AS, Rakitin VY: Ethylene synthesis in petunia stigma tissues governs the growth of pollen tubes in progamic phase of fertilization. Russ J Plant Physl. 2011, 58 (3): 402-408.
Villarino GH, Bombarely A, Giovannoni JJ, Scanlon MJ, Mattson NS: Transcriptomic analysis of Petunia hybrida in response to salt stress using high throughput RNA sequencing. PLoS One. 2014, 9 (4): e94651-
Soetaert SSA, Van Neste CMF, Vandewoestyne ML, Head SR, Goossens A, Van Nieuwerburgh FCW, Deforce DLD: Differential transcriptome analysis of glandular and filamentous trichomes in Artemisia annua . BMC Plant Biol. 2013, 13: UNSP 220-
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, 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 (7): 644-652.
O'Neil ST, Dzurisin JDK, Carmichael RD, Lobo NF, Emrich SJ, Hellmann JJ: Population-level transcriptome sequencing of nonmodel organisms Erynnis propertius and Papilio zelicaon . BMC Genomics. 2010, 11: 310-
Anders S, Huber W: Differential expression analysis for sequence count data. Genome Biol. 2010, 11: R106-
Langfelder P, Horvath S: Fast R functions for robust correlations and hierarchical clustering. J Stat Softw. 2012, 46 (11): 1-17.
Langfelder P, Horvath S: WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008, 9: 559-
Woltering EJ, Vrije T, Harren F, Hoekstra FA: Pollination and stigma wounding: same response, different signal?. J Exp Bot. 1997, 48 (5): 1027-1033.
Gilissen LJW: Style-controlled wilting of flower. Planta. 1977, 133 (3): 275-280.
Yan H, Zhang H, Chen M, Jian H, Baudino S, Caissard J, Bendahmane M, Li S, Zhang T, Zhou N, Qiu X, Wang Q, Tang K: Transcriptome and gene expression analysis during flower blooming in Rosa chinensis `Pallida'. Gene. 2014, 540 (1): 96-103.
Breeze E, Harrison E, McHattie S, Hughes L, Hickman R, Hill C, Kiddle S, Kim Y, Penfold CA, Jenkins D, Zhang C, Morris K, Jenner C, Jackson S, Thomas B, Tabrett A, Legaie R, Moore JD, Wild DL, Ott S, Rand D, Beynon J, Denby K, Mead A, Buchanan-Wollaston V: High-resolution temporal profiling of transcripts during Arabidopsis leaf senescence reveals a distinct chronology of processes and regulation. Plant Cell. 2011, 23 (3): 873-894.
Guo Y, Gan S: Convergence and divergence in gene expression profiles induced by leaf senescence and 27 senescence-promoting hormonal, pathological and environmental stress treatments. Plant Cell Environ. 2012, 35 (3): 644-655.
Jansen R, Greenbaum D, Gerstein M: Relating whole-genome expression data with protein-protein interactions. Genome Res. 2002, 12 (1): 37-46.
Verlinden S: Flower senescence. The Molecular Biology and Biotechnology of Flowering. Edited by: Jordan B. 2006, CABI Publishing, Cambridge, MA, 404-
Jones ML, Woodson WR: Pollination-induced ethylene in carnation: role of stylar ethylene in corolla senescence. Plant Physiol. 1997, 115 (1): 205-212.
Lovell PJ, Lovell PH, Nichols R: The control of flower senescence in petunia (Petunia hybrida). Ann Bot. 1987, 60 (1): 49-59.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G: Gene Ontology: tool for the unification of biology. Nat Genet. 2000, 25 (1): 25-29.
Avila-Ospina L, Moison M, Yoshimoto K, Masclaux-Daubresse C: Autophagy, plant senescence, and nutrient recycling. J Exp Bot. 2014, 65 (14): 3799-3811.
Altschul SF, Madden TL, Schaffer AA, Zhang JH, Zhang Z, Miller W, Lipman DJ: Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997, 25 (17): 3389-3402.
Fromm J, Hajirezaei M, Wilke I: The biochemical response of electrical signaling in the reproductive system of Hibiscus plants. Plant Physiol. 1995, 109 (2): 375-384.
Wędzony M, Filek M: Changes of electric potential in pistils of Petunia hybrida Hort. and Brassica napus L. during pollination. Acta Physiol Plant. 1998, 20 (3): 291-297.
Dumas C, Gaude T: Fertilization in plants: is calcium a key player?. Semin Cell Dev Biol. 2006, 17 (2): 244-253.
Lenartowska M, Bednarska E, Butowt R: Ca2+ in the pistil of Petunia hybrida Hort. during growth of the pollen tube: cytochemical and radiographic studies. Acta Biol Cracov Bot. 1997, 39: 79-89.
Tintor N, Ross A, Kanehara K, Yamada K, Fan L, Kemmerling B, Nuernberger T, Tsuda K, Saijo Y: Layered pattern receptor signaling via ethylene and endogenous elicitor peptides during Arabidopsis immunity to bacterial infection. Proc Natl Acad Sci U S A. 2013, 110 (15): 6211-6216.
Elleman CJ, Dickinson HG: Commonalties between pollen/stigma and host/pathogen interactions: calcium accumulation during stigmatic penetration by Brassica oleracea pollen tubes. Sex Plant Reprod. 1999, 12 (3): 194-202.
Boutrot F, Segonzac C, Chang KN, Qiao H, Ecker JR, Zipfel C, Rathjen JP: Direct transcriptional control of the Arabidopsis immune receptor FLS2 by the ethylene-dependent transcription factors EIN3 and EIL1. Proc Natl Acad Sci U S A. 2010, 107 (32): 14502-14507.
Waithaka K, Dodge LL, Reid MS: Carbohydrate traffic during opening of gladiolus florets. Acta Hortic. 2001, 543: 217-226.
Wiemken V, Wiemken A, Matile P: Physiology of flowers of Ipomoea tricolor (Cav): studies on excised flowers and recovery of a phloem exudate. Biochem Physiol Pfl. 1976, 169 (4): 363-376.
Bieleski RL, Reid MS: Physiological changes accompanying senescence in the ephemeral daylily flower. Plant Physiol. 1992, 98 (3): 1042-1049.
Nichols R, Ho LC: Effects of ethylene and sucrose on translocation of dry-matter and 14C-sucrose in cut flower of glasshouse carnation (Dianthus caryophyllus) during senescence. Ann Bot. 1975, 39 (160): 287-296.
Yanagisawa S, Yoo SD, Sheen J: Differential regulation of EIN3 stability by glucose and ethylene signalling in plants. Nature. 2003, 425 (6957): 521-525.
van der Meulen-Muisers JJM, van Oeveren JC, van der Plas LHW, van Tuyl JM: Postharvest flower development in Asiatic hybrid lilies as related to tepal carbohydrate status. Postharvest Biol Tec. 2001, 21 (2): 201-211.
Lenartowska M, Rodriguez-Garcia MI, Bednarska E: Immunocytochemical localization of esterified and unesterified pectins in unpollinated and pollinated styles of Petunia hybrida Hort. Planta. 2001, 213 (2): 182-191.
Lenartowska M, Krzeslowska M, Bednarska E: Pectin dynamic and distribution of exchangeable Ca2+ in Haemanthus albiflos hollow style during pollen-pistil interactions. Protoplasma. 2011, 248 (4): 695-705.
O'Donoghue EM, Somerfield SD, Watson LM, Brummell DA, Hunter DA: Galactose metabolism in cell walls of opening and senescing petunia petals. Planta. 2009, 229 (3): 709-721.
O'Donoghue EM, Somerfield SD, Heyes JA: Vase solutions containing sucrose result in changes to cell walls of sandersonia (Sandersonia aurantiaca) flowers. Postharvest Biol Tec. 2002, 26 (3): 285-294.
O'Donoghue EM, Somerfield SD, Heyes JA: Organization of cell walls in Sandersonia aurantiaca floral tissue. J Exp Bot. 2002, 53 (368): 513-523.
De Vetten NC, Huber DJ: Cell wall changes during the expansion and senescence of carnation Dianthus caryophyllus petals. Physiol Plantarum. 1990, 78 (3): 447-454.
Shibuya K, Niki T, Ichimura K: Pollination induces autophagy in petunia petals via ethylene. J Exp Bot. 2013, 64 (4): 1111-1120.
Eisinger W: Role of cytokinins in carnation flower senescence. Plant Physiol. 1977, 59 (4): 707-709.
Saks Y, Vanstaden J, Smith MT: Effect of gibberellic-acid on carnation flower senescence: evidence that the delay of carnation flower senescence by gibberellic acid depends on the stage of flower development. Plant Growth Regul. 1992, 11 (1): 45-51.
Porat R, Halevy AH: Enhancement of petunia and dendrobium flower senescence by jasmonic acid methyl ester is via the promotion of ethylene production. Plant Growth Regul. 1993, 13 (3): 297-301.
Jofuku KD, Denboer BGW, Vanmontagu M, Okamuro JK: Control of Arabidopsis flower and seed development by the homeotic gene APETALA2 . Plant Cell. 1994, 6 (9): 1211-1225.
Yan X, Zhang L, Chen B, Xiong Z, Chen C, Wang L, Yu J, Lu C, Wei W: Functional identification and characterization of the Brassica Napus transcription factor gene BnAP2, the ortholog of Arabidopsis thaliana APETALA2 . PLoS One. 2012, 7 (3): e33890-
Karlova R, Rosin FM, Busscher-Lange J, Parapunova V, Do PT, Fernie AR, Fraser PD, Baxter C, Angenent GC, de Maagd RA: Transcriptome and metabolite profiling show that APETALA2a is a major regulator of tomato fruit ripening. Plant Cell. 2011, 23 (3): 923-941.
Izumi M, Wada S, Makino A, Ishida H: The autophagic degradation of chloroplasts via Rubisco-containing bodies is specifically linked to leaf carbon status but not nitrogen status in Arabidopsis. Plant Physiol. 2010, 154 (3): 1196-1209.
Guiboileau A, Yoshimoto K, Soulay F, Bataille M, Avice J, Masclaux-Daubresse C: Autophagy machinery controls nitrogen remobilization at the whole-plant level under both limiting and ample nitrate conditions in Arabidopsis. New Phytol. 2012, 194 (3): 732-740.
Guiboileau A, Avila-Ospina L, Yoshimoto K, Soulay F, Azzopardi M, Marmagne A, Lothier J, Masclaux-Daubresse C: Physiological and metabolic consequences of autophagy deficiency for the management of nitrogen and protein resources in Arabidopsis leaves depending on nitrate availability. New Phytol. 2013, 199 (3): 683-694.
Ono Y, Wada S, Izumi M, Makino A, Ishida H: Evidence for contribution of autophagy to Rubisco degradation during leaf senescence in Arabidopsis thaliana . Plant Cell Environ. 2013, 36 (6): 1147-1159.
Wang Y, Yu B, Zhao J, Guo J, Li Y, Han S, Huang L, Du Y, Hong Y, Tang D, Liu Y: Autophagy contributes to leaf starch degradation. Plant Cell. 2013, 25 (4): 1383-1399.
van Doorn WG, Beers EP, Dangl JL, Franklin-Tong VE, Gallois P, Hara-Nishimura I, Jones AM, Kawai-Yamada M, Lam E, Mundy J, Mur LAJ, Petersen M, Smertenko A, Taliansky M, Van Breusegem F, Wolpert T, Woltering E, Zhivotovsky B, Bozhkov PV: Morphological classification of plant cell deaths. Cell Death Differ. 2011, 18 (8): 1241-1246.
Chen WQ, Provart NJ, Glazebrook J, Katagiri F, Chang HS, Eulgem T, Mauch F, Luan S, Zou GZ, Whitham SA, Budworth PR, Tao Y, Xie ZY, Chen X, Lam S, Kreps JA, Harper JF, Si-Ammour A, Mauch-Mani B, Heinlein M, Kobayashi K, Hohn T, Dangl JL, Wang X, Zhu T: Expression profile matrix of Arabidopsis transcription factor genes suggests their putative functions in response to environmental stresses. Plant Cell. 2002, 14 (3): 559-574.
Besseau S, Li J, Palva ET: WRKY54 and WRKY70 co-operate as negative regulators of leaf senescence in Arabidopsis thaliana . J Exp Bot. 2012, 63 (7): 2667-2679.
Miao Y, Zentgraf U: A HECT E3 ubiquitin ligase negatively regulates Arabidopsis leaf senescence through degradation of the transcription factor WRKY53. Plant J. 2010, 63 (2): 179-188.
Chang X, Donnelly L, Sun D, Rao J, Reid MS, Jiang C: A petunia homeodomain-leucine zipper protein, PhHD-Zip, plays an important role in flower senescence. PLoS One. 2014, 9 (2): e88320-
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 (2): 825-840.
Dietrich K, Weltmeier F, Ehlert A, Weiste C, Stahl M, Harter K, Droege-Laser W: Heterodimers of the Arabidopsis transcription factors bZIP1 and bZIP53 reprogram amino acid metabolism during low energy stress. Plant Cell. 2011, 23 (1): 381-395.
Hummel M, Rahmani F, Smeekens S, Hanson J: Sucrose-mediated translational control. Ann Bot. 2009, 104 (1): 1-7.
Yuasa T, Nagasawa Y, Osanai K, Kaneko A, Tajima D, Htwe NMPS, Nang MPS, Ishibashi Y, Iwaya-Inoue M: Induction of a bZIP type transcription factor and amino acid catabolism-related genes in soybean seedling in response to starvation stress. J Botany. 2013, 2013 (Article ID 935479): 1-8.
Gerats T, Vandenbussche M: A model system comparative for research: Petunia . Trends Plant Sci. 2005, 10 (5): 251-256.
Schneider CA, Rasband WS, Eliceiri KW: NIH Image to ImageJ: 25 years of image analysis. Nat Methods. 2012, 9 (7): 671-675.
Zhong S, Joung J, Zheng Y, Chen Y, Liu B, Shao Y, Xiang JZ, Fei Z, Giovannoni JJ: High-throughput illumina strand-specific RNA sequencing library preparation. Cold Spring Harb Protoc. 2011, 2011 (8): 940-949.
Kofler R, Orozco-terWengel P, De Maio N, Pandey RV, Nolte V, Futschik A, Kosiol C, Schloetterer C: PoPoolation: a toolbox for population genetic analysis of next generation sequencing data from pooled individuals. PLoS One. 2011, 6 (1): e15925-
Martin M: Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011, 17 (1): 10-12.
Fu L, Niu B, Zhu Z, Wu S, Li W: CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics. 2012, 28 (23): 3150-3152.
Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M: Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005, 21 (18): 3674-3676.
Activities at the Universal Protein Resource (UniProt). Nucleic Acids Res. 2014, 42 (D1): D191-D198.
Zdobnov EM, Apweiler R: InterProScan - an integration platform for the signature-recognition methods in InterPro. Bioinformatics. 2001, 17 (9): 847-848.
Li H, Durbin R: Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009, 25 (14): 1754-1760.
VENNY: An interactive tool for comparing lists with Venn Diagrams. In , [http://bioinfogp.cnb.csic.es/tools/venny/index.html]
Chapin LJ, Jones ML: Ethylene regulates phosphorus remobilization and expression of a phosphate transporter (PhPT1) during petunia corolla senescence. J Exp Bot. 2009, 60 (7): 2179-2190.
Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc B. 1995, 57 (1): 289-300.
Jensen LJ, Kuhn M, Stark M, Chaffron S, Creevey C, Muller J, Doerks T, Julien P, Roth A, Simonovic M, Bork P, von Mering C: STRING 8-a global view on proteins and their functional interactions in 630 organisms. Nucleic Acids Res. 2009, 37: D412-D416.
Guo AY, He K, Liu D, Bai SN, Gu XC, Wei LP, Luo JC: DATF: a database of Arabidopsis transcription factors. Bioinformatics. 2005, 21 (10): 2568-2569.
Salaries and research support were provided by State and Federal funds appropriated to the Ohio Agricultural Research and Development Center, The Ohio State University (Journal Article Number HCS 14-09). Research was also funded by the D. C. Kiplinger Endowment, The American Floral Endowment Gus Poesch Fund, and SEEDS: The OARDC Research Enhancement Competitive Grants Program. We thank Jason Van Houten and Esther van der Knaap for help with the library construction and experimental design.
The authors declare that they have no competing interests.
Electronic supplementary material
Additional file 1: Biological replicate correlation scatterplots. Pearson correlation coefficients (R) were calculated between the normalized count data from each (A) pollinated and (B) unpollinated biological replicate and graphed. The blue line represents the slope of the Pearson correlation. (TIFF 978 KB)
Additional file 2: Differential gene expression results and Module I and III annotations. Differentially expressed genes were identified using pairwise comparisons (12-P/U, 18-P/U, 24-P/U and P/U) in DESeq2. Using BLASTx, gene annotation was obtained from the UniProt A. thaliana proteome. Each spreadsheet, which can be accessed through the tabs, contains the DESeq2 expression results for a single comparison, and two additional tabs contain the annotation results used for the KEGG enrichment of Modules I and III. (XLSX 1 MB)
Additional file 3: WGCNA cluster dendrogram. Dendrogram of modules based on biweight midcorrelation calculations from variance stabilized count data generated by DESeq2. The colors of the dynamic tree cut correspond to the modules assigned for each gene. The merged dynamic colors display the assigned changes when the stringency threshold of 0.75 was used. Module I (red), Module II (cyan) and Module III (grey60) were identified as pollination-associated modules. (TIFF 5 MB)
Additional file 4: Gene ontology enrichment results of differentially expressed genes. Gene ontology enrichment was performed on each list of differentially expressed genes from the DESeq2 pairwise comparison analyses (12-P/U, 18-P/U, 24-P/U and P/U). The enriched GO terms are organized by Biological Process, Cellular Component, and Molecular Function. Each spreadsheet, which can be accessed through the tabs, contains the GO term enrichment results for a single comparison. (XLSX 55 KB)
Additional file 5: Pathway mapping results for the genes belonging to the enriched pollination-associated KEGGs. TAIR codes were obtained from the top A. thaliana BLASTx results for each DESeq2 differentially expressed contig and WGCNA pollination-associated module contig. Through hypergeometric testing, 13 enriched KEGG pathways were identified in the differentially expressed 18-P/U and P/U gene lists, and six enriched KEGG pathways were identified from the WGCNA pollination-associated modules. Each gene name corresponds with the adjacent A. thaliana TAIR code. The KEGG pathway codes and names have boldface type. Each spreadsheet, which can be accessed through the tabs, contains the enriched KEGG pathways belonging to the DESeq2 pairwise results (i.e. 18-PU and PU) or to the designated WGCNA module (i.e. Module I and Module III). (XLSX 26 KB)
Additional file 6: 18-hour down-regulated KEGG pathways. TAIR codes from the top A. thaliana BLASTx hits of the differentially expressed 18-P/U genes were mapped to the KEGG database and the (A) Ribosome KEGG pathway, (B) Carbon fixation in photosynthetic organisms KEGG pathway, and the (C) Ribosome biogenesis in eukaryotes KEGG pathway were found to be enriched. Red boxes represent genes that were up-regulated 18 hours after pollination in petunia corollas, while dark green boxes represent down-regulated genes. Light green boxes represent A. thaliana genes that have been previously identified, while white boxes represent genes that belong to the KEGG pathway, but have no currently identified A. thaliana ortholog. (TIFF 4 MB)
Additional file 7: Primers used for quantitative PCR. Primers used to confirm expression patterns of select sequences from the RNA-seq analysis by quantitative PCR. (PDF 13 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
About this article
Cite this article
Broderick, S.R., Wijeratne, S., Wijeratn, A.J. et al. RNA-sequencing reveals early, dynamic transcriptome changes in the corollas of pollinated petunias. BMC Plant Biol 14, 307 (2014). https://doi.org/10.1186/s12870-014-0307-2