Hypoxia-related signatures predicts survival, immunosuppression and PARP inhibitor resistance in HCC

Abstract

Background: Despite extensive research on hypoxia in hepatocellular carcinoma (HCC), previous studies have relied on pre-existing hypoxia gene sets, limiting their specificity. We developed a novel approach using direct comparison of hypoxic versus normoxic HCC cell lines to establish a more precise hypoxia signature.

Methods: Through differential gene expression analysis of HCC cells under controlled oxygen conditions (GSE185969) and integration with TCGA-LIHC data, we identified and validated a highly specific 29-gene hypoxia signature. We performed comprehensive immune profiling and genomic instability analyses using multi-omics approaches.

Results: Our HCC-specific hypoxia signature demonstrated superior prognostic value (AUC: 0.805, 0.805, 0.748 at 1/3/5 years) compared to conventional hypoxia markers. High-risk tumors showed distinct immunosuppressive features including reduced CD8 + T cells and elevated Th2 cells, along with significantly increased expression of immune checkpoints CD274 (PD-1, p < 0.05) and CD276 (B7-H3, r = 0.62, p < 0.001). Notably, we uncovered an unexpected inverse relationship between hypoxia-induced genomic instability and PARP inhibitor sensitivity, challenging current therapeutic paradigms.

Conclusion: Our methodology establishes a more precise hypoxia signature specific to HCC, advancing beyond traditional approaches. The paradoxical finding of reduced PARP inhibitor sensitivity in genomically unstable tumors reveals new complexities in hypoxia-driven treatment resistance, suggesting the need for alternative therapeutic strategies in hypoxic HCC.

Keywords: Genomic instability; Hepatocellular carcinoma; Hypoxia; Prognostic signature; Tumor immune microenvironment.

Fig. 1

Transcriptional and functional analysis of hypoxia-associated genes in hepatocellular carcinoma. A Volcano plot of differentially expressed genes (DEGs). Log2 fold change (log2FC) is plotted against – log10 (P-value). Significantly downregulated genes (blue dots) and upregulated genes (red dots) are shown. Dot size represents statistical significance. Key genes of interest are indicated by arrows. B Ridge plot visualization of the top 10 enriched Gene Ontology (GO) terms. Functional enrichment patterns are represented by Normalized Enrichment Scores (NES), with positive enrichment in red and negative in blue. C Ridge plot depicting the top 10 enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Color scheme follows (B), with NES values indicated by color intensity. D Venn diagram showing the overlap between prognostic genes (left, n = 68) identified through univariate Cox regression analysis and differentially expressed genes (right, n = 528). E Enrichment analysis of prognostic genes (n = 68) presented as a bar plot, with significance levels indicated by p-values. F Functional network analysis of enriched terms clustered by biological similarity. Nodes with identical cluster IDs are spatially grouped. G Network visualization of enriched terms with p-value overlay. Node size correlates with gene count, and color intensity reflects statistical significance

Fig. 2

LASSO-cox model analysis and risk score validation. A Cross-validation plot for LASSO coefficient selection. The x-axis shows log(lambda) values and y-axis displays partial likelihood deviance. Vertical dashed lines represent optimal lambda values. B Coefficient profile plot showing the path of each gene’s coefficient against log(lambda). Each curve represents an individual gene’s coefficient trajectory, demonstrating the regularization effect on feature selection. C Risk score distribution plot depicting patient stratification. Patients are arranged from low to high-risk scores, with blue and red dots representing low-risk and high-risk patients respectively. The dotted line indicates the median risk score threshold used for group stratification. D Patient survival status distribution. Blue dots indicate surviving patients while red dots represent deceased patients, arranged according to risk score ranking. E Kaplan-Meier survival analysis comparing high-risk (red) and low-risk (blue) groups. Survival curves demonstrate significant prognostic stratification between risk groups, with corresponding hazard ratio (HR) and 95% confidence interval (CI). F Time-dependent ROC curves showing predictive accuracy at 1-year (green), 3-year (blue), and 5-year (red) time points for the LASSO-Cox model

Fig. 3

Clinical validation of the prognostic model. A Nomogram incorporating clinical variables (age, gender, tumor staging [T, N, M stages], and risk score) for predicting 1-, 3-, and 5-year survival probabilities in cancer patients. Line length represents the relative contribution of each variable to survival outcomes. B–D Calibration curves for survival prediction at 1 year (B), 3 years (C), and 5 years (D). X-axis shows predicted survival probability from the nomogram; y-axis shows observed survival probability adjusted by Kaplan-Meier estimation. Dashed diagonal line represents ideal prediction

Fig. 4

Analysis of immune cell infiltration and TIDE features between risk groups. A Heatmap visualization of differential immune cell infiltration patterns between high- and low-risk groups, showing significantly different immune cell types (p < 0.05). Immune scores are displayed as top annotations. B Box plots comparing immune-related metrics (Immune Score, Stromal Score, and Microenvironment Score) between risk groups. C Correlation plot showing the relationship between significant immune cell infiltration levels and risk scores. D Violin plots depicting the distribution of MDSC Exclusion scores in high-risk (red) and low-risk (blue) groups, with embedded box plots showing quartile ranges and median values. MDSC demonstrated the most significant inter-group difference (p = 1.65e- − 04). E Correlation matrix illustrating the relationships between ICB-related gene expression and risk scores, where color intensity represents correlation strength (red: positive correlations; blue: negative correlations)

Fig. 5

Analysis of genomic features and DNA damage response characteristics between risk groups. A–F Distribution plots showing differences in key genomic and DNA damage response features between high- and low-risk groups: A tp53_score indicating TP53 pathway dysfunction, B PARPi7 score reflecting PARP inhibitor sensitivity, C eCARD score representing targeted therapy response potential, D CNA_frac_altered showing chromosomal instability, E CNA_n_segs representing genomic segment alterations, and F CNA_n_focal_amp_del indicating focal amplification and deletion events. Asterisks denote statistical significance (*P < 0.05, **P < 0.01, ***P < 0.001). G Bar plot illustrating the distribution of PARPi7-positive and PARPi7-negative cases across high-risk and low-risk groups, with percentages indicated for each category