Transcriptome analysis of Rafflesia cantleyi flower stages reveals insights into the regulation of senescence

Rafflesia is a unique plant species existing as a single flower and produces the largest flower in the world. While Rafflesia buds take up to 21 months to develop, its flowers bloom and wither within about a week. In this study, transcriptome analysis was carried out to shed light on the molecular mechanism of senescence in Rafflesia. A total of 53.3 million high quality reads were obtained from two Rafflesia cantleyi flower developmental stages and assembled to generate 64,152 unigenes. Analysis of this dataset showed that 5,166 unigenes were differentially expressed, in which 1,073 unigenes were identified as genes involved in flower senescence. Results revealed that as the flowers progress to senescence, more genes related to flower senescence were significantly over-represented compared to those related to plant growth and development. Senescence of the R. cantleyi flower activates senescence-associated genes in the transcription activity (members of the transcription factor families MYB, bHLH, NAC, and WRKY), nutrient remobilization (autophagy-related protein and transporter genes), and redox regulation (CATALASE). Most of the senescence-related genes were found to be differentially regulated, perhaps for the fine-tuning of various responses in the senescing R. cantleyi flower. Additionally, pathway analysis showed the activation of genes such as ETHYLENE RECEPTOR, ETHYLENE-INSENSITIVE 2, ETHYLENE-INSENSITIVE 3, and ETHYLENE-RESPONSIVE TRANSCRIPTION FACTOR, indicating the possible involvement of the ethylene hormone response pathway in the regulation of R. cantleyi senescence. Our results provide a model of the molecular mechanism underlying R. cantleyi flower senescence, and contribute essential information towards further understanding the biology of the Rafflesiaceae family.


Results
Transcriptome assembly and functional annotation. Transcriptome sequencing was carried out on the perigone lobe of the F1 and F2 flower stages of the R. cantleyi flower (Fig. 1). Rafflesia flowers have a short blooming period of up to seven days before they wither and die. From our observation in the field (data not shown), flowers after three days of blooming showed significant fungal infections, which will affect endogenous gene expression. Thus, to minimize this effect on the transcriptome analysis, we chose to focus on the F1 and F2 flower stages in this study.

Differential expression analysis and GO enrichment. Differential expression analysis was conducted
to identify unigenes with significant differential expression between F1 and F2. The results of this analysis showed that a total of 5,166 unigenes identified were differentially expressed between F1 and F2. Of these, 4,062 unigenes were significantly up-regulated (log 2 FC ≥ 2; FDR ≤ 0.05) and 1,104 unigenes were significantly downregulated (log 2 FC ≤ -2; FDR ≤ 0.05). Further analysis from LSD 2.0 showed that 1,073 unigenes were identified to be involved in senescence, of which 715 unigenes were up-regulated and 358 unigenes were down-regulated (Supplementary Dataset S2).
Gene ontology (GO) annotation carried out on the differentially expressed genes (DEGs) showed the highest match to cell (GO:0005623) and organelle (GO:0043226) for the cellular component category; catalytic activity (GO:0003824) and binding (GO:0005488) for the molecular function category; while the biological process category was dominated by metabolic (GO:0008152) and cellular process (GO:0009987) (Fig. 3). GO enrichment analysis of the DEGs was carried out to identify significantly over-represented DEGs (Supplementary Dataset S3). For up-regulated DEGs, 22 GO terms were significantly enriched, including transportation, lipid metabolism, cell cycle, response to stress, cellular protein modification process, nucleotide and protein binding, and enzyme regulator activity. In contrast, 27 GO terms for the down-regulated DEGs were significantly enriched, including meiotic nuclear division, response to biotic stimulus, post-embryonic development, and cell growth.
Other than that, the DEGs identified showed functions in the regulation of transcription factors, nutrient remobilization, redox regulation, and ethylene regulation ( www.nature.com/scientificreports/ (ATG ) and transporter genes involved in nutrient remobilization were up-regulated. A slight difference in the expression pattern was recorded for unigenes involved in redox regulation, in which most genes were constitutively expressed at both stages. This pattern of constitutive expression was similar to unigenes involved in the ethylene biosynthesis but unigenes involved in the ethylene signal transduction showed significant differential expression (Fig. 4).  between the two flower stages were selected for RT-qPCR validation analysis. Results showed that QUINOLI-NATE SYNTHASE (QS) has the highest increase in expression, followed by ETHYLENE RECEPTOR (ETR), ATG4, and CATALASE (CAT ) (Fig. 5a), while PECTIN ESTERASE (PE) has the highest decrease in expression, followed by MYB and AGL8 (Fig. 5b). Overall, the RT-qPCR results showed a high correlation with that of RNAseq (R 2 > 0.8742) (Fig. 5c).

Discussion
In this study, besides the more commonly used databases such as NR and Swiss-Prot, a BLAST search was also performed against the leaf senescence-specific LSD 2.0 database. In the absence of a specific database for flower senescence, the LSD 2.0 database was utilized since most of the leaf senescence genes can also be found in the flower. Furthermore, leaf and floral organs arise from the same primordium, thus it is expected that the senescence genes are conserved in these two organs, as observed in A. thaliana 22 and Erysimum linifolium 23 . BLAST results showed that the majority of the homologs are from A. thaliana (79.14%) and M. acuminata (16.98%), and this may be partly due to both of these species having the largest datasets in the LSD 2.0 database, 69.89% for A. thaliana and 16.46% for M. acuminata. GO annotation revealed that DEGs were mainly responsible for fundamental biological regulation and metabolism common in other plants. This suggests that R. cantleyi flower senescence is a coordinated biological process that could be regulated similarly to other plants 24 . The transcriptional control mechanism, which involves various transcription factors, plays a pivotal role in regulating gene expression during plant development and senescence 25 . Many transcription factors have been identified previously and found to be differentially expressed during the development and senescence of various flower systems 17 . In R. cantleyi, members of MYB, bHLH, NAC, and WRKY were found to be up-regulated at F2. All of these families were known to be involved in various kinds of biological processes, including responses to stresses and injuries [26][27][28][29] . MYB102 functions in response to biotic and abiotic stresses, such as defense against insects, injuries, and osmotic stresses [30][31][32] , whereas, MYB7 is a transcription factor involved in the flavonoid biosynthesis induced during salinity stress 33 . On the other hand, the function of bHLH family members MYC2 34 , bHLH092 27 , and ICE1 [35][36][37] in abiotic stress response has been well characterized. Moreover, expression data from genome-wide transcriptome analyses in many plants such as A. thaliana 38 , soybean 39 , rice 40 , and Chrysanthemum lavandulifolium 41 revealed that a significant proportion of NAC genes are related to stress responses 42 . Studies conducted on the expression of WRKY showed rapid induction when plants are exposed to a variety of stresses or defensive signals, including plant defense against attack and senescence 29,43,44 .
The up-regulation of these transcription factors is in line with other studies on flower senescence. The same phenomenon can be seen in Dianthus and Gardenia, which showed a high level of expression in MYB during flower senescence 45,46 . Other than that, the expression pattern of bHLH between the two R. cantleyi flower stages is similar to that in Gardenia, whereby this transcription factor family plays a major role in regulating flower senescence 46,47 . The expression pattern of NAC in the R. cantleyi flower was similar to that of A. thaliana leaf and flower senescence 22,48 .
The transcription factors that showed significant down-regulation in R. cantleyi were identified to be from the bZIP, MADS-box, and AUX/IAA families, which play important roles in flower development [49][50][51] . The transcription factors from the bZIP family plays important roles in regulating seed maturation and flower development 49 . These functions are also shared with the MADS-box transcription factors, which are known to be involved in the development of the reproductive structure and determination of meristems and organ identity 50 . Additionally, Flowers have a fixed lifespan determined by their status in sexual reproduction such that flowers that have been pollinated or are no longer receptive to pollination will be programmed to senensce 54 . In some plant species, the symptoms of flower senescence are visibly shown through wilting or withering of petals, whereas in other plants, the petals abscise while still turgid 55 . For plant species with petal wilting, such as in Rafflesia, the main function of flower senescence is proposed to allow the remobilization of nutrients to other parts of the plant 16,56,57 . One of the processes during the remobilization of nutrients is protein degradation, which involves multiple protease enzymes. These enzymes degrade proteins by hydrolyzing internal peptide bonds, or more commonly known as proteolysis, and this activity is irreversible.
Cysteine protease is an enzyme that has been reported to be involved in the remobilization of essential nutrients from senescing floral tissue. In some plant species such as Sandersonia 58 , Narcissus 59 , Alstroemeria 56 , and  www.nature.com/scientificreports/ Rosa hybrida 60 , genes encoding cysteine protease were found to be induced during flower senescence. In this study, a unigene encoding cysteine protease in the R. cantleyi flower was identified. Differential expression analysis showed that the cysteine protease gene ATG4 (UN020537) was up-regulated at F2 (Table 1; Supplementary Dataset S4). This gene was identified to be involved in autophagy, which is the main pathway for the degradation of proteins in vacuole activated by environmental stresses 61 .
Other than that, unigenes encoding other related autophagy proteins were also identified. Besides ATG4, the up-regulation of ATG7, ATG8, and ATG9 at F2 were also similar to those in A. thaliana during leaf senescence 62 . The expression patterns of ATG4 and ATG8 in R. cantleyi are similar to Ipomea nil and Petunia with increasing expression during petal senescence [63][64][65] . Furthermore, the activation of the autophagy pathway is negatively regulated by TARGET OF RAPAMYCIN (TOR). The down-regulation of TOR caused the autophagy pathway activation in A. thaliana 66 . In this study, a unigene that encodes TOR was down-regulated at F2, indicating that the autophagy mechanism in the R. cantleyi flower is probably activated by the regulation of TOR.
Besides that, unigenes involved in transporter activity, such as ABC transporter and transporters of sugar, nitrates, and polyamines were also identified as DEGs (Table 1; Supplementary Dataset S4). In A. thaliana, the expression of transporter protein was found to be up-regulated during leaf senescence, which reflects nutrient remobilization during senescence 67 . Overall, the differential expression of unigenes related to nutrient remobilization suggests transportation of nutrients from the senescing flower of R. cantleyi, perhaps for the development of seeds and fruits.
The reactive oxygen species (ROS) represent free radicals derived from oxygen molecules, such as singlet oxygen ( 1 O 2 ), superoxide (O 2 · − ), hydroxyl (OH·), hydroperoxyl (HO 2 · − ), and hydrogen peroxide (H 2 O 2 ) 68 . These radicals can oxidize various molecules in plants, which gives adverse effects on membrane integrity and protein stability leading to programmed cell death 69 . ROS are produced during normal plant development and the production increases when triggered by environmental stresses such as drought, salinity, and pathogen attack. ROS are extremely harmful in high concentrations, thus they need to be removed from the cell through various activities of antioxidative enzymes, such as superoxide dismutase (SOD), catalase (CAT), and ascorbate oxidase (ASO). In previous studies of flower senescence, such as in Gladious 70 , Hemerocallis 71,72 , and Rosa hybrida 73 , various free oxygen radicals were detected in flower petals. The increase of these free radicals caused the upregulation of antioxidative enzymes SOD, CAT, and ASO during flower senescence 55 . In this study, unigenes that are involved in redox regulation encoding antioxidative enzymes CAT (UN004263), SOD (UN012078), ASO (UN037470), and other peroxidases (UN001999, UN063210) were identified (Table 1; Supplementary Dataset S4). Our results revealed that these unigenes showed constitutive expressions at both flower stages, except CAT (UN004263), which was found to be up-regulated in F2. High expression of ROS-related unigenes at both flower stages indicates the regulation of redox activity in R. cantleyi.
Flower senescence involves highly coordinated events that are regulated by endogenous signals, such as phytohormones and environmental factors, including biotic and abiotic stresses 74 . Several plant hormones have been reported to be involved in senescence, of which ethylene and abscisic acid (ABA) function to induce flower senescence, while cytokinin, gibberellic acid (GA), and auxin delay the process 16,75 . In most plant species, ethylene is a phytohormone that initiates flower senescence 76 . Results from differential expression analysis were able to help identify genes involved in ethylene hormone regulation in R. cantleyi.
Ethylene is synthesized through cysteine and methionine metabolism pathways mediated by enzyme precursors ACO and ACS 77 . In this study, unigenes encoding ACO (UN064036) and ACS (UN023149) were constitutively expressed at both flower stages. However, unigenes involved in the ethylene signal transduction pathway were identified to be differentially expressed. Among them, positive regulators ETHYLENE-INSENSITIVE 2 (EIN2) (UN007860), ERF (UN012030), and ETR (UN024806), and transcription factor ETHYLENE-INSENSI-TIVE 3 (EIN3) (UN025327) were significantly up-regulated, whereas a negative regulator CTR1 (UN004001) was significantly down-regulated at F2 (Table 1; Supplementary Dataset S4). ETR serves as a receptor that binds to ethylene gases and the loss of ETR function delays senescence in flowers and leaves as recorded in Petunia and Nicotiana sylvestris [78][79][80] . On the other hand, in the ethylene signal transduction pathway, the binding of ethylene to the ETR receptor causes the CTR gene to be inactive. The deactivation of CTR activates other downstream components such as transcription factors EIN2, EIN3, and ERF, followed by ethylene response 17 . Constitutive expression of unigenes involved in the ethylene biosynthesis and the up-regulation of unigenes involved in the ethylene signal transduction pathway suggest the role of ethylene regulation in flower senescence of R. cantleyi (Fig. 4).
Results of the transcriptome data analysis allowed us to build a model for the regulation of senescence in R. cantleyi, based on a senescence regulation pathway proposed by Tripathi and Tuteja (2007) 81 (Fig. 6). Senescence involves death at the cellular level known as PCD 13 . In R. cantleyi, PCD might be induced by environmental factors, such as biotic and abiotic stresses, which hasten the process of flower senescence and lead to the regulation of ethylene hormone and various transcription factors. Various pathways involving phytohormones and the regulation of transcription and signal transduction are known to activate the expression of senescence-associated genes 17 . The regulation of these transcription factors suggests activation of senescence-associated genes that are involved in the process of nutrient remobilization and redox regulation.

Conclusion
Transcriptome generation and analysis focused on DEGs for two stages of R. cantleyi flowers yielded a model on the regulation of flower senescence in R. cantleyi. This model consists of the regulation of transcription factors and ethylene-related genes, together with nutrient remobilization, redox regulation, and the activation of senescence-associated genes, which ultimately lead to a coordinated death of a floral organ. The understanding  (Fig. 1). Tissue samples were acquired from the perigone lobes of two different flowers that were collected from the same area of sampling and under a relatively similar environment. These samples were dissected and surface-sterilized using 10% (v/v) Clorox and rinsed with distilled water on site. Following that, the tissue samples were flash-frozen in liquid nitrogen before being transported to the laboratory. All the samples collected were stored at − 80 °C until further use. Total RNA extraction was carried out  Data pre-processing and de novo assembly of transcriptome. Raw sequencing reads generated (deposited to the NCBI SRA database with the accession numbers SRR7544088 (F1) and SRR7544087 (F2)) were quality-filtered using SolexaQA (v.2.2) 83 . High-quality reads were acquired by removing adaptor sequences and low-quality reads with Phred score less than 20 and read length less than 50 bp. Furthermore, the reads were screened for contaminating sequences from other organisms by aligning them using Bowtie2 (v2.3.0) 84 to all genomes of bacteria (version 66, 7 July 2014), viruses, and fungi (version 73, 2 Nov 2015) from the NCBI database. Reads that aligned to these genomes were compared against the plant database Plaza 3.0 85 and reads without any alignment to the plant genomes were discarded. Finally, the paired-end clean reads were subjected to de novo assembly using Trinity (version v2.2.0) 86 and iAssembler (version v1.3.2) 87 . The R. cantleyi flower transcriptome assembly was deposited in the TSA database with the accession number GIQT00000000.
Functional annotation of R. cantleyi transcriptome. To predict the putative functions of assembled unigenes, functional annotation was performed using BLASTX against NR and Swiss-Prot protein databases with an E-value cutoff of 1E-6. Unigenes involved in senescence were identified by performing similarity searches against 44 species of dicots and monocots in the Leaf Senescence Database (LSD 2.0) 88 . Gene ontology (GO) annotation was implemented on the annotated sequences from NR using the Blast2GO program and was further classified into three categories: biological process, molecular function, and cell component. The summary of GO classification was visualized using WEGO 89 . GO enrichment analysis was performed using Fisher's exact test (FDR < 0.05) based on the up-regulated and down-regulated unigenes to predict the overrepresented GO function. KEGG pathway assignment was performed on unigenes using KEGG Automatic Annotation Server (KAAS) 90 to predict metabolic pathways involved in the R. cantleyi flower.
Differential gene expression analysis. Gene expression levels were estimated by mapping paired-end clean reads from each flower stage (F1 and F2) to the assembled transcriptome using the RSEM method 91 . Due to the great difficulties of on-site sampling of specific Rafflesia flower stages in the rainforest and the inability of cultivating Rafflesia for study under laboratory or controlled environment conditions, only one biological replicate was available for transcriptomic analysis. The gene expression was based on TPM (transcripts per kilobase million) that normalizes transcript length first before sequencing depth. Unigenes with TPM values greater than 0 were considered as expressed. By pairwise comparisons of the two libraries (F1 vs. F2), the DEGs were identified using the edgeR package 92 based on the read counts of unigenes in different libraries. Trimmed Mean of M-values (TMM) normalization with an automatic method to estimate dispersion without replicates allowed single-replicate DEG analysis in edgeR, assuming that the counts are not too small. The significance of differential gene expression was judged by the false discovery rate (FDR ≤ 0.05) and fold change (|Log 2 FC|≥ 2) as the cut-off threshold.
Gene expression analysis by RT-qPCR. Seven DEGs related to transcription factors, nutrient remobilization, redox regulation, and ethylene hormone pathway were selected from the transcriptome dataset and expression profiles were investigated by reverse transcription-quantitative polymerase chain reaction (RT-qPCR) analysis. Gene-specific primers (Supplementary Dataset S5) were designed using Primer-Blast. TUBU-LIN (TUB) and ELONGATED FACTOR 1A (EF1A) were selected as reference genes after they were identified as the most stable genes in the R. cantleyi flower by the Bestkeeper software 93 . For RT-qPCR analysis, 20 ng/µL DNAse-treated RNA was used as a template in a 20 µL total reaction with QuantiNova® SYBR® Green RT-PCR Kit (Qiagen, USA) according to the manufacturer's instructions. The reactions were incubated in the CFX96™ Real-Time Detection System (Bio-Rad, USA) with the following cycling conditions: 10 min at 50 °C and 2 min at 95 °C for the reverse transcription, 40 cycles of 5 s at 95 °C and 10 s at 60 °C. In this analysis, three perigone lobes of each flower stage were treated as biological replicates, and each PCR reaction included three technical replicates. The changes in gene expression were calculated based on the 2 −∆∆Ct method 94 . Statistical analysis of the data was performed using unpaired Student's t-test at a P < 0.05 significance level. Correlation analysis of relative expression values from RNA-seq and RT-qPCR was performed using Microsoft Excel.