Lactate metabolism-driven tumor heterogeneity and molecular signatures in intrahepatic cholangiocarcinoma

Abstract

BACKGROUND

Intrahepatic cholangiocarcinoma (ICC) is a highly aggressive liver malignancy with limited therapeutic options and poor prognosis. Recent evidence indicates that lactate metabolism (LM) plays a pivotal role in tumor metabolic reprogramming, immune evasion, and disease progression; however, the heterogeneity and regulatory mechanisms of LM activity within ICC remain largely undefined.

AIM

To systematically characterize LM-driven heterogeneity and its molecular and functional implications in ICC.

METHODS

Single-cell RNA sequencing and bulk transcriptomic datasets were integrated to characterize LM heterogeneity in ICC. High-dimensional weighted gene co-expression network analysis and multiple machine-learning algorithms (least absolute shrinkage and selection operator, random forest, gradient boosting machine, adaptive best subset selection, and decision tree) were employed to identify LM-associated feature genes. CytoTRACE and CellChat analyses were used to assess differentiation potential and intercellular communication among malignant epithelial subpopulations. Kyoto Encyclopedia of Genes and Genomes and Gene Ontology enrichment analyses were performed to elucidate biological functions. A random forest model combined with SHapley Additive exPlanation (SHAP) interpretability analysis identified the most predictive LM-related gene. Functional assays, including quantitative polymerase chain reaction, cell counting kit-8, colony formation, wound-healing, and transwell experiments, were conducted to validate CYC1 in ICC cell lines.

RESULTS

Malignant ICC cells were stratified into three LM-activity subtypes (high, intermediate, and low) exhibiting distinct transcriptional programs and differentiation trajectories. Twelve LM-associated feature genes GPX3, CYC1, NME1, GSTP1, MGST1, ALDH3A1, TALDO1, SNRPB, TKT, NAA20, G6PD, and RPL13A were identified as key molecular markers linked to aggressive phenotypes and poor prognosis. Among them, CYC1 showed the highest predictive accuracy (area under the curve = 0.844) and strongest model contribution (SHAP = 0.091), marking it as the principal LM-related driver gene. Functional experiments confirmed that CYC1 knockdown significantly suppressed ICC cell proliferation, migration, and invasion, validating its oncogenic role in promoting malignant progression.

CONCLUSION

This integrative single-cell and machine-learning study delineates the molecular heterogeneity of LM in ICC and identifies twelve feature genes linking LM with tumor aggressiveness. These findings provide novel insight into LM-driven oncogenic mechanisms and propose CYC1 and other LM-associated genes as potential biomarkers and therapeutic targets for ICC.

Key Words: Intrahepatic cholangiocarcinoma; Lactate metabolism; Single-cell RNA sequencing; Metabolic heterogeneity; Tumor microenvironment

Core Tip: This study integrates single-cell and bulk transcriptomic data with machine learning to uncover lactate metabolism-driven heterogeneity in intrahepatic cholangiocarcinoma (ICC). Twelve lactate metabolism-related genes were identified as prognostic biomarkers, and functional validation of CYC1 revealed its role in promoting tumor invasion. These findings provide new insights into ICC progression and offer promising targets for diagnosis and therapy.

INTRODUCTION

Intrahepatic cholangiocarcinoma (ICC) is a highly aggressive liver malignancy originating from the epithelial cells of the intrahepatic bile ducts. Its incidence has been steadily increasing in recent years, with China accounting for over half of all cases worldwide[1-3]. Due to the lack of early clinical symptoms and effective diagnostic biomarkers, most ICC patients are diagnosed at advanced stages. Despite recent advances in targeted therapy and immunotherapy, prognosis remains dismal, largely due to the high degree of tumor heterogeneity and the complex immune microenvironment[4,5]. Therefore, identifying reliable molecular biomarkers is essential for improving diagnostic precision, guiding therapeutic strategies, and understanding the underlying biology of ICC progression.

Metabolic reprogramming is a hallmark of cancer, and lactate metabolism (LM) has emerged as a key contributor to tumor growth and immune evasion. The Warburg effect, describing cancer cells’ preference for glycolysis over oxidative phosphorylation even under normoxia, leads to excessive lactate accumulation. This not only provides an alternative energy source but also acidifies the tumor microenvironment, thereby suppressing antitumor immunity and promoting invasion and metastasis[6,7]. Beyond its role in energy metabolism, lactate has emerged as a critical signaling molecule that influences various biological processes within the tumor microenvironment, regulating gene expression, protein modification, and cell differentiation, ultimately enhancing tumor aggressiveness[8-10].

Recent studies have further linked LM to therapy resistance and immune modulation. Li et al[11] showed that tumor-derived lactate upregulates RUBCNL to promote bevacizumab resistance in colorectal cancer. Qian et al[12] demonstrated that monocarboxylate transporter 4-mediated lactate secretion in STK11/LKB1-mutant lung adenocarcinoma induces an immunosuppressive environment via M2 macrophage polarization and T-cell inhibition. Lactate-induced post-translational modifications, such as lactylation, have also emerged as key regulators of tumor biology. In the study of ICC, the ribosomal protein nucleolin has been identified as a lactylation target, with modification at lysine 477 mediated by the acetyltransferase p300 promoting tumor proliferation and invasion[13]. Similarly, p-PCK1 has been shown to interact with LDHA to enhance lactate production, leading to KAT7-mediated lactylation of SPRING at lysine 82, thereby increasing resistance to ferroptosis in ICC cells and suggesting p-PCK1 as a potential biomarker[14]. These findings underscore the multifaceted role of lactate in shaping tumor biology and therapeutic response.

However, the heterogeneity of LM activity within ICC and its regulatory mechanisms remain largely unexplored. A comprehensive understanding of LM-driven molecular alterations at single-cell resolution could reveal new therapeutic vulnerabilities. In this study, we integrated single-cell RNA sequencing (scRNA-seq) and bulk transcriptomic data to systematically dissect LM heterogeneity in ICC. By combining high-dimensional weighted gene co-expression network analysis (hdWGCNA) with multiple machine learning algorithms [least absolute shrinkage and selection operator (LASSO), random forest, gradient boosting machine (GBM), adaptive best subset selection (ABESS), and decision tree], we identified twelve key LM-associated genes and constructed predictive models. Among them, CYC1, a mitochondrial electron transport gene, exhibited the highest predictive power and strongest association with LM activity. Functional validation confirmed that CYC1 knockdown significantly inhibited proliferation, migration, and invasion in ICC cells, supporting its role as a potential prognostic and therapeutic target.

Together, our findings provide a comprehensive framework for understanding LM-associated heterogeneity and molecular signatures in ICC, linking metabolic reprogramming to malignant progression.

MATERIALS AND METHODS

Data acquisition

ScRNA-seq datasets (GSE138709[15], GSE171899[16], GSE151530[17,18], GSE189903[19], and GSE201425[20]) comprising 42 ICC and 12 normal samples were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Bulk RNA-seq data (36 ICC and 9 normal samples) were retrieved from The Cancer Genome Atlas-cholangiocarcinoma (TCGA-CHOL) database (https://portal.gdc.cancer.gov/projects/TCGA-CHOL) (Supplementary Table 1). Based on a previous study[21], LM-related genes were collected from the Molecular Signatures Database (MSigDB) database (https://www.gsea-msigdb.org/gsea/msigdb) (Supplementary Table 2).

Data processing and scRNA-seq analysis

ScRNA-seq data were processed using the “Seurat” R package (4.3.0). High-quality cells expressing between 200-7000 genes with < 20% mitochondrial content were retained. Data were normalized with the “Log-normalization” method, scaled with “scale data”, and the top 3000 highly variable genes were identified using “find variable features”. Principal component analysis was performed for dimensionality reduction, and batch effects were corrected using the “Harmony” (0.1.0) algorithm[22]. Cell clustering was conducted via “FindClusters” with a resolution parameter of 0.8, and the uniform manifold approximation and projection (UMAP) algorithm was used for visualization. Cluster annotation was guided by canonical cell-type markers and distinct gene expression profiles to ensure accurate cell-type identification. To identify malignant cells exhibiting extensive chromosomal copy number variations (CNVs), we inferred CNV profiles using the “inferCNV” package (1.8.1). A CNV score was computed as the mean of squared CNV values across all chromosomes, representing the overall amplitude of chromosomal alterations. Cells were classified as malignant or non-malignant based on the bimodal distribution of CNV-derived malignancy scores relative to reference cells[23,24].

Evaluation of LM activity

To evaluate LM activity at the single-cell level, five computational algorithms were employed: AUCell (1.12.2)[25], UCell (1.3.1)[26], singscore (1.12.0)[27], single-sample Gene Set Enrichment Analysis (ssGSEA) (1.40.1)[28], and AddModuleScore (4.3.0)[29]. The resulting LM scores from each algorithm were normalized to a 0-1 scale and integrated to generate a composite LM activity score. Each method provided a distinct mechanism for evaluating gene set activity. Gene expression rankings were used by AUCell to determine LM activity, while UCell normalized rank scores to estimate enrichment. ssGSEA generated relative enrichment scores by comparing the expression levels of genes within and outside the LM set, and AddModuleScore computed the weighted average expression of gene sets, normalizing the results according to the singscore ranking framework. After transforming the raw score matrix to a zero-mean, unit-variance scale and rescaling values to the (0, 1) range for consistency, Z-score standardization and min-max normalization was applied to enhance comparability across algorithms. The normalized feature values were then integrated to produce a composite LM activity score. To classify cells into biologically meaningful categories, a tertile-based cutoff at the 33^rd and 66^th percentiles was applied, defining low-LM (LLM), intermediate-LM (ILM), and high-LM (HLM) groups. These thresholds corresponded to inflection points in the LM score distribution, indicating natural metabolic transitions. This integrative scoring and classification framework provided a robust and quantitative assessment of LM activity across single cells.

hdWGCNA

hdWGCNA was performed to identify gene modules associated with LM and transcriptional heterogeneity in ICC. Genes were grouped into co-expression modules based on pairwise correlation of expression profiles, resulting in the construction of a weighted gene co-expression network[30,31]. Module-trait association analyses were conducted to explore the relationship between LM activity and module clustering, with a particular focus on modules implicated in tumor progression. Modules exhibiting strong correlations with elevated LM activity were considered functionally significant. Intra-module hub genes were identified as key regulatory elements potentially driving malignant transformation and intratumoral heterogeneity in ICC.

Intercellular communication analysis

Intercellular communication within the tumor microenvironment was inferred from scRNA-seq data using the “CellChat” R package (1.6.1)[32]. Annotated cell populations were analyzed to construct comprehensive communication networks. Overall signaling patterns between cell types were visualized using netVisual circle, while netVisual bubble was applied to depict significant ligand-receptor interactions and their associated signaling pathways.

Functional enrichment analysis

Functional pathway activity across cell groups was assessed using the “GSVA” R package (1.40.1) with gene sets from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb), specifically the c2.cp.kegg.v7.4.symbols.gmt [Kyoto Encyclopedia of Genes and Genomes (KEGG) canonical pathways] and h.all.v2022.1.Hs.symbols.gmt (Feature gene sets) collections[33]. To further explore biological functions, KEGG and Gene Ontology (GO) enrichment analyses (biological processes, cellular components, and molecular functions), were conducted using the “clusterProfiler” R package (4.0.5), facilitating interpretation of molecular functions and pathway interactions.

Cellular trajectory inference

The developmental potential of malignant cells was assessed using the “CytoTRACE” (0.3.3) which assigns scores based on transcriptional diversity, with higher scores indicating lower differentiation and greater proliferative capacity[34]. To further investigate cellular differentiation trajectories, the “Monocle” algorithm (2.20.0) was applied by projecting high-dimensional scRNA-seq data into a low-dimensional space, enabling pseudotemporal ordering and identification of lineage progression among malignant cell populations[35].

Feature gene selection and LM activity prediction based on machine learning

LM-associated feature genes were identified using a consensus-based feature selection strategy, integrating five distinct algorithms: Random forest, GBM, decision tree, LASSO and ABESS with random forest (4.7.1), glmnet (4.1.7), gbm (2.1.8), rpart (4.1.16), and abbess (0.4.6) R packages[36-39]. Genes repeatedly selected by two or more algorithms were defined as LM core genes, ensuring stability and minimizing method-specific bias. To predict LM activity at the single-cell level, eight machine learning algorithms including k-nearest neighbor, naive Bayes, random forest, recursive partitioning and regression trees, support vector machines, linear discriminant analysis, and extreme gradient boosting, were implemented using the ‘mlr3’ R package (0.15.0). The dataset was split into training (80%) and testing (20%) subsets, models were tuned via five-fold internal cross-validation, and generalization was assessed using ten-fold external cross-validation. Model performance was evaluated by mean area under the curve (AUC), with the best-performing model selected. SHapley Additive exPlanation (SHAP) values were used to quantify each gene’s contribution to LM activity prediction[39].

Cell culture

Human ICC cell lines HCCC9810 and HuCCT1 were obtained from the Cell Resource Center of the Shanghai Institute of Life Sciences. HuCCT1 cells were cultured in RPMI-1640 medium (Gibco, United States) supplemented with 10% fetal bovine serum (FBS) (Gibco, United States), while HCCC9810 cells were maintained in Dulbecco’s Modified Eagle’s Medium (Gibco, United States) containing 10% FBS. All cells were incubated at 37 °C in a humidified atmosphere with 5% carbon dioxide.

Quantitative real-time polymerase chain reaction

Total RNA was extracted using TRIzol reagent (Invitrogen, United States). Reverse transcription was performed with 500 ng RNA using the HiScript RT Mix Kit (Vazyme, China), and quantitative polymerase chain reaction (PCR) was carried out using the SYBR Green Kit (Vazyme) on a real-time PCR system (Applied Biosystems). Glyceraldehyde-3-phosphate dehydrogenase was used as internal control. Primer sequences were synthesized by Qingdao Biotechnology Company (China) and are listed in Supplementary Table 3.

Cell proliferation

Cell proliferation was assessed using the cell counting kit-8 (CCK-8) (Vazyme, China). Cells (2 × 10³ per well) were seeded into 96-well plates, and 10 μL of CCK-8 reagent was added at the indicated time points (days 1-5). After incubation at 37 °C for 2 hours, absorbance at 450 nm was measured using a microplate reader (Thermo, United States).

Colony formation assays

For colony formation, 1 × 10³ cells were seeded into 6-well plates and cultured for approximately 14 days until visible colonies appeared. Colonies were fixed with 4% paraformaldehyde for 15 minutes, stained with 0.1% crystal violet (Solarbio, China) for 20 minutes, air-dried, and counted under a microscope.

Transwell assay

Cell migration and invasion were evaluated using Transwell assay. A total of 2 × 10⁴ cells suspended in 200 μL serum-free medium were added to the upper chamber, which was either uncoated (migration) or Matrigel-coated (invasion) (BD Biosciences, United States). The lower chamber contained 600 μL medium with 10% FBS as a chemoattractant. After incubation, cells that migrated or invaded were fixed with 4% paraformaldehyde, stained with 0.1% crystal violet, and counted under a light microscope.

Statistical analysis

All statistical analyses were conducted in R (version 4.1.3). Group comparisons were performed using the Student’s t-test or Mann-Whitney U test for continuous variables, and the χ² or Fisher’s exact test for categorical variables, as appropriate. For transcriptomic analyses, two-sided tests were applied, and false discovery rate correction was used to adjust for multiple comparisons. Pearson’s correlation assessed variable associations. Statistical significance was defined as two-tailed P < 0.05.

RESULTS

Single-cell transcriptomic landscape of ICC

To comprehensively delineate the cellular heterogeneity of ICC, we integrated five publicly available scRNA-seq datasets comprising a total of 42 ICC tumor samples and 12 matched normal controls. We assessed key quality control metrics including the number of detected genes per cell (nFeatureRNA), total transcript counts (nCountRNA), and the proportions of mitochondrial (percent.MT), hemoglobin (percent.HB), and ribosomal (percent.Ribosome) gene expression (Supplementary Figure 1A), to ensure data reliability and exclude low-quality cells. Correlation analysis revealed a strong association between nFeatureRNA and nCountRNA (r = 0.93), indicating that cells with higher transcript counts also tend to express more genes, while percent.MT and percent.Ribosome showed weak or negative correlations (Supplementary Figure 1B). After batch correction, UMAP embedding demonstrated consistent integration across all datasets and samples (Supplementary Figure 1C).

To characterize the cellular composition of ICC tissue, we applied UMAP for dimensionality reduction, identifying 31 transcriptionally distinct cell clusters (Figure 1A). Based on canonical marker expression, 11 major cell types were annotated: B cells, T/natural killer (NK) cells, plasma cells, monocytes, macrophages, dendritic cells (DCs), fibroblasts, endothelial cells, epithelial cells, mast cells, and hepatocytes (Figure 1B). Representative marker genes used for cell-type identification were as follows: B cells (MS4A1, CD79A, CD79B, CD19); T/NK cells (CD3D, CD3E, IL7R, CCL5, NKG7); Monocytes (S100A12, S100A9, FCN1, LYZ); Macrophages (CD68, CD163, C1QA, C1QB, MARCO); DCs (CD1C, CD1E, CLEC10A, CLEC9A); Fibroblasts (PDGFRA, COL1A1, DCN, RGS5, LUM); Mast cells (KIT, GATA2, CPA3); Endothelial cells (VWF, CD34, ICAM2); Epithelial cells (EPCAM, CDH1, KRT7, KRT18); Plasma cells (CD38, MZB1, XBP1, PRDM1); Hepatocytes (ALB, TTR, APOC3). Dot-plot visualization confirmed the specific expression patterns of these marker genes across distinct clusters (Figure 1C). UMAP feature plots further illustrated the expression density of representative marker genes across the identified cell subgroups, supporting accurate cell-type annotation and highlighting intratumoral heterogeneity (Figure 1D).

Figure 1 Single-cell landscape and cell-type annotation of intrahepatic cholangiocarcinoma. A: Uniform manifold approximation and projection (UMAP) projection identifying 31 transcriptionally distinct clusters; B: Cell-type annotation based on canonical marker gene expression, revealing 11 major populations; C: Dot-plot visualization of representative marker genes used for each annotated cell type and their relative expression level; D: UMAP feature plots showing spatial expression patterns of selected canonical markers. UMAP: Uniform manifold approximation and projection; NK: Natural killer; DC: Dendritic cell.

Evaluation of LM activity in ICC

To investigate the functional relevance of LM in ICC, we performed a multi-level analysis integrating bulk and single-cell transcriptomic data. In the TCGA-CHOL cohort, patients were stratified by LM-related gene expression, revealing that those with high LM activity exhibited significantly reduced overall survival (OS) compared to the low LM activity group (P = 0.0043, Figure 2A). At single-cell resolution, LM activity was quantified and compared across cellular subpopulations within the ICC microenvironment, revealing that tumor-derived bile duct epithelial cells displayed the highest LM activity (Figure 2B). To distinguish malignant from non-malignant epithelial cells, infer CNV analysis was applied, identifying six epithelial clusters based on chromosomal CNV profiles (Figure 2C). Clusters K1, K2, K4, K5, and K6 exhibited significantly higher CNV alterations compared to, indicating their malignant nature (P < 0.001, Figure 2D). Consistently, malignant epithelial cells demonstrated significantly elevated LM activity relative to normal epithelial cells (P < 0.001, Figure 2E). Five computational algorithms AUCell, UCell, singscore, ssGSEA, and AddModuleScore were employed to quantify LM activity across cell types. All methods consistently indicated heightened LM activity in malignant epithelial cells within the ICC microenvironment (Figure 2F and G). UMAP visualization illustrated the spatial distribution of LM activity throughout the ICC microenvironment, highlighting strong enrichment within malignant epithelial clusters (Figure 2H). Correspondingly, heatmap analysis confirmed that LM activity varied across cell types, with malignant epithelial cells displaying the highest LM activity levels (Figure 2I). These findings underscore the potential functional significance of LM activity in ICC pathogenesis.

Figure 2 Evaluation of lactate metabolism activity in intrahepatic cholangiocarcinoma. A: Kaplan-Meier analysis of the Cancer Genome Atlas-cholangiocarcinoma cohort showing that patients with high lactate metabolism (LM) activity exhibited significantly reduced overall survival compared to the low LM group (P = 0.0043); B: Violin plots comparing LM activity across major cell types within the intrahepatic cholangiocarcinoma (ICC) microenvironment; C: Heatmap of infer copy number variation (CNV) analysis identifying six epithelial clusters based on chromosomal CNV profiles; D: Boxplots illustrating significantly higher CNV scores in malignant clusters (K1, K2, K4, K5, and K6) compared with non-malignant cluster K3 (P < 0.001); E: Comparison of LM activity between normal and malignant epithelial cells, showing markedly elevated LM activity in the malignant population (P < 0.001); F and G: LM activity quantified by five computational algorithms (AUCell, UCell, singscore, single-sample Gene Set Enrichment Analysis, and AddModuleScore), consistently indicating higher LM activity in malignant epithelial cells; H: Uniform manifold approximation and projection visualization displaying the spatial distribution of LM activity scores across the ICC microenvironment; I: Heatmap of LM score distribution highlighting cell-type-specific variation. ^cP < 0.001. ^dP < 0.0001. NS: Not significant; TCGA-CHOL: The Cancer Genome Atlas-cholangiocarcinoma; LM: Lactate metabolism; NK: Natural killer; DC: Dendritic cell; CNV: Copy number variation; ssGSEA: Single-sample Gene Set Enrichment Analysis; UMAP: Uniform manifold approximation and projection.

Heterogeneity of LM activity in ICC malignant cells

To explore the heterogeneity of LM activity within malignant biliary epithelial cells, we analyzed the malignant cell population identified by infer CNV. LM activity scores displayed an approximately normal distribution (Figure 3A). Based on score thresholds, 4270 cells (66%) with LM scores > 2.83 were classified as the HLM activity group, 4145 cells (33%) with LM scores < 2.38 were designated as the LLM activity group, and 4144 cells with intermediate scores were categorized as the ILM activity group. These three subgroups exhibited distinct spatial expression characteristics and density distribution within the ICC microenvironment (Figure 3B and C). To evaluate the differentiation potential of these malignant subpopulations, we applied CytoTRACE analysis to single-cell transcriptomic data. CytoTRACE scores were significantly higher in HLM cells than in ILM or LLM cells (Figure 3D and E), consistent with a positive correlation between LM activity and CytoTRACE scores (r = 0.60; Figure 3F and G). This indicates that cells with elevated LM activity are transcriptionally less differentiated and possess enhanced developmental plasticity. Integrating CytoTRACE with pseudo time inference further delineates the differentiation trajectory of malignant epithelial cells. The trajectory originated from the HLM subpopulation and progressed toward ILM and LLM states, outlining a continuous differentiation landscape within the malignant compartment (Figure 3H). Collectively, these findings demonstrate marked intratumoral heterogeneity of LM activity, which is closely coupled with the differentiation hierarchy of malignant epithelial cells in ICC.

Figure 3 Lactate metabolism heterogeneity intrahepatic cholangiocarcinoma malignant cells. A: Distribution of lactate metabolism (LM) activity scores in malignant epithelial cells. Cells were stratified into low-LM (LLM) (< 2.38), intermediate-LM (ILM), and high-LM (HLM) (> 2.83) activity groups; B: Spatial distribution of LLM, ILM, and HLM subgroups; C: Density maps illustrating distinct spatial patterns of LM activity; D: Differentiation potential across malignant cells; E: CytoTRACE scores among subgroups; F: The co-expression relationship of LM scoring and CytoTRACE score; G: Correlation plot of LM activity and CytoTRACE scores (r = 0.60, Pearson), showing a strong positive association; H: The development trajectory of malignant epithelial cell subpopulation. ^cP < 0.001. LM: Lactate metabolism; LLM: Low-lactate metabolism; ILM: Intermediate-lactate metabolism; HLM: High-lactate metabolism; UMAP: Uniform manifold approximation and projection; CI: Confidence interval.

Intercellular communication associated with LM activity in ICC

To elucidate how LM activity influences tumor microenvironment communication, we compared intercellular signaling between HLM and LLM malignant epithelial subgroups. HLM cells exhibited a markedly higher number and overall strength of intercellular interactions than LLM cells (Figure 4A). Analysis of signaling patterns revealed that endothelin signaling, calcitonin receptor signaling, growth differentiation factor family, epidermal growth factor family, annexin family signaling (ANNEXIN), and secreted phosphoprotein 1 dominated outgoing signals, whereas insulin-like growth factor family, ANNEXIN, pleiotrophin, midkine (MDK), and macrophage migration inhibitory factor (MIF) were major contributors to incoming signaling (Figure 4B). Quantitative evaluation showed that HLM cells displayed stronger outgoing interactions, while LLM cells exhibited relatively enhanced incoming signaling activity (Figure 4C). Ligand–receptor network analysis identified MIF and MDK as pivotal mediators linking LM activity with microenvironmental communication, suggesting their potential role in ICC progression (Figure 4D). Consistent with this, HLM cells showed a substantially greater number of both ligand and receptor interactions across diverse cell types compared with LLM cells (Figure 4E and F). Pathway activation profiling using PROGENy revealed significant enrichment of oncogenic pathways in HLM cells, including epidermal growth factor receptor, mitogen-activated protein kinase pathway, Janus kinase-signal transducer and activator of transcription, p53, and androgen signaling (Figure 4G). This pattern indicates a strong association between elevated LM activity and tumor proliferation, metastasis, and immune evasion. Complementary GSVA analysis demonstrated that HLM cells were enriched in oxidative phosphorylation, MYC TARGETS V1, and adipogenesis pathways, whereas LLM cells were enriched in allograft rejection, KRAS signaling DN, and spermatogenesis (Figure 4H). Additionally, HLM cells were extensively involved in various metabolic pathways, including drug metabolism, fatty acid metabolism, and glycolysis, highlighting their metabolic reprogramming (Figure 4I).

Figure 4 Biological function analysis of lactate metabolism activity subgroups in intrahepatic cholangiocarcinoma. A: Number and strength of intercellular interactions in intrahepatic cholangiocarcinoma (ICC) tumor microenvironment; B: Outgoing and incoming signaling patterns in ICC microenvironment; C: Scatter plot comparing outgoing vs incoming interaction strength of cell types; D: Ligand-receptor interaction analysis for malignant lactate metabolism (LM) cells; E and F: Interaction pairs of cell types interacting with malignant LM activity cells in providing ligand and receptor; G: PROGENy pathway activity heatmap showing elevated oncogenic signaling in LM activity subgroups; H: GSVA enrichment of hallmark pathways distinguishing high-LM and low-LM; I: Metabolic pathway enrichment analysis in LM activity subgroups. NK: Natural killer; DC: Dendritic cell; LLM: Low-lactate metabolism; ILM: Intermediate-lactate metabolism; HLM: High-lactate metabolism; PI3K: Phosphatidylinositol 3-kinase; TGF: Transforming growth factor; VEGF: Vascular endothelial growth factor; NF-κB: Nuclear factor kappa-B; TNF: Tumor necrosis factor; JAK-STAT: Janus kinase-signal transducer and activator of transcription; MAPK: Mitogen-activated protein kinase; EGFR: Epidermal growth factor receptor.

Identification of signature genes driving LM activity

To delineate key transcriptional regulators underlying LM activity heterogeneity in ICC, we performed hdWGCNA on malignant epithelial cells. A scale-free topology fit index of 0.8 was achieved using a soft-thresholding power of β = 14 (Figure 5A), yielding eight distinct co-expression modules (Figure 5B). The top 10 hub genes from each module and their module eigengene correlations were visualized (Figure 5C). UMAP projections revealed spatially distinct expressions of module-specific genes (Figure 5D). Comparative module analysis indicated that the red, blue, brown, black, and green modules were predominantly enriched in the HLM group, whereas the yellow and turquoise modules were more associated with the LLM group (Figure 5E). Collectively, 500 hub genes identified across these modules (Supplementary Table 4) represent potential regulators of LM-associated transcriptional programs in ICC.

Figure 5 Identification of lactate metabolism-associated gene co-expression modules in malignant intrahepatic cholangiocarcinoma cells. A: Determination of soft-thresholding power (β = 14) achieving scale-free topology fit ≥ 0.8; B: Gene dendrogram and module color assignment by high-dimensional weighted gene co-expression network analysis; C: Module eigengene correlations and top hub genes per module; D: Uniform manifold approximation and projection visualization of module-specific gene expression across malignant cells; E: Module-trait relationship plot showing relative enrichment of modules across lactate metabolism activity subgroups. hdWGCNA: High-dimensional weighted gene co-expression network analysis; kME: Eigengene connectivity; LLM: Low-lactate metabolism; ILM: Intermediate-lactate metabolism; HLM: High-lactate metabolism.

To further characterize genes linked to elevated LM activity, differential expression analysis identified 627 upregulated genes in HLM malignant cells with average log2 fold change (average log₂ fold change) > 0.25 and adjusted P value < 0.05. (Supplementary Figure 2A and Supplementary Table 5). Intersection of these differentially expressed genes with hdWGCNA hub genes yielded 232 shared genes (Supplementary Figure 2B and Supplementary Table 6). Pearson correlation analysis subsequently identified 190 signature genes significantly associated with LM activity (Supplementary Figure 2C and Supplementary Table 7). Functional enrichment analyses revealed that these LM-associated genes were significantly significant associations with autoimmune diseases, metabolic processes, and tumor-related signaling pathways (Supplementary Figure 2D and Supplementary Table 8). GO enrichment analysis further indicated significant enrichment in biological processes such as ribonucleoprotein complex biogenesis (GO: 0022613), translational initiation (GO: 0006413), and viral gene expression (GO: 0019080). Cellular component analysis showed predominant localization to the ribosome (GO: 0005840) and mitochondrial inner membrane (GO: 0005743), while molecular function terms were enriched for structural constituents of the ribosome (GO: 0003735), cadherin binding (GO: 0045296), and oxidoreductase activity acting on nicotinamide adenine dinucleotide phosphate (NADPH) (GO: 0016655) (Supplementary Figure 2E and Supplementary Table 9). Together, these findings define a set of LM-associated signature genes driving metabolic and translational reprogramming in malignant ICC cells.

Selection and prognostic evaluation of genes characterizing LM activity

To identify robust marker genes associated with elevated LM activity in ICC malignant epithelial cells, we applied a multi-algorithm machine learning framework integrating five models: Random forest, GBM, LASSO, decision tree, and ABESS. Random forest and GBM analyses each identified the top 20 genes ranked by feature importance (Figure 6A and B; Supplementary Tables 10 and 11). LASSO regression, optimized through 10-fold cross-validation, selected 104 predictive genes (Figure 6C; Supplementary Table 12). Decision tree and ABESS models further contributed 20 and 60 genes, respectively (Figure 6D and E; Supplementary Tables 13 and 14). Integrating all model outputs revealed 12 consensus genes strongly associated with high LM activity: GPX3, CYC1, NME1, GSTP1, MGST1, ALDH3A1, TALDO1, SNRPB, TKT, NAA20, G6PD, and RPL13A (Figure 6F). These genes represent stable LM-related molecular signatures for downstream functional and prognostic analyses.

Figure 6 Multi-algorithm identification of lactate metabolism-related feature genes. A: Importance of signature genes evaluated in random forest; B: Importance of signature genes generated from gradient boosting machine; C: Optimal lambda tuning and cross-validation in least absolute shrinkage and selection operator; D: Importance of signature genes generated from decision tree; E: Candidate optimal signature genes obtained in adaptive best subset selection; F: Intersection diagram summarizing overlapping genes across the five algorithms, identifying 12 shared lactate metabolism-associated genes. GBM: Gradient boosting machine; DT: Decision tree; ABESS: Adaptive best subset selection; LASSO: Least absolute shrinkage and selection operator.

To assess the prognostic relevance of these feature genes, we performed receiver operating characteristic (ROC) analysis at the single-cell level. The AUC values demonstrated strong predictive performance for each gene: CYC1 (AUC = 0.844), NME1 (AUC = 0.831), TKT (AUC = 0.819), MGST1 (AUC = 0.816), G6PD (AUC = 0.806), ALDH3A1 (AUC = 0.796), NAA20 (AUC = 0.795), GPX3 (AUC = 0.794), SNRPB (AUC = 0.787), TALDO1 (AUC = 0.783), RPL13A (AUC = 0.777) and GSTP1 (AUC = 0.758) (Supplementary Figure 3A). Expression profiling revealed that these genes were predominantly expressed in HLM malignant cells, with minimal expression in LLM cells (Supplementary Figure 3B). UMAP-based visualization confirmed distinct spatial enrichment within the malignant compartment (Supplementary Figure 3C). Bulk RNA-seq data further validated the upregulation of most feature genes in ICC tumors compared with normal tissues, except ALDH3A1, GPX3, MGST1, and RPL13A, which displayed intermediate differences (Supplementary Figure 3D). Correlation analysis revealed strong positive associations between LM activity scores and the expression levels of individual feature genes, highlighting their potential regulatory roles (Supplementary Figure 3E). Kaplan-Meier survival analysis in the TCGA-CHOL cohort revealed that elevated expression of GSTP1, MGST1, NAA20, NME1, RPL13A, SNRPB, TKT, and TALDO1 was significantly associated with reduced OS (Supplementary Figure 3F). Collectively, these results highlight a 12-gene LM-related signature that not only reflects metabolic reprogramming within ICC malignant cells but also predicts patient prognosis.

Construction and interpretation of the optimal predictive model

Although 12 feature genes were identified through five machine-learning algorithms, their comparative predictive performance required further evaluation. To systematically benchmark model efficiency and stability, we constructed multiple classifiers and compared their performance. Among all algorithms tested, the random forest model achieved the highest mean AUC and exhibited the most consistent performance across resampling iterations (Figure 7A; Supplementary Table 15). ROC and precision-recall (PR) analyses confirmed the superior accuracy of the random forest classifier, with both AUC and PRAUC values reaching 0.97 in the test set, demonstrating excellent predictive capability (Figure 7B and C). Decision curve analysis further indicated that the random forest model provided the greatest net clinical benefit, highlighting its robustness and practical applicability (Figure 7D). To uncover the internal decision mechanisms of the model, we applied SHAP analysis to quantify each gene’s contribution to the final prediction. Based on SHAP value rankings, the 12 feature genes were prioritized as follows: CYC1, NME1, G6PD, NAA20, TKT, GSTP1, TALDO1, SNRPB, RPL13A, MGST1, GPX3, and ALDH3A1 (Figure 7E). Notably, CYC1 exhibited the highest SHAP value (0.091) among all features, indicating that it exerted the strongest positive influence on the model’s predicted LM activity. Scatter-plot visualization of individual SHAP values further illustrated strong positive correlations between feature genes and model output probability (Figure 7F). Together, these results confirm that the random forest classifier is the most accurate and interpretable model for LM prediction, with CYC1 emerging as the principal molecular driver within the predictive framework and a compelling candidate for downstream experimental validation.

Figure 7 Benchmarking and interpretation of the optimal lactate metabolism-prediction model. A: Benchmark analysis for model comparison; B: Receiver operating characteristic (ROC) and precision-recall (PR) curves examined the accuracy of different classifiers; C: ROC and PR curves confirming random forest model accuracy [area under the curve (AUC) = 0.97, PRAUC = 0.97]; D: Decision curve analysis for clinical applicability assessment; E: SHapley Additive exPlanation (SHAP) feature importance analysis of 12 lactate metabolism-related genes; F: Correlation between SHAP values and the expression of feature genes. AUC: Area under the curve; PR: Precision-recall; TPR: True positive rate; FPR: False positive rate; SHAP: SHapley Additive exPlanation.

CYC1 knockdown inhibits the malignant activity of ICC cells

To further investigate the functional role of CYC1 in ICC, we performed loss-of-function assays in two ICC cell lines, HCCC9810 and HuCCT1. Quantitative PCR confirmed efficient knockdown of CYC1 expression in both lines following small interfering RNA transfection (Figure 8A). Cell proliferation essays using CCK-8 demonstrated that CYC1 silencing significantly reduced the growth rate of ICC cells (Figure 8B and C). Similarly, colony formation assays revealed a marked decline in both colony number and size upon CYC1 knockdown (Figure 8D and E). Furthermore, wound-healing and transwell migration/invasion assays showed that loss of CYC1 markedly impaired cellular motility (Figure 8F-H). Quantitative analyses confirmed that both migration and invasion were significantly reduced in CYC1-knockdown groups (Figure 8I and J). Collectively, these findings confirm that CYC1 promotes ICC cell proliferation, migration, and invasion, consistent with its computational identification as a core lactate-metabolism driver gene and highlighting its critical role in enhancing tumor aggressiveness.

Figure 8 Functional role of CYC1 in intrahepatic cholangiocarcinoma cell lines. A: Quantitative polymerase chain reaction validation of CYC1 knockdown efficiency in HCCC9810 and HuCCT1 cells; B and C: Cell counting kit-8 assay showing reduced proliferation following CYC1 silencing; D and E: Colony-formation assay demonstrating decreased colony number and size after CYC1 knockdown; F and G: Wound-healing assay showing impaired migratory capacity in CYC1-deficient cells; H-J: Transwell migration and invasion assays confirming that CYC1 depletion significantly suppresses motility in both intrahepatic cholangiocarcinoma cell lines. ^aP < 0.05. ^bP < 0.01. ^cP < 0.001. mRNA: Messenger RNA; OD: Optical density; NC: Negative control; Si: Small interfering RNA targeting CYC1.

DISCUSSION

ICC is an aggressive malignancy with a rapidly rising incidence and a 5-year survival rate below 10%[40,41]. Aberrant LM has emerged as a feature of tumor progression, driving angiogenesis, metabolic reprogramming, and immune evasion[42]. However, previous studies based on bulk transcriptomic data have overlooked the cellular heterogeneity underlying LM activity in ICC.

In this study, we integrated scRNA-seq and bulk transcriptomics to comprehensively delineate LM heterogeneity in ICC and identify prognostic biomarkers with potential clinical relevance. Our results reveal that LM activity was unevenly distributed across the tumor microenvironment, with malignant epithelial cells exhibiting markedly elevated metabolic activity. By integrating hdWGCNA and multiple machine learning algorithms, we identified 12 robust LM-associated feature genes (GPX3, CYC1, NME1, GSTP1, MGST1, ALDH3A1, TALDO1, SNRPB, TKT, NAA20, G6PD, and RPL13A). Among these, CYC1 emerged as the top-ranked feature by SHAP analysis and was functionally validated in vitro, underscoring its central role in promoting ICC cell proliferation, migration, and invasion.

Mechanistically, CYC1 encodes a core component of the mitochondrial cytochrome b-c1 complex (complex III), which mediates electron transfer to cytochrome c and drives oxidative phosphorylation[43,44]. Dysregulated CYC1 expression may enhance mitochondrial respiratory efficiency and elevate reactive oxygen species production, thereby promoting biosynthetic flux and redox adaptation that support rapid tumor proliferation[44,45]. Pathogenic mutations in CYC1 cause mitochondrial complex III deficiency with clinical features such as lactic acidosis, which links CYC1 function to cellular lactate handling and redox homeostasis[46]. Notably, previous reports have shown that CYC1 overexpression correlates with tumor aggressiveness and poor survival in breast cancer[47]. Our findings extend this oncogenic role to ICC, where CYC1 knockdown significantly suppressed malignant phenotypes, indicating that CYC1 functions as a metabolic nexus linking mitochondrial bioenergetics to tumor progression and may represent a promising therapeutic target in ICC.

Beyond CYC1, other top-ranked genes also contribute to metabolic reprogramming and tumor aggressiveness. NME1, a nucleoside-diphosphate kinase and well-characterized metastasis suppressor, modulates nucleotide metabolism and cell motility[48,49]. G6PD sustains NADPH production via the oxidative pentose-phosphate pathway, maintaining redox balance and contributing to chemotherapy resistance[50,51]. NAA20, the catalytic subunit of N-terminal acetyltransferase B, promotes tumorigenesis by aberrant N-terminal acetylation (e.g., of liver kinase B1), thereby attenuating adenosine 5’-monophosphate-activated protein kinase and activating mechanistic target of rapamycin signaling[52]. The convergence of these genes on redox regulation, biosynthesis, and mitochondrial activity reveals the intricate metabolic network that contributes to ICC heterogeneity. Despite emerging evidence, their distinct contributions to ICC pathogenesis remain elusive. Systematic dissection of their functional and regulatory networks is essential to understanding the metabolic and molecular basis of ICC and to developing effective targeted therapies.

Our functional enrichment analysis revealed that HLM cells exhibit upregulation of oxidative phosphorylation and ribosome biogenesis, and mitochondrial translation pathways, suggesting that lactate serves not merely as a glycolytic byproduct but also as a metabolic substrate fueling oxidative metabolism. Additionally, cancer cells can oxidize lactate via mitochondrial lactate dehydrogenase to generate pyruvate, which enters the tricarboxylic acid cycle and supports oxidative phosphorylation[53-55]. Elevated LM activity in ICC may therefore reflect a hybrid metabolic phenotype, enabling tumor cells to flexibly switch between glycolysis and oxidative phosphorylation. This metabolic plasticity facilitates adaptation to hypoxic and nutrient-limited microenvironments, contributing to tumor resilience and therapy resistance[10,56]. Additionally, our CellChat analysis revealed metabolic crosstalk between LM activity and the ICC tumor microenvironment. MIF and MDK were identified as key mediators: MIF regulates glycolysis and immune suppression, promoting M2 macrophage polarization and T-cell exhaustion, while MDK enhances tumor-stroma interactions and angiogenesis via phosphoinositide 3-kinase-protein kinase B signaling and lactate exchange[57,58]. These findings suggest that LM not only supports cancer cell energetics but also shapes the immune-metabolic landscape of the tumor microenvironment, establishing a feedback loop between metabolic adaptation and immune evasion[59,60]. Lactate accumulation profoundly affects the immune microenvironment. Elevated extracellular lactate suppresses cytotoxic T-cell and NK-cell function, promotes regulatory T-cell expansion, and fosters the differentiation of tumor-associated macrophages toward an immunosuppressive phenotype[61,62]. Recent evidence highlights the interactions between metabolic and immune signaling in shaping immunotherapy outcomes[63-65]. Integrating these findings, we propose that LM represents a central node linking metabolic reprogramming and immune modulation in ICC. Therapeutically, combining metabolic inhibitors with immune checkpoint blockade could overcome metabolic immunosuppression and improve treatment efficacy in metabolically active ICC subtypes.

Although this study provides important insights involving LM in ICC, certain limitations must be considered when interpreting the findings. While single-cell transcriptomic and in vitro functional data highlight the oncogenic role of CYC1, further in vivo studies and metabolic flux assays are required to confirm its regulatory mechanism in LM and immune modulation. Additionally, external validation is needed for other feature genes (e.g., G6PD, NME1) to strengthen the robustness of the model. The LM score and prognostic model require validation in independent cohorts and correlation with metabolic phenotyping in future study. Spatial transcriptomics and patient-derived xenograft models would also help clarify the spatial interactions between metabolic and immune networks in ICC.

CONCLUSION

In summary, our study presents an integrative single-cell and computational framework to dissect LM-related heterogeneity and molecular signatures in ICC. By linking metabolic activity to tumor microenvironment remodeling and therapeutic responsiveness, our work highlights LM as a contributing factor to ICC pathogenesis and progression, and identifies potential metabolic vulnerabilities for precision therapy. Future studies should focus on targeting LM-associated pathways and integrating metabolic inhibitors with immunotherapy to develop more effective treatment strategies.

ACKNOWLEDGEMENTS

We thank the members of Dr. Wu’s and Dr. Wang’s laboratories for their valuable technical assistance and insightful discussions. We also acknowledge the support of the Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, as well as Mayo Clinic, for providing essential experimental resources. We further appreciate the use of publicly available datasets from TCGA, GEO, and other open-access databases, which greatly contributed to this study.