Comprehensive Bioinformatics Analysis and Experimental Verification RNF186 Is a Recurrence Signature Gene of Hepatocellular Carcinoma that Promotes Cell Proliferation

1  Introduction

Hepatocellular carcinoma (HCC) is the third leading cause of cancer-related deaths worldwide, accounting for 90% of primary liver cancer cases [1–3]. Despite advancements in treatment, including surgical resection and neoadjuvant therapy, approximately 70% of patients with resectable HCC experience relapse and metastasis, significantly limiting long-term survival [4–6]. Therefore, a deeper understanding of the molecular mechanisms underlying HCC progression and recurrence is essential for improving treatment strategies and developing personalized medicine.

Recurrence of HCC post-resection can be classified as early (within weeks or months) or late (more than two years post-resection) [7]. Early recurrence, which accounts for nearly 70% of all cases, is often attributed to micrometastases or incomplete tumor removal [8], whereas late recurrence is typically associated with newly formed primary tumors [4]. Ongoing clinical trials are exploring the efficacy of adjuvant immunotherapy to reduce recurrence risk after surgery [9]. However, molecular studies of HCC have largely focused on primary tumors, leaving the evolutionary relationship between primary and recurrent tumors poorly understood, with several biological and clinical questions remaining unanswered.

The emergence of next-generation sequencing (NGS) and single-cell sequencing technologies has led to the creation of extensive public tumor databases, enabling the integration and reanalysis of large datasets to uncover crucial insights into tumor pathogenesis [10,11]. Recently, prognostic models based on gene signatures have been increasingly adopted [12,13].

This study aimed to screen for prognostic features of recurrent hepatocellular carcinoma (HCC) by integrating multi-omics data and further strengthening the level of evidence through cellular function experiments.

2  Materials and Methods

2.1 Data Collection

The workflow is shown in Fig. 1. Single-cell transcriptome data from 12 primary hepatocellular carcinoma (HCC) tumors and 6 recurrent tumors were obtained from the study of Sun et al. [10]. 374 tumor samples and 50 adjacent normal samples were obtained from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/, accessed on 02 September 2024), supplemented with gene expression information and clinical information, including overall survival, survival status, patient age, sex, and TNM stage, from the GTEx database (https://gtexportal.org/home/aboutAdultGtex, accessed on 23 October 2025). Additionally, ICGC-LIRI-JP (https://dcc.icgc.org/, accessed on 23 October 2025) and GSE14520 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE14520, accessed on 23 October 2025) were collected as independent validation cohorts. Batch effects were corrected using ComBat (from the sva package) during data merging. Clinical annotations included overall survival, survival status, patient age, sex, and TNM stage.

Figure 1: Work model. Abb: PCA: principal components analysis; tSNE: t-distributed stochastic neighbor embedding; UMAP: Uniform Manifold Approximation and Projection; TCGA: The Cancer Genome Atlas; HCC: Hepatocellular carcinoma; Cox: Cox regression model; CCK-8: Cell Counting Kit-8; GRID2: Glutamate Ionotropic Receptor Delta Type Subunit 2; RNF186: Ring Finger Protein 186; SLC4A10: Solute Carrier Family 4 Member 10.

2.2 Single-Cell Sequencing Analysis

As the dataset had undergone prior preprocessing and quality control in the original publication, downstream analyses were performed directly using the Seurat (v5.0.2) R package [14]. Principal component analysis (PCA) followed by t-distributed stochastic neighbor embedding (t-SNE) was utilized for dimensionality reduction, and the Louvain algorithm was used for unsupervised clustering. Cellular identities were assigned with reference to the Human primary cell atlas (HPCA) in the SingleR (4.3.1) package [15].

2.3 Cell-Cell Interaction Network Analysis

After normalization and batch effect correction were performed using the NormalizeData function in the Seurat package, the CellChat v.1.6.1 R package was used to quantitatively infer and analyze intercellular communication networks based on scRNA-seq data, aiming to predict major cellular signal inputs and outputs and elucidate how cells and signaling coordinate various functions [16]. Signaling pathways were identified based on known ligand-receptor pairs (CellChatDB.human, https://github.com/sqjin/CellChat, accessed on 02 September 2024) and their expression in specific cell types. Adjust the key parameters of CellChat (such as the number of signal path groupings k value) to observe the stability of the results.

2.4 Risk Score Prognostic Model Construction

After standardizing and normalizing the gene expression data, the extracted differentially expressed genes associated with HCC recurrence were further screened using univariate Cox regression to enhance their prognostic relevance (p < 0.05). Subsequently, a multigene prognostic risk score was constructed using a minimum absolute contraction and selection operator (Lasso) regression model, and a nomogram was plotted by combining it with clinical parameters [17]. The reliability of the scores was assessed using calibration curves, Kaplan-Meier analysis, and the area under the time-dependent receiver operating characteristic (ROC) curve (AUC, 95% CI).

2.5 Functional Enrichment Analysis and Protein-Protein Interaction Network Construction

The functional enrichment analysis of genes was performed using the DAVID tool (version 6.8, https://davidbioinformatics.nih.gov/tools.jsp, accessed on 02 September 2024), with annotations for its terminology including Gene Ontology (GO) terminology and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp/kegg/, accessed on 02 September 2024) (FDR < 0.05, log2 |Fold Change (FC)| > 2) [18,19]. The protein-protein interaction (PPI) network was constructed by uploading the gene list to the STRING database (version 11.0, https://string-db.org/, accessed on 02 September 2024) [20].

2.6 Evaluation of Immune Status, Tumor Stemness, and Genomic Heterogeneity

Three deconvolution algorithms—CIBERSORT, EPIC, and ESTIMATE—were used to comprehensively characterize immune cell infiltration, with expression matrices normalized and batch-effect corrected prior to each. The results from all three methods were analyzed to assess the overall immune status of the sample, rather than being limited to any single method [21,22]. Data quality control in calculating microsatellite instability (MSI) and tumor mutational burden (TMB) scores included excluding samples with incomplete mutation annotation, sequencing artifacts, or mutation counts exceeding three standard deviations of the cohort mean, while missing values were directly removed [23]. The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/, accessed on 02 September 2024) was used to predict the immune escape potential and response to immune checkpoint inhibitors in each subgroup [24].

2.7 Protein Structure Prediction and Mutation Analysis

Protein structures of the genes were accessed from the GeneCards database (https://www.genecards.org/, accessed on 02 September 2024). 296 samples with detected mutations were obtained from cBioPortal (https://www.cbioportal.org/, accessed on 02 September 2024), of which 153 (51.7%) were used for localization. The chi-square test was used to assess mutation frequency and variant types, including single nucleotide polymorphisms (SNPs), oligonucleotide polymorphisms (ONPs), insertions, and deletions, between the high-expression and low-expression groups.

2.8 Cell Culture and Transfection

Huh7 (CC0102) and PLC/PRF/5 (CC0111) cell lines were purchased from Guangzhou Saiku Biotechnology Co., Ltd. (Guangzhou, China). Both cell lines were authenticated by short tandem repeat (STR) profiling and confirmed to be free of mycoplasma contamination prior to use. Cells were cultured in a humidified incubator at 37°C in a supplemental medium (DMEM, Cat. No. 12491015, Thermo Fisher Scientific, Shanghai, China) containing 10% fetal bovine serum (FBS, Cat. No. C2739, Beyotime, Shanghai, China) and two antibiotics (penicillin 100 U/mL, Cat. No. ST486-10g, Beyotime and streptomycin 100 U/mL, Cat. No. C0222, Beyotime). For gene overexpression, SESN2 (NM_031459.5) and RNF186 (NM_019062.2) were generated by cloning human cDNAs into the pcDNA3-Flag vector. For gene silencing experiments, targeted small interfering RNAs (siRNAs) (BGI, Shenzhen, China) were transfected into cells using Lipofectamine 2000 reagent (Cat. No. 11668500, Thermo Fisher Scientific), following the manufacturer’s instructions. The final concentration of siRNA in the medium was adjusted to 20 nm. The specific siRNA sequences utilized are detailed in Table 1.

2.9 Western Blot Analysis

Forty-eight hours after siRNA transfection, total protein was extracted for immunoblot analysis. Cells were lysed in buffer containing 1% sodium dodecyl sulfate (SDS, Cat. No. ST626, Beyotime), and the total protein concentration was quantified using the bicinchoninic acid (BCA) assay kit (Cat. No. A55860, Thermo Fisher Scientific) according to the manufacturer’s protocol.

Equal amounts of protein (18 μg) lysate were separated by SDS-PAGE (10%) and transferred onto PVDF membranes (Cat. No. FFP19, Beyotime), followed by blocking with a 5% solution of skim milk (Cat. No. P0216-300g, Beyotime) powder dissolved in TBST at room temperature for 2 h. The membranes were incubated with primary antibodies in 4°C overnight targeting GRID2 (Dilution ratio, 1:1000, Cat. No. ab198499, Abcam, Shanghai, China), RNF186 (Dilution ratio, 1:1000, Cat. No. H00054546-M01, Biotechne, Wuhan, China), SLC4A10 (Dilution ratio, 1:1000, Cat. No. 27197-1-AP, Proteintech, Wuhan, China), HA (Dilution ratio, 1:5000, Cat. No. ab236632, Abcam), Flag (Dilution ratio, 1:2000, Cat. No. 20543-1-AP, Proteintech), SESN2 (Dilution ratio, 1:1000, Cat. No. ab178518, Abcam), or β-Actin (Dilution ratio, 1:3000, Cat. No. 66009-1-Ig, Proteintech) as a loading control. Secondary antibody incubation was performed using goat anti-rabbit (Dilution ratio, 1:2000, Cat. No. 27197-1-AP, Beyotime) or goat anti-mouse (Dilution ratio, 1:3000, Cat. No. 27197-1-AP, Beyotime) and incubated for 2 h at room temperature. Detection was performed using enhanced chemiluminescence (ECL, Cat. No. P0018S, Beyotime) reagents (Smart-Lifesciences) to visualize protein bands. All experiments were repeated three times.

2.10 Ubiquitination Assays and Co-Immunoprecipitation

Cells were cultured in 6 cm culture dishes, and during cell passage, the seeding density was maintained at at least 30% (usually one culture dish was passaged to three culture dishes) for ubiquitination analysis and immunoprecipitation experiments. When cell growth reaches a density of 60%–80%, cells were transfected according to the indicated conditions for 48 h and lysed using guanidine denaturing buffer (6 M guanidine hydrochloride, 0.1 M Na2HPO4/NaH2PO4, 10 mm imidazole, pH 8.0). After sonication, whole cell extracts were incubated with 30 μl anti-Flag magnetic beads (Cat. No. P2115, Beyotime) and washed, and the supernatant was immunoprecipitated at 4°C for 4 h or overnight (ensuring total protein content is above 1 mg). The immune complexes were washed three times with native lysis buffer, and the pulled-down proteins were analyzed by immunoblotting. For immunoprecipitation, cells were lysed using RIPA buffer (Cat. No. P0013, Beyotime) (non-denaturing lysis buffer). All other steps were the same as for the ubiquitination assay. All experiments were repeated three times. All immunoblot images were subjected to exposure analysis using Image Lab (v.3.0, Bio-Rad, Shanghai, China).

2.11 Real-Time Quantitative Polymerase Chain Reaction (RT-qPCR)

Total RNA was isolated from Huh7 cells using TRIzol (Cat. No. R0016, Beyotime, Shanghai, China) reagent following the standard protocol and using Exonuclease I (1 μl, Cat. No. 70073Z2500UN, Thermo Fisher Scientific) to avoid genomic contamination. First-strand complementary DNA (cDNA) synthesis was carried out using the iScript cDNA Synthesis Kit (Cat. No. 1708890, Bio-Rad, Hercules, CA, USA) using 20 ng of total RNA per replicate well. Quantitative PCR amplification was performed with Fast SYBR Green Master Mix (Bio-Rad, Hercules, CA, USA) on a real-time PCR system, using a three-step cycle of 35–40 times (95°C (15 s) + 60°C (15 s) + 72°C (45 s)). All experiments were biologically replicated three times.

Actin served as the internal control gene for normalization. Relative gene expression levels were calculated using the comparative 2−ΔΔCt method. The primer sequences used in the analysis were as follows:

GRID2-F: 5′-TTTGTCCGTCTGGTGGTCTC-3′

GRID2-R: 5′-CAGTGCGAAATACCTCATCATCC-3′

RNF186-F: 5′-GGCTGGTGGATCCTGCTGAC-3′

RNF186-R: 5′-CTGCTCTCTGCAGAAGAGAG-3′

SLC4A10-F: 5′-GTCACAGGCATCGTGGTCATA-3′

SLC4A10-R: 5′-CCTTCACGCCAACAAATCTCAT-3′

SESN2-F: 5′-AGGACTACCTGCGGTTCGCC-3′

SESN2-R: 5′-GTAGTGCAGGCGCCAGAAGCTGG-3′

Actin-F: 5′-GACGGCCAGGTCATCACTATTG-3′

Actin-R: 5′-CCACAGGATTCCATACCCAAGA-3′

2.12 Colony Formation Assay

Approximately 1000 cells per well were plated into 6-well plates and cultured under standard conditions for a period of two weeks to allow colony development. Cells were fixed in 4% paraformaldehyde for 15 min at room temperature. Fixed cells were stained with 0.5% (w/v) crystal violet solution (Cat. No. C0775, Sigma-Aldrich, Shanghai, China) to visualize colonies. Images were acquired using the native camera of an iPhone 14 Pro. Clone counting was performed using ImageJ (v.1.53e, National Institutes of Health, Bethesda, MD, USA) software, following the official standard procedure.

2.13 Cell Viability Assay

Cell viability was assessed 48 h following siRNA transfection utilizing the Cell Counting Kit-8 (CCK-8, (Cat. No. C0037, Beyotime)). Briefly, cells were seeded in 96-well plates at a density of 2 × 104 cells per well and allowed to adhere for 12 h. Each well then received 100 μl of complete DMEM medium along with 10 μl of the CCK-8 reagent, followed by incubation for 1 h at 37°C. Absorbance was measured at 450 nm using a Multiskan™ FC microplate reader (Cat. No. 1410101, Thermo Scientific, Waltham, MA, USA). All viability measurements were performed in triplicate to ensure reproducibility. Absorbance standardization: Five technical replicates were performed for each sample well, and the average value was taken as the final result. Standardization calculation: Sample value/Maximum sample value in this group.

2.14 Statistical Analysis

Comparative analyses between groups were conducted using the Wilcoxon rank-sum test, while Spearman’s rank correlation was applied to examine associations between variables. All statistical tests were two-tailed, with significance defined as p < 0.05. Analyses were conducted using R software (version 4.2.2).

3  Results

3.1 Single-Cell Analysis of Hepatocellular Carcinoma

We obtained single-cell sequencing data from 18 hepatocellular carcinoma (HCC) samples, comprising 12 primary and 6 early-recurrent cases, as described by Sun et al. [10]. Using the “Seurat” package, we identified the top 3000 highly variable genes, which exhibited significant intercellular expression differences essential for cell-type identification (Fig. 2A). Principal component analysis (PCA) and t-SNE were employed for dimensionality reduction, clustering the samples into 18 distinct cell groups (Fig. 2B). Key differentially expressed genes were visualized via volcano plots (Fig. 2C). The top 10 heterogeneous genes were mapped in a t-SNE plot, showing cluster-specific expression, with APOA2 notably enriched (Fig. A1). Further annotation revealed seven major cell types: macrophages, T cells, NK cells, monocytes, endothelial cells, fibroblasts, and adipocytes (Fig. 2D,E). Comparison of t-SNE plots between primary and recurrent samples revealed minimal differences in overall cell distribution (Fig. 2F–H). However, recurrent samples displayed a notable reduction in macrophages and a significant increase in endothelial cells compared to primary samples (Fig. 2I).

Figure 2: Single-cell sequencing analysis identifies relapse-related genes. (A) Volcano plot illustrating the top 3000 highly variable genes. (B) Heatmap showing the expression of key differential genes across different clusters. (C) Volcano plot displaying the top 10 key genes with the greatest differences across different cell types. (D) Principal component analysis (PCA) and t-SNE dimensionality reduction separating cells into 18 distinct clusters. (E) Cell annotation using the SingleR package, categorizing cells into 7 major groups. (F) t-SNE 2D plots colored by primary or recurrent samples. (G,H) t-SNE 2D plots of 12 primary samples (G) or 6 recurrent samples (I). (I) Proportion of cell content in 12 primary samples or 6 recurrent samples.

3.2 Cell Communication Analysis and Identification of Relapse-Specific Genes

Cell Chat analysis revealed an increase in both the number and intensity of intercellular communications in recurrent samples relative to primary samples (Fig. 3A,B). Endothelial cells exhibited the highest number of communications, while macrophages showed the fewest (Fig. 3C). In terms of communication intensity, interactions between fibroblasts and NK cells were the strongest (Fig. 3C). Comparative analysis of communication-related genes between primary and recurrent samples identified 126 differentially expressed genes (Fig. 3D, Table A1), most of which were enriched in pathways such as cell adhesion molecules, cytokine-cytokine receptor interactions, and the PI3K-Akt pathway (Fig. 3E). Pathway enrichment analysis of the top 10 differentially expressed genes within each cell cluster revealed significant enrichment in cytokine-cytokine receptor interactions, natural killer cell-mediated cytotoxicity, and chemokine signaling pathways (Fig. 3F). Notably, chemoattractant cytokine ligand (CCL) signaling exhibited the most substantial changes in recurrent samples, particularly in endothelial cells compared to primary samples (Fig. A2).

Figure 3: Extraction of relapse-specific genes and cell-cell communication. (A,B) Intercellular communication between 7 cell types. The left side is the communication network, the thickness represents the number of communications (A) or strength (B), and the color represents the direction of communication. The right side is the quantitative graph. (C) The number and strength of intercellular communication are represented by a heat map. The color depth represents the relative value. The left side is the number of communications, and the right side is the communication strength. The observed cell types are highlighted in the red box. (D) Gene expression heat map showing 126 differentially expressed genes between primary and relapse samples. (E) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of 126 differentially expressed genes. (F) KEGG analysis of 180 recurrent signature genes.

3.3 Risk Score Construction Based on Relapse-Specific Genes Using the LASSO-Cox Algorithm

To identify relapse-specific genes associated with prognosis, we conducted univariate survival analysis on 180 differentially expressed genes, identifying 34 genes with significant associations (p < 0.05) (Fig. 4A). LASSO regression was used to select 9 key genes—CDCA2, RNF186, GTSE1, CDCA5, GRID2, AMIGO1, CERS1, ANKRD66, and SLC4A10—to construct a prognostic risk score model (Fig. 4B,C).

Figure 4: Development of a risk score model using the Least absolute shrinkage and selection operator (LASSO) algorithm and integration with clinical parameters. (A) Univariate Cox analysis identifying 34 prognostic genes. (B,C) Gene selection using the LASSO-Cox algorithm for constructing the optimal model. (D) Expression patterns and survival time relationships of 9 risk score genes. (E) Kaplan-Meier analysis demonstrating poorer prognosis in the high-risk group compared to the low-risk group. (F) Receiver operating characteristic curve (ROC) curve analysis showing Area Under Curve (AUC) values for the risk score model at 1, 3, and 5 years. (G) Nomogram integrating risk scores and clinical factors to predict 1-, 3-, and 5-year survival rates. (H) Calibration curve demonstrating the model’s accuracy in predicting patient survival rates. (I) Kaplan-Meier analysis indicating worse prognosis for patients with higher nomogram scores. (J) ROC curve showing AUC values for the nomogram scores at 1, 3, and 5 years. (K) Decision curve analysis to verify model effectiveness.

RiskScore = 0.0644CDCA2 + 0.0625RNF186 + 0.0035GTSE1 + 0.0001CDCA5 + 1.0301GRID2 + 0.0093AMIGO1 + 0.0148CERS1 + 1.5702ANKRD66 − 0.1633SLC4A10

Heatmap analysis indicated consistent expression patterns across most selected genes, except for CDCA5 (Fig. 4D). Kaplan-Meier survival analysis showed that the high-risk group had a significantly lower mortality rate than the low-risk group (Fig. 4E). The prognostic accuracy of the model was assessed using ROC curves, yielding AUC values of 0.76, 0.70, and 0.71 for 1-year, 3-year, and 5-year survival, respectively (Fig. 4F). A nomogram combining clinical parameters and the risk score was developed to predict 1-year, 3-year, and 5-year survival outcomes (Fig. 4G). Calibration curves confirmed the robustness of the model’s predictive capabilities (Fig. 4H). Kaplan-Meier analysis demonstrated that patients with higher nomogram scores had worse prognoses (Fig. 4I), and further ROC analysis yielded AUC values of 0.93, 0.96, and 0.94 for 1-year, 3-year, and 5-year survival, respectively (Fig. 4J). The DCA decision curve showed that the prediction of 1-year survival was good, but the performance was distorted for 3- and 5-year survival, which may be due to the low proportion of patients with survival of more than 3 and 5 years in the entire cohort (Fig. 4K).

3.4 Immune Landscape, Tumor Stemness, and Genomic Heterogeneity of Risk Scores

To assess the association between the risk score and immune status, we divided the TCGA-LIHC cohort into high-risk and low-risk groups. The CIBERSORT algorithm, analyzing 22 immune cell populations, revealed differences in immune cell infiltration between the groups, particularly in T cells, macrophages, and monocytes (Fig. 5A). Similarly, the EPIC algorithm showed higher macrophage infiltration in the low-risk group (Fig. 5B), while the Estimate algorithm also indicated increased immune infiltration in low-risk patients (Fig. 5C). Additionally, we extracted three categories of immune pathway genes (chemokines, receptors, and MHC) from the TCGA-LIHC cohort for comparison based on gene communication and signaling alterations. Most immune genes exhibited differential expression between the two groups. Specifically, chemokine-related genes were predominantly upregulated in the high-risk group, with a few exceptions showing downregulation (Fig. 5D), while receptor and MHC-related genes were consistently upregulated (Fig. 5E,F). Furthermore, the TIDE algorithm predicted higher risk scores in patients expected to respond to immunotherapy (Fig. 5G). High-risk patients exhibited lower scores for IFNG, Merck18, CD8, CD274, TIDE, Dysfunction, and TAM M2, but higher Exclusion and MDSC scores (Fig. 5H,I). Regarding tumor stemness, two algorithms produced conflicting results for risk groups, while other algorithms showed no significant differences (Fig. 5J). Genomic heterogeneity analysis demonstrated higher MSI and LOH scores in the high-risk group, suggesting a poor prognosis (Fig. 5K).

Figure 5: Analysis of immune status, tumor stemness, and genomic heterogeneity under risk score grouping. (A–C) Cibersort, Estimation of Proportions of Immune and Cancer cells (EPIC), and Estimate algorithms respectively calculated 22 immune cell infiltration scores (A), 8 immune cell infiltration scores (B), and overall immune scores (C) of high-risk and low-risk groups in the TCGA cohort samples and performed differential tests. (D–F) Differential expression of genes in three immune pathways between high-risk and low-risk groups, including chemokine (D), receptor (E), and major histocompatibility complex (MHC) (F). (G) Tumor Immune Dysfunction and Exclusion (TIDE) algorithm evaluates the difference in risk scores between the expected immunotherapy response group and the non-response group. (H) TIDE algorithm evaluates Interferon-γ (IFNG), Merck 18, CD8 Subunit Alpha (CD8), and CD274 scores in high-risk and low-risk groups. (I) TIDE algorithm evaluates TIDE, Dysfunction, Exclusion, Myeloid-derived suppressor cells (MDSC), cancer associated fibroblasts (CAF), and Tumor-Associated Macrophages (TAM) M2 scores in high-risk and low-risk groups. (J) Six scoring methods evaluate tumor stemness in high-risk and low-risk groups. (K) Seven scoring methods evaluate tumor genomic heterogeneity in high-risk and low-risk groups. *indicates p < 0.05, **indicates p < 0.01, ***indicates p < 0.001 and ****indicates p < 0.0001, -indicates no significance.

3.5 Expression and Prognostic Analysis of Risk Score Genes

To elucidate the role of the 9 risk score genes in hepatocellular carcinoma (HCC), we performed Kaplan-Meier (KM) survival analysis to evaluate their association with patient outcomes. In the TCGA cohort, the results showed that most genes alone did not predict prognosis, with the exception of AMIGO1 (Fig. 6A). In GSE14520, the transcript levels of CERS1, GTSE1, and SLC4A10 could distinguish between different lineages, while GRID2 and RNF186 did not (Fig. A3A). In the ICGC-LIRI-JP cohort, the transcript levels of CDCA2, CDCA5, and GTSE1 were closely associated with patient prognosis, while AMIGO1, CERS1, RNF186, and SLC4A10 were not significant (Fig. A3B). Multivariate COX regression can be used to evaluate the simultaneous contribution of all included variables, while controlling for covariates was evaluated to more fully assess the prognostic role of each gene in the risk score model. Inclusion of the nine genes in multivariate COX regression analysis identified GRID2, RNF186, and SLC4A10 as independent prognostic factors (Fig. 6B). Their inclusion in a prognostic model effectively distinguished different prognostic subgroups (Fig. 6C), with ROC curve areas under the curve (AUC) for 1-, 3-, and 5-year survival predictions at 0.77, 0.71, and 0.72, respectively (Fig. 6D). A single-cell t-SNE expression heatmap revealed distinct expression patterns of the 9 genes, with SLC4A10 showing higher expression levels, while CDCA2 and CDCA5 exhibited similar expression patterns (Fig. 6E).

Figure 6: Prognostic significance and single-cell analysis of the 9 risk score genes. (A) Kaplan-Meier survival analysis of the 9 risk score genes. (B) Multivariate Cox regression analysis of the 9 risk score genes. (C) Kaplan-Meier survival analysis of the multivariate Cox regression risk model. (D) ROC curve of the multivariate Cox regression risk model. (E) Expression distribution of the 9 risk score genes in the single-cell t-SNE plot.

3.6 Inhibition of RNF186, but Not GRID2 or SLC4A10, Suppresses Huh7 Cell Proliferation

Using Alphafold AI, we predicted the spatial conformations of GRID2, RNF186, and SLC4A10, which corresponded to proteins with molecular weights of 113 kDa, 24 kDa, and 125 kDa, respectively (Fig. A4A–C). In the TCGA cohort, RNF186 and SLC2A10 were lower expressed in tumor tissues, while GRID2 was nonsense (Fig. A4D–F). We used siRNA to suppress the expression of GRID2, RNF186, and SLC4A10, which were confirmed by real-time quantitative PCR and immunoblotting (Fig. 7A–F). Subsequent assays evaluated the proliferative capacity of Huh7 and PLC/PRF/5 hepatocellular carcinoma cells. CCK-8 and colony formation assays showed that GRID2 inhibition did not affect Huh7 and PLC/PRF/5 proliferation (Fig. 7G–J), while RNF186 inhibition significantly reduced cell proliferation (Fig. 7K–N). The results of the colony formation assay showed that si-SLC4A10#1 reduced the viability of PLC/PRF/5 cells, although this finding was not supported by Huh7 cells (Fig. 7O–R).

Figure 7: Functions of glutamate ionotropic receptor delta type subunit 2 (GRID2), ring finger protein 186 (RNF186), and solute carrier family 4 member 10 (SLC4A10) in hepatocellular carcinoma (HCC). (A–C) Real-time fluorescence quantitative polymerase chain reaction (RT-qPCR) detects the efficiency of si-RNA, including GRID2 (A), RNF186 (B), and SLC4A10 (C). (D–F) Western immunoblotting was used to detect the protein knockdown efficiency of si-RNA, including GRID2 (D), RNF186 (E), and SLC4A10 (F). (G,H) Cell counting kit-8 (CCK-8) (G) and colony formation assay (H) to detect the effect of si-GRID2 on the proliferation activity of Huh7 cells. (I,J) CCK-8 (I) and colony formation assay (J) to detect the effect of si-GRID2 on the proliferation activity of PLC/PRF/5 cells. (K,L) CCK-8 (K) and colony formation assay (L) to detect the effect of si-RNF186 on the proliferation activity of Huh7 cells. (M,N) CCK-8 (M) and colony formation assay (N) to detect the effect of si-RNF186 on the proliferation activity of PLC/PRF/5 cells. (O,P) CCK-8 (O) and colony formation assay (P) to detect the effect of si-SLC4A10 on the proliferation activity of Huh7 cells. (Q,R) CCK-8 (Q) and colony formation assay (R) to detect the effect of si-SLC4A10 on the proliferation activity of PLC/PRF/5 cells. ns indicates no significance, *indicates p < 0.05, **indicates p < 0.01, ***indicates p < 0.001.

3.7 RNF186 Promotes Huh7 Cell Proliferation via SESN2

Immunohistochemistry data from the HPA database showed that there was no significant difference in RNF186 in tumor tissues compared with normal liver tissues (Fig. A5A). STRING analysis identified interactions between RNF186 and the JAK family, inflammatory bowel disease-related proteins, NOD-like receptor signaling proteins, and components of cytokine-cytokine receptor interactions (Fig. A5B,C). Mutation profiling showed that low expression of RNF186 was associated with higher mutation frequencies of genes like CTNNB1, BIRC6, and RGPD4, while high expression correlated with increased mutation frequencies of MUC16, PCDH9, and SORCS3 (Fig. 8A). Furthermore, correlation analysis suggested that RNF186 was not correlated with various immune checkpoint genes at the transcriptional level (Fig. A5D), but the transcriptional level of RNF186 was higher in the immunotherapy response group evaluated by TIDE (Fig. 8B). In addition, in the sample immunotherapy scores evaluated by the “TIDE” algorithm, Dysfunction, CAF, and TIDE scores were closely related to the expression of RNF186 (Fig. A5E,F). Not only that, but RNF186 expression was also closely related to the stage of HCC and whether it was primary (Fig. 8C,D), although it was not related to tumor lymph node metastasis, whether it was distant metastasis, and tumor grade (Fig. A5G–I).

Figure 8: RNF186 promotes Huh7 cell proliferation in a manner dependent on SESN2. (A) Mutation analysis of the TCGA-LIHC cohort grouped by RNF186 transcript levels. (B) Evaluation of the relationship between RNF186 expression and immune response predicted by the TIDE algorithm. (C) Evaluation of the relationship between RNF186 expression and HCC clinical stage. (D) Evaluation of the relationship between RNF186 expression and whether HCC is primary or not. (E) KEGG enrichment analysis of the interacting protein profile of RNF186. (F) Huh7 cells were transfected with plasmids as indicated in the figure, lysed using the modified lysis buffer (Guan-HCl), and incubated with anti-Flag magnetic beads for 4–8 h at 4°C. Western blotting was performed using the antibodies indicated in the figure. (G) Huh7 cells were transfected with plasmids as indicated in the figure, lysed using the non-modified lysis buffer (RIPA), and incubated with anti-Flag magnetic beads for 4–8 h at 4°C. Western blotting was performed using the antibodies indicated in the figure. (H) Western blot analysis of Huh7 cells with or without overexpression of RNF186 using Flag or β-Actin antibodies. (I) Real-time fluorescence quantitative PCR to detect the inhibitory efficiency of si-SESN2. (J) Western immunoblotting detects the protein knockdown efficiency of si-SESN2. (K,L) CCK-8 (H) and colony formation assay (I) to detect the effect of si-SESN2 on the proliferation activity of Huh7 cells. (M,N) CCK-8 (J) and colony formation assay (K) to detect the effect of si-SESN2 on the proliferation activity of PLC/PRF/5 cells. (O,P) CCK-8 (L) and colony formation assay (M) were used to detect the effects of si-SESN2 transfection on the proliferation activity of Huh7 cells with or without overexpression of RNF186. (Q,R) CCK-8 (N) and colony formation assay (O) were used to detect the effects of si-SESN2 transfection on the proliferation activity of PLC/PRF/5 cells with or without overexpression of RNF186. ns and -indicate no significance, *indicates p < 0.05, **indicates p < 0.01, ***indicates p < 0.001.

To investigate potential downstream targets of RNF186 in cell proliferation regulation, we used BioGrid to identify interacting proteins. KEGG pathway enrichment analysis of these interacting proteins revealed significant enrichment in proteasomal degradation and ubiquitination pathways. Notably, among the enriched pathways, the crossover protein between the p53 and mTOR signaling pathways is SESN2 (Fig. 8E), suggesting that RNF186 may regulate cell proliferation in an SESN2-dependent manner [25,26]. To test this hypothesis, ubiquitination experiments in Huh7 cells showed that inhibition of RNF186 reduced the ubiquitination modification of SESN2 (Fig. 8F). Furthermore, exogenous SESN2 protein successfully co-immunoprecipitated endogenous RNF186 from Huh7 cells (Fig. 8G). We modulated SESN2 and RNF186 expression using siRNA and Flag-RNF186 plasmids, respectively, with RT-qPCR and immunoblotting confirming the efficacy of these interventions (Fig. 8H–J). Cell proliferation was assessed in Huh7 and PCL/PRF/5 cells using CCK-8 and colony formation assays. SESN2 knockdown significantly enhanced the proliferation of both cell lines (Fig. 8K–N). Furthermore, under SESN2 knockdown conditions, Flag-RNF186 overexpression failed to promote cell proliferation in either Huh7 or PCL/PRF/5 cells (Fig. 8O–R), indicating that SESN2 suppression effectively abrogates the proliferative effects of RNF186 in Huh7 cells.

4  Discussion

In this study, we conducted an in-depth analysis of single-cell sequencing data from liver cancer samples, focusing on the identification and functional characterization of genes associated with recurrence. We selected high-quality cells, clustered them into 18 distinct populations, and annotated seven major cell types. Notably, endothelial cells and macrophages emerged as key cell clusters linked to recurrence. We further developed a prognostic model based on macrophage-specific genes using the LASSO-Cox algorithm, demonstrating significant correlations between risk scores, immune landscapes, and patient outcomes. Multivariate Cox regression and functional assays identified GRID2, RNF186, and SLC4A10 as pivotal genes, with RNF186 shown to promote Huh7 cell proliferation in an SESN2-dependent manner.

While previous studies have primarily focused on various immune subsets, our work expands this knowledge by identifying specific cell clusters and prognostic genes associated with recurrence in liver cancer [11,27,28]. Macrophages and endothelial cells were identified as hallmark clusters in recurrent patients. Indeed, macrophages play a crucial role in liver cancer progression and recurrence [29–31]. Specifically, a reduction in macrophage populations in recurrent samples could reflect a diminished immune surveillance capacity, favoring tumor immune evasion [32,33]. It is particularly important to note that the subtype composition of macrophages may also shift. M1-like macrophages, typically associated with anti-tumor activity, may decrease, while M2-like macrophages, which promote tumor progression through immunosuppressive cytokine secretion and tissue remodeling, may be enriched [34].

Endothelial cell dysfunction, on the other hand, is a key driver of tumor angiogenesis and vascular remodeling, which are critical processes for sustaining tumor growth and metastasis [11,35]. An increase in endothelial cells in recurrent samples may reflect enhanced angiogenic activity, contributing to a pro-recurrence microenvironment [30]. Moreover, endothelial cells also influence the tumor immune microenvironment by modulating immune cell trafficking and response to immunotherapy [36]. Their enrichment in recurrent samples underscores their role in shaping tumor progression and resistance to therapeutic interventions [30,37]. Additionally, chemokine-mediated cell communication was notably altered in these clusters in recurrent samples, suggesting immune microenvironment changes in recurrent tumors. Although chemokines are implicated in tumor progression and immune evasion [37,38], their role in HCC recurrence remains poorly understood [39].

Multiple immune deconvolution algorithms were used to ensure the robustness of immune infiltration assessment. Although CIBERSORT, EPIC, and ESTIMATE rely on distinct reference matrices and computational principles, the major immune infiltration patterns identified by these methods were largely consistent. Minor discrepancies likely reflect methodological differences rather than biological contradictions. Considering that each algorithm captures unique aspects of the tumor microenvironment, we adopted a complementary strategy to interpret convergent trends rather than combine absolute values. The findings reveal a complex immune and genomic landscape in high-risk patients. Despite lower scores for key immunotherapy-related markers such as IFNG, CD8, CD274, and TIDE, these patients demonstrate a paradoxical susceptibility to immunotherapy. This discrepancy may be attributed to the interplay between tumor microenvironment features and genomic heterogeneity [40,41]. High-risk patients showed elevated Exclusion and MDSC scores, indicating an immune-suppressive microenvironment, which could be offset by other factors such as increased tumor mutational burden (TMB), microsatellite instability (MSI), or loss of heterozygosity (LOH) [42,43]. These genomic alterations may enhance the presentation of neoantigens, thereby improving immune recognition and response to immune checkpoint blockade, despite low baseline expression of classical markers like CD274 [44,45]. Enhanced tumor stemness can lead to enhanced tumor immune evasion [46,47]. In our analysis, however, different tumor stemness algorithms produced inconsistent results, with four out of six algorithms showing no significant differences. These discrepancies may arise from the distinct methodological principles of each algorithm. For instance, the RNAss and EREG.EXPss scores are calculated based on different molecular profiles, including RNA expression and methylation signatures, which may capture different aspects of tumor biology. The conflicting results from tumor stemness algorithms also highlight the complexity of the tumor’s biological behavior and its interaction with the immune system. Further exploration of these mechanisms, particularly the relationship between genomic instability and immunotherapy response, is warranted to clarify these findings and optimize therapeutic strategies for high-risk patients.

RNF186, an E3 ubiquitin ligase, is implicated in various diseases, such as ulcerative colitis and bladder cancer, by modulating substrate ubiquitination [48,49]. This study identified RNF186 as a signature gene for HCC recurrence and confirmed its pro-proliferative effect through cell function experiments, which is dependent on the ubiquitination regulation of SESN2. Regarding HCC recurrence, we tend to support the view that RNF186 maintains the proliferative capacity of HCC cells, which supports the recurrence of HCC, because in other diseases, genomic mutations in RNF186 itself often lead to the onset of the disease [48,50]. Conversely, there is also evidence that RNF186 is regulated by other molecules to promote cancer progression [49]. Human macrophages with the RNF186-A64T variant exhibit impaired NOD2 induction, reducing innate immune responses [51], suggesting that RNF186 dysregulation may affect macrophage function and contribute to HCC recurrence. Low RNF186 expression correlates with higher mutation frequencies in CTNNB1, BIRC6, and RGPD4, suggesting its potential role in suppressing Wnt/β-catenin signaling, maintaining apoptosis regulation, and preserving genomic stability [52,53]. Conversely, high RNF186 expression associates with increased mutations in MUC16, PCDH9, and SORCS3, which may indicate distinct tumor subtypes with unique molecular dependencies, such as immune evasion, disrupted tumor suppressor functions, and altered neurotrophin signaling [54,55]. These findings highlight the dual role of RNF186 in influencing HCC recurrence and immunotherapy relevance, underscoring its potential as a biomarker or therapeutic target for HCC.

Despite these significant findings, several limitations must be acknowledged. For example, although 18 samples (12 primary samples and 6 recurrent samples) were included in this study, the sample number was still small for hepatocellular carcinoma, which may lead to abnormal results that are caused by the heterogeneity of hepatocellular carcinoma. The recurrence prediction model in this study is based on scRNA-seq data, but the high cost and complex analysis limit its clinical application. To improve feasibility, key genes can be validated by more economical detection methods (such as qPCR or immunohistochemistry) in the future, and the stability of the model can be further verified in independent cohorts. In addition, combining other clinical parameters to optimize the model will help to achieve its wide application in personalized treatment, and its clinical utility requires evaluation through prospective studies. Moreover, additional in vivo and in vitro mechanistic studies are necessary to elucidate the precise roles of GRID2, RNF186, and SLC4A10 in liver cancer.

5  Conclusion

This study provides a comprehensive analysis of recurrence-related cell clusters and genes in liver cancer, revealing their prognostic significance and therapeutic potential. Functional experiments demonstrate that RNF186 promotes cell proliferation in an SESN2-dependent manner, laying a foundation for future research into the tumor microenvironment and the discovery of novel therapeutic targets for liver cancer.

Acknowledgement: Not applicable.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: Shanbao Ke: Conceptualization, Data curation, Software, Writing—original draft. Junya Yan: Data curation, Methodology, Writing—original draft. Xiao Feng: Data curation, Writing—original draft. Baiyu Li: Conceptualization, Project administration, Funding acquisition, Writing—review & editing. All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The datasets presented in this study can be found in online repositories, which deposited into the CNGB Sequence Archive (CNSA) of the China National GeneBank DataBase (CNGBdb, https://db.cngb.org/cnsa/) with accession number CNP0000650. Bulk RNA data comes from the TCGA project. All public data are fully described in the methods section.

Ethic Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

Appendix A

Figure A1: Expression distribution of the top 10 heterogeneous genes in single-cell t-SNE plot.

Figure A2: Heat map showing the differences in cell communication between chemoattractant cytokine ligand (CCL) signaling in primary and relapse samples.

Figure A3: The prognostic efficacy of the risk score genes was validated in an external cohort. (A) Kaplan-Meier survival analysis of risk score genes in the external cohort GSE14520. (B) Kaplan-Meier survival analysis of risk score genes in the external cohort ICGC-LIRI-JP.

Figure A4: Functions of GRID2, RNF186, and SLC4A10 in hepatocellular carcinoma (HCC). (A) GRID2 protein structure predicted by alphafold. (B) RNF186 protein structure predicted by alphafold. (C) SLC4A10 protein structure predicted by alphafold. (D) Difference in GRID2 transcription levels between normal and tumor groups in the TCGA-LIHC cohort. (E) Difference in RNF186 transcription levels between normal and tumor groups in the TCGA-LIHC cohort. (F) Difference in SLC4A10 transcription levels between normal and tumor groups in the TCGA-LIHC cohort. **p < 0.01, ****p < 0.0001, -no significance.

Figure A5: RNF186 promotes Huh7 cell proliferation in a manner dependent on SESN2. (A) Immunohistochemical detection of RNF186 expression in normal and HCC tissues in the HPA database. (B) Protein interaction network of RNF186 in the STRING database. (C) KEGG analysis of 40 RNF186 interacting proteins. (D) Correlation analysis between RNF186 expression and immune checkpoint genes. (E) Evaluation of RNF186 expression and tumor grade. (F,G) Differential analysis of immunotherapy scores calculated by the “TIDE” algorithm grouped by RNF186 expression. (H) Evaluation of the relationship between RNF186 expression and tumor lymph node metastasis. (I) Evaluation of the relationship between RNF186 expression and distant tumor metastasis. *p < 0.05, -no significance.

References

1. Bray F, Laversanne M, Sung H, Ferlay J, Siegel RL, Soerjomataram I, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2024;74(3):229–63. doi:10.3322/caac.21834. [Google Scholar] [PubMed] [CrossRef]

2. Wang Y, Deng B. Hepatocellular carcinoma: molecular mechanism, targeted therapy, and biomarkers. Cancer Metastasis Rev. 2023;42(3):629–52. doi:10.1007/s10555-023-10084-4. [Google Scholar] [PubMed] [CrossRef]

3. Llovet JM, Kelley RK, Villanueva A, Singal AG, Pikarsky E, Roayaie S, et al. Hepatocellular carcinoma. Nat Rev Dis Primers. 2021;7:6. doi:10.1038/s41572-020-00240-3. [Google Scholar] [PubMed] [CrossRef]

4. Wu J, Chan YT, Lu Y, Wang N, Feng Y. The tumor microenvironment in the postsurgical liver: mechanisms and potential targets of postoperative recurrence in human hepatocellular carcinoma. Med Res Rev. 2023;43(6):1946–73. doi:10.1002/med.21967. [Google Scholar] [PubMed] [CrossRef]

5. Tabrizian P, Abdelrahim M, Schwartz M. Immunotherapy and transplantation for hepatocellular carcinoma. J Hepatol. 2024;80(5):822–5. doi:10.1016/j.jhep.2024.01.011. [Google Scholar] [PubMed] [CrossRef]

6. Llovet JM, Pinyol R, Yarchoan M, Singal AG, Marron TU, Schwartz M, et al. Adjuvant and neoadjuvant immunotherapies in hepatocellular carcinoma. Nat Rev Clin Oncol. 2024;21(4):294–311. doi:10.1038/s41571-024-00868-0. [Google Scholar] [PubMed] [CrossRef]

7. Sherman M. Recurrence of hepatocellular carcinoma. N Engl J Med. 2008;359(19):2045–7. doi:10.1056/nejme0807581. [Google Scholar] [PubMed] [CrossRef]

8. Huang X, Song C, Zhang J, Zhu L, Tang H. Circular RNAs in breast cancer diagnosis, treatment and prognosis. Oncol Res. 2023;32(2):241–9. doi:10.32604/or.2023.046582. [Google Scholar] [PubMed] [CrossRef]

9. Kudo M. A paradigm change in the treatment strategy for hepatocellular carcinoma. Liver Cancer. 2020;9(4):367–77. doi:10.1159/000507934. [Google Scholar] [PubMed] [CrossRef]

10. Sun Y, Wu L, Zhong Y, Zhou K, Hou Y, Wang Z, et al. Single-cell landscape of the ecosystem in early-relapse hepatocellular carcinoma. Cell. 2021;184(2):404–21.e16. doi:10.1016/j.cell.2020.11.041. [Google Scholar] [PubMed] [CrossRef]

11. Li K, Zhang R, Wen F, Zhao Y, Meng F, Li Q, et al. Single-cell dissection of the multicellular ecosystem and molecular features underlying microvascular invasion in HCC. Hepatology. 2024;79(6):1293–309. doi:10.1097/HEP.0000000000000673. [Google Scholar] [PubMed] [CrossRef]

12. Li D, Shi Z, Liu X, Jin S, Chen P, Zhang Y, et al. Identification and development of a novel risk model based on cuproptosis-associated RNA methylation regulators for predicting prognosis and characterizing immune status in hepatocellular carcinoma. Hepatol Int. 2023;17(1):112–30. doi:10.1007/s12072-022-10460-2. [Google Scholar] [PubMed] [CrossRef]

13. Zhang X, Xiao Z, Zhang X, Li N, Sun T, Zhang J, et al. Signature construction and molecular subtype identification based on liver-specific genes for prediction of prognosis, immune activity, and anti-cancer drug sensitivity in hepatocellular carcinoma. Cancer Cell Int. 2024;24(1):78. doi:10.1186/s12935-024-03242-3. [Google Scholar] [PubMed] [CrossRef]

14. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36(5):411–20. doi:10.1038/nbt.4096. [Google Scholar] [PubMed] [CrossRef]

15. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. 2019;20(2):163–72. doi:10.1038/s41590-018-0276-y. [Google Scholar] [PubMed] [CrossRef]

16. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. 2021;12(1):1088. doi:10.1038/s41467-021-21246-9. [Google Scholar] [PubMed] [CrossRef]

17. Alhamzawi R, Ali HTM. The Bayesian adaptive lasso regression. Math Biosci. 2018;303:75–82. doi:10.1016/j.mbs.2018.06.004. [Google Scholar] [PubMed] [CrossRef]

18. Dennis G, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, et al. DAVID: database for annotation, visualization, and integrated discovery. Genome Biol. 2003;4(9):R60. doi:10.1186/gb-2003-4-9-r60. [Google Scholar] [CrossRef]

19. Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587–92. doi:10.1093/nar/gkac963. [Google Scholar] [PubMed] [CrossRef]

20. Szklarczyk D, Kirsch R, Koutrouli M, Nastou K, Mehryary F, Hachilif R, et al. The STRING database in 2023: protein-protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res. 2023;51(D1):D638–46. doi:10.1093/nar/gkac1000. [Google Scholar] [PubMed] [CrossRef]

21. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. 2018;1711:243–59. doi:10.1007/978-1-4939-7493-1_12. [Google Scholar] [PubMed] [CrossRef]

22. Racle J, Gfeller D. EPIC a tool to estimate the proportions of different cell types from bulk gene expression data. Methods Mol Biol. 2020;2120:233–48. doi:10.1007/978-1-0716-0327-7_17. [Google Scholar] [PubMed] [CrossRef]

23. Palmeri M, Mehnert J, Silk AW, Jabbour SK, Ganesan S, Popli P, et al. Real-world application of tumor mutational burden-high (TMB-high) and microsatellite instability (MSI) confirms their utility as immunotherapy biomarkers. ESMO Open. 2022;7(1):100336. doi:10.1016/j.esmoop.2021.100336. [Google Scholar] [PubMed] [CrossRef]

24. Brinkman EK, Van Steensel B. Rapid quantitative evaluation of CRISPR genome editing by TIDE and TIDER. Methods Mol Biol. 2019;1961:29–44. doi:10.1007/978-1-4939-9170-9_3. [Google Scholar] [PubMed] [CrossRef]

25. Wang Z, Huang Y, Zhang Y, Zhu H, Amin MR, Chen R, et al. Knockdown of SESN2 exacerbates cerebral ischemia-reperfusion injury through enhancing glycolysis via the mTOR/HIF-1α pathway. CNS Neurosci Ther. 2025;31(3):e70314. doi:10.1111/cns.70314. [Google Scholar] [PubMed] [CrossRef]

26. Haidurov A, Zheltukhin AO, Snezhkina AV, Krasnov GS, Kudryavtseva AV, Budanov AV. p53-regulated SESN1 and SESN2 regulate cell proliferation and cell death through control of STAT3. Cell Commun Signal. 2025;23(1):105. doi:10.1186/s12964-025-02104-3. [Google Scholar] [PubMed] [CrossRef]

27. Zheng C, Zheng L, Yoo JK, Guo H, Zhang Y, Guo X, et al. Landscape of infiltrating T cells in liver cancer revealed by single-cell sequencing. Cell. 2017;169(7):1342–56.e16. doi:10.1016/j.cell.2017.05.035. [Google Scholar] [PubMed] [CrossRef]

28. Ma L, Wang L, Khatib SA, Chang CW, Heinrich S, Dominguez DA, et al. Single-cell atlas of tumor cell evolution in response to therapy in hepatocellular carcinoma and intrahepatic cholangiocarcinoma. J Hepatol. 2021;75(6):1397–408. doi:10.1016/j.jhep.2021.06.028. [Google Scholar] [PubMed] [CrossRef]

29. Zhang W, Zhangyuan G, Wang F, Jin K, Shen H, Zhang L, et al. The zinc finger protein Miz1 suppresses liver tumorigenesis by restricting hepatocyte-driven macrophage activation and inflammation. Immunity. 2021;54(6):1168–85.e8. doi:10.1016/j.immuni.2021.04.027. [Google Scholar] [PubMed] [CrossRef]

30. Li Z, Pai R, Gupta S, Currenti J, Guo W, Di Bartolomeo A, et al. Presence of onco-fetal neighborhoods in hepatocellular carcinoma is associated with relapse and response to immunotherapy. Nat Cancer. 2024;5(1):167–86. doi:10.1038/s43018-023-00672-2. [Google Scholar] [PubMed] [CrossRef]

31. Hossain M, Shim R, Lee WY, Sharpe AH, Kubes P. Gata6+ resident peritoneal macrophages promote the growth of liver metastasis. Nat Commun. 2022;13(1):4406. doi:10.1038/s41467-022-32080-y. [Google Scholar] [PubMed] [CrossRef]

32. Cohen Shvefel S, Pai JA, Cao Y, Pal LR, Bartok O, Levy R, et al. Temporal genomic analysis of homogeneous tumor models reveals key regulators of immune evasion in melanoma. Cancer Discov. 2024:553–77. doi:10.1158/2159-8290.CD-23-1422. [Google Scholar] [CrossRef]

33. Lemaitre L, Adeniji N, Suresh A, Reguram R, Zhang J, Park J, et al. Spatial analysis reveals targetable macrophage-mediated mechanisms of immune evasion in hepatocellular carcinoma minimal residual disease. Nat Cancer. 2024;5(10):1534–56. doi:10.1038/s43018-024-00828-8. [Google Scholar] [PubMed] [CrossRef]

34. Ding J, Liu H, Zhang X, Zhao N, Peng Y, Shi J, et al. Integrative multiomic analysis identifies distinct molecular subtypes of NAFLD in a Chinese population. Sci Transl Med. 2024;16(772):eadh9940. doi:10.1126/scitranslmed.adh9940. [Google Scholar] [PubMed] [CrossRef]

35. Tan J, Fan W, Liu T, Zhu B, Liu Y, Wang S, et al. TREM2+ macrophages suppress CD8+ T-cell infiltration after transarterial chemoembolisation in hepatocellular carcinoma. J Hepatol. 2023;79(1):126–40. doi:10.1016/j.jhep.2023.02.032. [Google Scholar] [PubMed] [CrossRef]

36. Young KZ, Plazyo O, Gudjonsson JE. Targeting immune cell trafficking and vascular endothelial cells in psoriasis. J Clin Invest. 2023;133(9):e169450. doi:10.1172/JCI169450. [Google Scholar] [PubMed] [CrossRef]

37. Chen RY, Yen CJ, Liu YW, Guo CG, Weng CY, Lai CH, et al. CPAP promotes angiogenesis and metastasis by enhancing STAT3 activity. Cell Death Differ. 2020;27(4):1259–73. doi:10.1038/s41418-019-0413-7. [Google Scholar] [PubMed] [CrossRef]

38. Xu X, Ye L, Zhang Q, Shen H, Li S, Zhang X, et al. Group-2 innate lymphoid cells promote HCC progression through CXCL2-neutrophil-induced immunosuppression. Hepatology. 2021;74(5):2526–43. doi:10.1002/hep.31855. [Google Scholar] [PubMed] [CrossRef]

39. Liu H, Ling CC, Yeung WHO, Pang L, Liu J, Zhou J, et al. Monocytic MDSC mobilization promotes tumor recurrence after liver transplantation via CXCL10/TLR4/MMP14 signaling. Cell Death Dis. 2021;12(5):489. doi:10.1038/s41419-021-03788-4. [Google Scholar] [PubMed] [CrossRef]

40. Vitale I, Shema E, Loi S, Galluzzi L. Intratumoral heterogeneity in cancer progression and response to immunotherapy. Nat Med. 2021;27(2):212–24. doi:10.1038/s41591-021-01233-9. [Google Scholar] [PubMed] [CrossRef]

41. Serrano A, Berthelet J, Naik SH, Merino D. Mastering the use of cellular barcoding to explore cancer heterogeneity. Nat Rev Cancer. 2022;22(11):609–24. doi:10.1038/s41568-022-00500-2. [Google Scholar] [PubMed] [CrossRef]

42. Klümper N, Grünwald V, Hartmann A, Hölzel M, Eckstein M. The role of microsatellite instability/DNA mismatch repair deficiency and tumor mutational burden as biomarkers in predicting response to immunotherapy in castration-resistant prostate cancer. Eur Urol. 2024;86(5):388–90. doi:10.1016/j.eururo.2024.04.026. [Google Scholar] [PubMed] [CrossRef]

43. Lenis AT, Ravichandran V, Brown S, Alam SM, Katims A, Truong H, et al. Microsatellite instability, tumor mutational burden, and response to immune checkpoint blockade in patients with prostate cancer. Clin Cancer Res. 2024;30(17):3894–903. doi:10.1158/1078-0432.CCR-23-3403. [Google Scholar] [PubMed] [CrossRef]

44. Campbell KM, Amouzgar M, Pfeiffer SM, Howes TR, Medina E, Travers M, et al. Prior anti-CTLA-4 therapy impacts molecular characteristics associated with anti-PD-1 response in advanced melanoma. Cancer Cell. 2023;41(4):791–806.e4. doi:10.1016/j.ccell.2023.03.010. [Google Scholar] [PubMed] [CrossRef]

45. Stadler ZK, Maio A, Chakravarty D, Kemel Y, Sheehan M, Salo-Mullen E, et al. Therapeutic implications of germline testing in patients with advanced cancers. J Clin Oncol. 2021;39(24):2698–709. doi:10.1200/JCO.20.03661. [Google Scholar] [PubMed] [CrossRef]

46. Agudo J, Miao Y. Stemness in solid malignancies: coping with immune attack. Nat Rev Cancer. 2025;25(1):27–40. doi:10.1038/s41568-024-00760-0. [Google Scholar] [PubMed] [CrossRef]

47. Wu B, Shi X, Jiang M, Liu H. Cross-talk between cancer stem cells and immune cells: potential therapeutic targets in the tumor immune microenvironment. Mol Cancer. 2023;22(1):38. doi:10.1186/s12943-023-01748-4. [Google Scholar] [PubMed] [CrossRef]

48. Rivas MA, Graham D, Sulem P, Stevens C, Desch AN, Goyette P, et al. A protein-truncating R179X variant in RNF186 confers protection against ulcerative colitis. Nat Commun. 2016;7:12342. doi:10.1038/ncomms12342. [Google Scholar] [PubMed] [CrossRef]

49. Hou T, Dan W, Liu T, Liu B, Wei Y, Yue C, et al. Deubiquitinase OTUD5 modulates mTORC1 signaling to promote bladder cancer progression. Cell Death Dis. 2022;13(9):778. doi:10.1038/s41419-022-05128-6. [Google Scholar] [PubMed] [CrossRef]

50. Park SH, Won SY, Kim YC. Identification of a disease-associated germline mutation (A64T) of the ring finger protein 186 gene (RNF186) in Korean patients with ulcerative colitis. Genes Dis. 2024;12(4):101438. doi:10.1016/j.gendis.2024.101438. [Google Scholar] [PubMed] [CrossRef]

51. Ranjan K, Hedl M, Sinha S, Zhang X, Abraham C. Ubiquitination of ATF6 by disease-associated RNF186 promotes the innate receptor-induced unfolded protein response. J Clin Invest. 2021;131(17):e145472. doi:10.1172/JCI145472. [Google Scholar] [PubMed] [CrossRef]

52. Song P, Gao Z, Bao Y, Chen L, Huang Y, Liu Y, et al. Wnt/β-catenin signaling pathway in carcinogenesis and cancer therapy. J Hematol Oncol. 2024;17(1):46. doi:10.1186/s13045-024-01563-4. [Google Scholar] [PubMed] [CrossRef]

53. Carson MD, Nejak-Bowen K. Wnt/β-catenin signaling in liver pathobiology. Annu Rev Pathol Mech Dis. 2025;20:59–86. doi:10.1146/annurev-pathmechdis-111523-023535. [Google Scholar] [PubMed] [CrossRef]

54. Cervantes-Villagrana RD, Albores-García D, Cervantes-Villagrana AR, García-Acevez SJ. Tumor-induced neurogenesis and immune evasion as targets of innovative anti-cancer therapies. Signal Transduct Target Ther. 2020;5(1):99. doi:10.1038/s41392-020-0205-z. [Google Scholar] [PubMed] [CrossRef]

55. Pandey S, Anang V, Schumacher MM. Mitochondria driven innate immune signaling and inflammation in cancer growth, immune evasion, and therapeutic resistance. Int Rev Cell Mol Biol. 2024;386:223–47. doi:10.1016/bs.ircmb.2024.01.006. [Google Scholar] [PubMed] [CrossRef]