Deciphering cholangiocarcinoma heterogeneity and specific progenitor cell niche of extrahepatic cholangiocarcinoma at single-cell resolution

Abstract

Background: Cholangiocarcinoma (CCA) is a highly heterogeneous malignancy, primarily comprising intrahepatic (iCCA) and extrahepatic (eCCA) subtypes. Reconciling the variability between iCCAs and eCCAs in clinical trials remains a challenge, largely due to the inadequate understanding of their shared and subtype-specific cellular heterogeneity. We aim to address this issue using single-cell and spatially resolved transcriptomic approaches.

Methods: We performed comprehensive single-cell RNA sequencing (scRNA-seq) by profiling 109,071 single cells from 28 samples, including chronic biliary inflammatory conditions (n = 7) and CCAs from different anatomical sites (n = 21). Findings were validated using external multi-omics datasets, tissue microarray cohort, spatial RNA in situ sequencing, CCA patient-derived organoids (PDOs), and mouse models.

Results: iCCAs and eCCAs exhibited distinct tumor ecosystems, with notable differences in cellular composition, diversity, and abundance across various cell types. Non-malignant epithelial cells displayed divergent precancer hallmarks from different biliary sites, with inflammatory extrahepatic bile ducts exhibiting early hijacking of the gastrointestinal metaplastic process. We identified seven meta-programs within cancer cells, mapped into four major subtypes. This subtyping was validated using external CCA cohorts and PDO models, distinguishing patients based on clinical outcomes and drug vulnerabilities. Specifically, iCCAs were associated with a senescent program, while eCCAs were enriched in an IFN-responsive program linked to adverse clinical outcomes and increased drug resistance. We identified a basal-like LY6D⁺ cancer cell subpopulation specific to eCCAs, which displayed significant stemness, drug resistance, and IFN-responsive features. This subpopulation was closely associated with an interferon-stimulated gene 15 (ISG15)-enriched mesenchymal and immune microenvironment. Functional assays demonstrated that ISG15 stimulation significantly boosted stemness, basal-like features, and drug resistance in eCCA cells, highlighting its pivotal role in sustaining the LY6D⁺ progenitor niches.

Conclusion: We present a comprehensive single-cell landscape of CCAs, uncovering the molecular heterogeneity between iCCA and eCCA subtypes. Transcriptomic subtyping of CCA cancer cells offers implications for clinical stratification and functional precision oncology. We identify basal-like epithelial progenitors and characterize their associated ISG15-enriched microenvironment in eCCAs. These findings hold significant promise for the development of novel prognostic biomarkers, therapeutic targets, and treatment strategies for CCAs.

Fig. 1

The single-cell landscape of included patients. (A) Experimental design and workflow of our study. (B) UMAP visualization of all cells from included patients colored by cell types. (C) Violin plot showing marker genes expression levels of each cell type. Numbers on the left side represent the cell types shown in (B). (D) UMAP plots showing representative markers. (E)The proportions of nine cell types in each patient. Cell types are colored the same as (B). (F) Sample type preference of various cell types analyzed by Ro/e. Numbers on the right side represent cell types

Fig. 2

Decoding shared and distinct malignant meta-programs between CCA subtypes. (A) NMF analysis of mEPCs. Heatmap is colored by Jaccard similarity. Numbers on the top represent seven meta-programs. (B) Heatmap showing gene signatures in seven meta-programs normalized by z-score. Meta-programs are labeled the same as (A). (C) Bubble plot showing enrichment of GO biological processes across seven meta-programs. (D) Violin plots comparing meta-program enrichment scores between different CCA subtypes. (E) Heatmap showing four subtypes of in-house CCA scRNA-seq cohort based on meta-program scores. LN, lymph node; MVI, microvascular invasion. (F-I) Representative IHC staining images of markers for four malignant subtypes of our CCA tissue microarray. (J-M) Kaplan-Meier curves comparing overall survival of four subtypes markers in our CCA tissue microarray cohort. The p values are from Wilcoxon test. *, p < 0.05; **, p < 0.01; ***, p < 0.001

Fig. 3

CCA subtypes recapitulated in PDO models linking to drug response. (A) Representative bright field and H&E staining images of 16 CCA PDOs. (B) Heatmap showing subtyping of 16 PDOs based on NMF meta-program scores. (C) Heatmaps displaying high-throughput screening of 175 drugs across four subtypes of PDOs, with the right and left side heatmap showing cell viability and normalized viability, respectively. Barplots on the left side of heatmaps showing scaled sums of z-scores representing the overall drug response. Drugs are colored by related pathways and are clustered by rows. (D-G) Barplot showing top 30 effective drugs in proliferative-subtype, senescent-subtype, IFN responsive-subtype and glandular-subtype PDOs

Fig. 4

Identification of eCCA-specific basal-like progenitor populations. (A) UMAP plot of mEPCs colored by clusters. (B) Pie chart showing the patient composition of cluster 7. (C) Scatter plot displaying upregulated marker genes of cluster 7. (D-F) UMAP density plot showing LY6D, KLK7 and CDH3 expression across mEPCs. (G) UMAP plot showing differential transcriptional entropy of mEPCs calculated by CytoTRACE. (H-I) Pseudotime analysis of mEPCs from eCCA-P and eCCA-D patients. Cells are colored by Seurat clusters and sample types. (J) The dynamic expression of marker genes of cluster 7 along the pseudotime. (K) Bubble plot showing enrichment of GO biological processes using cluster 7 DEGs

Fig. 5

eCCA-specific basal-like progenitors showing basal cell and stemness features, with increased drug resistance. (A) Multiplex IHC of LY6D⁺ EPCs in three eCCA patients, indicted by white arrows. DAPI: deep blue. Pan-CK: light yellow. LY6D: pink. Scale bar: 10 μm. (B) RNA in situ sequencing of three eCCA patients. LY6D⁺ and LY6D⁻ EPCs are labeled by red and blue, respectively. Cell nuclei are labeled by white. (C) RNA in situ sequencing revealing differential expression of basal, stemness and proliferative marker genes between LY6D⁺ and LY6D⁻ EPCs. (D) Sorting of LY6D⁺ and LY6D⁻ subpopulations in QBC-939 eCCA cell line by flow cytometry. (E) Bar plot comparing RT-qPCR results of marker genes between LY6D⁺ and LY6D⁻ subpopulations of the QBC-939 eCCA cell line. (F) Barplot comparing spheroid formation capacity between LY6D⁺ and LY6D⁻ subpopulations of the QBC-939 eCCA cell line. (G) Dose-response curves showing differential responses to gemcitabine, cisplatin, and 5-fluorouracil between LY6D⁺ and LY6D⁻ subpopulations of QBC-939 eCCA cell line. (H) Photos comparing subcutaneous tumor growth after injection of LY6D⁺ and LY6D⁻ populations from the QBC-939 eCCA cell line in nude mice. (I) Bar plot comparing average tumor volumes after injecting LY6D⁺ or LY6D⁻ cells. (J) Western blot analyses comparing protein expression of LY6D and stemness markers between LY6D⁺ and LY6D⁻ tumors. Statistical analyses are from t-test. *, p < 0.05; **, p < 0.01; ***, p < 0.001

Fig. 6

scRNA-seq analysis of CCA mesenchymal microenvironment. (A) UMAP showing the mesenchymal cells colored by cell subtypes. iFib, inflammatory fibroblast; senCAF, senescent CAF; iPer, inflammatory pericyte; myoPer, myofibroblast pericyte. (B) UMAP showing mesenchymal cells colored by sample types. (C) Sample preference of the mesenchymal subtypes analyzed by Ro/e. (D) Heatmap showing average expression of marker genes across eight mesenchymal subtypes. (E) Barplots showing enrichment of GO biological processes across eight mesenchymal subtypes. (F) Pie chart displaying the percentage contribution of sample sources in ISG15 Per. (G) Bubble plot showing cell-cell interaction of EDN signaling pathways between eCCA mEPC subclusters and Per_ISG15

Fig. 7

scRNA-seq revealing IFN-signaling activated immune subpopulations. (A) UMAP of myeloid cells colored by subtypes. Mono, monocytes; Macro, macrophages; DC, dendritic cells. (B) Bubble plots showing the expression of marker genes across myeloid cell subtypes. (C) Sample preference of myeloid cells subtypes analyzed by Ro/e. (D) Heatmap showing functional enrichment scores across monocyte and macrophage subtypes. TAM, tumor-associated macrophage. (E) Heatmap showing functional enrichment scores across different sample types. (F) Pie chart displaying the percentage contribution of sample sources in Macro_ISG15. (G) Kaplan-Meier curves of overall survival in ICGC-eCCA cohort patients grouped by Macro_ISG15 signature score. (H) UMAP visualization of lymphocytes colored by subtypes. (I) Bubble plot showing the expression of marker genes of lymphocyte subtypes. (J) Sample preference of lymphocyte subtypes analyzed by Ro/e. (K, L) Heatmap showing functional enrichment scores of CD4⁺ T cell and CD8⁺ T cell subtypes. (M) Pie chart displaying the percentage contribution of sample sources in CD8T_ISG15. (N) Kaplan-Meier curves of overall survival in ICGC-eCCA cohort patients grouped by CD8T_ISG15 signature score

Fig. 8

Specific mesenchymal and immune niches sustaining basal-like progenitors. (A) Circular plot displaying the distribution of sample sources for Macro_ISG15, CD8T_ISG15, Per_ISG15 and mEPC_c7: patient sources (left) and sample type origins (right). (B) Heatmap showing average expression of ISG15 and representative upstream genes across mesenchymal and immune cell subtypes. (C) Boxplot comparing ISG15 expression levels between iCCAs and eCCAs in the ICGC cohort. (D) Correlation heatmap analyzing expression of ISG15 and representative DEGs of c7 mEPC clusters in ICGC eCCA patients. (E) Kaplan-Meier curves of overall survival in eCCA patients stratified by ISG15 expression in the ICGC-eCCA dataset. (F) Spatial RNA ISS visualization of LY6D⁺ EPCs and ISG15-expressing pericytes, macrophages and T cells. LY6D⁺ EPCs: red. Pericytes: light brown, Macrophages: light blue. T cells: grey. ISG15: green cross. (G, H) Spheroid formation capacity analysis between vehicle- and ISG15 (100 ng/mL)-treated QBC-939 eCCA cell line. (I) Bar plot comparing RT-qPCR results of marker genes between vehicle- and ISG15-treated QBC-939 eCCA cell line. (J) WB analysis comparing proteins levels of LY6D, KLK7 and stemness markers between vehicle- and ISG15 (100 ng/mL)-treated QBC-939 eCCA cell line. Statistical analyses are obtained from t-test. *, p < 0.05; **, p < 0.01; ***, p < 0.001