Abstract Background Hepatocellular carcinoma (HCC) is a deadly malignancy known for its ability to evade immune surveillance. NOP2/Sun RNA methyltransferase family member 2 (NSUN2), an RNA methyltransferase involved in carcinogenesis, has been associated with immune evasion and energy metabolism reprogramming. This study aimed to examine the molecular mechanisms underlying the involvement of NSUN2 in immune evasion and metabolic reprogramming of HCC. Methods Single‐cell transcriptomic sequencing was applied to examine cellular composition changes, particularly immune cell dynamics, in HCC and adjacent normal tissues. Bulk RNA‐seq and proteomics identified key genes and proteins. Methylation sequencing and methylated RNA immunoprecipitation (MeRIP) were carried out to characterize the role of NSUN2 in 5‐methylcytosine (m5C) modification of sterol O‐acyltransferase 2 (SOAT2). Clinical samples from 30 HCC patients were analyzed using reverse transcription‐quantitative polymerase chain reaction and Western blotting. Gene expression was manipulated using CRISPR/Cas9 and lentiviral vectors. In vitro co‐culture models and metabolomics were used to study HCC cell‐T cell interactions, energy metabolism, and immune evasion. Tumor growth in an orthotopic mouse model was monitored by bioluminescence imaging, with subsequent measurements of tumor weight, volume, and immunohistochemical staining. Results Single‐cell transcriptomic analysis identified a marked increase in malignant cells in HCC tissues. Cell communication analysis indicated that tumor cells might promote cancer progression by evading immune clearance. Multi‐omics analyses identified NSUN2 as a key regulator in HCC development. MeRIP confirmed that NSUN2 facilitated the m5C modification of SOAT2. Analysis of human HCC tissue samples demonstrated pronounced upregulation of NSUN2 and SOAT2, along with elevated m5C levels in HCC tissues. In vitro experiments uncovered that NSUN2 augmented the reprogramming of energy metabolism and repressed the activity and cytotoxicity of CD8^+ T cells, contributing to immune evasion. In vivo studies further substantiated the role of NSUN2 in fostering immune evasion and tumor formation of HCC by modulating the m5C modification of SOAT2. Conclusions The findings highlight the critical role of NSUN2 in driving HCC progression through the regulation of m5C modification on SOAT2. These findings present potential molecular markers for HCC diagnosis and therapeutic targets for its treatment. Keywords: 5‐methylcytosine modification, energy metabolism reprogramming, Hepatocellular carcinoma, immune escape, NSUN2, SOAT2 __________________________________________________________________ Abbreviations HCC hepatocellular carcinoma cDNA complementary DNA PCA principal component analysis CNV copy number variation GO Gene Ontology KEGG Kyoto Encyclopedia of Genes and Genomes PPI protein‐protein interaction RT room temperature shRNA short hairpin RNA RRA Robust rank aggregation GEO Gene Expression Omnibus DEG differentially expressed genes EtOH ethanol OCR oxygen consumption rate IHC Immunohistochemistry ELISA Enzyme‐linked immunosorbent assay LASSO Least Absolute Shrinkage and Selection Operator SVM‐RFE Support Vector Machine‐Recursive Feature Elimination WB Western blotting RT‐qPCR reverse transcription quantitative polymerase chain reaction PBS phosphate‐buffered saline TEAB triethylammonium bicarbonate MS2 Mass Spectrometry Level 2 AGC Automatic Gain Control TME tumor microenvironment FBS fetal bovine serum UMAP Uniform manifold approximation and projection GAPDH glyceraldehyde 3‐phosphate dehydrogenase m5C 5‐methylcytosine MeRIP methylated RNA immunoprecipitation NSUN2 NOP2/Sun RNA methyltransferase family member 2 SOAT2 sterol O‐acyltransferase 2 scRNA‐seq single‐cell RNA sequencing FC fold change LC‐MS/MS liquid chromatography‐tandem mass spectrometry 1. BACKGROUND Hepatocellular carcinoma (HCC) is one of the most common malignancies worldwide, characterized by high incidence and mortality rates [[36]1, [37]2, [38]3]. Despite advancements in HCC treatment and management, prominent challenges remain, particularly in elucidating the mechanisms of tumor immune escape [[39]4, [40]5]. Tumor immune escape, a critical process enabling tumor cells to escape immune detection and subsequent clearance, contributes markedly to tumor progression and metastasis [[41]4, [42]6, [43]7]. Understanding these mechanisms is critical for the development of novel therapeutic approaches tailored to counteract immune escape in HCC. Recent advancements in high‐throughput sequencing (HTS), particularly single‐cell transcriptomics, have provided valuable insights into the tumor microenvironment (TME) of HCC [[44]8, [45]9, [46]10, [47]11]. Detailed analysis of single‐cell transcriptomics data from HCC and adjacent normal tissues has demonstrated prominent alterations in cellular composition, particularly in immune cell populations, which assume pivotal roles in immune escape in the context of HCC [[48]12]. Furthermore, the integration of HTS with proteomics and meta‐analysis can pinpoint key genes and proteins involved in HCC progression and offer promising targets for therapies against immune escape [[49]8, [50]13, [51]14]. Of note, NOP2/Sun RNA methyltransferase family member 2 (NSUN2), an RNA methyltransferase, has been increasingly recognized for its crucial involvement in the oncogenic phenotypes of various cancers, including HCC [[52]15, [53]16, [54]17]. Given its pronounced role in cancer biology, the function of NSUN2 in HCC, particularly concerning immune escape, has gained growing attention. NSUN2‐mediated 5‐methylcytosine (m5C) modification of mRNA has been proposed to enhance the methylation of sterol o‐acyltransferase 2 (SOAT2) to trigger metabolic reprogramming and allow tumor cells to evade immune surveillance [[55]18]. Furthermore, the aberrant m5C modification of long non‐coding RNA H19 by NSUN2 has been associated with poor differentiation in HCC, highlighting its pivotal role in tumor development [[56]19]. Based on existing evidence on the significance of NSUN2 in tumor progression, this study was undertaken to clarify its involvement in immune escape in HCC, with a specific focus on its impacts on the m5C modification of SOAT2 RNA. Multi‐omics analyses, incorporating single‐cell transcriptomics and proteomics, and targeted In vivo and In vitro experimental investigations were utilized to elucidate the mechanistic underpinnings of NSUN2‐mediated m5C modification of SOAT2 in HCC. The study sought to uncover the molecular mechanisms contributing to immune escape and to pinpoint potential therapeutic targets. This study integrated multi‐omics analyses with experimental investigations to provide insights into the significance of NSUN2 in HCC. 2. MATERIALS AND METHODS 2.1. Construction of the mouse HCC model A total of 44 male C57BL/6J mice (aged 12 weeks, weighing 21‐25 g) were purchased from Hunan Slac Jingda Laboratory Animal Co., Ltd. (Changsha, Hunan, China). The mice were individually housed in the specific pathogen‐free animal facility, where conditions were controlled at 60%‐65% humidity and a temperature of 25 ± 2°C. Ad libitum access to food and water was provided, and the mice were maintained on a 12‐h light‐dark cycle. Following a one‐week acclimatization, the health of each mouse was assessed before initiating the experiment. All experimental procedures and protocols for animal use were reviewed and approved by the Animal Ethics Committee of Renji Hospital, School of Medicine, Shanghai Jiaotong University (Shanghai, China). In this study, an orthotopic HCC model was developed in C57BL/6J mice with underlying liver cirrhosis. The induction of liver cirrhosis was achieved through intraperitoneal injections of carbon tetrachloride (CCl[4], 0.5 mL/kg, 488488, Sigma‐Aldrich, St. Louis, MO, USA) twice a week for 6 weeks. Subsequently, Hepa1‐6 cells (murine HCC cell line, CRL‐1830, American Type Culture Collection [ATCC], Manassas, VA, USA; 1 × 10^8 cells/mL) were thoroughly mixed with Matrigel (354262, Corning, New York, NY, USA) at a 4:1 ratio and subsequently injected into the left hepatic lobe. The Matrigel mixture was maintained on ice until injection and gently mixed prior to use. Although Matrigel provides a supportive scaffold conducive to tumor growth, it is essential to acknowledge that Matrigel does not fully capture the complexity of the TME, particularly in terms of extracellular matrix components, stromal cells, immune cells, and signaling molecules. This limitation may affect the representativeness of the In vivo conditions observed. The mice were anesthetized using inhaled isoflurane (2% in 3 L/min oxygen carrier gas), and the abdominal area was prepared by shaving and sterilization. A longitudinal incision was made along the midline using ophthalmic scissors to expose the xiphoid, with the incision held open by curved metal hooks threaded with sutures. A small pad (2 cm in diameter) was placed on the mouse's back to expose the liver, and gentle pressure was applied to the ribs on both sides to extrude a portion of the left hepatic lobe. The focal distance of the stereomicroscope was adjusted to clearly expose the left hepatic lobe. The cell and matrix mixture was gently agitated to ensure even distribution. A 10 µL microinjection needle was used to slowly aspirate 5 µL of the cell suspension. The needle was inserted into the liver surface at a 20°C angle, leveled parallel to the table, and the suspension was injected gradually. The liver tissue at the injection site displayed a rapid color change, becoming lighter or forming cloudy translucent bubbles. The needle remained in position for 1 min post‐injection, then was slowly withdrawn from the affected area, pausing for 2‐4 min before complete removal. The needle entry point was carefully cauterized with a needle‐type electric soldering iron to prevent tumor cell leakage. The liver surface was then wiped with a cotton swab soaked in 75% ethanol (EtOH), and the left hepatic lobe was repositioned into the abdominal cavity. The peritoneum and skin were sutured, and the suture area was disinfected with 75% EtOH. The mice were positioned on a heating pad set to 35°C until full recovery, after which they were returned to their cages for regular care. Tumor growth was assessed 4 days after injection using In vivo live imaging, high‐frequency ultrasound imaging and bioluminescence imaging, with tumor volume measurements recorded every 7 days [[57]20]. Tumor volume was recorded at weeks 1, 2, and 4 post‐injection. The average tumor size reached approximately 10 mol/L in diameter by the end of the study. The In vivo live imaging was performed using the IVIS Spectrum Imaging System (PerkinElmer, Waltham, MA, USA), equipped with a cooled CCD camera to detect bioluminescent and fluorescent signals. 2.2. Histological analysis Upon completion of the experiment, the mice were euthanized, and the entire liver was harvested for histological analysis. The liver tissues were immobilized in 10% formalin (Formal‐Fixx, [58]R21923, Thermo Fisher Scientific Inc., Madison, WI, USA) and subsequently processed for histopathological evaluation. Hematoxylin and eosin (H&E, H3136, Sigma‐Aldrich) staining was carried out to visualize the tissue architecture. The tumor regions exhibited characteristic features of HCC, while the adjacent regions appeared normal, indicating that the tumor was localized and did not spread to adjacent normal tissues. 2.3. Single‐cell RNA sequencing (scRNA‐seq) of liver tissues from HCC mouse model Two randomly selected mice with successful modeling and 12‐week feeding were used in this study. Their tumor tissues and adjacent normal tissues were disrupted and digested with 1 mg/mL collagenase (C2674, Sigma‐Aldrich) at 37°C for 10 min. This process was followed by a brief incubation with ethylenediaminetetraacetic acid (EDTA) (25200072, Gibco, Waltham, MA, USA) at 37°C for an additional 5 min to create single‐cell suspensions. Individual cells were captured using the C1 Single‐Cell Auto Prep System (Fluidigm, Inc., South San Francisco, CA, USA). mRNA was released from the lysed cells on the chip and reverse‐transcribed into complementary DNA (cDNA). This cDNA was pre‐amplified on a microfluidic chip, preparing it for subsequent sequencing steps. Library construction was then completed using the amplified cDNA, and scRNA‐seq was performed on the HiSeq 4000 Illumina platform (Illumina, San Diego, CA, USA) with paired‐end reads, a read length of 2 × 75 bp, and estimated 20,000 reads per cell. The scRNA‐seq data were analyzed to substantiate the overexpression of NSUN2 in tumor cells, rather than in other immune cell types. The “Seurat” package (Seurat v4.3.0.1, Satija Lab, New York, NY, USA; available at [59]https://satijalab.org/seurat/) was utilized to segregate the data into distinct cell clusters and pinpoint cell types based on known lineage‐specific marker genes. Distinct clusters corresponding to epithelial cells, immune cells, stromal cells, and other cell populations were identified. Quality control of the data was ensured by applying thresholds where the number of detected genes (nFeature_RNA) ranged between 200 and 5,000, and the mitochondrial gene content (percent.mt) was below 25%. For subsequent analyses, the top 2,000 highly variable genes were selected based on variance [[60]21]. Following cell type annotation, the “FeaturePlot” function in “Seurat” was applied to visualize NSUN2 expression across different cell types, offering a direct view of its distribution, particularly in tumor cells versus immune cells. The potential impact of NSUN2 on the activity of endothelial progenitor cells (EPCs) in the TME was examined by reanalyzing scRNA‐seq data from the Hepa1‐6 HCC model. The levels of EPC‐related marker genes, specifically CD34 and VEGFR2, were determined to assess changes in response to NSUN2 overexpression or knockout (KO). Clustering analysis of the scRNA‐seq data was executed using the “Seurat” package, and EPC‐related genes significantly expressed in different cell populations were identified through the “FindMarkers” function. The “CellChat” package was also utilized to explore the signaling pathway activities of these genes between different cell types, providing further insights into the role of EPCs in the NSUN2‐mediated TME. 2.4. Uniform manifold approximation and projection (UMAP) clustering and cell annotation Dimensionality reduction of the scRNA‐Seq dataset was achieved through principal component analysis (PCA) on the top 2,000 highly variable genes. The initial 30 principal components (PCs) were selected for further downstream investigation, as determined by the “Elbowplot” function of the “Seurat” package. The identification of major cell types was accomplished using the “FindClusters” function in “Seurat”. The UMAP algorithm was subsequently applied for non‐linear dimensionality reduction of the scRNA‐seq data. Cell types were annotated by comparing them with known lineage‐specific marker genes and utilizing the “CellMarker” database available online ([61]http://biobigdata.hrbmu.edu.cn/CellMarker/). The “CellCycleScoring” function was employed to assess the cell cycle phases of the sampled cells. Detailed marker genes corresponding to the annotations of each cell type are provided in Supplementary Table [62]S1. Batch correction of the sample data was executed using the “Harmony” package (v0.1.1, [63]https://github.com/immunogenomics/harmony), followed by the application of the “ElbowPlot” function to prioritize the standard deviations of the PCs. The R package “CellChat” (v1.6.1, Jin Lab, [64]https://github.com/sqjin/CellChat) was employed to explore pathway activity across different cell populations and further assess the functional differences between cell types in tumor versus normal groups. 2.5. Malignant cell extraction Malignant cells were identified by first isolating epithelial cells with the “subset” function, followed by copy number variation (CNV) analysis using the “inferCNV” package (v1.14.1, [65]https://github.com/broadinstitute/inferCNV). The “K‐means” algorithm was implemented using the “stats” package in R (v4.2.1, R Core Team, Vienna, Austria) to further exclude non‐malignant cells. Temporal analysis of these cells was performed using the “Monocle2” package (v2.26.0, Trapnell Lab, [66]https://cole‐trapnell‐lab.github.io/monocle‐release/). During the cell sorting process, data dimensionality reduction was performed using the “DDRTree” method with default parameter settings: the maximum dimensionality was set to 2 (max_components = 2), and the resolution for pseudotime state inference was set to high resolution (reduction_method = DDRTree). In addition, highly variable genes were selected as ordering genes (ordering_genes), which were screened based on significant differences in gene expression (q‐value < 0.05). Trajectory construction was based on the expression patterns of the ordering genes, and the evolutionary pathways of the cells were inferred through pseudotime analysis. During the inference process, a log‐transformed expression matrix was used, and branch points and developmental directions were determined based on clustered cells. [[67]22, [68]23]. 2.6. RNA‐seq and analysis of liver tissues from HCC mouse model Three mice, successfully modeled and maintained for 12 weeks, were randomly selected following the construction of the HCC mouse model. Criteria for humane endpoints included a weight reduction surpassing 20% of body weight, severe lethargy, an inability to access nourishment or hydration, or indications of severe, unmanageable pain or distress. Euthanasia was performed by cervical dislocation under deep anesthesia induced by 2% isoflurane (R510‐22‐10, Ruivode, Shenzhen, Guangdong, China). Total RNA was isolated from HCC and adjacent normal tissues of the mice using the Total RNA Isolation Kit (AM1912, Thermo Fisher Scientific). The optical density (OD) of the RNA was quantified at 450 nm using the BioSpectrometer basic UV‐vis spectrophotometer (Eppendorf, Hamburg, Germany). RNA integrity was verified through agarose gel electrophoresis. RNA with an OD260/OD280 ratio ranging from 1.8 to 2.0 was considered high‐quality RNA. The high‐quality RNA underwent reverse transcription to generate cDNA, which was then utilized for constructing the RNA library. Sequencing was carried out on the Illumina NextSeq 500 platform, and the raw image data obtained from sequencing were transformed into raw reads through the base‐calling process. Quality control was performed using the “Cutadapt” software (v3.5, [69]https://github.com/marcelm/cutadapt), which removed adapter sequences and filtered out low‐quality reads. Specifically, reads containing adapter sequences, shorter than 30 bp, and with bases having a Q‐value less than 20 were considered low‐quality and removed. The remaining high‐quality reads after filtering were designated as “clean reads”. The processed clean reads were subsequently aligned to the mouse reference genome (mm10, GRCm38) utilizing the “Hisat2” package (v2.2.1, [70]https://daehwankimlab.github.io/hisat2/) [[71]24]. The identification of differentially expressed genes (DEGs) from the bulk RNA‐seq data was performed using the “limma” package (v3.54.2) based on the screening criteria set at |log[2] fold change (FC)| > 6 and P < 0.001. Visualization of the DEGs was achieved by generating volcano plots with the “ggplot2” package (v3.4.2) and heatmaps with the “pheatmap” package (v1.0.12). Enrichment analysis of the target gene set was performed using the Disease Ontology (DO) database ([72]http://www.disease‐ontology.org/). DO pathway enrichment analysis was conducted with the R packages “enrichplot (v4.6.2),” “org.Hs.eg.db (v3.16.0),” and “clusterProfiler (v1.18.3)” packages. Visualization of the results was generated using the “ggplot2” package. The SangerBox database (available at: [73]http://sangerbox.com/home.html) was utilized for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the identified DEGs. The Gene Expression Profiling Interactive Analysis (GEPIA) database ([74]http://gepia.cancer‐pku.cn/detail.php) was used to analyze the gene expression patterns in The Cancer Genome Atlas‐Liver Hepatocellular Carcinoma (TCGA‐LIHC). Furthermore, the Xiantaozi academic database ([75]https://www.xiantaozi.com/) was employed to examine the impact of individual gene expression on immune cell infiltration within TCGA‐LIHC. 2.7. Preparation and determination of proteomics samples Three HCC‐bearing mice were sacrificed after 12 weeks, and HCC and adjacent normal tissues were collected. These samples were ground into fine powder using liquid nitrogen in a mortar. The resulting powder was then placed in a 5 cm^3 centrifuge tube and subjected to ultrasonic disruption with an ultrasonic cell disruptor (SCIENTZ‐IID, Scientz, Ningbo, Zhejiang, China) while kept in an ice bath. The treatment medium comprised phenol (100206, Sigma‐Aldrich) extraction buffer containing 10 mmol/L dithiothreitol (DTT) (R0861, Solarbio, Beijing, China), 1% protease inhibitor mixture (P6731, Solarbio), and 2 mmol/L EDTA (E1170, Solarbio). An equivalent volume of Tris‐saturated phenol (HC1380, Solarbio) at pH 8.0 was added to the mixture, which was then vortexed for 4 min. The phenol phase, located at the top, was carefully transferred to a fresh centrifuge tube. Following this, an equal volume of saturated ammonium sulfate (59845, Sigma‐Aldrich) in methanol (34860, Sigma‐Aldrich) at a 1:5 ratio relative to the phenol solution was added. The solution was left undisturbed overnight to promote protein precipitation. Next, centrifugation was conducted at 4°C and 5,000 ×g for 10 min, and the remaining pellet was rinsed once with cold methanol and three times with cold acetone, all under 4°C conditions. The purified protein was re‐dissolved in 8 mol/L urea (U8020, Solarbio), and the protein concentration was quantified using a BCA assay kit (P0012, Beyotime, Shanghai, China) [[76]25]. 2.8. Enzymatic digestion, peptide labeling gradient, and nano‐liquid chromatography‐tandem mass spectrometry (nano‐LC‐MS/MS) analysis Each sample underwent enzymatic digestion with 50 µg of protein. The protein solution was mixed with 5 mol/L DTT (DTT‐RO, Sigma‐Aldrich) for a 30‐min incubation at 56°C. Subsequently, acetamide was introduced to achieve a final concentration of 11 mmol/L, followed by a 15‐min incubation at room temperature (RT). Urea concentration in the solution was then diminished to below 2 mol/L, and trypsin (25200056, Thermo Fisher Scientific) was added at a 1:50 (w/w) ratio of trypsin to protein for overnight digestion at 37°C. The digestion was extended for additional 4 h after adding more trypsin at a 1:100 (trypsin:protein) ratio. Following enzymatic digestion, peptides were purified and desalted using a HyperSep™ C18 column (60108‐302, Thermo Fisher Scientific), after which they were vacuum‐dried. The dried peptides were then reconstituted in 0.5 mol/L triethylammonium bicarbonate (TEAB) (90114, Thermo Fisher Scientific) and labeled through the Tandem Mass Tag (TMT) Reagent Kit (90064CH, Thermo Fisher Scientific). The TMT reagent was thawed and reconstituted in acetonitrile (113212, Sigma‐Aldrich), followed by a 2‐h incubation at RT. Post‐labeling, peptides were desalted, vacuum centrifuged (700 ×g, 10 min) at RT to remove cell debris and nuclei, followed by additional centrifugation (3,000 ×g, 25 min) to isolate the peptides. Labeled peptides were pooled, reconstituted with the Pierce™ High pH Reversed‐Phase Peptide Separation Kit (84868, Thermo Fisher Scientific), and divided into 15 fractions. Each fraction was dried and reconstituted in 0.1% formic acid (159002, Sigma‐Aldrich). For nano‐LC separation, 2 µg of peptides from each sample were loaded onto the Easy nLC 1200 nano‐UPLC system (Thermo Fisher Scientific). The separated peptides were then analyzed using a Q‐Exactive HCX mass spectrometer (Thermo Fisher Scientific) over 60 min, with an electrospray voltage of 2.1 kV in positive ion mode. The mass scan range for precursor ions was set at 350‐1,200 m/z, with a primary mass resolution of 60,000 at m/z 200. Higher‐energy collisional dissociation (HCD) served as the MS2 Activation Type, with an isolation window of 20 Th (Thomson) and a normalized collision energy of 32 eV [[77]26]. 2.9. Protein‐protein interaction (PPI) network analysis Quantitative proteomics analysis was employed to pinpoint differentially expressed proteins (DEPs) in HCC tissue samples. Protein levels between HCC tissues and adjacent normal tissues were compared using MS and bioinformatics tools. The identified DEPs were used to construct PPI networks. The DEPs were imported into the STRING database (available at: [78]https://cn.string‐db.org/) for network analysis, and the PPI network was visualized using Cytoscape v3.10.0 software [[79]26]. 2.10. Clinical sample collection The study included 30 patients with HCC who underwent surgical treatment at Renji Hospital, School of Medicine, Shanghai Jiaotong University from January 2022 to January 2023. HCC tissues were collected from these patients, and adjacent normal tissues located 2 cm away from the cancerous lesion were selected as controls. The inclusion criteria specified that patients were aged between 18 and 80 years, had undergone surgical treatment without prior anti‐cancer therapies such as radiation or chemotherapy, and had no history of other diseases. The tissues were flash‐frozen immediately using liquid nitrogen and then stored at ‐80°C. Informed consent was obtained from all patients, and the study was approved by the Ethics Committee of Renji Hospital, School of Medicine, Shanghai Jiaotong University, in compliance with the Helsinki Declaration. 2.11. Isolation of primary cells of mice Three healthy C57BL/6J male mice (aged 6‐8 weeks, weighing 21‐24 g) were obtained from Hunan Slac Jingda Laboratory Animal Co., Ltd. The mice were euthanized by cervical dislocation under deep anesthesia induced by 2% isoflurane. CD8^+ T cells were isolated from the spleen and lymph nodes of the mice using the Mouse CD8^+ T Cell Isolation Kit (8804‐6822‐74, Invitrogen, Carlsbad, CA, USA). Following a 10‐min centrifugation (3,000 ×g) at 4°C, the supernatant was retrieved and passed through a 0.2 µm filter [[80]27]. Prior to experimentation, isolated CD8^+ T cells (2 × 10^6 cells) were stimulated for 24 h with 5 µg/mL mouse CD3/CD28 activator (11452D, Gibco). The cells were maintained in a humidified Heracell™ Vios 160i CO[2] Incubator (Thermo Fisher Scientific) under conditions of 37°C and 5% CO[2]. The culture medium consisted of RPMI‐1640 (11875093, Gibco) enriched with 10% fetal bovine serum (FBS) (A5670701, Gibco) and 1% penicillin/streptomycin (15140148, Gibco) [[81]28, [82]29, [83]30, [84]31]. 2.12. Construction and culture of In vitro model The mouse HCC cell line Hepa1‐6 was cultured in Dulbecco's modified eagle medium (DMEM) (11965092, Gibco) containing 10% FBS, 10 µg/mL streptomycin, and 100 U/mL penicillin. The HCC cell line H22 (CL‐0341, Procell, Wuhan, Hubei, China) was maintained in RPMI‐1640 medium enriched with 10% FBS and 1% penicillin/streptomycin. Before the experiments, the identity and purity of both cell lines were confirmed through short tandem repeat (STR) genotyping to ensure proper cell line identification and avoid any cross‐contamination. Additional control experiments were executed to examine the independent effect of NSUN2 overexpression on CD8^+ T cell activity. The number of CD8^+ T cells was kept constant, while the number of tumor cells was varied in the co‐culture setup. The co‐culture of HCC cells and CD8^+ T cells at a 1:5 ratio for 48 h was performed by incubating the HCC cell lines with pre‐processed CD8^+ T cells in a co‐culture system at 37°C with 5% CO[2]. The HCC cell lines were placed in the upper insert of Transwell plates (CLS3412, Corning, New York, NY, USA), while CD8^+ T cells were inoculated in the lower insert. The two cell types were separated by a 0.4 µm membrane, allowing only soluble factors to pass through. After 3 days, the supernatant and cultured cells were collected for further analysis. The 293T cell line (CRL‐3216) was acquired from ATCC and cultured in DMEM enriched with 10% FBS, 10 µg/mL streptomycin, and 100 U/mL penicillin. All cell lines were maintained in a humidified incubator at 37°C with 5% CO[2], and subculturing was carried out upon reaching 80%–90% confluency [[85]32]. The cell co‐culture groups were established as follows: (1) control group: co‐culture with varying numbers of tumor cells without modification of NSUN2 expression; (2) NSUN2‐wild‐type (WT) group: co‐culture of Hepa1‐6/H22 cells expressing NSUN2‐WT and CD8^+ T cells; (3) NSUN2‐KO group: co‐culture of NSUN2‐KO Hepa1‐6/H22 cells and CD8^+ T cells; (4) negative control (NC)‐overexpression (oe) group: co‐culture of Hepa1‐6/H22 cells transduced with an NC‐oe vector and CD8^+ T cells; (5) NSUN2‐WT + NC‐oe group: co‐culture of Hepa1‐6/H22 cells expressing NSUN2‐WT, transduced with an NC‐oe vector, along with CD8^+ T cells; (6) NSUN2‐KO + NC‐oe group: co‐culture of NSUN2‐KO Hepa1‐6/H22 cells transduced with an NC‐oe vector, along with CD8^+ T cells; (7) NSUN2‐KO + SOAT2‐oe group: co‐culture of NSUN2‐KO Hepa1‐6/H22 cells transduced with an SOAT2‐oe vector, along with CD8^+ T cells. 2.13. Construction of NSUN2‐KO cells using CRISPR/Cas9 NSUN2‐KO cells were generated using the CRISPR/Cas9 technique. The single guide RNA (sgRNA) sequences were as follows: NSUN2‐sgRNA1: 5'‐TGGCTATCCCGAGATCGTAA‐3' (PAM: AGG); NSUN2‐sgRNA2: 5'‐CTTGAGGAAGTCCCCGTTGT‐3' (PAM: TGG). The sgRNA sequences were integrated into the Lenti‐CRISPR v2 plasmid, designed to express the Streptococcus pyogenes Cas9 nuclease gene (Hanbio Biotechnology, Shanghai, China). The CRISPR/Cas9 system, delivered via the Lenti‐CRISPR v2 vector, was employed to transduce cells and facilitate the creation of NSUN2‐KO cells. The strategy ensured that the genomic region between sgRNA1 and sgRNA2 was completely excised, leading to a full knockout of NSUN2 at the DNA level and preventing its mRNA expression. Selection of transduced cells was performed using 4 µg/mL puromycin (A1113803, Gibco) to screen for sgRNA plasmid and donor sequence. Surviving clones, isolated through limiting dilution, were screened via reverse transcription‐quantitative polymerase chain reaction (RT‐qPCR) and Western blotting (WB) to confirm NSUN2‐KO. DNA sequencing was subsequently conducted to validate the NSUN2‐KO cells [[86]33, [87]34]. 2.14. Development of lentiviral vectors for silencing and overexpression Potential short hairpin RNA (shRNA) target sequences of SOAT2 were identified based on the mouse cDNA sequence from NCBI GeneBank (available at: [88]https://www.ncbi.nlm.nih.gov/genbank/). Three specific target sequences against SOAT2 were designed, along with a non‐targeting sequence used as a sh‐NC. The shRNA target sequences used are listed in Supplementary Table [89]S2. Oligonucleotides were synthesized by GenePharma^® (Shanghai, China). The lentiviral packaging system was established using the pLKO.1 vector, which serves as a lentiviral gene silencing vector. Co‐transfection of 293T cells with the lentiviral constructs was performed using Lipofectamine 2000 (11668500, Invitrogen) upon reaching 80%–90% confluence. After 48 h of culture, the supernatant was collected, filtered, and centrifuged at 1,500 ×g for 10 min at 4°C to harvest the viral particles. The viruses were harvested from exponentially growing cells, and viral titration was determined. Overexpression lentiviral vectors were developed and packaged by Genechem (Shanghai, China) using the lentiviral vector‐platelet‐derived growth factor receptor alpha (LV‐PDGFRα) [[90]35, [91]36, [92]37, [93]38]. 2.15. Cell transfection Prior to the construction of the In vitro model, Hepa1‐6/H22 cells at the exponential phase were detached with trypsin to create a cell suspension (5 × 10^4 cells/mL). A volume of 2 mL of this suspension was seeded into each well of a 6‐well plate. Lentivirus, at a multiplicity of infection of 10 and a viral titer of 1 × 10^8 TU/mL, was introduced into the culture medium, followed by a 48‐h incubation. Stable cell lines were established by treating the cultures with 2 µg/mL puromycin (Sigma‐Aldrich) and maintained in culture for 2 weeks [[94]39]. The transfection of Hepa1‐6/H22 cells was categorized into the following groups: (1) sh‐NC group: transduced with lentivirus harboring NC vector; (2) shSOAT2 #1 group: transduced with lentivirus harboring shSOAT2 vector #1; (3) shSOAT2 #2 group: transduced with lentivirus harboring shSOAT2 vector #2; (4) shSOAT2 #3 group: transduced with lentivirus harboring shSOAT2 vector #3; (5) NC‐oe group: transduced with lentivirus harboring NC‐oe vector; (6) NSUN2‐oe group: transduced with lentivirus harboring NSUN2‐oe vector; (7) SOAT2‐oe group: transduced with lentivirus harboring SOAT2‐oe vector. The design and synthesis of all the mentioned lentivirus vectors were performed by Guangzhou Raybo Biotech (Guangzhou, Guangdong, China). 2.16. RNA stability analysis The stability of RNA of HCC cells was assessed by treating the cells with 5 µg/mL actinomycin D ([95]SBR00013, Sigma‐Aldrich). Samples were harvested at 2, 4, 6, and 8 h post‐incubation, and RNA was subsequently extracted for RT‐qPCR analysis [[96]40]. The primers used for RT‐qPCR are as follows: glyceraldehyde 3‐phosphate dehydrogenase (GAPDH; housekeeping gene) forward: 5'‐GGAGAGTGTTTCCTCGTCCC‐3'; GAPDH reverse: 5'‐ATGAAGGGGTCGTTGATGGC‐3'; NSUN2 forward: 5'‐TGTGCGCTCCGCGTG‐3'; NSUN2 reverse: 5'‐CCCAGCCCGCCTGG‐3'; SOAT2 forward: 5'‐TGGAGGTGCAACATTTCCGA‐3; SOAT2 reverse: 5'‐AGCATCAACCTGCCCTCATC‐3'. The length of the PCR products for GAPDH, NSUN2, and SOAT2 was approximately 150‐200 bp. The reaction conditions for RT‐qPCR were set as follows: an initial denaturation at 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. Relative expression was quantified using the 2^−ΔΔCt method. 2.17. Extraction of lipid metabolites from tissues and cells The extraction of metabolites followed a precise method. After reaching confluency, co‐cultured HCC cells were digested with Trypsin‐EDTA, and pre‐treated HCC tissues were freeze‐dried. Samples were then resuspended in a pre‐cooled EtOH/acetonitrile/water solution (2:2:1, v/v; 100 µL), vortexed, and ultrasonicated at 4°C for 30 min. Following a 10‐min incubation at ‐20°C, the samples underwent centrifugation (14,000 ×g) for 20 min at 4°C. Thereafter, the supernatant was vacuum‐dried and subsequently reconstituted in acetone/water (1:1, v/v; 100 µL). After vortexing, the solution was centrifuged under the same conditions again for 15 min. The resulting supernatant was reserved for further analysis. The high‐performance liquid chromatography (HPLC) procedure was conducted under specific conditions to ensure optimal results. A hydrophilic interaction liquid chromatography (HILIC) column was used in conjunction with an ultra‐high‐performance liquid chromatography system (UHPLC‐Q‐Exactive Orbitrap HRMS, Thermo Fisher Scientific). The column was kept at a stable temperature of 40°C, while the flow rate was set at 0.4 mL/min, with a sample injection volume of 2 µL. The Quadrupole Time‐of‐Flight (Q‐TOF) MS conditions were established as follows. The electrospray ionization source was configured with ion source gases I and II, both set at 50, with a curtain gas of 30 and a source temperature of 500°C. The ion spray voltage was configured at ±5,500 V for both ion modes. Precursor ion scanning was conducted in high‐sensitivity mode utilizing information‐dependent acquisition, with collision energy fixed at 35 ± 15 eV and declustering potential set at ±80 V for both negative and positive ions [[97]41, [98]42, [99]43]. 2.18. Cholesterol level measurement The raw MS data in wiff.scan format were first converted to. mzXML files with ProteoWizard MSConvert. Peak calibration and retention time alignment were performed through Cross‐platform Chromatogram Management System (OpenChrom, [100]http://www.openchrom.net). The extracted peak areas were then subjected to isotope and adduct annotation using the Collection of Algorithms of MEtabolite pRofile Annotation (CAMERA) algorithm set. Cholesterol compounds were identified by comparing the m/z (<10‰) and Tandem Mass Spectrometry (MS/MS) spectra against the mzCloud (available at: [101]www.mzcloud.org/), mzVault, and Masslist databases [[102]41, [103]44]. Cholesterol distribution and its possible involvement in intercellular signaling were assessed by measuring cholesterol concentrations within the cells and in the surrounding culture medium. After 48 h of standard cell culture, the supernatant was collected for analysis. Cholesterol levels in the medium were determined using the Amplex™ Red Cholesterol Assay Kit (A12216, Thermo Fisher Scientific), where the samples were combined with assay reagents and subjected to an enzymatic reaction. The absorbance was then recorded at 570 nm, with a standard curve applied for accurate cholesterol quantification. 2.19. Meta‐analysis Nine datasets relevant to HCC were retrieved from nine human Gene Expression Omnibus (GEO) datasets ([104]GSE25599, [105]GSE105130, [106]GSE112705, [107]GSE135631, [108]GSE164359, [109]GSE202069, [110]GSE202853, [111]GSE214846, and [112]GSE216613) (Supplementary Table [113]S3). In this meta‐analysis, HCC tumor samples from different patient populations were analyzed to capture the variation in NSUN2 expression across different geographic regions and etiological contexts. For Asian patients, the samples were predominantly from HCC cases linked to hepatitis B virus (HBV) or hepatitis C virus (HCV) infections. In contrast, the American patient samples were primarily associated with HCC cases driven by metabolic‐associated steatohepatitis. Furthermore, NSUN2 expression in these samples was measured using RT‐qPCR, and the expression was compared between HCC and adjacent normal tissues. The “meta” package (v6.5‐0, [114]https://github.com/guido‐s/meta) in R software (The R Foundation for Statistical Computing, Vienna, Austria; [115]https://www.r‐project.org) was utilized for meta‐analysis, with continuous variables as the data type. Due to variations in detection techniques and measurement units among the datasets, the standard mean difference (SMD) and 95% confidence interval (CI) were utilized to estimate the effect size. A heterogeneity test was performed to examine the suitability of combining the studies, using the Q test and I^2 test to determine the presence and extent of heterogeneity based on τ^2 value, P value, and I^2 statistics. A fixed‐effects model was applied in the absence of statistical heterogeneity (P > 0.05, I^2 < 50%), suggesting that differences among studies were attributable to random sampling error. Conversely, when significant heterogeneity was found (P < 0.05, I^2 > 50%), a random‐effects model was adopted for meta‐analysis. Subgroup analysis was employed to further investigate data heterogeneity, while leave‐one‐out sensitivity analysis was conducted to assess the robustness of the results [[116]45]. The sensitivity analysis involved sequentially removing each study to observe any significant changes in the summary effect size. Consistency in the SMD values confirmed the stability of the results, indicating robust meta‐analysis findings. Egger's test was performed to statistically evaluate publication bias, where a P value greater than 0.05 indicated no prominent publication bias. 2.20. Observation of mitochondrial morphology Mitochondria were visualized by transfecting cells with the mitochondrial green fluorescent protein (mito‐GFP) plasmid ([117]C10600, Invitrogen) or by staining with 100 nmol/L Mito‐Tracker Red dye (M7512, Invitrogen) at 37°C. The cells were immobilized with 4% paraformaldehyde (PFA) and permeabilized with 0.1% Triton X‐100 (P0096, Beyotime). Samples underwent incubation in a solution containing 0.3 mol/L glycine (50046, Sigma‐Aldrich) and PBS with 5% normal goat serum (16210072, Gibco) for 1 h. If required, cells were further treated with 500 n mol/L AlexaFluor660‐phalloidin (A22285, Invitrogen) or 100 nmol/L tetramethylrhodamine isothiocyanate (TRITC)‐phalloidin (P1951, Sigma‐Aldrich) and 10 µmol/L 4',6‐diamidino‐2‐phenylindole (DAPI) (C1002, Beyotime) for 30 min at 22°C. Mitochondrial length was determined by creating z‐axis maximum intensity projections from the red channel (Mitotracker or mito‐dsRed) at 0.2 µm intervals. Cell regions with clearly visible mitochondria were identified, and the line tool in Nikon Elements software (v5.30, Nikon Corporation, Tokyo, Japan) was employed to measure 25‐30 mitochondria per cell. Mitochondria shorter than 1.0 µm were classified as fragmented. Cells in which over 80% of the mitochondria displayed fragmentation were deemed to be undergoing mitochondrial fission [[118]46]. 2.21. Seahorse analysis The oxygen consumption rate (OCR) was assessed using the Extracellular Flux System 24 (Agilent, Palo Alto, CA, USA). Cells underwent 24‐h EtOH treatment in a medium containing 5% charcoal‐stripped serum. Following trypsinization, the cells were plated (3 × 10^4 cells/100 µL per well) and subjected to incubation for 4 h to facilitate cell attachment. An unbuffered XF assay medium (103680‐100, Agilent, Palo Alto, CA, USA) was prepared and supplemented with 2 mol/L GlutaMAX (35050079, Gibco), 1 mol/L sodium pyruvate (11360070, Gibco), and 10 mol/L glucose (pH 7.4 at 37°C; 15023021, Gibco). The medium was equilibrated at 37°C and 0.04% CO[2] for 30 min before the experiment. OCR was measured using an extracellular flux analyzer, which was programmed to operate in three cycles: mix (150 s), wait (120 s), and measure (210 s) [[119]46]. 2.22. Glucose consumption assay Glucose levels were quantified with the glucose fluorescence assay kit (CA‐G005, eEnzyme, Gaithersburg, MD, USA). Cell lysates were diluted to a protein concentration of 10 µg, mixed with reaction buffer, and incubated for 10 min at 37°C. Fluorescence was then measured using a microplate reader supplied by Molecular Devices (Sunnyvale, CA, USA) [[120]46]. 2.23. Lactate production measurement Lactate secretion was quantified using an L‐lactate assay kit (ab65331, Abcam, Cambridge, UK). The collected cell culture medium was diluted with assay buffer and subjected to a 20‐min incubation with the reaction mixture at RT. Absorbance was then recorded at 450 nm using the Molecular Devices^® microplate reader to determine lactate levels [[121]46]. 2.24. ATP production evaluation ATP production was quantified using the ATP assay kit (ab113849, Abcam). Cells were lysed in the provided lysis buffer, and the resulting lysate was transferred into a 96‐well plate. Following the addition of the reaction mixture, the plate was subjected to a 20‐min incubation at RT. Luminescence was measured using the Molecular Devices^® microplate reader. ATP measurements were normalized based on the cell number per well, which was quantified using the cell counting kit‐8 (CCK‐8, CK04, Dojindo, Kumamoto Prefecture, Kyushu Island, Japan). 2.25. Mitochondrial isolation Mitochondria were isolated using the mitochondrial isolation kit (K256‐25, BioVision, San Francisco, CA, USA). Cells were homogenized 70 times using a Dounce homogenizer and suspended in buffer A containing a protease inhibitor. The homogenate underwent an initial centrifugation at 700 ×g for 10 min to eliminate cell debris and nuclei. The pellet obtained from this step was then further subjected to centrifugation at 3,000 ×g for 25 min to yield the mitochondria. For enhanced mitochondrial purity, the pellet was resuspended in buffer C and centrifuged (12,000 ×g) for 10 min. This final centrifugation step yielded the isolated mitochondria [[122]46]. 2.26. Measurement of mitochondrial membrane potential (MMP) The JC‐1 MMP assay kit (40706ES60, Yeasen, Shanghai, China) was employed to examine MMP. Cells were plated (1 × 10^5 cells per well) into a 6‐well culture plate, followed by a 20‐min incubation with JC‐1 dye working solution (1 mL per well) at 37°C. Carbonyl cyanide m‐chlorophenyl hydrazone, provided in the kit, was diluted to 50 µmol/L and added to the cell culture in a 1:1000 ratio as a positive control, with a subsequent 20‐min treatment. After the incubation period, the supernatant was discarded, and the cells were subjected to two washes in JC‐1 staining buffer (1×). Cells were then resuspended in 2 mL of serum‐ and phenol red‐containing culture medium. The detection of JC‐1 monomers was performed using a fluorescence or laser confocal microscope, set to excitation at 490 nm and emission at 530 nm, allowing for clear visualization of the monomers [[123]47]. 2.27. 3‐(4,5‐dimethylthiazol‐2‐yl)‐2,5‐diphenyltetrazolium bromide (MTT)‐based viability assay Cells were plated into 96‐well culture plates (3‐5 × 10^3 cells/well) and allowed to grow for 48 h. Afterward, MTT solution (10 mg/mL, ST316, Beyotime) was introduced to the cell suspension, followed by a further 4‐h incubation. Following dimethyl sulfoxide (DMSO) treatment, absorbance was measured using a UV‐visible spectrophotometer (Eppendorf, Hamburg, Germany) at OD 490 nm for tumor cells and OD 570 nm for CD8^+ T cells [[124]48]. 2.28. 5‐ethynyl‐2'‐deoxyuridine (EdU)‐based proliferation assay Cells were seeded in 24‐well plates, after which EdU ([125]C10337, Invitrogen) was introduced into the culture medium at 10 µmol/L. After a 2‐h incubation, the medium was aspirated, and the cells were immobilized with 4% PFA for 15 min. Cells were then permeabilized with 0.5% Triton X‐100, followed by staining with 50 µmol/L EdU (100 µL per well). DAPI was used for nuclear staining. Fluorescence microscopy (BX63, Olympus, Tokyo, Japan) was employed to visualize EdU‐positive cells [[126]49]. 2.29. Transwell assays for migration and invasion Transwell inserts (8 µm) were coated with Matrigel (E1270, Sigma‐Aldrich) for invasion assays, while migration assays were performed without the coating. Following the 48‐h incubation, cells were resuspended in a serum‐free medium to achieve 1 × 10^5 cells/mL. From this suspension, 200 µL was seeded (2 × 10^4 cells/well) and added to the upper compartment of each insert, while 800 µL of 20% FBS‐containing medium was placed in the lower insert. Following a 24‐h incubation at 37°C, the cells were immobilized in formaldehyde, stained with crystal violet, and any non‐migrated cells were removed. Invasive cells were imaged using an inverted microscope (CKX53, Olympus, Tokyo, Japan) and further analyzed using ImageJ software (v1.53, National Institutes of Health, Bethesda, MD, USA) [[127]50, [128]51]. 2.30. Terminal deoxynucleotidyl transferase dUTP nick end labeling (TUNEL)‐based apoptosis assay Samples were immobilized with 4% PFA (60536ES60, Zhongke Yasen Biotechnology, Shanghai, China) for 15 min, followed by permeabilization with 0.25% Triton X‐100 at RT for 20 min. Subsequently, the cells were cultured with 5% bovine serum albumin (BSA, 36101ES25, Zhongke Yasen Biotechnology) and stained using TUNEL staining reagent (C1086, Bestbio Biotechnology, Shanghai, China). DAPI was used to stain nuclei (blue fluorescence). Images of TUNEL‐positive cells, indicating apoptosis through green fluorescence, were captured using a fluorescence microscope (LSM 880, Carl Zeiss AG, Germany). The percentage of TUNEL‐positive cells was then determined from five randomly chosen fields to calculate the apoptotic rate [[129]52]. 2.31. Co‐culture of tumor cells and CD8^+ T cells A co‐culture experiment was conducted to explore the effect of different tumor cell densities on CD8^+ T cell activity, with NSUN2 expression levels remaining unaltered in the tumor cells. WT Hepa1‐6 cells were divided into low (1 × 10^4), medium (5 × 10^4), and high (1 × 10^5) densities. CD8^+ T cells, isolated from healthy C57BL/6J mice, were stimulated with CD3/CD28 activators (1 µg/mL, 11453D, Gibco) for 24 h prior to co‐culturing with the Hepa1‐6 cells at a 1:1 ratio. The co‐cultures were incubated at 37°C in a humidified environment with 5% CO[2] for 48 h. 2.32. Flow cytometric analysis of CD8^+ T cell function Following the co‐culture period, CD8^+ T cell function and apoptosis were examined via flow cytometry. The cells were labeled using an Annexin V‐FITC/PI double staining kit (556547, BD Bioscience, San Diego, CA, USA) as per the manufacturer's instructions. CD8^+ T cell viability and apoptosis were then evaluated using a flow cytometer (FACSCalibur, BD Bioscience, San Diego, CA, USA). 2.33. In vivo animal experiments The procedure for modeling with the Hepa1‐6 cell line and lentiviral vector adhered to the description provided earlier. For the H22 cell line, cells at a density of 1 × 10^8 cells/mL were directly injected into the left liver lobe of the mice. After 2 weeks, the engraftment of pretreated (i.e., pre‐KO) HCC cell lines, Hepa1‐6 and H22, was assessed in mice using the CRi Maestro In vivo imaging system (Cambridge Research & Instrumentation, Woburn, MA, USA) to examine the bioluminescent signal of firefly luciferase [[130]53]. For imaging of excised tissues, mice were euthanized via cervical dislocation under deep anesthesia using 2% isoflurane. After euthanasia, tumors were isolated and either processed for formalin‐immobilized, paraffin‐embedded tissue or rapidly immersed in liquid nitrogen to ensure thorough coverage by the freezing medium, minimizing bubble formation for subsequent analysis [[131]54]. The mice were randomly assigned to six groups, each comprising six mice: the NSUN2‐WT + NC‐oe group (injection of NSUN2‐WT cells infected with NC‐oe), NSUN2‐KO + NC‐oe group (injection of NSUN2‐KO cells infected with NC‐oe), and NSUN2‐KO + SOAT2‐oe group (injection of NSUN2‐KO cells infected with SOAT2 overexpression vector). Each group was subsequently transplanted with either Hepa1‐6 or H22 cells. 2.34. Robust rank aggregation (RRA) examination Six GEO transcriptomic datasets related to mouse HCC ([132]GSE39401, [133]GSE89689, [134]GSE101818, [135]GSE148379, [136]GSE159518, and [137]GSE200708) were sourced, as outlined in Supplementary Table [138]S4. The RRA method was employed to pinpoint reliable DEGs by integrating data from multiple microarray studies, thereby reducing inconsistencies. Prior to conducting the RRA analysis, gene lists indicating upregulation and downregulation were generated for each dataset based on the FC differences between control and HCC samples. The “Robust Rank Aggregation” R package was utilized to merge the DEG lists, with a significance threshold of P‐value < 0.05 and |logFC| > 1 used to filter for significant DEGs [[139]55]. 2.35. Immune cell infiltration analysis The CIBERSORT algorithm (R Script v1.03) was employed to estimate the relative proportions of 25 subtypes of immune cell in the samples. The interactions among different immune cell populations in HCC samples were explored using the CIBERSORT R script, which was downloaded from the CIBERSORT website [[140]56, [141]57]. 2.36. Identification of disease‐specific feature cells using the least absolute shrinkage and selection operator (LASSO) algorithm The glmnet R package was used to perform LASSO logistic regression, with the optimal number of feature cells determined by minimizing the lambda value [[142]58]. In the Oncomine database ([143]https://www.oncomine.org), we utilized the Gene Expression Module and the Immune Cell Analysis Module. First, the Gene Expression Module was used to extract NSUN2 and SOAT2 expression data in TCGA‐LIHC samples. Subsequently, the Immune Cell Analysis Module was employed to assess their correlations with CD8^+ T cell infiltration levels. Based on the median gene expression values in the TCGA‐LIHC dataset, samples were grouped into high and low expression categories, and the immune infiltration scores of CD8^+ T cells corresponding to high and low expression levels of NSUN2 and SOAT2 were presented [[144]59]. To further explore cell types influencing the tumor microenvironment (TME) during HCC progression, transcriptome sequencing data from the six mouse GEO datasets related to HCC were obtained. The datasets underwent batch correction and integration using the “limma” R package. The accuracy of batch correction was validated through density distribution and UMAP distribution analyses. 2.37. mRNA bisulfite sequencing (BS‐seq) and analysis of liver tissues in NSUN2‐KO mice Liver tissue samples of mice from the NSUN2‐WT and NSUN2‐KO groups (two from each) using Hepa1‐6 cells, were rapidly frozen in liquid nitrogen and stored at ‐80°C. These tissues were ground into powder, and total RNA was extracted from four samples using the Total RNA Isolation Reagent Kit (12183555, Invitrogen). Subsequently, mRNA was then captured with VAHTS mRNA Capture Beads (N401, Vazyme, Nanjing, Jiangsu, China), and strand‐specific mRNA libraries were generated using the VAHTS Stranded mRNA‐seq Library Prep Kit for Illumina V2 (NR612, Vazyme, Nanjing, Jiangsu, China). For methylation analysis, 200 ng of mRNA underwent bisulfite conversion with the EZ RNA Methylation Kit (R5002, Zymo Research, Orange County, CA, USA) with slight modifications to the standard protocol, including incubation at 70°C for 30 min and subsequent incubation at 60°C for 1 h. RNA integrity was confirmed using the RNA 6000 Pico Kit (5067‐1513, Agilent Technologies, Santa Clara, CA, USA), followed by next‐generation sequencing library construction with the VAHTS Stranded mRNA‐seq Library Prep Kit (NR612‐01/02, Vazyme, Nanjing, Jiangsu, China). Paired‐end sequencing (2 × 150 bp) was conducted by HaploX Genomics Center (Shenzhen, Guangdong, China) [[145]60]. The “limma” package in R software was utilized to pinpoint genes with differentially altered m5C levels in the BS‐seq data, with selection criteria set at |log2FC| > 1 and P‐value < 0.05. Data visualization was conducted using the “ggplot2” package for volcano plots and the“pheatmap” package for heatmaps. The m5C‐Atlas database ([146]http://180.208.58.19/m5c‐atlas/) was employed to predict potential m5C methylation sites in SOAT2 [[147]61, [148]62, [149]63, [150]64]. 2.38. Methylated RNA immunoprecipitation‐PCR (MeRIP‐PCR) analysis MeRIP was performed using the m5C MeRIP Kit (GS‐ET‐003, CLOUDSEQ, Shanghai, China) to detect m5C modifications. Total RNA (300 µg) was first fragmented into approximately 100 nucleotide‐long sequences. These fragments underwent incubation with protein A/G magnetic beads (88802, Thermo Fisher Scientific) pre‐coated with m5C antibodies (2.5 µg; ab10805, Abcam). The RNA fragments, bound to the m5C antibody‐coated beads, were mixed and incubated with gentle rotation for 2 h. The immunoprecipitated, methylated RNA fragments were then eluted and prepared for PCR analysis [[151]65, [152]66]. 2.39. Immunohistochemistry (IHC) Tissue samples from tumors were preserved overnight in 4% PFA, sectioned at 4 µm, and deparaffinized in xylene. The sections were rehydrated through a gradient of EtOH concentrations, and antigen retrieval was achieved by boiling them in 0.01 mol/L citrate buffer for 15‐20 min. Endogenous peroxidases were quenched with 3% H[2]O[2] for 30 min, followed by a 20‐min blocking with goat serum. The sections were incubated with primary Ki67 antibody at RT for 1 h, and subsequently treated with goat anti‐rabbit secondary antibody IgG at 37°C for 20 min. Please refer to Supplementary Table [153]S5 for the antibody list. Streptavidin‐peroxidase (SP) was introduced and allowed to incubate the samples at 37°C for 30 min. Next, diaminobenzidine (DAB; P0202, Beyotime) was applied for chromogenic detection, and sections were counterstained with hematoxylin (C0107, Beyotime). The slides were dehydrated and mounted for microscopic analysis. The Ki67 positivity rate was assessed by counting the number of positive cells in five randomly selected high‐power fields on each slide [[154]67]. 2.40. Flow cytometry analysis Target cells were stained with rat anti‐mouse antibodies (Supplementary Table [155]S6), incubated at 4°C for 30 min without light exposure, and centrifuged at 1,500 ×g and 4°C for 10 min to discard the supernatant. The cells were then immobilized with 2% PFA (30525‐89‐4, Sigma‐Aldrich) and analyzed within 24 h using a FACS Aria II flow cytometer (BD Bioscience, San Diego, CA, USA). Data were processed using Flowjo CE and analyzed with FlowJo software [[156]68]. 2.41. Enzyme‐linked immunosorbent assay (ELISA) The levels of interferon‐gamma (IFN‐γ) and granzyme B (GZMB) were quantified in mouse tumor tissue and In vitro cell models using ELISA kits (Invitrogen) as per the manufacturer's protocols: IFN‐γ (BMS606‐2) and GZMB (BMS6029) [[157]69, [158]70]. 2.42. Gene expression quantification Total RNA was extracted from tissues or cells using Trizol reagent (15596026, Invitrogen), followed by concentration and purity assessment using the NanoDrop LITE spectrophotometer (ND‐LITE‐PR, Thermo Fisher Scientific). RNA was reverse‐transcribed into cDNA using the PrimeScript RT reagent Kit (RR047Q, TaKaRa, Kyoto, Japan). Next, RT‐qPCR was conducted with SYBR Green PCR Master Mix (4309155, Applied Biosystems, Madison, WI, USA) on the ABI PRISM 7500 system. Primers were custom‐designed and synthesized by TaKaRa (Supplementary Table [159]S7), with GAPDH as the housekeeping gene. Gene expression was calculated using the 2^−ΔΔCt method [[160]71, [161]72, [162]73]. 2.43. Protein expression quantification Lysis of tissues or cells was performed using enhanced radioimmunoprecipitation assay (RIPA) buffer containing protease inhibitors (P0013B, Beyotime). Protein concentration was measured with a BCA kit (P0012, Beyotime). Proteins were separated via 10% sodium dodecyl sulfate‐polyacrylamide gel electrophoresis (SDS‐PAGE) and transferred onto polyvinylidene difluoride (PVDF) membranes. Membranes were blocked with 5% BSA at RT for 2 h. Samples were incubated overnight with diluted rabbit anti‐mouse primary antibodies (Supplementary Table [163]S8) and then with horseradish peroxidase (HRP)‐conjugated goat anti‐rabbit secondary antibody (ab6721, Abcam) at RT for 1 h. Blots were visualized through Pierce™ ECL substrate (32209,Thermo Fisher Scientific) and captured using the Bio‐Rad image analysis system (Bio‐Rad, Hercules, CA, USA). Band intensity was quantified using Image J, with GAPDH as the loading control [[164]72]. 2.44. Statistical analysis R language (version 4.2.1) and RStudio (version 2022.12.0‐353) were used for statistical analysis, with Perl (version 5.30.0) for file processing and GraphPad Prism (version 8.0) for visualization. Data are reported as mean ± standard deviation. Pearson correlation analysis was conducted to assess the relationships between sequencing depth and other metrics. Student's t‐tests were used for two‐group comparison, while one‐way analysis of variance (ANOVA) was employed for multigroup comparisons. Two‐way ANOVA was applied for analyzing data across different time points, with post‐hoc Bonferroni correction applied. The significance threshold was set at P < 0.05. Multivariate Cox regression analysis was conducted to pinpoint prognostic factors for HCC (P < 0.05). LASSO regression was used using the glmnet R package (v4.1‐6, [165]https://glmnet.stanford.edu/) to reduce data dimensionality and select key variables, with 10‐fold cross‐validation (cv.glmnet function) applied to minimize model error and determine the optimal lambda. Additionally, the support vector machine recursive feature elimination (SVM‐RFE) method was used to identify HCC‐specific feature genes via the e1071 R package (v1.7‐13, [166]https://cran.r‐project.org/web/packages/e1071/). Cross‐validation determined the number of features, and top‐ranked genes were selected as key variables. 3. RESULTS 3.1. scRNA‐seq analysis for differences in cell types and quantities in HCC tissues To explore the progression of HCC and alterations in the TME, an HCC mouse model was established (Figure [167]1A), and liver tissues were analyzed. The results confirmed the success of the HCC mouse model establishment (Figure [168]1B, C). The histological examination further validated the model, with HE staining revealing intact liver tissue structure adjacent to the tumor and no pronounced pathological changes (Figure [169]1D). These results confirmed that the model successfully replicated HCC development with cirrhosis, while maintaining liver health in other regions. FIGURE 1. FIGURE 1 [170]Open in a new tab Cell clustering and annotation of scRNA‐seq data. (A) Schematic diagram illustrating the establishment of the HCC mouse model. (B) Representative histological images of the liver from normal mice and HCC model mice. Black arrows indicate tumor lesions. (C) Tumor volume changes measured at different time points. (D) HE‐stained images showing the normal liver tissue structure adjacent to the tumor. Black arrows indicate tumor lesions. (E) Highly variable genes were selected through variance analysis, and the gene names of the top 10 highly variable genes were annotated. (F) PCA results showing the distribution of cells on PC1 and PC2. Each point represents a cell. C1 and C2 represent two adjacent normal tissue samples, while T1 and T2 represent two HCC tissue samples. (G) Distribution of standard deviations of PCs. Important PCs have larger standard deviations. (H) Visualization of clustering results using UMAP, showing the clustering and distribution of cells. (I) Expression pattern of known cell lineage‐specific marker genes in different cell types. (J) Three‐dimensional UMAP clustering results showing the clustering and distribution of cells. (K) The difference in cell abundance between normal adjacent tissue samples and HCC samples. (A‐D) HCC and adjacent normal tissues from 6 mice; (E‐K) HCC and adjacent normal tissues randomly selected from 2 of the 6 mice used in (A‐D). Abbreviations: HCC, hepatocellular carcinoma; scRNA‐seq, Single‐cell RNA sequencing; HE, hematoxylin‐eosin; PCA, principal component analysis; UMAP, Uniform manifold approximation and projection; PCs, principal components; TME, tumor microenvironment. Tumor tissues and adjacent normal tissues from two HCC mice were then subjected to scRNA‐seq. After applying quality control criteria, low‐quality cells and redundant genes were excluded, yielding an expression matrix of 14,899 genes and 17,898 cells (Supplementary Figure [171]S1A). nCount_RNA was negatively correlated with percent.mt and positively correlated with nFeature_RNA, but not correlated with percent.HB (Supplementary Figure [172]S1B), confirming the robustness of the filtered data. For further analysis, the top 2,000 highly variable genes were selected for downstream processing (Figure [173]1E). The CellCycleScoring function was utilized to determine the cell cycle phases of the sample cells, and the results showed significant distribution in the G1, S, and G2/M phases (Supplementary Figure [174]S1C). Following preliminary data normalization, linear dimensionality reduction was then performed using PCA, successfully separating the cell populations based on gene expression variability (Figure [175]1F). The gene expression heatmap for PC1–PC6 demonstrated significant differences in key gene sets across different principal components, reflecting the biological diversity present in the dataset (Supplementary Figure [176]S1D). After batch correction using the Harmony package, the clustering of cells across samples was significantly improved (Supplementary Figure [177]S1E, F). Additionally, ElbowPlot was used to rank the standard deviations of the PCs, with the results indicating that PC1–PC30 effectively captured the main variance of highly variable genes (Figure [178]1G), providing strong support for subsequent analysis. The results indicated that the preprocessed data effectively reflected the biological diversity of the samples while successfully addressing technical variation. Subsequently, clustering analysis was conducted, resulting in the identification of 20 distinct clusters, each characterized by specific marker gene expression profiles (Figure [179]1H). Cell lineage‐specific marker genes were identified through a literature search, and the clusters were annotated using the CellMarker online database (Figure [180]1I, Supplementary Figure [181]S1G). This analysis revealed eight distinct cell categories (Figure [182]1J). Relative to adjacent normal tissues, HCC tissues exhibited a pronounced decrease in neutrophils and stromal cells, whereas the proportion of epithelial cells or cancer cells was markedly elevated (Figure [183]1K). Non‐linear dimensionality reduction was performed using the UMAP algorithm, and the analysis revealed that the identified clusters were clearly separated in spatial distribution, with minimal overlap. Clustering at a resolution of 0.4 further revealed the hierarchical structure and specific transcriptional characteristics of the cell subpopulations. Each cluster exhibited unique connectivity and proportions, reflecting the well‐differentiated state of the cell subgroups (Supplementary Figure [184]S2). These findings demonstrate that HCC samples can be divided into 20 clusters representing eight distinct cell subpopulations. Notably, HCC tissues showed a pronounced reduction in neutrophils and stromal cells, alongside a marked increase in epithelial cells or cancer cells. 3.2. Characterization of cell subtypes and immune escape of tumor cells in HCC samples through CNV and trajectory analysis Cell subtype alterations in HCC samples were investigated by validating the cell annotation, which confirmed their accuracy (Supplementary Figure [185]S3). The “inferCNV” package was then utilized to detect large‐scale CNVs, with dendritic cells serving as controls. The majority of epithelial cells exhibited CNVs (Figure [186]2A). Normal cells, epithelial cells, or cancer cells were clustered into 12 classes based on CNV profiles, referred to as “Observation” in the inferCNV package. After excluding Class 6, 11,965 malignant cells were identified (Figure [187]2B, C). FIGURE 2. FIGURE 2 [188]Open in a new tab Results of malignant cell extraction and cell communication analysis. (A) CNV analysis of scRNA‐seq sequencing data, with dendritic cells as the reference. The yellow color in the observation represents normal dendritic cells, while the cyan color represents epithelial cells or cancer cells. Blue indicates DNA copy loss, while red indicates DNA copy gain. (B) K‐means clustering algorithm for clustering the observation cells in Figure A. (C) CNV scores after K‐means clustering. (D) UMAP plot of malignant cell clusters. (E) Differentiation trajectory plot of malignant epithelial cells analyzed by Monocle2. (F) Distribution of cells in different states along the pseudotime. (G) Reordered gene ranking after filtering out genes with a q‐value greater than 0.05 in malignant cells. The plot shows the relationship between gene expression variability and mean expression, with the significant genes above the red line. (H) Cell‐cell communication differences between adjacent normal tissue samples (n = 2) and tumor samples (n = 2). The thickness of the lines on the left represents the number of pathways, while the thickness of the lines on the right represents the interaction strength. Red lines indicate increased signals in the tumor group compared to the normal group, and blue lines indicate decreased signals in the tumor group. (I) Heatmap of differences in interaction quantity or interaction strength. The colored bars at the top represent the sum of the column values (incoming signals), while the colored bars on the right represent the sum of the row values (outgoing signals). In the color bar, red (or blue) indicates increased (or decreased) signals in the Tumor group compared to the Normal group. “Relative values” represent the relative interaction signals between different cell types, comparing the quantity and strength of the interaction signals using a normalization method. (J) Total count of information flow between adjacent normal and tumor tissues. (K) Calculation of differential information flow between adjacent normal and tumor tissues, defined as the sum of communication probabilities between all cell pairs in the inferred network (i.e., total weight in the network). Red signal pathways are enriched in the Normal group, while green signal pathways are enriched in the Tumor group. Abbreviations: scRNA‐seq, Single‐cell RNA sequencing; CNV, copy number variation; UMAP, Uniform manifold approximation and projection. Furthermore, these malignant cells were further clustered into eight classes (Figure [189]2D, Supplementary Figure [190]S4). Trajectory analysis using the “Monocle2” package highlighted the dynamic progression and heterogeneity of these malignant cells, with State 1 cells at the starting point of the pseudotime trajectory axis (Figure [191]2E, F). Differential expression gene analysis along the pseudotime axis identified key genes with q values below 0.05 (Figure [192]2G). The top 50 significantly altered genes (P < 0.05) were extracted, with 42 of these genes specifically associated with State 1, forming the core of gene set A (Supplementary Figure [193]S5). These results highlight the unique transcriptional features of State 1 in the malignant cell process, suggesting its key role in the initiation of tumor progression and the early molecular events mediating cell transformation. To further investigate the functional differences between cell types in the tumor and normal groups, pathway activities were explored using the R package “CellChat”. Figure [194]2H illustrates the number of cell‐cell communications and interaction strengths among the eight cell subtypes. Notably, interactions between epithelial cells/cancer cells and T cells were relatively high in terms of number and strength. However, a comparison between the groups revealed a prominent weakening of these interactions in the tumor group versus the normal group (Figure [195]2I). Further analysis focused on the total number and differences in information flow between the groups. The tumor group demonstrated a pronounced reduction in total information flow relative to the normal group (Figure [196]2J). Specifically, there was a marked decrease in information flow related to the CD45 pathway among immune cells, whereas information flow related to the junctional adhesion molecule (JAM) pathway, which plays a role in immune escape by cancer cells, markedly increased (Figure [197]2K). JAMs are essential for maintaining cell‐cell junctions and have been shown to contribute to the immune escape of cancer cells, as supported by recent studies in HCC [[198]74, [199]75, [200]76]. These results suggest that tumor cells may drive cancer progression through immune escape. The results indicate that the malignancy of cells in HCC samples exhibits significant heterogeneity, with distinct transcriptional states and intercellular interactions compared to normal samples. The reduction in CD45‐related immune interactions and the increase in JAM‐related adhesion pathways suggest that tumor cells may evade immune responses and enhance invasive characteristics by remodeling the tumor microenvironment. 3.3. Identification and validation of key genes and proteins associated with HCC progression Key genes associated with HCC were identified through bulk RNA‐seq with HCC tissues and adjacent normal tissues from three HCC‐bearing mice. Differential analysis of the sequencing data revealed 89 genes that were differentially expressed in HCC tissues versus adjacent normal tissues (Figure [201]3A). Disease Ontology Enrichment analysis linked these 89 genes to several liver diseases, such as hepatitis and hepatoblastoma (Supplementary Figure [202]S6A). GO enrichment analysis showed that the DEGs were primarily enriched in molecular functions such as proteoglycan binding, RNA transmembrane transport activity, and RNA methyltransferase activity (Supplementary Figure [203]S6B). Additionally, KEGG pathway analysis indicated that DEGs were mainly involved in pathways associated with fatty acid metabolism, HCC cells, and glycolysis/gluconeogenesis (Supplementary Figure [204]S6C). These results indicate the strong connection between these DEGs and liver diseases as well as various metabolic processes. FIGURE 3. FIGURE 3 [205]Open in a new tab Integration of bulk RNA‐seq and proteomics in a meta‐analysis to identify key factors associated with HCC. (A) Volcano plot depicting differential expression of mRNA between three samples of HCC tissue and their matched normal tissue from mice using bulk RNA‐seq data. (B) The trend of gene regression coefficients under different lambda values in the LASSO regression process. The x‐axis (L1 Norm) represents the sum of the absolute values of the regression coefficients and increases as lambda decreases. As lambda increases, more genes retain nonzero coefficients, indicating their importance in the model. The gene set with the smallest lambda value is selected as the candidate gene set. (C) The process of determining the optimal lambda value for the LASSO model through cross‐validation. The red dashed line marks the point where lambda is smallest (i.e., the point of minimal cross‐validation error), and the corresponding gene set is considered to have the best predictive ability. The numbers above the plot denote the number of nonzero coefficients (i.e., selected variables) for each lambda value. As lambda decreases (moving to the right), the number of selected variables gradually increases, indicating that the model includes more features. (D) The feature genes were selected by the SVM‐RFE algorithm. As the number of feature genes decreases, RMSE initially fluctuates but then decreases, reaching its lowest value when 36 key genes remain (marked in blue). The final 36 key genes with optimal predictive ability were retained as feature genes. (E) Venn diagram showing the intersection between HCC‐relevant mRNAs selected by LASSO regression and SVM‐RFE and DEGs obtained from bulk RNA‐seq of tumor tissues and paired normal tissues from 3 HCC‐bearing mice. (F) Heatmap revealing the expression levels of Nsun2 and Cdkn3 in the sequencing data. (G) Proportion of DEPs between three samples of HCC tissue and their matched normal tissue in mice. (H) Comparison of the number of upregulated and downregulated DEPs in HCC tissues. (I) Candidate DEP interactome network, with lines between circles represent PPIs. Blue nodes represent downregulated proteins, while red and orange nodes represent upregulated proteins. (J) Core protein adjacency node counts in the gene interaction network. (K) Heatmap illustrating the differential expression of 25 DEPs in the proteomics data between HCC and adjacent normal tissues from 3 mice. (L) Venn diagram showing the intersection between Gene set A, Gene set B, and Protein set C. (M) Forest plot comparing NSUN2 expression between tumor and normal groups. (N) Forest plot comparing NSUN2 expression between ethnicity subgroups within the tumor and normal groups. (O) Funnel plot for assessing publication bias. (P) Sensitivity analysis conducted using the one‐by‐one exclusion method. In the volcano plot, blue dots represent markedly downregulated mRNA in HCC tissue, red dots represent markedly upregulated mRNA in HCC tissue, and grey dots represent mRNA with no significant difference in expression. Abbreviations: HCC, hepatocellular carcinoma; SMD, standardized mean difference; CI, confidence interval; LASSO, Least Absolute Shrinkage and Selection Operator; SVM‐RFE, Support Vector Machine‐Recursive Feature Elimination; DEPs: differentially expressed proteins; RMSE, root mean square error. Furthermore, reanalysis of the scRNA‐seq data revealed a notable upregulation of EPC‐related marker genes, CD34 and VEGFR2, in NSUN2‐overexpressing HCC samples (Supplementary Figure [206]S6D, E). This finding suggests that NSUN2 might promote the recruitment of EPCs, thereby facilitating tumor angiogenesis and remodeling of the TME. Additionally, pathway analysis uncovered that under NSUN2 overexpression, EPC‐related genes were markedly associated with multiple immune escape escape‐related signaling pathways (Supplementary Figure [207]S6F). These results indicate that NSUN2 not only promotes HCC progression by regulating the metabolic activity of tumor cells but also potentially supports tumor growth and immune escape escape by modulating EPC activity. Subsequently, multivariate Cox regression analysis was conducted using LASSO regression. As the penalty parameter lambda increased, the regression coefficients of the genes gradually approached zero (Figure [208]3B). The optimal lambda value was selected through cross‐validation (Figure [209]3C). This approach effectively reduced the data dimensionality and identified a set of genes with the best predictive ability. Meanwhile, the SVM‐RFE analysis method was used to extract HCC‐specific feature genes. The optimal model was achieved when 36 key genes remained (Figure [210]3D). Finally, the intersection of DEGs identified through bulk RNA‐seq and those selected by LASSO regression and SVM‐RFE analysis resulted in gene set B, which included two key mRNAs: NSUN2 and Cyclin‐dependent kinase inhibitor 3 (CDKN3) (Figure [211]3E). The expression patterns of NSUN2 and CDKN3 (Figure [212]3F) highlighted these genes as potential biomarkers due to their marked differential expression between tumor and normal tissues, suggesting their important roles in HCC progression. To further pinpoint key proteins involved in the occurrence and progression of HCC, TMT‐labeled quantitative proteomic analysis was employed to compare HCC tissues with adjacent normal tissues in the three HCC‐bearing mice and identified 10,249 proteins in the HCC tissue samples, among which 466 DEPs were detected. Of these, 147 proteins were upregulated, while 319 were downregulated (Figure [213]3G, H). The top 50 DEPs, ranked by the absolute value of log2FC, were selected, and unconnected nodes were excluded. A PPI network was then constructed, highlighting 19 protein nodes with high degrees of connectivity (Figure [214]3I). The network analysis involved calculating the number of adjacent nodes for each gene, followed by the removal of proteins with zero adjacent nodes, which resulted in a refined protein set C, comprising 29 proteins (Figure [215]3J). Among them, only TOP2A and NSUN2 were upregulated, while the other proteins exhibited downregulated expression (Figure [216]3K). According to the intersection of gene sets A and B with the selected protein set C, one key factor, NSUN2, was pinpointed (Figure [217]3L). The scRNA‐seq data analysis indicated a marked upregulation of NSUN2 in tumor cells relative to immune cells (Supplementary Figure [218]S6G). Specifically, the FeaturePlot revealed that NSUN2 expression was predominantly localized in tumor cells, with notably low levels detected in immune cells (Supplementary Figure [219]S6H). These results highlight the selective high expression of NSUN2 in HCC tumor cells, further implicating its role in facilitating tumor proliferation and immune escape. To validate the role of NSUN2 in HCC, nine transcriptome mRNA sequencing data were retrieved from the nine human GEO datasets, followed by a meta‐analysis to confirm the differential expression of NSUN2. Given the statistical heterogeneity detected among the studies (I^2 = 52%, P = 0.030), a random‐effects model was employed. The meta‐analysis revealed that NSUN2 expression was markedly elevated in the tumor group relative to the normal group (SMD = 1.12, 95% CI = 0.79‐1.45) (Figure [220]3M). Subgroup analysis based on the ethnicity of the study subjects revealed that in the Asian population, NSUN2 expression was markedly higher in the tumor group relative to the normal group (SMD = 1.09, 95% CI = 0.74‐1.44). Conversely, in the American population, the differential expression of NSUN2 between the tumor and normal groups was not statistically significant (Figure [221]3N). A funnel plot was utilized to evaluate publication bias, which showed that two studies did not align with the symmetry of the funnel plot. Nevertheless, Egger's test revealed no statistically significant difference (P > 0.05), indicating the absence of publication bias in this analysis (Figure [222]3O). Moreover, sensitivity analysis was conducted using the leave‐one‐out sensitivity analysis to ensure the reliability of the included study results, with the SMD values showing no significant change (Figure [223]3P), thus confirming the robustness of the meta‐analysis results. In summary, NSUN2 was identified as a critical factor, playing a substantial role in the onset of HCC. 3.4. NSUN2 facilitated the malignant phenotypes of HCC cells NSUN2, a widely expressed protein in mammalian cells, belongs to the RNA methyltransferase family. This protein, present in both the cytoplasm and nucleus, plays an essential role in transcription and RNA methylation modification [[224]77]. Existing research has shown that NSUN2 is overexpressed in various tumors, including breast and colorectal cancers, and is linked to a range of malignant phenotypes [[225]78, [226]79]. To confirm the association between NSUN2 and HCC development, RT‐qPCR analysis was conducted on HCC and adjacent normal tissues from 30 HCC patients. The clinicopathologic characteristics of these patients are summarized in Supplementary Table [227]S9. NSUN2 was found to be markedly higher in HCC tissues relative to adjacent normal tissues (Figure [228]4A). FIGURE 4. FIGURE 4 [229]Open in a new tab Effects of NSUN2 on the biological functions of HCC. (A) RT‐qPCR was performed to evaluate the gene expression levels of NSUN2 in HCC and adjacent normal tissues collected from 30 HCC patients. (B) Schematic diagram of CRISPR/Cas9 gene‐editing technology used to construct NSUN2‐KO HCC cells. (C‐D) Expression levels of NSUN2 in CRISPR/Cas9 gene‐edited NSUN2‐KO cells detected by RT‐qPCR and Western blotting. (E) Schematic diagram of HCC cells overexpressing NSUN2 through lentivirus transfection. (F‐G) Expression levels of NSUN2 in cells overexpressing NSUN2 plasmid detected by RT‐qPCR (F) and Western blotting after lentivirus transfection (G). (H) Cell viability of Hepa1‐6 cells in each group measured by MTT assay. (I) Proliferation ability of Hepa1‐6 cells in each group evaluated by EdU assay. (J‐K) Migration and invasion capacities of Hepa1‐6 cells in each group determined by Transwell experiment. (L) Apoptosis of Hepa1‐6 cells in each group analyzed by TUNEL assay. Cell experiments were performed in triplicate. Abbreviations: HCC, hepatocellular carcinoma; RT‐qPCR, reverse transcription quantitative polymerase chain reaction; KO: knock out; MTT: Methylthiazolyldiphenyl‐tetrazolium bromide; EdU: 5‐ethynyl‐2`‐deoxyuridine; TUNEL: TdT‐mediated dUTP Nick‐End Labeling; NSUN2, NOP2/Sun RNA methyltransferase family member 2. CRISPR/Cas9 gene‐editing technology was employed to create NSUN2‐KO Hepa1‐6/H22 cells (Figure [230]4B). RT‐qPCR and WB assays were employed to measure NSUN2 expression levels in monoclonal cells, and clones without detectable NSUN2 expression were selected for subsequent expansion and experiments (Figure [231]4C, D, Supplementary Figure [232]S7A, B). Concurrently, NSUN2‐oe Hepa1‐6/H22 cells were generated using lentiviral constructs (Figure [233]4E). The transfection efficiency of NSUN2 was confirmed through RT‐qPCR and WB assay, demonstrating a pronounced increase in NSUN2 expression after transfection (Figure [234]4F, G, Supplementary Figure [235]S7C, D). Malignant phenotypes of HCC cells were assessed in response to NSUN2‐oe and NSUN2‐KO. Relative to the NC‐oe group, the NSUN2‐oe group exhibited markedly increased cell viability, proliferation, migration, and invasion abilities, along with a prominent reduction in apoptosis. Conversely, when compared to the NSUN2‐WT group, the NSUN2‐KO group exhibited markedly diminished cell viability, proliferation, migration, and invasion capacities, along with a notable increase in apoptosis (Figure [236]4H–L, Supplementary Figure [237]S7E–I). Taken together, these findings demonstrate that NSUN2 enhances the malignant phenotypes of HCC cells. 3.5. Impact of NSUN2 on immune cell interactions and tumor immune escape in HCC As illustrated in Figure [238]2H–J, a notable reduction was noted in interactions between epithelial cells or cancer cells and immune cells in HCC tissues relative to adjacent normal tissues, suggesting a possible association with immune escape by HCC cells. Existing literature consistently demonstrates the pivotal role of the TME in HCC progression, wherein cancer cells frequently adopt strategies such as immune escape, immunotolerance, or immune suppression to evade immune cell attacks [[239]80, [240]81, [241]82]. To delve deeper into the cell types that markedly influence the TME during the progression of HCC, transcriptome sequencing data were obtained from the six mouse GEO datasets (Figure [242]5A). The data were subjected to batch correction and integrated using the “limma” package in R (Supplementary Figure [243]S8A), with the accuracy of the batch‐corrected results validated through density distribution and UMAP distribution analyses (Supplementary Figure [244]S8B, C). FIGURE 5. FIGURE 5 [245]Open in a new tab Analysis of immune infiltration and experimental validation results. (A) Presentation of mouse sample information for 6 GEO datasets. A total of 77 samples (23 adjacent normal samples and 54 HCC samples) were collected from six mouse GEO datasets. (B) Distribution of immune cell proportions in 23 adjacent normal samples and 54 HCC samples from the 6 mouse GEO datasets. (C) LASSO analysis results. The dashed line indicates the log(λ) value corresponding to the optimal Binomial Deviance and the number of cells retained. The red dashed line indicates the optimal Binomial Deviance value corresponding to log(λ). The gray dashed lines represent the upper and lower confidence intervals. (D) LASSO analysis results of different cell types. (E) Comparison of the levels of 22 immune cells between the adjacent normal tissue and HCC samples. (F) Venn diagram showing the feature immune cells calculated by the LASSO regression model and the different cells in the adjacent normal tissue and HCC samples. (G) Measurement of Hepa1‐6 cell activity by MTT assay. (H) Mechanism diagram showing the CD8^+ T cell‐mediated killing of tumor cells. (I) Detection of GZMB and IFN‐γ content in Hepa1‐6 cell culture medium by ELISA.(J‐K) Quantification of IFN‐γ and GZMB positive CD8^+ T cells by flow cytometry. Cell experiments were performed in triplicate. Abbreviations: HCC, hepatocellular carcinoma; GEO, Gene Expression Omnibus; LASSO, Least Absolute Shrinkage and Selection Operator; MTT: Methylthiazolyldiphenyl‐tetrazolium bromide; ELISA: Enzyme linked immunosorbent assay. Subsequently, the CIBERSORT algorithm was used to calculate the proportions of 25 immune cell types infiltrating HCC samples in an integrated dataset composed of the six mouse GEO datasets. The relative distribution of immune cell proportions between the normal group and the HCC group were compared. Overall, M0 macrophages, M1 macrophages, Th1 cells, and Th17 cells were significantly increased in HCC tissues, while CD8+ T cells and monocytes also exhibited notable changes (Figure [246]5B), indicating a significant remodeling of the immune cell composition within the TME during HCC progression. Based on the CIBERSORT results, we used the LASSO regression model to perform LASSO regression on 25 immune cell subtypes from CIBERSORT results (Figure [247]5C), identifying 9 key immune cell subtypes associated with HCC (Figure [248]5D). Additionally, we compared the differences in immune cell infiltration between HCC and adjacent normal tissues, resulting in eight immune cell subtypes with significant differences in infiltration levels, including activated CD8^+ T cells, M0 macrophages, M1 macrophages, naïve CD4^+ T cells, Th1 cells, Th17 cells, monocytes, and immature DCs (Figure [249]5E). Integrating the findings from the LASSO regression model, five cell subtypes (activated CD8^+ T cells, M1 macrophages, Th17 cells, monocytes, and immature DCs) were identified in the intersection (Figure [250]5F). Synthesizing these observations with the data from Figure [251]2 and earlier experiments led to the hypothesis that NSUN2, an HCC‐related gene, might facilitate immune escape by impairing the efficacy of CD8^+ T cells against epithelial cells or cancer cells, thereby promoting HCC development. This hypothesis is supported by previous research indicating the overexpression of NSUN2 in HCC tissues and its involvement in RNA methylation [[252]15, [253]16, [254]17], which may play a role in immune escape. To investigate the impact of NSUN2 expression in cancer cells on CD8^+ T cells, the MTT assay results demonstrated that CD8^+ T cell viability was markedly increased in the NSUN2‐KO group versus the NSUN2‐WT group, while it was markedly diminished in the NSUN2‐oe group compared to the NC‐oe group (Figure [255]5G, Supplementary Figure [256]S9A). Previous evidence has confirmed that the secretion of GZMB and IFN‐γ by CD8^+ T cells is crucial for their capacity to recognize and eliminate tumor cells (Figure [257]5H) [[258]83]. In the present study, ELISA results indicated that GZMB and IFN‐γ levels were increased in the NSUN2‐KO group, whereas a marked decrease was observed in the NSUN2‐oe group (Figure [259]5I, Supplementary Figure [260]S9B). Flow cytometry analysis further corroborated these findings, revealing a prominent increase in IFN‐γ‐ and GZMB‐positive cells in the NSUN2‐KO group and a notable reduction in the NSUN2‐oe group (Figure [261]5J, K, Supplementary Figure [262]S9C, D). In conclusion, these preliminary findings suggest that NSUN2‐KO in HCC cells amplifies the activity and cytotoxicity of CD8^+ T cells, while NSUN2 overexpression diminishes these effects. These results indicate that NSUN2 may play a pivotal role in modulating the immune response within the HCC microenvironment, potentially facilitating immune escape by the tumor. 3.6. Effect of different tumor cell density on CD8^+ T cell activity and apoptosis The study shifted to verify whether different tumor cell numbers could independently affect CD8^+ T cell activity regardless of NSUN2 expression. We conducted co‐culture experiments using Hepa1‐6 cells and CD8^+ T cells without altering NSUN2 expression. Flow cytometry analysis revealed a pronounced reduction in CD8^+ T cell viability and elevation in CD8^+ T cell apoptosis with increasing tumor cell densities (Supplementary Figure [263]S10A, B), indicating that tumor cell density could affect CD8^+ T cell activity independently of NSUN2 levels. In the NSUN2‐oe and NSUN2‐KO groups, NSUN2 overexpression markedly repressed CD8^+ T cell activity relative to the control group, while NSUN2‐KO enhanced the cytotoxicity of CD8^+ T cells (Supplementary Figure [264]S10C, D). These findings, coupled with the new experimental results, further confirm the critical role of NSUN2 in immune escape in HCC, not only through the regulation of tumor cell metabolism but also by directly modulating CD8^+ T cell activity to support immune escape. 3.7. Functional study of NSUN2‐mediated m5C modification in enhancing mRNA transcription of SOAT2 NSUN2, a key m5C RNA methyltransferase, has been extensively studied for its capacity to enhance the transcription of oncogenic genes through m5C modifications, thereby facilitating cancer progression [[265]15, [266]84]. To elucidate the downstream targets of NSUN2, mRNA BS‐seq was conducted on HCC tissues derived from an in situ HCC mouse model harboring NSUN2‐WT and NSUN2‐KO Hepa1‐6 cells. This approach was employed to assess total m5C content in NSUN2‐WT and NSUN2‐KO tumor tissues. The results demonstrated a notable reduction in total m5C levels after NSUN2‐KO (Figure [267]6A). The m5C modification rates of specific RNA sequences were examined in HCC tissues. It was uncovered that the m5C levels of SOAT2 in NSUN2‐KO mouse HCC tissues were markedly decreased versus the NSUN2‐WT group (Figure [268]6B), suggesting that SOAT2 is a likely downstream target of NSUN2‐mediated m5C modification. Furthermore, predictions from the m5C‐Atlas database identified two m5C sites on SOAT2 (Figure [269]6C), indicating that SOAT2 can undergo m5C methylation mediated by m5C‐modifying enzymes. FIGURE 6. FIGURE 6 [270]Open in a new tab m5C‐seq analysis results and experimental validation. (A) Detection of relative m5C levels in tumor tissues from different groups of mice (2 mice/group). (B) Heatmap showing the relative m5C levels of representative mRNAs with significant expression changes in tumor tissues from NSUN2‐KO and NSUN2‐WT mice. (C) Prediction results from the m5C‐Atlas database for potential m5C methylation sites on SOAT2. (D) Immune infiltration scores of CD8^+ T cells corresponding to the expression levels of NSUN2 and SOAT2, grouped according to the median gene expression values from the TCGA‐LIHC dataset composed of 419 samples. (E) Scatter plot analysis of the correlation between the expression of NSUN2 and SOAT2 in m5C‐sequencing data from different groups of mice (2 mice/group). (F‐G) MeRIP‐qPCR and RT‐qPCR were utilized to detect the m5C and mRNA levels of SOAT2 in HCC and adjacent normal tissues collected from 30 HCC patients. (H‐K) Expression levels of NSUN2 and SOAT2 in different groups of cells detected by RT‐qPCR and Western blotting. (L) Stability of SOAT2 in Hepa1‐6 cells of different groups tested by Actinomycin D experiment. (M) Impact of NSUN2 knock out or overexpression on the m5C levels of SOAT2 in Hepa1‐6 cells analyzed using MeRIP‐qPCR. (N‐O) Cholesterol levels in various groups of Hepa1‐6 cells detected by metabolomics. (P) Quantitative analysis of cholesterol levels in the cell culture medium. Cell experiments were repeated 3 times. Abbreviations: m5C, 5‐methylcytosine; MeRIP, methylated RNA immunoprecipitation; RT‐qPCR, reverse transcription quantitative polymerase chain reaction; NSUN2, NOP2/Sun RNA methyltransferase family member 2; KO: knock out; SOAT2, sterol O‐acyltransferase 2. Furthermore, the impact of NSUN2 and SOAT2 expression on the immune infiltration of CD8^+ T cells in TCGA‐LIHC was analyzed using the Oncomine database. The results indicated that high expression of NSUN2 and SOAT2 was associated with a pronounced reduction in CD8^+ T cell enrichment scores (Figure [271]6D). Additionally, a correlation analysis of NSUN2 and SOAT2 expression in transcriptome sequencing data revealed a positive correlation between these two genes (Figure [272]6E). SOAT2, also known as acetyl‐CoA acetyltransferase 2 (ACAT2), is involved in lipid metabolism by facilitating the esterification of free cholesterol into cholesterol esters, which are stored in cytoplasmic lipid droplets and serve as intracellular cholesterol reserves [[273]85, [274]86]. The accumulation of cholesterol and its derivatives has been strongly linked to cancer and liver diseases [[275]87, [276]88, [277]89, [278]90]. Consistent with these findings, clinical tissue experiments revealed a notable upregulation of SOAT2 mRNA and m5C levels in HCC tissues versus adjacent normal tissues (Figure [279]6F, G). Further, SOAT2 was knocked down in HCC cells to confirm that NSUN2 directly mediated the m5C modification of SOAT2 mRNA and affected its transcription activity. RT‐qPCR analysis confirmed that sh‐SOAT2 #2 exhibited the optimal knockdown efficiency (Supplementary Figure [280]S11A), leading to its selection for subsequent experiments. The mRNA and protein expression levels were evaluated using WB assay and RT‐qPCR following NSUN2‐KO or SOAT2 knockdown. The data revealed a notable reduction in the expression of SOAT2 mRNA and protein in NSUN2‐KO HCC cells relative to the NSUN2‐WT group, whereas NSUN2 overexpression markedly increased the expression of SOAT2 versus the NC‐oe group (Figure [281]6H, I, Supplementary Figure [282]S11B, C). Notably, SOAT2 knockdown resulted in a substantial decrease in SOAT2 mRNA and protein expression in HCC cells, with no observable changes in NSUN2 mRNA and protein expression (Figure [283]6J, K, Supplementary Figure [284]S11D, E). In addition, NSUN2‐KO accelerated the degradation of SOAT2 mRNA, whereas NSUN2 overexpression prolonged the half‐life of SOAT2 mRNA (Figure [285]6L, Supplementary Figure [286]S11F). MeRIP‐qPCR experiments demonstrated that NSUN2‐KO diminished the m5C levels of SOAT2 mRNA in HCC cells, while NSUN2 overexpression reversed the results (Figure [287]6M, Supplementary Figure [288]S11G). Furthermore, metabolomics analysis was conducted to assess the cholesterol content in the supernatant of HCC Hepa1‐6 and H22 cells. The findings revealed a pronounced reduction in cholesterol content in NSUN2‐KO/sh‐SOAT2 cells relative to NSUN2‐WT/sh‐NC groups, while overexpression of NSUN2 markedly increased cholesterol content compared to the NC‐oe group (Figure [289]6N, O, Supplementary Figure [290]S11H, I). We next moved to investigate the role of cholesterol between tumor cells and CD8^+ T cells by quantifying the cholesterol levels in the cell culture medium. The results showed that cholesterol levels in the medium were markedly higher in the NSUN2‐oe group compared to the control group (P < 0.05) (Figure [291]6P), a trend consistent with the increase in intracellular cholesterol levels. These results suggest that cholesterol may serve a dual function in tumor cells, acting both as a metabolic substrate and as a signaling molecule that facilitates immune escape. In summary, our findings demonstrate that NSUN2 mediates m5C modification of SOAT2 in HCC, thereby enhancing its mRNA transcription and contributing to cholesterol‐mediated immune escape. 3.8. NSUN2‐mediated m5C modification of SOAT2 facilitated immune escape in HCC cells A recent study had highlighted that cholesterol within the TME drove the excessive uptake of fatty acids by CD8^+ T cells, leading to lipid oxidation damage and ferroptosis [[292]91]. This process impaired the cytotoxic function of CD8^+ T cells, thereby facilitating immune escape by cancer cells. Pathway enrichment analysis revealed that DEGs were intricately linked to reprogramming pathways such as lipid metabolism and glycolysis (Supplementary Figure [293]S6). Based on these findings, it was hypothesized that NSUN2 could modulate the lipid metabolism reprogramming of cancer cells through m5C modification of SOAT2, thereby promoting immune escape in HCC cells. To investigate this hypothesis, the mRNA and protein expression levels of genes were analyzed following the KO of NSUN2 and (or) overexpression of SOAT2 using WB assay and RT‐qPCR. The results showed that relative to the NSUN2‐WT + NC‐oe group, both NSUN2 and SOAT2 mRNA expression levels were markedly diminished in the NSUN2‐KO + NC‐oe group; compared to the NSUN2‐KO + NC‐oe group, both mRNA and protein expression levels of SOAT2 were markedly increased in the NSUN2‐KO + SOAT2‐oe group, while the expression level of NSUN2 remained unchanged (Figure [294]7A, B). NSUN2‐KO in H22 cells showed the similar results (Supplementary Figure [295]S12A, B). FIGURE 7. FIGURE 7 [296]Open in a new tab NSUN2 mediates m5C modification of SOAT2 to regulate energy metabolism reprogramming and promote immune escape. (A‐B) Expression levels of NSUN2 and SOAT2 in different groups of Hepa1‐6 cells detected by RT‐qPCR and Western blotting. (C) MMP of different groups of cells detected by JC‐1 staining. (D) Changes in mitochondrial morphology observed by fluorescence microscopy after Mito‐tracker staining. (E) Glucose uptake in different groups of cells. (F‐G) Lactate and ATP production in different groups of cells; (H) OCR in different groups of cells. (I) Cell viability measured by MTT assay. (J) Levels of GZMB and IFN‐γ in cell culture medium detected by ELISA. (K) Cholesterol levels in different groups of cells measured by metabolomics. (L) Cell viability measured by MTT assay. (M) Proliferation capacity of different groups of cells detected by EdU assay. (N‐O) Migration and invasion capabilities of different groups of cells measured by Transwell assay. (P) Apoptosis rate of different groups of cells detected by TUNEL assay. Cell experiments were repeated 3 times. Abbreviations: RT‐qPCR, reverse transcription quantitative polymerase chain reaction; KO: knock out; MTT: Methylthiazolyldiphenyl‐tetrazolium bromide; EdU: 5‐ethynyl‐2`‐deoxyuridine;TUNEL: TdT‐mediated dUTP Nick‐End Labeling; NSUN2, NOP2/Sun RNA methyltransferase family member 2; SOAT2, sterol O‐acyltransferase 2; ELISA: Enzyme linked immunosorbent assay. The MMP in cancer cells was assessed using JC‐1 dye. We observed a decrease in MMP in the NSUN2‐KO + NC‐oe group compared to the NSUN2‐WT + NC‐oe group, while an increase in MMP was observed in the NSUN2‐KO + SOAT2‐oe group compared to the NSUN2‐KO + NC‐oe group (Figure [297]7C, Supplementary Figure [298]S12C). Mitochondrial morphology alterations were visualized through Mito‐tracker fluorescence staining and microscopy. The results unveiled that compared to the NSUN2‐WT + NC‐oe group, the mitochondrial network in the NSUN2‐KO + NC‐oe group was fragmented and scattered, whereas the mitochondria were more complex and developed in the NSUN2‐KO + SOAT2‐oe group (Figure [299]7D, Supplementary Figure [300]S12D). Further investigation into glycolysis‐related parameters indicated a notable decrease in glucose consumption and a notable increase in lactate production, along with a reduction in ATP generation in the NSUN2‐KO + NC‐oe group compared to the NSUN2‐WT + NC‐oe group; conversely, the NSUN2‐KO + SOAT2‐oe group showed enhanced glucose consumption and diminished lactate output, alongside increased ATP production when compared to the NSUN2‐KO + NC‐oe group (Figure [301]7E–G, Supplementary Figure [302]S12E–G). Moreover, Seahorse analysis highlighted a pronounced decrease in the OCR in the NSUN2‐KO + NC‐oe group compared to the NSUN2‐WT + NC‐oe group, while a prominent increase in OCR was observed in the NSUN2‐KO + SOAT2‐oe group relative to the NSUN2‐KO + NC‐oe group (Figure [303]7H, Supplementary Figure [304]S12H). The viability of CD8^+ T cells in an In vitro co‐culture cell model was measured using MTT assay. It was observed that CD8^+ T cell viability was markedly decreased in the NSUN2‐KO + NC‐oe group compared to the NSUN2‐WT + NC‐oe group but was markedly increased in the NSUN2‐KO + SOAT2‐oe group compared to the NSUN2‐KO + NC‐oe group (Figure [305]7I, Supplementary Figure [306]S12I). ELISA results found that the levels of GZMB and IFN‐γ were markedly diminished in the cell supernatant of the NSUN2‐KO + NC‐oe group compared to the NSUN2‐WT + NC‐oe group but were markedly increased in the NSUN2‐KO + SOAT2‐oe group compared to the NSUN2‐KO + NC‐oe group (Figure [307]7J, Supplementary Figure [308]S12J). Metabolomic analysis was employed to measure cholesterol content within tumor cells. The results demonstrated a prominent decrease in cholesterol content in the NSUN2‐KO + NC‐oe group compared to the NSUN2‐WT + NC‐oe group; in contrast, the cholesterol content was markedly increased in the NSUN2‐KO + SOAT2‐oe group compared to the NSUN2‐KO + NC‐oe group (Figure [309]7K, Supplementary Figure [310]S12K). The biological functions of HCC cells were examined using the MTT, EdU, Transwell, and TUNEL assays. The results indicated that cell viability, proliferation, migration, and invasion abilities were markedly diminished in the NSUN2‐KO + NC‐oe group compared to the NSUN2‐WT + NC‐oe group, and cell apoptosis was markedly increased; in contrast, these malignant phenotypes were markedly increased in the NSUN2‐KO + SOAT2‐oe group compared to the NSUN2‐KO + NC‐oe group, and cell apoptosis was markedly restricted (Figure [311]7L, P, Supplementary Figure [312]S12L–P). These results suggest that NSUN2 regulates energy metabolism reprogramming through m5C modification of SOAT2, thereby facilitating immune escape in cancer cells. 3.9. NSUN2 promoted tumor growth by regulating SOAT2 m5C modification‐mediated immune escape in H22 cells The aforementioned In vitro cell experiments confirmed that NSUN2‐KO could facilitate the reprogramming of energy metabolism through the mediation of SOAT2 m5C modification. This process enabled cancer cells to evade immune surveillance, thereby promoting the initiation and progression of HCC. RT‐qPCR and WB assay analyses were performed to verify whether this mechanism affected the In vivo tumorigenic ability of HCC cells. The findings revealed a pronounced reduction in the expression levels of NSUN2 and SOAT2 in the tumor tissues of the NSUN2‐KO + NC‐oe group compared to the NSUN2‐WT + NC‐oe group; conversely, compared to the NSUN2‐KO + NC‐oe group, the NSUN2‐KO + SOAT2‐oe group exhibited a pronounced increase in the expression of SOAT2, while the expression of NSUN2 remained unchanged (Figure [313]8A, B, Supplementary Figure [314]S13A, B). After 2 weeks, tumor growth in the mouse model was monitored using bioluminescent imaging. The results demonstrated that tumor growth was suppressed in the NSUN2‐KO + NC‐oe group relative to the NSUN2‐WT + NC‐oe group, while the NSUN2‐KO + SOAT2‐oe group exhibited enhanced tumor growth (Figure [315]8C, Supplementary Figure [316]S13C). FIGURE 8. FIGURE 8 [317]Open in a new tab Influence of NSUN2 regulation on the expression of SOAT2 in HCC and its impact on tumor growth. (A) Gene expression levels of NSUN2 and SOAT2 in tumor tissues of mice (6 mice/group) were measured using RT‐qPCR; (B) Protein expression levels of NSUN2 and SOAT2 in tumor tissues of mice were detected using Western blotting. (C) Tumor growth was monitored at different time points using In vivo live imaging; (D) Cholesterol content in tumor tissues of mice was measured using metabolomics analysis; (E) Morphology of tumor tissues in different groups of mice; (F) Tumor growth status in different groups of mice; (G) Tumor tissue weights in different groups of mice; (H‐K) Protein expression levels of CD8α (H‐I) and Ki67 (J‐K) in tumor tissues of mice were detected using immunohistochemical staining; (L) Contents of GZMB and IFN‐γ in tumor tissues were measured using ELISA. Abbreviations: RT‐qPCR, reverse transcription quantitative polymerase chain reaction; KO: knock out; NSUN2, NOP2/Sun RNA methyltransferase family member 2; SOAT2, sterol O‐acyltransferase 2. Moreover, lipid metabolism metabolomics analysis was undertaken for tumor tissues from each group of mice, accompanied by assessments of tumor volume and weight. The data indicated that the NSUN2‐KO + NC‐oe group displayed marked decreases in cholesterol content, tumor volume, and weight compared to the NSUN2‐WT + NC‐oe group, whereas the NSUN2‐KO + SOAT2‐oe group exhibited markedly increased cholesterol content, tumor volume, and weight (Figure [318]8D–G, Supplementary Figure [319]S13D–G). IHC of mouse tumor tissues further showed that compared to the NSUN2‐WT + NC‐oe group, the NSUN2‐KO + NC‐oe group displayed a notable increase in CD8α positivity rate and a corresponding decrease in Ki67 protein positivity rate, whereas the NSUN2‐KO + SOAT2‐oe group exhibited a notable decrease in CD8α positivity rate and a reverse trend in the Ki67 protein positivity rate (Figure [320]8H–K, Supplementary Figure [321]S13H–K). ELISA findings revealed that compared to the NSUN2‐WT + NC‐oe group, the NSUN2‐KO + NC‐oe group showed a pronounced decrease in GZMB and IFN‐γ in tumor tissues; conversely, compared to the NSUN2‐KO + NC‐oe group, the NSUN2‐KO + SOAT2‐oe group showed a prominent increase in GZMB and IFN‐γ in tumor tissues (Figure [322]8L, Supplementary Figure [323]S13L). These findings collectively indicate that NSUN2 promotes immune escape in H22 cells by regulating SOAT2 m5C modification, ultimately promoting tumor growth In vivo. 4. DISCUSSION Given the critical need to understand the molecular mechanisms driving HCC, this study revealed that NSUN2 regulated SOAT2 via m5C modification, contributing to energy metabolism reprogramming and immune escape. Based on the scRNA‐seq, bulk RNA‐seq, and meta‐analysis data, we identified key DEGs and DEPs involved in HCC progression. The overexpression of NSUN2 in HCC cells was found to be associated with diminished activity and cytotoxicity of CD8^+ T cells, as evidenced by curbed levels of key cytokines such as IFN‐γ and GZMB. These observations demonstrate the significance of NSUN2 in altering the TME to favor immune escape. This study provided a detailed elucidation of the mechanisms underlying the effects of NSUN2 on modulating energy metabolism reprogramming via m5C modification of SOAT2. scRNA‐seq identified a notable increase in malignant cell populations in HCC tissue, and subsequent cell communication analysis revealed that tumor cells may augment cancer development by evading immune clearance, which is a critical mechanism in tumor development and metastasis [[324]92, [325]93, [326]94]. These findings corroborated previous studies on the role of NSUN2 in oncogenesis but also offered a more refined understanding of how HCC cells exploited metabolic pathways to escape immune detection [[327]95, [328]96]. The study also demonstrated that NSUN2‐mediated m5C modification of SOAT2 enhanced cholesterol biosynthesis, thereby contributing to the reprogramming of energy metabolism in HCC cells. This metabolic reprogramming supported the rapid proliferation and survival of tumor cells by providing essential resources [[329]97, [330]98], such as increased lipid reserves, which are critical for membrane synthesis and energy production [[331]99, [332]100]. Furthermore, the elevated levels of intracellular and extracellular cholesterol observed upon NSUN2 overexpression suggested a dual role for cholesterol in both cellular metabolism and modulation of the TME [[333]88, [334]101]. NSUN2 may contribute to the suppression of immune cell function by altering the TME, particularly by reducing the activity of CD8^+ T cells, which are crucial for anti‐tumor immunity [[335]16]. Our meta‐analysis highlighted significant differential expression of NSUN2 between HCC and adjacent normal tissues, with a particularly strong association observed in the Asian population. This population‐specific expression pattern may be linked to the etiology of HCC, where chronic hepatitis infections are more prevalent and are known contributors to the development of HCC [[336]102, [337]103, [338]104]. The combination of meta‐analysis and advanced methylation sequencing techniques established NSUN2 as a critical factor in the progression of HCC, which coincided with existing studies [[339]15, [340]105, [341]106, [342]107, [343]108, [344]109]. These results were consistent with existing literature on the regulatory role of NSUN2 across various tumor types, thereby underscoring its importance in tumorigenesis and highlighting its potential as a therapeutic target in HCC. Furthermore, our research confirmed that NSUN2 facilitated the m5C modification of SOAT2, thereby enhancing its stability and transcriptional activity. This epigenetic modification likely confered a survival advantage to HCC cells by modulating metabolic pathways that are critical for rapid cell proliferation and immune escape [[345]110, [346]111]. Based on the CRISPR/Cas9 technology, we constructed KO models for key genes and generated cell lines with either ablated or overexpressed NSUN2 levels. The findings from these models demonstrated that NSUN2 conferred a pivotal role in energy metabolism reprogramming by mediating the m5C modification of SOAT2, a process that is essential for maintaining the metabolic flexibility of cancer cells [[347]112, [348]113, [349]114]. These results align with previous studies and provide more direct and conclusive evidence of the involvement of NSUN2 and SOAT2 in the metabolic adaptations that support HCC progression. The identification of NSUN2 and SOAT2 as key players in HCC progression holds promising scientific and clinical significance. SOAT2 is a key enzyme involved in cholesterol metabolism [[350]115, [351]116, [352]117]. The increased expression and stability of SOAT2 in response to NSUN2 overexpression suggested that SOAT2 imparted critical impacts on the metabolic adaptations that enable HCC cells to thrive under hostile conditions, including immune surveillance [[353]118, [354]119]. These findings highlight SOAT2 as a crucial effector in the NSUN2‐driven metabolic pathway, contributing to the aggressive nature of HCC. This study demonstrated the potential of NSUN2 or SOAT2 as novel indicators for the early detection and prognostic evaluation of HCC. The identification of these molecules offered new avenues for therapeutic intervention, as targeting NSUN2 or SOAT2 could disrupt the metabolic reprogramming and immune escape mechanisms that are vital for tumor survival [[355]17, [356]120–[357]122]. The strategic inhibition of these targets may present innovative therapeutic strategies for HCC treatment, potentially improving patient outcomes. The potential impact of NSUN2 on the activity of EPCs in the TME was examined. EPC‐related genes significantly expressed in different cell populations, providing further insights into the role of EPCs in the NSUN2‐mediated TME. Hepatic‐derived POSTN induced elevated CCL2 expression and secretion in EPCs, and CCL2 promoted prometastatic traits in HCC [[358]123]. The immunomodulatory function of epithelial cells has emerged as a pivotal determinant of the immune response, endothelial cells can express immune checkpoint molecules, such as PD‐L1, which can modulate immune cell activity [[359]124]. In the future, we will further focus on impact of NSUN2 on the activity of EPCs in the TME. However, this study hads several limitations that must be acknowledged. A noticeable limitation is the reliance on In vitro models, which may not fully reflect the complexity of the TME In vivo. The use of Matrigel as a matrix component, while helpful in mimicking some aspects of the TME, fails to capture the diverse cell populations and intricate signaling networks present in a living organism. Consequently, the generalizability of our findings is restricted by the absence of human clinical data, which limits the applicability of these results to broader patient populations. Furthermore, this study is primarily observational, validating the functionality of NSUN2 and SOAT2 through cell and animal experiments, yet it remains unclear whether these molecules can reliably serve as prognostic indicators for HCC in human clinical settings. Another limitation of our research is the incomplete exploration of the potential downstream effects of NSUN2‐mediated m5C modifications. This leaves a gap in the understanding of the broader impact of these modifications on HCC progression. The study did not delve deeply into the signaling pathways and regulatory mechanisms associated with NSUN2 and SOAT2, which are critical for understanding their roles in HCC progression. Additionally, our research did not explore all potential downstream effects of NSUN2‐mediated m5C modifications, leaving a gap in the understanding of the broader impact of these modifications on HCC progression. Further investigations are necessary to pinpoint additional therapeutic targets and to elucidate the complex signaling pathways and regulatory mechanisms associated with NSUN2 and SOAT2. Such studies could provide a more thorough understanding of HCC biology and reveal novel therapeutic opportunities. Considering the findings of this study, several avenues for further exploration are suggested. More clinical research is needed to verify the role and clinical significance of NSUN2 in HCC, as well as its potential as a therapeutic target. Further investigations can explore the interaction of NSUN2 with other related molecules in the immune escape of HCC to uncover additional regulatory mechanisms and pinpoint new treatment targets. Moreover, potential treatment strategies, such as targeted drug design against NSUN2 and the combined application of immunotherapy and NSUN2‐targeted treatment, should be explored to enhance the effectiveness of HCC therapies. 5. CONCLUSIONS In conclusion, this study highlights the pivotal role of NSUN2 in the progression of HCC by regulating SOAT2 through m5C modification. A particularly notable finding is the ability of NSUN2 to attenuate the activity and cytotoxicity of CD8^+ T cells, thereby promoting immune escape in HCC cells. This enhanced understanding of the NSUN2‐SOAT2 interaction demonstrates its potential as a therapeutic target, offering new possibilities for disrupting these mechanisms to improve immune response and inhibit tumor progression in HCC. AUTHOR CONTRIBUTIONS JHJ, DC, FL, FG, JCC and THY conceived and designed the study. CXX, FG, JCC performed the experiments. THY analyzed the data. JHJ, FG, JCC and THY wrote the manuscript. All authors reviewed and approved the final version of the manuscript. CONFLICT OF INTEREST STATEMENT The authors declare that they have no competing interests. ETHICS APPROVAL AND CONSENT TO PARTICIPATE All experimental procedures and protocols for animal use were reviewed and approved by the Animal Ethics Committee of Renji Hospital, School of Medicine, Shanghai Jiaotong University. Informed consent was obtained from all patients, and the study was approved by the Ethics Committee of Renji Hospital, School of Medicine, Shanghai Jiaotong University, in compliance with the Helsinki Declaration. CONSENT FOR PUBLICATION Not applicable. Supporting information Supporting Information [360]CAC2-45-846-s001.docx^ (2.5MB, docx) ACKNOWLEDGMENTS