Comparative transcriptomic analysis reveals the cold acclimation during chilling stress in sensitive and resistant passion fruit (Passiflora edulis) cultivars

Plant materials

HJG was introduced in the Bannahuangguo of Xishuangbanna Botanical Garden in Yunnan, and it is a cold-sensitive accession. TN1 comes from Taiwan, and it is a cold-resistant purple passion fruit. The cutting seedling heights ranged from 29 to 38 cm, and the seedlings were transplanted in Nanning experimental field (Guangxi, China, 22.85°N, 108.26°E) on May 25, 2019. The first sampling time was November 25, 2019, at 10 a.m., and the temperature was 25 °C. The second sampling time was January 18, 2020, at 10 a.m., and the temperature was 7 °C. The fresh leaves of passion fruit were snap frozen in liquid nitrogen, and then stored in −80 °C freezer. Each sample had three biological replicates. Under NT condition, the three biological replicates of HJG, HJGA1, HJGA2 and HJGA3, were recorded as A1; the three biological replicates of TN1, TN1A1, TN1A2 and TN1A3 were recorded as A2. Under CS condition, three biological replicates of HJG, HJGB1, HJGB2 and HJGB3, were denoted as B1; three biological replicates of TN1, TN1B1, TN1B2 and TN1B3, were denoted as B2.

RNA extraction, sequencing, assembly and annotation

Total RNA was extracted with RNAprep Pure kit (Tiangen, Beijing, China) according to the manufacturer’s instructions. Nanodrop2000 (Shimadzu, Japan) was used to detect the concentration and purity of the extracted RNA. Agarose gel electrophoresis was used to detect the integrity of the RNA, and Agilent 2100 (Agilent, America) was used to determine the RIN value. A single library requires 1 µg of RNA, with a concentration of ≥ 50 ng/µL, and OD260/280 between 1.8 and 2.2. Magnetic beads (Invitrogen, America) with Oligo (dT) was used to pair with the 3′ poly A tail of eukaryotic mRNA, thus isolating mRNA from total RNA. Subsequent, reverse synthesis of cDNA was performed. These libraries above were sequenced using an Illumina NovaSeq 6000 sequencer (Illumina Inc., USA) and 150-bp paired-end reads were generated. In order to ensure the accuracy of subsequent analysis, the original sequencing data were filtered first to obtain clean data.

We used Trinity v2.8.6 (Haas et al., 2013) to splice the transcript fragments to obtain transcripts, and then used CD-HIT to cluster the transcript sequences to remove redundant sequences and get all the unigene sequence sets for the subsequent analysis. Bowtie 2 (Salzberg & Langmead, 2012) was used to align the sequencing data to the reconstructed unigene sequence set, and the alignment file was mainly used for subsequent unigene quantification and differential expression analysis. The unigene sequences were compared with the NCBI non-redundant protein sequences (NR), Swiss-Prot, TrEMBL, KEGG, GO, Pfam, and EuKaryotic Orthologous Groups (KOG) databases using Basic Local Alignment Search Tool (BLAST). Finally, HMMER3 (Finn, Clements & Eddy, 2011) was used to align the amino acid sequence of unigene with the Pfam database to obtain the annotation information of unigenes.

Enrichment analysis of DEGs

The read counts and transcripts per million reads (TPM ) were calculated using RSEM (Li & Dewey, 2011) and bowtie2 (Salzberg & Langmead, 2012). The DEGs were identified through the software packages of Bioconductor 3.11-DESeq2 (Love, Huber & Anders, 2014). The screening threshold is false discovery rate (FDR ) <0.05, and log2 fold change (FC (condition 2/condition 1) for a gene) >1 or log2FC <−1. The DEGs were classified, and GO and KEGG enrichment analysis were subsequently performed.

Weighted gene co-expression network analysis

We followed these steps below for weighted gene co-expression network analysis (WGCNA): (i) screening DEGs for WGCNA cluster analysis; (ii) calling the R package to cluster the DEGs; (iii) calling ggplot2 in the R package to draw the clustering heat map and histogram of each module; (iv) using the top GO terms to perform GO enrichment analysis on each module; (v) calling Fisher-test function in R for KEGG enrichment analysis; and (vi) using Cytoscape3.8.0 (Su et al., 2014) to draw network diagram. A signed network was constructed using the blockwiseModules function, with the following parameters: power = 14, minModuleSize = 30, mergeCutHeight = 0.25, corType = pearson. When co-expressed genes are defined according to the above-mentioned standards, each gene is assigned a module number and corresponding module color; otherwise the ‘gray’ module was used.

Validation of the CS-related genes using RT-qPCR

We selected 11 genes related to plant hormone signaling, fatty acid metabolism and plant-pathogen interaction using GO and KEGG databases, and 4 hub genes in WGCNA for validation. The primers were designed using Primer3 (Table S1). Using HIS as the reference gene (Liu et al., 2017A), RT-qPCR was used to analyze the expression levels of 15 genes in B1 and B2.

The identical RNA samples as RNA-seq experiments were used for RT-qPCR analysis. The detailed experimental method refers to Wu et al. (2020). The relative gene expression level was calculated by reference to the 2^−ΔΔCt method (Livak & Schmittgen, 2001). All unigenes expression analysis were performed in triplicates. The values represented arithmetic averages of three replicates, and the data were expressed as a mean plus and minus standard deviation (mean ± SD).

Statistical analysis

CASAVA was used for base calling. Subsequently, we used SeqPrep for quality control of raw sequencing data. Pearson correlation coefficient is used to measure the correlation between samples. The package heatmap of R was used to prepare the correlation between samples and DEGs expression pattern clustering. Data of RT-qPCR was analyzed using Excel 2007. The figures were prepared using Origin 9.65.

Quality control and assembly of passion fruit transcriptome sequences

To compare gene expression profiles of the two passion fruit cultivars under NT and CS, transcriptome sequencing and analysis were performed. After decontamination and adaptor removal, 533,935,574 raw reads were obtained from 12 samples (NCBI accession number: PRJNA634206), a total of 80.09 Gb clean reads and 6.67 Gb per sample. The Q30 base percentage was 93.22% and GC content was 44.64% (Table 1).

The clean reads were assembled into transcripts using the Trinity in paired-end method, and 211,874 transcripts were obtained (https://doi.org/10.6084/m9.figshare.13489869.v2). The CD-HIT was then used to cluster the transcripts, which yielded 47,353 unigenes (https://doi.org/10.6084/m9.figshare.13489863.v4) with a mean length of 1,211 bp, N50 length of 2,368 bp, and N90 of 450 bp. Afterwards, Bowtie2 was used to align the sequences of each sample to the unigene sequence set, with an average alignment ratio of 77.89% (Table 2).

Unigene function annotation

The assembled unigenes were annotated to databases including the NR, Swiss-Prot, TrEMBL, KEGG, GO, Pfam, and KOG, to which 97.92%, 70.40%, 97.82%, 33.97%, 36.16%, 61.43%, and 48.92% of unigenes were mapped, respectively. A total of 47,353 unigenes acquired annotation information (Table 3). The number of annotated unigenes in NR and TrEMBL was the largest, which were 46,369 and 46,323, respectively.

In the GO database, 17,123 unigenes were annotated and matched to three major categories: biological process (BP), cellular component (CC) and molecular function (MF). Enriched BP terms were mainly about “metabolic process” (4,350) and “cellular process” (2,191). Enriched CC terms were mainly about “membrane part” (1,270) and “cell part” (890). Enriched MF terms were mainly about “binding” (7,367) and “catalytic activity” (5,715) (Fig. 1A).

In KOG database, 23,164 unigenes were annotated, which were clustered into 25 categories. The unigenes were mainly about “signal transduction mechanisms” (2,439) and “posttranslational modification, protein turnover, chaperones” (2,138) (Fig. 1B).

In KEGG database, 16,086 unigenes were annotated. According to the functions, these unigenes were enriched in 9 pathways. The enriched pathways were mainly about “metabolism” (10,045) and “organismal systems” (4,505) (Fig. 1C).

Comparative analysis of DEGs in two cultivars at two stages

In order to gain insights on the adaptation of passion fruit to CS, the TPM method was used to analyze the gene expression levels in the two stages (Fig. S1). The correlation coefficient between the three biological replicates was 0.87 in HJGA, 0.98 in TN1A, 0.96 in HJGB, 0.99 in TN1B, and the average correlation coefficient value was 0.95 (Fig. S2 ), indicating that the reproducibility of this study was good and the experimental results were reliable.

The software package DESeq2 was used to perform differential expression analysis of unigenes. There were 3,248 and 4,340 DEGs at two stages, respectively. After, analysis of the DEGs for the two stages, we found that the DEGs between HJG and TN1 were increased by 33.6% under CS condition (1,092), and 87.5% (955) were up-regulated (Fig. 2A).

Cluster analysis of gene expression can intuitively reflect the level of gene expression and expression patterns in multiple samples. We used the DEGs to perform cluster analysis on A1 vs A2 (Fig. 2B) and B1 vs B2 (Fig. 2C). The results showed that the difference between the three biological replicates of each group was small, which again confirmed the rationality of sample selection.

GO and KEGG pathway enrichment analysis of DEGs

There were 1,182 up-regulated unigenes, and 2,066 down-regulated unigenes at stage A; and there were 2,137 up-regulated unigenes and 2,203 down-regulated unigenes at stage B.

GO enrichment analysis indicated that “metabolic process” (542), “oxidation–reduction process” (156), “protein phosphorylation” (92), “carbohydrate metabolic process”(73), “organic substance catabolic process” (40), and “catabolic process” (40), “extracellular region” (10), “apoplast” (8), “cell wall” (8), “external encapsulating structure” (8), “catalytic activity” (634), “transferase activity” (228), “oxidoreductase activity” (167), “metal ion binding” (146), “cation binding” (146), and “transition metal ion binding” were enriched at stage A (110) (Table S2). But “oxidation–reduction process” (187), “phosphate-containing compound metabolic process” (171), “phosphorus metabolic process” (171), “macromolecule modification” (170), “cellular protein modification process” (169), and “protein modification process” (169), “membrane” (213), “intrinsic component of membrane” (99), and “integral component of membrane” (97),“catalytic activity” (837), “transferase activity” (326), “cation binding” (202), “metal ion binding” (201), “oxidoreductase activity” (198), “phosphotransferase activity”, and “alcohol group as acceptor” (165) were enriched at stage B (Table S3).

The GO terms in A (P > 0.05) were compared to B (P ≤ 0.05), and the unigenes were mainly about “protein phosphorylation” (GO:0006468 , 61) “phosphorylation” (GO:0016310, 61), “response to stimulus” (35), “lipid metabolic process” (19), “response to chemical” (13), “membrane” (73), “intrinsic component of membrane” (29), “integral component of membrane” (28) “catalytic activity”, “acting on a protein” (77), “transferase activity”, “transferring phosphorus-containing groups” (70), “kinase activity” (67), “phosphotransferase activity”, “alcohol group as acceptor” (67), and “protein kinase activity” (62) (Table S4).

The KEGG pathway enrichment analysis can reveal the main metabolic pathways and signal transduction pathways in which the DEGs were involved, and the prevailing pathways were as follows: “ribosome” (42), “carbon metabolism” (39), “biosynthesis of amino acids” (30), “starch and sucrose metabolism” (21), “and cysteine and methionine metabolism” (20) at stage A; “plant hormone signal transduction” (31), “plant-pathogen interaction” (27), “fatty acid metabolism” (21), “cysteine and methionine metabolism” (20) (Fig. 3A). The KEGG pathway in A (P > 0.05) were compared to B (P ≤ 0.05), and the unigenes were mainly about “plant-pathogen interaction” (17), “plant hormone signal transduction” (14), and “fatty acid metabolism” (8) (Fig. 3B).

WGCNA analysis

After background correction and normalization of the unigenes expression data, we filtered out the abnormal and minor changed unigenes. Finally, we obtained 12,471 highly expressed unigenes (Table S5). In this study, when the soft threshold was 16 (Fig. S3), the gene topology matrix expression network was closest to the scale-free distribution. A gene cluster tree was constructed based on the correlation between genes, and each branch corresponded to a cluster of gene sets with highly correlated expression levels (Fig. S4A).

According to the standard of mixed dynamic shear, the gene modules were classified and the eigenvector of each module was calculated. The modules close to each other were merged, and 8 co-expression modules were obtained (Fig. S4B). Each module used different colors to represent the clustered genes. The turquoise module had the most clustered genes (4,171), the red module contained the fewest (81), and the grey module contained the unigenes that couldn’t be included in any module.

The DEGs were used to draw the heat map of each module in the four sample groups. The brown and yellow modules showed less changes in differential unigenes between the early and late HJG, but showed larger changes in differential unigenes between early and late TN1 (Fig. 4), which is consistent with the chilling resistance feature of TN1. Therefore, we selected the unigenes of these two modules for in-depth GO and KEGG pathway analysis.

In the brown module, the significant GO terms were “cellular macromolecule metabolic process”, “phosphate-containing compound metabolic process”, “phosphorus metabolic process”, “protein phosphorylation”, “stimulus”, “transferase complex”, “riboflavin synthase complex”, “photosystem I reaction center”, “photosystem I” “binding”, “metal ion binding”, “cation binding”, “phosphotransferase activity”, “alcohol group as acceptor”, “kinase activity” (Table S6). In the KEGG pathway analysis, the prevailing pathways were “plant hormone signal transduction”, “MAPK signaling pathway”, and “starch and sucrose metabolism” (Fig. 4A).

In the yellow module, the significant GO terms were “cellular process”, “macromolecule modification”, “phosphorus metabolic process”, “cellular protein modification process”, “protein modification process”, “cell periphery”, “photosystem”, “photosynthetic membrane”, “thylakoid”, “extracellular region”, “3-deoxy-7-phosphoheptulonate synthase activity”, “alkylbase DNA N-glycosylase activity”, “DNA-3-methyladenine glycosylase activity”, “DNA N-glycosylase activity”, and “method adenosyltransferase activity” (Table S7). The significantly enriched pathways included “biosynthesis of amino acids”, “plant hormone signal transduction”, “ABC transporters”, “starch and sucrose metabolism”, “folate biosynthesis” and “other pathways” that might be related to CS (Fig. 5B).

The constructed network was visualized with the Cytoscape in the brown and yellow modules, and got 19 hub unigenes which mainly related to “MAPK signaling pathway”, “plant hormone signal transduction”, “starch and sucrose metabolism”, “fatty acid biosynthesis” and “photosynthesis in the brown module” (Fig. 6A). In the yellow module, 13 hub unigenes were mainly related to “plant hormone signal transduction”, “MAPK signaling pathway”, “starch and sucrose metabolism”, and “fatty acid degradation” (Fig. 6B).

Validation of gene expression changes during chilling acclimation

We used the RT-qPCR method to validate the expression levels of 15 unigenes. The results showed that the RT-qPCR expression patterns of the 15 unigenes were consistent with RNA-seq analysis (Fig. 7, Table S8). RT-qPCR analysis showed that the 9 unigenes were ≥2 or ≤0.5 fold-change. Comparison with B1 and B2, TPM value of 12 unigenes were ≥2 or ≤0.5 fold-change. The results showed that seventy-five percent DEGs could be validated using RT-qPCR, and DEGs analysis were highly reliable.