De Novo Transcriptome Sequencing in Passiflora edulis Sims to Identify Genes and Signaling Pathways Involved in Cold Tolerance

The passion fruit (Passiflora edulis Sims), also known as the purple granadilla, is widely cultivated as the new darling of the fruit market throughout southern China. This exotic and perennial climber is adapted to warm and humid climates, and thus is generally intolerant of cold. There is limited information about gene regulation and signaling pathways related to the cold stress response in this species. In this study, two transcriptome libraries (KEDU_AP vs. GX_AP) were constructed from the aerial parts of cold-tolerant and cold-susceptible varieties of P. edulis, respectively. Overall, 126,284,018 clean reads were obtained, and 86,880 unigenes with a mean size of 1449 bp were assembled. Of these, there were 64,067 (73.74%) unigenes with significant similarity to publicly available plant protein sequences. Expression profiles were generated, and 3045 genes were found to be significantly differentially expressed between the KEDU_AP and GX_AP libraries, including 1075 (35.3%) up-regulated and 1970 (64.7%) down-regulated. These included 36 genes in enriched pathways of plant hormone signal transduction, and 56 genes encoding putative transcription factors. Six genes involved in the ICE1–CBF–COR pathway were induced in the cold-tolerant variety, and their expression levels were further verified using quantitative real-time PCR. This report is the first to identify genes and signaling pathways involved in cold tolerance using high-throughput transcriptome sequencing in P. edulis. These findings may provide useful insights into the molecular mechanisms regulating cold tolerance and genetic breeding in Passiflora spp.


Introduction
Passiflora is the largest genus of the Passifloraceae family, with more than 500 species [1]. Passiflora species are distributed throughout Latin America, and Brazil and Colombia are the centers of diversity for this genus [2], and many of these species are widely cultivated for their edible fruit, medicinal efficacy, and ornamental properties. In the early 20th century, Passiflora edulis Sims as an edible fruit was introduced to China, mainly in Taiwan, Guangdong, Guangxi and Fujian. P. edulis is known for its taste, is used in Brazilian traditional folk medicine and is included in pharmacopoeias of several countries [3,4]. Leaf extracts of P. edulis are considered to treat alcoholism, anxiety and insomnia [5]. The flower has been used for the treatment of cough and bronchitis, and the seed oil as a lubricant and massage oil [6].

Transcriptome Assembly and Function Annotation
Transcriptome assembly was achieved using Trinity version r20140413p1 [23] based on the left.fq and right.fq, with the min_kmer_cov set as 2 by default and all other parameters set as their defaults. For function annotation, the longest transcript of each gene was defined as the 'unigene'. All assembled unigenes were BLASTed in Nr, Nt, Pfam, KOG/COG, Swiss-Prot, KEGG ortholog database (KO) and Gene ontology (GO) databases using BLAST2GO of version 2.5 with a cut-off E-value of 10 −6 [24].

Differential Expression Analysis
Gene expression level of all samples was estimated by mapping clean reads to the Trinity transcripts assembly using RSEM version 1.2.15 [25], with the bowtie2 parameter set at mismatch 0. Then, read counts for each gene were obtained from the mapping results. Prior to Differentially expressed genes (DEG) analysis, the read counts were normalized using the edgeR program package with the Trimmed Mean of M-values method [26,27]. The DEGseq R package was used to analyze differential expression of two P. edulis samples. The p-value was adjusted using q-value [28]. The significant differential expression threshold was set as q-value < 0.005 and |log 2 (foldchange)| > 1 [29]. The identified DEGs were used for GO enrichment analyses, which were performed using the GOseq R package (version 1.10.0), based on the Wallenius non-central hypergeometric distribution [30]. Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis was performed using KOBAS version 2.0.12 [31].

Quantitative Real-Time PCR (qRT-PCR)
Total RNA was extracted from various samples with the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) and treated with RNAse-free DNase I (Ambion, Austin, Texas, USA). Of the DNase-treated RNA, 1 µg was reverse transcribed using random hexamer primers. The resulting cDNA was diluted, and 2 µL of the diluted cDNA was used. Specific primers (Tm, 58-61 • C) were designed to generate PCR products of 70-150 bp (Table S1). The qRT-PCR was performed on an ABI ViiA 7 Real-time PCR platform. FastStart Universal SYBR Green Master (Rox) was used for qRT-PCR assays according to the manufacturer's protocol. For each sample, three replicates were performed in a final volume of 20 µL containing 10 µL of SYBR Premix Ex Taq (2×), 0.4 µL of 50 × ROX Reference Dye II, 0.4 µL (10 µM) of each primer, 2 µL of cDNA, and 6.8 µL of dH2O. Thermo-cycling conditions were as follows: initial denaturation at 95 • C for 30 s, followed by 40 cycles of 95 • C for 5 s, and 60 • C annealing and extension for 34 s. All reactions were performed in triplicate. The specificity of the PCR amplification procedures was checked with a heat dissociation protocol after the final cycle of the PCR to ensure that each amplicon was a single product. Relative expression was calculated as the difference in delta-Ct between the target gene and the internal control, histone H3.3 (HIS) gene.

Transcriptome Sequencing and De Novo Assembly
In total, there were 60,881,198 raw reads generated from KEDU_AP and 71,719,162 from GX_AP. The sequencing raw data were submitted to the Short Reads Archive (SRA) database under the accession number SRP106510. Of the raw reads from KEDU_AP, more than 95.87% bases had a q value ≥ 20, and for GX_AP 95.15%. The GC-contents were 44.84% and 45.20% for KEDU_AP and GX_AP, respectively. After removing low-quality reads, 57,840,324 clean reads from KEDU_AP and 68,443,694 from GX_AP were obtained. These were used for de novo assembly.
The Trinity software generated 127,474 transcripts with an average length of 1077 bp and an N50 of 2057 bp (

Function Annotation and Classification
All the 86,880 assembled unigenes were annotated to the seven databases using the BLAST algorithm (Table 2). In total, 64,067 unigenes were annotated, accounting for 73.74% (Table 2). There were 10,457 (12.03%) unigenes successfully annotated in all seven databases. Analyses showed that 60,028 (69.09%) unigenes had high homology with sequences in the Nr database and 46,093 (53.05%) in the Nt database. The details of other database proportions are shown in Table 2. There were 46,754 unigenes divided into three functional GO categories: biological process (BP), cellular component (CC) and molecular function (MF) (Figure 1). In the BP category, these matched unigenes were annotated to 25 GO terms, with the three top terms being 'cellular process' (27,925), 'metabolic process' (26,692) and 'single-organism process' (21,336). For CC and MF, these unigenes were clustered into 20 and 10 GO terms, respectively, with 'cell' (16,101) and 'binding' (27,248) as the largest subcategories. There were 17,742 unigenes divided into 26 groups for KOG  (post-translational modification, protein turnover and chaperon) and then J (translation, ribosomal structure and biogenesis).  A total of 18,366 unigenes were assigned to 19 metabolic pathways in the KEGG database ( Figure 3). These 19 KEGG pathways were clustered into five branches: cellular processes (A) of 1164 unigenes, environmental information processing (B) of 999 unigenes, genetic information processing (C) of 4935 unigenes, metabolism (D) of 10,430 unigenes and organismal systems (E) of 838 unigenes. The metabolic pathways with the largest number of unigenes were 'carbohydrate metabolism' (2101), followed by 'translation' (2087) and 'folding, sorting and degradation' (1528).  (post-translational modification, protein turnover and chaperon) and then J (translation, ribosomal structure and biogenesis).

DEG Identification, GO and KEGG Enrichment Analysis
A total of 80,810 unigenes (Fragments per kilobase of exon per million fragments mapped (FPKM) > 0.3) in two groups were identified, with 57,213 unigenes in common. There were 3045 significant DEGs identified between the KEDU_AP and GX_AP libraries. For these DEGs, if log2Foldchange >1, the DEG was considered as up-regulated but if log2Fold change <−1, it was considered as down-regulated. The 3045 DEGs included 1075 up-regulated and 1970 down-regulated DEGs ( Figure 4). Between the two libraries of enriched GO terms, the analysis showed that 'catalytic activity' and 'single-organism metabolic process' had the highest degree of enrichment ( Figure 5). In the KEGG enrichment analysis (Figure 6), the most significantly enriched pathway was 'flavonoid biosynthesis', which mainly consisted of up-regulated DEGs. The second most significantly enriched pathway was 'phenylpropanoid biosynthesis', which also mainly consisted of up-regulated DEGs. Down-regulated DEGs were dominant in the 'photosynthesis-antenna proteins' and 'amino sugar and nucleotide sugar metabolism' (Figure S1). The greatest numbers of DEGs were involved in 'starch and sucrose metabolism', followed by 'phenylpropanoid biosynthesis'.  Between the two libraries of enriched GO terms, the analysis showed that 'catalytic activity' and 'single-organism metabolic process' had the highest degree of enrichment ( Figure 5). In the KEGG enrichment analysis (Figure 6), the most significantly enriched pathway was 'flavonoid biosynthesis', which mainly consisted of up-regulated DEGs. The second most significantly enriched pathway was 'phenylpropanoid biosynthesis', which also mainly consisted of up-regulated DEGs. Down-regulated DEGs were dominant in the 'photosynthesis-antenna proteins' and 'amino sugar and nucleotide sugar metabolism' (Figure S1). The greatest numbers of DEGs were involved in 'starch and sucrose metabolism', followed by 'phenylpropanoid biosynthesis'.

DEGs Involved in the Plant Hormone Signal Transduction Pathways
Plant hormones play an important role in plant cold tolerance. A total of 36 DEGs and 677 background unigenes were in the enriched pathway of plant hormone signal transduction ( Figure  7). Plant hormone signal transduction contained eight sub-paths. In the auxin signal transduction sub-pathway, there were two up-regulated and nine down-regulated genes. In the cytokinin signal   Between the two libraries of enriched GO terms, the analysis showed that 'catalytic activity' and 'single-organism metabolic process' had the highest degree of enrichment ( Figure 5). In the KEGG enrichment analysis (Figure 6), the most significantly enriched pathway was 'flavonoid biosynthesis', which mainly consisted of up-regulated DEGs. The second most significantly enriched pathway was 'phenylpropanoid biosynthesis', which also mainly consisted of up-regulated DEGs. Down-regulated DEGs were dominant in the 'photosynthesis-antenna proteins' and 'amino sugar and nucleotide sugar metabolism' (Figure S1). The greatest numbers of DEGs were involved in 'starch and sucrose metabolism', followed by 'phenylpropanoid biosynthesis'.

DEGs Involved in the Plant Hormone Signal Transduction Pathways
Plant hormones play an important role in plant cold tolerance. A total of 36 DEGs and 677 background unigenes were in the enriched pathway of plant hormone signal transduction ( Figure  7). Plant hormone signal transduction contained eight sub-paths. In the auxin signal transduction sub-pathway, there were two up-regulated and nine down-regulated genes. In the cytokinin signal

DEGs Involved in the Plant Hormone Signal Transduction Pathways
Plant hormones play an important role in plant cold tolerance. A total of 36 DEGs and 677 background unigenes were in the enriched pathway of plant hormone signal transduction (Figure 7). Plant hormone signal transduction contained eight sub-paths. In the auxin signal transduction sub-pathway, there were two up-regulated and nine down-regulated genes. In the cytokinin signal transduction sub-pathway, there were three up-regulated and four down-regulated genes. There was just one up-regulated gene in the gibberellin signal transduction sub-pathway. Only one down-regulated gene was related to jasmonic acid signal transduction. There were up-regulated and down-regulated genes in the brassinolide and salicylic acid signal transduction pathways, which of them all contained three up-regulated genes and one down-regulated gene. Ethylene and ABA have vital functions in cold-stress signaling. In the ethylene signal transduction sub-pathway, there were two up-regulated genes and one down-regulated gene. In the ABA signal transduction sub-pathway, there were two up-regulated genes and one down-regulated gene for PYR/PYL, two down-regulated genes for PP2C, and one up-regulated and two down-regulated genes for SnRK2. transduction sub-pathway, there were three up-regulated and four down-regulated genes. There was just one up-regulated gene in the gibberellin signal transduction sub-pathway. Only one down-regulated gene was related to jasmonic acid signal transduction. There were up-regulated and down-regulated genes in the brassinolide and salicylic acid signal transduction pathways, which of them all contained three up-regulated genes and one down-regulated gene. Ethylene and ABA have vital functions in cold-stress signaling. In the ethylene signal transduction sub-pathway, there were two up-regulated genes and one down-regulated gene. In the ABA signal transduction sub-pathway, there were two up-regulated genes and one down-regulated gene for PYR/PYL, two down-regulated genes for PP2C, and one up-regulated and two down-regulated genes for SnRK2.

Verification by qRT-PCR
In the ICE1-CBF-COR pathway, we identified one ICE1, one COR, two ICE1-like and two CBF genes. These genes were up-regulated in cold-tolerant variety 'Pingtang 1' (KEDU_AP) by qRT-PCR verification (Figure 8). To verify the reliability of the transcriptome data, 11 DEGs were randomly selected and examined using qRT-PCR at the transcriptional level ( Figure S2). The expression These results indicate that our transcriptomic analysis was highly reproducible and reliable.
In the ICE1-CBF-COR pathway, we identified one ICE1, one COR, two ICE1-like and two CBF genes. These genes were up-regulated in cold-tolerant variety 'Pingtang 1' (KEDU_AP) by qRT-PCR verification (Figure 8). To verify the reliability of the transcriptome data, 11 DEGs were randomly selected and examined using qRT-PCR at the transcriptional level ( Figure S2). The expression patterns of the 11 DEGs were consistent with the transcriptome data (R 2 = 0.81321, p-value = 0.002334). These results indicate that our transcriptomic analysis was highly reproducible and reliable.

Discussion
As a key messenger, the cytosolic Ca 2+ occupies an important role in response to cold stress [32,33]. Research shows that cold-induced Ca 2+ influx is positively correlated with accumulation of cold-induced transcripts in Arabidopsis [34] and alfalfa (Medicago sativa L.) [35]. Cold stress increases the cytosolic Ca 2+ level, and this is then sensed by Ca 2+ sensor proteins such as calmodulin (CaM) and calcineurin B-like proteins (CBLs) [36,37]. The increased levels of Ca 2+ in plant cells affect the expression level of CBF and COR genes in the cold signaling pathway [7,38]. In this study, 332 unigenes were annotated as CaMs and CBLs, and included 14 DEGs. Among these DEGs, seven (four CaM and three CBL) genes were up-regulated in KEDU_AP, and seven (six CaM and one CBL) were down-regulated ( Figure S3 and Table 3). Plant protein kinases, such as Mitogen-activated protein kinases (MAPKs), play a central role in cellular signaling. MKK2 is a member of the MAPK family, and is upstream of MPK4 and MPK6, which are activated by low temperatures [14]. Previous studies showed that the influx of calcium is

Discussion
As a key messenger, the cytosolic Ca 2+ occupies an important role in response to cold stress [32,33]. Research shows that cold-induced Ca 2+ influx is positively correlated with accumulation of cold-induced transcripts in Arabidopsis [34] and alfalfa (Medicago sativa L.) [35]. Cold stress increases the cytosolic Ca 2+ level, and this is then sensed by Ca 2+ sensor proteins such as calmodulin (CaM) and calcineurin B-like proteins (CBLs) [36,37]. The increased levels of Ca 2+ in plant cells affect the expression level of CBF and COR genes in the cold signaling pathway [7,38]. In this study, 332 unigenes were annotated as CaMs and CBLs, and included 14 DEGs. Among these DEGs, seven (four CaM and three CBL) genes were up-regulated in KEDU_AP, and seven (six CaM and one CBL) were down-regulated ( Figure S3 and Table 3). Plant protein kinases, such as Mitogen-activated protein kinases (MAPKs), play a central role in cellular signaling. MKK2 is a member of the MAPK family, and is upstream of MPK4 and MPK6, which are activated by low temperatures [14]. Previous studies showed that the influx of calcium is involved in cold-stress regulation of MAPKs [39]. The receptor-like kinase CRLK1 may associate calcium signaling with the MAPK cascade by binding to calcium, and CaM interacts with MEK kinase 1 (MEKK1) [40] (Figure 9). We identified 93 MAPK unigenes, which included 11 DEGs: nine up-regulated and two down-regulated ( Figure S3 and Table 3). involved in cold-stress regulation of MAPKs [39]. The receptor-like kinase CRLK1 may associate calcium signaling with the MAPK cascade by binding to calcium, and CaM interacts with MEK kinase 1 (MEKK1) [40] (Figure 9). We identified 93 MAPK unigenes, which included 11 DEGs: nine up-regulated and two down-regulated ( Figure S3 and Table 3). In this work, 56 differentially expressed TFs were identified, and were divided into six TF families: AP2/EREBP, WRKY, bZIP, MYB, NAC and B-ARR ( Figure S4 and Table 3). Most of them have been reported to be associated with cold tolerance in plants [41][42][43][44]. The TFs of CBFs belong to the AP2/EREBP family [45]. Three CBF genes are known in Arabidopsis, and are induced with plant exposure to cold, and regulate COR gene accumulation [46,47] (Figure 9). Overexpression of OsMYB4 could improve freezing tolerance in Arabidopsis by increasing COR gene expression, and OsMYBS3 and OsMYB2 enhanced cold resistance in rice [48][49][50]. We identified six genes in the ICE1-CBF-COR pathway of P. edulis, and their expression levels in the cold-tolerant variety 'Pingtang 1' were higher than that in cold-susceptible 'Purple Fragrance 1'.
In addition to the TFs, other important genes involved in plant cold tolerance include late-embryogenesis-abundant proteins (LEAs), antifreeze proteins (AFPs), superoxide dismutase (SOD) and proline [51][52][53]. LEAs are highly hydrophilic and provide protection for plants during cold stress. We also found four differentially expressed LEAs ( Figure S3 and Table 3), with a high expression level in KEDU_AP.
In this work, 56 differentially expressed TFs were identified, and were divided into six TF families: AP2/EREBP, WRKY, bZIP, MYB, NAC and B-ARR ( Figure S4 and Table 3). Most of them have been reported to be associated with cold tolerance in plants [41][42][43][44]. The TFs of CBFs belong to the AP2/EREBP family [45]. Three CBF genes are known in Arabidopsis, and are induced with plant exposure to cold, and regulate COR gene accumulation [46,47] (Figure 9). Overexpression of OsMYB4 could improve freezing tolerance in Arabidopsis by increasing COR gene expression, and OsMYBS3 and OsMYB2 enhanced cold resistance in rice [48][49][50]. We identified six genes in the ICE1-CBF-COR pathway of P. edulis, and their expression levels in the cold-tolerant variety 'Pingtang 1' were higher than that in cold-susceptible 'Purple Fragrance 1'.
In addition to the TFs, other important genes involved in plant cold tolerance include late-embryogenesis-abundant proteins (LEAs), antifreeze proteins (AFPs), superoxide dismutase (SOD) and proline [51][52][53]. LEAs are highly hydrophilic and provide protection for plants during cold stress. We also found four differentially expressed LEAs ( Figure S3 and Table 3), with a high expression level in KEDU_AP.
In this study, we identified candidate genes and signaling pathways associated with cold tolerance by transcriptome sequencing of cold-tolerant variety 'Pingtang 1' and cold-susceptible 'Purple Fragrance 1'. The heatmap of DEGs showed that many genes related to cold-tolerance had higher expression levels in the cold-tolerant variety, indicating that they played a role in the chilling stress response. This report is the first to identify genes and signaling pathways involved in cold tolerance using high-throughput transcriptome sequencing in P. edulis. These findings may provide useful insights into the molecular mechanisms regulating cold tolerance and genetic breeding in Passiflora spp.

Conclusions
P. edulis is mainly distributed in tropical and subtropical regions and has difficulty surviving at low temperatures. To elucidate the molecular mechanisms of cold tolerance, we collected the aboveground parts of cold-tolerant variety 'Pingtang 1' and cold-susceptible 'Purple Fragrance 1' and subjected them to high-throughput sequencing. In total, 26,284,018 clean reads were obtained, and 86,880 unigenes with a mean size of 1449 bp were assembled. There were 3045 significant differentially expressed genes (DEGs) identified between 'Pingtang 1' (KEDU_AP) and 'Purple Fragrance 1' (GX_AP). To provide reference for the cold-tolerance breeding of P. edulis, we screened many DEGs, constructed their expression profiles and analyzed their functions; some potentially vital cold-tolerance genes and transcription factors were identified, and the expression levels of the two varieties were compared and analyzed. As the first report on the high-throughput sequencing of cold-tolerant P. edulis, this study should provide novel insights into cold-tolerance genes for P. edulis and be a valuable molecular basis for study in Passiflora spp.