Integrative Multi-Omics and Functional Validation Reveal the Role of the TACE Refractoriness-Associated Gene ATP1B3 in Hepatocellular Carcinoma

Introduction

Hepatocellular carcinoma (HCC) is the most common type of liver cancer, accounting for approximately 90% of liver cancer cases.¹ Despite the recent plateau in its incidence, HCC remains the third leading cause of cancer-related mortality worldwide.^2,3 The poor overall survival (OS) rate of HCC patients (5-year OS < 20%) is primarily attributed to the fact that the majority of cases are diagnosed at intermediate or advanced stages, where curative treatment options are no longer viable.^4,5

Transarterial chemoembolization (TACE) is the most widely used treatment for unresectable HCC and is recommended as the standard therapy for Barcelona Clinic Liver Cancer (BCLC) stage B HCC.^6,7 It is also utilized as a palliative treatment for advanced disease, a bridging therapy to liver transplantation, and as a tumor downstaging strategy prior to surgical resection.^1,4,8 TACE delivers cytotoxic agents through intra-arterial injection, followed by embolization of tumor-feeding blood vessels, inducing cytotoxic and ischemic effects within the tumor.⁹ Despite its widespread use, the efficacy of TACE varies significantly among patients. The objective response rate of TACE is approximately 52%, with OS rates of 70%, 52%, and 32% at 1, 2, and 5 years, respectively, and a median OS of 19.4 months.¹⁰ These findings reveal contrasting patient outcomes, as 30% of patients die within the first year due to poor responses, while 32% survive beyond five years. This suggests the existence of a subgroup of patients who derive substantial benefit from TACE. The phenomenon of patients showing poor response to TACE is termed “TACE refractoriness” or “TACE failure”, first introduced by the Japan Society of Hepatology in 2010.¹¹ In addition, emerging evidence indicates that TACE refractoriness is strongly associated with an unfavorable prognosis in HCC patients.^12,13 Therefore, identifying the specific genes and pathways underlying TACE refractoriness is critical for advancing HCC research and improving therapeutic outcomes.

Recent advancements in sequencing methods and bioinformatic algorithms have transformed multi-omics technology into an indispensable tool for HCC research.^14–16 Nevertheless, the inherent complexity and immense scale of sequencing data present substantial challenges in pinpointing critical genes. Machine learning, a branch of artificial intelligence, provides powerful tools for analyzing complex datasets by uncovering patterns and predicting outcomes from unknown data. Additionally, machine learning systems can continuously update and refine their performance by integrating new knowledge.¹⁷ With ongoing advancements in machine learning, more accurate and sophisticated algorithms enable reducing data dimensionality and identifying key disease-related genes with greater precision.¹⁸ For example, Zhang et al identified macrophage-specific PLA2G7 as a key biomarker associated with prognosis and immunotherapy response in HCC through single-cell RNA sequencing and machine learning analysis, demonstrating the potential of such approaches to uncover novel therapeutic targets.¹⁹

In this study, we conducted unsupervised clustering analysis to stratify HCC patients treated with TACE into two distinct subgroups based on their gene expression profiles. These subgroups demonstrated significant differences in clinical outcomes and molecular characteristics, highlighting the heterogeneity underlying TACE refractoriness. Among the differentially expressed genes identified between the two clusters, ATP1B3 emerged as a key gene associated with TACE refractoriness, as determined through univariate Cox regression, lasso regression, differential expression, and survival analyses. To further elucidate the role of ATP1B3 in HCC, we integrated transcriptomics, genomics, and scRNA sequencing data to perform multi-omics analyses, including functional, mutational, immune, cell subset composition, and cell communication analyses. The functional significance of ATP1B3 in HCC progression was validated through in vitro assays, including differential expression, migration, invasion, cell cycle, chemotherapy sensitivity, and apoptosis experiments. Our findings suggest that the TACE refractoriness-associated gene ATP1B3 may serve as a promising therapeutic target for improving treatment outcomes in HCC patients.

Materials and Methods

Data Collection

A total of seven public datasets were analyzed in this study, covering scRNA-seq, high-throughput RNA sequencing (RNA-seq), microarray data, and clinical treatment cohorts. The HCC scRNA-seq datasets GSE149614²⁰ and GSE242889²¹ were obtained from the Gene Expression Omnibus (GEO). Gene expression data and complete clinical information for the TCGA-LIHC and GSE14520²² cohorts were downloaded from TCGA and GEO, respectively, while somatic mutation and copy number variation (CNV) data for TCGA-LIHC were retrieved from cBioPortal. Additional microarray datasets, GSE104580 and GSE109211,²³ were retrieved from GEO. GSE104580 contains data on TACE treatment, while GSE109211 focuses on sorafenib therapy. Information regarding patient response to treatment in both datasets was provided by the original authors. Lastly, the Imvigor210 cohort,²⁴ an immune therapy dataset, was included for further analysis.

Unsupervised Clustering and Functional Profiling in TACE-Treated Cohort

Unsupervised clustering analysis was performed on 104 TACE-treated patients from the GSE14520 cohort using the ConsensusClusterPlus package.²⁵ To identify the most variable genes, the median absolute deviation (MAD) was calculated for each gene, and the top 150 genes with the highest variability were selected for clustering. The optimal number of clusters was determined based on cumulative distribution function (CDF) curves, changes in the area under the CDF, and the consensus heatmap. This was further validated using the proportion of ambiguous clustering (PAC) method.²⁶ Kaplan-Meier survival analysis was conducted to assess the ability of the identified subgroups to predict OS and recurrence-free survival (RFS) in TACE-treated patients. Differentially expressed genes (DEGs) between the two subgroups were identified using the limma package, with thresholds of adjusted P value < 0.05 and |logFC| > 1. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the DEGs were performed using the ClusterProfiler package.²⁷ To identify potential therapeutic compounds and target pathways, the Connectivity Map (CMap) database was utilized. For CMap analysis, 150 upregulated and 150 downregulated genes were submitted simultaneously to identify candidate compounds based on similarities in gene expression profiles.

Workflow of Gene Selection

First, univariate Cox regression analysis was performed to identify TACE prognosis-related genes among the upregulated genes associated with TACE refractoriness, as determined through differential expression analysis. Subsequently, univariate Cox regression was also applied to the TCGA cohort to further identify genes associated with HCC prognosis. To refine the selection, key genes were identified using lasso regression. Finally, core gene for further investigation were comprehensively determined by evaluating three critical factors: (1) differences in gene expression between HCC tissues and adjacent non-tumor tissues, (2) their prognostic significance in HCC patients, and (3) their associations with HCC malignancy progression.

Identification of ATP1B3 Co-Expressed Genes and Enrichment Analysis

The TCGA-LIHC cohort was analyzed to identify genes co-expressed with ATP1B3. The Spearman correlation coefficient was used to evaluate the correlation between ATP1B3 expression and other genes. The top 150 positively and negatively correlated genes were selected for GO and KEGG pathway enrichment analysis. In addition, the Spearman correlation analysis and gene set enrichment analysis (GSEA) were further performed to investigate the association between ATP1B3 expression and hallmark gene sets from the Molecular Signatures Database (MSigDB).

Somatic Mutation and Copy Number Variation Analysis

The genomic variation landscape was examined by analyzing somatic mutations and CNV data. Differences in somatic mutations between ATP1B3-high and ATP1B3-low groups were assessed by evaluating single nucleotide polymorphisms (SNPs), insertions and deletions (INDELs), tumor mutation burden (TMB), and mutation frequencies. Visualizations were created using the maftools package.²⁸ The top 20 frequently mutated genes (FMGs), identified based on mutation frequency, were considered potential drivers of tumorigenesis.²⁹ CNV profiles were compared between the two groups, emphasizing the top 10 genes with the highest frequencies of amplification (AMP) or homozygous deletion (HOMDEL) in each group.

Comprehensive Analysis of Immune-Omics Molecular Characterization and Immunotherapy Response

The relative infiltration levels of 28 types of immune cells were quantified using the ssGSEA method,^30,31 and their correlation with ATP1B3 expression was calculated using Spearman correlation analysis. Using the IOBR package,³² we collected previously published gene signatures associated with immunotherapy response, immune suppression, and immune exclusion. A unified scoring method was applied to calculate the enrichment scores for each sample, enabling a comprehensive analysis of immunological differences between patients with high or low ATP1B3 expression. To evaluate immunotherapy response, the TIDE algorithm³³ was used, and the findings were further validated in the Imvigor210 cohort.

scRNA-Seq Analysis

The Seurat R package³⁴ was used for processing scRNA-seq data. Quality control standards for scRNA-seq data included the following criteria: (1) the number of detected features (nFeature_RNA) > 500, and (2) the percentage of mitochondrial genes (percent_mito) < 20%. Subsequently, the Harmony R package³⁵ was used to integrate data from fifteen samples. Dimensionality reduction and visualization were performed using t-SNE. Cell clustering was conducted using the “FindNeighbors” and “FindClusters” functions in Seurat. Finally, cell clusters were annotated to ten known cell types based on canonical marker genes. CNVs were assessed by calculating a moving average of relative expression using InferCNV.³⁶

Cell-Cell Communication Analysis

Cellular communication analysis was conducted using the R package CellChat (v1.6.1). The normalized expression matrix was first imported into a CellChat object through the createCellChat function. Default preprocessing steps were then performed, including identifyOverExpressedGenes, identifyOverExpressedInteractions, and projectData. Potential ligand-receptor interactions were identified using the computeCommunProb, filterCommunication, and computeCommunProbPathway functions. The results were visualized with the ktplots package (v2.3.0).

Cell Culture and Transfection

HCC-LM3 and Huh-7 cell lines were obtained from the National Collection of Authenticated Cell Cultures (Shanghai, China). THLE-2 cell line was obtained from Meisen Chinese Tissue Culture Collections (Zhejiang, China). HCC-LM3 cells were cultured in DMEM medium supplemented with 10% fetal bovine serum (FBS, Vivacell, C04001-050) and 1% penicillin-streptomycin. Huh-7 cells were maintained in MEM medium containing 10% FBS and 1% penicillin-streptomycin. THLE-2 cells were cultured in BEGM medium supplemented with 10% FBS and 5 ng/mL EGF.

Small interfering RNA (siRNA) targeting ATP1B3 was synthesized by TsingkeBio (Beijing, China). Transient transfection was performed using Opti-MEM™ reduced serum medium (Thermo Fisher Scientific), siRNA (5 µL), and Lipofectamine™ RNAiMAX Transfection Reagent (5 µL, Invitrogen, 13778030). The components were thoroughly mixed and incubated at room temperature for 10 minutes. The resulting mixture was then added to 6-well plate cultures and incubated for 48 hours before further analysis.

Real-Time Quantitative RT–PCR

Total cellular RNA was extracted using the RNA-Quick Purification Kit (ES Science, RN001) following the manufacturer’s protocol. RNA concentration and purity were determined using the NanoQ spectrophotometer (NanoQ, China). Complementary DNA (cDNA) was synthesized using the Hifair^® AdvanceFast 1st Strand cDNA Synthesis Kit (YisBio, 11150ES60) according to the manufacturer’s instructions. Relative quantitative PCR was performed using the Applied Biosystems ABI 7500 system through the measurement of real-time SYBR Green fluorescence. Relative gene expression levels were calculated using the 2^−ΔΔCt method to determine fold changes, normalized to GAPDH expression levels as an internal control.

Western Blotting

Total proteins were extracted using RIPA lysis buffer on ice and quantified using the Bradford assay. Proteins were separated by SDS-PAGE and transferred onto a PVDF membrane using electroblotting. The membrane was blocked with 5% skim milk dissolved in TBST for 1 hour at room temperature and incubated overnight at 4°C with a specific primary antibody. After washing three times with TBST, the membrane was incubated with HRP-linked secondary antibodies (1:10,000 dilution) for 1 hour at room temperature. The enhanced chemiluminescent detection kit was employed to visualize the proteins bands. The primary antibodies used in this study included anti-ATP1B3 (Proteintech, 11142-1-AP, 1:1000) and anti-β-tubulin (Cell Signaling Technology, 2146S, 1:1000).

Cell Apoptosis Analysis

After 48 hours of siRNA transfection, the culture medium was replaced with fresh complete medium, and the cells were cultured for an additional 48 hours. Both floating and adherent cells were then collected and stained with the Annexin V-FITC/PI Staining Kit (Beyotime, C1062L) according to the manufacturer’s protocol. The samples were incubated in the dark at room temperature for 15 minutes and then analyzed using a flow cytometer (BD FACSCanto, BD Biosciences, USA). Results were analyzed using FlowJo software.

Cell Cycle Analysis

HCC cells were collected after 48 hours of siRNA transfection and were stained with the cell cycle kit (MULTI SCIENCES, CCS012) according to the manufacturer’s protocol. Cells were resuspended in PI staining solution and the mixture was vortexed for 10 seconds to ensure homogeneity and incubated in the dark at room temperature for 30 minutes. After staining, the samples were analyzed using a flow cytometer (BD FACSCanto, BD Biosciences, USA).

Wound Healing Assay

When the cell density reached over 90% confluence, a sterile 200-μL pipette tip was used to create a cross-shaped scratch at the bottom of the well. The wells were then washed three times with PBS to remove suspended cells. The cells were cultured in serum-free medium under standard culture conditions (37°C, 5% CO₂), and images of the scratch area were captured every 12 hours using an inverted microscope. The relative wound areas were quantified using ImageJ software.

Invasion Assay

Matrigel (Corning, 354234) was diluted with serum-free medium at a ratio of 1:3, and the mixture was evenly spread onto the membrane of a Transwell insert (Corning, 3422). The inserts were incubated at 37°C with 5% CO₂ for 1 hour to allow the Matrigel to polymerize. Then, 100 µL of serum-free medium containing 5 × 10⁴ cells was added to the upper chamber, and 500 µL of complete medium with 10% FBS was added to the lower chamber as a chemoattractant. The plates were incubated at 37°C with 5% CO₂ for 18–24 hours. After incubation, non-invaded cells on the upper surface of the membrane were removed using a cotton swab. The invaded cells on the lower surface of the membrane were fixed with 4% paraformaldehyde for a period of time, stained with crystal violet, and rinsed gently with distilled water. The membranes were air-dried, and images were captured using an inverted microscope.

CCK8 Assay

Drug sensitivity was determined using the CCK8 assay. Cells were seeded in a 96-well plate (10,000 cells/well) and cultured at 37°C with 5% CO₂ for 24 hours. Cisplatin (MCE, HY-17394) was then added to the wells and continued incubation for 48 hours. Following treatment, CCK-8 solution (beyotime, C2015L) was added to each well to evaluate cell viability, and the plate was incubated at 37°C for an additional 3 hours. The absorbance was measured at 450 nm using a microplate reader (Tecan, Tecan infinite m200).

Statistical Analyses

Statistical analyses were performed using R (version 4.3.1) and GraphPad Prism (version 9.5). Spearman correlation coefficient was used to evaluate the correlation between two continuous variables. The chi-square test was applied to compare categorical variables, while continuous variables were analyzed using the Wilcoxon rank-sum test or the t-test based on the distribution of the data. Survival analyses were conducted using the Kaplan-Meier method, and differences between survival curves were assessed with the Log rank test. All P values were two-sided, and statistical significance was defined as P < 0.05.

Results

Unsupervised Clustering Identifies Two Subtypes with Distinct Responses to TACE Treatment

To explore the molecular heterogeneity of HCC patients treated with TACE, we performed unsupervised clustering analysis on the gene expression profiles of 104 patients. The consensus heatmap, CDF curves and the PAC statistic indicated that the optimal number of clusters was two (Figures 1A and S1A–C). Principal component analysis (PCA) further supported theses clustering results (Figure 1B). A Log rank test revealed significant differences between the two groups in OS (P < 0.001) and RFS (P = 0.002), indicating better prognoses for patients in cluster 1 for both OS and RFS (Figure 1C and D). We then divided the TACE-treated patients into two subgroups, namely adjuvant TACE and post-recurrence TACE. Results showed that cluster 1 continued to have superior OS and RFS compared to cluster 2 (Figure 1E–H). Additionally, we analyzed HCC patients who underwent only surgical resection. Cluster 1 showed a greater survival benefit compared to surgery alone, whereas cluster 2 did not exhibit significant survival advantages (Figure 1I and Table S1). This suggests heterogeneity in the effectiveness of TACE treatment between clusters 1 and 2.

Figure 1 Unsupervised clustering reveals distinct prognostic clusters of HCC patients treated with TACE. (A) Unsupervised consensus clustering illustrates two clusters in 104 HCC patients treated with TACE. (B) PCA illustrates the distribution of two clusters. (C) OS and (D) RFS of patients receiving TACE treatment. (E) OS and (F) RFS of patients receiving adjuvant TACE treatment. (G) OS and (H) RFS of patients receiving post-recurrence TACE treatment. (I) Comparison of OS among cluster1, cluster2 and patients who underwent surgical resection without additional therapy. All data are from the GSE14520 cohort.

Comparative Clinical and Molecular Analysis Between Cluster 1 and Cluster 2

To further explore the differences between cluster 1 and cluster 2, we performed a comprehensive comparative analysis. Clinical characteristic evaluation revealed significant differences between the two clusters, particularly in BCLC stage and patient status (all P < 0.05). Cluster 2 patients experienced more aggressive disease progression and higher mortality rates compared to cluster 1 (Figure 2A). Differential gene expression analysis identified 95 upregulated and 255 downregulated genes between the two clusters (Figure 2B). GO and KEGG functional enrichment analyses highlighted the involvement of several metabolism-related pathways (Figure 2C and D). Immune cell infiltration analysis revealed higher levels of immunosuppressive cells, such as MDSCs and Tregs, in cluster 2 (Figure 2E). Using the TIDE algorithm, we found that cluster 2 had higher TIDE scores, indicating a poorer response to immunotherapy (Figure 2F). Furthermore, we explored potential therapeutic strategies using CMap database, which matches gene expression profiles with drug responses. We identified 20 compounds with potential therapeutic effects for patients in cluster 2 (Figure 2G), targeting pathways that could be utilized to develop novel therapeutic options (Figure 2H). These findings provide a comprehensive characterization of the clinical and molecular differences between cluster 1 and cluster 2, suggesting distinct underlying biological mechanisms.

Figure 2 Different clinicopathological characteristics, functional enrichment, immune infiltration, and CMap analysis between the two clusters. (A) Pie charts showing the differences in clinicopathological factors between the two clusters in the GSE14520 cohort. Status refers to the survival status of patients who received TACE treatment during the follow-up period (up to 67.3 months). (B) Volcano plot of differential genes. (C) GO and (D) KEGG enrichment analysis. (E) Differences in immune cell infiltration in patients of two clusters. (F) TIDE analysis between the two clusters. (G) Heatmap of enrichment score generated from potential therapeutic compounds. (H) Description of mode of action of compounds targeting corresponding molecular pathways. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. Abbreviation: ns, no significance.

Identification of Key Genes Associated with TACE Refractoriness in Cluster 2

To identify genes driving TACE refractoriness in cluster 2, we analyzed 75 highly expressed genes using univariate Cox regression. Forty-four genes were significantly associated with OS (Figure S2). Validation in TCGA HCC cohort, using both univariate and lasso regression analyses, further narrowed this list to 11 key genes: MAFF, MARCKSL1, TRAF5, CXCL8, NT5DC2, CDC20, G6PD, ATP1B3, BICC1, NCAPG, SPP1 (Figure 3A–C). We then systematically evaluated the expression patterns and prognostic significance of these 11 key genes. The results showed that five genes, namely NT5DC2, CDC20, ATP1B3, NCAPG, and SPP1, were significantly overexpressed in tumor tissues and associated with poor prognosis (Figures S3A–J, S4A–J, S5A–J, S6A–J, S7A–J, and S8A–J). Among these, ATP1B3 stood out due to its pronounced overexpression in both non-paired and paired tumor tissues, and its strong correlation with poor prognosis across multiple survival endpoints, including OS, PFI, DFI, and DSS (all P < 0.05) (Figure 3D–I). Furthermore, ATP1B3 expression was significantly upregulated in advanced HCC, including higher pathological grades (Grade 3 and 4) and advanced clinical stages (stage III and IV) (Figure 3J and K), and was associated with HCC diagnosis (Figure 3L). Additionally, ATP1B3 had the highest hazard ratio in the univariate Cox analysis. Therefore, we selected ATP1B3 as the gene of interest for further investigation into its role in TACE refractoriness.

Figure 3 Identification of key genes associated with TACE refractoriness. (A) Univariate Cox regression analyses. (B) Lasso analysis of 11 genes, lambda min = 0.024. (C) Cross-validation of results of lasso analysis (lambda min = 0.024, number = 11). Comparison of ATP1B3 expression levels between (D) tumor and tumor-adjacent tissues, (E) as well as between tumor and paired paracancerous tissues in HCC patients. The impact of ATP1B3 on (F) OS, (G) DSS, (H) DFI, and (I) PFI in HCC analyzed using Kaplan-Meier analysis. Box plots showing ATP1B3 expression between different groups stratified by (J) histological grade and (K) clinical stage. (L) The diagnostic value of ATP1B3 for HCC evaluated using ROC curve analysis (AUC = 0.832, 95% CI = 0.776–0.889).

Functional Characterization of ATP1B3 in HCC

To further investigate the role of ATP1B3 in HCC, we performed a correlation analysis between ATP1B3 and all other genes (Figure 4A). The top 150 positively and 150 negatively correlated genes were identified and subjected to GO and KEGG enrichment analyses. GO analysis revealed significant associations with metabolic catabolic processes, regulation of coagulation, and enzyme regulatory activities (Figure 4B), while KEGG analysis linked ATP1B3 to key metabolism-related pathways, including cholesterol and pyruvate metabolism (Figure 4C). We also performed a correlation analysis between ATP1B3 expression and 50 hallmark gene sets, identifying significant correlations with Myc targets V1, glycolysis, and inflammatory response pathways (Figure 4D). GSEA further demonstrated that ATP1B3 high-expression groups were enriched in pathways related to Myc targets V1, TNF-α signaling via NF-κB, and epithelial-mesenchymal transition, suggesting potential involvement in cancer-related biological processes (Figure 4E). In addition, we depicted the landscape of genomic variations between the high and low ATP1B3 expression groups. Analysis of the top 20 FMGs revealed differences in mutation frequencies, SNPs, INDELs, and TMB (Figure 4F). Specifically, TP53 mutations were more frequent in the high ATP1B3 group (Figure S9). CNV analysis, including the top 10 amplified and homozygous deleted genes, further highlighted distinct genomic landscapes between the two groups (Figure 4G). These findings underscore the diverse genomic features associated with varying ATP1B3 expression levels.

Figure 4 Functional roles and genomic alternations of ATP1B3 in HCC. (A) Volcano plot depicting the distribution of genes co-expressed with ATP1B3. (B) GO and (C) KEGG enrichment analyses of co-expressed genes. (D) Correlation analysis between hallmark genesets and ATP1B3 expression. (E) GSEA analysis of ATP1B3. (F) The waterfall plot illustrates the differences in frequently mutated genes in HCC between the two groups. (G) Amplified and homozygously deleted genes identified between the two groups.

Association of ATP1B3 with Sorafenib Response, TACE Response, and Immunotherapy Outcomes

Given ATP1B3’s association with poor survival and tumor prognosis in HCC, we next investigated its role in modulating response to therapies such as sorafenib and TACE, as well as its impact on the tumor immune microenvironment. First, we analyzed the correlation between ATP1B3 expression and sorafenib response. ATP1B3 expression was significantly higher in sorafenib non-responders (Figure 5A), and the high ATP1B3 group had a larger proportion of non-responders, suggesting that elevated ATP1B3 contributes to sorafenib resistance (Figure 5B). Moreover, sorafenib-related targets, including BRAF, RAF1, FLT1, KIT, FLT3, and FLT4, were upregulated in the ATP1B3 high-expression group (Figure 5C). Next, we examined ATP1B3 expression in TACE responders and non-responders. Similar to sorafenib, ATP1B3 expression was significantly higher in TACE non-responders (Figure 5D), and a larger proportion of patients with high ATP1B3 expression failed to respond to TACE (Figure 5E). To further understand the potential mechanisms behind these associations, we analyzed the relationship between ATP1B3 expression and immune cell infiltration. ATP1B3 expression was positively correlated with immunosuppressive cells, such as Tregs and MDSCs (Figures 5F and S10), and associated with higher immune exhaustion, exclusion, and suppression scores, suggesting its role in promoting an immunosuppressive tumor microenvironment (Figure 5G–I). Furthermore, ATP1B3 was positively linked to immune checkpoint molecules, including PDCD1, CD274, CTLA4, IDO1, and LAG3 (Figure 5J), indicating its potential role in immune evasion. Consistently, the ATP1B3 high-expression group exhibited higher TIDE scores, suggesting a possible association between ATP1B3 expression and immunotherapy response (Figure 5K). To further explore this, we analyzed the Imvigor210 cohort, which includes patients treated with immune checkpoint inhibitors. High ATP1B3 expression was associated with worse prognosis (Figure 5L), and a larger proportion of these patients experienced SD or PD following immunotherapy (Figure 5M).

Figure 5 ATP1B3 expression correlates with treatment responses and immune features in HCC. (A) ATP1B3 expression in responders and non-responders to sorafenib in the GSE109211 cohort. (B) Proportion of responders in ATP1B3 low- and high-expression groups for sorafenib in the GSE109211 cohort. (C) Expression of sorafenib targets in ATP1B3 high- and low-expression groups. (D) ATP1B3 expression in responders and non-responders to TACE in the GSE104580 cohort. (E) Proportion of responders in ATP1B3 low- and high-expression groups for TACE in the GSE104580 cohort. (F) Heatmap illustrates the relationship between ATP1B3 expression and immune cell infiltration. Red indicates a positive correlation, while blue indicates a negative correlation. Blank cells in the heatmap represent correlations that are not statistically significant (P > 0.05). (G) Immune suppression, (H) exhaustion, and (I) exclusion signatures in ATP1B3 high- and low-expression groups. (J) Chord plot showing correlations between immune checkpoint genes and ATP1B3 expression. (K) TIDE scores in ATP1B3 high- and low-expression groups. (L) Kaplan-Meier survival curves for OS in ATP1B3 high- and low-expression groups in the Imvigor210 cohort. (M) Objective response rates (CR, PR, SD, PD) in ATP1B3 high- and low-expression groups. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. Abbreviation: ns, no significance.

scRNA Analysis of ATP1B3 Expression and Tumor Microenvironment Interactions

To gain a deeper understanding of the role of ATP1B3 in HCC at the single-cell level, we integrated 15 HCC samples from two single-cell RNA sequencing datasets. Batch effects from different datasets and samples were corrected (Figures 6A, B and S11A, B). After dimensionality reduction and unsupervised cell clustering, ten distinct cell subsets were identified based on marker gene expression, including HCC cells (ALB, HP, FGA, FGB, AFP, SPINK1, GPC3, AKR1C1, TF), endothelial cells (VWF, PECAM1), fibroblasts (COL1A1, COL3A1, COL1A2), macrophages (C1QA, C1QB), monocytes (S100A8, LYZ), dendritic cells (CLEC10A, CD1C), T cells (CD3D, CD3E), natural killer (NK) cells (NKG7, GNLY, KLRD1), B cells (CD79A, MS4A1), and plasma cells (MZB1, IGHG1) (Figure 6C and D). In addition, inferCNV analysis confirmed the malignant nature of HCC cells (Figure S12). We then visualized ATP1B3 expression across these cell subsets and found that ATP1B3 was expressed in all of them (Figure 6E). To further explore the functional differences between ATP1B3-positive and ATP1B3-negative HCC cell subpopulations, we performed GO and KEGG enrichment analyses. GO analysis showed that differentially expressed genes were enriched in translation, metabolic processes, immune response, and enzyme regulation (Figure 6F). KEGG analysis revealed that ATP1B3-positive subpopulations were associated with pathways such as antigen processing and presentation, oxidative phosphorylation, glycolysis/gluconeogenesis, Th17 cell differentiation, and glutathione metabolism. In contrast, ATP1B3-negative subpopulations were enriched in pathways including retinol metabolism, drug metabolism-cytochrome P450, fatty acid degradation, xenobiotic metabolism, and tryptophan metabolism (Figure 6G). Next, we conducted a cell-cell communication analysis to explore the interactions among various cell subsets. Interestingly, ATP1B3-positive HCC cell subpopulations exhibited stronger communication with macrophages, DCs, and monocytes compared to ATP1B3-negative HCC cell subpopulations (Figure 6H). Specific ligand-receptor (LR) pairs mediating these interactions included CCL5-CCR1 and LGALS9-CD45 (Figure 6I). Taken together, these findings suggest that while ATP1B3 is widely expressed across various cell subsets, ATP1B3-positive HCC cell subpopulations have distinct functional properties compared to ATP1B3-negative cell subpopulations.

Figure 6 Single-cell RNA sequencing reveals ATP1B3-associated cellular heterogeneity and intercellular communications in HCC. t-SNE plots of two HCC datasets (A) before and (B) after batch correction. (C) t-SNE plot of 10 identified cell subsets. (D) Bubble plot of signature gene expression across cell subsets. (E) ATP1B3 expression across cell subsets. (F) GO and (G) KEGG enrichment analyses of differentially expressed genes between ATP1B3-positive and ATP1B3-negative HCC cell subpopulations. (H) Heatmap of interaction profiles between cell subsets. (I) Dot plot depicting LR interactions between ATP1B3-positive or ATP1B3-negative HCC cells and immune cells (macrophages, DCs, and monocytes).

Experimental Validation of ATP1B3 Function in HCC Cell Lines

To investigate the expression of ATP1B3 in HCC, we used RT-qPCR to assess its levels in HCC cell lines (HCC-LM3, Huh-7) and the normal hepatic cell line THLE-2, revealing higher expression in HCC cells compared to normal cells (Figure 7A). To further elucidate the role of ATP1B3 in HCC initiation and progression, we knocked down ATP1B3 in HCC-LM3 and Huh-7 cells, confirming knockdown efficiency using RT-qPCR and WB (Figure 7B–D). Subsequently, functional assays were performed to evaluate the biological effects of ATP1B3 knockdown. A wound healing assay demonstrated that ATP1B3 knockdown reduced the migratory capacity of HCC cells. (Figure 7E). Consistently, transwell invasion assays showed a significant reduction in the invasive potential of ATP1B3-silenced cells (Figure 7F and G). Flow cytometry analysis further revealed that ATP1B3 knockdown induced significant G1 phase arrest in HCC-LM3 and Huh-7 cells (Figure 8A and B) and increased apoptotic cell populations (Figure 8C and D). Additionally, CCK8 assay demonstrated enhanced cisplatin sensitivity in ATP1B3-silenced cells (Figure 8E). Collectively, these findings confirm that ATP1B3 silencing inhibits HCC cell migration and invasion, induces G1 phase arrest, promotes apoptosis, and enhances chemosensitivity, highlighting its critical role in HCC progression and chemoresistance.

Figure 7 In vitro validation of ATP1B3 function in HCC. (A) qRT-PCR analysis of ATP1B3 expression in HCC-LM3, Huh-7, and THLE-2 cells. Knockdown efficiency of ATP1B3 confirmed by qRT-PCR in HCC-LM3 and Huh-7 cells (B), and by Western blot in (C) HCC-LM3 and (D) Huh-7 cells. (E) Scratch wound-healing assay showing reduced migration of Huh-7 cells after ATP1B3 knockdown. Transwell invasion assay demonstrating decreased invasiveness of ATP1B3-silenced (F) HCC-LM3 and (G) Huh-7 cells. **P < 0.01; ***P < 0.001; ****P < 0.0001.

Figure 8 Effects of ATP1B3 knockdown on cell cycle, apoptosis, and chemosensitivity in HCC cells. Flow cytometry analysis of cell cycle in HCC-LM3 (A) and Huh7 (B) cells. Flow cytometry analysis of apoptosis in HCC-LM3 (C) and Huh7 (D) cells. (E) CCK8 assay for cell viability. ****P < 0.0001.

Discussion

The heterogeneous response to TACE among HCC patients underscores the complexity of tumor biology and the multifactorial nature of treatment resistance. While TACE is globally recognized as the standard of care for intermediate-stage HCC and is also commonly utilized in advanced cases, approximately 30% of patients exhibit poor responses, resulting in significantly shortened survival.^6,7,10 This phenomenon, termed TACE refractoriness, has been strongly associated with adverse clinical outcomes and underscores the need to address the underlying mechanisms driving this refractoriness. Given the molecular and microenvironmental heterogeneity of HCC, identifying the genetic and molecular pathways that contribute to TACE refractoriness is critical for advancing our understanding of HCC pathophysiology and for developing more effective, personalized therapeutic strategies. Among these pathways, ATP1B3, which encodes the β3 regulatory subunit of Na⁺/K⁺-ATPase, has emerged as a potential candidate. Beyond its canonical role in maintaining cellular ion homeostasis, ATP1B3 influences disease progression through various mechanisms, including promotion of gastric cancer progression via PI3K/AKT signaling,³⁷ modulation of HBV replication through NF-κB activation,³⁸ and association with poor prognosis and immune infiltration in HCC.³⁹ However, its specific role in HCC progression and treatment resistance remains poorly understood.

In our study, we used an unsupervised clustering method to identify molecular subtypes based on the gene expression profiles of TACE-treated HCC patients.⁴⁰ This analysis identified two distinct clusters with significant differences in prognosis, clinical characteristics, immune microenvironment composition, and therapeutic vulnerabilities. Notably, Cluster 2 was associated with worse OS and RFS compared to Cluster 1 in both adjuvant and post-recurrence TACE procedure. Interestingly, when we compared these clusters to surgically resected HCC patients, we observed that the OS of cluster 2 tended to be even worse than that of surgical resection patients, suggesting that TACE may not provide survival benefits for cluster 2. The poor prognosis of cluster 2 may be attributed to its immunosuppressive tumor microenvironment, as we observed higher infiltration of MDSCs and Tregs in this cluster, consistent with previous studies linking MDSCs and Tregs infiltration to poor outcomes across multiple cancers.^41–43 TIDE analysis further indicated that cluster 2 was less likely to respond to immunotherapy. There findings underscore the need to further identify the key genes and pathways driving the aggressive behavior and treatment resistance of cluster 2, which could provide valuable insights into its underlying biology and reveal novel therapeutic targets.

The rapid advancement of machine learning has facilitated its widespread application in predicting disease states, uncovering biomarkers, identifying key prognostic genes, and more.^44–46 Using a combination of the lasso algorithm and univariate Cox regression analysis, 11 prognostic genes associated with TACE refractoriness were identified. To further narrow down the candidate genes, their expression patterns in tumor tissues, associations with prognostic outcomes (OS, DFS, DFI, and PFI), and relevance to HCC progression were systematically analyzed. Among these candidates, ATP1B3 emerged as a pivotal gene, characterized by its markedly elevated expression in HCC tissues, strong association with malignant progression, and the highest HR among the 11 genes. Moreover, elevated ATP1B3 expression was significantly correlated with poor prognosis and demonstrated diagnostic potential for HCC. These findings highlight ATP1B3 as a key gene potentially driving TACE refractoriness in cluster 2, providing a rationale for further investigation into its functional role in HCC.

To elucidate the functional role of ATP1B3 in HCC progression, a correlation analysis was conducted to identify ATP1B3 co-expressed genes, followed by enrichment analysis of these genes. The results revealed that ATP1B3 was associated with metabolism-related pathways, including cholesterol and pyruvate metabolism, both of which have been implicated in tumor development.^47,48 GSEA further demonstrated that high ATP1B3 activated tumor-associated pathway, such as Myc targets V1 and epithelial-mesenchymal transition, which are known to promote tumor progression.^49,50 In addition, TP53 mutation is associated with poor prognosis in HCC.⁵¹ Notably, we observed that the TP53 mutation rate was higher in the ATP1B3-high group, further supporting the association between high ATP1B3 expression and adverse clinical outcomes. Regarding treatment responsiveness, ATP1B3 expression was negatively correlated with the efficacy of TACE and sorafenib, with patients in the ATP1B3-high group exhibiting lower response rates and poorer therapeutic outcomes. Moreover, the immune microenvironment is closely linked to prognosis and the success of immunotherapy.⁵² Our results indicated that tumors with high ATP1B3 expression exhibited an immunosuppressive state, accompanied by increased expression of immune checkpoint-related genes. Consistently, the TIDE algorithm predicted that patients with high ATP1B3 expression were less likely to respond to immunotherapy. Supporting this, analysis of the Imvigor210 cohort further demonstrated that high ATP1B3 expression was associated with worse prognosis and reduced immune response rates following immunotherapy. Together, these findings support the role of ATP1B3 as a potential driver of HCC progression and treatment resistance.

To further elucidate the role of ATP1B3 within the tumor microenvironment, single-cell RNA sequencing was performed to analyze the microenvironment at single-cell resolution.⁵³ ATP1B3 was detected in both malignant tumor cells and immune cells. Comparative analysis of ATP1B3-positive and ATP1B3-negative malignant cell subpopulations revealed enrichment of metabolic and immune-related pathways, consistent with our earlier findings. Furthermore, cell-cell communication analysis demonstrated stronger interactions between ATP1B3-positive malignant cells and myeloid cells, mediated by CCL5-CCR1 and LGALS9-CD45 ligand-receptor interactions, both of which are known to promote tumor progression and immune evasion.^54,55 Finally, in vitro validation experiments confirmed that ATP1B3 is highly expressed in HCC cell lines. Knocking down ATP1B3 inhibited tumor cell proliferation, migration, and invasion, while inducing cell cycle arrest and apoptosis. These results further underscore the role of ATP1B3 as a key driver of HCC progression and treatment resistance. Despite these findings, several limitations should be noted. First, while we explored ATP1B3-associated pathways and functions, the precise molecular mechanisms underlying its role in HCC progression and treatment resistance remain to be fully elucidated. In addition, the clinical data on treatment responsiveness, particularly for TACE and sorafenib, were primarily retrospective and require prospective validation. Furthermore, due to the lack of publicly available cohorts for immunotherapy in HCC patients, we utilized the Imvigor210 cohort as an external cohort to validate the predictive power of ATP1B3 on immune response. Moreover, scRNA-seq data revealed that ATP1B3 expression is higher in immune cells, such as monocytes and macrophages, than in HCC cells, suggesting a potential role in immune regulation. While this study focused on tumor cells due to their strong association with HCC progression and diagnostic potential, the role of ATP1B3 in immune cells remains unclear and warrants further investigation. Addressing these gaps will provide deeper insights into the biological and therapeutic relevance of ATP1B3.

Conclusion

In this study, we identified two molecularly distinct HCC subtypes with divergent responses to TACE, showing differences in survival and immune microenvironment. Through univariate and lasso analyses, ATP1B3 was identified as a key gene associated with TACE refractoriness. Its overexpression correlated with metabolic reprogramming, immune suppression, and resistance to TACE, sorafenib, and immunotherapy. Functional assays demonstrated that ATP1B3 promotes HCC migration, invasion, and chemoresistance, while its knockdown inhibited these phenotypes and enhanced cisplatin sensitivity. These findings position ATP1B3 as both a biomarker for patient stratification and a promising therapeutic target in HCC.