C. utilis was cultured at 30°C for 24 h, and maintained in the stationary phase using the YEPD medium containing 2% peptone, 2% glucose, and 1% yeast extract. The genomic DNA was extracted by phenol-chloroform extraction and ethanol precipitation according to the standard procedures. The DNA samples were then subjected to sequencing using nanopore sequencing technology (Kono and Arakawa, 2019). The third-generation nanopore sequencing technology offers several distinct advantages. Its long reading length facilitates the assembly of large genomes and significantly improves genome integrity. Additionally, it provides real-time sequence information, but it does have a drawback: the error rate of individual reads is high, necessitating repeated sequencing for error correction. In comparison, Illumina HiSeq sequencing is more accurate than nanopore sequencing, but it has much shorter read lengths. However, by using Illumina HiSeq data to correct nanopore sequencing data, the error rate of the third-generation nanopore sequencing can be reduced to a level similar to second-generation sequencing data. Combining nanopore sequencing technology with Illumina HiSeq sequencing technology leads to more accurate and in-depth data mining results than using second-generation sequencing alone. The experimental procedures, including sample quality detection, library construction, library quality detection and library sequencing, were conducted according to the standard protocol provided by ONT (Karamitros and Magiorkinis, 2018). For genome assembly, the sequencing reads obtained from nanopore sequencing were analyzed with Canu v1.5, and the merged genomes were error-corrected by mapping them to the Illumina Hi-seq reads using Pilon (Walker et al., 2014). To assess the quality of the assembled genome, the assembled genome data was mapped using BWA. Finally, a BUSCO test was applied to evaluate the gene content of the assembled genome (Seppey et al., 2019).

Genome component prediction of C. utilis involved the identification and classification of various DNA sequences, including repeating DNA sequences, non-coding RNA, pseudogenes, and coding genes. To identify and classify repeating DNA sequences, Repeat-Masker (version 4.0.6) was utilized. To determine the tRNA genes in the genome of C. utilis, the tRNAscan-SE tool was employed. The microRNA and rRNA genes were identified using Infernal 1.1 and Rfam databases, respectively (Nawrocki and Eddy, 2013). Pseudogenes, which are non-functional copies of genes, were predicted using GeneWise based on characteristics such as the absence of a promoter, premature stop codons, or frameshift mutations. Coding genes were retrieved using a combination of five predictors and homology-based gene prediction with GeMoMa (version 1.3.1) (Keilwagen et al., 2019). To generate consensus gene models, de novo gene predictions and protein alignments were integrated using the Evidence Modeler (version 1.1.1) tool.

Genome annotation and gene function prediction were executed by aligning to the sequences in seven databases as follows: Pfam, TrEMBL, KOG, Nr, KEGG, Swiss-Prot, and GO (Wood et al., 2020). The membrane transport proteins were predicted by TCDB. The factors associated with pathogen-host interactions were analyzed with PHI. The identification of carbohydrate-active enzymes was carried out through the CAZymes database (Tuveng et al., 2019). To ensure biological significance, the best alignment result was selected as annotation. The E-values were calculated, and a threshold with an acceptable E-value (E-value < 10⁵) was then selected. The annotated proteins were then analyzed for GO annotations (biological process, cellular component and molecular function) using WEGO (Web Gene Ontology Annotation Plot) tool (Ye et al., 2018).

With an understanding of the intrinsic genomic diversity of the strain, transcriptome analysis was performed on both wild-type C. utilis (WT) and δ-zein gene genetically engineered C. utilis (RCT). After incubation with YEPD for 24 h, total RNA was extracted and purified using TRIzol reagent according to the manufacturer’s protocol. The RNA integrity of RNA samples was assessed using the Bioanalyzer 2,100 system, and then confirmed by denaturing agarose gel electrophoresis. The transcriptome library was constructed using the NR603-VAHTSTM Total RNA-seq (H/M/R) Library Prep Kit for Illumina. Then, a 2 × 150 bp paired-end (PE150) sequencing was conducted on the Illumina Novaseq^™ 6,000 platform according to the manufacturer’s protocol (Modi et al., 2021). After removing low-quality bases and undetermined bases from the sequencing reads, the reads were mapped to the C. utilis genome using the HISAT2 software. To reconstruct a comprehensive transcriptome, all transcriptomes from the different samples were merged using the gffcompare software. Once the final transcriptomes were generated, the expression levels of all transcripts were estimated using StringTie and ballgown. These tools allow for the quantification of gene expression by calculating Fragments Per Kilobase of transcript per Million mapped reads (FPKM). After normalization to a reference gene, the expression levels of genes in the experimental group RCT were compared with those in the control group WT. For experiments with three biological replicates, differential expression analysis was performed using the DESeq2 software package (Shahjaman et al., 2020). To identify significantly differentially expressed mRNAs, a criterion of fold change (FC) ≥ 2 and a false discovery rate (FDR) < 0.01 were applied. Once the DEGs were identified, they were subjected to GO and KEGG pathway enrichment analyses using the clusterProfiler package in R. Additionally, GSEA was performed using the clusterProfiler package in R (Fornes et al., 2020).

Following the RNA-Seq data analysis, gene co-expression network analysis was performed using the WGCNA package in R (Liang et al., 2023). As a result, specific gene modules of interest were identified and further investigated. To gain insights into the functional significance of the gene modules identified by WGCNA, GO and KEGG analyses were conducted.

To assess phenotypic divergence between the WT and RCT C. utilis, an Omnilog™ phenotype microarray (PM) system was employed. OmniLog phenotype microarray technology, developed by Biolog, was utilized to monitor and compare the differences in phenotype under various growth conditions (Karlsson et al., 2023). The carbohydrate/carbon source utilization profiles of WT and RCT C. utilis were examined by PM1 and PM2 phenotype microarrays containing 190 carbon sources. The phenotype microarray analyses were conducted following the manufacturer’s instructions, with two biological replicates for each growth condition (Luque-Sastre et al., 2021). In all the growth experiments, both WT and RCT C. utilis were precultured with YPDA agar medium for 24 h. The monoclonal strains cells were removed from the YPDA agar plate using a sterile swab, and then transferred into a sterile capped tube containing 20 mL of nutrient supplement (Table 1). To achieve a uniform suspension of the cell sample, the cell suspension was gently stirred in the capped tube using the swab. After stirring, the turbidity of the suspension was examined and the cell concentration was adjusted to achieve 62% T (transmittance) in the biolog turbidimeter. In PM1 and PM2 inoculating fluids, 0.50 mL of 62% T cell suspension was added to the vial containing 23.5 mL of IFY-0 (1.2×) + Dye mix D (75×) + D-glucose (24×). In PM4 inoculating fluids, 1.50 mL of 62% T cell suspension was mixed with 70.5 mL of IFY-0 (1.2×) + Dye mix D (75×) + D-glucose (24×). Then, 100 μL of the suspension was added to each well (Morrison et al., 2017). PM1–2 and PM4 were cultured at 30°C for 72 h in an OmniLog automated incubator/reader. The digital imagery of this instrument was used to track and monitor changes in the respiration of the strains cultures in each well over the incubation period. Meanwhile, the reduction of tetrazolium dye by respiring cells was measured in each well every 15 min using the OmniLog system (Cruz et al., 2021). Cellular respiration activity was evaluated based on the area of a region bounded by a color development time-series. The cells grown without carbon source were used as a control group. The results were analyzed using the Omnilog PM software (Biolog) according to the manufacturer’s instructions (Vaillant et al., 2022). Only the most significant (p-value < 0.05) differences in metabolic activities under the experimental conditions of this study were shown.

The RNA-seq data were validated by qRT-PCR. Four genes were chosen from the DEGs for data confirmation. The concentration and purity of RNA samples were assessed using the NanoDrop 2000 system. To quantify gene expression in real time, specific oligo nucleotide primer pairs were designed against target sequences in the C. utilis genome. The designed primer pairs, which had nearly identical annealing temperatures, were subjected to SYBR Green-based qRT-PCR analysis. The thermal cycling conditions were as follows: an initial heat-denaturing step at 95°C for 15 s, followed by 40 cycles of amplification at 95°C for 30 s and 60°C for 1 min. After amplification, the melting curves of the PCR products were determined by gradually increasing the temperature from 60 to 95°C. To ensure the reliability of the results, all samples were run in triplicate. The housekeeping gene glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an endogenous control. The relative expression levels of target genes were calculated using the 2−ΔΔCT method (Yasin et al., 2023). The qRT-PCR data was compared with the transcriptome study. All data are presented as mean ± SD for three replicates. The PCR primers used for qRT-PCR analysis are listed in Supplementary Table S1.

The genome of WT C. utilis was sequenced by Oxford Nanopore Technologies (ONT) using the Illumina Hi-seq 2,500 platform. This study generated 7,439,441,967 bp clean data with a depth of 568.63×. The final genome assembly was 13,082,961 bp in size, which consisted of scaffolds (≥1 kb) and contigs (≥1 kb) with an N50 of 2,124,471 bp. The GC content of the assembled C. utilis genome was 44.53% (Table 2). The reference genome employed for genome sequencing is that of Cyberlindnera jadinii. The statistical analysis was performed based on the length of the reads and the accumulated amount of data, which allowed for meaningful interpretation and representation of the data through statistical graphs and charts.

The sequencing reads were individually assembled into the contigs with Canu (version 1.5) assembler (Yao et al., 2020). The quality of the assembly was evaluated using two methods. Firstly, the de novo assembled genome was aligned to the Illumina Hi-seq data using BWA software, and the analysis revealed that over 92.94% of the sequences were successfully mapped (Supplementary Table S2). Secondly, the assembly was examined for the presence of conserved core eukaryotic genes using the Benchmarking Universal Single-Copy Orthologs (BUSCO) datasets, and it was found to contain 279 of these genes (Simão et al., 2015), with a BUSCO completeness score of 96.21%. The results obtained from the aforementioned analyses provided strong evidence that the genome assembly of C. utilis was of high quality.

Further comprehensive analysis identified a total of 162,803 bp of repetitive sequences within the assembled genome of C. utilis. These repetitive sequences accounted for approximately 1.24% of the total assembled genome. Among them, Class II had the largest number of repeat elements with a length of 98,478 bp, accounting for 0.75% of total assembly. These results indicated that Class II repeat elements, which encompassed a total of 168 genes, were the most abundant among the identified repetitive sequences. Additionally, the genome of C. utilis contained four pseudogenes, characterized by frame shift and/or premature stop codon mutations, occupying a total of 336 bp. In addition to the identified repetitive sequences, a total of 213 noncoding RNAs (ncRNAs) were also annotated in the C. utilis genome, including 3 ribosomal RNAs (rRNAs), 167 transfer RNAs (tRNAs), and 43 other ncRNAs (Table 2). After masking all repeat sequences, 5,732 protein-coding genes (PCGs) were predicted from the C. utilis genome. The prediction of these genes was based on a combination of ab initio and homology-based gene prediction methods. The number of genes predicted by each approach is listed in Supplementary Table S3.

Out of the 5,732 predicted genes, 127 genes did not yield any results for gene functional annotation in the general database. A total of 5,605 gene annotations were categorized by different databases: Protein family database of alignments and hidden Markov models (Pfam; 4,889 genes), Translation of the EMBL nucleotide sequence (TrEMBL; 5,593 genes), Karyotic Orthology Groups (KOG; 3,949 genes), National Center for Biotechnology Information’s non-redundant protein sequences (Nr; 5,593 genes), GO (2,351 genes), The Swiss-Prot section of the Universal Protein Knowledgebase (Swiss-Prot; 5,013 genes), KEGG (3,131 genes), Transporter Classification Database (TCDB; 187 genes), Pathogen Host Interactions Factors database (PHI; 1,706 genes), and Carbohydrate-Active EnZymes Database (CAZymes; 182 genes). Among the 3,949 genes annotated in the KOG functional classification group, 583, 406, 333, 301 and 275 genes were assigned to five categories: R, “general function prediction only”; O, “Posttranslational modification, protein turnover, chaperones”; J, “Translation, ribosomal structure and biogenesis”; U, “Intracellular trafficking, secretion, and vesicular transport”; and T; “Signal transduction mechanism,” respectively.

GO assignments were used to classify the functions of the predicted genes. A total of 2,351 genes from C. utilis were assigned with different GO terms: biological processes, cellular components and molecular functions (Figure 1A). Biological processes included cellular process (1,251 genes), metabolic process (1,167 genes), single-organism process (804 genes), etc. Cellular components included cell part (951 genes), cell (937 genes), organelle (617 genes), etc. Molecular functions included catalytic activity (1,241 genes), binding (889 genes), transporter activity (144 genes), etc. Moreover, KEGG analysis of the predicted genes revealed that 3,131 genes were categorized into the following KEGG categories: Ribosome (126 genes), Biosynthesis of amino acids (121 genes), Carbon metabolism (99 genes), etc. (Figure 1B). These results indicated that a significant proportion of genes in C. utilis were involved in translation and metabolism.

Comparing the transcriptomes of genetically engineered C. utilis strain (RCT) with the wild-type strain (WT) using RNA-seq technology can provide valuable insights into the genome effect of the δ-zein gene in RCT. The clean reads from each sample were subjected to sequencing and then aligned with the genome obtained in this study, and the alignment efficiency ranged from 93.26 to 94.13%. Among the aligned reads, 89.32 and 89.60% were found to originate from mature mRNA and were mapped to exon regions, respectively (Supplementary Table S4). Alternative splicing of mRNA derived from precursor mRNA could result in the translation of different proteins. The alternative splicing events identified in this study were categorized into 12 types (Figure 2A). To evaluate the adequacy of sequencing data and ensure that it is sufficient for follow-up analysis, the saturation of the number of genes was detected (Figures 2B,C). All samples demonstrated high sequencing quality and correlated well with biological replicates (Supplementary Figure S1). In total, 252 DEGs were identified, including 138 up-regulated genes and 114 down-regulated genes (Table 3). To assess the suitability of the strain as a host for recombinant protein production, we leveraged these transcriptome data, in conjunction with targeted assays. To further depict the overall distribution of gene expression and the fold change in expression levels between the two samples, the MA plot of DEGs was generated (Figure 2D).

The DEGs identified in the two comparison groups. Were subjected to GO enrichment analysis to elucidate their functional annotations in terms of biological processes, molecular functions and cellular components. In cellular components (Figure 3A), DEGs were associated with “integral component of membrane” (GO: 0016021, 57 genes), “ribosome” (GO: 0005840, 14 genes), “mitochondrial large ribosomal subunit” (GO: 0005762, 13 genes), “plasma membrane” (GO: 0005886, 9 genes), and “mitochondrial inner membrane” (GO: 0005743, 8 genes). In molecular functions (Figure 3B), DEGs were mainly involved in “structural constituent of ribosome” (GO: 0003735, 37 genes), “transmembrane transporter activity” (GO: 0022857, 15 genes), “pyridoxal phosphate binding” (GO: 0030170, 7 genes), “rRNA binding” (GO: 0019843, 4 genes) and “NADP binding” (GO: 0050661, 4 genes). In biological processes (Figure 3C), the most highly enriched annotations were “translation” (GO: 0006412, 29 genes), “mitochondrial translation” (GO: 0032543, 9 genes), “oxidation–reduction process” (GO: 0055114, 7 genes), “transmembrane transport” (GO: 0055085, 7 genes), “lipid metabolic process” (GO: 0006629, 5 genes), “sulfate assimilation (GO: 0000103, 4 genes),” “amino acid transport” (GO: 0006865, 4 genes), and “carbohydrate catabolic process” (GO: 0016052, 3 genes).

In addition, the functions of DEGs were further analyzed by KEGG pathway enrichment. The enriched pathways of DEGs between the WT and RCT strains were mainly related to “Ribosome” (ko03010; 25 regulator genes). Moreover, as shown in Figure 4, the distribution pathways of DEGs included “Peroxisome” (ko04146; 11 regulator genes), “Arginine and proline metabolism” (ko00330; 9 regulator genes), “Carbon metabolism” (ko01200; 9 regulator genes) and “Fatty acid degradation” (ko00071; 8 down-regulated genes). The three main KEGG pathways of the up-regulated genes were “Sulfur metabolism” (ko00920), “Aminoacyl-tRNA biosynthesis” (ko00970), and “Cysteine and methionine metabolism” (ko00270), whereas the largest KEGG group of the down-regulated genes was “Fatty acid degradation” (ko00071). Based on these results, we speculate that these pathways may be associated with the effects of δ-zein gene on protein synthesis and lipid synthesis.

To explore the effect of δ-zein gene on C. utilis metabolism, gene set enrichment analysis (GSEA) was performed on the sets of DEGs. GSEA was used to examine several pathways that exhibited significant enrichment among the identified DEGs (Choi et al., 2022). In the analysis, the gene sets of interest were derived from GO terms and KEGG pathways. These gene sets were used to assess the enrichment of up-or down-regulated genes within these specific pathways. The up-regulated gene sets were significantly enriched in various pathways, including “Ribosome” (ko03010; 101 genes), “Biosynthesis of amino acids” (ko01230; 51 genes), and “Ribosome biogenesis in eukaryotes” (ko03008; 48 genes), while down-regulated gene sets were enriched in “MAPK signaling pathway-yeast” (ko04011; 55 genes).

The RNA-seq results were experimentally verified by quantitative Real-Time PCR (qRT-PCR) assays. The expression levels of 4 DEGs (MET17, CHA1, SAMS1, and ARO8) were measured by qRT-PCR, and the results demonstrated a concordant expression of these genes between RNA-seq data and qRT-PCR analysis. The linear correlation analyses between RNA-Seq and qRT-PCR yielded the following correlation coefficients (R²) values: 0.86 (Supplementary Figure S2). As shown in Supplementary Figure S3, the expression levels of MET17, CHA1, SAMS1, and ARO8 genes were observed to be higher in the RCT strain compared to the WT strain. Notably, the expression of MET17 was approximately 1.9 times higher in RCT than in WT, CHA1 showed a 6.2-fold increase, SAMS1 exhibited a 1.6-fold increase, and ARO8 displayed a 3.4-fold increase. The high similarity observed between the qRT-PCR and RNA-seq data of DEGs provides strong evidence for the reliability and accuracy of the transcriptome data analysis. The concordance between these two independent experimental techniques indicates that the conclusions drawn from the transcriptome data are biologically reliable.

WGCNA was performed to gain a deeper understanding of the relationship among the key module genes (Nangraj et al., 2020). To further analyze the features of the module genes, GO and KEGG analyses were conducted. The most significant GO terms for biological process, cellular component, and molecular function (Figure 5A) are demonstrated in the network plot. Genes in the turquoise module were mainly involved in “metabolic process,” “single-organism process,” “cellular process,” “membrane,” “membrane part,” “cell,” “catalytic activity,” “binding,” and “transporter activity.” The KEGG analysis revealed that genes in the turquoise module were mainly associated with “ribosome,” “biosynthesis of amino acids,” “MAPK signaling pathway,” “meiosis,” and “carbon metabolism” (Figure 5B).

To reveal the differential metabolic and genetic information processing between the WT and RCT strains, phenotype microarrays analysis was conducted on the two groups. The metabolic capacity of WT and RCT in utilizing carbohydrate/carbon sources was assessed and validated using a combination of PM1-PM2, offering a novel approach to comprehensively analyze the system-wide response of genetically modified cells. This innovative technique provides valuable insights, complementing data obtained from molecular methods such as RNA-seq analysis. Multiple growth conditions, encompassing a diverse range of 190 carbon sources, were tested by inoculating PM1-PM2 plates with WT and RCT strains. The PM system, a 96-well plate-based assay, was utilized to monitor the variances in cell respiration activity between WT and RCT. Each well of the plate contained a distinct medium and an equivalent amount of tetrazolium dye that develops a purple color in the reduced form. Throughout the growth and maintenance of the strains, the utilization of carbon sources induced distinct phenotypic switching in both WT and RCT. The data collected for each strain line in the OmniLog System were processed by the PM suite of software. This software can automatically compare two strain lines over hundreds of phenotypes. The output from Omnilog was color-coded, facilitating the visual comparison of data between the two strains (Li C. et al., 2022; Li Z. et al., 2022). The metabolic profile curves of WT and RCT could be overlaid, allowing for a clear visualization of the differences between them. The phenotypic responses exhibited by the WT strain were represented by the color red in the strain line, while the RCT strain was represented by the color green. When both strain lines demonstrated similar phenotypic responses, the corresponding kinetic curve appeared yellow. This color-coding scheme allowed for easy identification and differentiation of phenotypic patterns between the two strains (Figure 6; Supplementary Figure S4). The color-coded graphic output facilitated a quick review and visual identification of changes between the two strain lines. By examining the output, significant differences between the strain lines could be easily identified. The software automatically highlighted these differences, allowing for a characterization of the phenotypic variations observed between the two strains. This feature streamlined the analysis process and helped to pinpoint and characterize the specific phenotypic differences detected (Filannino et al., 2016).

The biolog phenotype microarrays revealed positive results for a range of substances. Both the WT and RCT strains demonstrated the ability to metabolize a significant portion of the tested carbon sources, with a success rate of 97.3%. Notably, in plate PM1, both strains were able to metabolize 95 out of 95 carbon sources, while in plate PM2, they metabolized 90 out of 95 carbon sources. However, there were slight differences in carbon metabolism between the WT and RCT strains. In the WT strain, the following substances were highly utilized for carbon metabolism: D-Galactonic Acid-γ-Lactone, D-Tagatose, D-Glucose-6-Phosphate, Pectin, and 5-Keto-D-Gluconic Acid. These particular carbon sources exhibited significant utilization rates in the WT strain, indicating their metabolic significance and potential importance in the cellular processes of the strain (Figure 4A). The observed increases in the utilization of substances such as D-Galactonic Acid-γ-Lactone, D-Tagatose, D-Glucose-6-Phosphate, Pectin, and 5-Keto-D-Gluconic Acid in the WT strain are indicative of enhanced Galactose metabolism and carbohydrate digestion and absorption. Moreover, the higher utilization of D-Arabinose and 2-Deoxy-D-ribose in the WT strain compared to the RCT strain indicated significant differences in pentose and glucuronate interconversions, as well as the pentose phosphate pathway. Additionally, the increase in D-Serine utilization in the WT strain reflected the changes in amino acid metabolism. By simplifying this intricate framework, we were able to identify and describe the phenotypic differences in the WT and RCT strains. In the phenotypic screening using PM4, we specifically focused on 35 different sulfur sources to investigate whether the expression of genes in the RCT strain varied during growth and maintenance, and whether this variability had phenotypic consequences. Our results indicated that the metabolic rates of sulfur sources and phosphorus sources were 100% (35/35 in plate PM4, wells F02-H12) and 89.8% (53/59 in plate PM4, wells A02-E12), respectively (Figure 6; Supplementary Figure S4). The RCT strain exhibited a significant increase in the utilization of O-Phospho-L-Tyrosine and L-Methionine Sulfone, indicating an elevated activity in sulfur metabolism (Figure 4A).