VDAC1 and VDAC2 as prognostic biomarkers and therapeutic targets in hepatocellular carcinoma

Introduction

Background

Metabolic reprogramming and immune dysregulation are core features driving tumor progression across cancers (1). Mitochondria, as central hubs for metabolism and signaling, sustain energy supply and regulate immune responses within the tumor microenvironment (2). Voltage-dependent anion channels (VDACs), located on the mitochondrial outer membrane, mediate transmembrane transport of molecules such as ATP and calcium, thereby influencing metabolism, apoptosis, and immune signaling (3). Among them, VDAC1 is widely overexpressed in various cancers and is associated with a poor prognosis, whereas VDAC2 exhibits more context-dependent effects. Despite overlapping functions, VDAC1 and VDAC2 display distinct characteristics, suggesting potential complementary roles in tumor biology (4,5).

Rationale and knowledge gap

Liver hepatocellular carcinoma (LIHC) is the most common primary liver cancer, characterized by high molecular and immunological heterogeneity, which complicates diagnosis and treatment (6). Current clinical indicators often fail to capture the tumor’s full molecular and immune features. While some studies have explored VDAC1/2 in liver cancer, comprehensive evaluations of their mechanistic roles, regulation, and clinical value—especially their joint contribution—remain limited (7,8). Understanding their complementary functions could improve prognostic biomarker systems and inform therapeutic strategies.

Objective

This study aims to systematically evaluate the expression patterns, prognostic value, and mechanistic roles of VDAC1 and VDAC2 in hepatocellular carcinoma (HCC) through integrated multi-omics analyses and functional experiments, and to explore their potential as complementary prognostic biomarkers and therapeutic targets. We present this article in accordance with the TRIPOD and MDAR reporting checklists (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1887/rc).

Methods

Analysis of VDAC1 and VDAC2 expression profiles

RNA-seq count data for 33 cancer types were obtained from The Cancer Genome Atlas (TCGA) via the GDC portal (https://portal.gdc.cancer.gov). The full names of tumor types along with their corresponding abbreviations are detailed in Table S1. Transcript abundance was normalized to transcripts per million (TPM) values, followed by log2(TPM+1) transformation. Samples with incomplete clinical annotations were excluded. Differential expression of VDAC1 and VDAC2 was assessed between tumor and matched normal tissues, and further analyzed across pathological stages.

Validation was performed using Gene Expression Omnibus (GEO) datasets: GSE14520 (225 LIHC tumors and 219 non-tumor tissues). Protein expression was confirmed using immunohistochemistry data from the Human Protein Atlas (HPA) Pathology module. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Prognostic and diagnostic evaluation of VDAC1 and VDAC2

We downloaded STAR-counts data and clinical information for 33 cancer types from the TCGA database (https://portal.gdc.cancer.gov), normalized the data in TPM format using log2(TPM+1), and retained samples with both RNA-seq and clinical data for analysis.

First, we performed receiver operating characteristic (ROC) curve analysis to assess the diagnostic performance of VDAC1 and VDAC2 in distinguishing tumors from normal tissues in LIHC, using the area under the curve (AUC) for accuracy.

Next, we conducted a pan-cancer Kaplan-Meier survival analysis to evaluate the prognostic value of VDAC1 and VDAC2. For LIHC, subgroup analyses assessed the impact of VDAC1/2 expression on clinical variables like T stage, N stage, M stage, and grade. Log-rank tests and hazard ratios (HRs) with 95% confidence intervals (CIs) were used to assess survival differences.

Univariate and multivariate Cox regression analyses were then performed, with forest plots generated to show HRs, p-values, and 95% CIs. A Nomogram was constructed using multivariate results to predict the 1-, 3-, and 5-year recurrence rates for LIHC.

Genetic and epigenetic alteration analysis

Genetic alterations of VDAC1 and VDAC2, including copy number variations (CNVs), were analyzed using TCGA datasets via the cBioPortal platform (https://www.cbioportal.org). Mutation frequencies across 33 cancer types were obtained from the “Cancer Types Summary” module, while detailed mutation sites and variant classifications were examined in the “Mutations” module. The prognostic significance of these alterations was evaluated using the “Survival” module.

To validate and expand these analyses, the Gene Set Cancer Analysis (GSCA) platform (http://bioinfo.life.hust.edu.cn/GSCA/) was used. The “Mutation” module was applied to assess CNV frequencies and survival differences between mutant and wild-type groups. The prognostic impact of CNV status was further evaluated based on overall survival (OS).

For epigenetic analysis, promoter region CpG methylation levels were retrieved using UALCAN (http://ualcan.path.uab.edu/). Associations between methylation status (high vs. low) and clinical outcomes were analysed through GSCA.

Tumor immune microenvironment and immunoregulatory analysis

To evaluate the immunological roles of VDAC1 and VDAC2, multi-dimensional analyses were conducted using publicly available datasets. Immune cell infiltration correlations were assessed via TIMER2.0 (http://timer.cistrome.org/). The “Immune Gene” module was used to quantify associations between gene expression and six immune cell types (B cells, CD4⁺ T cells, CD8⁺ T cells, macrophages, neutrophils, and dendritic cells). The “Immune Outcome” module assessed the combined effect of gene expression and immune infiltration on OS in the LIHC cohorts.

CIBERSORTx (https://cibersortx.stanford.edu/) was applied using its LM22 signature matrix to estimate the relative abundance of 22 immune cell subtypes. Correlation analyses were performed between VDAC1/2 expression and individual immune cell fractions.

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was applied to estimate immune escape potential and to predict likely responsiveness to immune checkpoint blockade (ICB). Additionally, correlations between VDAC1/2 expression and ten representative immune checkpoint genes (including PDCD1, CD274, CTLA4, HAVCR2, TIGIT, LAG3, PDCD1LG2, SIGLEC15, IGSF8, and IT-PRIPL1) were examined based on TCGA transcriptomic data.

Cell culture and small interfering RNA (siRNA) interference

Human HCC cell lines (Hep3B, PLC) were obtained from the Chinese Academy of Medical Sciences (Beijing, China). Cells were maintained in RPMI-1640 medium (Gibco) supplemented with 10% fetal bovine serum (FBS; Gibco), 100 U/mL penicillin, and 100 µg/mL streptomycin (Solarbio) at 37 ℃ under 5% CO₂ humidified atmosphere.

Gene-specific silencing was performed using small interfering RNAs (siRNAs) targeting VDAC1 and VDAC2, with a non-targeting siRNA serving as a negative control (NC; all from GenePharma, Shanghai, China). siRNA sequences were:

VDAC1-1: 5'-GAAUAGCAGCCAAGUAUCAG-3'
VDAC1-2: 5'-ACACUAGGCACCGAGAUUA-3'
VDAC2-1: 5'-GGAGUGUAUAAACCUUGGUUGUGAU-3'
VDAC2-2: 5'-ACUCUGAGGCCUGGUGUGAAGCUUA-3'

Transfection was conducted using HiPerFect reagent (QIAGEN) according to the manufacturer’s protocol. Briefly, cells seeded in 6-well plates were transfected with 10 nM siRNA in antibiotic-free complete medium. Cells were harvested 72 hours post-transfection for downstream analyses.

Colony formation assay

To assess clonogenic potential post-siRNA transfection, cells were plated in 24-well plates at a density of 1×10³ cells per well and maintained under standard culture conditions for 10–14 days. Cultured medium was replaced every 2–3 days. Upon completion, cell colonies were rinsed with phosphate-buffered saline (PBS), fixed using 4% paraformaldehyde for 15 minutes, and stained with 0.1% crystal violet at room temperature for 30 minutes. After washing and air-drying, colonies containing over 50 cells were quantified either manually or with ImageJ. Each experiment was independently repeated three times.

Cell proliferation assay

Cell proliferation was measured using the Cell Counting Kit-8 (CCK-8; KeyGEN BioTECH, Nanjing, China) following the manufacturer’s protocol. In brief, 1.0×10³ cells were seeded into each well of a 96-well plate containing 100 µL of complete medium. At 0, 24, 48, and 72 hours post-seeding, 10 µL of CCK-8 solution was added to each well and incubated at 37 ℃ for 1–2 hours. Absorbance at 450 nm was then recorded using a microplate reader (Infinite M200 PRO; TECAN, Switzerland) to determine cell viability. All assays were carried out in triplicate.

Apoptosis assay

For apoptosis analysis, cells were seeded into 6-well plates at a density of 2×10⁵ cells per well and treated with the indicated reagents for 48 hours. After treatment, cells were harvested, washed twice with cold phosphate-buffered saline (PBS), and stained using the Annexin V Alexa Fluor 647 Apoptosis Detection Kit (KeyGEN BioTECH, Nanjing, China) following the manufacturer’s protocol. The stained cells were then analyzed by flow cytometry using a Gallios flow cytometer (Beckman Coulter Inc., Brea, CA, USA). Data were processed using CytExpert software (Beckman Coulter Inc., Brea, CA, USA). Live cells were defined as Annexin V⁻/PI⁻, while apoptotic cells included both early apoptotic (Annexin V⁺/PI⁻) and late apoptotic/necrotic populations (Annexin V⁺/PI⁺).

Reactive oxygen species (ROS) detection

Cells were seeded in six-well plates to 70–80% confluence, then washed with PBS and incubated with 2 µM DCFH-DA (ROS probe) for 30 minutes at 37 ℃ according to the Beijing Yuntian Biotech ROS Detection Kit instructions. After incubation, cells were washed three times with PBS to remove excess probe. Intracellular ROS levels were analyzed by flow cytometry (BD Accuri C6 Plus), and data were processed using FlowJo software. Fluorescence intensity was directly correlated with ROS levels, allowing comparison of different treatment conditions.

Transwell assay

Transwell migration assays were performed using 24-well plates with 8 µm pore inserts. Hep3B and PLC cells (2×10⁴ cells/well) were seeded in the upper chamber with serum-free medium, and the lower chamber contained 10% FBS as a chemoattractant. After 24 hours of incubation at 37 ℃, non-migrated cells were removed, and the migrated cells on the lower surface were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. Migration was assessed by counting the number of migrated cells under a light microscope (Leica) at 200× magnification.

Quantitative real-time PCR (qRT-PCR)

qRT-PCR analysis was performed to assess the efficiency of VDAC1/2 knockdown and to evaluate their impact on genes associated with epithelial-mesenchymal transition (EMT). Total RNA was isolated 48 hours post-transfection using TRIzol^® Reagent (Thermo Fisher Scientific), followed by quantification with a NanoDrop 2000 spectrophotometer (A260/A280 ratios ≥1.8). We reverse transcribed 1 µg RNA to cDNA using the PrimeScript RT^® Kit (Takara). We amplified target genes with SYBR^® Pre-mix Ex Taq™ II (Vazyme) on a CFX96 Touch™ system (Bio-Rad).

Relative gene expression was calculated via the 2^−ΔΔCt method using GAPDH as endogenous control. Primer sequences for VDAC1, VDAC2, and EMT markers (E-cadherin, N-cadherin, Vimentin, Snail) are listed in Table S2.

Western blot analysis

To validate VDAC1/2 knockdown efficiency and assess EMT pathway modulation, total cellular protein was extracted using RIPA lysis buffer (Biosharp) supplemented with 1 mM phenylmethylsulfonyl fluoride (PMSF). Protein concentration was quantified via the BCA Protein Assay Kit (Beyotime Biotechnology). Equal amounts of protein (10 µg) were denatured in loading buffer at 100 ℃ for 5 min, separated by 10% SDS-PAGE, and electrophoretically transferred to 0.22-µm PVDF membranes (pre-activated in methanol for 30 sec) at 300 mA for 90 min. Membranes were blocked with 5%non-fat milk for 1 h at room temperature, then incubated overnight at 4 ℃ with primary antibodies (detailed in Table S3), followed by IRDye 800CW-conjugated secondary antibodies (1:10,000 dilution) for 1 h at room temperature. Protein signals were detected using an Odyssey Dual-Color Infrared Imaging System (LI-COR Biosciences) with GAPDH as a loading control.

Statistical analysis

Statistical analyses in this study were conducted using the above online databases and R software (46, R version 3.6.3), as previously described. Experimental data were analyzed using GraphPad Prism 9.0 (GraphPad Software, La Jolla, CA, USA). Differences between the two groups were compared using Student’s t-test, and results are expressed as mean ± standard deviation. For apoptosis assays involving multiple groups and conditions, statistical significance was determined using two-way analysis of variance (two-way ANOVA) followed by Tukey’s post hoc test. Statistical significance was reported at *P<0.05, **P<0.01, ***P<0.001.

Results

The expression of VDAC1 and VDAC2 are upregulated in LIHC tissue

To comprehensively characterize the expression patterns of VDAC1 and VDAC2 across human cancers, we first analyzed TCGA pan-cancer transcriptomic data from 33 tumor types (Figure 1A-1D). VDAC1 was significantly upregulated in the majority of tumors, with particularly high expression in COAD, LIHC, BRCA, BLCA, CHOL, KIRC, LUAD, PRAD, READ, THCA, UCEC, HNSC, and STAD as well as in CESC, ESCA, and KICH. In contrast, VDAC2 displayed greater inter-cancer variability. Elevated VDAC2 levels were observed in LUSC, LIHC, and CHOL, whereas downregulation was noted in BLCA, KICH, READ, THCA, and COAD.

Figure 1 Expression profiles of VDAC1 and VDAC2 in pan-cancer. (A-D) The expression levels of VDAC1 and VDAC2 across 33 cancer types and 24 adjacent normal tissues as analyzed using the TCGA database. (E) Expression patterns of VDAC1 and VDAC2 in the GSE14520 cohorts. (F) Representative IHC images showing VDAC1 (https://www.proteinatlas.org/ENSG00000213585-VDAC1/tissue) and VDAC2 (https://www.proteinatlas.org/ENSG00000165637-VDAC2/tissue) expression in human cancer tissue samples obtained from the HPA database (Original magnification, ×200). (G) VDAC1 and VDAC2 expression levels across different cancer stages, as determined by pathological staging. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001. LIHC, liver hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; TPM, transcripts per million.

To validate these patterns, we examined two independent GEO datasets (GSE14520). Both confirmed the TCGA results, showing elevated VDAC1 and VDAC2 expression in LIHC (Figure 1E).

Protein-level validation using the HPA further supported these findings (Figure 1F). Both VDAC1 and VDAC2 showed strong cytoplasmic staining in tumor tissues from LIHC, consistent with mRNA upregulation.

To explore the relationship between gene expression and tumor progression, we next examined VDAC1 and VDAC2 expression across clinical stages (I–IV) in LIHC (Figure 1G). VDAC1 expression showed a stepwise increase from stage I to stage IV, with a statistically significant elevation observed in stage III compared to stage I (P=0.02). Although the differences among all stages did not consistently reach statistical significance—likely due to the limited sample size in stage IV—the overall upward trend suggests that VDAC1 may be gradually activated during tumor progression. A similar pattern was observed for VDAC2. Compared with normal liver tissues, VDAC2 expression was significantly elevated across all tumor stages (P<0.001), and its levels also increased progressively with clinical stage. Notably, expression in stage II (P=0.007) and stage III (P<0.001) was significantly higher than in stage I, reinforcing the possible stage-dependent role of VDAC2 in LIHC progression.

VDAC1 and VDAC2 serve as robust diagnostic markers and unfavorable prognostic factors in LIHC

To evaluate the clinical relevance of VDAC1 and VDAC2, we first assessed their diagnostic performance. ROC analyses revealed that both VDAC1 and VDAC2 exhibited strong discriminatory power between tumor and normal samples, with AUC values of 0.917 (95% CI: 0.887–0.947) and 0.902 (95% CI: 0.869–0.936), respectively (Figure 2A). Moreover, univariate and multivariate Cox regression analyses confirmed the independent prognostic significance of both genes (Figure 2B,2C). After adjusting for clinical variables, high expression of VDAC1 (adjusted HR =1.678, 95% CI: 1.182–2.382, P=0.004) and VDAC2 (adjusted HR =1.477, 95% CI: 1.038–2.100, P=0.03) were both independently associated with poorer OS. Finally, time-dependent calibration curves demonstrated good predictive performance of nomogram models incorporating VDAC1 and VDAC2, with strong concordance between predicted and observed 1-, 3-, and 5-year survival (Figure 2D).

Figure 2 Diagnostic, Cox regression and predictive performance analyses of VDAC1 and VDAC2 in LIHC. (A) ROC curves of VDAC1 and VDAC2 showing their predictive accuracy, with AUC values and 95% CIs. (B,C) Univariate and multivariate Cox regression analyses evaluating the prognostic significance of VDAC1 (B) and VDAC2 (C) along with clinicopathological variables. HRs with 95% CIs are presented. (D) Calibration curves of nomograms for 1-, 3-, and 5-year overall survival prediction based on VDAC1 and VDAC2. AUC, area under the curve; CI, confidence interval; FPR, false positive rate; HR, hazard ratio; LIHC, liver hepatocellular carcinoma; M, metastasis; N, node; ROC, receiver operating characteristic; T, tumor; TPR, true positive rate.

We next investigated their prognostic significance across multiple cancer types (Figure S1). VDAC1 overexpression was associated with reduced OS in several cancers, including BRCA, LUAD, HNSC, LIHC, PAAD, and SKCM, while showing improved OS in COAD and KIRC (Figure S1). Similarly, VDAC2 displayed heterogeneous associations, predicting worse outcomes in LIHC, BLCA, LUAD, and PAAD, but improved survival in KIRC (Figure S1). Focusing on LIHC, Kaplan-Meier analyses demonstrated that high expression levels of both VDAC1 (HR =1.612, P=0.007) and VDAC2 (HR =1.585, P=0.009) were significantly correlated with poorer OS, suggesting their roles as unfavorable prognostic biomarkers (Figure 3A-3D). Stratified subgroup analyses further highlighted stage-specific effects. Elevated VDAC1 expression predicted reduced survival in patients with T3&T4 tumors (P=0.02), N0 status (P=0.005), and M0 status (P=0.005), whereas no significant associations were observed in the T1&T2 or Grade 1&2 subgroups (Figure 3A-3D). By contrast, VDAC2 exhibited a broader prognostic impact across early-stage subgroups, where high expression significantly correlated with worse survival in patients with T1&T2 tumors (P=0.03), N0 (P=0.02), M0 (P=0.043), and Grade 1&2 (P=0.02), indicating its potential as a sensitive indicator of early disease progression (Figure 3A-3D).

Figure 3 Survival analyses of VDAC1 and VDAC2 in LIHC. (A) Correlations between VDAC1 expression levels and patient survival in LIHC. (B) Kaplan-Meier survival curves and boxplots showing VDAC1 expression in LIHC patients at different stages (T1&T2, T3&T4, N0, M0, Grade 1&2). (C) Correlations between VDAC1 expression levels and patient survival in LIHC. (D) Kaplan-Meier survival curves and boxplots showing VDAC1 expression in LIHC patients at different stages. CI, confidence interval; HR, hazard ratio; LIHC, liver hepatocellular carcinoma; M, metastasis; N, node; T, tumor.

Distinct genomic and epigenetic alterations underlie VDAC1/2 expression and prognosis in LIHC

Using the cBioPortal platform, we characterized the somatic alterations of VDAC1 and VDAC2 across TCGA cohorts (Figure 4A-4D). Both genes exhibited distinct mutation spectra across cancer types. VDAC1 was primarily affected by copy number amplification and deletion, whereas VDAC2 displayed amplifications as well as a diverse range of mutation types, including missense and truncating mutations. The highest alteration frequencies of VDAC1 occurred in KIRC and CHOL, while VDAC2 showed a relatively high mutation rate in UCEC and ESCA. CNVs were frequently observed as heterozygous gains and losses, with VDAC1 amplification being particularly pronounced in KIRC and CHOL, while VDAC2 alterations were more widespread, especially in LIHC, THCA, and KIRC.

Figure 4 Genetic alteration analysis in pan-cancer. (A,B) Genetic alteration in VDAC1 and VDAC2 in pan-cancer tissues. (C,D) The alteration frequency with the mutation type of VDAC1 and VDAC2 in different cancers. Structural variant¹ indicates structural variants are displayed preferentially over copy number alterations when both alterations are present in the same sample. VUS, variants of uncertain significance.

Integration of CNV and methylation profiles further revealed regulatory mechanisms underlying VDAC1/2 expression (Figure S2) (Figure 5A,5B). Structural variations, amplifications, and deep deletions constituted the predominant alteration types of both genes across cancers. In LIHC, VDAC1 was mainly characterized by copy number amplification, whereas VDAC2 tended to undergo copy number loss. Importantly, VDAC1 amplification was significantly correlated with poorer OS (P<0.001; Figure 5C), while CNV of VDAC2 showed no significant prognostic impact (P=0.26; Figure 5C).

Figure 5 Genetic and epigenetic features of VDAC1 and VDAC2 in pan-cancer. (A) Heterozygous CNVs of VDAC1 and VDAC2 across cancer types (data from GSCA). (B) CNV profiles of VDAC1 and VDAC2 in each cancer type, including heterozygous and homozygous alterations. (C) Survival differences between CNV and WT VDAC1 and VDAC2 in LIHC (data from GSCA). (D) Statistically significant differences exist in the promoter methylation levels of VDAC1 and VDAC2 between normal tissues and tumor tissues in LIHC. (E) Survival differences between high and low methylation of VDAC1 and VDAC2 in LIHC (data from GSCA). CNVs, copy number variations; GSCA, Gene Set Cancer Analysis; LIHC, liver hepatocellular carcinoma; OS, overall survival; TCGA, The Cancer Genome Atlas.

The promoter methylation level of VDAC1 was significantly reduced in primary LIHC tissues compared with normal liver tissues (P<0.001; Figure 5D), indicating a hypomethylation status in tumors.

Kaplan-Meier survival analysis further demonstrated that patients with lower VDAC1 promoter methylation exhibited significantly poorer OS compared to those with higher methylation levels (log-rank P=0.02; Figure 5E). By contrast, methylation changes in VDAC2 were not statistically significant (P=0.29; Figure 5D,5E). Collectively, these results indicate that VDAC1 and VDAC2 undergo distinct genomic and epigenetic regulation across cancers. In particular, VDAC1 exhibits prominent alterations—especially copy number amplification and promoter hypomethylation—which are strongly associated with poor prognosis in LIHC.

VDAC1 and VDAC2 exhibit differential immune regulatory roles with potential immunotherapeutic value in LIHC

We conducted an extensive immunogenomic analysis using TCGA datasets to evaluate the immunological roles of VDAC1 and VDAC2 in LIHC (Figure S3). This analysis revealed significant cancer-type-specific heterogeneity in the immune regulatory functions of these genes, prompting a focused investigation on LIHC. Our findings highlighted distinct immune infiltration patterns associated with VDAC expression levels. Specifically, high VDAC1 expression was linked to increased infiltration of immunosuppressive cells such as M0 macrophages, regulatory T cells (Tregs) (9), and resting dendritic cells, while concurrent reductions in M1 macrophages, mast cells, and CD4⁺ memory T cells were observed (Figure 6A) (10). These results suggest that VDAC1 may contribute to immune escape mechanisms in LIHC. Conversely, high VDAC2 expression was associated with increased infiltration of M0 macrophages and CD4⁺ memory T cells, but reduced levels of M1 macrophages, γδ T cells, and B cells (Figure 6B) (11). Notably, the depletion of γδ T cells, which play a critical role in cytotoxic immune responses, suggests that VDAC2 overexpression may impair antitumor immunity (11).

Figure 6 Correlation of VDAC1 and VDAC2 expression with immune cell infiltration in LIHC. Boxplots showing the differential infiltration levels of immune cell subtypes between VDAC1 (A) and VDAC2 (B) high and low expression groups in LIHC. (C) Lollipop plots illustrating the correlation between VDAC1 and VDAC2 expression and 24 immune cell types in LIHC. Kaplan-Meier survival analysis of LIHC patients stratified by VDAC1 (D) and VDAC2 (E) expression and immune cell infiltration levels. (F) TIDE scores comparing high and low expression groups of VDAC1 and VDAC2 in LIHC. (G) Heatmap showing the correlation between VDAC1 and VDAC2 expression and immune checkpoint gene expression across multiple cancer types. *, P<0.05; **, P<0.01; ***, P<0.001. HR, hazard ratio; LIHC, liver hepatocellular carcinoma; TIDE, Tumor Immune Dysfunction and Exclusion.

Survival analysis integrating VDAC expression with immune cell infiltration levels revealed that immune-dependent prognostic effects were more pronounced in LIHC (Figure 6C-6E). Low VDAC1 expression combined with high CD4⁺ T cell infiltration was associated with significantly improved survival outcomes, while high VDAC2 expression and elevated CD4⁺ T cell infiltration also predicted a favorable prognosis. On the other hand, high M0 macrophage infiltration in the low-expression contexts of both VDAC1 and VDAC2 was linked to worse survival outcomes, indicating that downregulation of VDACs may exacerbate the negative impact of M0 macrophages in immunosuppressive environments.

To assess the potential for immunotherapy, we performed TIDE analysis, which revealed that both high VDAC1 and VDAC2 expression in LIHC were associated with elevated TIDE scores, reflecting increased immune escape and reduced sensitivity to immunotherapy (Figure 6F). Additionally, immune checkpoint analysis showed that VDAC1 expression correlated positively with key immune checkpoint molecules, such as PD-1 and CTLA-4, suggesting that VDAC1 may facilitate immunosuppressive signaling (Figure S4) (Figure 6G). Similarly, VDAC2 expression showed strong positive correlations with HAVCR2, LAG3, PDCD1, and IT-PRIPL1, which were upregulated in the high-expression group, further supporting the involvement of VDAC2 in immune evasion (Figure 6G).

VDAC1 and VDAC2 regulate proliferation, apoptosis, ferroptosis, and migration in HCC cells

To explore the biological roles of VDAC1 and VDAC2 in LIHC, we selected Hep3B and PLC cell lines based on their cancer-specific expression and prognostic patterns identified through TCGA multi-omics analyses. siRNA-mediated knockdown of VDAC1 and VDAC2 in both cell lines was confirmed by qRT-PCR and Western blotting. Functional assays revealed that depletion of eithergene significantly reduced colony formation in both Hep3B and PLC cells (Figure 7A,7B), with CCK-8 assays showing that VDAC1 silencing notably suppressed proliferation (Figure 7C). Apoptosis assays demonstrated differential effects: VDAC1 or VDAC2 knockdown significantly increased apoptosis in Hep3B cells, whereas in PLC cells, VDAC1 depletion had only a modest effect, but VDAC2 knockdown markedly elevated apoptosis, indicating a stronger anti-apoptotic role for VDAC2 in this context (Figure 7C,7D). We further shift our focus to the mechanism of ferroptosis—a form of iron-dependent, lipid peroxidation-driven cell death closely associated with mitochondrial oxidative stress and the VDAC channel (12,13). In the PLC cell line, knockdown of VDAC1 significantly reduces intracellular ROS levels, whereas no significant change is observed in the Hep3B cell line (Figure 7E,7F). Meanwhile, knockdown of VDAC2 increases ROS levels in both cell lines, which Fer-1, a typical ferroptosis inhibitor, can partially or fully reverse (Figure 7E,7F).

Figure 7 Effects of VDAC1 and VDAC2 on proliferation, apoptosis, and ROS levels in hepatocellular carcinoma cell lines. We selected hepatocellular carcinoma cell lines (Hep3B and PLC) and performed siRNA-mediated knockdown of VDAC1 and VDAC2, with a non-targeting siRNA (NC) as a control. (A) Representative images of colony formation assays stained with 0.1% crystal violet for 30 min (original magnification, ×1). (B) Cell viability was assessed using CCK-8 assays at the indicated time points. (C,D) Apoptosis analysis by flow cytometry using Annexin V-FITC/PI staining and quantification of apoptotic cell percentages. (E,F) ROS levels were measured by flow cytometry after staining with DCFDA, and the quantification of ROS-positive cells is shown. DCFDA, 2'-7'-dichlorodihydrofluorescein diacetate.

To examine the roles of VDAC1 and VDAC2 in tumor cell migration, Transwell assays were conducted in both cell lines, revealing that silencing eithergene significantly inhibited migration (Figure 8A). To further elucidate the underlying mechanism, RT-qPCR analysis of EMT-associated markers was conducted (Figure 8B,8C). Knockdown of VDAC1 or VDAC2 led to a marked increase in the epithelial marker E-cadherin, accompanied by a reduction in mesenchymal markers (N-cadherin, Vimentin) and the EMT-inducing transcription factor ZEB1, suggesting that loss of either gene suppresses EMT (14). These results were confirmed at the protein level through Western blotting, showing a robust reduction of N-cadherin along with efficient knockdown of VDAC1 and VDAC2 (Figure 8D). These findings collectively suggest that both VDAC1 and VDAC2 regulate cell proliferation, apoptosis, and migration in LIHC cells, with potential roles in EMT activation and tumor invasion (Figure 8D).

Figure 8 Effects of VDAC1 and VDAC2 on cell migration and EMT process in hepatocellular carcinoma cell lines. (A) Transwell migration assays of Hep3B and PLC cells after transfection with si-NC, si-VDAC1, or si-VDAC2. (B,C) Relative mRNA expression of EMT markers (e.g., E-cadherin, N-cadherin, Vi-mentin) detected by qPCR. (D) Protein expression of EMT-related markers was analyzed by Western blotting. EMT, epithelial-mesenchymal transition; NC, negative control; qPCR, quantitative polymerase chain reaction.

Discussion

Key findings

This study systematically evaluated the expression patterns, stage-specific prognostic value, and potential biological functions of VDAC1 and VDAC2 in LIHC. Both proteins were found to be significantly upregulated in LIHC and independently predict patient survival. VDAC1 has a stronger prognostic impact in advanced tumors, whereas VDAC2 remains predictive in early stages. Multi-omics analysis revealed that VDAC1 upregulation is primarily driven by copy number amplification and promoter hypomethylation, while VDAC2 is more influenced by post-transcriptional regulation and the tumor microenvironment. These differences explain their complementary roles across tumor stages. Further pathway enrichment analyses highlighted that both VDAC1 and VDAC2 are closely associated with ferroptosis and EMT pathways, suggesting that mitochondrial metabolism and redox homeostasis are crucial for tumor progression and immune regulation. VDAC1 is more involved in metabolic activation linked to EMT, whereas VDAC2 regulates ferroptosis and oxidative stress responses, demonstrating their complementary functions in different stages of tumor progression.

Strengths and limitations

This study integrates multi-omics analyses, immune profiling, and functional validation, providing a comprehensive evaluation of VDAC1 and VDAC2 in LIHC. The use of TCGA data strengthens the generalizability of the findings, particularly in elucidating the roles of these proteins in prognosis and therapy. However, the study relies primarily on retrospective data, with limited representation of stage IV HCC, which may impact its applicability to advanced disease stages. While multi-omics data offer valuable insights, further in vivo and clinical validation is needed. Additionally, the study does not account for tumor heterogeneity across different etiologies (e.g., viral or alcohol-related HCC), which could influence VDAC expression and function. Prospective studies are essential to confirm the clinical relevance and therapeutic potential of VDAC1 and VDAC2 in liver cancer.

Comparison with similar research

Previous studies have reported the high expression of VDAC1 and its association with prognosis in various cancers, while the biological functions of VDAC2 display clear context-dependence across different tumor types (15-18). This study further expands these observations in HCC, revealing their complementary prognostic value in different stages and their relationship with immune infiltration. Notably, a recent high-impact study has indicated that both VDAC1 and VDAC2 are involved in mitochondrial-mediated cell fate processes (18). VDAC1 has been found to affect sensitivity to ferroptosis inducers (e.g., erastin) and plays a dual role in energy metabolism and lipid peroxidation balance (19,20); VDAC2 plays a key role in inflammatory signaling and ferroptosis regulation, linking mitochondrial function with oxidative stress and immune responses (21,22). Our results systematically validate and extend these mechanisms in liver cancer, particularly emphasizing the complementary functional features of both in ferroptosis and EMT pathways. This provides new evidence for understanding the role of the VDAC family in tumor metabolism, redox homeostasis, and immune regulation.

Explanations of findings

The prognostic role of VDAC1 is primarily observed in advanced tumors, likely due to its high expression driven by copy number amplification and promoter hypomethylation. It regulates mitochondrial membrane permeability and metabolite exchange, promoting energy metabolism and lipid peroxidation, thus enhancing tumor growth. VDAC1 may also increase sensitivity to ferroptosis under stress, reflecting its context-dependent dual functions. In contrast, VDAC2 maintains mitochondrial homeostasis, inhibits excessive ferroptosis, and regulates inflammation. Knockdown of VDAC2 in PLC and Hep3B cells increased ROS levels, which were partially reversed by the ferroptosis inhibitor Fer-1, suggesting VDAC2 suppresses ferroptosis-related oxidative stress. Knockdown of VDAC1 reduced ROS in PLC cells but had no effect in Hep3B cells, highlighting differences in redox regulation across cell types.

The ferroptosis-inhibiting effect of VDAC2 aligns with its association with poor prognosis. By stabilizing mitochondrial function, reducing ROS, and suppressing ferroptosis, VDAC2 promotes tumor cell survival and immune evasion. This explains the link between high VDAC2 expression and poor prognosis, supporting its role in early-stage tumorigenesis. Overall, VDAC1 drives tumor progression and invasion, while VDAC2 regulates ferroptosis and immune homeostasis, with complementary roles in different cancer stages.

Implications and actions needed

The complementary features of VDAC1 and VDAC2 make them potential stage-specific biomarkers and therapeutic targets in HCC. The metabolic and epigenetic activation of VDAC1 positions it as a target for metabolic intervention and energy-blocking therapies. Its association with EMT suggests that combining VDAC1-targeted therapies with anti-migration or anti-invasion drugs could improve prognosis. The role of VDAC2 in ferroptosis regulation opens new directions for ferroptosis-related therapies. Reducing VDAC2 expression or inhibiting its protective functions may enhance tumor cell sensitivity to ferroptosis inducers and immune checkpoint inhibitors. Future studies should focus on elucidating the specific roles of VDAC1/2 in ferroptosis and EMT pathways, exploring their context-dependent effects in different molecular subtypes and therapeutic settings, and evaluating their feasibility as predictive or therapeutic biomarkers through prospective clinical trials.

Conclusions

VDAC1 and VDAC2 exhibit complementary prognostic and functional features in HCC. VDAC1 is driven by copy number amplification and promoter hypomethylation, promoting tumor progression through metabolic reprogramming and EMT activation. VDAC2 stabilizes mitochondrial function, reduces ROS, and inhibits ferroptosis signaling, protecting tumor cells from oxidative damage and promoting immune evasion. Together, they form a critical axis linking mitochondrial metabolism, oxidative stress, and immune regulation, providing a novel theoretical foundation and potential targets for stage-specific risk assessment and precision therapy in HCC.

Acknowledgments

The authors would like to express their appreciation to GEO and TCGA for providing open access to their valuable datasets, which significantly contributed to the success of this study.

Footnote

Reporting Checklist: The authors have completed the TRIPOD and MDAR reporting checklists. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1887/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1887/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1887/prf

Funding: This study was funded by the foundation of Chengdu Municipal Science and Technology Bureau (No. 2024-YF05-00853-SN), the Third People’s Hospital of Chengdu Scientific Research Project (No. 2023PI19), and the Medical Research Project of Chengdu Municipal Health Commission (No. 2025455).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1887/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Laversanne M, Sung H, 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:229-63. [Crossref] [PubMed]
2. Du H, Xu T, Yu S, et al. Mitochondrial metabolism and cancer therapeutic innovation. Signal Transduct Target Ther 2025;10:245. [Crossref] [PubMed]
3. Shoshan-Barmatz V, De S, Meir A. The Mitochondrial Voltage-Dependent Anion Channel 1, Ca(2+) Transport, Apoptosis, and Their Regulation. Front Oncol 2017;7:60. [Crossref] [PubMed]
4. Shoshan-Barmatz V, Maldonado EN, Krelin Y. VDAC1 at the crossroads of cell metabolism, apoptosis and cell stress. Cell Stress 2017;1:11-36. [Crossref] [PubMed]
5. Shoshan-Barmatz V, Shteinfer-Kuzmine A, Verma A. VDAC1 at the Intersection of Cell Metabolism, Apoptosis, and Diseases. Biomolecules 2020;10:1485. [Crossref] [PubMed]
6. Bao M, Wu A. Understanding the heterogeneity in liver hepatocellular carcinoma with a special focus on malignant cell through single-cell analysis. Discov Oncol 2024;15:241. [Crossref] [PubMed]
7. Binnewies M, Roberts EW, Kersten K, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med 2018;24:541-50. [Crossref] [PubMed]
8. Passaniti A, Kim MS, Polster BM, et al. Targeting mitochondrial metabolism for metastatic cancer therapy. Mol Carcinog 2022;61:827-38. [Crossref] [PubMed]
9. Zong Y, Deng K, Chong WP. Regulation of Treg cells by cytokine signaling and co-stimulatory molecules. Front Immunol 2024;15:1387975. [Crossref] [PubMed]
10. Boutilier AJ, Elsawa SF. Macrophage Polarization States in the Tumor Microenvironment. Int J Mol Sci 2021;22:6995. [Crossref] [PubMed]
11. Liu Y, Zhang C. The Role of Human γδ T Cells in Anti-Tumor Immunity and Their Potential for Cancer Immunotherapy. Cells 2020;9:1206. [Crossref] [PubMed]
12. Jin X, Tang J, Qiu X, et al. Ferroptosis: Emerging mechanisms, biological function, and therapeutic potential in cancer and inflammation. Cell Death Discov 2024;10:45. [Crossref] [PubMed]
13. Oh SJ, Ikeda M, Ide T, et al. Mitochondrial event as an ultimate step in ferroptosis. Cell Death Discov 2022;8:414. [Crossref] [PubMed]
14. Kalluri R, Weinberg RA. The basics of epithelial-mesenchymal transition. J Clin Invest 2009;119:1420-8. [Crossref] [PubMed]
15. Fang Y, Liu J, Zhang Q, et al. Overexpressed VDAC1 in breast cancer as a novel prognostic biomarker and correlates with immune infiltrates. World J Surg Oncol 2022;20:211. [Crossref] [PubMed]
16. Magrì A, Reina S, De Pinto V. VDAC1 as Pharmacological Target in Cancer and Neurodegeneration: Focus on Its Role in Apoptosis. Front Chem 2018;6:108. [Crossref] [PubMed]
17. Zhou K, Yao YL, He ZC, et al. VDAC2 interacts with PFKP to regulate glucose metabolism and phenotypic reprogramming of glioma stem cells. Cell Death Dis 2018;9:988. [Crossref] [PubMed]
18. Yuan S, Sun R, Shi H, et al. VDAC2 loss elicits tumour destruction and inflammation for cancer therapy. Nature 2025;640:1062-71. [Crossref] [PubMed]
19. Garcia-Gimenez A, Ditcham JE, Azazi DMA, et al. CREBBP inactivation sensitizes B cell acute lymphoblastic leukemia to ferroptotic cell death upon BCL2 inhibition. Nat Commun 2025;16:4274. [Crossref] [PubMed]
20. Yang J, Lu X, Hao JL, et al. VSTM2L protects prostate cancer cells against ferroptosis via inhibiting VDAC1 oligomerization and maintaining mitochondria homeostasis. Nat Commun 2025;16:1160. [Crossref] [PubMed]
21. Luo P, Zhang Q, Shen S, et al. Mechanistic engineering of celastrol liposomes induces ferroptosis and apoptosis by directly targeting VDAC2 in hepatocellular carcinoma. Asian J Pharm Sci 2023;18:100874. [Crossref] [PubMed]
22. She H, Tan L, Du Y, et al. VDAC2 malonylation participates in sepsis-induced myocardial dysfunction via mitochondrial-related ferroptosis. Int J Biol Sci 2023;19:3143-58. [Crossref] [PubMed]