De novo sequencing and analysis of the transcriptome during the browning of fresh-cut Luffa cylindrica 'Fusi-3' fruits

Fresh-cut luffa (Luffa cylindrica) fruits commonly undergo browning. However, little is known about the molecular mechanisms regulating this process. We used the RNA-seq technique to analyze the transcriptomic changes occurring during the browning of fresh-cut fruits from luffa cultivar ‘Fusi-3’. Over 90 million high-quality reads were assembled into 58,073 Unigenes, and 60.86% of these were annotated based on sequences in four public databases. We detected 35,282 Unigenes with significant hits to sequences in the NCBInr database, and 24,427 Unigenes encoded proteins with sequences that were similar to those of known proteins in the Swiss-Prot database. Additionally, 20,546 and 13,021 Unigenes were similar to existing sequences in the Eukaryotic Orthologous Groups of proteins and Kyoto Encyclopedia of Genes and Genomes databases, respectively. Furthermore, 27,301 Unigenes were differentially expressed during the browning of fresh-cut luffa fruits (i.e., after 1–6 h). Moreover, 11 genes from five gene families (i.e., PPO, PAL, POD, CAT, and SOD) identified as potentially associated with enzymatic browning as well as four WRKY transcription factors were observed to be differentially regulated in fresh-cut luffa fruits. With the assistance of rapid amplification of cDNA ends technology, we obtained the full-length sequences of the 15 Unigenes. We also confirmed these Unigenes were expressed by quantitative real-time polymerase chain reaction analysis. This study provides a comprehensive transcriptome sequence resource, and may facilitate further studies aimed at identifying genes affecting luffa fruit browning for the exploitation of the underlying mechanism.


Introduction
Luffa cylindrica (i.e., luffa) of the family Cucurbitaceae is one of the most important vegetables and widely used medicinal plants in China. People are increasingly starting to pay more attention to diet, nutrition, and health. Because luffa is nutritious, delicious, and has beneficial effects on human health, the land area on which it is cultivated continues to increase [1][2][3]. There are currently two luffa species, namely towel gourd (Luffa cylindrica Roem.) and sinkwa a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 The total RNA used for a quantitative real-time polymerase chain reaction (qRT-PCR) assay was extracted from each sample using the EZNA Plant RNA Kit (Bio-tek, Beijing, China). The quantity and quality of the RNA samples were assessed using the NanoDrop ND-1000 spectrophotometer (Thermo Scientific, CA, USA) and a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The libraries were then sequenced using the HiSeq™ 2000 system (Illumina Inc.). Equal amounts of RNA from three samples collected at the same browning time point were combined.

Determination of total phenol content
The total phenol content was determined with the Folin-Ciocalteu assay, using a modified version of the method described by Dewanto et al. (2002) [24]. The absorbance (750 nm) of a 25-μL aliquot of extract was recorded, and the total phenol concentration was expressed as gallic acid equivalents (mg GAE g −1 fresh weight) based on a calibration curve (1-6 μg mL −1 gallic acid). Each biological replicate was measured three times.

De novo transcriptome assembly and functional annotation
The RNA sequencing (RNA-seq) raw data were filtered by removing adapter sequences with in-house Perl scripts to obtain high-quality reads. Reads with low-quality regions (reads with a base quality < 20) were excluded using the sliding window trimming approach. The remaining reads were de novo assembled using Trinity software (http://trinityrnaseq.github.io/) [25], with min_kmer_cov set to 2 by default and all other parameters also set to their default values.
We used the luffa Unigene sequences to search the National Center for Biotechnology Information nonredundant (NCBInr) protein, Gene Ontology (GO), Eukaryotic Orthologous Groups of proteins (KOG), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Swiss-Prot databases. Protein functional annotations were based on a BLASTX search (E < 10 −5 ), which was used to identify the most similar proteins. The GO terms were assigned to Unigenes using Blast2GO [26]. The distribution of the GO functional classifications for the Unigenes was plotted using WEGO software [27]. We identified the genes potentially important for luffa browning using the aforementioned databases.

Analysis of Unigene expression levels
Reads from four cDNA libraries were mapped to the assembled Unigenes using Bowtie [28]. A fragments per kilobase of exon model per million mapped reads (FPKM) fold-change > 2 and a false discovery rate (FDR) < 0.05 were used as the thresholds to identify significant differences in gene expression [29]. Expression-level differences between two samples were examined using the DEGseq R package [30]. A fold-change ! 2 (in both directions; log 2 ratio ! 1) and FDR 0.05 were used as standards to determine the significance of gene expression level differences among samples. A time series cluster analysis, which is based on the short timeseries expression miner (STEM) method (http://www.cs.cmu.edu/~jernst/stem/) [31], was used to identify global trends and similar temporal model expression patterns among all differentially expressed genes (DEGs).

Quantitative real-time polymerase chain reaction analysis
The RNA-seq data were validated by qRT-PCR using gene-specific primers. Total RNA was extracted from fruits, and a 2-μg sample was used to prepare cDNA with the PrimeScript™ II 1st Strand cDNA Synthesis kit (Takara). The resulting cDNA was used as the template for qRT-PCR, which was conducted using a real-time PCR plate (Applied Biosystems, Foster City, CA, USA), SYBR green as the fluorescent reagent, and an ABI 7500 Fast Real-time PCR system (Applied Biosystems). A standard curve for each target gene was constructed based on qRT-PCR data from a series of diluted cDNA samples. The 20-μL qRT-PCR solutions consisted of 2 μl cDNA, 10 μL SYBR Premix Ex Taq™ II, 6.8 μL EASY Dilution buffer, and 0.6 μL 10 μM forward and reverse primers. The PCR was conducted with the following program: 95˚C for 30 s; 40 cycles of 95˚C for 10 s and 60˚C for 34 s. The 18S rRNA gene (GenBank accession: KM656452), which was stably expressed according to the RNA-seq data, was used to normalize expression levels [32]. The 2 −ΔΔCt method was used to calculate relative gene expression levels with Microsoft Excel software. The qRT-PCR experiment consisted of three independent biological replicates, and the expression levels calculated for each sample were based on three technical replicates.

Determination of total phenol content
Previous studies indicated that among Luffa cylindrica cultivars, 'Fusi-3' produces fruits that are the most susceptible to browning [33]. Additionally, the total phenol content of browning fresh-cut vegetables under constant temperatures tends to increase over time [13,34]. For a comprehensive transcriptome analysis of luffa fruits, total RNA was extracted from peeled 'Fusi-3' fruit slices. The total phenol content exhibited an increasing trend during the browning treatment period (i.e., 0-6 h) (Fig 1).

RNA sequencing and assembly of reads
To analyze the transcription profiles of luffa fruits undergoing browning, four cDNA libraries were created from 'Fusi-3' fruits at 18 days after pollination, and analyzed using the Illumina HiSeq™ 2000 system. We generated about 92 million sequence reads. After removing adapter sequences and low-quality reads, we obtained 91 million valid sequencing reads, with an average length of 90.2 bp. After a de novo assembly analysis, 58,073 Unigenes were obtained, with an average length of 896.93 bp and an N50 of 1,510 bp. Additionally, 29.86% of the Unigenes were > 1,000 bp long ( Table 1).

Sequence annotation
The annotation of Unigenes provided information regarding functions, protein sequence similarities, and KOG, GO, and KEGG pathway details. The data were screened against the following four public databases: NCBInr, Swiss-Prot, KOG, and KEGG. Unigenes were annotated according to the best BLASTX matches (E < 10 −5 ) and sequence identities > 30%. A total of 35,345 (60.86%) Unigenes were annotated. The assembled sequences included 37,321 (60.75%), 24,427 (30.87%), 20,546 (30.52%), and 13,021 (18.83%) Unigenes with matches in the NCBInr, Swiss-Prot, KOG, and KEGG databases, respectively (S1 Table and Fig 2).  (Fig 3). The fact that the luffa Unigenes had the most matches with C. melo sequences suggested that C. melo may be a useful reference material for further analyses of luffa functional genes. A total of 27,858 Unigenes were annotated by one or more databases (i.e., 78.81% of all assembled Unigenes), implying they have relatively well conserved functions.
The luffa Unigenes were searched against the KOG database to predict and classify their possible functions. Of 37,321 NCBInr matches, 31,567 sequences had KOG classifications among 25 KOG categories (Fig 4). The most common categories were "general function prediction only", with 7,857 Unigenes, and "posttranslational modification, protein turnover, and chaperones" (3,626 Unigenes). "Cell motility" had the fewest Unigenes, with only 17. Our data revealed that annotated Unigenes were associated with most of the activities related to growth and development.
The GO assignments were used to classify the functions of the luffa transcripts. Based on sequence homologies, the 24,021 annotated sequences with BLAST hits to PlantGDB proteins were categorized into 45 functional groups (Fig 5). "Metabolic process", "cellular process", and "binding" were the dominant terms in the "biological process", "cellular component", and "molecular function" main categories, respectively. We also identified a relatively large number of genes associated with "catalytic activity", "single-organism process", and "cell", with only a few genes related to "extracellular region part", "nucleoid", "extracellular matrix", and "locomotion".
Pathway-based analyses may be useful for clarifying the biological functions of genes. The KEGG pathway database contains information regarding networks of intracellular molecular interactions, as well as their organism-specific variations [35]. To identify the biological pathways that are active during the browning of fresh-cut luffa fruits, we mapped the 58,073 Unigenes to the reference canonical pathways in the KEGG database. A total of 13,021 Unigenes were annotated (S2 Table). The top 10 pathways (i.e., "ribosome", "carbon metabolism", "plant hormone signal transduction" "biosynthesis of amino acids", "starch and sucrose metabolism", "protein processing in endoplasmic reticulum", "endocytosis", "spliceosome", "purine metabolism", and "RNA transport") and three potential browning-related pathways (i.e., "phenylpropanoid biosynthesis", "peroxisome", "tyrosine metabolism") are presented in  Moreover, there were 47, 131, and 232 Unigenes associated with the "tyrosine metabolism", "peroxisome" and "phenylpropanoid biosynthesis" pathways, respectively. These results indicated that diverse metabolic processes are active in luffa fruit flesh, resulting in the synthesis of diverse metabolites. These annotations may serve as a valuable resource for investigating specific processes, functions, and pathways, and facilitate the identification of novel genes involved in the browning of luffa fruit flesh.

Analysis of differentially expressed genes
Based on a published protocol, we developed an algorithm to identify genes that were differentially expressed between two samples [36]. The FDR was used to determine the P-value threshold in multiple analyses for calculating the expression-level differences between two samples [37]. We calculated the number of expressed sequence tags corresponding to each gene in each library to estimate gene expression levels. A total of 27,301 Unigenes were differentially expressed during the browning of luffa fruits (S3 Table). The highest number of DEGs was detected during the comparison between the 0 and 6 h time points [i.e., 16,509 (3,489 up-and 13,020 down-regulated at 6 h) (log 2 ratio ! 1; FDR 0.05; P 0.01)] (Fig 7). We used the STEM method to refine the sets of genes that were differently expressed at a minimum of four time points (i.e., 0, 1, 3, and 6 h). This method is commonly used for gene clustering in transcriptomic studies [31]. A total of 27,301 DEGs were clustered into 26 possible model profiles based on expression patterns (Fig 8).
The DEGs were then subjected to GO and KEGG pathway enrichment analyses.
processes of interest or metabolic pathways. We conducted pathway enrichment analyses to identify significantly enriched metabolic or signal transduction pathways associated with the DEGs [38]. According to GO functional annotations and the screening of 58,073 Unigenes, we  identified 11 Unigenes with potential roles in the browning of fresh-cut luffa fruits, including three PPO genes (Unigene0043875, Unigene0050171, and Unigene0015096), two PAL genes (Unigene0044748 and Unigene0024978), one POD gene (Unigene0011770), two CAT genes (Unigene 0036262 and Unigene0033876), and three SOD genes (Unigene0021782, Uni-gene0008835, and Unigene0015506), as well as four WRKY transcription factors (S1 Table). The Unigenes were annotated based on the NCBInr protein database and divided into the following categories: oxidoreductase, ammonia lyase, antioxidant, catalytic activity, and DNAbinding function related. Furthermore, a KEGG analysis revealed that 15 genes may affect tyrosine metabolism, phenylpropanoid biosynthesis, or the peroxisome pathway (S4 Table). Additionally, with the aid of rapid amplification of cDNA ends technology, we obtained the full-length CAT (Unigene0036262) sequence, which has been deposited in the GenBank database ( Table 2).

Validation of differentially expressed genes based on RNA-seq data by quantitative real-time PCR
The heat maps for the 11 browning-related genes and four WRKY transcription factors indicated that the expression levels varied at different time points (i.e., 0, 1, 3, and 6 h) (Fig 9). We used qRT-PCR with gene-specific primers ( Table 3) to analyze the expression profiles of genes likely involved in the browning of fresh-cut luffa fruits. The expression levels were mainly upregulated relative to the 0 h expression levels, except for SOD (Unigene0021782), which was down-regulated at the 1 h time point, and PPO (Unigene0015096), CAT (Unigene0033876), and SOD (Unigene0008835 and Unigene0015506), which were down-regulated at the 1, 3, and 6 h time points (Fig 10). A linear regression analysis revealed a strong positive correlation (R 2 = 0.9694; P 0.01) between the fold-change values based on qRT-PCR and RNA-seq data (Fig 11). These results confirm the reliability of the RNA-seq data generated in this study.

Discussion
Luffa is an important vegetable crop worldwide. It is rich in nutrients and medicinal compounds, including carotenoids, β-carotene, and α-carotene, which maintain the immune system [39]. To identify genes involved in the browning of luffa fruits, total RNA was extracted from 'Fusi-3' fruits and used for mRNA preparation, fragmentation, and cDNA synthesis. The cDNA was sequenced using the Illumina Genome Analyzer IIx platform, and the resulting sequencing data were subjected to bioinformatic analysis. We generated 92,010,328 raw reads from a 200-bp insert (mean size) cDNA library. After trimming adapter sequences and removing low quality reads, 90,226,276 clean reads were obtained. Of these clean reads, 94.19% were considered high-quality reads, which were analyzed further. We annotated 35,345 Unigenes (i.e., 60.86% of the assembled Unigenes), suggesting their functions were relatively conserved. In addition to identifying genes or proteins, we also analyzed the associated GO terms and potential roles in metabolic pathways. Combining our data regarding Unigenes expressed in fresh-cut luffa fruits with the results of a previous transcriptomic study of the same luffa variety reveals that luffa fruits may be a useful source of abundant genes [40]. The current study represents the first attempt at applying Illumina paired-end sequencing technology to investigate the changes in the transcriptome during the browning of fresh-cut luffa fruits. Enzymatic processes are generally recognized as the main determinants of browning. The accumulation of transcripts for genes associated with the oxidation of phenolic compounds is a hallmark of enzymatic browning. Many studies have confirmed that browning-related genes are differentially expressed in fresh-cut vegetables [41][42][43]. The PPO and POD genes have been associated with the browning of fruit and vegetable tissues, and their encoded enzymes function together to induce browning [44,45]. The PPO enzyme encoded by members of the PPO multigene family is considered to be a major factor responsible for enzymatic browning of several fresh-cut vegetables [17,46]. Zhou et al. (2015) [47] reported that the enzymatic browning of fresh-cut Chinese chestnut can be inhibited by applying salicylic acid to decrease PPO activity. Their findings provide indirect evidence that PPO influences the browning of fresh-cut Chinese chestnut. Moreover, a previous study revealed that a 3-day treatment with 1 mg L −1 ozone decreased POD activity in fresh-cut lettuce, resulting in obviously inhibited browning [48]. Additionally, Chisari et al. (2007) [46] reported that the POD enzyme also affected the browning of fresh-cut melon. In this study, the observed total phenol contents and the results of a KEGG analysis revealed that three PPO genes (Unigene0015096, Unigene0043875, and Unigene0015096) and one POD gene (Unigene0011770) related to the tyrosine metabolism pathway were differentially expressed while the phenolic acid (i.e., gallic acid) content increased in fresh-cut luffa fruit slices (S4 Table) [49,50]. The production of phenolic acids, including gallic acid, L-tyrosine, chlorogenic acid, and gentisic acid, catalyzed by browningrelated enzymes (i.e., PPO and POD) can be affected by the tyrosine metabolism pathway [51][52][53][54]. However, the exact roles of PPO and POD during responses to the browning of fresh-cut luffa fruits remains to be determined.
In addition to PPO and POD, the PAL, CAT, and SOD genes are also considered to have important roles in the enzymatic browning of fruits and vegetables [55,56]. The PAL enzyme is reportedly involved in the browning of fresh-cut potatoes [13]. Additionally, a lettuce PAL homolog was observed to be responsible for regulating the browning of fresh-cut edges [57]. Our results indicated that the relative expression levels of two luffa PAL genes (Uni-gene0044748 and Unigene0024978) involved in the phenylpropanoid biosynthesis pathway were up-regulated (S4 Table). It is widely accepted that phenylpropanoid biosynthesis pathways are important for controlling the browning of cut vegetables and fruits (S4 Table) [58][59][60][61]. Additionally, the CAT enzyme has a pivotal role in the browning of fresh-cut fruits and vegetables [62][63][64]. Dong et al. (2015) [44] proposed that the activity of CAT can be enhanced by a short-term carbon dioxide treatment, leading to inhibited browning of fresh-cut burdock during storage at low temperatures. Meanwhile, Fan et al. (2016) [65] reported that significantly higher activities of antioxidant enzymes (i.e., SOD and CAT) were maintained in browning 'Laiyang' pear fruits, although PPO activity was inhibited. Another important enzyme related to the browning of fresh-cut fruits and vegetables is SOD, which is essential for enhancing antioxidant activities in fresh-cut Chinese water chestnut [66]. In this study, two CAT genes (Unigene0036262 and Unigene0033876) and three SOD genes (Unigene0021782, Uni-gene0008835, and Unigene0015506) related to the peroxisome pathway were also differentially expressed during the browning of fresh-cut luffa fruits (S4 Table).
The browning of fruits and vegetables is caused by the overproduction of reactive oxygen species [67][68][69]. Members of the WRKY transcription factor superfamily have been detected exclusively in plants, and affect diverse physiological responses and metabolic pathways [70,71]. Previous studies confirmed that WRKY transcription factors are important inducers of the expression of genes associated with the reactive oxygen species pathway during plant responses to biotic, abiotic, and oxidative stresses [72,73]. In this study, we observed that the expression levels of four Unigenes (Unigene0018509, Unigene002 1412, Unigene0025291, and Unigene0034271) annotated as WRKY transcription factors were significantly up-regulated at various time points (i.e., 0, 1, 3 and 6 h) during the browning of fresh-cut luffa fruits. Analyses revealed that these four transcription factors (GenBank ID: KY621843-KY621846) were members of the Group II WRKY superfamily, and contained a conserved WRKY domain and C 2 H 2 -type zinc fingers [74]. Furthermore, the current study is the first to reveal a relationship between these four WRKY transcription factors and the browning of luffa fruits. Further studies involving structural and functional characterizations are required to determine whether the identified candidate genes influence the browning of fresh-cut luffa fruits.
Transcriptome profiles of different luffa cultivars showed that the expression levels of PPO POD and PAL genes related with browning in XTR05 (browning sensitive) were significantly higher than that in YLB05 (browning resistant) [40]. In this study, we also identified three PPO genes, one POD gene and two PAL genes were differentially expressed during the browning of fresh-cut luffa fruits (i.e., after 1-6 h). Moreover, we identified two CAT genes, three SOD genes and four WRKY genes, which were differentially expressed and may play an important role in browning process of fresh-cut luffa fruits. We also observed the differential expression of other Unigenes, such as Unigene0019122, Unigene0003760, and Unigene0031055, which were predicted to encode an ethylene-responsive transcription factor ERF109-like protein, an abscisic acid receptor PYL8-like protein, and gibberellin 2-oxidase, respectively ( Table 2). These Unigenes may contribute to the browning of fruits and vegetables [40,75,76]. However, their potential effects on the browning of fresh-cut luffa fruits will need to be confirmed in future studies.

Conclusions
We herein describe the first use of Illumina sequencing technology to investigate the changes in the transcriptome during the browning of fresh-cut luffa fruit slices. A total of 27,301 differentially expressed Unigenes were detected in four sequenced libraries. Additionally, 11 browningrelated genes from five gene families (i.e., PPO, PAL, POD, CAT, and SOD) as well as four WRKY transcription factors were observed to be differentially regulated in fresh-cut luffa fruits. These genetic resources and putative signaling pathways related to luffa defense responses against browning may be useful for future molecular studies of L. cylindrica.
Supporting information S1