To investigate the genetic relationships and phenotypic diversity among these citrus varieties, phylogenetic analysis and morphological observations were performed. Phylogenetic analysis was conducted using the TPM1uf+I+G4 model parameters based on 222,094 genome-wide single-nucleotide polymorphism (SNP) sites. Yuzu (Citrus junos) originated in the upper reaches of the Yangtze River in China and is considered a hybrid between Ichang papeda and mangshanyeju mandarin. It was later introduced to Japan and Korea during the Tang Dynasty (Rahman et al., 2001; Wu et al., 2021). In this study, ‘Ziyangxiangcheng’ and ‘Xiecheng’ are indigenous Chinese cultivars; ‘Kitou’, ‘Tada Nishiki’, and ‘Zairai’ are native Japanese cultivars; ‘Line131’ is a regional Korean cultivar. Phylogenetic results indicate that yuzu cultivars exhibit closer genetic relationships among themselves compared to other citrus types. Chinese yuzu cultivars exhibited a more distant genetic relationship with Japanese and Korean yuzu varieties, consistent with their historical dispersal routes (Figure 1A). Kabosu (Citrus sphaerocarpa) originated in China. It is believed that the fruit was introduced to Japan during the Edo period. Legend has it that a Japanese physician purchased the fruit from a monk and subsequently cultivated it in Oita Prefecture. It is speculated to be a hybrid of ichang papeda and bitter orange, but its origin is unknown (Tomiyama et al., 2011). Xiang yuan (Citrus wilsonii) is a wild citrus that originated from China. Recent research suggests that Xiang yuan may be a hybrid variety between pummelo and yuzu, with pummelo likely serving as its maternal parent (Demarcq et al., 2021). Kabosu and Xiang yuan are more closely related to Chinese yuzu, possibly due to their shared origin in China and similar parental relatives (Figure 1A). Citrumelo and citrange are hybrids of trifoliate with ‘Duncan’ grapefruit and sweet orange, respectively (Stover et al., 2021). They share common trifoliate parents and thus cluster together phylogenetically (Figure 1A).

Consistent with these genetic relationships, yuzu cultivars with closer evolutionary proximity displayed greater similarity in fruit morphology, particularly in size and color, whereas substantial phenotypic variation was observed among the other citrus types, such as the thicker albedo found in the peel of the Xiang yuan fruit (Figure 1A).

Aroma is a critical quality trait of aromatic acid citrus. To evaluate aroma characteristics among the selected citrus varieties, electronic nose (E-nose) analysis was performed using peel samples from all 11 cultivars. Most citrus varieties exhibited similar aroma response patterns, whereas Kabosu and citrumelo showed distinct aroma profiles (Figure 1B). Principal component analysis (PCA) of the E-nose data further confirmed that Kabosu and citrumelo were clearly separated from the other citrus varieties along the principal components, indicating substantial differences in aroma characteristics (Figure 1C). Among the ten metal oxide sensors, W5S, W1S, W1W, W2S, and W2W showed the strongest responses to citrus aromas (Figure 1B). Consistently, loading analysis demonstrated that these five sensors contributed most strongly to the variation captured by PC1 and PC2 (Figure 1D), suggesting that they are key sensors for discriminating aroma differences among aromatic acid varieties.

To further elucidate the chemical basis of citrus aroma formation, volatile organic compounds (VOCs) were analyzed using headspace solid-phase microextraction coupled with gas chromatography–mass spectrometry (HS-SPME–GC–MS). The results revealed a total of 972 VOCs detected across all citrus samples, comprising 26.23% terpenoid, 18.00% ester, 10.08% ketone, 9.67% heterocyclic compound, 8.64% alcohol, 5.45% hydrocarbons, 5.35% aldehyde, 4.73% acid, 3.09% amine, 2.78% phenol, 2.16% aromatics, 2.06% ether, 1.03% nitrogen compounds, 0.62% sulfur compounds, 0.10% halogenated hydrocarbons. Of these, terpenoids represented the most abundant and diverse class of VOCs (Figure 2A; Supplementary Table S1). The E-nose sensor responses were consistent with VOC composition. Sensor W1W is primarily sensitive to inorganic sulfides and many terpenes, with its elevated response values correlating with high levels of terpenoid VOCs (Figures 1B, 2A). Similarly, W2S is primarily sensitive to most alcohols, aldehydes, and ketones; W2W, primarily sensitive to aromatics and organic sulfides. The elevated response values of W2S and W2W are consistent with the high abundance of the corresponding VOCs (Figures 1B, 2A).

Multivariate statistical analyses were performed to further evaluate differences in the VOCs profiles. In PCA analysis, all varieties of yuzu exhibited relatively clustered distributions. We observed relative differences between the China-native varieties ‘Ziyangxiangcheng’ (ZYXC), ‘Xiecheng’ (XC) and, other yuzu, including ‘Kitou’ (KT), ‘Tada Nishiki’ (TN), ‘yamane’ (YM), ‘Zairai’ (ZR), and ‘Line131’ (L131) (Figure 2B). VOCs of Xiang yuan (XY) and yuzu exhibit similar traits, particularly in relation to KT, TN, YM, ZR, and L131 (Figure 2B). Furthermore, the results of orthogonal partial least squares discriminant analysis (OPLS-DA) were consistent with PCA (Figure 2C). Heatmap and hierarchical cluster analysis of VOCs across all samples revealed that KT, TN, ZR, and L131 exhibited the most similar VOC profiles. Compared to ZYXC and XC, XY exhibited greater similarity to KT, TN, ZR, and L131. Furthermore, the aromatic compounds of yuzu differed from those of Kabosu (KBS), citrumelo (SCM), and Citrange (CA) (Figure 2D).

The sensory impact of VOCs depends on both their concentration and odor threshold, expressed as relative odor activity value (rOAV). VOCs with rOAV > 1 were considered to contribute significantly to citrus aroma. After excluding compounds lacking odor threshold data, 308 key VOCs were identified, consisting of 77 terpenoids, 46 esters, 38 heterocyclic compounds, 33 alcohols, 27 aldehydes, 23 ketones, 21 phenols, 12 aromatics, 8 ethers, 6 acids, 6 hydrocarbons, 4 nitrogen compounds, 4 amines, 3 sulfur compounds (Supplementary Table S2). Based on rOAV values, the top 10 VOCs from each citrus variety were selected, yielding 18 VOCs with consistently high aroma contributions (Figure 3). These compounds impart diverse sensory notes such as fruity, mushroom-like, sulfurous, green, floral, and sweet notes to citrus fruits (Supplementary Figure S1). Among them, 1-p-menthen-8-thiol (rOAV: 4895599−169411622), 3−mercapto−3−methylbutyl formate (rOAV: 1083440−7624089), and 3−octen−2−one (rOAV: 224445−3031934) ranked among the top contributors across all citrus varieties and were particularly enriched in yuzu cultivars (Figure 3). These compounds make a major contribution to the formation of citrus aromas. In addition, we found that Damascone, beta- and 2,4-Undecadienal, 1-Octen-3-one, and 2-Furanmethanethiol, 5-methyl- exhibited highly similar accumulation patterns across the different varieties (Figure 3).

To identify characteristic aroma compounds distinguishing citrus varieties, OPLS-DA based on rOAV values of key VOCs was performed. OPLS-DA is a supervised pattern regression analysis method based on orthogonal partial least squares regression. Through interactive residual validation ANOVA, the explanatory variables (R²) and predictive capability (Q²) of the OPLS-DA model were determined to be 0.969 and 0.787, respectively. Results indicate this model exhibits high classification accuracy and strong predictive capability for characteristic aroma components across different citrus varieties. The variable importance in projection (VIP) values from OPLS-DA can be used to evaluate each variable’s contribution to classification. Compounds with VIP > 1 were considered characteristic compounds. Higher VIP values indicate greater compositional differences among citrus varieties and a stronger contribution to sample differentiation. Results reveal a total of 115 characteristic aroma compounds with VIP values >1 (Supplementary Table S2). Among these, the differential metabolites safranal, cis−7−decen−1−al, borneol, endo−borneol, terpinen−4−ol, L-4-terpineol, and benzyl methyl sulfide exhibited particularly high VIP values (VIP > 2) (Figure 4). These compounds exhibit stronger aromas in ZYXC, XC, and CA, moderate aromas in YM, TN, KT, ZR, L131, and XY, and weaker aromas in KPS and SCM. In addition, 4−heptanol, 3−methyl−, citral, bis(2−chloroethyl) ether are characteristic aroma compounds of KBS and SCM (Figure 4; Supplementary Figure S2). α−farnesene, undecanol−5, and 1−undecanol were more abundant in ZYXC. Camphene, benzene, tertbutyl−, L-camphene, β-thujene, 1-methylformylpyrrole, propanoic acid, pentyl ester, hexanoic acid, ethyl ester, and butanoic acid, 1−methylbutyl ester accumulated more significantly in the CA sample (Figure 4; Supplementary Figure S2). These characteristic VOCs possess numerous flavors, resulting in the diverse aromatic profiles exhibited by different citrus fruits. For instance, the top 30 characteristic VOCs account for the variations in scents such as fresh, herbal, citrus, woody, camphor, fruity, sweet, and green among these citrus varieties (Supplementary Figure S3).

Key VOCs in citrus possess high rOAV values and play a major role in the formation of citrus aroma. Characteristic VOCs contribute to the aromatic differences among various citrus varieties. To elucidate the molecular regulatory mechanisms underlying citrus aroma formation and varietal divergence, we selected three VOCs with the strongest olfactory impact (highest rOAVs) and one VOC showing the greatest variation among cultivars (highest VIP value) for further analysis. Among the VOCs with the strongest aroma intensity, 1-p-menthene-8-thioli, also known as grapefruit mercaptan, is the primary aroma compound found in grapefruit. 3-mercapto-3-methylbutyl formate is a carboxylic acid ester widely used as a flavoring agent, characterized by intense catty, caramel, onion, roasted coffee, and sulfurous notes. 3-octen-2-one is a ketone with a potent earthy and mushroom-like odor, commonly applied in flavor and fragrance formulations. Among the VOCs exhibiting the greatest aromatic variation across citrus varieties, safranal is a monoterpene compound presenting as a pale straw-colored oily liquid with fresh, herbal, and rosemary-like aromas. To elucidate the gene regulatory mechanisms underlying the synthesis pathways of these important VOCs, transcriptome analysis was conducted. Weighted Gene Co-expression Network Analysis (WGCNA) was employed as a systems biology approach to identify modules of highly correlated genes across large-scale gene expression datasets. This method constructs a weighted gene co-expression network, clusters genes into distinct modules, and relates module eigengenes to VOC abundance, thereby enabling the identification of gene networks associated with specific VOC. Moreover, WGCNA facilitates the detection of hub genes within each module, which may serve as core regulators of VOC biosynthesis. To construct the gene co-expression network based on the appropriate scale-free topology criterion, set the soft-thresholding power of β = 14, MinModuleSize = 25, and MergeCutHeight = 0.75 (Supplementary Figure S4). Correlation analysis between module eigengenes and VOC levels revealed that the green module exhibited a strong positive correlation with 1-p-menthene-8-thiol (r = 0.85). The tan module showed high correlations with both 3-mercapto-3-methylbutyl formate (r = 0.90) and 3-octen-2-one (r = 0.89), while the darkgreen module was correlated with safranal (r = 0.70) (Supplementary Figure S5). These results indicate that these gene co-expression modules are closely involved in regulating the biosynthesis of their corresponding VOCs.

Genes with the highest degree of connectivity in essential modules are known as hub genes, and they may have important biological functions related to traits. The co-expression network was visualized using Cytoscape, and genes with the highest connectivity within each key module were identified as hub genes (Supplementary Table S4). As shown in the Figure 5, Cs_ont_5g010460 was identified as the hub gene in the green module and was closely associated with the regulation of 1-p-menthene-8-thiol biosynthesis. The tan module’s hub gene, Cs_ont_5g018210, is highly associated with the synthesis of 3-Mercapto-3-methylbutyl formate and 3-Octen-2-one. The darkgreen module’s hub gene, Cs_ont_3g018770, plays a crucial role in safranal synthesis. To further demonstrate the close relationship between hub genes and their corresponding key VOCs, Pearson correlation analysis was conducted. The results showed a significant positive correlation between the expression level of Cs_ont_5g010460 and the accumulation of 1-p-menthene-8-thiol (r = 0.736). Similarly, Cs_ont_5g018210 exhibited strong positive correlations with 3-mercapto-3-methylbutyl formate (r = 0.865) and 3-octen-2-one (r = 0.916). In addition, Cs_ont_3g018770 showed a significant correlation with safranal content (r = 0.691) (Figure 6). These results indicate that the expression levels of these hub genes are closely associated with the accumulation patterns of their respective VOCs, supporting their potential regulatory roles in citrus aroma formation.

To further validate the expression patterns of the hub genes, RT-PCR analysis was conducted. The results showed that the expression levels of Cs_ont_5g010460, Cs_ont_5g018210, and Cs_ont_3g018770 determined by RT-PCR were highly consistent with the RNA-seq data, confirming the reliability of the transcriptome analysis for quantifying gene expression (Figure 7). Cs_ont_5g010460 exhibited significantly higher expression in XC than in the other citrus varieties, which was consistent with the elevated accumulation of 1-p-menthene-8-thiol in XC. Cs_ont_5g018210 showed relatively low expression levels in CA and SCM, corresponding to the reduced accumulation of the key VOCs 3-mercapto-3-methylbutyl formate and 3-octen-2-one in these varieties. Cs_ont_3g018770 was highly expressed in ZYXC, XC, and CA, displaying an expression pattern similar to the accumulation trend of the characteristic VOC safranal (Supplementary Figures S2, S3, S7).

Functional annotation further revealed the putative biological roles of these hub genes. Cs_ont_5g010460 encodes ovate family protein 16 (OFP16), a plant-specific transcriptional regulator that functions as a transcription factor or repressor and plays an important role in modulating downstream gene expression. Cs_ont_5g018210, encoding SKP1-interacting partner 16 (SKIP16), is an F-box protein and a core component of the SCF (SKP1–Cullin–F-box) ubiquitin ligase complex. This protein mediates the selective recognition and ubiquitination of target substrates, thereby regulating protein turnover through the 26S proteasome pathway. Cs_ont_3g018770 encodes a zeta-carotene desaturase-like (ZDS-like) protein, whose homologs catalyze the introduction of double bonds into colorless zeta-carotene, leading to the formation of colored carotenoids such as lycopene (Chen et al., 2017). This enzyme is essential for the biosynthesis of metabolites, plant development, and the enhancement of crop nutritional quality. Based on their potential biological functions, these genes may control the biosynthesis of these key VOCs by regulating gene transcription, protein levels, and catalytic reactions.

This study employed the ‘Ziyangxiangcheng’ (ZYXC), ‘Xiecheng’(XC), ‘Kitou’ (KT), ‘Tada Nishiki’ (TN), ‘Zairai’(ZR), ‘Yamane’ (YM), ‘Line131’ (L131) yuzu (Citrus junos) cultivars. Additionally, ‘Kabosu’ (KBS) (Citrus spaerocarpa), ‘Swingle’ citrumelo (SCM), Xiang yuan (XY) (Citrus wilsonii), and citrange (CA) were included. All citrus fruits were cultivated at the experimental field of Zhejiang Citrus Research Institute (E: 121.9, N: 28.38), with uniform cultivation practices applied throughout. By 1st December, all fruits had reached maturity with consistent ripeness levels. Citrus peel samples were randomly collected, immediately flash-frozen using liquid nitrogen, and stored at -80 °C for subsequent experiments. Three biological replicates were established for each cultivar.

The PNE3 electronic nose system (Airsense Analytics, Schwerin, Germany) was employed to analyse volatile compounds released from citrus peel samples. This olfactory system incorporates 10 different semi-conductive metal oxide semiconductor sensors responsive to distinct chemical classes, enabling comprehensive characterization of aroma profiles. For analysis, approximately 3g of citrus peel is placed within a sealed vial. The vial is then sealed and positioned within the sensor chamber of the PNE3 electronic nose system. The system operates at 25 °C with a carrier gas ambient air flow rate of 200 mL/min. Measurement parameters were based on the methodology described by Chen et al (Chen et al., 2020). Steady-state response values were recorded for 57 to 60 s during each test, with data analysis performed using WinMuster software (version 1.6.2).

VOCs were extracted using headspace solid-phase microextraction (HS-SPME). Approximately 0.5 g of ground sample was transferred into a 20 mL glass headspace vial (Agilent Technologies, Palo Alto, CA, USA) and sealed with a polytetrafluoroethylene (PTFE)/silicone septum. Sodium chloride was added to enhance VOC release by salting-out effects. The vial was equilibrated at 60 °C for 5 min with agitation. VOCs were extracted by exposing a 120 μm divinylbenzene/carboxen/polydimethylsiloxane SPME fiber to the headspace of the sample for 15 min at 60 °C. After extraction, the SPME fiber was immediately inserted into the injection port of the gas chromatograph for thermal desorption.

VOC analysis was performed using a gas chromatography–mass spectrometry (GC–MS) system (Agilent 7890 GC coupled with a 5977 MSD). Separation was achieved on a DB-5MS capillary column (30 m × 0.25 mm × 0.25 μm, 5% phenyl-polymethylsiloxane) (Agilent Technologies, Palo Alto, CA, USA). The GC oven temperature program was set as follows: initial temperature at 40 °C for 3.5 min, increased to 100 °C at 10 °C min⁻¹, and finally to 280 °C at 25 °C/min (held for 5 min). Helium was used as the carrier gas at a constant flow rate of 1.2 mL min⁻¹. The MS was operated in electron impact (EI) mode at 70 eV, with a mass scan range of m/z 35–500. The ion source and quadrupole temperatures were set to 230 °C and 150 °C, respectively. VOCs were detected using selected ion monitoring (SIM) mode.

Relative quantification of VOCs was performed based on peak area normalization, and the results were expressed as the relative percentage of total ion chromatogram (TIC) peak areas. The normalized VOC data were subjected to multivariate statistical analyses, including principal component analysis (PCA) and Orthogonal Partial Least Squares Discriminant Analysis (OPLS-DA), using MetaboAnalyst package in R (version 4.5.1) (Pang et al., 2024). Metabolites with a variable importance in projection (VIP) value > 1.0 were considered as key VOCs.

Total RNA was extracted from citrus peels using the TRIzol reagent (Invitrogen, USA) following the manufacturer’s instructions. RNA concentration and purity were assessed using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). RNA samples with an RNA integrity number (RIN) ≥ 7.0 were used for library construction. Each experimental group consisted of three independent biological replicates. Libraries were amplified by PCR and sequenced on a NovaSeq X Plus platform (Illumina, Inc., San Diego, CA, USA) to generate 150-bp paired-end reads by Metware Biotechnology Co., Ltd. (Wuhan, China). Raw sequencing reads were processed to remove adapter sequences, low-quality reads, and reads containing poly-N using fastp (version 0.23.4) (Chen et al., 2018). Clean reads were mapped to the Citrus sinensis reference genome using HISAT2 with default parameters (version 2.2.0) (Kim et al., 2015). Uniquely mapped reads were retained for calculating read counts. Gene expression levels were quantified as fragments per kilobase of transcript per million mapped reads (FPKM) using featureCounts (Liao et al., 2013).

Weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA package in R to identify gene modules associated with key VOCs (Langfelder and Horvath, 2008). Genes with low expression levels (FPKM < 1) were filtered out prior to network construction to reduce noise. An expression matrix of the remaining genes was used to calculate pairwise Pearson correlation coefficients. A soft-thresholding power (β) was selected based on the scale-free topology criterion using the pickSoftThreshold function. The adjacency matrix was constructed and subsequently transformed into a topological overlap matrix (TOM), which measures the network connectivity between genes. Hierarchical clustering was performed based on the TOM dissimilarity measure, and gene modules were identified using the dynamic tree cut method with a minimum module size of 25 genes. Highly similar modules were merged using a threshold of 0.25. Each module was summarized by its module eigengene (ME), representing the first principal component of the gene expression matrix within a module. Module–trait relationships were evaluated by correlating ME with sample traits (VOC levels) using Pearson correlation analysis. Modules with high correlation coefficients and significant p-values were considered biologically relevant. Hub genes within key modules were identified based on highest degree of connectivity.

Residual gDNA was removed and cDNA synthesized using the Hifair^® III 1st Strand cDNA Synthesis SuperMix for qPCR (YEASEN, Shanghai, China). 2×HS Taq Universal SYBR Green Qpcr Master Mix (SAIPUBIO, Shenzhen, China) and BIOER FQD-96A (BIOER, Hangzhou, China) were used for qPCR experiments. Actin3 was used as internal control gene in citrus. Three biological replicates were performed for each experiment and relative gene expression was calculated following the 2 ^-ΔΔCT method. Values are presented as mean values. Primer information is provided in Supplementary Table S5.