Saccharomyces cerevisiae EC1118a used in this study was provided and preserved by the State Key Laboratory of Microbiology, School of Bioengineering, Qilu University of Technology. The plasmids and primers used in this study are listed in Supplementary Tables S1, S2, respectively. Enzyme preparations used in this study were purchased from Vazyme Biotechnology Co., Ltd. (Nanjing, China).

Saccharomyces cerevisiae was cultured in YPD medium (containing 1% yeast extract, 2% peptone, and 2% glucose); 200 μg/mL of G418 (geneticin, Coolaber, Beijing, China) was added to YPD medium to screen for KanMX-resistant transformants in recombinant strains. YPG medium (1% yeast extract, 2% peptone, and 2% galactose) was used to select strains with KanMX resistance gene eliminated. The grape juice culture medium was prepared from Chardonnay grapes that were destemmed, crushed, and pressed to extract the juice. After centrifugation, the juice was sterilized and set aside for later use (total sugar, 206 g/L; total acid: 5.8 g/L; pH 3.6). SO₂ (in the form of K₂S₂O₅) at 30, 60, and 90 mg/L was added to the grape juice culture medium. Each concentration has three parallel groups, and equal amounts of yeast cells were inoculated into grape juice culture medium. At the end of fermentation, the concentrations of isobutanol, isoamyl alcohol, and 2-phenylethyl alcohol were determined. A UV spectrophotometer was used to measure optical density at 600 nm to determine the biomass concentration and draw a standard curve.

Yeast cells cultured in grape juice medium containing various concentrations of SO₂ were harvested during the logarithmic phase of growth, centrifuged at a low temperature (8,000 × g, 4°C), and immediately frozen in liquid nitrogen. Total RNA was extracted using the TRIzol reagent. The quality and integrity of RNA nucleic acid samples were checked using agarose gel electrophoresis and Agilent 2,100 Bioanalyzer (Agilent Technologies Inc., CA, United States). Measuring the purity and concentration of RNA nucleic acid samples using the NanoDrop™ One/OneCsystem (Thermo Fisher Scientific, MA, United States).

After enriching the eukaryotic mRNA with a poly A tail using magnetic beads with oligo (dT), the mRNA is fragmented using ultrasound and reverse transcribed into the first strand of cDNA (Cui et al., 2023). The RNA chain was degraded with RNase H, and the second strand of cDNA was synthesized using dNTPs as the raw material in the DNA polymerase I system. Purified double-stranded cDNA was end-repaired, A-tailed, and connected to sequencing adapters (Peng et al., 2023). AMPure XP beads were used to screen approximately 200 bp of cDNA. PCR amplification was performed. Finally, AMPure XP beads were used to purify the PCR products and obtain a library (Duan et al., 2022). The sequencing strategy was PE150, and the sequencing platform was Illumina Novaseq 6,000. Sequencing of cDNA libraries was done using an Illumina sequencing platform by GeneDenovo Biotechnology Co., Ltd. (Guangzhou, China).

Use Fastp (v0.18.0) to perform quality control on the raw reads off the machine, and filter low-quality data to obtain clean reads. Use the short read alignment tool Bowtie2 (v2.2.8) to align the clean reads to the ribosomal RNA database, remove the aligned ribosomal reads without allowing mismatches, and use the remaining unmapped reads for subsequent transcriptomes analyze. HISAT2 (v2.2.4) (Kim et al., 2015) was used to align the sequences obtained by paired-end sequencing to the reference genome (Ensembl_release47).¹ Stringtie was used to reconstruct transcripts based on the comparison results, and RSEM (v1.3.1) was used to quantify expression abundance and variation through FPKM values (number of transcript fragments per kilobase/million mapped reads).we used R language gmodels (v2.18.1)² to carry out principal component analysis (PCA) to study the distance relationship between samples through dimensionality reduction. Differential gene expression analysis was performed using the DESeq2 (Love et al., 2014) software. Genes meeting the criteria FDR1 < 0.05 and |log2FC| > 1 were considered as significant DEGs.

After slowly thawing 6 sets of parallel samples and quality control samples (QC Dunn et al., 2011) at 4°C, take an appropriate amount of samples and add them to the pre-cooled solution (methanol: acetonitrile: aqueous solution = 2:2:1, v/v) Mix well, sonicate at low temperature for 30 min, then let stand at −20°C for 10 min, and centrifuge at low temperature for 20 min (14,000 × g, 4°C). Take the supernatant and dry it under vacuum. During mass spectrometry analysis, add 100 μL acetonitrile aqueous solution (acetonitrile: aqueous solution = 1:1, v/v) to reconstitute. After centrifugation for 15 min (14,000 × g, 4°C), the supernatant was taken as the loading solution. Equal volumes from the samples were used to determine the instrument status and equilibrium chromatography prior to injection. Mass spectrometry was used to evaluate system stability during the experimental process. LC–MS/MS (Triple quadrupole: Q-TOF) analysis was performed to quantify and characterize metabolites, followed by data quality evaluation and analysis.

Preprocess the original data and use ProteoWizard software to convert the original data into. MzML format. The XCMS program was then used to extract the data and perform data integrity checks and normalization of the total peak area of the data to ensure parallel comparisons between samples and metabolites. Two ionization methods, positive ion mode (POS) and negative ion mode (NEG), are used for data quality control, and the data can be analyzed with higher coverage and better detection results. The R language pheatmap (Wu et al., 2022) was used to standardize the data z-score and perform cluster analysis, combined with the VI p value of the multivariate statistical analysis PLS-DA (Bylesjö et al., 2006) and the univariate statistical analysis T-test p value to screen for significantly different metabolites between different comparison groups (Saccenti et al., 2014). Utilizing the rich pathway data from differential metabolites, we constructed a metabolic relationship network diagram to elucidate the interactions between different pathways, aiding in the identification of core metabolic pathways. The threshold for significant differences was VIP ≥ 1 and a t-test <0.05 in the PLS-DA model.

The Gene Ontology (GO) database³ was used for the GO functional analysis. GO functional classification annotations of DEGs help determine the functions of these genes, while GO functional significance enrichment analysis determines the main biological functions of DEGs. The Kyoto Encyclopedia of Genes and Genomes (KEGG)⁴ was used to analyze large-scale molecular datasets generated from molecular-level information, especially through high-throughput experimental techniques such as genome sequencing, to understand the advanced functions and utility of biological systems (Wang et al., 2021).

To screen and obtain a set of correlated genes and metabolites that influence sample grouping and analyze correlation characteristics, three model analyses were performed based on gene expression and metabolite abundance data. First, a pathway functional model was used to query the KEGG metabolic pathways shared by genes and metabolites and analyze their correlation characteristics (Cho et al., 2016). Second, gene expression and metabolite abundance data were used to construct a bidirectional orthogonal projections to latent structures (O2PLS) model (el Bouhaddani et al., 2016); through model prediction, we obtained a combined analysis of related gene and metabolite sets. Finally, when the sample group was ≥3, a correlation coefficient model was used to calculate the Pearson correlation coefficient between gene expression and metabolite abundance (Bartel et al., 2015; Hamanishi et al., 2015). The model analysis was visualized as a heat map and network diagram.

Using the yeast genome as a template, the upstream and downstream homologous fragments of the target gene were PCR amplified. Additionally, the plasmid pUG6 was used as a template to PCR amplify the Loxp-KanMX-Loxp resistance fragment to construct the knockout component. These three knockout component fragments were transformed into the parental strain using the lithium acetate/polyethylene glycol (PEG) method (Gietz and Schiestl, 2007) to achieve homologous recombination between the transformed fragments and the yeast genome (Supplementary Figure 1). Transformants were initially screened on YPD agar containing 200 μg/mL G418, and PCR was used to confirm the accurate integration of the KanMX knockout cassette from a single colony, which was subjected to two strain activations. Subsequently, lithium acetate/ PEG method was used to transform the plasmid pSH69 into the abovementioned colonies containing the correct integration. Plasmid transformation was verified. Colonies containing the plasmids were induced and cultured on YPG agar medium. KanMX-resistant fragments were removed through the expression of Cre recombinase induced by galactose in the medium. Subsequently, the plasmid pSH69 was lost through subculture, and the recombinant strain with the desired knockout was obtained.

A GC–MS method was used to determine the content of higher alcohols. An 8 mL sample was transferred into a 15 mL sample bottle with 1.5 g NaCl and 100 μL butyl acetate (internal standard) and preheated to 45°C for 10 min using a magnetic stirrer/heating stage. The material was then extracted with a DVA/CAR/PDMS extraction head, and analyzed using GC–MS. Each sample was analyzed three times. This was done on an InertCap capillary column (30 mm × 0.25 mm × 0.25 μm) using the temperature protocol: hold at 40°C for 3 min, heat to 130°C at 2°C/min, and heat to 220°C at 10°C/min retention for 4 min; the injector and detector temperature 250°C, and no split injection; electron impact ion source, electron energy 70 eV, ion source temperature 230°C, and full scan mode. Based on the chromatographic retention time and mass spectrometry information of the standard product, quantitative analysis was performed by comparing the qualitative and internal standard methods with reference to the standard spectral library (NIST17).

To perform fluorescence reverse transcriptase real-time quantitative PCR (qPCR), RNA was extracted using the SPARK Easy Yeast RNA, and Script II RT Plus Kits (Sparkjade Biotechnology Co., Ltd., Shandong, China) was used to reverse transcribe the RNA. Taq Pro Universal SYBR qPCR Master Mix was obtained from Vazyme Biotechnology (Nanjing, China). The PCR procedure consisted of an initial pre-denaturation at 95°C for 30 s, followed by amplification consisting of 40 cycles of denaturation at 95°C for 5 s, annealing and polymerization at 60°C for 30 s, and a melting curve stage at 95°C for 15 s and at 60°C for 1 min. The ubiquitin-conjugating enzyme E26 gene UBC6 was used as the internal reference gene, and the 2^−ΔΔCt method was used for quantitative analysis (Teste et al., 2009).

All experiments were conducted in at least three parallel groups, and statistical analysis was performed in GraphPad Prism 9.5.1 (San Diego, California, United States). The results were expressed as mean and standard deviation (SD), and were analyzed using one-way analysis of variance (ANOVA). There was a significant difference, and the statistical significance level was set at p value <0.05.

The generation of higher alcohols was correlated with SO₂ concentration (Figure 2A). At 30, 60, and 90 mg/L SO₂, the total alcohol content was 91.513, 113.677, and 98.649 mg/L, respectively. Compared with the 30 mg/L SO₂, total high alcohol content increased by 16.892 and 7.797% at the 60 and 90 mg/L SO₂, respectively; this involved an increase in the amounts of isobutanol by 31.064 and 5.964%, isoamyl alcohol by 16.473 and 7.163%, and 2-phenylethyl alcohol by 16.824 and 11.064%, respectively. The yield of higher alcohols showed an overall trend of being highest at 60 mg/L SO₂, and showed a decline at 90 mg/L SO₂. As the initial SO₂ concentration increased, yeast growth was inhibited (Figure 2B). The inhibitory effect was evident at 90 mg/L, and the duration of the lag phase was significantly prolonged with increasing SO₂ concentrations.

Based on the transcriptomic analysis results, PCA showed reliable repeatability in the three parallel samples and identified 2,766 DEGs, of which 937 were upregulated and 1829 were downregulated (Figures 3A,B). GO and KEGG identified a total of 6,510 GO terms that were enriched in biological process (BP; 34.29%), molecular function (MF; 30.38%), and cellular component [CC; 35.33%] (Figure 3C). According to the results of the BP analysis, with the increase in SO₂ concentration, the number of upregulated and downregulated DEGs also increased. Metabolic processes related to the synthesis of higher alcohols include pyruvate metabolism that is catabolized into alcohol through the Ehrlich pathway (GO:0006090), glycolysis (GO:0006096), amino acid metabolism (GO:1901605), oxidation–reduction process (GO:0055114), redox coenzyme metabolism (GO:0006733), and oxidative phosphorylation (GO:0006119). The results of the MF analysis relating to the synthesis of higher alcohols included phosphoenolpyruvate carboxykinase (ATP) activity (GO:0004612), aldehyde dehydrogenase (NAD) activity (GO:0004029), pyruvate decarboxylase activity (GO:0004737), and transaminase activity (GO:0008483). KEGG enrichment analysis revealed that the pathway classes with significantly enriched differential genes mainly included metabolic pathways, carbon metabolism, amino acid biosynthesis, glycolysis, and gluconeogenesis (Figures 4A,B). The pathways related to the synthesis of higher alcohols mainly involve carbon metabolism (ko01200), pentose phosphate pathway (ko00030), pyruvate metabolism (ko00620), glycolysis/gluconeogenesis (ko00010), and the valine, leucine, and isoleucine anabolic pathways (ko00290). In addition, the regulation of DEGs examined by RT-qPCR was consistent with that observed using RNA-seq data (Supplementary Figures 2A–D).

We conducted a trend analysis on three groups of gradient samples with varying initial SO₂ concentrations (30, 60, and 90 mg/L). Clustering was based on the gene expression patterns of the samples (the shape of the expression curve in multiple stages), in order to facilitate the selection of representative gene sets that met specific biological characteristics in the clustering result (Zeng et al., 2019). The overall findings of the trend analysis are shown in Supplementary Figure 3. A total of seven change trends was observed. Based on the change trend observed in the cumulative concentration of higher alcohols (total higher alcohol) and the gene expression function, the following three gene expression trends were analyzed. Twenty-eight genes exhibited a gene expression trend that initially increased and subsequently decreased in relation to the SO₂ content, including ADH4, ADH5, GDH2, FBA1, LEU2, TIR1, and TPI1. Eighty-nine genes showed an initial downward and a subsequent upward trend, these included ALD6 and SOL4. GPD1 and PCK1 showed a continuous downward trend.

The metabolites identified by metabolomics and related to the synthesis of higher alcohols are shown in Supplementary Table S3. Metabolomic analysis identified 946 differential metabolites in POS mode, of which 416 were upregulated and 530 were downregulated (Figure 5A). In total, 943 differential metabolites were identified in NEG mode, of which 485 were upregulated and 458 were downregulated (Figure 5B). KEGG analysis was performed on different metabolites in the POS and NEG modes. The results showed that the different metabolic pathways in POS mode (Figure 5C) mainly included glycerophospholipid metabolism (ko00564), glutathione metabolism (ko00480), and the biosynthesis of alkaloids derived from the shikimate pathway (ko01063). The main differential metabolites were acetylcholine (C01996), acetyl coenzyme A (C00024), and glutathione (C00051). In NEG mode (Figure 5D), the affected metabolic pathways included the pentose phosphate pathway (ko00030), amino acid metabolism (ko01230), metabolic pathway (ko01100), and carbon metabolism (ko01200). The main differential metabolites include S-adenosylmethionine (C00019), pyruvate (C00022), α-ketoglutarate (C00026), glycerol (C00116), threonine (C00188), isoleucine (C00407), and glutamine acid (C00302).

In addition, the core pathways in POS mode (Figure 5E) were the glycerophospholipid metabolism pathway, whose main related pathways are ether lipid metabolism and glycolysis/gluconeogenesis, and metabolic pathways; the glutathione metabolism pathway, whose main related pathways are arginine biosynthesis, cysteine and methionine metabolism; biosynthesis of alkaloids derived from the shikimate pathway, whose main related pathways are phenylalanine metabolism, citrate cycle, phenylalanine, tyrosine and tryptophan biosynthesis; thiamine metabolism, whose main related pathways are glycine, serine and threonine metabolism, cysteine and methionine metabolism. The core pathway in NEG mode (Figure 5F) was the metabolic pathway, with its main related pathways being the pentose and glucuronate interconversions, the pentose phosphate pathway, fatty acid biosynthesis, glycine, serine, and threonine metabolism, glyoxylate and dicarboxylate metabolism, and valine, leucine, and isoleucine biosynthesis.

The DEGs obtained from transcriptomics and differential metabolites obtained from metabolomics were subjected to KEGG pathway analysis. The most common metabolic pathways include glycolysis/gluconeogenesis, carbon metabolism, amino acid biosynthesis, tryptophan metabolism, and the TCA cycle. The O2PLS model was used to perform a correlation analysis on the two sets of scientific data. The O2PLS load diagram (Supplementary Figures 4A–C) of the joint part (the part that is highly related to the transcriptome and metabolome) and the two sets of scientific correlation load diagrams (Supplementary Figures 5A–C) were used to obtain a combined analysis of the correlated gene and metabolite sets. The correlation between the top 250 DEGs and differential metabolites with correlation coefficients and absolute values greater than 0.5 is represented by a heat map (Figure 6A) and a network diagram (Figure 6B). In the correlation analysis, differential metabolites related to the synthesis of higher alcohols include pyruvate, valine, leucine, isoleucine, tyrosol, and glycerol, and related genes include CDC19 (pyruvate kinase), CYB2 (L-lactate dehydrogenase, cytochrome b2), ENO2 (phosphopyruvate hydratase), ALD4 (mitochondrial aldehyde dehydrogenase), ALD6 (cytosolic aldehyde dehydrogenase), ARO10 (phenylpyruvate decarboxylase), PDCs (pyruvate decarboxylase), ADH4 (alcohol dehydrogenase isoenzyme type IV), ADH5 (Alcohol dehydrogenase isoenzyme V), BDH2 (alcohol dehydrogenase), TPI1 (triose phosphate isomerase), HXT2 (glucose transporter), FBA1 (fructose 1,6-bisphosphate aldolase), and GDH2 (NAD⁺-dependent glutamate dehydrogenase).

First, DEGs and differential metabolites were obtained through single-omics analyses such as PCA and differential enrichment analysis. Based on the changing characteristics of gene expression at different omics levels, DEGs and different metabolites were subjected to O2PLS analysis and various correlation analyses to explore the potential molecular connections between different omics levels. Starting from the target differential metabolites related to the anabolism of higher alcohols, the associated differential genes were obtained and combined using existing research and trend analysis. Ten target DEGs (ADH4, GDH2, SER33, PDC6, LEU2, CHA1, ARO10, TIR1, BDH2, and ADH5) were chosen for the knockout experiments. Single-gene knockout strains were constructed to explore the effect of their knockout on the production of higher alcohols. The total alcohol synthesis in each knockout strain was lower than that in the original strain (Figures 7A,B). Among them, the EC1118a-ADH4, EC1118a-SER33, and EC1118a-ADH5 knockout strains showed a more significant reduction effect, with the total higher alcohol content reduced by 16.479, 15.620, and 14.771%, respectively. For isobutanol, only the EC1118a-SER33 knockout strain showed a decreasing trend and the effect was significant with a reduction of 71.346%; whereas for isoamyl alcohol, the EC1118a-ADH4 and EC1118a-TIR1 knockout strains showed a significant reduction of 21.003 and 17.334%, respectively. Regarding 2-phenylethanol, the effects of EC1118a-GDH2, EC1118a-CHA1, and EC1118a-ARO10 knockout were more significant, with reductions of 25.198, 18.234, and 18.016%, respectively.