Comprehensive Bioinformatics Analysis and Experimental Verification RNF186 Is a Recurrence Signature Gene of Hepatocellular Carcinoma that Promotes Cell Proliferation
Shanbao Ke, Junya Yan, Xiao Feng, Baiyu Li
TL;DR
This study identifies RNF186 as a gene linked to liver cancer recurrence and shows it promotes cancer cell growth, offering new insights for precision treatment.
Contribution
Novel identification of RNF186 as a recurrence signature gene in HCC with experimental validation of its role in cell proliferation.
Findings
A risk model stratified HCC patients into high- and low-risk groups with distinct survival outcomes.
RNF186 was identified as an independent prognostic gene promoting HCC cell proliferation.
High-risk patients showed altered immune profiles and reduced immunotherapy response.
Abstract
Tumor recurrence is a major determinant of poor prognosis in hepatocellular carcinoma (HCC), yet its cellular and molecular basis remains incompletely understood. This study aimed to identify recurrence-associated genes at single-cell resolution and to develop a prognostic model for predicting survival outcomes and immunotherapy responsiveness in HCC. Single-cell RNA sequencing data from 12 primary and 6 recurrent HCC samples were integrated and analyzed to identify genes characteristic of recurrence. After quality control, principal component analysis, and t-SNE-based clustering were used to identify highly variable genes for cell clustering and annotation. Based on macrophage characteristic genes, a recurrence-related risk score was constructed using a LASSO-Cox regression model, and a nomogram integrating clinical variables was developed. Prognostic performance was assessed using…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8| Gene | Forward (5′-3′) | Reverse (5′-3′) |
|---|---|---|
|
| UAAGAUCUCCUCAUUCUGGUU | CCAGAAUGAGGAGAUCUUACA |
|
| ACAAACGUCACUGAAAAUGUG | CAUUUUCAGUGACGUUUGUUG |
|
| UCUUGAACUGCUUGGAAAGGG | CUUUCCAAGCAGUUCAAGAAG |
|
| UGUUGCAGGGUCUUGGUGCAG | GCACCAAGACCCUGCAACAGU |
|
| UCACAUUCUGUGGAGCCAGAA | CUGGCUCCACAGAAUGUGACC |
|
| AUUCUGUGGAGCCAGAAUGAC | CAUUCUGGCUCCACAGAAUGU |
|
| UGUGUUUUGAGAAUAGAACGA | GUUCUAUUCUCAAAACACACU |
|
| UCUCAAAGUGUGUUUUGAGAA | CUCAAAACACACUUUGAGAAA |
|
| UUUCUCAAAGUGUGUUUUGAG | CAAAACACACUUUGAGAAAGA |
|
| UUUCACGGCCUCGGAAGUCCG | GACUUCCGAGGCCGUGAAAAC |
|
| AGCUUCUAGAACCCCGGGGCG | CCCCGGGGUUCUAGAAGCUCC |
|
| UCCAGAACGCUCCGGAAUCCG | GAUUCCGGAGCGUUCUGGAGC |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsFerroptosis and cancer prognosis · Single-cell and spatial transcriptomics · Cancer Immunotherapy and Biomarkers
Introduction
1
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.
Materials and Methods
2
Data Collection
2.1
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.
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.
Single-Cell Sequencing Analysis
2.2
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].
Cell-Cell Interaction Network Analysis
2.3
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.
Risk Score Prognostic Model Construction
2.4
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).
Functional Enrichment Analysis and Protein-Protein Interaction Network Construction
2.5
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, log_2_ |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].
Evaluation of Immune Status, Tumor Stemness, and Genomic Heterogeneity
2.6
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].
Protein Structure Prediction and Mutation Analysis
2.7
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.
Cell Culture and Transfection
2.8
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.
Western Blot Analysis
2.9
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.
Ubiquitination Assays and Co-Immunoprecipitation
2.10
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 Na_2_HPO_4_/NaH_2_PO_4_, 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).
Real-Time Quantitative Polymerase Chain Reaction (RT-qPCR)
2.11
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^′^
Colony Formation Assay
2.12
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.
Cell Viability Assay
2.13
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 × 10^4^ 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.
Statistical Analysis
2.14
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).
Results
3
Single-Cell Analysis of Hepatocellular Carcinoma
3.1
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).
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.
Cell Communication Analysis and Identification of Relapse-Specific Genes
3.2
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).
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.
Risk Score Construction Based on Relapse-Specific Genes Using the LASSO-Cox Algorithm
3.3
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).
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).
Immune Landscape, Tumor Stemness, and Genomic Heterogeneity of Risk Scores
3.4
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).
*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.
Expression and Prognostic Analysis of Risk Score Genes
3.5
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).
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.
Inhibition of RNF186, but Not GRID2 or SLC4A10, Suppresses Huh7 Cell Proliferation
3.6
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).
*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.
RNF186 Promotes Huh7 Cell Proliferation via SESN2
3.7
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).
*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.
Discussion
4
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.
Conclusion
5
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.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Bray 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; 38572751 · doi ↗ · pubmed ↗
- 2Wang Y, Deng B. Hepatocellular carcinoma: molecular mechanism, targeted therapy, and biomarkers. Cancer Metastasis Rev. 2023;42(3):629–52. doi:10.1007/s 10555-023-10084-4; 36729264 · doi ↗ · pubmed ↗
- 3Llovet 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/s 41572-020-00240-3; 33479224 · doi ↗ · pubmed ↗
- 4Wu 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; 37102365 · doi ↗ · pubmed ↗
- 5Tabrizian 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; 38253289 · doi ↗ · pubmed ↗
- 6Llovet 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/s 41571-024-00868-0; 38424197 PMC 11984461 · doi ↗ · pubmed ↗
- 7Sherman M. Recurrence of hepatocellular carcinoma. N Engl J Med. 2008;359(19):2045–7. doi:10.1056/nejme 0807581; 18923166 · doi ↗ · pubmed ↗
- 8Huang X, Song C, Zhang J, Zhu L, Tang H. Circular RN As in breast cancer diagnosis, treatment and prognosis. Oncol Res. 2023;32(2):241–9. doi:10.32604/or.2023.046582; 38186573 PMC 10765117 · doi ↗ · pubmed ↗