Analysis of Volatile Aroma Components and Regulatory Genes in Different Kinds and Development Stages of Pepper Fruits Based on Non-Targeted Metabolome Combined with Transcriptome

Aroma is a crucial attribute affecting the quality of pepper and its processed products, which has significant commercial value. However, little is known about the composition of volatile aroma compounds (VACs) in pepper fruits and their potential molecular regulatory mechanisms. In this study, HS-SPME-GC-MS combined with transcriptome sequencing is used to analyze the composition and formation mechanism of VACs in different kinds and development stages of pepper fruits. The results showed that 149 VACs, such as esters, alcohols, aldehydes, and terpenoids, were identified from 4 varieties and 3 development stages, and there were significant quantitative differences among different samples. Volatile esters were the most important aroma components in pepper fruits. PCA analysis showed that pepper fruits of different developmental stages had significantly different marker aroma compounds, which may be an important provider of pepper’s characteristic aroma. Transcriptome analysis showed that many differential genes (DEGs) were enriched in the metabolic pathways related to the synthesis of VACs, such as fatty acids, amino acids, MVA, and MEP in pepper fruits. In addition, we identified a large number of differential transcription factors (TFs) that may regulate the synthesis of VACs. Combined analysis of differential aroma metabolites and DEGs identified two co-expression network modules highly correlated with the relative content of VACs in pepper fruit. This study confirmed the basic information on the changes of VACs in the fruits of several Chinese spicy peppers at different stages of development, screened out the characteristic aroma components of different varieties, and revealed the molecular mechanism of aroma formation, providing a valuable reference for the quality breeding of pepper.


Introduction
Peppers are one of the most important vegetables and originate from the tropical and subtropical regions of Central and South America, where more than thirty closely-related species of the genus Capsicum are found. There are six economically cultivated species: Capsicum chinense, C. annuum, C. frutescens, C. baccatum, C. assamicum, and C. pubescens [1], with different fruit shapes, sizes, spiciness, aromas, and tastes. Pepper fruits are rich in capsaicin, minerals, vitamin C, etc., with properties such as nutritional supplement, antioxidant, antibacterial, anti-inflammatory, thrombus dissolution, blood pressure reduction, and cholesterol reduction. As a result, they are widely used as an application in medical, antiseptic, chemical, food, and other fields [2,3].
In this paper, the composition and relative content of VACs in pepper fruits from different varieties and developmental stages were determined by headspace solid-phase microextraction (HS-SPME) combined with gas chromatography-mass spectrometry (GC-MS). Moreover, RNA sequencing (RNA-seq) technology was used to analyze the gene expression of fruits in corresponding varieties and development stages, and the structural genes and co-expression networks related to the synthesis path of aroma volatiles were identified by metabolome-transcriptome association analysis. The results will be helpful in explaining the reasons for the different aromas of pepper fruits from different varieties and developmental stages and analyzing the metabolic pathways and molecular regulation mechanisms of VACs in pepper fruits from the omics level, which will lay a foundation for further improving the aroma of pepper fruits and promoting the improvement of pepper fruit flavor through gene and metabolic engineering methods.

Volatile Aroma Characteristics of Pepper Fruit from Different Samples
The specific sample phenotype information is shown in Figure 1 and Table S1, the total ion chromatograms for the GC-MS detection are shown in Figure S1, and the QC sample quality test is shown in Figure S2. A total of 149 VACs were identified in 36 samples (four varieties × three developmental stages × three biological replicates) and were composed of 53 esters (35.57% of the total), 27 aldehydes (18.12%), 10 alcohols, 20 hydrocarbons (13.42%), 20 terpenoids (13.42%), and 19 other VACs ( Table 1). As the dominant volatiles, esters mainly include Z17 (trans-2-hexenyl isovalerate), Z11 (methyl salicylate), Z7 (namyl isovalerate), Z29 (6-methylhept-4-en-1-yl 3-methyl-butanoate), Z12 (4-methylpentyl 2-methylbutanoate), Z28 (6-Methyl-4-heptenyl 2-methylbutanoate), Z39 (4-methylpentyl 8-methylnonanoate) components, followed by aldehydes such as Q3 (Hexanal), Q4 (e-2hexenal), Q16 (9-undecenal, 2,10-dimethyl-), Q21 (pentadecanal), etc. The third largest number are terpenes and hydrocarbons, the former mainly including T3 (Z-β-Farnesene), T4 (himach-ala-2,4-diene), T1 (α-Cubenene), T20 (α-Ionone), the latter comprising TH15 (hexade-cane), TH4 (2-methyl-hexadecane), TH12 (tetradecane), etc. In this experiment, only three ketones and four phenolic components were detected. We found that the types of volatile compounds detected in different samples were basically the same, while the total relative content of volatiles showed obvious recognizable differences (Figure 2a and Table 1). With the development of fruits, the relative content of total volatiles in GJ and HDL1 gradually increased from 663.67 µg/kg and 811.76 µg/kg to 931.69 µg/kg and 1210.62 µg/kg, respectively. CTJ showed a trend of increasing first and then decreasing, with a maximum value of 761.66 µg/kg, while HDL2 showed a trend of decreasing from 1082.97 µg/kg at the green stage to 796.88 µg/kg at the mature stage. In addition, the concentration of total VACs compounds between different varieties at the  Note: values are means ± standard deviation; "-": indicates not detected; "The total aroma content of volatiles" was calculated as the sum of all aroma metabolites; RT: retention time. We found that the types of volatile compounds detected in different samples were basically the same, while the total relative content of volatiles showed obvious recognizable differences ( Figure 2a and Table 1). With the development of fruits, the relative content of total volatiles in GJ and HDL1 gradually increased from 663.67 µg/kg and 811.76 µg/kg to 931.69 µg/kg and 1210.62 µg/kg, respectively. CTJ showed a trend of increasing first and then decreasing, with a maximum value of 761.66 µg/kg, while HDL2 showed a trend of decreasing from 1082.97 µg/kg at the green stage to 796.88 µg/kg at the mature stage. In addition, the concentration of total VACs compounds between different varieties at the same stage showed noticeable differences, and the concentration of total VACs in HDL2 (1082.97 µg/kg) was significantly higher than that of other varieties at the green stage, while HDL1 was significantly higher than other varieties at the breaking stage (1070.67 µg/kg) and mature stage (1210.62 µg/kg).  16.33% at the green stage to 3.17% at the maturation stage. Besides CTJ, terpenoids had a high content in at least one developmental stage in the other three pepper varieties. The proportion of terpenoids in GJ increased significantly at the breaking and maturation stages (about 20%), which became the second only to esters in pepper fruit. However, there was an opposite phenomenon in HDL2. Terpenoids were mainly accumulated at the green stage (15.12%) and significantly decreased to 2.18% at the maturation stage of HDL2. In CTJ and HDL1, the relative content of terpenoids maintained a stable proportion and was no significant change at three stages, ranging from 4.64 to 6.01% (CTJ) and 11.73 to 14.69% (HDL1), respectively. In addition, the total volatile esters content in HDL1 (712.61 µg/kg) was significantly higher than that in other varieties at maturity, while aldehydes (38.42 µg/kg) and alcohols (47.00 µg/kg) were significantly lower than other varieties (Table S2).  The relative content of VACs (esters, aldehydes, alcohols, terpenoids, etc.) from different samples exhibited a dynamic change in different groups ( Figure 2b and Table S2). In the four varieties, the content ratio of volatile esters in HDL1 and HDL2 during the whole stage reached 50-60%, which was much higher than that of other groups. The relative content of esters in GJ in the green stage was the highest (56.62%), while the proportion of esters in the breaking and maturation stages decreased slightly. Although the proportion of total volatile esters in CTJ was relatively low, it was still above 30.64%. The aldehydes accounted for 16.33% (HDL1), 15.96% (HDL2), 22.84% (GJ), and 24.54% (CTJ) at the green stage, respectively, and decreased gradually with the development of the fruit. The proportion of alcohols in CTJ, HDL2, and GJ increased with fruit development, from 5.41%, 4.43%, and 5.54% to 18.24%, 10.41%, and 9.47%, respectively, while HDL1 showed the opposite accumulation pattern, with 16.33% at the green stage to 3.17% at the maturation stage. Besides CTJ, terpenoids had a high content in at least one developmental stage in the other three pepper varieties. The proportion of terpenoids in GJ increased significantly at the breaking and maturation stages (about 20%), which became the second only to esters in pepper fruit. However, there was an opposite phenomenon in HDL2. Terpenoids were mainly accumulated at the green stage (15.12%) and significantly decreased to 2.18% at the maturation stage of HDL2. In CTJ and HDL1, the relative content of terpenoids maintained a stable proportion and was no significant change at three stages, ranging from 4.64 to 6.01% (CTJ) and 11.73 to 14.69% (HDL1), respectively. In addition, the total volatile esters content in HDL1 (712.61 µg/kg) was significantly higher than that in other varieties at maturity, while aldehydes (38.42 µg/kg) and alcohols (47.00 µg/kg) were significantly lower than other varieties (Table S2).

Analysis of Differential VACs in Pepper Fruits from Different Samples
The difference in VACs from different pepper fruits was analyzed using the relative content of each of the metabolites. Principal component analysis (PCA) showed that the samples were significantly different (Figure 3a). Hierarchical cluster analysis (HCA) indicated that the experimental data were reliable (Figure 3b). In the cluster heat map (Figure 3c), the VACs of different pepper fruit showed significant differences. Most of them were significantly accumulated in the middle and late stages (breaking and maturation stages), clearly separated from the green stage, while metabolites of specific biochemical pathways were without obvious clustering. Moreover, the accumulation patterns of VACs were different between different varieties at the same stage. The differential accumulation of these volatile aroma substances may be the reason for elucidating the aroma formation mechanism of pepper fruit.

Analysis of Differential VACs in Pepper Fruits from Different Samples
The difference in VACs from different pepper fruits was analyzed using the relative content of each of the metabolites. Principal component analysis (PCA) showed that the samples were significantly different (Figure 3a). Hierarchical cluster analysis (HCA) indicated that the experimental data were reliable (Figure 3b). In the cluster heat map (Figure  3c), the VACs of different pepper fruit showed significant differences. Most of them were significantly accumulated in the middle and late stages (breaking and maturation stages), clearly separated from the green stage, while metabolites of specific biochemical pathways were without obvious clustering. Moreover, the accumulation patterns of VACs were different between different varieties at the same stage. The differential accumulation of these volatile aroma substances may be the reason for elucidating the aroma formation mechanism of pepper fruit.  Moreover, we constructed PLS-DA models of VACs to analyze the similarities and differences of VACs in pepper fruits at different development stages, and in different varieties of pepper fruits at the mature stage, based on the fact that fresh peppers are mostly eaten at the mature stage. The quality test results showed that the PLS-DA model was reliable (Figures S3 and S4). Finally, according to the screening conditions VIP > 1, p < 0.05, and |log2FC| ≥ 1, 70 significantly different metabolites (Table S3 and Figure S5) from all the detected VACs were screened. There were 24 different VACs at three development stages fruits of CTJ (CTJ-a vs. CTJ-b and CTJ-a vs. CTJ-c), among which 10 VACs were significantly up-regulated at the breaking and maturation stages compared with the green stage 30 different VACs in GJ (GJ-a vs. GJ-b and GJ-a vs. GJ-c), 12 VACs were significantly up-regulated and 9 VACs were significantly down-regulated at the breaking and maturation stages compared with the green period, Among them, up-regulated VACs included Z8 (4-Methyl-pentyl isobutyrate), Z28, TH4,TH12, T20, and the down-regulated VACs mainly included O4 (2-pentyl-furan), C8 (e-2-hexadecacen-1-ol), Q3, Q4, etc. 32 significantly different VACs were found in the HDL1 comparison group at different developmental stages (HDL1-a vs. HDL1-b and HDL1-a-V HDL1-c), 12 VACs such as Z7, Z29, Z31 (valeric acid, benzyl ester), and TH5 (heptane) were significantly upregulated compared to the green stage in the middle and late stages, while Z11, Z12, Z15 (butyric acid, 2-methyl, hexyl ester), C3, C8, and the other 9 VACs were significantly downregulated. In addition, 30 VACs were identified as significant differences in HDL2 (HDL2-a vs. HDL2-b and HDL2-a vs. HDL2-c), including 6 up-regulated and 12 down-regulated VACs in the middle and late stages. Among them, up-regulated VACs included Z39, Z22 (4-methyl-amyl tiglate), C9, TH19 (z-3-heptadecene), etc., and down-regulated VACs mainly included Z35 (Isopentyl 8-methylnon-6-enoate), Z36 (pentadecanoic acid, 3-methylbutyl ester), Z23 (5-methyl-hexyl 3-methylbutanoate), Z11, Z12, etc. Forty VACs showed significant differences in mature pepper fruits of different varieties (GJ-c vs. CTJ-c, HDL1-c vs. CTJ-c, and HDL2-c vs. CTJ-c), among which volatile esters, aldehydes, alcohols, terpenoids, and hydrocarbons were the main VACs, such as Z12, Z13, Q25 (z-9,17-octadecadienal), C8, T3, Q16, etc. The VACs in pepper fruits of different varieties showed different accumulation patterns compared with CTJ mature fruits, HDL1 and GJ had more up-regulated metabolites, while HDL2 had more down-regulated metabolites. Notably, 31 of all 53 esters (58.49% of the total) showed significant differences in different comparison groups, while 22 esters were branched-chain esters, indicating that esters (especially branched-chain esters) may play a more important role in the change of pepper fruit aroma.

Screening of Marker Aroma Compounds in Pepper Fruits from Different Samples
Interestingly, most of VACs with high content at the green stage decreased gradually with fruit ripening, significantly different from the differential metabolites with high abundance in chili fruit at the breaking and maturation stages. Therefore, to screen out the marker aroma compounds in pepper fruits at different stages, we performed PCA analysis on the relative expression of 70 selected differential metabolites. In the PCA loading plot ( Figure S6), the variable metabolites in different quadrants were the most important VACs describing the specific stage of pepper in this study. The results showed that Q3 and Z11 were the marker volatiles of CTJ at the green stage, Z13 and TH13 (tetradecane, 2-methyl-) were the typical aroma components of CTJ at the breaking stage, and F1 and Z12 were the marker aroma compounds at the mature stage. In GJ, Z11, Z12, and Q3 were significantly expressed at the green stage, and Z17, T4, and T3 were the characteristic volatiles at the breaking stage, and Z29, Z7, Z22, and C2 were the marker aroma compounds of mature fruit. In HDL1 pepper fruit, Z11 and Z13 were the marker aroma compounds at the green stage, Z29 and Z38 (4-methyl-pentyl 8-methylnon-6-enoate) were the marker volatiles at the breaking stage, and Z17 and T3 were significantly expressed at the mature stage. At the green stage of HDL2, the main marker volatiles were T4 (Himachal-2,4-diene) and Z11, while Z17 was significantly expressed at the breaking stage. Z32 (2-methyl butyl 8-methyl-6-enoate) is a marker aroma compound in mature pepper fruit.

Transcriptome Analysis of Pepper Fruits from Different Samples
To further explore the expression changes of genes related to VACs in pepper fruit, we sequenced the transcriptome of 36 samples. Sequencing quality and data statistics are shown in Table S4, and pepper reference genome mapping data statistics are shown in Table S5. We then analyzed the gene expression patterns in different samples of pepper fruits based on the FPKM values of genes ( Figure S7). Cluster analysis (Figure 4a) showed that the biological repetition of pepper fruit samples from different samples had good classification. According to the PCA analysis results (PC1: 60.86% and PC2: 25.4%), the samples differed significantly in different varieties and stages (Figure 4b). Then, we used p < 0.05 and |log2FC| ≥ 1 as the screening conditions to screen the DEGs between different comparison groups. As shown in (Figure 4c,d), the number of DEGs in the four pepper varieties showed an increasing trend first and then decreased (Figure 4c), indicating that the gene expression change is the most active at the breaking stage, which may be the critical period for pepper to synthesize and metabolize VACs. The comparison of different varieties (Figure 4d) showed that there were more DEGs at the breaking stage and at the mature stage, which might be the critical period for the formation of differential VACs in different varieties. Furthermore, DEGs in different comparison groups showed that the number of down-regulated genes was much higher than that of up-regulated genes. the mature stage. At the green stage of HDL2, the main marker volatiles were T4 (Himachal-2,4-diene) and Z11, while Z17 was significantly expressed at the breaking stage. Z32 (2-methyl butyl 8-methyl-6-enoate) is a marker aroma compound in mature pepper fruit.

Transcriptome Analysis of Pepper Fruits from Different Samples
To further explore the expression changes of genes related to VACs in pepper fruit, we sequenced the transcriptome of 36 samples. Sequencing quality and data statistics are shown in Table S4, and pepper reference genome mapping data statistics are shown in Table S5. We then analyzed the gene expression patterns in different samples of pepper fruits based on the FPKM values of genes ( Figure S7). Cluster analysis (Figure 4a) showed that the biological repetition of pepper fruit samples from different samples had good classification. According to the PCA analysis results (PC1: 60.86% and PC2: 25.4%), the samples differed significantly in different varieties and stages (Figure 4b). Then, we used p < 0.05 and |log2FC| ≥ 1 as the screening conditions to screen the DEGs between different comparison groups. As shown in Figure 4c,d), the number of DEGs in the four pepper varieties showed an increasing trend first and then decreased (Figure 4c), indicating that the gene expression change is the most active at the breaking stage, which may be the critical period for pepper to synthesize and metabolize VACs. The comparison of different varieties (Figure 4d) showed that there were more DEGs at the breaking stage and at the mature stage, which might be the critical period for the formation of differential VACs in different varieties. Furthermore, DEGs in different comparison groups showed that the number of down-regulated genes was much higher than that of up-regulated genes. To further study the change of genes at the transcriptional level, we performed Venn analysis on the DEGs (Figures 5a-c and S8A-C). There were 5282, 4790, and 371 common DEGs (Figure 5a-c) at three stage comparison groups of the four varieties (a vs. b, a vs. c, and b vs. c), respectively. While 3273, 4807, and 4548 common DEGs ( Figure S8A-C) were in four variety comparison groups at three stages (GJ vs. CTJ, HDL1 vs. CTJ, and HDL2 vs. CTJ), respectively. Furthermore, the results of GO enrichment analysis showed that at level 2, a total of 45 functional subclasses were enriched between the comparison groups at stages (Figure 5d-f), including 22 biological processes, 13 cellular components, and 10 molecular functions. A total of 46 functional subclasses, including 21 biological processes, 13 cellular components, and 12 molecular functions, were enriched among different varieties ( Figure S8D-F). KEGG enrichment analysis showed that 117, 114, and 57 metabolic pathways were enriched in three comparison groups at three stages. The pathways related to the synthesis and metabolism of VACs mainly included carotenoid biosynthesis (cann00906), fatty acid elongation (cann00062), alpha-Linolenic acid metabolism (cann00592), phenylpropanoid biosynthesis (cann00940), and valine, leucine and isoleucine degradation (cann00280), fatty acid degradation (cann00071), and fatty acid biosynthesis (cann00061) (Figure 5g-i). These metabolic pathways were significantly enriched (p < 0.05) in the comparison group of green and maturity stage, green stage and breaking stage, but were generally not significantly enriched in the comparison group of the maturity and breaking stages, indicating that the gene expression related to the synthesis of VACs mainly occurred in the period before the breaking stage, among which the breaking stage is the critical turning point. Common DEGs among different varieties were enriched to 34 metabolic pathways at the green stage, but most pathways were not significant ( Figure S8G). At the breaking stage, 116 metabolic pathways ( Figure S8H) were enriched, among which the pathways related to the synthesis of VACs were significantly enriched, including carotenoid biosynthesis, phenylalanine metabolism (cann00360), and fatty acid degradation. The maturation stage of different pepper variety comparison groups enriched to 115 metabolic pathways ( Figure S8I). Phenylpropanoid biosynthesis is the most significant enriched metabolic pathway, followed by ubiquinone and other terpenoid-quinone biosynthesis (cann00130) and phenylalanine metabolism, and the volatile ester synthesis of linoleic acid metabolism (cann00591) pathway is significantly enriched.

The Expression Changes of Key Genes in the Biosynthesis Pathway of VACs
Volatile esters (especially branched-chain esters) and aldehydes, which are the most abundant and changing substances in pepper fruits, are mainly produced by the fatty acid pathway (FAP) and branched-chain amino acid pathway (BCAAP). According to transcriptome annotation and homologous sequence alignment, 90 genes involved in FAP and 35 genes involved in BCAAP were identified with differential expression ( Figure 6 and Table S6). In the FAP (Figure 6 (left)), we identified 13 genes encoding LOX, 1 HPL, 11 ADH, and 3 AAT, among which most of the DEGs were significantly downregulated in pepper fruit development (p < 0.05), which is consistent with the change of VACs of straight-chain fatty acids such as hexanal and (e)-2-hexenal, indicating that they may play an important role in this process. Similarly, the genes encoding ACC, MAT, FAD, AOS, and ACX showed similar expression patterns. In addition, we found the LOC107850471 gene encoding SAD, the LOC107839570 gene encoding 9S-LOX, and the LOC107867711 gene encoding ADH had significantly high expression at the breaking and maturation stages, indicating that these genes may be associated with the significant accumulation of VACs of straight-chain fatty acids such as hexyl hexanoate in middle and late pepper fruits. Compared with different varieties, most of the DEGs showed lower expression levels compared with CTJ in GJ, HDL1, and HDL2 pepper fruits of the same period.

The Expression Changes of Key Genes in the Biosynthesis Pathway of VACs
Volatile esters (especially branched-chain esters) and aldehydes, which are the most abundant and changing substances in pepper fruits, are mainly produced by the fatty acid pathway (FAP) and branched-chain amino acid pathway (BCAAP). According to transcriptome annotation and homologous sequence alignment, 90 genes involved in FAP and 35 genes involved in BCAAP were identified with differential expression (Figure 6 and Table S6). In the FAP (Figure 6 (left)), we identified 13 genes encoding LOX, 1 HPL, 11 ADH, and 3 AAT, among which most of the DEGs were significantly downregulated in pepper fruit development (p < 0.05), which is consistent with the change of VACs of straight-chain fatty acids such as hexanal and (e)-2-hexenal, indicating that they may play an important role in this process. Similarly, the genes encoding ACC, MAT, FAD, AOS, Branched-chain aminotransferase (BCATs) perform the last step of synthesis of Lleucine, L-isoleucine, and L-valine and the first step of catabolism to produce volatile branched-chain esters in the BCAAP (Figure 6 (right) and Table S6). Three BCAT genes showed significant differences in different samples, among which LOC107851037 was significantly downregulated during pepper fruit development, while LOC107867296, showed higher expression at the breaking and maturation stages, indicating that the gene of this family may have some control in the generation of volatile esters in BCAAP comprehensively and coordinately. Interestingly, the expression level of LOC107867296 was significantly higher in HDL2, which may be related to the significant accumulation of 4-methylpentyl 8-methylnonanoate (Z39) in pepper fruit at the breaking stage of HDL2, but was almost undetectable in the other three pepper varieties. As the key enzyme that catalyzes the straight-chain aldehydes to produce acids in BCAAP, seven genes encoding Aldehyde dehydrogenase (ALDH) were identified, and most of these genes' expression are significantly up-regulated in the fruit development, among which LOC107862407 has more significant expression, indicating that it may have a more critical role in the synthesis of branched ester. In our study, 10 Carboxylesterase (CXE) genes were identified, among which LOC107872419 and LOC107859022 gene expression were significantly reduced in the middle and late stages of pepper fruit development, which may be related to the significant increase of most branched-chain esters in pepper fruit at breaking and maturation stages. In addition, we found that the LOC107867571 and LOC107859022 genes encoding CXE were more highly expressed in HDL2 than the other three varieties, which may associate with 4-methylpentyl 8-methylnonanoate (Z39) significant accumulation at the mature stage.  A total of 20 volatile terpene compounds were detected in this experiment, mainly consisting of monoterpenes and sesquiterpenes. As shown in Figure 7 and Table S6, our transcriptome data identified 49 differentially expressed genes involved in MVA and MEP pathways. 1-Deoxy-D-xylulose 5-phosphate synthase (DXS), 3-hydroxy-3-methyl-glutaryl-CoA synthase (HMGS), and 3-hydroxy-3-methylglutaryl-CoA reductase (HMGCR) are considered to be the critical rate-limiting enzymes in the MVA and MEP pathways. In Methyl salicylate (Z11, MeSA) is widely found in various fruits and vegetables, often exhibiting green and minty aromas. In this experiment, MeSA mainly accumulated significantly and was the primary characteristic VACs of pepper fruits at the green stage. We detected 12 DEGs involved in MeSA synthesis, including 3 genes encoding DAHP synthase, 3 genes encoding AroE, 1 gene encoding MtSK, 1 gene encoding ICS and 4 genes encoding SAMT, most of these genes had higher expression levels ( Figure S9 and Table S6) at green stages. Moreover, comparing CTJ with the other three pepper varieties, we found that one gene encoding DAHP synthase (LOC107867463) and two genes encoding SAMT (LOC107864268 and LOC107849647), were significantly more expressed at breaking and mature stages, which was consistent with the accumulation pattern of MeSA, indicating that these genes have a significant regulatory role in the MeSA biosynthesis in pepper fruit.
A total of 20 volatile terpene compounds were detected in this experiment, mainly consisting of monoterpenes and sesquiterpenes. As shown in Figure 7 and Table S6, our transcriptome data identified 49 differentially expressed genes involved in MVA and MEP pathways. 1-Deoxy-D-xylulose 5-phosphate synthase (DXS), 3-hydroxy-3-methyl-glutaryl-CoA synthase (HMGS), and 3-hydroxy-3-methylglutaryl-CoA reductase (HMGCR) are considered to be the critical rate-limiting enzymes in the MVA and MEP pathways. In this study, we identified three differentially expressed DXS genes, three HMGs, and five HMGCRs. Although the number of these gene families was small, their expression levels have changed significantly in pepper fruits of different samples, which may play an important role in accumulating terpene compounds in pepper fruits. As the last key enzyme for the synthesis of terpenes, 16 TPSs were identified, most of which were significantly expressed at the green stage, and the expression of LOC107868331, LOC107846955, LOC107877838, and LOC107840471 was more prominent, indicating that they play a more critical role in the synthesis of terpene aromatic compounds. The expression level of LOC107878220 (TPS) was significantly higher in the fruits of GJ at the breaking and maturation stages than in the other three varieties, which was consistent with the accumulation pattern of T4 (himachali-2,4-diene) in the four pepper varieties. In this experiment, five genes encoding carotenoid cleavage dioxygenase (CCD) were identified, which may be involved in regulating the formation of α-ionone in pepper fruit.

Gene Expression Verified by Real-Time Quantitative PCR (qRT-PCR)
Fifteen key DEGs were selected for qRT-PCR verification analysis, among which SAD, ADH, SAD, HPL, AOS, AOC, FAD, LOX, ACX, BCAT, PD, and FAD in FAP and BCAAP, and TPS, and DXS genes were related to volatile terpenoid synthesis. As shown in Figure S10, the expression profiles of 15 genes were consistent with the expression results of RNA-Seq, indicating that the transcriptome data in this experiment were accurate and reliable.

Correlation Analysis of TFs and Key DEGs in the Synthesis Pathway of Aroma Components
TFs play an important role in the synthesis of plant VACs. A total of 2755 genes were identified as TFs, of which 1965 genes showed significant differences between at least two samples and were further divided into 58 families. Figure 8a shows the cluster heatmap of FPKM values of all DEGs in the same TF family in each group. The results showed that these TF families had significant differential expression in different pepper fruit samples. Furthermore, we analyzed the correlation between the key DEGs in the FAP, BCAAP, and terpenoid biosynthesis pathways and the differential TFs. As shown in Figure 8b, a large number of key DEGs involved in the synthesis of aroma components have a significant correlation with TFs (|r| > 0.9), p-value < 0.05), indicating that these transcription factors may be involved in the volatile biosynthesis by regulating the expression of key genes in the VACs synthesis pathway. LOC107877838, and LOC107840471 was more prominent, indicating that they play a more critical role in the synthesis of terpene aromatic compounds. The expression level of LOC107878220 (TPS) was significantly higher in the fruits of GJ at the breaking and maturation stages than in the other three varieties, which was consistent with the accumulation pattern of T4 (himachali-2,4-diene) in the four pepper varieties. In this experiment, five genes encoding carotenoid cleavage dioxygenase (CCD) were identified, which may be involved in regulating the formation of α-ionone in pepper fruit.

Co-Expression Network Module Analysis Identifies DEGs Modules Involved in the Biosynthesis of VACs
To gain insight into the molecular mechanisms regulating the formation of VACs in pepper fruits, we constructed gene networks and co-expression analysis using nonredundant DEGs. The results showed that these DEGs were clustered into ten major modules, and genes in the same module showed higher co-expression and correlation relationships. The correlation analysis between the gene module and trait module with 70 differential VACs as phenotypic traits showed that ( Figure 9a) the dark-olive-green module was positively correlated with the accumulation of most aroma components, while the light-yellow module was negatively correlated with the most VACs. The genes in these two modules were used for further cluster analysis, and the results showed that the genes in the dark-olive-green module (Figure 9b) had lower expression levels in CTJ than in other samples, while the genes in the light-yellow module showed the opposite expression patterns, higher expression of genes was mainly in CTJ (Figure 9c). In order to screen for hub genes associated with the synthesis of VACs, we investigated key genes in the volatile biosynthesis pathway in the dark-olive-green and light-yellow modules and TFs with significant correlation with these essential genes (r > 0.8 or <−0.8, p-value < 0.05). A total of 8 TF genes and 3 VACs synthesis-related genes were identified in the dark-olive-green module, and 22 TF genes and 4 VACs synthesis-related genes were identified in the light-yellow module (Table S7). We then selected these TFs and aroma synthesis-related genes as essential hub genes to construct a co-expression network of the dark-olive-green and light-yellow modules. The results showed that (Figure 9d,e) most of these essential hub genes are closely related to other genes in the plate and are closely linked. It can be speculated that they play an essential role in regulating the synthesis of aroma substances. Meanwhile, many TFs were identified as crucial hub genes, indicating that they play an essential role in forming VACs in pepper fruits. Notably, the genes directly involved in the aroma volatile synthesis pathway did not largely cluster in the dark-olive-green and light-yellow modules, possibly due to the low expression level of most genes involved in the volatile synthesis and the diversity of volatile aroma products. two modules were used for further cluster analysis, and the results showed that the genes in the dark-olive-green module (Figure 9b) had lower expression levels in CTJ than in other samples, while the genes in the light-yellow module showed the opposite expression patterns, higher expression of genes was mainly in CTJ (Figure 9c). In order to screen for hub genes associated with the synthesis of VACs, we investigated key genes in the volatile biosynthesis pathway in the dark-olive-green and lightyellow modules and TFs with significant correlation with these essential genes (r > 0.8 or

Aroma Perception from Dynamic Changes of VACs in Pepper Fruits
The volatile components and relative content are important reasons for people's unique aroma perception of different types of fruits and vegetables [29]. In this study, we detected more than 140 VACs with different contents in a variety of pepper fruits, indicating that the aroma differences were mainly affected by the relative content of VACs. We found that HDL1 (Hainan Huangdenglong pepper) and HDL2 (hybrid Huangdenglong pepper) are both processed varieties of chili sauce, and their genetic relationship is similar. However, HDL1 mature fruit has a higher content of flavor components, and the accumulation pattern of HDL2 flavor components is the opposite, which may be the reason why its chili sauce is more popular with consumers. Studies have shown that straight and branched-chain esters can give Capsicum chinense pepper fruits a solid fruity/exotic aroma [23,30]. In our study, esters were the most critical aroma components in four pepper varieties, similar to the results of previous studies on Cachucha peppers [23], Habanero pepper [31], and Tabasco pepper [32] and accumulated in the late stage of most fruit development in HDL1 and GJ (ghost pepper), which is consistent with the accumulation pattern of fruity esters [33,34], further explaining the cause of the processing aroma of HDL1 mature fruit. However, the accumulation of esters gradually decreased during HDL2 and increased first and then decreased during CTJ (Chaotian pepper) fruit development, indicating the accumulation pattern difference of esters in different cultivars and developmental stages. The content of esters in HDL 1 mature fruit was significantly higher than other varieties at the same stage, which may be an important reason for for its rich aroma. Aldehydes were the main aroma components at the green stage of pepper fruits and showed a decreasing expression trend with fruit development, which may be related to the gradual transformation of aldehydes into esters with fruit ripening. Terpenes, showing floral and fruity aroma in tea, were detected in pepper fruit in this experiment, and the accumulation pattern was affected by varieties and development stages as well as esters. Although their concentration is usually low, these compounds usually have a low aroma threshold [35], which has an important impact on the aroma quality of pepper.

Analysis of Marker Aroma Compounds in Pepper Fruits from Different Samples
Generally, the special flavor of vegetables, fruits, and tea comes from the marker aroma compounds. (E)-2-hexenal and (Z)-3-hexenal are the most important aromatic compounds with strong green grass and green leaf odor in cherry tomato [36], while linalool oxide and linalool are the main components of citrus flower fragrance [37], and hexyl 2-methylbutyrate is the main reason for the "honey crisp" apple flavor [34]. The characteristic components (E, E)-2,4-hexagon make Jingshan green tea have a floral and fatty flavor [38]. In Brazil's tabasco pepper, α-Pinene, hexyl 2-methylbutyrate, hexyl 3methylbutanoate, and hexyl butanoate are the main contributors to the typical sweet, herbal, and pepper-like flavors [32]. Our research confirmed methyl salicylate (Z11) was the main marker aroma compound at the green stage of four pepper varieties and described as green and sweet in the Brazilian tabasco pepper [32], hexanal (Q3) described in cherry tomatoes as having the smell of freshly cut grass [36], are another characteristic VACs of CTJ and GJ at green stage, however 4-methylpentyl-2-methylbutyrate (Z12) not only accumulated significantly at the green phase of GJ, and was the characteristic volatile of CTJ mature fruit with herbal flavor [32]. Moreover, as one of the main marker volatiles of CTJ and HDL1, 4-methylpentyl 3-methylbutanoate (Z13) may show the same odor as Z12. The study found that the characteristic component trans-2-hexenylisovalerate (Z17) with green and fruit-like odor accumulated significantly at the green stage of HDL2 and breaking stage of GJ, while the hexyl hexanoate (Z27) with an apple fruit flavor is the main marker volatile component in the mature fruit of HDL1, and is usually described as fresh cut grass and vegetable flavor [34]. The above showed that the marker aroma components are different in pepper fruits of different varieties and development stages, which can be used as important indicators to study the flavor quality of pepper. However, how they affect the specific aroma quality of these pepper fruits still needs to be determined. Further sensory evaluation and gas-chromatography-olfactometry (GC-O) methods are needed to study the specific effects of these volatile aroma substances on the aroma of pepper fruit [33,39].

Many DEGs and TFs Are Involved in the Synthesis of VACs in Pepper Fruits
The expression of key genes in amino acids, fatty acids, MVA, and MEP metabolism are an essential factor affecting the formation of aldehydes, alcohols, and esters producing VACs [9,14]. The Lipoxygenase pathway is the key pathway for the biosynthesis of straightchain aldehydes, alcohols, and esters, including four-step enzymatic reactions, in which LOX, as the first step in the synthesis of straight-chain volatile fatty acid derivatives, has been proven to be able to catalyze the oxidation of polyunsaturated fatty acids in apples to produce volatile C6 aldehydes [40]. HPL catalyzes the conversion of LOX reaction products into aldehydes, and its silencing significantly reduces the production of C6 volatiles in olive [41], such as hexanol, hexanal, and (E)-2-hexanol. ADH catalyzes the reduction of straight-chain aldehydes in fruits to corresponding straight-chain alcohols, which has been proven to be closely related to the formation of watermelon C9 alcohols [42]. AAT catalyzes linear alcohols and acyl carrier proteins (ACPs) to form volatile straight-chain esters in many horticultural plants, such as flowers, fruits, and vegetables [6]. Our study identified 13 LOX, 1 HPL, 11 ADH, and 3 AAT genes with different expression patterns in pepper fruits of different varieties and developmental stages with similar trends to some volatile compounds, such as hexanal, (e)-2-hexenal, hexanoic acid-hexyl ester. This suggests that these DEGs may regulate the production of straight-chain volatiles in pepper fruit.
Branched-chain esters are the most critical VACs in pepper fruits in this study, mainly produced by BCAAP. As the last enzyme in the BCAAP, the catalytic product branchedchain amino acid of BCAT provides a precursor for the synthesis of volatile branched-chain compounds [15]. Studies have shown that BCAT is related to the formation of branchedchain esters in fruits such as bananas [43], melons [16], and apples [29]. Another study on peppers showed that the accumulation of most branched esters was significantly positively correlated with the CcBcat4 gene [33]. Our study confirmed three differentially expressed BCAT, of which the expression pattern of LOC107867296 was consistent with the changing trend of 4-methylpentyl-8-methyl-nonanoate (Z39) content, indicating that it may play a key regulatory role in Z39. CXE is a key enzyme for the synthesis of branched-chain esters, which has been shown to hydrolyze flavor esters in peach [17], strawberry [44], and tomato [45]. In this study, 10 genes encoding CXE were identified, of which LOC107867571 and LOC107859022 encoding CXE gene decreased significantly in the breaking and maturation stages of pepper fruit, which may be related to the accumulation of volatile branchedchain esters in pepper fruit in the middle and late growth. Methyl salicylate (Z11) produced by SAMT-catalyzed SA is the most important aromatic compound in pepper fruit [46]. We identified four SAMT-coding genes; however, their expression pattern is not completely consistent with the accumulation pattern of methyl salicylate (MeSA), which indicates that the synthesis of MeSA in pepper fruit may be regulated by other factors that need further study. As the most important enzyme-producing terpene compounds and the gene family with diverse functions in plants [14], TPS directly participated in the synthesis of monoterpene, sesquiterpene, and carotenoid derivatives, the main VACs of HDL1 and GJ fruits. In our experiment, 16 genes encoding TPS were identified as having different expression patterns, which may be related to their functions and the diversity of substrates produced. In addition, we found many other key DEGs related to the synthesis pathway of volatile aromatic compounds, which may play an important role in regulating the synthesis of volatile aromatic compounds in pepper fruits.
TFs play an essential role in the metabolism of aroma volatiles [47]. It has been reported that MYB regulates the production of volatile compounds in strawberries [48], PAP1 can enhance the synthesis of phenylpropanoids and terpenoids in roses [49], the synthesis of volatile esters and monoterpenoids in kiwifruit is closely related to NAC transcription factors [50], AP2/ERF transcription factors have been proven related to the production of volatile terpenes in sweet orange fruits [20], and in tomatoes, RIN transcription factors can regulate the production of volatile aldehydes and alcohols in fruits by regulating the expression of key genes in the LOX pathway [51]. In our experiment, we found that many TFs (C2H2, C3H, and MYB, etc.) were significantly correlated with the critical DEGs and differential metabolites involved in the synthesis of aroma components, indicating that these TFs may be involved in the synthesis of aroma components by regulating the expression of critical genes in the synthesis pathway of VACs in pepper fruits. Therefore, future studies will focus on the specific functions of these TFs in regulating the synthesis of volatile aroma compounds in pepper.

Plant Materials
The four pepper varieties used in this experiment were preserved by the Pepper Research Group of Hainan University respectively, including "HDL1" (Capsicum chinense L., local varieties unique to Hainan Province, China), "HDL2" (Capsicum chinense L., an exotic hybrid whose fruit shape is different from HDL1 is being popularized in Hainan Province), "GJ" (Capsicum chinense L., a highly hot variety originating from India), "CTJ" (Capsicum annuum L., an annual dried pepper variety), specific sample phenotype information is shown in Table S1 and Figure 1. The four pepper varieties are representative in terms of their close affinity, economic species, and varieties, which can reflect the cultivation of pepper in China. All peppers were planted in plant growth rooms in October 2021 (temperature: 23 ± 0.5 • C; light duration: 16 h; light intensity: 6500lx), with ten plants for each cultivar, and cultivation and management measures were consistent. From December 2021 to February 2022, fruits of four pepper varieties were collected at the green, breaking, and maturation stages based on fruit size and color. Thirty fruits were randomly selected from each stage and divided into three biological replicates, totaling 36 fruit samples. Each fruit sample was homogenized by a mixing method, quickly frozen in liquid nitrogen, and stored in an ultra-low temperature refrigerator at −80 • C for subsequent transcriptome and metabolomics analysis.

Sample Preparation and Extraction
Shanghai OE Biotech Co., Ltd. completed untargeted metabolome sequencing. HS-SPME extracted VACs from pepper fruit samples: Fruit samples are frozen fresh fruit and were immediately powdered with a grinder (MM400, Retsch, Dusseldorf, GER), then defrosted at room temperature prior to analysis. Accurately weigh the 1 g of sample powder and add to an ethanol solution containing 20 µg/mL n-alkanes (C7-C40) as an internal standard (configuration method: 980 µL ethanol solution (99.7%) + 20 µL n-alkanes mother liquor (1 mg/mL)) and then transfer to a sealed headspace injection bottle (20 mL, Agilent, Santa Clara, CA, USA) for the release of volatiles. The headspace bottle was equilibrated at a temperature of 60 • C and a shaking speed of 450 rpm for 10 min. Then the extraction head (50/30 µm DVB/CAR/PDMS, Sigma, Shanghai, China) was inserted into the headspace part of the sample and extracted for 60 min [52,53].

GC-MS Analysis of the VACs
Agilent 7890B chromatography and a 5977B mass spectrometer equipped with a DB-5MS capillary column (30 m × 0.25 mm × 0.25 µm, Agilent J&W Scientific, Folsom, CA, USA) carried out the qualitative and quantitative analysis of volatile aromatic compounds. The extracted sample was directly injected into the injection port of the gas chromatographmass spectrometer, desorbed at 250 • C for 5 min, and the injection port temperature was set at 230 • C. The high-purity helium (Purity ≥ 99.999%) was used as the carrier gas. The initial oven temperature was set at 40 • C for 1 min, increased to 230 • C at 4 • C/min and held for 1.5 min, then increased to 250 • C at a rate of 10 • C/min and kept for 2 min. Mass spectrometry was recorded in the ionization mode of electron impact ion source (EI) energy of 70 eV. The ion source temperature was set to 230 • C, and the quadrupole temperature was 150 • C. Mass spectrometry data were extracted by full scan mode (SCAN) with a mass scanning range of 40-500 m/z. In the process of mass spectrometry analysis, all samples were mixed as quality control (QC) samples to test the stability of the system mass spectrometry platform during the whole experiment. MS-DIAL software analyzed the GC-MS raw data for peak detection, peak recognition, MS2 Dec deconvolution, and peak alignment. They then compared with MS in the NIST database to determine the species of VACs. VACs with MS matching scores greater than 80% were retained, and the final substance type was determined after manual identification and comparison. The relative content of volatiles was calculated by the standard internal method. The formula is as follows: f[relative content of volatiles (µg/kg) ] = peak area of target × 0.1 µg peak area of internal standard × 1000 g

Screening of Differential VACs and Marker Aroma Compounds
Multivariate statistical analysis was used to screen the differential VACs in samples: The relative content of VACs was used as the analysis data, and the PLS-DA model of VACs in different pepper fruit samples was constructed by (PLS-DA). The VACs in the model were screened by variable weight value (VIP) > 1, and a T-test was performed to verify its significance.
Reference: Liu et al.'s method for screening marker aroma compounds [24].

RNA Extraction and Library Preparation
Total RNA was extracted from 36 samples (4 pepper varieties × 3 developmental stages × 3 biological replicates) using FastPure ® Plant Total RNA Isolation Kit (Vazyme, China). The purity and quantitative analysis of extracted RNA were evaluated using the NanoDrop 2000 ultraviolet spectrophotometer (Science Corporation, USA). The integrity of RNA was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The qualified total RNA was used to construct the cDNA library.

RNA Sequencing Analysis
The transcriptome sequencing and DEGs expression analysis were performed by Shanghai OE Biotech Co., Ltd. (Shanghai, China). Pre-process the raw data to remove low-quality reads containing Ploy-N. The obtained clean reads were mapped to the pepper genome Pepper Zunla_1_Ref_v1.0 (https://www.ncbi.nlm.nih.gov/data-hub/genome/ GCF_000710875.1/ (accessed on 30 March 2022).) using the default parameter HISAT2 software. Then the transcript assembly was performed, and the transcript abundance was quantified. The calculation formula of transcript abundance FPKM value is as follows: f(FPKM) = cDNA Fragments Mapped Fragments (Millions) × Transcript Length (kb)

DEGs Analysis
DEGs expression analysis was performed using DESeq (2012) R software and identified as DEGs under the conditions of p-value < 0.05 and |log2FoldChange| > 1. GO (Gene Ontology) enrichment analysis was performed using GOseq R package software, and KEGG (Kyoto Encyclopedia of Genes and Genomes) signaling pathway enrichment analysis was performed using KOBAS 2.0 software.

qRT-PCR Validation of Transcriptome Data
Total RNA was extracted from 36 pepper samples using FastPure ® Plant Total RNA Isolation Kit (Vazyme, China), and reverse transcription was performed using HiScript ® III All-in-one RT SuperMix (Vazyme, China). 15 DEGs related to the synthesis of VACs were selected for qPCR, and gene-specific primers were designed using Primer3 software (Table S8). The qRT-PCR reaction was performed using Quant Studio 3 real-time fluorescence quantitative PCR instrument (ABI, Shanghai, China) and ChamQ Universal SYBR qPCR Master Mix. Quantitative Expression Data of Genes Calculated by 2 −∆∆Ct Method with Actin Reference Gene.

Construction of Co-Expression Analysis by WGCNA
The co-expression network was constructed with 6204 DEGs (p-value < 0.05 and |log2FoldChange| > 1 with FPKM ≥ 1 and standard deviation ≤ 0.8), and the weighted gene co-expression network analysis was performed using the R software package. In order to obtain genes related to the synthesis of aroma volatiles, the relative content data of 70 differential metabolites were used as phenotypic traits to correlate with modules, and the significance of DEGs analyzed the correlation between trait data and gene expression data. We used Cytoscape3.9.1 software to visualize the co-expression module.

Statistical Analysis
Trials were conducted entirely randomized, with three replicates per trial and data expressed as mean ± standard deviation (SD). SPSS v19.0 software was used to perform statistical analysis and one-way analysis of variance (ANOVA) of the data, and the Duncan test was used to evaluate the significant difference between samples (p < 0.05). Graphics were mainly drawn using Tbtools, Origin 2021, Excel 2021, Cytoscape 3.9.1, prism8, and other drawing software.

Conclusions
In summary, we found significant quantitative differences in VACs from several Chinese spicy peppers and one India pepper, and volatile esters were the most crucial aroma components in different variety and developmental stage pepper fruits. The higher accumulation of total VACs and volatile esters in Hainan Huangdenglong pepper during fruit development may be the reason for better processing aroma quality. We found that the marker aroma compounds in different pepper fruit samples were different, which may have a crucial effect on the aroma of pepper fruit. Meanwhile, many DEGs and TFs related to volatile compound biosynthesis were identified in pepper fruit. Moreover, the breaking stage may be critical for forming aroma compounds. At last, the comprehensive analysis of transcriptome and metabolomics identified potential candidate genes and transcription factors that may be involved in volatile compound biosynthesis in pepper, which can be essential candidate genes for subsequent functional verification studies. This paper is the first attempt to systematically study the metabolic pathways and molecular regulation mechanisms of VACs in pepper, which provides a theoretical basis for the molecular breeding of high-quality pepper.

Data Availability Statement:
The datasets presented in this study can be found in online repositories of the National Center for Biotechnology Information (NCBI). The accession number is PRJNA932246.