Abstract Background Hepatocellular carcinoma (HCC), a prevalent and highly lethal malignancy, is notorious for its aggressive nature and inherent tendency to metastasize, posing significant challenges in clinical management and prognosis. Hypoxia, a pivotal characteristic of the tumor microenvironment (TME) in hepatocellular HCC, is intimately linked to disease progression and unfavorable patient outcomes, underscoring its critical role in shaping the malignant behavior of this cancer. Methods Our research leveraged single-cell RNA sequencing technology to dissect the heterogeneity of the HCC TME, focusing on hypoxia-related genes. To probe the effects of hypoxia on HCC invasion, we developed the Cell Hypoxia-Related Prognostic Feature (CHPF). We analyzed transcriptome data from the The Cancer Genome Atlas (TCGA) database and the [28]GSE149614 dataset, employing computational methods such as UMAP, Weighted correlation network analysis (WGCNA), and CellChat to identify hypoxia cells, characterize cell subsets, and elucidate intercellular communications. Results Our analysis revealed significant heterogeneity in hypoxia cell populations within the HCC TME, with distinct expression patterns of hypoxia-related genes in neoplastic and immune cells. Our analysis revealed distinct hypoxia subpopulations within HCC, with significant overexpression of genes like MEG3, KLF6 and JUN in hypoxia cells. We identified a unique hypoxia subpopulation with high invasive potential and constructed a prognostic model based on H2-specific transcription factors including LRP10、MED8、NOL10、NOP58 and REXO4. The model demonstrated significant predictive value for lifespan of patients as verified in the TCGA dataset and an external validation group. Conclusion Key transcription factors like NOP58, MED8 play pivotal roles in hypoxia-induced HCC invasion and metastasis, and a predictive model based on these factors forecasts HCC survival. Our findings provide novel molecular markers and therapeutic targets for HCC, highlighting the importance of considering the hypoxia TME in diagnostic and treatment strategies. Supplementary Information The online version contains supplementary material available at 10.1007/s12672-025-03201-y. Keywords: Hepatocellular carcinoma, Hypoxia, Single-cell RNA sequencing, Tumor microenvironment, Prognostic model Background Globally, HCC stands as a prevalent and formidable malignancy, posing a significant health challenge. According to the latest epidemiological surveys, primary liver cancer ranks second in mortality in China, posing a grave threat to human life and health due to its stealthy nature, proneness to metastasis, and aggressive malignancy [[29]1, [30]2]. Despite surgery being the primary therapeutic modality for primary HCC, postoperative recurrence and metastasis remain significant clinical challenges, significantly impacting patient prognosis. From 2005 to 2020, the annual mortality rate of HCC decreased by only 1–3% [[31]1], boasting a 5-year survival rate nearing 10% [[32]3]. Remarkably, more than 90% of postoperative fatalities among HCC patients can be traced back to tumor recurrence and metastasis, highlighting their critical role in compromising patient outcomes. The overwhelming majority of liver cancer patients receive their diagnosis at an advanced stage of the disease, often accompanied by distant metastases to multiple organs such as the lungs, bones, and adrenal glands [[33]4–[34]6]. The “premetastatic niche” (PMN) concept explains the immunosuppressive changes in tissues that prepare for tumor metastasis [[35]7]. Exosomes from the tumor microenvironment (TME) can influence PMN formation, aiding in tumor cell colonization and immune evasion [[36]8]. Hypoxia is a crucial negative element influencing the outlook for individuals with liver cancer. Studies show that the PaO[2] level in normal liver tissue is approximately 33 mmHg, whereas it significantly decreases in HCC tissues, averaging only 6 mmHg [[37]9]. The imbalance between insufficient oxygen supply and high oxygen consumption by neoplastic cells is pivotal in inducing hypoxia [[38]10]. Furthermore, liver cancer treatments such as TACE and TKI may also exacerbate hypoxia, thereby facilitating lung metastasis [[39]11–[40]13]. Hypoxia-inducible factor 1 (HIF-1), particularly its HIF-1α subunit, is crucial in mediating the response to low oxygen levels, promoting tumor growth, metastasis, and resistance to therapies [[41]14, [42]15]. The role of HIF-1α in HCC prognosis is complex, with some studies linking its overexpression to poor outcomes [[43]16–[44]20], while others find no significant association in non-cirrhotic patients [[45]21]. HIF-1α can enhance tumor cell proliferation [[46]22], promote EMT [[47]23], and increase angiogenesis through the regulation of VEGF [[48]24]. Since the induction of oxidative stress is a crucial mechanism for the antitumor effects of TKIs [[49]25], metabolic reprogramming under hypoxia contributes to drug resistance by reducing ROS levels in TKI-treated HCC [[50]26]. Furthermore, alterations in metabolite levels, including glucose [[51]27], lactate [[52]28], and adenosine [[53]29] within the TME, collectively foster an immunosuppressive milieu that significantly hinders the efficacy of ICIs in HCC. Hypoxia triggers a range of metabolic changes through HIFs, which enhance drug resistance and result in poor treatment outcomes for patients with HCC. Studying the hypoxia microenvironment in HCC faces challenges: (1) The mechanisms of how TAMs respond to hypoxia are not well understood, despite knowing their role in immunosuppression and resistance to anti-PD-L1 therapy [[54]30]. (2) HCC involves many mutations, with 30–40 per cell and 5–8 being driver genes, making targeted treatments comple [[55]31, [56]32]. (3) The diverse immune microenvironment in liver cancer contributes to drug resistance and poor outcomes but is not well-studied [[57]33]. (4) Limited efficacy of FDA-approved drugs and ineligibility for surgery in most patients [[58]34–[59]36], highlighting limitations in current treatment strategies. Single-cell RNA sequencing has advanced our understanding of HCC’s tumor heterogeneity and immune microenvironment, identifying potential biomarkers and therapeutic targets [[60]37–[61]40]. Within this examination, we harness the liver malignancy genomic data from the TCGA archive, combined with the [62]GSE149614 dataset, to systematically investigate the phenotypic expression and its correlation with clinical outcomes of hypoxia-related genes in HCC. Our objectives are to pinpoint genes related to hypoxia that are linked to the prognosis of HCC, explore their mechanisms within the hypoxia TME, and provide novel molecular indicators and therapeutic objectives for HCC detection and management. Methods Acquisition and processing of transcriptomic data In this study, RNA molecular profiles and medical records from HCC patients (n = 368) were retrieved from TCGA database for model construction [[63]41]. To validate the model’s stability and accuracy, the LIRI-JP dataset from the ICGC database, comprising microarray data from 232 liver cancer samples, was employed as a validation set. All data underwent Transcripts Per Million (TPM) format conversion followed by log2 transformation to facilitate subsequent analysis. Acquisition and processing of single-cell RNA sequencing data The single-cell dataset used in this study was sourced from [64]GSE149614 in the Gene Expression Omnibus (GEO) database [[65]42], encompassing six liver cancer samples: three primary tumors (HCC07T, HCC08T, HCC10T), two portal vein tumor thrombi (PVTT), and one metastatic lymph node (MLN). We conducted the data analysis utilizing the Seurat package in the R programming environment, specifically R version 4.1.3. Quality control criteria included mitochondrial content below 20%, hemocyte content below 3%, and UMI counts and gene counts ranging from 200−50,000 and 200–8000, respectively. The data underwent normalization, followed by the selection of the top 2,000 most variable genes. Subsequently, we performed data transformation utilizing the NormalizeData, FindVariableFeatures, and ScaleData functions from the Seurat suite, aiming to reduce the effects of cell cycle variability on our analysis. Harmony algorithm was applied for batch effect correction. We employed UMAP and t-SNE algorithms to reduce the dimensionality of the data, while clustering analysis was performed using the Louvain algorithm. Differential gene analysis employed the FindAllMarkers function, using selection criteria established at a p-value of less than 0.05, log2 fold change greater than 0.25, and an expression proportion exceeding 0.1. Identification of hypoxia cells An open-source software, CHPF (available at github.com/yihan1221/CHPF), was utilized in this study to predict cellular hypoxia conditions by integrating single-cell transcriptomic profiles with low-oxygen-induced gene clusters. Cell annotation analysis Cell annotation was performed using a panel of cell markers, including neoplastic cells (CDH1, EPCAM, KRT18, KRT19,), fibroblasts (SLRR1B, CD90, COL1A1, COL1A2), endothelial cells (CD31, CLDN2, VEGFR-1, RAMP2), T-cells (CD3D/E/G, IMD7), NK-cells (GIG1, NKG5, CD56, CD94), B-cells (CD79A, AGM1, IgG3, IGHA2), and myeloid cells (AMYLD5, SCARA2, CD16, CD68). Individual cluster analyses were performed on the cancer cells to investigate the diversity within the tumor, with results depicted through UMAP and t-SNE plots, as well as bar graphs and heatmaps. WGCNA analysis and enrichment analysis We used the WGCNA package to examine gene modules correlated with the H-group cells, and performed gene enrichment analysis with the clusterProfiler package [[66]43], taking advantage of the GOBP and KEGG databases. The enrichment outcomes for both the H-group and N-group cells were graphically represented using the EnrichmentMap and AutoAnnotate plugins within the Cytoscape platform. Analysis of intercellular communication CellChat package was used to assess potential intercellular communication [[67]44]. The gene expression matrix, once normalized, was fed into the CellChat framework to establish a CellChat entity. Preliminary data handling encompassed the application of functions such as identifyLIHCerExpressedGenes, identifyLIHCerExpressedInteraction, and ProjectData. The likelihood of ligand-receptor interactions was deciphered through the execution of computeCommunProb, filterCommunication, and computeCommunProbPathway utilities, culminating in the assembly of a network mapping cellular communications utilizing the aggregateNet function. Cell development trajectory analysis Monocle was used to construct pseudo-temporal trajectories of cells, and analysis was performed using Monocle2 (version 2.22.0) using the DDRTree dimensionality reduction method with input as a highly variant gene set in time direction according to the starting cell cluster [[68]45]. Single-cell CNV analysis Employing the infercnv package, we estimated the copy number variations within the neoplastic cells, benchmarking against endothelial cells as a comparative baseline. For every tumor cell, a Copy Number Variation (CNV) score was derived to quantify these genetic alterations. Single-cell transcription factor analysis The SCENIC package was employed to predict transcription factors in H2 and N2 cell populations, with GRNboost2 software used for gene co-expression analysis to construct gene regulatory networks [[69]46]. Important nodes in the network were assessed by degree, and the top 1% of genes or transcription factors were selected for further analysis. Construction of prognostic model We focused on the transcription factors and genes that were distinctive to the H2 cell subset. Initially, a univariate Cox regression analysis was applied to identify genes associated with patient prognosis. Subsequently, the LassoCox method was employed to develop a predictive model of patient outcomes. Using this model, we assigned risk scores to each patient, and then categorized them into either a high-risk or low-risk group, with the division made at the median risk score value. Gene Set Variation Analysis (GSVA) was used for pathway enrichment analysis, scoring KEGG and HALLMARK gene sets using GSVA R package using default parameters [[70]47]. Cell differentiation potential analysis CytoTRACE was used to assess the differentiation potential of cells and run to calculate cell differentiation status scores on a normalized expression matrix based on default parameters provided by the authors (R package CytoTRACE, version 0.3.3) [[71]48]. Statistical analysis All steps of data manipulation, statistical computations, and visualization were executed utilizing R software, specifically version 4.1.3. The correlation between two continuous variables was assessed by calculating the Pearson’s correlation coefficient. For the comparison of categorical variables, a Chi-squared test was employed, whereas for continuous variables, either the Wilcoxon rank-sum test or the t-test was selected as appropriate. The survival package in R was utilized to perform Cox regression analysis and Kaplan-Meier survival analysis. Results Single-cell atlas of HCC Cell clustering and dimensionality reduction were achieved using the UMAP algorithm, cell identities were determined by referencing markers that are distinctive to each cell type. The UMAP plot showcased distinct cell types including neoplastic cells, fibroblasts, endothelial cells, B cells, NK T cells, and myeloid cells (Fig. [72]1A). CHPF software was utilized to classify cells based on their hypoxia status, revealing that neoplastic cells, fibroblasts, endothelial cells, and myeloid cells were primarily hypoxia, whereas B cells and NK T cells were mostly normoxia (Fig. [73]1B). Further analysis of hypoxia cell proportions across different samples and tissue types demonstrated the highest hypoxia prevalence in neoplastic cells, with relatively high proportions also observed in fibroblasts, endothelial cells, and myeloid cells (Fig. [74]1C). A Sankey diagram illuminated the associations between cell type, sample type, and hypoxia status, indicating that hypoxia cells were concentrated in crucial regions of tumor metastasis, such as primary tumors, portal vein tumor thrombi, and metastatic lymph nodes (Fig. [75]1D). Heatmap analysis showcased the intensity of hypoxia-triggered marker genes in neoplastic and immune cells (Fig. [76]1E). Fig. 1. [77]Fig. 1 [78]Open in a new tab HCC single-cell map analysis. A UMAP dimensionality reduction plot displaying the distribution of various cell types within single-cell data. B Application of the CHPF software for the identification of hypoxia status within single-cell data. C A bar chart is presented to show the distribution of hypoxia cells among various samples and tissue categories. D A Sankey diagram is provided, illustrating the relationship between different cell types and their hypoxia conditions. E A heatmap is displayed, depicting the expression profiles of genes related to hypoxia in both tumor and immune cells Hypoxia analysis of immune cell subpopulations This study delved into the characteristics of immune cells in the hypoxia TME. Dimensionality reduction and clustering analysis depicted the distribution of immune cell subpopulations, distinguishing hypoxia from normoxia states. Hypoxia cells were predominantly found in macrophages, myeloid dendritic cells (mDCs), and monocytes (Fig. [79]2A–C). Sankey diagram analysis revealed that these hypoxia immune cells primarily originated from primary tumors and portal vein tumor thrombi (Fig. [80]2D). CellChat software was employed to analyze intercellular communication between immune and neoplastic cells, uncovering potential communication networks and presenting communication maps along with ligand-receptor interaction diagrams (Fig. [81]2E, F). Notably, high communication levels were observed between SPP1 and CD44, prompting further analysis of their expression patterns in the UMAP plot (Fig. [82]2G, H). Correspondingly, a heatmap of MIF signaling pathway communication among various cell types was also presented (Supplementary Fig. [83]1B). Fig. 2. [84]Fig. 2 [85]Open in a new tab Hypoxia analysis of immune cell subsets. A–C UMAP plots illustrating the clustering of immune cell subsets and their hypoxia status. D Sankey diagram analyzing the association between immune cell subsets and hypoxia status. E, F Presentation of the communication network between immune and neoplastic cells. G,H Single-cell data showing the expression patterns of SPP1 and CD44 Hypoxia analysis of neoplastic cell subpopulations Cluster analysis of neoplastic cells based on hypoxia status identified four hypoxia clusters (H1-H4) and three normoxia clusters (N1-N3) (Fig. [86]3A and B). Hypoxia scoring revealed that the H2 subpopulation exhibited the highest degree of hypoxia (Fig. [87]3C). Heatmap analysis showcased the expression of marker genes in each cluster, with significant upregulation of MEG3, KLF6, JUN, among others, in the H2 subpopulation (Fig. [88]3D). GOBP enrichment analysis and Cytoscape application illuminated functional enrichments of hypoxia-specific marker genes, involving biological processes such as blood coagulation, fibrin clot formation, neutrophil migration and chemotaxis, and cadmium ion and copper ion homeostasis (Fig. [89]3E). WGCNA analysis identified key gene modules associated with each subpopulation, such as the high correlation of the H2 subpopulation with the MEred module and the H3 subpopulation with the MEyellow module (Fig. [90]3F). GOBP enrichment analysis of these gene modules revealed processes including cytoplasmic translation, protein ubiquitination, T-cell activation, leukocyte cell-cell adhesion, and others (Fig. [91]3G). CytoTRACE analysis demonstrated that the H2 subset exhibited the most pronounced differentiation metric (Fig. [92]3H). Pseudo-temporal analysis using monocle3 software showed that neoplastic cells with severe hypoxia were located in the middle and late stages of the pseudo-temporal trajectory (Fig. [93]3I). In addition, we presented the distribution of tumor cell subsets, hypoxia scores, and hypoxia-specific subset scores by UMAP plots and boxplots (Supplementary Fig. [94]2A–C). Fig. 3. [95]Fig. 3 [96]Open in a new tab Hypoxia and functional analysis of tumor cell subsets. A–C UMAP plots and ridge plots depicting the distribution of hypoxia and normoxia tumor cell subsets and their hypoxia scores. D A heatmap is utilized to illustrate the expression levels of signature genes across various cellular subpopulations. E GOBP enrichment analysis revealing the functional characteristics of hypoxia subsets. F WGCNA analysis demonstrating the correlation between gene modules and cell subsets. G GOBP enrichment analysis of specific gene modules, left panel is MEred module gene enrichment results for H2 subset and right panel is MEyellow module gene enrichment results for H3 subset. H CytoTRACE analysis assessing the differentiation degree of cell subsets. I Pseudotime analysis using Monocle3 software revealing the developmental trajectory of cell subsets Analysis of tumor-related pathways, copy number variations, and transcription factors To gain insights into the relationship between hypoxia and HCC aggressiveness, signature genes related to oxygen deprivation, infiltration, cell death, blood vessel formation, and epithelial transformation. were retrieved from the CancerSEA database, and activity scores for each cell subpopulation were calculated using GSVA (Supplementary Fig. [97]3). Correlation analysis revealed positive correlations between hypoxia scores and the other four signature scores (Fig. [98]4A). Invasion analysis showed that H2 subpopulation displayed the utmost invasive potential, whereas N2 subpopulation obtained the lowest (Fig. [99]4B). CNV analysis revealed CNV profiles of each tumor subpopulation (Supplementary Fig. [100]4), and CNV conditions of neoplastic cells was assessed utilizing endothelial cells served as a benchmark (Fig. [101]4C). Transcription factor analysis further identified key genes and transcription factors in the H2 and N2 subpopulations, extracting the top 1% of associated nodes and presenting Venn diagrams (Fig. [102]4D, E). A cumulative total of 119 key genes were pinpointed within the H2 subset, while 70 were discovered in the N2 subset. Correlation analysis of H2-specific transcription factors and genes with Hallmark pathways revealed a heatmap of positive correlations with markers and P-values < 0.05 (Fig. [103]4F). The gene sets from the MSigDB hallmark collection were utilized to evaluate the cancer-related pathways that are associated with transcription factors and their corresponding target genes. Some important genes, such as NOP58, were positively associated with certain markers (Fig. [104]4F). Fig. 4. [105]Fig. 4 [106]Open in a new tab Correlation analysis between hypoxia and tumor biological processes. A A scatter plot is employed to scrutinize the relationship between hypoxia scores and additional biological functions. B Boxplot comparing the scores of biological processes among different cell subsets. C The infercnv software predicted the CNV status of neoplastic cells, with endothelial cells serving as a comparative baseline. D Network diagram illustrating the interactions of key transcription factors in H2 and N2 subsets. E Venn diagram showing the overlap of key transcription factors in H2 and N2 subsets. F Heatmap analysis of the correlation between H2-specific transcription factors and Hallmark pathways Construction of a prognostic model based on H2-specific transcription factors and genes We employed univariate Cox analysis to sift through genes to identify those linked to patient prognosis. Subsequently, a stepwise regression analysis was conducted to pinpoint genes that have an association with overall survival outcomes. Ultimately, five transcription factors—LRP10, MED8, NOL10, NOP58, and REXO4—were identified for constructing the prognostic model (Fig. [107]5A; Supplementary Fig. [108]5C, D). The LassoCox algorithm was applied to optimize the model construction (Supplementary Fig. [109]5A, B), patients were categorized into two distinct groups—high risk and low risk—based on the median of their risk scores. A correlation analysis was performed between the gene expression of the model and various clinical traits, revealing that the survival rate for patients classified in the high-risk group was markedly inferior to that of the low-risk group, with a statistically significant difference (P < 0.0001) as depicted in Fig. [110]5B. Multivariate Cox regression analysis validated that the model is an independent predictor of patient outcomes for those with HCC, with an odds ratio of 4.63 (95% confidence interval: 2.85 to 7.50), and a P value less than 0.001 (as shown in Fig. [111]5C). The prognostic accuracy of the five transcription factor signatures was confirmed using an external dataset named LCGC-LIRI-JP. This validation demonstrated that individuals with HCC who were categorized in the high-risk group had significantly lower survival rates compared to those in the low-risk group. The survival disparity between the two groups was found to be statistically significant (P < 0.0001), as depicted in Fig. [112]5D. Fig. 5. [113]Fig. 5 [114]Open in a new tab Prognostic model built with H2-specific transcription factors. A A heatmap is constructed to integrate clinical indicators for analyzing the expression patterns of genes within our model. B Survival analysis is conducted to assess the efficacy of the prognostic model using data from the TCGA dataset. C Forest plot displaying the risk scores from the multifactorial Cox analysis. D The predictive accuracy of the prognostic model is confirmed using an external dataset Drug analysis Based on information from the GDSCv2 database, this study calculated correlations between H2-specific transcription factors and genes and drugs, screening transcription factors with|Cor| >0.3 and P < 0.05, and presenting them in a bar chart (Fig. [115]6A). These drugs were involved in multiple biological processes, including DNA replication, apoptosis regulation, genomic integrity, mitosis, and PI3K/MTOR signaling pathways (Fig. [116]6B). NOP58 was found to be significantly associated with multiple drugs targeting different biological pathways, such as PRIMA-1MET targeting p53 signaling, Obatoclax Mesylate regulating apoptosis, AZD7762 acting on cell cycle, Mirin affecting genomic integrity, and Pyridostatin and Epirubicin targeting DNA replication (Fig. [117]6C). To forecast possible pharmaceuticals for supplementary treatment aimed at hypoxia, the online analysis tool of the CMap database (QUERY [clue.io]) was used, with H2-specific transcription factors and genes as Up genes and N2-specific transcription factors and genes as Down genes. A Volcano plots show scores for some potential compounds in different cell lines that may have some inhibitory effect on neoplastic cells (Fig. [118]6D). The top 5 compounds in terms of|Score| were analyzed in each cell line, totaling 44 compounds (Fig. [119]6E). Finally, statistical analysis of the pathways targeted by these compounds was performed (Fig. [120]6F), revealing potential candidates for future hypoxia-targeted therapies. Fig. 6. [121]Fig. 6 [122]Open in a new tab Drug response analysis. A Bar graph analyzing the correlation between H2-specific transcription factors and drug responsiveness. B Bar graph showcasing the biological pathways targeted by drugs. C Network diagram revealing the pathways of drug interactions with the model gene NOP58. D Volcano plot presenting the distribution of compound scores from the CMap database. E Bubble chart showcasing the top 5 compounds scored in each cell line. F Bubble chart analyzing the target pathways of these compounds Discussion Our comprehensive investigation conducted an elaborate examination of the TME of HCC using single-cell RNA sequencing, with a particular concentrate on the effects of hypoxia upon neoplastic and immune cell subpopulations. Our results exposed considerable diversity within the HCC TME under hypoxia conditions. The major hypoxia cell populations comprised neoplastic cells, fibroblasts, endothelial cells, and myeloid cells. Neoplastic cells undergo metabolic reprogramming, such as the Warburg effect, shifting towards anaerobic glycolysis to produce metabolites like lactate, which not only provides energy but also promotes further hypoxia by altering the pH of the microenvironment and increasing lactate concentrations [[123]49]. Additionally, neoplastic cells enhance antioxidant mechanisms, like increased expression of glutathione, Nrf2, and HIF-1α, to combat cellular necrosis caused by increased ROS [[124]50]. Fibroblasts and endothelial cells also exhibit heterogeneity under hypoxia. Fibroblasts release growth factors, cytokines, and chemokines that stimulate the development of cancer-associated fibroblasts. These fibroblasts restructure the neighboring stroma, making it easier for neoplastic cells to migrate and invade [[125]51]. Endothelial cells promote angiogenesis under hypoxia, enhancing tumor oxygenation and nutrient supply, thereby supporting tumor growth and metastasis [[126]52]. Myeloid cells, especially tumor-associated macrophages (TAMs) and MDSCs, have significant functions in the hypoxia TME. TAMs can switch from an M1 to an M2 phenotype, promoting tumor growth and immunosuppression [[127]53]. MDSCs inhibit T-cell activity through secretion of immunosuppressive factors like TGF-β and IL-10, aiding tumor escape from immunological monitoring [[128]54]. Moreover, hypoxia is closely associated with tumor aggressiveness and therapeutic resistance, as it activates HIFs and VEGF signaling, promoting angiogenesis and metastasis [[129]52]. It also fosters tumor stem cell maintenance and phenotypic switching, increasing tumor heterogeneity and therapeutic resistance [[130]55]. These discoveries highlight the intricate function of the oxygen-deprived milieu in the progression of HCC, potentially influencing tumor aggressiveness, therapeutic resistance, and prognosis. The clustering of hypoxia cells in key metastatic regions highlights their potential driving force in tumor progression. The hypoxia microenvironment significantly impacts immune cell subpopulations. Our analysis showed enhanced communication between macrophages, dendritic cells, and monocytes under hypoxia, likely promoting immune evasion and tumor immune editing. Hypoxia, a characteristic feature of the majority of malignant solid neoplasms, fosters immune resistance and suppression/tolerance, adversely affecting antitumor effector cell function [[131]56]. Hypoxia is associated with adverse prognosis, heightened tumor diversity, the rise of chemotherapy-resistant strains, and the bypassing of immunosurveillance [[132]57]. Under low-oxygen conditions, an immunosuppressive state is induced through HIF-1-mediated pathways, which result in the modification of molecular signatures, the movement of immune cells, and the stimulation of new blood vessel formation [[133]58]. Moreover, oxidative stress is known to induce lipid peroxidation, stress in the endoplasmic reticulum, and dysfunction of Tregs, all factors that contribute to the disruption of normal immune responses [[134]58]. In addition, it has been observed that hypoxia conditions can modulate the phenotype of MDSCs in the local environment, with HIF-1α playing a role in this transformation, causing a switch in the suppressive activity of splenic MDSCs from a targeted to a broad-spectrum approach [[135]59]. Macrophages and neutrophils, key innate immune components, display distinct characteristics under hypoxia, with hypoxia influencing morphology, migration, chemotaxis, endothelial cell adhesion, bacterial killing, differentiation/polarization, and protumor activity [[136]60]. Furthermore, hypoxia-induced conditions can lead to an increase in PD-L1 expression on myeloid-derived suppressor cells through the action of HIF-1α, thereby hindering the activation of T-cells [[137]61]. Monocytes’ roles in the TME are also modulated by hypoxia. Bridging innate and adaptive immunity, monocytes can induce Immune passivity, vasculogenesis, and malignant cell propagation. Nevertheless, they are also capable of producing antitumor effector cells and stimulating the activity of antigen-presenting cells [[138]62]. Thus, these results emphasize considering the hypoxia microenvironment in immunotherapy strategies. Therapeutic interventions targeting hypoxia adaptation in specific immune cell subpopulations may be effective. For instance, improving TME oxygenation could weaken protective immunosuppressive signals, unleashing restrained tumor-reactive T cells and NK cells [[139]63]. Additionally, inhibiting HIF-1α may alter MDSC function and differentiation, enhancing antitumor immune responses [[140]59]. The impact of the hypoxia microenvironment on immune cell subpopulations is significant, particularly in promoting immune evasion and tumor immune editing. Therefore, when developing new immunotherapy strategies, it is crucial to consider the influence of the hypoxia microenvironment and implement therapeutic interventions targeting the hypoxia adaptation of specific immune cell subpopulations. These results emphasize the importance of incorporating the hypoxia microenvironment into immunotherapy strategies, potentially necessitating therapeutic interventions tailored to the hypoxia adaptation of specific immune cell subpopulations. Our hypoxia analysis of neoplastic cell subpopulations uncovered an affirmative link between hypoxia and the virulence of the tumor. Notably, the H2 subpopulation, with high hypoxia scores, was associated with high neoplastic cell aggressiveness, potentially linked to hypoxia-induced gene expression changes promoting angiogenesis and extracellular matrix remodeling. Genes like MEG3, KLF6, JUN, upregulated in hypoxia subpopulations, contribute to cell multiplication, movement, infiltration, and vasculogenesis [[141]64–[142]67]. Hypoxia-associated genes enrich in processes like cytoplasmic translation, protein ubiquitination, T-cell activation, and leukocyte cell adhesion, informing future therapeutic strategies targeting the hypoxia microenvironment. Our results showed positive correlations between hypoxia scores and HCC aggressiveness, apoptosis, angiogenesis, and EMT, especially in the H2 subpopulation. This suggests hypoxia may affect HCC prognosis by enhancing neoplastic cell aggressiveness. Hypoxia enhances neoplastic cell invasion and metastasis, correlating with increased malignancy [[143]68, [144]69]. Transcriptional factor analysis identified key genes and transcription factors in H2 and N2 subpopulations, like NOP58 and MED8, highly expressed in H2. These unique transcription factors and genes informed a prognostic model validated in clinical data, offering personalized treatment possibilities. Furthermore, based on information from the GDSCv2 database, we examined the relationship between unique transcription factors and genes in the H2 subpopulation with drugs, particularly the significant association between NOP58 and multiple drugs, providing new targets for future treatment strategies targeting these drugs. NOP58, a core boxC/D snoRNP protein, regulates RNA processing and modification [[145]70], potentially impacting drug metabolism and resistance-related gene expression. NOP58’s role in snoRNP complexes and as a novel candidate marker regulating neoplastic self-renewal and maturation processes [[146]71, [147]72] underscores its therapeutic potential. Similarly, MED8, part of the evolutionarily conserved Mediator complex involved in transcriptional regulation [[148]73], is mutated in colorectal cancer [[149]74] and associated with poor prognosis in renal carcinoma [[150]75]. In HCC, elevated MED19, MED15, MED23 expression correlates with prognosis [[151]76–[152]78]. As HCC often arises in chronic inflammation [[153]79], complement activation, enhancing tumor growth and metastasis [[154]80, [155]81], may drive inflammatory processes. Thus, MED8 may influence HCC prognosis by modulating complement activation pathways. This study primarily relies on comprehensive computational analyses and bioinformatic modeling to explore hypoxia-related heterogeneity in HCC. Although we identified key hypoxia-associated genes such as LRP10, NOP58, and NOL10 with potential roles in tumor progression, their biological functions under hypoxia conditions remain to be experimentally validated. Functional assays including gene knockdown or overexpression in HCC cell lines exposed to hypoxia would provide critical evidence, but due to current limitations in our laboratory resources and experimental conditions, these investigations were not feasible within the scope of this work. In addition, our hypoxia scoring model (CHPF) has not yet been systematically benchmarked against established hypoxia signatures such as Buffa, Ragnum, and Winter scores, nor validated across independent single-cell or bulk RNA-seq datasets. Similarly, experimental validation of drug responses predicted by our model, for example through treatment of HCC cell lines or patient-derived organoids under hypoxia conditions, remains to be performed. These represent important next steps to confirm the clinical relevance and therapeutic potential of our findings. Moreover, the single-cell RNA sequencing dataset used comprised a limited number of samples. While the original dataset included about ten tumor samples, only three were selected for our analysis based on stringent criteria including data quality, clinical annotation completeness, and sequencing consistency. This small sample size may introduce sampling bias and limit the generalizability of our conclusions. Future studies incorporating larger, more diverse patient cohorts are necessary to validate the robustness of the hypoxia-associated cellular subpopulations and prognostic gene signatures identified here, and to better capture the inter-patient heterogeneity characteristic of HCC. Despite these limitations, our study offers valuable insights into the multifaceted role of hypoxia in shaping the tumor microenvironment and highlights promising molecular targets. We hope these findings will motivate further experimental and clinical research to deepen understanding and improve therapeutic strategies for HCC. Conclusion In summary, this study used CHPF to identify cellular hypoxia states, elucidating intratumoral heterogeneity in HCC. We unveiled hypoxia’s multifaceted roles in tumor development and constructed a prognostic model. By isolating a unique hypoxia subpopulation with high invasive potential, we discovered distinct transcriptomic and genomic alterations. Key transcription factors like LRP10, MED8, NOL10, NOP58, and REXO4 play pivotal roles in hypoxia-induced HCC invasion and metastasis, and a predictive model based on these factors forecasts HCC survival. Our discoveries offer new molecular indicators and therapeutic objectives for the detection and management of HCC, guiding upcoming studies and medical procedures. Electronic supplementary material Below is the link to the electronic supplementary material. [156]12672_2025_3201_MOESM1_ESM.tif^ (2MB, tif) Supplementary Material 1: Figure 1 Cell communication analysis. A Heatmap displaying the patterns of cell-cell communication in single-cell data. B Heatmap analyzing the communication of the MIF signaling pathway among different cell types [157]12672_2025_3201_MOESM2_ESM.tif^ (5.5MB, tif) Supplementary Material 2: Figure 2 Features of neoplasm cellular subpopulations. A UMAP plot depicting the distribution of tumor cell subsets. B UMAP plot presenting the hypoxia scores of neoplastic cells. C Boxplot analyzing the hypoxia scores of hypoxia subsets [158]12672_2025_3201_MOESM3_ESM.tif^ (2.4MB, tif) Supplementary Material 3: Figure 3 Hallmark pathway analysis. Heatmap displaying the GSVA scores of Hallmark pathways in various tumor subsets [159]12672_2025_3201_MOESM4_ESM.tif^ (903.4KB, tif) Supplementary Material 4: Figure 4 CNV analysis. Boxplot presenting the CNV scores of tumor subsets [160]12672_2025_3201_MOESM5_ESM.tif^ (5.2MB, tif) Supplementary Material 5: Figure 5 Prognostic model analysis. A, B Presentation of modeling results using the LassoCox algorithm. C Forest plot analyzing the univariate Cox risks of model genes. D Bar graph displaying the coefficients of genes in the prognostic model Acknowledgements