Construction of a Prognostic Model for Hepatocellular Carcinoma Based on Macrophage Polarization-Related Genes

Introduction

Hepatocellular carcinoma (HCC) is a highly aggressive malignancy with a significant global impact, resulting in a substantial number of deaths worldwide.^1,2 The limited availability of effective pharmacotherapeutic interventions contributes to the high mortality rate in patients.³ HCC is often diagnosed at an advanced stage because of the absence of noticeable clinical symptoms in the early stages.⁴ Consequently, the overall survival (OS) rates for HCC treatment remain suboptimal, highlighting the need for improved therapeutic strategies.^5,6 To address these challenges, it is crucial to gain a comprehensive understanding of the molecular mechanisms involved in HCC development and to identify novel therapeutic targets. Moreover, the identification of novel prognostic markers and the development of predictive models for HCC prognosis have significant implications for guiding therapeutic interventions in patients with HCC. By identifying reliable prognostic markers and implementing accurate predictive models, clinicians can make informed decisions regarding personalized treatment approaches, ultimately improving patient outcomes.

HCC is a malignancy closely associated with inflammation. The stromal component of HCC primarily consists of fibroblasts, endothelial cells, tumor-infiltrating inflammatory cells, and soluble molecules.^7–10 The tumor microenvironment (TME), which encompasses these cellular components and molecular factors, plays a vital role in regulating neoplastic cell behavior, tumorigenesis, and therapeutic responses.^11,12 Tumor-associated macrophages (TAMs) are the major contributors.¹³ TAMs can be classified into two subtypes: M1 (classical activated macrophages) and M2 (alternatively activated macrophages), each with distinct receptors, cytokines, and chemokines.^14,15 M1-type TAMs express inflammatory factors such as tumor necrosis factor α (TNF-α), IL-1β, IL-6, IL-12, and IL-23, which promote a Th1 (cytotoxic) type response to inhibit tumor growth.¹⁶ Conversely, M2-type TAMs express cytokines such as IL-4, IL-13, IL-8, and IL-10 and growth factors such as epidermal growth factor (EGF), platelet-derived growth factor (PDGF), vascular endothelial growth factor A (VEGFA), and basic fibroblast growth factor (bFGF).¹⁷ The polarization of M1/M2 macrophages and the secretion of these cytokines play critical roles in regulating tumor cell proliferation, tumor cell metabolism, and tumor-related angiogenesis.^18,19 Therefore, it is essential to identify genes and biomarkers related to tumor-associated macrophage polarization for the effective treatment and prognosis of HCC. Such markers will aid in monitoring HCC immunotherapy responses and furthering our understanding of immune infiltration mechanisms.²⁰ Therefore, in this study, we employed bioinformatic methodologies to identify differentially expressed genes in HCC, establish a prognostic model related to macrophage polarization to predict the prognosis of patients with HCC, and validate the robust performance of our model. The findings from this study provide novel insights into the influence of macrophages on HCC and serve as a theoretical foundation for enhancing prognosis prediction and treatment strategies for patients with HCC.

Materials and Methods

Data Acquiring

HCC-related datasets were obtained from the The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) databases. The TCGA-HCC dataset consists of 50 normal tissue samples and 374 cancer tissue samples, with survival information available for 368 of these samples. The ICGC-LIHC dataset served as the validation set and included 232 cancer tissue samples with survival time records available for 215 samples. To identify the relevant genes associated with macrophage polarization (MP), we mined 35 MP-related genes (MPRGs) from published papers. These datasets and gene sets formed the basis of our study to investigate the relationship between MPRGs and HCC prognosis.²¹

Weighted Gene Co-Expression Network Analysis (WGCNA)

MPRGs score was was calculated using the “GSVA” package through ssGSEA analysis.²² To identify gene modules associated with the MPRG score in the TCGA-HCC dataset, we used the “WGCNA” package.²³ Initially, clustering analysis was performed on the samples and outliers were excluded from further analysis. Subsequently, an appropriate soft threshold is determined for the data. A scale-free network was then constructed and genes were divided into several modules using a minimum module size of 30, as determined by the dynamic tree-cutting algorithm. Similar modules were subsequently merged, with the merge CutHeight set to 0.45. The correlation between the MPRG score and each module was assessed. Finally, modules showing the strongest association with the MPRG score (based on a significance threshold of p < 0.05 and |correlation| ≥ 0.3) were identified, and key genes within these modules were selected for subsequent analysis.

Recognition and Enrichment Analysis of Sifferentially Expressed MPRGs (DEMPRGs)

To identify differentially expressed genes (DEGs), we used the “DESeq2” package²⁴ with a significance threshold of adjusted p < 0.05 and |Log2FC| > 2 on the TCGA-HCC dataset.²⁴ DEMPRGs were identified by overlapping them with the key module genes related to the MPRGs score. Volcano Plot was plotted by package “ggplot2” (H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.). A heat map was drawn using the ‘pheatmap package (Kolde R (2019). Pheatmap: Pretty Heatmaps. R package version 1.0.12). GO and KEGG enrichment analyses of DEMPRGs were carried out using the “clusterProfiler” package (adjusted p < 0.05).²⁵

Construction of Risk Model

A total of 368 HCC patients with survival information in the TCGA-HCC dataset (training set) and 215 HCC patients with survival records from the ICGC-LIHC dataset (validation set) were used to establish and validate the MP-related gene signature. The univariate Cox (survival package) (Terry M. Therneau, Patricia M. Grambsch (2000). Modeling Survival Data: Extending the Cox Model Springer, New York. ISBN 0–387-98784-3.) and LASSO regression (glmnet package)²⁶ analyses were initially applied to screen for MP-related prognostic genes in the TCGA-HCC dataset.

Next, we conducted multivariate Cox analysis to identify prognostic genes. We constructed an MP-related prognostic model using multivariate Cox regression analysis in combination with stepwise regression. To categorize patients with HCC into high- and low-risk groups, we utilized the median risk score as a threshold. To evaluate the performance of our prognostic model, we employed several measures. First, Kaplan-Meier (K-M) survival analysis was performed to assess the differences in overall survival between the high- and low-risk groups. Additionally, receiver operating characteristic (ROC) curve analysis was conducted to determine the predictive accuracy of our model. Finally, risk curves were generated to visualize the relationship between the risk score and patient outcomes. Collectively, these analyses provided insights into the effectiveness of our MP-related prognostic model in predicting patient prognosis.

Correlation and Independent Prognostic Analysis Between Clinical Features and Risk Score

First, the differences in the risk score in clinical information (M, N, T, tumor stage, grade, sex, age, and extent of adjacent hepatic tissue inflammation) were analyzed using the Wilcoxon test. Next, the independent prognostic ability of clinical features and risk scores was explored using univariate and multivariate Cox models. Finally, we developed a nomogram model based on independent prognostic factors.

GSEA and TME Analysis

First, we downloaded ‘h.all.v7.5.1.symbols.gmt gene sets in MSigDB and GSEA were performed. The enriched pathways in the two risk subgroups are presented. The differences in the ESTIMATE, stromal, and immune scores between the two risk subgroups were analyzed using the Wilcoxon test. In addition, the correlation between biomarkers and tumor purity, ESTIMATE, stromal, and immune scores was analyzed.

Immunotherapy and Chemotherapy Analyses

First, the expression of immune checkpoints was analyzed using the Wilcoxon test in high-risk groups. Next, TIDE and immunophenoscore (IPS) scores were used to evaluate sensitivity to immunotherapy. The correlation between the biomarkers and the TIDE and IPS scores was further evaluated.

We downloaded invasion-related genes from the cancer SEA database and epithelial-mesenchymal transition (EMT)-related genes from the dbEMT2 database. We then analyzed the correlation between the risk score and invasion-related and EMT-related genes using the Spearman method. Next, the ssGSEA algorithm was used to calculate invasion and EMT scores based on invasion-related and EMT-related genes. Correlations between biomarkers, invasion, and EMT scores were analyzed using the Spearman method. Finally, the “pRRophetic” package (Geeleher (2014) was used. _pRRophetic: Predicts the clinical chemotherapeutic response from pre-treatment tumor gene expression levels. _. R package version 0.5.) software was used to calculate drug sensitivity of each sample. IC50 values of the two risk groups are presented.

Clinical Samples

A total of 64 para-tumor HCC tissues were obtained after informed consent was obtained from patients who underwent liver resection at the Western China Hospital of Sichuan University from 2018 to 2020. All subjects provided informed consent for inclusion before participating in the study, and the tissue was used for future studies.

IHC Staining

Eight tissue specimens from patients with HCC were subjected for immunohistochemistry (IHC) staining. The primary antibodies used for IHC staining were anti-SCN4A polyclonal antibody and DF4513 (Affbiotech). IHC staining was performed according to the manufacturer’s instructions using the Envision+ kit (Dako, Carpinteria, CA, USA). Immunohistochemical images were captured using an Olympus BX63 microscope (Olympus, Japan), allowing the visualization and analysis of the staining patterns and distribution of the target protein in the HCC tissue samples.

Cell Lines and Reagents

The HCC cell line HepG2 and the monocyte cell line THP-1 were obtained from Procell Life Science and Technology (Wuhan, China). HepG2 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM; HyClone, Cytiva, USA) supplemented with 10% fetal bovine serum (FBS; Gibco) at 37 °C. THP-1 monocytes were maintained in RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS; Gibco) and 0.05nM β-mercaptoethanol. THP-1 monocytes were incubated with 100 ng/mL phorbol 12-myristate 13-acetate (PMA, Sigma, P8139) for 24 h in RPMI 1640 medium to differentiate them into macrophages. Macrophage M2 polarization was achieved by incubating PMA-differentiated macrophages with 20 ng/mL of interleukin 4 (GMP-CD03, Novoprotein, Suzhou, China) and 20 ng/mL of interleukin 13 (CC89, Novoprotein, Suzhou, China) for 24 h in the presence of PMA. This protocol allowed us to generate M2-polarized macrophages from THP-1 monocytes for further experiments.

RNA Interference Transfection

SCN4A-specific and negative control siRNAs were designed and obtained from RiboBio (Guangzhou, China). Prior to siRNA transfection, the cells were cultured in medium without fetal bovine serum (FBS) for 24 h. The siRNA constructs were transfected into the cells using Lipofectamine RNAiMAX (Invitrogen), following the manufacturer’s protocol. At 4–6 hours post-transfection, the cells were washed with phosphate-buffered saline (PBS) to ensure the complete removal of the siRNA constructs from the medium. Subsequently, the cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% FBS, which provides the growth factors and nutrients necessary for cell viability and proliferation.

Real-Time Quantitative PCR

Total RNA was extracted from the samples using a High Pure RNA Isolation Kit (RNeasy Mini Kit, 74104, QIAGEN) following the manufacturer’s instructions. A PrimeScript™ RT Reagent Kit (Perfect Real Time, Takara, Japan) was used to obtain complementary DNA (cDNA). Specific primers were used for the amplification of target genes. The expression levels of the selected genes were measured in triplicate using a LightCycler 96 system (Roche). To ensure accurate normalization, the expression of each sample was normalized to that of the reference gene ACTB.

Western Blotting

Cell lysates were prepared using RIPA buffer for protein extraction. For Western blotting, the following primary antibodies were used: anti-SCN4A polyclonal antibody DF4513 (AffBiotech) and anti-β-actin monoclonal antibody TA-09 (ZSBIO, Beijing, China). To detect immune complexes, enhanced chemiluminescence detection reagents (4A BIOTECH, Beijing, China) were used according to the manufacturer’s instructions. This allowed for the visualization of the protein bands on the Western blot, aiding in the analysis and interpretation of the results.

Results

Identification of MPRGs (Macrophage Polarization Related Genes) Score-Related Module Genes (Scale-Free R² = 0.9)

A flowchart of the study is presented in Figure 1. The ssGSEA algorithm was used to calculate the MPRGs score for each sample. Subsequently, WGCNA was used to identify module genes associated with the MPRG score. Initially, we removed five outlier samples, resulting in 419 samples that were used for the subsequent analysis (Figure 2A). Next, we created a scale-free network using a soft threshold of β = 8, which achieved a scale-free R² value of 0.9, indicating an appropriate fit for the data (Figure 2B). We then performed a cluster dendrogram analysis using the dynamic tree cut algorithm, resulting in the division of the network into 40 modules. After merging similar modules, we obtained a final set of 16 modules (Figure 2C). We further investigated the correlation between the modules and MPRGs score. Analysis revealed that the cyan module exhibited a strong correlation with the MPRGs score (Figure 2D). Finally, we selected 3510 genes from the cyan module for subsequent analysis considering their potential relevance to the MPRGs score.

Figure 1 The brief flowchart of this study.

Figure 2 Identification of MPRGs score-related module genes (A) Sample clustering to detect the outliers in the TCGA-HCC dataset. The red line represents the height of the screening outliers. (B) Analyses of the appropriate soft threshold power and minimum mean connectivity to construct topological overlap matrix T(OM). (C) Clustering dendrograms of the co-expression network modules. (D) The correlation analysis between the TME score and module, correlation coefficient, and p-values are shown.

Identification of DEMPRGs and Enrichment Analysis

We identified 2942 differentially expressed genes (DEGs) in the HCC groups of the TCGA-HCC dataset, consisting of 361 downregulated genes and 2581 up-regulated genes (Figure 3A). To further investigate the genes that were differentially expressed and relevant to the MPRG scores, we performed an intersection analysis between the 2942 DEGs and 3510 module genes. This analysis identified 283 differentially expressed MPRGs (DEM-PRGs). These 283 DEM-PRGs represent the overlap between the DEGs and the module genes, signifying their potential significance in relation to both differential expression and MPRG scores. Moreover, we conducted GO and KEGG enrichment analyses to gain insights into the biological functions and signaling pathways associated with the differentially expressed MPRGs (DEMPRGs). These analyses revealed various significant biological processes including epithelial cell proliferation, endothelial cell proliferation, regulation of opsonization, axonogenesis, positive chemotaxis, cell recognition, cell chemotaxis, and sprouting angiogenesis. Additionally, several signaling pathways were found to be related to DEMPRGs, such as cytokine-cytokine receptor interactions, systemic lupus erythematosus, malaria, nicotine addiction, and amphetamine addiction (Supplementary Table 1, Figure 3B and C).

Figure 3 Identification of DEMPRGs and enrichment analysis (A) Volcano plot of DEGs between two clusters with different MP-related gene expression. The red and green points represent up- and down-regulated genes with statistical significance, respectively. (B) The bar plot of enriched Gene Ontology (GO) terms. BP indicates biological process, CC indicates cellular component, and MF indicates molecular function. The length of each column represents the count of genes; the shade of color represents the p value. (C) KEGG analysis of DEGs that were upregulated of DEGs.

Creation and Validation of a MP-Related Risk Model

To screen for genes associated with HCC prognosis in the TCGA-HCC dataset (training set), 49 genes were mined using a univariate Cox model (p<0.05) (Table 1 and Figure 4A). 12 genes (SCN4A, CD5L, EBF3, ADGRB2, HOXD9, HOXD10, KDM8, CTHRC1, CLEC1B, MSC, KCNQ3, and MAPT) were further acquired by LASSO regression analysis (Figure 4B and C). Six prognostic genes (SCN4A, EBF3, ADGRB2, HOXD9, CLEC1B, and MSC) were identified using multivariate COX and stepwise regression analyses to create a risk model (Table 2 and Figure 4D). Meanwhile, a risk score model was developed to calculate a patient’s risk score (Riskscore = h0 (t) * exp ((−0.162 SCN4A expression) + (−0.109 * EBF3 expression) + (0.108 * ADGRB2 expression) + (0.117 * HOXD9 expression) + (−0.129 * CLEC1B expression) + (0.078 * MSC expression))), and according to median of risk score, we classified patients into high and low risk groups. The patient survival time gradually decreased as the risk score increased. The results of the heat map revealed that ADGRB2, HOXD9 and MSC were highly expressed in high-risk patients, while the high expression of CLEC1B, SCN4A and EBF3 was associated with a low risk of HCC (Figure 4E and F). Kaplan-Meier survival analysis revealed that high-risk HCC patients had significantly poorer survival than those with low risk (Figure 4G). Finally, the results of the ROC curves illustrated that the AUC values for OS were 0.737 (1 year), 0.708 (3 years), and 0.742 (5 years) in the TCGA-HCC dataset, indicating the decent performance of the prognostic model (Figure 4H).

Table 1 Result of Univariate Cox Model

Table 2 Prognostic Genes Identified by Multivariate COX

Figure 4 Creation and validation of a MP-related risk model (A) Forest plot of the univariate Cox regression analysis showing the prognostic metabolic genes. Genes with HR < 1 have better overall survival outcomes, while genes with HR > 1 have worse overall survival outcomes. (B) Lasso regression. Finally, 12 MP-related genes were obtained. (C) Lasso filters variables. (D) Forest plot of HR of six genes. (E and F) Expression heat map, risk score distribution and survival status. (G) Survival curve of low and high risk subgroup. (H) Time-dependent ROC curves for the efficacy evaluation of the 6-gene MP-related signature. * P < 0.05, ** P ≤ 0.01.

To further validate the predictive applicability of the risk model, the aforementioned analyses were performed using the ICGC-LIHC dataset (the validation set). Consistent trends were observed in the ICGC-LIHC dataset (Supplementary Figure 1).

Clinical Significance and Independent Prognostic Analysis of the Risk Model

To explore the link between clinical characteristics and risk scores in patients with HCC, the differences in risk scores among different clinical characteristics were evaluated using the Wilcoxon test. The results showed that the risk score was correlated with the HCC grade, HCC stage, and T stage (p<0.05) (Figure 5A).

Figure 5 Clinical significance and independent prognostic analysis of the risk model (A) The correlation between risk score and patient parameters. (B) Risk score of different stages in HCC. (C) Forest maps of the univariate and multivariate Cox regression analysis, including risk score, grade, and stage. (D) ROC curves presented AUC value of risk score. (E) Nomogram based on risk score and clinical factors, including a prediction of 1-, 3-, and 5-year OS. (F) Calibration curves were predicted and observed 1-, 3-, and 5-year survival nomogram consistency.

To determine the prognostic ability of clinicopathological features and risk score, univariate and multivariate Cox analyses were performed. We found that risk score, age, and grade could be used as independent prognostic factors (p<0.05) (Tables 3 and 4, Figure 5B and C). The ROC curves presented an AUC value (0.708) for the risk score (Figure 5D). We then established a nomogram model with risk score, stage, and grade to predict the survival rate of patients at different times (Figure 5E). The nomogram model had a significant predictive effect for patients with HCC (Figure 5F).

Table 3 Result of Independent Prognostic Univariate COX Analysis

Table 4 Result of Independent Prognostic Multivariate COX Analysis

GSEA and TME Analysis in Two Risk Subgroups

GSEA was performed to analyze the differential GO enrichment pathways in the two risk groups. The results revealed that pathways of the G2M checkpoint, KRAS signaling down regulation (DN), and E2F targets were enriched in the high-risk groups, whereas fatty acid metabolism, xenobiotic metabolism, bile acid metabolism, and adipogenesis were related to the low-risk group (Figure 6A).

Figure 6 GSEA analysis showing the enriched pathways of the high-and low-risk groups. (A) Left chart: Multiple GSEA showing enriched pathways in the high-risk group. Right chart: Multiple GSEA showing enriched pathways in the low-risk group. (B) Violin plots of stromal scores, immune scores and ESTIMATE scores between high and low risk groups. (C) The correlation between MPRGs and stromal scores, immune scores, ESTIMATE scores.* P < 0.05, ** P < 0.01.

Next, we analyzed the differences in TME between the two risk subgroups. The results revealed that the immune and stromal scores were significantly difference between the two risk subgroups (Figure 6B). The correlation between each biomarker and tumor purity, ESTIMATE, and stromal and immune scores is presented in (Figure 6C). MSC and CLEC1B were positively correlated with ESTIMATE and stromal and immune scores, whereas MSC and CLEC1B were negatively correlated with tumor purity.

Correlation Analysis of Metastasis/Invasion Indicators in HCC

To better assess the effect of MPRGs on HCC development, a correlation between the expression of invasion-related and EMT-related genes was performed. The results demonstrated a positive correlation between the risk score and the majority of invasive-related genes (Figure 7A), while also revealing diverse associations between the risk score and EMT-related genes (Figure 7B). The correlations of biomarkers with invasion and EMT scores are presented in Figure 7C, in which ADGRB2, MSC and HOXD9 were positively correlated with EMT and invasion scores, indicating that the biomarkers were related to tumor invasion.

Figure 7 Correlation analysis of metastasis/invasion indicators in HCC. (A) The correlation between risk scores and invasion-related related genes. (B) The correlation between risk scores and EMT related genes. (C) The correlation between MPRGs and invasion score, EMT score. ** P < 0.01.

Immunotherapy and Chemotherapy Response in Two Risk Subgroups

Previous studies have shown that macrophages play a dominant role in immunotherapy and are related to the response to immune checkpoint blockade.^27,28 Based on this, we analyzed the expression of many immune checkpoints in the low- and high-risk groups. The results indicated that most of the immune checkpoints were significantly different between the low- and high-risk groups (Figure 8A). Drug sensitivity tests suggested that high-risk patients had lower IC50 values for mitomycin C, cisplatin, gemcitabine, rapamycin, and paclitaxel, indicating higher chemotherapeutic sensitivity to these drugs. Patients in the low-risk group had a lower IC50 value for Doxorubicin, suggesting that they had a higher response to chemotherapy of Doxorubicin (Figure 8B).

Figure 8 Immune checkpoint expressions and chemotherapy response in two risk subgroups (A) Immune checkpoint expressions in high risk group and low risk group. (B) IC50 value of chemotherapeutic drugs in high risk group and low risk group. ns=not significant, * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001.

SCN4A is Highly Expressed in HCC and Inhibition of SCN4A Promotes Macrophage M2 Polarization in HCC Tumor Micro-Environment

As the MP-related risk model was composed of six genes (SCN4A, EBF3, ADGRB2, HOXD9, CLEC1B, and MSC), their expression levels were extracted from the TCGA-HCC dataset. Combined with the sample grouping information of the data set, the expression of biomarkers in the data set was plotted by the wilcox.test method using the R package “ggplot2” (Figure 9A). In the subsequent analysis, we investigating the expression of SCN4A in 64 HCC tissues using qPCR. These results indicate that SCN4A exhibited higher expression levels in HCC tissues than in para-carcinoma tissues. Of the 64 HCC tissues analyzed, 59 showed elevated expression of SCN4A, while 6 displayed decreased expression (Figure 9B). To further validate these findings, IHC staining for SCN4A was performed on 8 HCC specimens. Among these samples, six patients exhibited higher expression levels of SCN4A in HCC tissues than in para-carcinoma tissues (Figure 9C), whereas two patients demonstrated lower expression levels (data not shown). These findings provided additional support for the upregulation of SCN4A in HCC tissues, as observed by qPCR analysis.

Figure 9 SCN4A is highly expressed in HCC tissues (A) MPRGs expressions in HCC of TCGA-HCC dataset. (B) qPCR result of SCN4A expressions in HCC of patients. (C) Representative IHC staining of SCN4A expressions in HCC of patients. Bars=100 μm. * P < 0.05, *** P < 0.001, **** P < 0.0001.

Given the strong association between M2 polarization and tumor development, we investigated whether SCN4A influenced M2 polarization in the HCC tumor microenvironment. To achieve this, we differentiated THP-1 monocytes into M2-type macrophages as previously described (Figure 10A). Subsequently, we employed siRNA to downregulate the expression of SCN4A in HepG2 cells and utilized a plasmid to overexpress SCN4A in HepG2 cells. To establish an indirect co-culture model, supernatants from HepG2 cells transfected with either siSCN4A/siCtrl or SCN4A plasmid/vector for 48 h were transferred to M2-type THP-1 cells (Figure 10B). After stimulation with HepG2 supernatants for 24 h, M2-type THP-1 cells were harvested for qPCR analysis. qPCR results demonstrated that, in comparison to THP-1 cells cultured in the supernatant of HepG2 cells transfected with siCtrl, the M2-type macrophage biomarker CD163 exhibited upregulation upon exposure to the supernatant of HepG2 cells transfected with siSCN4A. Conversely, when compared to THP-1 cells cultured in the supernatant of HepG2 cells transfected with the vector, the M2-type macrophage biomarker CD163 showed downregulation upon exposure to the supernatant of HepG2 cells transfected with the SCN4A plasmid. These findings imply that SCN4A may play a role in inhibiting the M2 polarization of macrophages within the HCC tumor microenvironment (Figure 10C). Next, in order to explore whether SCN4A acts as an oncogene or a tumor suppressor gene, we further investigated the impact of SCN4A on liver cancer cell proliferation. Results from CCK-8 assays and colony formation experiments indicated that downregulation of SCN4A expression in HepG2 cells promotes their proliferation, while overexpression of SCN4A inhibits their proliferation (Figure 10D and E). These findings suggest that SCN4A may function as a tumor-suppressive gene in liver cancer. To further confirm that SCN4A is a tumor-suppressive gene in liver cancer, we explored the changes in the expression levels of several molecules associated with liver cancer progression in HepG2 cells following knockdown/overexpression of SCN4A. We found that upon knocking down SCN4A, the expression of molecules closely linked to hepatocellular carcinoma progression, such as β-catenin, c-myc, CyclinD1, NOTCH-1, HIF-1α, VEGFA, and N-cadherin, increased, while the expression of the liver cancer inhibitory gene E-cadherin decreased. Conversely, in HepG2 cells overexpressing SCN4A, the expression of β-catenin, c-myc, CyclinD1, NOTCH-1, HIF-1α, VEGFA, and N-cadherin decreased, while the expression of E-cadherin, the liver cancer inhibitory gene, increased (Figure 10F). These research findings indicate that SCN4A is a gene that suppresses the progression of liver cancer.

Figure 10 SCN4A expression in HCC cells could regulate M2 macrophage polarization in HCC tumor microenvironment and HCC progression (A) THP-1 cells were induced to differentiate into macrophages displaying M2 polarization by treatment with IL-4 (20 ng/mL) and IL-13 (20 ng/mL), in the presence of PMA (100 ng/mL). (B) Schematic of the indirect co-culture model. (C) qPCR result of CD163 expression in MP macrophages of the indirect co-culture model. (D) CCK-8 experiments showed the growth curve of HepG2 cells. (E) Colony formation assay showed the growth of HepG2 cells. (F)Western blot analysis of SCN4A, β-catenin, c-Myc, CyclinD1, NOTCH-1, HIF-1α, VEGFA, Vimentin, N-Cadherin and E-Cadherin in HepG2 cells. **** P < 0.0001.

Discussion

HCC is a highly malignant tumor with poor prognosis worldwide. TAMs have been found to play a crucial role in all stages of HCC.²⁹ Recent studies have revealed the significance of M1 and M2 macrophage polarization in the development of HCC. M1 TAMs primarily exhibit antitumor functions, whereas M2 TAMs contribute to cancer cell proliferation, aggressiveness, and metastasis through various mechanisms.^14,30 Moreover, M2 TAMs are involved in promoting angiogenesis in liver cancers. This effect is achieved by upregulating the expression of tyrosine kinase receptor Tie-2, human macrophage metalloelastase (HME), and VEGF. These factors collectively enhance angiogenesis within the tumor and facilitate HCC.^31–34

Based on these findings, we integrated MPRGs with clinical parameters to develop a prognostic model for HCC. GO and KEGG pathway analyses revealed that these MPRGs were primarily associated with axon development, epithelial cell function, and played a significant role in HCC prognosis through growth factor activity, extracellular matrix structural constituents, and binding to extracellular components. These results align with those of previous studies, suggesting that morphological changes, such as axon sprouting, occur during macrophage polarization. Consequently, macrophages may interact with epithelial cells, whereas numerous extracellular components can bind to macrophages. Using univariate Cox regression, LASSO regression, and multivariate Cox regression, we selected six MPRGs, ADGRB2, HOXD9, MSC, CLEC1B, SCN4A, and EBF3, to construct the risk model. Among these genes, ADGRB2, HOXD9, and MSC exhibited higher expression levels in patients at higher risk, whereas CLEC1B, EBF3, and SCN4A were highly expressed in patients at a lower risk of HCC. These findings provide valuable insights into the potential utility of MPRGs as a prognostic marker for HCC.

To gain a deeper understanding of MPRGs included in the model, we conducted a literature review to explore their known functions. The first gene, ADGRB2 (Adhesion G Protein-Coupled Receptor B2), also referred to as BAI2 (brain angiogenesis inhibitor 2), is a p53-target gene that specifically encodes a brain-specific angiogenesis inhibitor.³⁵ To date, no study has specifically explored the role of ADGRB2 in cancer have been published. Therefore, the function of ADGRB2 in cancer remains unclear. Another MPRG, HOXD9 (Homeobox D9), belongs to a superfamily of genes that regulates the development and control of various cellular processes. In HCC, HOXD9 promotes metastasis by regulating ZEB1.³⁶ Furthermore, HOXD9 has been frequently used in risk prediction models for HCC using bioinformatics methods, indicating its potential significance in prognostic evaluation.^37–40 The MSC gene encodes musculin, which acts as a transcriptional repressor, impeding myogenesis and activating E-box-dependent muscle genes. Its involvement in the TME suggests that it may play a role in HCC prognosis by influencing the TME dynamics. Further investigations are needed to elucidate the precise mechanisms underlying HCC.⁴¹

These results are consistent with our findings that ADGRB2, HOXD9 and MSC are highly expressed in high-risk patients. CLEC1B (C-Type Lectin Domain Family 1 Member B), also known as C-type lectin-like receptor 2 (CLEC-2), is expressed on platelets, Kupffer cells and other immune cells, and binds to various ligands.⁴² Notably, the expression of CLEC1B, and Programmed death-ligand 1 (PD-L1) has been identified as a predictor of clinical outcomes in HCC cases associated with tumor hemorrhage.⁴³ Furthermore, low expression of CLEC1B has been observed in HCC, highlighting its potential role as a biomarker for suppressing the microenvironment, stemness, and hypoxia.^44,45 EBF3 (Early B cell Factor 3) is a member of the early B-cell factor (EBF) family, which comprises DNA binding transcription factors. EBF proteins play critical roles in B cell differentiation, bone development, and neurogenesis. Although there are reports suggesting that EBF3 may promote metastasis in certain tumors,^46,47 it is generally considered to function as a tumor suppressor.^48–51 SCN4A (Sodium Voltage-Gated Channel Alpha Subunit 4) is predominantly expressed in skeletal muscles, and its mutations have been associated with various myotonia and periodic paralysis disorders.⁵² Interestingly, a correlation between SCN4A and HCC prognosis was initially identified using bioinformatics.^53,54 This highlights the potential significance of SCN4A in predicting the outcomes of previously unknown HCC cases. Our findings are consistent with these results, indicating that CLEC1B, EBF3, and SCN4A are highly expressed in patients at a lower risk.

Following the modeling of these six MPRGs, we proceeded with relevant verification steps. First, we established a nomogram model incorporating the risk score, stage, and grade, which demonstrated a favorable predictive effect for patients with HCC. Subsequently, we analyzed the correlation between tumor purity, ESTIMATE scores, stromal scores, immune scores, and the six MPRGs. Additionally, we examined the relationship between each biomarker and tumor purity, ESTIMATE scores, stromal scores, and immune scores. The results revealed that MSC and CLEC1B were positively correlated with ESTIMATE and stromal and immune scores, while demonstrating a negative correlation with tumor purity. Furthermore, we analyzed the association of each biomarker with invasion scores, which indicated a positive correlation between ADGRB2 and MSC with tumor invasion.

According to previous studies, it was the third time that SCN4A was identified to be strongly correlated with HCC prognosis by our research. However, the nature of this relationship remains unclear. To delve deeper into its role, we conducted qPCR and IHC staining of SCN4A and discovered that it exhibited higher expression levels than adjacent para-tumor tissues. To further investigate the effect of SCN4A, we induced THP-1 cells into an M2-like subtype and established an indirect co-culture model using the supernatants from HepG2 cells. Interestingly, When the SCN4A gene was knocked down in HepG2 cells, the supernatants obtained induced M2 polarization in THP-1 cells. Conversely, when the SCN4A gene was overexpressed in HepG2 cells, the resulting supernatants reduced M2 polarization in THP-1 cells. This observation aligns with our risk model, which implicates SCN4A in patients at lower risk of HCC. Based on these findings, we can infer that SCN4A deficiency promotes polarization of M2 macrophages in HCC. Thus, SCN4A may act as a negative regulator of M2 macrophages in HCC. Furthermore, we confirmed through cellular experiments that SCN4A is a gene that inhibits HCC.

We also analyzed the enrichment pathways of the high and low-risk groups. The E2F targets signaling, G2M checkpoint signaling and KRAS signaling were enriched in high-risk group and signalings related to adipogenesis, bile acid metabolism, fatty acid metabolism were enriched in low risk group. E2F targets are a group of genes that encode transcription factors in higher eukaryotes, participating in the cell cycle regulation and synthesis of DNA in mammalian cells. Activation of E2F targets promotes cell cycle progression. G2M checkpoint signaling is also a signaling pathway involved in the cell cycle. KRAS is one of the important members of the RAS family, and its main downstream signaling pathways include important signaling transduction pathways such as PI3K/AKT. It plays a crucial regulatory role in signaling pathways involved in tumor cell growth and angiogenesis processes. However, the enriched pathways in the low-risk group are primarily associated with metabolism, particularly lipid metabolism. The above research findings indicate that in the future development of anti-tumor drugs, efforts should be focused on targeting the E2F, G2M checkpoint, and KRAS signaling pathways. Additionally, exploring the development of drugs that regulate tumor lipid metabolism could also be beneficial.

In addition, we compared immunotherapy and chemotherapy responses between high- and low-risk subgroups. This analysis revealed that the patients in the high-risk subgroup exhibited greater sensitivity to mitomycin C, cisplatin, gemcitabine, rapamycin, and paclitaxel. Conversely, patients in the low-risk subgroup showed a better response to doxorubicin-based chemotherapy. These findings have crucial implications for guiding clinical decision making in terms of chemotherapy selection for patients with HCC.

This study has several limitations that should be acknowledged. First, our experimental focus was primarily on SCN4A, with only eight sections of IHC staining conducted to examine its effect on M2 polarization. Further investigations are needed to determine whether SCN4A overexpression can reduce CD163 expression in M2-like THP-1 cells, explore the impact of SCN4A on M1 polarization in HCC, and elucidate the precise mechanism by which SCN4A regulates M1/M2 polarization. Additionally, it is important to note that the expression of the other five MPRGs in the HCC samples was not thoroughly examined in this study. Future studies should aim to clarify their expression and role in HCC. Moving forward, we plan to expand our sample size to validate and optimize our model, specifically focusing on elucidating the mechanisms by which these genes modulate M1 or M2 polarization in HCC. By addressing these limitations, we can gain a deeper understanding of the complex interplay between the tumor microenvironment and HCC progression. Despite these limitations, our study provides a TME score-related prognostic model for HCC, contributing valuable insights that may aid the discovery of new treatment strategies for this disease.