The status of CA19-9 is associated with the aggressive phenotype and poor prognosis of ICC patients.

To further evaluate the clinical significance of CA19-9 in ICC, this study collected data from seven medical centers comprising 1062 ICC patients and conducted a retrospective analysis on the relationship between CA19-9 status and clinical phenotypes. Preliminary analysis revealed a specific association between CA19-9 levels, aggressive phenotype and poor prognosis in ICC patients across multiple clinical cohorts (Supplementary Fig. 1A-AK). Although some cohorts did not achieve statistical significance due to limited sample sizes, the overall trend was evident. Subsequently, we merged all patient data into a combined cohort for further analysis. The results showed that the overall survival (OS) and recurrence-free survival (RFS) of CA19-9 positive patients were significantly lower than those of negative patients (Fig. 1A, B). Moreover, compared to CA19-9 negative patients, the CA19-9 positive had significantly increased proportions of perineural invasion, vascular invasion, lymph node metastasis, and distal metastasis (Fig. 1C–F). Multivariate COX analysis further confirmed that the positive status of CA19-9 is an independent risk factor on both OS and RFS (Fig. 1G, H, Supplementary Table 5, 6). These findings support the notion that the positive status of CA19-9 is closely related to the aggressive phenotype and poorer prognosis in ICC patients.

Fig. 1

Survival analysis and aggressive characteristics of CA19-9 positive ICC patients in combined cohort. A, B The KM curves show that patients in the CA19-9 positive group exhibit poorer OS (A) and RFS (B) in the combined ICC cohort. C–F The bar plots show the proportion differences of perineural invasion (C), vascular invasion (D), lymph node metastasis (E), and distal metastasis (F) between CA19-9 positive and negative groups in the combined ICC cohort. G, H The forest plots display the multivariate Cox regression result for OS (G) and RFS (H) in the combined ICC cohort. P-values were calculated using the Log-rank test in (A, B), the Chi-square test in (C–F), and the multivariate COX regression analysis in (G, H)

The genomic characteristics between CA19-9 positive and negative ICC patients.

We ffirst compared the top 50 mutated genes in the ZS-ICC and MSK-ICC cohorts (Fig. 2A, Supplementary Fig. 2A). We found that only ten genes were shared between the two cohorts (highlighted in red), suggesting differences in the mutational landscape of different ICC populations. In both datasets, the critical driver mutation genes, including TP53, KRAS, IDH1/2, and BAP1, showed high mutation frequencies in the rankings, further confirming their essential role in ICC development. When comparing the mutation frequencies of 10 genes shared by two cohorts in the two groups, we observed that IDH1/2, PBRM1, and BAP1 mutation frequencies were generally higher in the CA19-9 negative group. In contrast, the KRAS gene is generally higher in the CA19-9 positive group. The mutation frequency distribution of the remaining six genes in the two cohorts is inconsistent (Fig. 2B, Supplementary Fig. 2B).

Fig. 2

Gene mutation pattern and tumor mutation burden of CA19-9 positive patients in the ZS-ICC cohort. A The waterfall plot shows the top 50 mutated genes in the ZS-ICC cohort. The genes marked in red are common to both cohorts. B The bar plot shows the mutation frequency differences of the ten shared high-ranked genes between CA19-9 positive and negative groups in the ZS-ICC cohort. ns: no significance. *P < 0.05. C The bar plot shows the genes with significant mutation frequency among all genes between CA19-9 positive and negative groups in the ZS-ICC cohort. Red represents a higher frequency in the positive group, while blue represents a higher frequency in the negative group. D, E The KM curves of OS between patients with and without KRAS (D) and IDH1/2 (E) mutations in the ZS-ICC cohort. F The violin plot shows the difference in tumor mutation burden between the CA19-9 positive and negative groups in the ZS-ICC cohort. P-values were calculated Using the Chi-square test in (B, C), the Log-rank test in (D, E), and the Student’s t-test in (F)

Next, we compared the differences of all mutated genes between the two groups and found that only KRAS and IDH1/2 were significant in both cohorts (Fig. 2C, Supplementary Fig. 2C). Furthermore, we investigated the potential impact of KRAS and IDH1/2 mutations on OS and RFS. The results indicate that KRAS mutations are associated with poorer OS and RFS (Fig. 2D, Supplementary Fig. 2E, F). ICC patients with IDH1/2 mutations exhibit better OS and RFS trends (Fig. 2E, Supplementary Fig. 2G,H). Considering that mutation burden is usually associated with poor prognosis of tumors, we further compared the tumor mutation burden (TMB) between the two groups. However, the differences in TMB in both cohorts were not statistically significant (Fig. 2F, Supplementary Fig. 2I). Collectively, our findings suggest that there are distribution differences of crucial driver genes between CA19-9 positive and negative ICC patients, and these phenomena may underlie the molecular basis for the clinical prognostic differences.

Transcriptome and proteomic analysis unveiled an enhanced glycolytic activity in CA19-9 positive patients

Through differential expression analysis, we identified 836, 343, and 710 significant DEGs (P < 0.05 and |logFC|> 0.5) between the CA19-9 positive and negative groups in ZS-ICC, TCGA-ICC, and GSE107943-ICC cohorts (Fig. 3A–C), respectively. Subsequently, we used GSEA to reveal the differences in Hallmark pathways between the two groups (Fig. 3D–F). An intersection analysis was conducted to pinpoint robust differentially activated pathways. The results indicated that the Hallmark glycolysis, Hallmark P53, Hallmark estrogen response late, and Hallmark allograft rejection pathways were significantly activated in all cohorts (Fig. 3G, Supplementary Table 7). In contrast, no shared downregulated pathways were obtained in the intersection analysis (Supplementary Fig. 3A). Of the four upregulated pathways, the Hallmark glycolysis pathway was associated with poorer prognosis across all three cohorts (Fig. 3H–J). In comparison, the prognostic impact of the other three pathways was less consistent (Supplementary Fig. 3B-D). Hallmark glycolysis scores were positively correlated with serum CA19-9 levels in all three cohorts (Fig. 3K–M). Considering that the transcriptomic data may not reflect the actual changes in protein level, we further analyzed the proteomic data from the ZS-ICC cohort. The enrichment analysis on the upregulated proteins showed significant enrichment in multiple glycolysis-related pathways (Fig. 3N, O), which confirmed the upregulation of the glycolytic pathway in CA19-9 positive ICC patients. In the ZS-ICC cohort, the multivariate COX analysis indicated that the Hallmark glycolysis score was an independent prognostic risk factor (Supplementary Fig. 3E, Supplementary Table 8). Furthermore, the Hallmark glycolysis scores showed a significant tumor elevation compared to normal adjacent tissues in five independent cohorts, confirming its association with malignant characteristics (Supplementary Fig. 4A-E). Notably, patients with KRAS mutations exhibited higher Hallmark glycolysis scores (Supplementary Fig. 4F), suggesting that glycolysis may serve as an intermediate factor in KRAS-mediated malignant phenotypes of ICC. In conclusion, we confirm that the glycolytic pathway was excessive activation in CA19-9 positive patients, which may be responsible for their poorer prognosis.

Fig. 3

The analysis of the transcriptome and proteome between CA19-9 positive and negative patients. A–C The volcano plots show the significantly differentially expressed genes (logFC > 0.5 & P < 0.05) between CA19-9 positive and negative groups in ZS-ICC (A), TCGA-ICC (B) and GSE107943-ICC (C) cohorts. D–F The bar plots show the differences of GSEA on hallmark pathways between CA19-9 positive and negative groups in ZS-ICC (D), TCGA-ICC (E), and GSE107943-ICC (F) cohorts. G The Venn diagram shows upregulated hallmark pathways overlapping across different ICC cohorts. H–J The KM curves show that patients in the CA19-9 positive group exhibit poorer OS in ZS-ICC (H), TCGA-ICC (I), and GSE107943-ICC (J) cohorts. K–M The scatter plots show the correlation between the serum CA19-9 levels and hypoxia glycolysis scores in ZS-ICC (K), TCGA-ICC (L) and GSE107943-ICC (M) cohorts. N The volcano plots show the significantly differentially expressed proteins (logFC > 0.5 & P < 0.05) between CA19-9 positive and negative groups in the ZS-ICC cohort. O The bubble plot shows the enrichment results on glycolysis-related pathways of the upregulated proteins (logFC > 0 & P < 0.05) in the ZS-ICC cohort. P-values were calculated using the Log-rank test in (H–J) and the Spearman correlation test in (K–M)

ScRNA-seq successfully identified three cell clusters associated with enhanced glycolytic activity in CA19-9 positive ICC patients

We turned to single-cell analysis to investigate the critical cell subpopulations responsible for the rise in glycolytic activity. After strict quality control, we obtained 97,768 high-quality cells from 20 ICC samples (Fig. 4A, Supplementary Fig. 5A, B). All cells were further divided into nine major cell clusters (Fig. 4A, B), and each cluster was defined by marker genes described in our methods section (Supplementary Fig. 5C, Supplementary Table 9). By evaluating the glycolytic score of single cells, we found that epithelial cells, stromal cells, and Mph/DC are the main clusters of glycolytic active cells (Fig. 4D, E). The scMetabolism analysis also indicated that these three cell types exhibited the highest glycolytic metabolic activity (Supplementary Fig. 5D, E). In terms of cell count, these three cell types were also numerically abundant, further supporting their significant impact on overall glycolytic activity (Fig. 4C). Notably, we observed an increase in both the number and strength of interactions between these three cell types and epithelial cells in the CA19-9 positive ICC (Fig. 4F, Supplementary Fig. 5F). Given the pivotal role of epithelial cells in tumor invasion and prognostic impact, these cell types with predominant glycolytic activities may play a vital role in the malignant progression of tumors.

Fig. 4

Single-cell analysis for all cells and epithelial cells. A The UMAP plot shows all cells colored by major cell types. B The dot plot shows the expression of the canonical marker genes across the eight major cell types. C The bar plot shows the numbers of each major cell type. D The UMAP feature plot shows the hallmark glycolysis scores in all cells. E The violin plot shows the difference in hallmark glycolysis scores in eight major cell types. F The heatmap shows eight major cell types’ interaction number differences between CA19-9 positive and negative groups. G The UMAP plot shows epithelial cells colored by seven epithelial subclusters. H The dot plot shows the expression of the marker genes across the seven epithelial subclusters. I, J The bar plot (I) and heatmap (J) show the proportion differences and Ro/e values of seven epithelial subclusters between CA19-9 positive and negative groups. *: Ro/e > 1. K The UMAP feature plots show the hallmark glycolysis scores (left) and glycolysis/gluconeogenesis scores (right) in epithelial cells. L The violin plots compare hallmark glycolysis scores (left) and glycolysis/gluconeogenesis scores (right) between Epi_SLC2A1 and other epithelial cells. M The GSEA result of hallmark glycolysis pathway between Epi_SLC2A1 and other epithelial cells. N The dot plot shows the expression of the LDHA, HK2, and PKM across the seven epithelial subclusters. O The violin plot shows the differences in GRCC scores between the CA19-9 positive and negative groups in the ZS-ICC cohort (up left panel). The KM curves compare the high and low GRCC score groups on OS in the ZS-ICC cohort (upright panel). The scatter plots show the correlation between GRCC scores, hallmark glycolysis scores (bottom left panel), and serum CA19-9 levels (bottom left panel) in the ZS-ICC cohort. P-values were calculated using the Student’s t-test in (L) and (O) (up left panel), the Log-rank test in (O) (upright panel), and the Spearman correlation test in (O) (bottom panel)

Based on the above findings, we conducted a more in-depth analysis of these three cell types. Specifically, we subdivided the epithelial cells into seven subclusters, designated as Epi_SLC2A1, Epi_APOE, Epi_IGS15, Epi_RPS27, Epi_FOS, Epi_FCER1G, and Epi_Cycling (Fig. 4G, H, Supplementary Table 10). The distribution of these subclusters varies between CA19-9 positive and negative ICC patients, with the Epi_SLC2A1, Epi_IGS15, Epi_FCER1G, and Epi_Cycling being more prevalent in CA19-9 positive ICC patients (Fig. 4I, J). Through InferCNV analysis, we identified a high level of CNV in these seven clusters, supporting their identity as malignant cells (Supplementary Fig. 5G, H). Furthermore, by calculating the glycolytic pathway score and scMetabolism analysis, we noted that the Epi_SLC2A1 exhibits the highest glycolytic activity compared to all other epithelial cell subclusters (Fig. 4K, L). Additionally, GSEA analysis revealed marked upregulation of the Hallmark glycolytic pathway in the Epi_SLC2A1 (Fig. 4M). Critically, key enzymes in the glycolytic pathway, such as LDHA, HK2, and PKM, were also expressed at significantly higher levels in the Epi_SLC2A1 compared to other epithelial subclusters (Fig. 4N), further confirming its role as the primary executor of glycolytic activity among epithelial cells. In the ZS-ICC cohort, we observed that the Epi_SLC2A1 top 50 gene signature scores were significantly higher in the CA19-9 positive group compared to the negative group. Furthermore, groups with higher Epi_SLC2A1 signature scores were correlated with poorer OS. Additionally, we found a positive correlation between the Epi_SLC2A1 signature scores, the Hallmark glycolysis scores, and the serum CA19-9 levels (Fig. 4O).

We conducted a more in-depth analysis of these three cell types based on the above findings. The stromal cells were categorized into eight subclusters (Supplementary Fig. 6A, B, Supplementary Table 11), Mph/DC cells were further divided into seven macrophage subclusters and two DC subclusters (Supplementary Fig. 6H, I, Supplementary Table 12). Through Ro/e analysis, we observed differences in the distribution of these cell subclusters among different groups (Supplementary Fig. 6C, J). Using glycolytic pathway scores and scMetabolism analysis, we screened the CAF_VEGFA and Mph-SPP1 subgroups, which exhibited significant glycolytic activity in their respective cell populations (Supplementary Fig. 6D, E, K, L), and these two clusters were more likely to be distributed in CA19-9 positive ICC patients. The results of the GSEA analysis further confirmed the significant upregulation of glycolytic pathways in these two clusters (Supplementary Fig. 6F, M). In addition, both cell clusters highly express critical genes involved in glycolysis, including HK2, PKM, and SLC2A1 (Supplementary Fig. 6G, N). In the ZS-ICC cohort, the top 50 gene signature scores for these two cell clusters were significantly higher in CA19-9 positive patients, and the group with high signature scores was significantly associated with poorer prognosis. We also found a positive correlation between the signature scores of the two clusters, glycolysis pathway scores, and serum CA19-9 levels (Supplementary Fig. 6O, P).

We performed mIF on 26 ICC samples collected at our center to validate our analytical findings. The staining results indicated a higher cellular density of the Epi_SLC2A1, CAF_VEGFA, and Mph_SPP1 cell subclusters in the CA19-9 positive patient group (Fig. 5A–C). Statistical analysis further confirmed that the positive group’s cell count per unit area was significantly increased compared to the negative group (Fig. 5D), which is consistent with the in silico results. Collectively, we successfully identified cell subgroups with enhanced glycolytic activity in epithelial cells, stromal cells, and Mph/DC cells. These subgroups are more prevalent in CA19-9 positive ICC patients and are associated with poorer prognosis.

Fig. 5

Glycolysis-related cell subclusters were verified using mIF and spatial transcriptomics analysis. A–C The representative mIF images of Epi_SLC2A1 (A, blue: DAPI, green: GLUT1, red: CK19), CAF_VEGFA (B, blue: DAPI, green: VEGFA, red: αSMA), and Mph_SPP1 (C, blue: DAPI, green: SPP1, red: CD68) in CA19-9 positive (bottom panel) and negative (up panel) groups. Scale bar = 50um. D The box plots show the differences in counts/unit area in Epi_SLC2A1 (left panel), CAF_VEGFA (mid panel), and Mph_SPP1 (right panel) between CA19-9 positive and negative groups (positive: n = 14, negative: n = 12). E The UMAP (up left panel) and spatial (up right panel) plots show all spots colored by tissue types. The spatial feature plots show the expression of KRT19 (bottom left panel) and ALB in all spots (bottom left panel). F The spatial feature plots show the Epi_SLC2A1 (left panel panel), CAF_VEGFA (mid panel panel), and Mph_SPP1 (right panel panel) signature scores in tumor region spots. G The correlation analysis between the Epi_SLC2A1 and CAF_VEGFA signature scores, between the Epi_SLC2A1 and Mph_SPP1 signature scores, and between the CAF_VEGFA and Mph_SPP1 signature scores in tumor region spots. H The violin plots show that the CA19-9 positive group has higher GRCC scores in the ZS-ICC cohort (up left panel). The KM curves show that the high GRCC scores group has a poorer OS in ZS-ICC cohorts (upright panel). The scatter plots show that the GRCC scores significantly correlate with hallmark glycolysis scores (bottom left panel) and serum CA19-9 levels (bottom right panel) in the ZS-ICC cohort. I The forest plot shows the multivariate Cox regression result for OS in the ZS-ICC cohort. P-values were calculated using the Mann-Whitney U test in (D), the Pearson correlation test in (G), the Student’s t-test in (H) (up left panel), the Log-rank test in (H) (upright panel), and the Spearman correlation test in (H) (bottom panel).

Glycolysis-related cellular communities are associated with the aggressive phenotype and adverse prognosis in CA19-9 positive ICC

We utilized spatial transcriptome technology to explore the spatial distribution and interactions of the three cell subclusters within the tumor microenvironment. Firstly, we divided the spots into two main clusters and classified them as tumor and pera-tumor regions based on the expression patterns of KRT19 and ALB (Fig. 5E). We separately assessed the infiltration scores of three glycolysis-related cell subclusters within the tumor region spots. The result revealed a pattern of co-localization for these subclusters within the tumor area (Fig. 5F). Furthermore, by analyzing the correlation of cell scores in each spot, we observed significant positive correlations among these cell subclusters (Fig. 5G). This finding suggests that these subclusters may be co-recruited within the tumor microenvironment and form an interactive cellular network. To further validate the association of these cellular communities with tumor progression, we constructed a GRCC signature by combining the top 50 genes from each of the three cell subclusters. In the ZS-ICC cohort, we found that the GRCC signature score was higher in CA19-9 positive patients and was significantly correlated with poorer prognosis. Additionally, the GRCC score showed a significant positive correlation with glycolysis pathway scores and serum CA19-9 levels (Fig. 5H). In multivariate COX regression analysis, a high GRCC signature score was identified as an independent risk factor affecting OS (Fig. 5I). Encouragingly, we observed broadly consistent results in both the TCGA-ICC cohort and the GSE179423-ICC cohort, indicating that our findings are reproducible and generalizable (Supplementary Fig. 7A, B). In groups with high GRCC signature scores, the incidence of aggressive phenotypes such as intrahepatic metastasis, neural invasion, vascular invasion, lymph node metastasis, and distal metastasis was also higher (Supplementary Fig. 7C), further emphasizing the close link between these glycolytic cellular communities and tumor aggressiveness. Finally, in five independent cohorts, we confirmed that the GRCC signature score in the tumor area was significantly higher than in the peritumoral area (Supplementary Fig. 7D). In summary, our results reveal that CA19-9 positive ICC patients tend to form a greater abundance of glycolysis-related cellular communities, which may lead to a more aggressive tumor phenotype and a poorer prognosis.

Hypoxia induces the formation of glycolysis-related cellular communities by promoting metabolic reprogramming

To investigate the origin and developmental trajectories of these three glycolytic-related cell clusters, we constructed evolutionary trajectories in each major cell type using the Monocle2 algorithm. In the epithelial cells, we identified five different cell states that form two distinct evolutionary branches (Lineage1 and Lineage2, Supplementary Fig. 8A). Similarly, we identified five and seven cellular states in the stromal cells and macrophages, forming three evolutionary branches (Lineage1, Lineage2, and Lineage3, Supplementary Fig. 8B, C). We observed that the three glycolysis-related cell clusters were mainly located at the end of the pseudotime trajectory. To further clarify the specific branches of these three clusters during evolution, we assessed the proportionate changes of each subcluster across different states. In epithelial cells, the Epi_SLC2A1 showed an increasing trend across both Lineage1 and Lineage2 branches (Supplementary Fig. 8D). In stromal cells, the CAF_VEGFA exhibited an increasing trend only in Lineage3, while in macrophages, the Mph_SPP1 showed an increasing trend solely in Lineage2 (Supplementary Fig. 8E, F). These results suggest that all two branches of epithelial cells represent the evolutionary trajectories of Epi_SLC2A1, while CAF_VEGFA and Mph_SPP1 need to follow specific evolutionary branches in their major cell types. We further assessed the hallmark pathways along their evolutionary trajectories. As expected, the glycolysis pathway was significantly enriched in the terminal states of all trajectories (Supplementary Fig. 9A-C). Concurrently, the oxidative phosphorylation (OXPHOS) pathway was more prevalent in the initial states, indicating a metabolic shift from OXPHOS to glycolysis—a common feature across these three cellular evolutionary trajectories. In trajectory plots, we distinctly tracked a decreasing trend in OXPHOS pathway scores and an increasing trend in glycolytic pathway scores along the trajectories, exhibiting a markedly opposite trend (Fig. 6A–C). Correlation analysis further confirmed significant associations between pseudotime and the signature scores of these metabolic pathways, particularly pronounced in the CA19-9 positive group (Fig. 6D–F). In GSEA conducted across three cohorts, we similarly observed opposing activation patterns between the OXPHOS and glycolysis pathways (Supplementary Fig. 8G). Integrating the above findings, our study reveals that these three glycolysis-related cell clusters may have evolved through metabolic reprogramming.

Fig. 6

The analysis of evolutionary trajectories and cellular interaction network. A-C. Trajectory plots of epithelial cells (A), stromal cells (B), and macrophages (C) colored by hallmark OXPHOS (left panel) and glycolysis (right panel) scores. D–F The correlation analysis between pseudotime and hallmark OXPHOS (left panel) and glycolysis (right panel) scores in epithelial cells (D), lineage 3 of stromal cells (E), and lineage 2 of macrophages (F). G The correlation analysis between pseudotime and hallmark hypoxia scores in epithelial cells (left panel), lineage 3 of stromal cells (mid panel), and lineage 2 of macrophages (right panel). H The scatter plots show the correlation between the hallmark hypoxia scores of each sample and the proportion of Epi_SLC2A1 (left panel), CAF_VEGFA (mid panel), and Mph_SPP1 (right panel). I The representative IHC images of the HIF1αstaining in CA19-9 positive (up panel) and negative group (bottom panel). Scale bar = 100um. J The box plot shows the differences in the fraction of HIF1αstaining between CA19-9 positive and negative groups (positive: n = 31, negative: n = 22). K The heatmaps of Nichenet analysis show regulatory patterns between CAF_VEGFA to Epi_SLC2A1 (up panel) and Mph_ SPP1 to Epi_SLC2A1 (bottom panel). L The bar plots show the GO and KEGG pathways enrichment results of the predicted target genes in Epi_SLC2A1. P-values were calculated using the Student’s t-test in (J) and the Spearman correlation test in (D–H)

ICC is characterized by its scarcity of blood vessels, which frequently results in hypoxic conditions within the tumor microenvironment. In the previous section, we observed an increasing trend of hypoxia signature in all three trajectories (Supplementary Fig. 9A-C). Given the close relationship between OXPHOS and glycolysis pathways with hypoxia, we hypothesized that hypoxia might be a potential factor driving metabolic shifts and ultimately promoting the formation of glycolysis-related cell clusters. To prove this hypothesis, we performed a correlation analysis between hypoxia pathway scores and pseudotime, showing a positive correlation with three glycolysis-related cellular evolutionary trajectories (Fig. 6G). Furthermore, the sample's hypoxia score was positively correlated with the proportions of Epi_SLC2A1, CAF_VEGFA, and Mph_SPP1 (Fig. 6H). Notably, samples with higher hypoxia scores were more common in the CA19-9 positive group, although this difference did not reach statistical significance due to the limited sample size (Supplementary Fig. 8H). In the ZS-ICC cohort, we also observed a positive correlation between hypoxia, glycolysis, and GRCC scores, and a negative correlation with OXPHOS scores (Supplementary Fig. 8I-K). Additionally, we performed IHC staining for the hypoxia marker HIF1a on tissue samples from 54 ICC patients. The results indicate that the CA19-9 positive group exhibits a higher fraction of HIF1a staining, and the difference was statistically significant (Fig. 6I, J). With increased HIF1a staining fraction, Epi_SLC2A1, CAF_VEGFA, and Mph_SPP1 were also increased (Supplementary Fig. 10A, B). Collectively, we observed a more hypoxic tumor microenvironment in CA19-9 positive ICC patients, which triggered more active metabolic reprogramming, leading to an increased accumulation of glycolysis-related cell communities.

Deciphering the cellular interaction network among glycolysis-related cellular communities

To further explore the interactions between glycolysis-related cell clusters, we used the NicheNet algorithm to analyze the ligand-receptor (LR) interaction networks between these cells. Intriguingly, the enrichment analysis on the predicted targeted genes of the three cell clusters all showed significant enrichment in oxygen level responses and glycolysis (Fig. 6L, Supplementary Fig. 11B, D). This discovery suggests that the interactions between the three clusters play an essential role in their development and functional maintenance. In the LR interaction between CAF_VEGFA and Mph_SPP1 and Epi_SLC2A1, we found several molecules among the high activity ranking ligands with crucial roles in cell proliferation and migration, such as TGFB1, CXCL12, FGF11, and GDF15 (Fig. 6K). Additionally, the enrichment analysis on the predicted targeted genes of Epi_SLC2A1 focused on biological processes like epithelial cell migration, cell growth, and cell adhesion (Fig. 6L), indicating that the other two cell clusters may regulate the EMT of Epi_SLC2A1. Further analysis showed that Epi_SLC2A1 has a marked upregulation of the EMT pathway and a higher EMT score compared to other epithelial cells (Fig. 7A, B), and the EMT score of epithelial cells is also highest at the end of the evolutionary trajectories, and is positively correlating with pseudotime (Fig. 7C). Using the CellChat analysis, we observed an increase in the number and intensity of communication between the three glycolysis-related cell types in CA19-9 positive ICC, and signaling pathways closely associated with EMT, such as the TGFb, ncWNT, and EGF pathways, exhibited significant upregulation (Fig. 7D, Supplementary Fig. 12A, B). Notably, when using Epi_SLC2A1 as the receptor, the TGFB1 ligand signal was significantly enhanced in the three cell types of the CA19-9 positive group (Fig. 7E). In addition, the activity of TGFB1 also ranks high in LR interaction analysis with CAF_VEGFA and Mph-SPP1 (Supplementary Fig. 11A, C), indicating that TGFB1 plays a crucial role in regulating the interaction network between glycolysis-related cell clusters. For CAF_VEGFA, the predicted targeted genes were mainly enriched in tissue remodeling, cell–matrix adhesion, and mesenchymal cell differentiation (Supplementary Fig. 11B), essential for tumor metastasis [49,50,51]. For Mph_SPP1, the enriched terms were related to lipid transport, lipid storage, and myeloid cell activation (Supplementary Fig. 11D), which have been implicated in M2 polarization of macrophages in previous studies [52,53,54]. In summary, our findings indicate that the crosstalk among the three glycolysis-related cell clusters constitutes a vicious cycle that promotes their own formation and maintenance plays a crucial role in tumor metastasis and invasion.

Fig. 7

The effects of glycolysis-related cellular communities on EMT, angiogenesis, and drug screening strategy. A The GSEA result of hallmark EMT pathway between Epi_SLC2A1 and other epithelial cells. B The UMAP feature (up panel) and violin (bottom panel) plots show the Epi_SLC2A1 has higher hallmark EMT scores compared to other epithelial cells. C The trajectory plot (up panel) shows the hallmark EMT scores and the correlation analysis (bottom panel) between pseudotime and hallmark EMT scores in epithelial cells. D The circle plots show the differences in the number (up panel) and intensity (bottom panel) of interactions within the three glycolysis-related cell subclusters between CA19-9 positive and negative groups. Red represents increasing, and blue represents decreasing in the positive group. E The dot plot shows the differences in TGFB1 signaling interaction of three glycolysis-related cell subclusters to Epi_SLC2A1 between CA19-9 positive and negative groups. F The dot plots show the expression of the VEGFA across the subclusters of stromal cells (left panel), epithelial cells (mid panel), and macrophages (right panel). G The dot plot shows the differences in VEGFA signaling interaction of three glycolysis-related cell types to endothelial cells between CA19-9 positive and negative groups. H The UMAP plot shows endothelial cells colored by seven endothelial subclusters. I The dot plot shows the differences in VEGFA signaling interaction between the three glycolysis-related cell subclusters and seven endothelial subclusters between CA19-9 positive and negative groups. J The overall workflow for screening candidate drugs targeting GRCC and CA19-9 positive ICC patients. P-values were calculated using the Student’s t-test in (B) (bottom panel) and the Pearson correlation test in C (bottom panel)

Glycolysis-related cellular communities promote tumor angiogenesis

The enrichment analysis on the predicted target genes of the three glycolysis-related cell clusters all showed positive regulation on vascular development (Fig. 6L, Supplementary Fig. 11B, D), suggesting that these cells may have the function of promoting tumor angiogenesis. VEGFA, a factor crucial for angiogenesis, activates endothelial cell proliferation by binding to its receptors on the surface of endothelial cells. In our data, we noted generally higher expression levels of VEGFA in these glycolysis-related cell clusters (Fig. 7F). In the ZS-ICC cohort, samples with high GRCC scores exhibited higher VEGFA expression, with a significant positive correlation between the two (Supplementary Fig. 12C). Moreover, CellChat analysis demonstrated a significant increase in VEGFA ligand signing between these three glycolysis-related cells and endothelial cells compared to other cells (Supplementary Fig. 12D). As expected, the communication probability between these three cell clusters and endothelial cells showed a general increase in the CA19-9 positive ICC patients (Fig. 7G). Subsequently, we subclassified endothelial cells into seven distinct subclusters (Fig. 6H, Supplementary Fig. 11E, F, Supplementary Table 13). We found that the EC_Tip_ESM1 had the closest interactions in VEGFA signing with the three glycolysis-related cells (Fig. 6I), suggesting it may be the critical subcluster influenced by these cells. In three independent cohorts, the top 50 signature score of the EC_Tip_ESM1 showed a positive correlation with GRCC scores, and KM analysis also indicated that the higher End_Tip_ESM1 score group tended to have poorer OS (Supplementary Fig. 11G, H). In summary, we found that high expression levels of VEGFA often accompany glycolysis-related cell clusters. VEGFA promotes endothelial cell proliferation and angiogenesis by activating receptors on endothelial cells. This process provides more blood supply to the tumor, which may enhance its invasiveness and metastatic potential of the tumor.

Therapeutic drug screening for GRCC and CA19-9 positive ICC patients

This study employed a two-stage strategy to screen therapeutic drugs for GRCC and CA19-9 positive ICC patients (Fig. 7I). Through this method, we identified 28 CTRP-derived drugs targeting GRCC and 12 drugs specifically suitable for CA19-9 positive ICC patients. After intersection analysis, we selected six drugs that yielded significantly lower AUC values in GRCC and CA19-9 positive ICC patients and showed negative correlations with the GRCC scores (Supplementary Fig. 13A–F). Among these candidate drugs, Selumetinib and Trametinib were the main components, with the remaining four drugs being combined applications of these two agents. By querying the CTRP database, we discovered that both are MEK inhibitors. As a pivotal regulatory factor in the MAPK signaling pathway, MEK plays a crucial role in controlling essential biological processes such as cell proliferation, differentiation, migration, and survival. Interestingly, we observed a significant enrichment of the MAPK signaling pathway in the enrichment analysis on the predicted target genes of Epi_SLC2A1 (Fig. 6L). This result further corroborates the potential application value of these candidate drugs.