Role of cancer stem cell heterogeneity in intrahepatic cholangiocarcinoma

Introduction

Intrahepatic cholangiocarcinoma (ICC), originated from epithelial cells of intrahepatic and extrahepatic bile ducts, is the second leading hepatic malignancy, characterized by high invasiveness and refractory to chemotherapy (1). Over the past decade, ICC has become a major global health concern due to its increasing morbidity and mortality, compounded by high recurrence and metastasis rates after treatment (2-4). According to statistics, only 5% of ICC patients survive more than 5 years, and satisfactory solutions for this low survival rate have not been established (5).

A major factor in ICC’s aggressive behavior is the presence of cancer stem cells (CSCs), a unique subset of tumor cells capable of self-renewal and differentiation. Unlike normal stem cells, CSCs evade regulatory mechanisms, allowing unchecked expansion and continuous production of tumorigenic progeny, which drive recurrence and metastasis in ICC (6,7). Critical signaling pathways play essential roles in maintaining CSC characteristics within ICC. The Notch pathway, for instance, is known to regulate CSC self-renewal and differentiation in liver cancers, where its activation supports tumorigenic properties and contributes to ICC progression and therapy resistance (8). Similarly, Wnt/β-catenin signaling is crucial for maintaining CSC stemness, and its abnormal activation in ICC is associated with increased invasiveness and therapeutic resistance, making it a promising target for disrupting CSC-driven tumor growth (9), positioning CSCs as a crucial focus in ICC therapy. And it is reported that targeting CSC can significantly reduce metastasis and recurrence rates in other cancers, including gastric cancer and breast cancer (10,11).

Recent advances in systemic therapies, including immunotherapy and targeted molecular treatments, offer promising options for addressing CSC-driven ICC progression. For instance, novel immunotherapies are being explored to enhance the immune response against ICC by targeting immune evasion mechanisms (12). Targeted therapies directed at specific mutations, like fibroblast growth factor receptor (FGFR) and isocitrate dehydrogenase (IDH), are also promising as they focus on molecular alterations unique to ICC tumors (13). However, these treatments come with risks; for example, immune checkpoint inhibitors (ICIs) have been linked to liver-specific toxicities, including elevated transaminase levels, which require careful monitoring, especially in ICC patients (14). This highlights the need for a refined therapeutic approach that balances efficacy with safety. However, ICC CSC populations are highly heterogeneous, with distinct subtypes exhibiting variable molecular characteristics and responses to treatment, making therapeutic targeting complex.

The advent of single-cell RNA sequencing (scRNA-seq) has transformed our understanding of CSC heterogeneity by providing an unprecedented resolution of the cellular composition within tumors. Unlike conventional bulk RNA sequencing, which aggregates gene expression across all cell types, scRNA-seq allows for a cell-by-cell analysis, uncovering the diversity of CSC subpopulations within ICC and their unique molecular profiles. This technology has already been applied in several cancers, such as esophageal squamous cell carcinoma and ovarian cancer, where it has enabled the identification of prognostic markers and therapeutic targets within CSC populations (15,16). In ICC, scRNA-seq has shed light on interactions between CSCs and cancer-associated fibroblasts (CAFs) (17,18), revealing the immune status of the ICC tumor microenvironment, and providing new biomarkers for ICC prognosis (19). However, existing studies on ICC have largely overlooked the detailed heterogeneity of CSCs, often lacking the resolution to identify distinct subpopulations and their unique roles in immune evasion and drug resistance, and no research has yet explored the relationship between CSC heterogeneity and patient prognosis in ICC. Here, we provide a comprehensive analysis of CSC subgroups (CD13⁺, CD90⁺ and CD133⁺), their signaling pathways, and interactions within the tumor microenvironment based on scRNA-seq data. By integrating this heterogeneity into a prognostic model, we aim to enhance outcome prediction and identify targeted therapeutic strategies, advancing personalized treatment approaches for ICC. We present this article in accordance with the REMARK reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1286/rc).

Methods

Data collection

Single-cell transcriptome data for CSCs in ICC and adjacent tissue samples were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), accession number GSE142784 (17), including four treatment-naive ICC samples, one recurrent ICC sample, and three adjacent tissues; clinical sample information and gene expression data for ICC were obtained from The Cancer Genome Atlas (TCGA) database (https://www.ncbi.nlm.nih.gov/geo/).

scRNA-seq data processing

R language scripts were written for analysis of the obtained scRNA-seq data. First, sequencing data were merged via the “Seurat” R package, and then log-normalization of the merged datasets was performed using the “NormalizeData” function. Next, the 2,000 most highly variable genes were selected after quality control using “FindVariableFeatures”. Subsequently, the “ScaleData” function was employed to further correct the data; the t-distributed stochastic neighbor embedding (t-SNE) algorithm was used to cluster cell groups; the “RunTSNE” function was used to generate clusters; the “FindAllMarker” function was used to annotate each cluster; and the “FindMarkers” function was used to identify differentially expressed genes between two clusters. In addition, the “Seurat” and “DoHeatmap” functions were used to assess features and generate heatmaps of gene expression, respectively, and SCENIC was used to further analyze changes in the expression of transcription factors. Clusters consisting of CSCs were extracted and reprocessed in the same way as described above, and each CSC type was further divided into subclusters. Marker genes for each stem cell subset were identified by comparison of the CSC subsets with the corresponding normal subsets, and P<0.05 was considered to indicate a significant difference. Finally, marker genes for each CSC subgroup were regarded as differentially expressed genes.

Pseudotime analysis of single-cell gene expression

Pseudotime analysis was performed using Monocle2 (version 2.18.0), and the assessed time was scaled from 0 to 1 (20). The indicated channel acted as the input dimension. Then, hub genes in each cluster were identified based on the “diffialGeneTest” function in the “Monocle2” package. According to the q value (q<0.01), the top 2,000 genes were filtered out. The expression profiles were reduced to two dimensions using the DDRTree method through the “reduceDimension” function in Monocle2 (max_components =2). Then, the “orderCells” function was applied to order malignant cells and assign a “time” value.

Calculation of cell signaling pathway activity

Hallmark gene sets of signaling pathways, including the epithelial mesenchymal transition (EMT) pathway, P53 pathway, Notch signaling pathway and Wnt/β-catenin pathway, were obtained from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/). Briefly, key genes involved in metastasis were obtained using MSigDB. Subsequently, CellChat was used to assess the viability of individual cells according to each gene set.

Cell-cell communication

To investigate cell-cell communication and interactions and identify mechanisms of communication at single-cell resolution, the R package “CellChat” (version 1.0.0) was applied to cells in five cell groups, namely, CD13 cells, CD133 cells, CD90 cells, dual leucine zipper kinase (DLK)1 + ICC cells, and ICC cells. Specifically, CellChat was used to construct cell-cell communication networks and visualize the signals of each cell group. Additionally, the “computeNetSimilarity” function was adopted for identifying signal groups with similar functions or structures, and the “netAnalysis_signalingRole_heatmap” function was used for visualizing signaling pathways.

Immunohistochemical staining

Specimens from 10 treatment-naive ICC patients were collected from the Third Affiliated Hospital of Sun Yat-sen University. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by the Institute Research Ethics Committee of the Third Affiliated Hospital of Sun Yat-sen University (No. [2023]02-167) and informed consent was taken from all individual participants.

The EnVision immunohistochemistry (IHC) system was used to analyze the expression of CD13, CD133, CD90, HSF1 and STAT1. The slides were incubated with rabbit-derived monoclonal antibodies against CK19 (ab76539, Abcam, Cambridge, USA, 1:1,000), CD90 (ab92574, Abcam, 1:200), CD13 (GB11356, Servicebio, Wuhan, China, 1:500), CD133 (ab222782, Abcam, 1:1,000), HSF1 (GB11932, Servicebio, 1:500) and STAT1 (GB111363, Servicebio, 1:1,000) as the primary antibodies and then incubated with a horseradish peroxidase (HRP)-labeled secondary antibody (Dako Envision, Glostrup, Denmark). Phosphate buffer saline (PBS) was used as the negative control.

Prediction of prognostic features

After collecting cholangiocarcinoma-related datasets from the TCGA, bulk RNA-seq data from the ssGSEA tool was used for data scoring, ConsensusClusterPlus was used for unsupervised clustering, and the entire dataset was used for observing the predictive ability of cell subgroups. In addition, Kaplan-Meier curves were generated to analyze recurrence-free survival (RFS) outcomes in different groups of patients.

Cell subgroup and immune-related analyses

After grouping the TCGA patients, CIBERSORT was used to observe immune-related cell scores and analyze the correlations among these scores, and ESTIMATE was used to calculate and observe the overall immune score and stromal score of tumor samples.

Prediction of the sensitivity of cell subgroups to drugs

Drug response prediction was performed using the “pRRophetic” package in R. Specifically, ridge regression was used to estimate the half-maximal inhibitory concentration (IC₅₀) for each patient, and the accuracy of prediction was determined by 10-fold cross-validation based on the Genomics of Drug Sensitivity in Cancer database.

Statistical analysis

Single-cell sequencing data were analyzed through the R package, single-cell analysis plots were generated via the “ggplot2” package, and Kaplan-Meier analysis was performed with the “glmnet” and “survival” packages. P<0.05 was considered to indicate statistical significance, and all the statistical analyses were performed using R language version 3.6.1.

Results

Sorting of ICC cell groups

Via the t-SNE algorithm, the heterogeneity of CSCs in ICC patients was sorted to observe. After dimensionality reduction and unsupervised cell clustering, we identified malignant (48.30%; EPCAM, KRT7, KRT19), biliary epithelia (30.18%; KRT7, KRT19), macrophages (8.74%; CD14), T (5.33%; CD2, CD3D, CD3E), endothelial (2.26%; ENG, VWF), dendritic (1.49%; CD1C), fibroblasts (2.13%; ACTA2, COL1A2), B (1.33%; CD79α), and hepatocytes (0.23%; APOC3, FABP1, APOA1) cells as nine distinct lineages based on marker gene expression (Figure 1A,1B). And EPCAM were highly expressed in malignant cells than the other cells (Figure 1C,1D). As EPCAM, an epithelial cell adhesion molecule and one of Wnt/beta-catenin signaling components, was identified as a marker for stem/progenitor cells, and serves as a prognostic marker, a therapeutic target, and an anchor molecule on circulating and disseminated tumor cells (CTCs/DTCs), which are considered the major source for metastatic cancer cells (21,22).

Figure 1 Sorting of intrahepatic cholangiocarcinoma cell groups. (A) Cancer stem cells in intrahepatic cholangiocarcinoma patients were grouped by the t-SNE algorithm according to surface receptor markers; (B) pie chart for proportions of each cell group; (C) distribution and expression levels of EPCAM, KRT7 and KRT19 antigen positive cell groups; (D) surface antigen expression for cell groups. t-SNE, t-distributed stochastic neighbor embedding.

Sorting of ICC stem cell groups

Subsequently, the CSCs in ICC were assessed using CD13, CD90, and CD133 (Figure 2). And Grassi et al. recently demonstrated that DLK1, a non-canonical NOTCH ligand, was discovered as a regulator of stem cell pools and tissue differentiation in liver (23). According to the t-SNE algorithm, the cholangiocarcinoma cells were divided into the following five groups: CD13-positive CSCs (14.74%), CD133-positive CSCs (18.72%), CD90-positive CSCs (0.41%), DLK1-positive cholangiocarcinoma cells (18.01%), and other cholangiocarcinoma cells (48.12%), respectively (Figure 2A,2C). Then, dot plot exhibiting the gene expression of different cell groups: ANPEP, SOX9, and ICAM1 were highly expressed in CD13-positive CSCs; PCAM, PROM1 and NANOG were highly expressed in CD133-positive CSCs; THY1 was most highly expressed in CD90-positive CSCs; and CD44, DLK1, SOX2 and APOC3 were most highly expressed in DLK1-positive cholangiocarcinoma cell groups (Figure 2B).

Figure 2 Sorting of intrahepatic cholangiocarcinoma stem cell groups. (A) Clustering of cancer stem cells in intrahepatic cholangiocarcinoma patients through the t-SNE algorithm based on surface receptor markers; (B) expression of tumor-associated genes in different cell groups; (C) pie chart showing the proportions of each cell group among all cells in tumor samples; (D) expression levels of different genes in different cell groups; (E) pseudotime analysis for paths of cell growth and development. CCA, cholangiocarcinoma; CSC, cancer stem cell; t-SNE, t-distributed stochastic neighbor embedding.

Further, heatmap of highly variable genes for the five major lineages presented that: FDCSP was highly expressed in CD13-positive CSCs; ITIH5, SPINK1, FXYD2, HSD17B11, and SERPINA1 were significantly upregulated in CD133-positive CSCs; SAA1 was highly expressed in DLK1-positive cholangiocarcinoma; and PSCA, LY6D, KRT6A, ALDH3A4, and AKR1C2 were significantly upregulated in cholangiocarcinoma cells vs. other cell types. These findings indicated that different expression levels of mRNAs in different cell groups predict their different developmental trajectories (Figure 2D). Moreover, the pseudotime trajectory axis derived from Monocle 2 showed the dynamic characteristics and heterogeneity of malignant epithelial cells: the cells in the Naive group were mainly CD13-positive and CD133-positive cells, while cells in the Recurrent group were mainly CD133-positive cells (Figure 2E).

Signaling pathway analysis in different ICC cell groups

To further investigate the differences in signaling pathway enrichment among various cell groups, we analyzed the expression patterns of ANPEP, PROM1, THY1, CEBPD, HSF1, and STAT1 using pseudotime analysis. This analysis revealed that the expression distribution of these genes corresponded well with the distribution of the respective cell groups (Figure 3A). Consequently, the differentially expressed genes were subjected to Hallmark gene set enrichment analysis to identify key signaling pathways. In summary, compared to cholangiocarcinoma cells, CSCs exhibited significant upregulation of EMT-related genes and downregulation of P53 pathway-related genes across the three subgroups. Additionally, genes associated with Notch signaling, Wnt/β-catenin signaling, and tumor necrosis factor-alpha (TNF-α) signaling were notably upregulated in the CD13, CD133, and CD90 subgroups, respectively (Figure 3B,3C).

Figure 3 Signaling pathway analysis for different intrahepatic cholangiocarcinoma cell groups. (A) Pseudotime analysis of the expression and distribution of ANPEP, PROM1, THY1, CEBPD, HSF1 and STAT1. Pseudo-time reflects the relative progression of cells along a trajectory and is dimensionless. (B) Hallmark analysis of differentially enriched signaling pathways in cancer stem cells; (C) expression of related signaling pathways in different cell groups according to t-SNE. CCA, cholangiocarcinoma; CSC, cancer stem cell; NES, Normalized Enrichment Score; GSEA, Gene Set Enrichment Analysis; AUC, area under the curve; t-SNE, t-distributed stochastic neighbor embedding.

Interactions between ICC cell groups and differences in transcription factors

Interactions between cell groups were analyzed using CellChat, which revealed the intercellular communication network among different cell clusters, indicating significant changes in interaction strength and cell type (Figure 4A). To clearly illustrate the differences in interactions, heatmaps were generated showing the relative expression of components of interaction networks, including genes in the EGF, HGF, Notch, WNT, and TGF-β pathways, across different ICC cell groups. Specifically, EGF signals may be transferred from cholangiocarcinoma cells to other subgroups, such as the CD13 or CD90 subgroups; HGF signals may be transferred from DLK1 subgroup cells to CD13 subgroup cells; Notch signals may be transferred from CD13 subgroup cells to other cells; and WNT signals may be transferred from CD90 subgroup cells to other cell groups (Figure 4B).

Figure 4 Interactions between intrahepatic cholangiocarcinoma cell groups and differences in transcription factors. (A) Cell-cell interaction; line segment thickness indicates the weight (upper) or number (lower) of interactions between groups, and loops indicate autocrine interactions. (B,C) The export (left) and input (right) signal patterns (B) and differences in transcription factors in different cell groups (C) shown in heatmaps. The bars represent relative signaling strength, which is dimensionless and derived from normalized interaction scores. (D) Violin plot showing the ex-pression of CEBPD, HSF1, and STAT1. (E) Representative IHC images for expression of CD13, CD133, CD90, HSF1 and STAT1 in ICC samples. For CK19, CD90, CD13 and CD133 staining, positive staining can be seen in the membrane of the tumor cells, while positive staining of HSF1 and STAT1 can be seen in the nuclear of the tumor cells. (Original magnification ×200; inset shows an enlarged area at ×400). CCA, cholangiocarcinoma; CSC, cancer stem cell; ICC, intrahepatic cholangiocarcinoma; IHC, immunohistochemistry.

Furthermore, the expression of transcription factors in different cell groups was analyzed using SCENIC, which indicated that CEBPD was highly expressed in the CD13 and CD133 subgroups, while HSF1 and STAT1 were highly expressed in the CD133 subgroup (Figure 4C). The expression distributions of these genes were also examined using t-SNE analysis, and the results were consistent with those of the differential ex-pression analysis (Figure 4D).

For further validation of the above results, we performed IHC to evaluate expression of cluster marker proteins and related transcription factors in 10 ICC samples. The out-comes showed that CD90 expression was restricted to a few cells within ICC, while CD13 and CD133 were partly expressed by ICC cells. In addition, approximately 50% of tumor epithelial cells were positive for HSF1 or STAT1, which overlapped the cells positive expressing CD13 or CD133 (Figure 4E).

Clinical prognostic ability of CSC subgroups

To assess the relationship between CSC subgroup enrichment and patient prognosis, we selected the top 20 differentially expressed genes with elevated expression in the three CSC populations (Figure 5A). First, ICC-related data were obtained from the TCGA database. Patients with cholangiocarcinoma in the TCGA cohort were then clustered and scored using unsupervised clustering with the ssGSEA tool. This process identified three CSC groups (clusters A, B, and C) (Figure 5B).

Figure 5 Clinical prognostic ability of cancer stem cell subgroups. (A) Upregulated differentially expressed genes in different cancer stem cell subgroups. (B) Clustering and scoring of intra-hepatic cholangiocarcinoma-related data in the TCGA database were performed with the ssGSEA tool. (C,D) Kaplan-Meier curves showing relapse-free survival in patients with high and low CD13/CD90/CD133 expression (C) and in different patient clusters (D). (E,F) Pattern recognition analysis for the correlations of different patient clusters with sex, pathological tissue, disease duration, and recurrence. CCA, cholangiocarcinoma; CSC, cancer stem cell; FC, fold change; CDF, cumulative distribution function; TCGA, the Cancer Genome Atlas.

Next, we evaluated the predictive value of CSC cluster membership and marker gene expression for RFS. The patients with high CSC133 expression had significantly shorter RFS compared to patients with low CSC133 expression [P=0.02, hazard ratio (HR) =0.7568, 95% confidence interval (CI): 0.3038–1.886]. Although the RFS was also shorter in the subgroup with high CD90 expression, this difference was not statistically significant (Figure 5C). Moreover, cluster A exhibited the shortest RFS among the highly enriched clusters (Figure 5D).

Additionally, pattern recognition analysis revealed the relationship between the expression levels of different genes and patient prognosis (Figure 5E). In cluster A, most patients were female, with tumors primarily located in the liver and a few around the hilum, and a high recurrence rate. For cluster B, none of the patients experienced recurrence; however, none of the clusters showed distinct characteristics in the disease course (Figure 5F). In summary, CSC heterogeneity is closely associated with ICC prognosis.

Correlation between CSC subgroup and immunity

To determine whether the CSC subgroup membership is associated with immune features, we assessed the levels of immune cells in clusters A, B, and C using CIBERSORT analysis. Notably, cluster A, which exhibited a high recurrence rate, had significantly reduced levels of activated B cells, effector memory CD4 T cells, and neutrophils. In contrast, cluster B, which had no recurrence events, showed significantly increased levels of neutrophils (Figure 6A).

Figure 6 Correlations between cancer stem cell subgroups and immune features. (A) CIBERSORT analysis method for observing immune cells in cluster A/B/C; (B) ESTIMATE algorithm for calculating immune and stromal scores in the 3 groups of patients; (C) observation of the expression of genes related to the immunotherapy response in the three groups; (D) CIBERSORT for calculating the correlation of cholangiocarcinoma immune cells with the CD133, CD13, and CD90 subgroups in the TCGA cohort. The bars represent the number of samples. *, P<0.05; **, P<0.01; ns, not significant. PRG, prognostic-related genes; CSC, cancer stem cell; MDSC, myeloid-derived suppressor cells; TCGA, the Cancer Genome Atlas.

We then calculated the immune score and stromal score for all patients using the ESTIMATE algorithm. The results indicated a significant correlation between the stromal score and subtype classification (Figure 6B). Further analysis of immunotherapy-related gene expression revealed high levels of CD276 and NRP1 in clusters A and C, while TNFRSF25 expression was markedly decreased in these clusters (Figure 6C).

Using CIBERSORT to calculate immune cell scores for ICC patients in the TCGA database, we found that CD133 expression negatively correlated with effector memory CD4 T-cell levels. Conversely, CD90 and CD13 expressions positively correlated with the levels of most cell types (Figure 6D). Overall, the heterogeneity of CSCs may be associated with variations in immune cell profiles, suggesting potential implications for immuno-therapy in ICC patients.

Drug sensitivity in CSC subgroups

To assess the sensitivity of different CSC subgroups to drugs, we predicted the sensitivity of different CSC subgroups to drugs via the pRRophetic tool. According to the prediction results, among the groups, the CD13-positive subgroup had significantly lower IC₅₀ values for etoposide (Figure 7A), GDC0941 (Figure 7B), and pazopanib (Figure 7C); the CD133-positive subgroup had significantly lower IC₅₀ values for OSI.906 (Figure 7D) and lapatinib (Figure 7E); and the CD90 subgroup had notably decreased sensitivity to midostaurin (Figure 7F), CMK (Figure 7G), AZD7762 (Figure 7H), and BX.795 (Figure 7I).

Figure 7 Drug sensitivity in cancer stem cell subgroups. The pRRophetic tool was applied to predict the IC₅₀ of different subgroups of CSCs for different drugs. (A-C) Sensitivity of the CD13-positive subgroup to etoposide, GDC0941, and pazopanib; (D,E) sensitivity of the CD133-positive subgroup to OSI.906 and lapatinib; (F-I) sensitivity of the CD90-positive subgroup to midostaurin, CMK, AZD7762, and BX.795. CSC, cancer stem cell; IC₅₀, half maximal inhibitory concentration.

Discussion

In this study, we employed scRNA-seq to comprehensively delineate the transcriptomic landscape of human ICCs. Five clusters (DLK⁺, CD13⁺, CD90⁺, CD133⁺, and other cholangiocarcinoma cells) were observed in ICC, which presented different signaling pathway activities, such as HSF1 and STAT1 were highly expressed in the CD133 cluster, and consistent with the results of IHC tests. And signaling, like Notch and Wnt/β-catenin, transferred between above subgroups. Further, subgroups favored varied immune response and drug sensitivity, with CD133⁺ patients experiencing significantly shortened RFS. These findings underscore the clinical potential of delineating ICC subgroups to predict prognosis and drug resistance, paving the way for personalized therapeutic strategies in ICC.

The initiation of cancer is stochastic, and cancer development and progression do not follow a fixed procedure. In other words, cancer development is dynamic and continuous, which causes tumor heterogeneity. Cancer cells have distinct molecular profiles that play a major role in disease progression and treatment failure. Tumor heterogeneity confers different proliferation and metastatic abilities to cancer cells, as well as different sensitivities to drugs (24). Notably, scRNA-seq, has emerged as an effective and accurate tool for directly investigating intratumorally heterogeneity in samples from patients with a variety of cancers, such as lung and prostate cancers, and has facilitated great advancements in cancer treatment and prognosis evaluation (25). However, there are few studies on CSC heterogeneity, especially CSC heterogeneity in ICC. Yet, scRNA-seq data on ICCs are now available.

By analyzing scRNA-seq data on ICCs, CSCs accounted for 33.87%, aligning with prior findings by Cardinale et al. (26). Further analysis indicated that CD13, CD90 and CD133 could divide CSCs into three subgroups. The hepatocyte marker genes of CD90⁺ subgroup cells are APOA1 and APOC3 (27), indicating their bidirectional differentiation potential. The expression levels of the proliferation-associated gene NANOG were significantly increased in the CD133⁺ subgroup (28), and the expression levels of the metastasis- and differentiation-associated gene SOX9 were notably increased in the CD13⁺ subgroup (29), indicating differences in the proliferation, differentiation, and metastasis abilities of these subgroup cells. Therefore, CSCs can be distinguished using CD13, CD90 and CD133, and we confirmed the heterogeneity of these tumor cells in clinical samples.

After further analysis, the EMT pathway was enriched in all three subgroups of CSCs relative to cholangiocarcinoma cells, while the P53 pathway was significantly downregulated, suggesting that CSCs have increased metastatic ability and decreased apoptosis (30). In addition, Notch signaling was significantly enriched in the CD13⁺ subgroup, Wnt/β-catenin signaling was enriched in the CD133⁺ subset, and TNF-α signaling was enriched in the CD90⁺ subset. Notch signaling is associated with cell differentiation and metastasis (31), Wnt/β-catenin signaling is correlated with cancer cell proliferation (32), and TNF-α signaling is involved in regulating immune regulation (33). According to the expression profile of transcription factors, CEBPD and STAT1 were most highly expressed in the CD13⁺ subgroup. Previous studies have shown that CEBPD can initiate the expression of proteins involved in cell differentiation (34), and STAT1 is involved in apoptosis (35). Moreover, HSF1 was most highly expressed in the CD133⁺ subgroup. HSF1 is a transcription factor that resists cellular stress and is utilized to drive cell proliferation in a variety of cancer cells (36). The analyses of transcription factor levels and signaling pathway enrichment suggested that the CD13⁺ subgroup may have greater metastatic and differentiation abilities, the CD133⁺ subgroup may have greater proliferative ability, and the CD90⁺ subgroup may include tumor cells that escape immunity. Interestingly, according to the pseudotime analysis, we observed the reverse differentiation of cholangiocarcinoma cells into CD13⁺, CD90⁺ and CD133⁺ cells in ICC. Choi et al. have also revealed this phenomenon and suggested that it may be caused by tumor cell evasion of the immune response or hypoxic environmental stress (37).

Based on pseudotime analysis, the impact of CSC subtypes on lymph node metastasis and recurrence varied. Consequently, we performed unsupervised clustering of online data using our defined hallmarks and found that CD133⁺ cell levels were strongly associated with RFS in patients. Yamashita et al. previously reported that CD133 expression was a poor prognostic factor in human lung adenocarcinoma, correlating with significantly shorter RFS (38). Similarly, Tachikawa et al. observed that positive CD133 expression was strongly linked to shortened RFS and high recurrence rates in colorectal cancer patients following chemoradiotherapy (39). Razmi et al. also have shown that in ICC patients, higher CD133 expression is associated with increased metastasis and recurrence rates, with patients displaying elevated CD133 levels experiencing shorter overall survival (OS) and RFS (40). Overall, CD133 is a prognostic risk marker across multiple cancers, supporting the credibility of our predictive model in this study. Furthermore, immunotherapy-related gene expression varied among the three clusters, suggesting that CSC heterogeneity may influence immunotherapy response. Using the pRRophetic tool, we predicted drug sensitivities across CSC subgroups (41,42), with results indicating that the three clusters also differed in their sensitivity to specific drugs.

This study holds considerable potential to advance the field by offering a foundational understanding of CSC-driven mechanisms in ICC. By precisely identifying the molecular characteristics and unique signaling pathways of specific CSC subtypes, we provide critical data for designing targeted therapies that address these CSC subpopulations, potentially leading to more effective and individualized treatments. Additionally, our approach in identifying prognostic biomarkers within CSC subpopulations enhances our ability to predict ICC outcomes and tailor treatments accordingly. Despite these advancements, several knowledge gaps remain. While our study has demonstrated the existence and roles of distinct CSC subpopulations in tumor progression and resistance, the mechanisms by which these CSC subtypes interact with the immune microenvironment and evade immune surveillance are still unclear. Additionally, our understanding of how CSC subtypes dynamically respond to treatments in vivo is limited. To address these gaps, future studies could incorporate multi-omics approaches, such as spatial transcriptomics and proteomics, to gain a comprehensive view of CSC interactions within the tumor microenvironment. Moreover, in vivo and in vitro models will be essential for elucidating the specific contributions of CSC subtypes to ICC progression, treatment response, and resistance. Looking ahead, we foresee significant advancements in CSC-targeted research within ICC over the next 5 years. As scRNA-seq technology and computational methods evolve, they are likely to provide even deeper insights into the molecular dynamics of CSCs, uncovering novel therapeutic targets and refining our understanding of the tumor microenvironment. Integrating scRNA-seq data with clinical datasets could also enable the development of more accurate prognostic models and improved prediction of treatment responses, facilitating personalized medicine for ICC patients. We also anticipate that novel CSC-specific therapeutic targets will emerge, paving the way for more precise and effective interventions against CSC-driven ICC progression.

This study has certain constraints that should be acknowledged. First, while scRNA-seq provides detailed snapshots of gene expression, it does not capture the temporal dynamics of CSC populations, which may change in response to treatment or disease progression. Furthermore, our study focused on specific CSC markers (CD13⁺, CD90⁺, CD133⁺), which may not encompass the full spectrum of CSC diversity within ICC. Although we validated the expression of these CSC subgroups in human ICC specimens through IHC, further in vivo and in vitro studies are needed to confirm the differential roles of CSC subtypes in tumor progression. Expanding the marker panel and conducting longitudinal studies would also strengthen these findings.

Conclusions

In summary, we identify CD13, CD90, and CD133 as key biomarkers for distinguishing CSC heterogeneity in ICC and suggest that the distribution of these subgroups correlates with differences in tumor progression, drug sensitivity, and RFS. We believe this study provides a timely contribution to ICC research; by advancing our understanding of CSC subpopulations and their roles in ICC progression, we are one step closer to developing therapies that can effectively target this aggressive, treatment-resistant cancer, ultimately improving patient outcomes.

Acknowledgments

We sincerely thank the public databases mentioned in this study and the patients participating in the study.

Footnote

Reporting Checklist: The authors have completed the REMARK reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1286/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1286/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1286/prf

Funding: This study was supported by grants from the Guangdong Basic and Applied Basic Research Foundation (No. 2021A1515111142), and the National Natural Science Foundation of China (NSFC) cultivating grant of the Third Affiliated Hospital of Sun Yat-sen University (No. 2024GZRPYQN06).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1286/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by the Institute Research Ethics Committee of the Third Affiliated Hospital of Sun Yat-sen University (No. [2023]02-167) and informed consent was taken from all individual participants.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Lee AJ, Chun YS. Intrahepatic cholangiocarcinoma: the AJCC/UICC 8th edition updates. Chin Clin Oncol 2018;7:52.
2. Kelley RK, Bridgewater J, Gores GJ, et al. Systemic therapies for intrahepatic cholangiocarcinoma. J Hepatol 2020;72:353-63. [Crossref] [PubMed]
3. Becht R, Wasilewicz MP. New Options for Systemic Therapies in Intrahepatic Cholangiocarcinoma (iCCA). Medicina (Kaunas) 2023;59:1174. [Crossref] [PubMed]
4. Moris D, Palta M, Kim C, et al. Advances in the treatment of intrahepatic cholangiocarcinoma: An overview of the current and future therapeutic landscape for clinicians. CA Cancer J Clin 2023;73:198-222. [Crossref] [PubMed]
5. Sirica AE, Gores GJ, Groopman JD, et al. Intrahepatic Cholangiocarcinoma: Continuing Challenges and Translational Advances. Hepatology 2019;69:1803-15. [Crossref] [PubMed]
6. Clara JA, Monge C, Yang Y, et al. Targeting signalling pathways and the immune microenvironment of cancer stem cells - a clinical update. Nat Rev Clin Oncol 2020;17:204-32. [Crossref] [PubMed]
7. Babaei G, Aziz SG, Jaghi NZZ. EMT, cancer stem cells and autophagy; The three main axes of metastasis. Biomed Pharmacother 2021;133:110909. [Crossref] [PubMed]
8. Ricci AD, Rizzo A, Brandi G. The DNA damage repair (DDR) pathway in biliary tract cancer (BTC): a new Pandora's box? ESMO Open 2020;5:e001042. [Crossref] [PubMed]
9. Ricci AD, Rizzo A, Brandi G. Immunotherapy in Biliary Tract Cancer: Worthy of a Second Look. Cancer Control 2020;27:1073274820948047. [Crossref] [PubMed]
10. Wang L, Jin Z, Master RP, et al. Breast Cancer Stem Cells: Signaling Pathways, Cellular Interactions, and Therapeutic Implications. Cancers (Basel) 2022;14:3287. [Crossref] [PubMed]
11. Li Y, Zhong X, Zhang Y, et al. Mesenchymal Stem Cells in Gastric Cancer: Vicious but Hopeful. Front Oncol 2021;11:617677. [Crossref] [PubMed]
12. Sahin TK, Rizzo A, Aksoy S, et al. Prognostic Significance of the Royal Marsden Hospital (RMH) Score in Patients with Cancer: A Systematic Review and Meta-Analysis. Cancers (Basel) 2024;16:1835. [Crossref] [PubMed]
13. Guven DC, Erul E, Kaygusuz Y, et al. Immune checkpoint inhibitor-related hearing loss: a systematic review and analysis of individual patient data. Support Care Cancer 2023;31:624. [Crossref] [PubMed]
14. Rizzo A, Santoni M, Mollica V, et al. Peripheral neuropathy and headache in cancer patients treated with immunotherapy and immuno-oncology combinations: the MOUSEION-02 study. Expert Opin Drug Metab Toxicol 2021;17:1455-66. [Crossref] [PubMed]
15. Liang L, Yu J, Li J, et al. Integration of scRNA-Seq and Bulk RNA-Seq to Analyse the Heterogeneity of Ovarian Cancer Immune Cells and Establish a Molecular Risk Model. Front Oncol 2021;11:711020. [Crossref] [PubMed]
16. Zheng L, Li L, Xie J, et al. Six Novel Biomarkers for Diagnosis and Prognosis of Esophageal squamous cell carcinoma: validated by scRNA-seq and qPCR. J Cancer 2021;12:899-911. [Crossref] [PubMed]
17. Zhang M, Yang H, Wan L, et al. Single-cell transcriptomic architecture and intercellular crosstalk of human intrahepatic cholangiocarcinoma. J Hepatol 2020;73:1118-30. [Crossref] [PubMed]
18. Zhong YJ, Luo XM, Liu F, et al. Integrative analyses of bulk and single-cell transcriptomics reveals the infiltration and crosstalk of cancer-associated fibroblasts as a novel predictor for prognosis and microenvironment remodeling in intrahepatic cholangiocarcinoma. J Transl Med 2024;22:422. [Crossref] [PubMed]
19. Su M, Qiao KY, Xie XL, et al. Development of a Prognostic Signature Based on Single-Cell RNA Sequencing Data of Immune Cells in Intrahepatic Cholangiocarcinoma. Front Genet 2020;11:615680. [Crossref] [PubMed]
20. Qiu X, Hill A, Packer J, et al. Single-cell mRNA quantification and differential analysis with Census. Nat Methods 2017;14:309-15. [Crossref] [PubMed]
21. Terris B, Cavard C, Perret C. EpCAM, a new marker for cancer stem cells in hepatocellular carcinoma. J Hepatol 2010;52:280-1. [Crossref] [PubMed]
22. Jeng KS, Chang CF, Sheen IS, et al. Cellular and Molecular Biology of Cancer Stem Cells of Hepatocellular Carcinoma. Int J Mol Sci 2023;24:1417. [Crossref] [PubMed]
23. Grassi ES, Pietras A. Emerging Roles of DLK1 in the Stem Cell Niche and Cancer Stemness. J Histochem Cytochem 2022;70:17-28. [Crossref] [PubMed]
24. Prasetyanti PR, Medema JP. Intra-tumor heterogeneity from a cancer stem cell perspective. Mol Cancer 2017;16:41. [Crossref] [PubMed]
25. Gao Y, Shang S, Guo S, et al. Lnc2Cancer 3.0: an updated resource for experimentally supported lncRNA/circRNA cancer associations and web tools based on RNA-seq and scRNA-seq data. Nucleic Acids Res 2021;49:D1251-8. [Crossref] [PubMed]
26. Cardinale V, Renzi A, Carpino G, et al. Profiles of cancer stem cell subpopulations in cholangiocarcinomas. Am J Pathol 2015;185:1724-39. [Crossref] [PubMed]
27. Du X, Lai S, Zhao W, et al. Single-cell RNA sequencing revealed the liver heterogeneity between egg-laying duck and ceased-laying duck. BMC Genomics 2022;23:857. [Crossref] [PubMed]
28. Vasefifar P, Motafakkerazad R, Maleki LA, et al. Nanog, as a key cancer stem cell marker in tumor progression. Gene 2022;827:146448. [Crossref] [PubMed]
29. Panda M, Tripathi SK, Biswal BK. SOX9: An emerging driving factor from cancer progression to drug resistance. Biochim Biophys Acta Rev Cancer 2021;1875:188517. [Crossref] [PubMed]
30. Lamouille S, Xu J, Derynck R. Molecular mechanisms of epithelial-mesenchymal transition. Nat Rev Mol Cell Biol 2014;15:178-96. [Crossref] [PubMed]
31. Misiorek JO, Przybyszewska-Podstawka A, Kałafut J, et al. Context Matters: NOTCH Signatures and Pathway in Cancer Progression and Metastasis. Cells 2021;10:94. [Crossref] [PubMed]
32. Zhang Y, Wang X. Targeting the Wnt/β-catenin signaling pathway in cancer. J Hematol Oncol 2020;13:165. [Crossref] [PubMed]
33. Cruceriu D, Baldasici O, Balacescu O, et al. The dual role of tumor necrosis factor-alpha (TNF-α) in breast cancer: molecular insights and therapeutic approaches. Cell Oncol (Dordr) 2020;43:1-18. [Crossref] [PubMed]
34. Sun X, Jefferson P, Zhou Q, et al. Dominant-Negative ATF5 Compromises Cancer Cell Survival by Targeting CEBPB and CEBPD. Mol Cancer Res 2020;18:216-28. [Crossref] [PubMed]
35. Wang W, Lopez McDonald MC, Kim C, et al. The complementary roles of STAT3 and STAT1 in cancer biology: insights into tumor pathogenesis and therapeutic strategies. Front Immunol 2023;14:1265818. [Crossref] [PubMed]
36. Chin Y, Gumilar KE, Li XG, et al. Targeting HSF1 for cancer treatment: mechanisms and inhibitor development. Theranostics 2023;13:2281-300. [Crossref] [PubMed]
37. Choi YJ, Ingram PN, Yang K, et al. Identifying an ovarian cancer cell hierarchy regulated by bone morphogenetic protein 2. Proc Natl Acad Sci U S A 2015;112:E6882-8. [Crossref] [PubMed]
38. Yamashita N, Oyama T, So T, et al. Association Between CD133 Expression and Prognosis in Human Lung Adenocarcinoma. Anticancer Res 2021;41:905-10. [Crossref] [PubMed]
39. Tachikawa Y, Kawai K, Ozaki K, et al. CD133(+)HIF-1α(-) Expression After Chemoradiotherapy Predicts Poor Prognosis in Rectal Cancer. Anticancer Res 2022;42:2033-43. [Crossref] [PubMed]
40. Razmi M, Ghods R, Vafaei S, et al. Clinical and prognostic significances of cancer stem cell markers in gastric cancer patients: a systematic review and meta-analysis. Cancer Cell Int 2021;21:139. [Crossref] [PubMed]
41. Chen J, Wu S, Peng Y, et al. Constructing a cancer stem cell related prognostic model for predicting immune landscape and drug sensitivity in colorectal cancer. Front Pharmacol 2023;14:1200017. [Crossref] [PubMed]
42. Lin G, Gao Z, Wu S, et al. scRNA-seq revealed high stemness epithelial malignant cell clusters and prognostic models of lung adenocarcinoma. Sci Rep 2024;14:3709. [Crossref] [PubMed]