Two passion fruit cultivar utilized in this study: purple passion fruit (Passiflora edulis ‘Tainong No. 1’), and yellow passion fruit (Passiflora edulis ‘Qinmi 9’), were both collected from the demonstration orchard in the Sujia Village, located in the Wanhe Zhelmu Town, Guilin City of Guangxi region. From each cultivar, a total of 100 fruits was meticulously selected from the fruit trees, divided into four groups, and then randomly sampled. The fruit samples were then classified into two maturity categories: ripe and unripe. The peels of unripe purple passion fruit were designated as TPU, and those of ripe fruits as TPR. For the yellow passion fruit, the peels of unripe fruits were labeled as TYU, and those of ripe fruit as TYR. The fruit samples were then washed with tap water and air-dried. The pectin layer was removed from the fruit peels, followed by a thorough washing of the pericarp layer and its subsequent drying on a bench. Subsequently, the peel samples were flash-frozen using liquid nitrogen, pulverized, and stored at -80 °C for further analysis. The fruit peels from the yellow variety were used for comparison.

Transcriptome sequencing was conducted in the UW Genetics Ltd. (Shenzhen, China). The total RNA was extracted from the peel samples of the purple passion fruit (TPU and TPR) and yellow passion fruit (TYU and TYR) using an RNA extraction kit (Cheng et al., 2023). Subsequently, Total RNA concentration and purity were measured by spectrophotometry. RNA integrity was assessed using the Agilent 2100 Bioanalyzer. After passing quality control, paired-end sequencing was performed on the Illumina HiSeq platform.

The library was subjected to sequencing using a second-generation sequencing platform (Illumina HiSeq 2000). An initial evaluation of the data was then conducted using the FastQC bioinformatics tool (Ward et al., 2020). Subsequently, the results were integrated and presented in a visual format using a MultiQC software (Ewels et al., 2016). Non-genomic sequences were trimmed using NGSQC Toolkit and Trimmomatic softwares with default parameters based on the evaluation results (Bolger et al., 2014). The softwares were employed for strict control of raw data, including removals of non-genomic sequences, reads with connectors (adapters). The presence of N bases in excess of this threshold has been demonstrated to affect the subsequent comparison results. Furthermore, low-quality bases (Q ≤ 20) were removed. Next, PCR duplicates generated during library construction were removed using FastUniq software to obtain clean reads for the subsequent comparisons (Xu et al., 2012). After quality control steps, the clean reads were compared to the reference genome in order to identify the source genes that were transcribed by these fragments. The clean reads were then aligned aganist the Passiflora reference genome using HISAT2, which yielded information regarding the genomic location, gene features, and sequence characteristics that were specific to the sequenced sample.

The truly differentially expressed genes (DEGs) were addressed based on their biological variability due to considerable variation in gene expression among individuals. To facilitate a comprehensive analysis of the entire transcriptome, biological replicates were established. Pearson’s correlation coefficient (R) was employed to assess the correlation between biological replicates in the transcriptome sequencing data. Pearson’s correlation coefficient (R) was also utilized to quantify the degree of correlation between biological replicates (Li et al., 2015). As R approaches 1, the correlation between the replicates increases.

In the RNA-Seq, transcript abundance functioned as an indicator of gene expression levels. FPKM (fragments per kilobase of transcript per million mapped reads) was utilized to quantify transcript abundance. This approach quantifies the expression level of each gene as a fraction of the total length of the transcript in millions of mapped reads (Mizushima et al., 2020). Principal component analysis (PCA) was a statistical method that does not require supervision. PCA was used to identify patterns in the data. The principal component 1 (PC1) represented the principal component that captured most of the variance in the data matrix, while the principal component 2 (PC2) captured additional and distinct variance. The employment of PCA was instrumental in extracting the maximum amount of information from the original variables, thereby ensuring their independence (Liu et al., 2003).

The DESeq2 method was employed to analyze the differential expression among the samples studied, thereby identifying the DEG sets between groups. The unstandardized data were utilized as the input for the differential analysis. Subsequent to the differential analysis, the Benjamini-Hochberg method was employed to correct the P-values for multiple testing, thereby yielding the false discovery rate (FDR). The identification of differential genes was facilitated by establishing a criteria of [log₂ fold change] ≥ 1 and FDR < 0.05. Subsequently, hierarchical clustering analysis was performed on FPKM expression data following differential gene centralization and standardization, and a corresponding heatmap was generated (Love et al., 2014; Varet et al., 2016).

To investigate the molecular basis of anthocyanin accumulation in passion fruit peel, RNA-seq analysis was performed on purple and yellow varieties at two developmental stages. Samples were collected from mature and immature fruits, with three biological replicates for each group, totaling 12 groups. RNA-Seq was performed on the Illumina-HiSeq platform, generating 559 million raw reads. After filtering, 535 million clean reads with an average length of 150 bp were obtained. The raw sequence data was evaluated using the Fastqc, and the summary reports were generated using the MultiQC. Based on these evaluations, the data was further filtered or trimmed to remove reads including sequences with substandard adapter sequences, excessive N bases, low-quality bases, and repetitive sequences caused by non-random fragmentation. The filtered clean reads was then aligned to the passion fruit reference genome using HISAT2 ( Table 1 ). The alignment results indicated that the alignment rate for each group was higher than 73%, indicating high data quality, good mapping efficiency and no contamination occurred during sequencing. The transcriptome data obtained in this sequencing can be used for subsequent analysis.

The correlation analysis of biological replicate samples requires that the R² value between replicates should be at least 0.8 or higher. Statistical analysis (seen in Figure 1A ) showed that, the R² value among biological replicates in each group was greater than 0.9, which met the requirements of quality criteria and confirmed the reliability of the data for subsequent analysis. The FPKM density distribution of the 12 groups was compared ( Figure 1B ), and the results showed that there was a large overlap among samples, indicating that the overall gene expression levels of the samples were similar. The non-overlapping parts might be attributed to variations in the different varieties and maturity stages of the passion fruit, which caused changes in gene expression levels.

In the principal component analysis ( Figure 1C ), the scattered points corresponding to the four sample groups displayed a trend of mutual aggregation within each groups, suggesting good internal repeatability and high similarity among the sample data; In contrast, the scattered points between the groups exhibited a pattern of mutual dispersion, indicating distinct discrimination among the groups, thus suitable for downstream differential expression and functional analyses.

The R package DESeq2 was employed to analyze the differential expression between samples, resulting in the identification of differentially expressed gene sets across various group. After the differential analysis, the Benjamini-Hochberg method was applied to adjust the P-values for multiple hypothesis testing. By setting the criteria of |log₂Fold Change| ≥ 1 and FDR < 0.05 as significance threshhold, the differentially expressed genes (DEGs) were screened, and the differences among genes between different combinations were also compared.

The Venn diagram ( Figure 2A ) illustrates the overlap of DEGs among different comparison groups, providing critical insights into the transcriptional regulation network underlying passion fruit peel coloration. Comparative analysis between mature purple-skin and yellow-skin passion fruit (TPR_vs_TYR) identified 7,528 DEGs, with 3,504 upregulated and 4,024 downregulated genes. This substantial change suggests variety-specific metabolic pathways are most active in mature fruit. Addtionally, fewer DEGs (5,039) were detected between immature purple-skin and yellow-skin fruit (TPU_vs_TYU), comprising 2,272 upregulated and 2,767 downregulated genes, indicating variety differences are already established early in development. This suggests that transcriptomic divergence increases along with fruit ripening, potentially corresponding to anthocyanin accumulation. Furthermore, a comparison between the immature stage of purple-skin passion fruit and the mature stage (TPU_vs_TPR) revealed 3,976 differentially expressed genes, including 1,894 up-regulated genes and 2,082 down-regulated genes, likely playing crucial roles in pigment accumulation during maturation.

Hierarchical clustering analysis was performed using normalized FPKM values of differentially expressed genes, and the resulting heatmap ( Figure 2B ) revealed distinctive expression patterns across samples. Notably, TPU and TPR samples clustered more closely together, as did TYU and TYR samples, indicating that genotype (purple-skin versus yellow-skin) exerts a stronger influence on global gene expression patterns than developmental stage. This genotype-dependent clustering suggests that the fundamental transcriptional programs governing fruit development are largely conserved within each variety across maturation stages. In contrast, there were marked differences in the color patterns between TPU and TYU, as well as between TPR and TYR. These differences in gene expression patterns between purple-skin and yellow-skin genotypes likely reflect divergent regulation of pigment biosynthesis pathways, particularly those involved in anthocyanin accumulation. The hierarchical clustering results provide a foundation for identifying the key regulatory modules that differentiate anthocyanin-accumulating and non-accumulating passion fruit genotypes.

To elucidate the biological functions of DEGs, KEGG pathway enrichment analysis was implemented ( Figure 3 ). In the comparison of TPR_vs_TYR ( Figure 3A ), 2628 of the 7528 DEGs were successfully annotated, accounting for 34.91% of the total DEGs, and they were involved in 113 metabolic pathways of KEGG. The largest number of annotated genes was found in the category of metabolism, with the metabolic pathway pathways and the biosynthesis of secondary metabolites having the most annotated DEGs, totaling 1215 and 724, respectively.

For the comparison of TPU_vs_TPR ( Figure 3B ), KEGG annotation was conducted on 3,976 DEGs, revealing that 1,478 of these genes were annotated, which accounted for 37.17% of the total DEGs, and they were involved in 72 KEGG metabolic pathways. Again, the majority of annotated genes were in the metabolism category, with 735 DEGs related to metabolic pathways and 455 to the biosynthesis of secondary metabolites.

Additionally, KEGG annotation was performed on 5,039 DEGs from the TPU_vs_TYU comparison ( Figure 3C ), where 1,832 of the DEGs were annotated, representing 36.36% of the total DEGs, and these genes were involved in 89 KEGG metabolic pathways. As with the previous comparisons, the largest number of annotated genes was in the metabolism category, with 884 DEGs associated with metabolic pathways and 536 with the biosynthesis of secondary metabolites.

After the genes have been annotated, the number of DEGs in each pathway was calculated, and a KEGG enrichment scatter plot of DEGs was drawn ( Figure 4 ). The KEGG enrichment scatter plot for the DEGs from the TPR_vs_TYR compariso ( Figure 4A ) display the 20 most significantly enriched pathway entries, with the DEGs annotated to the biosynthetic pathways of secondary metabolites being the most numerous, far exceeding other pathways. Therefore, the Rich factor of this pathway is the smallest and the q-value is the lowest, indicating significant enrichment.

Similarly, the KEGG enrichment scatter plot for the DEGs from the TPU_vs_TPR comparison ( Figure 4B ) also shows the 20 most significantly enriched pathway entries, where the DEGs associated with the biosynthetic pathways of secondary metabolites are again the most numerous. This pathway exhibits the smallest Rich factor and the lowest q-value, confirming its significant enrichment.

Furthermore, the KEGG enrichment scatter plot for the DEGs from the TPU_vs_TYU comparison ( Figure 4C ) presents the 20 most significantly enriched pathway entries, with the DEGs annotated to the biosynthetic pathways of secondary metabolites remaining the most numerous. As with the previous comparisons, this pathway has the smallest Rich factor and the lowest q-value, indicating significant enrichment.

In the GO term enrichment analysis, 7528 DEGs from the TPR_vs_TYR comparison were categorized into 58 subclasses across the three broad categories ( Figure 5A ). Within the BPs category, the categories were “Cellular Processes” (45.1%), “Metabolic Processes” (38.2%), and “Response to Stimulus” (23%). The subclasses annotated under the CCs category included “Cells” (55.6%), “Organelles” (43%), and “Membranes” (29.4%). In the MFs category, the annotations were predominantly “Binding” (45%) and “Catalytic Activity” (39.1%).

In the TPU_vs_TPR comparison, 3976 DEGs were categorized into 56 subclasses across the three broad categories ( Figure 5B ). For the BPs category, the distribution was as follows: “Cellular Processes” (47.2%), “Metabolic Processes” (40.1%), and “Response to Stimulus” (26.8%). In the CCs category, the DEGs were annotated to “Cells” (58.6%), “Organelles” (43.4%), and “Membranes” (32.8%). For the MFs category, the predominant subclasses were “Binding” (47.5%) and “Catalytic Activity” (41.5%).

In the TPU_vs_TYU comparison, 5039 DEGs were classified into 57 subclasses within the three broad categories ( Figure 5C ). Within the category of BPs, the most prevalent subclasses were “Cellular Processes” (45.1%), “Metabolic Processes” (38.3%), and “Response to Stimulus” (25.1%). For the CCs category, the annotations were “Cells” (55.6%), “Organelles” (42.6%), and “Membranes” (30.1%). Finally, the MFs category was included “Binding” (45.6%) and “Catalytic Activity” (41.5%).

In the TPR_vs_TYR comparison, 4321 DEGs were annotated to KOGs. As illustrated in Figure 6A , the five most prevalent KOG classifications are “Function Prediction” (934, 12.41%), “Transduction Machinery” (432, 5.74%), “Post-Translational Modifications, Protein Turnover, Molecular” (396, 4.90%), “Biosynthesis of Secondary Metabolites, Transport and Metabolism” (308, 4.09%), and “Carbohydrate Transport and Metabolism” (260, 3.45%).

In the TPU_vs_TPR comparison, 2333 DEGs were annotated to the KOG database. As illustrated in Figure 6B , the top five KOG classifications are as follows: “Function Prediction” (504, 12.68%), “Transduction Mechanism” (229, 5.76%), “Secondary Metabolite Biosynthesis, Transport and Metabolism” (222, 5.58%), “Post-Translational Modifications, Protein Turnover, Molecular” (176, 4.43%), and “Carbohydrate Transport and Metabolism” (150, 3.77%).

In the TPU_vs_TYU comparison, 2880 DEGs were annotated to KOGs. As illustrated in Figure 6C , the five most prevalent KOG classifications are: “Function Prediction” (639, 12.68%), “Transduction Mechanism” (282, 5.60%), “Post-Translational Modifications, Protein Turnover, Molecular” (230, 4.56%), “Biosynthesis of Secondary Metabolites, Transport and Metabolism” (217, 4.31%), and “Carbohydrate Transport and Metabolism” (199, 3.95%).

This study employed weighted gene co-expression network analysis (WGCNA) to cluster the DEGs in the peel samples of purple and yellow passion fruit at varying maturities. The genes were organized into multiple modules, which were found to be associated with anthocyanin content. As demonstrated in Figure 7A , the clustering of DEGs was executed using anthocyanin content as an external trait. A heat map was generated to illustrate the sample clusters for this trait, revealing that TPR had the highest correlation. The construction of the co-expression network was facilitated by employing an optimal soft threshold, a process that enabled the aggregation of genes into discrete modules. The construction of the gene clustering tree was performed to cluster the distance between these modules, which also illustrates their relationships.

As illustrated in Figure 7B , the upper portion of the figure presents the hierarchical clustering tree of the genes, while the lower portion illustrates the gene modules. The genes that exhibited strong relatedness were observed to cluster together and be assigned to the same module. The heat map illustrates the modules as color blocks on the left and the corresponding ranges as a color bar on the right. In the central section, darker colors indicate stronger correlations, with red signifying positive correlations and the blue denoting negative correlations. The numerical values within each cell denote the degree of correlation and its statistical significance. The results indicated a positive correlation between the turquoise module and anthocyanin content, while a negative correlation was observed between the blue module and the same variable. Subsequent analyses were conducted on the anthocyanin content related to the turquoise module. To identify the genes with a high membership in the turquoise module, the “Gene Significance” (GS) and “Module Membership” (MM) metrics were employed.

As illustrated in Figure 7C , the gene that demonstrates a high degree of correlation with a specific trait exhibits a robust correlation between “Gene Significance” (GS) and “Module Membership” (MM) within its designated module. This finding suggests that these genes play crucial roles within the key modules.

The construction of the gene co-expression network for the two related modules was facilitated by the Cytoscape software. The turquoise module identified 5 hub genes ( Figure 8A ), and the blue module similarly identified 5 hub genes ( Figure 8B ) using the MCC algorithm.

In order to verify the reliability of the transcriptional sequencing data, 10 genes related to anthocyanin content were selected for a real-time fluorescence quantitative PCR verification ( Table 2 ). The β-actin gene was utilized as an internal reference for the test samples by SYBR@Green I chimeric fluorescence method to detect the relative expression levels of each factor at the two maturity levels ( Figure 9 ). The relative expression levels of RT-qPCR were calculated using the 2-ΔΔCT method, while the RNA-seq data were represented by the log2-transformed FPKM values. We conducted a correlation analysis on these two datasets. The results are shown in Figure 10 , with a Pearson correlation coefficient (r = 0.968 > 0.9, p < 0.001), indicating that the two methods have extremely strong consistency in detecting gene expression changes. This high consistency confirms the reliability and accuracy of our transcriptome data.

The results demonstrated the up-regulation of PeCA (4-coumarate-CoA ligase activity), PeMYC2 (transcription factor MYC2), and PeCHI (chalcone isomerase) gene expression and the down-regulation of PeKIA (kinase inhibitor activity) and PeEIA (enzyme inhibitor activity) expressions. Consequently, the expression of these genes was found to be instrumental in the promotion of anthocyanin accumulation and synthesis.