In the present study, clinical datasets were utilized to evaluate the scRNA-seq results obtained from The National Center for Biotechnology Information Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/; GSE138709) (1). Samples in this dataset (GSE138709) had been collected from 5 patients with ICC, including 3 male patients and 2 female patients. The original databases did not include the age for all patients; therefore, the present study is unable to provide an overall age range or median age for the patient cohort.

All sequencing reads were aligned to the human genome (reference no. GRCh38) using Cell Ranger software (version 3.0.2; 10× Genomics). Subsequently, the refined feature barcode matrix was utilized as the basis for data analysis.

The Seurat package (version 5.0.2; http://satijalab.org/seurat/articles/install.html) was used for processing and analyzing the scRNA-seq data in R language (version 4.3.1, 2023-06-16). The individual libraries were transformed into Seurat objects. The data were filtered to focus on genes expressed across ≥10 cells and each cell exhibited gene expression levels within the range of 500–20,000 genes, with mitochondrial transcripts accounting for <15% of the total. The cutoffs were set after performing a quality control evaluation of each library. Subsequently, the data underwent logarithmic transformation for normalization and variable genes were identified through an in-depth mean and variance analysis. The genes were scaled and principal components were derived. To determine the most informative principal components, an elbow plot was used and the leading components were applied for dimensionality reduction using the Uniform Manifold Approximation and Projection (UMAP) dimensionality algorithm (19).

The main cell populations present in both the stromal tissue and immune compartments were identified by analyzing the average expression levels of particular gene markers in each cluster: i) Malignant, epithelial cell adhesion molecule; ii) B cells, B cell scaffold protein with ankyrin repeats 1 and CD79A; iii) CD4 cells, CD3D and CD4; iv) CD8 cells, CD3E and CD8A; v) macrophages, CD14 and apolipoprotein E; vi) myeloid cells, SP100A9 and lysozyme (LYZ); vii) dendritic cells (DCs), CD1C and C-type lectin domain containing 9A (CLEC9A); viii) natural killer (NK) cells, neural cell adhesion molecule 1 and NK cell granule protein 7; ix) hepatocytes, apolipoprotein C1 and fatty acid binding protein 1; and x) fibroblasts, decorin and actin α2, smooth muscle. Initially, the dataset was filtered to include only these specific cells. There were only 568 cholangiocytes in the original literature, resulting in cholangiocyte-specific markers showing limited representation in our results (GSE138709) (1). Following the initial filtering, analysis was performed to identify subpopulations of T cells and myeloid cells, as previously described (19).

For each patient, the subpopulations of T cells and myeloid cells were examined by isolating the corresponding subset based on patient identity and processing the raw data according to the aforementioned methodology. Patient data with <100 cells in either the myeloid or T-cell categories for the individual patient analysis were excluded (19). The marker identification function within the Seurat package was used to identify the signature genes associated with each cluster (19).

The visualizations, including UMAPs, heat maps, violin plots and box plots, displayed gene expression levels that were normalized using the Scran method (version 1.10.2; http://bioconductor.org/packages/release/bioc/vignettes/scran/inst/doc/scran.html) (20). The calculated expression values were derived from the raw expression data.

To analyze the myeloid and T-cell populations that exhibited clustering specific to individual patients, the Harmony algorithm (version 1.2.3; http://github.com/immunogenomics/harmony) was used to mitigate the influence of patient-related variations, also known as batch effects (19). A cross-verification of the analysis outcomes for the clusters and their associated markers with those from a single-library cluster was conducted to ensure that the Harmony algorithm focused solely on technical variations, without accounting for biological factors (21).

The inferCNV software (Bioconductor; version 3.10; http://bioconductor.riken.jp/packages/3.10/bioc/html/infercnv.html) was used to estimate the initial copy number variants (CNVs) (22). CNVs for each region were calculated based on cellular expression levels, using a cutoff value of 0.1 to refine the data. All cells, with the exception of epithelial cells, were categorized as healthy cells. Subsequently, a denoising procedure was used to generate refined CNV profiles. The diffusion pseudotime analysis was performed using Monocle3 software (Bioconductor; version 3.10; http://bioconductor.jp/packages/3.10/bioc/html/monocle.html) on an anndata object that had been processed with Harmony correction and the results were converted into a Seurat-compatible format (23).

The established clusters of immune cell subpopulations, specifically myeloid and T cells, were labeled based on the average expression levels of the subsequent genetic indicators: i) Macrophages, CD14 and apolipoprotein E; ii) myeloid cells, SP100A9 and LYZ; iii) DCs, CD1C and CLEC9A; iv) CD8 cells, CD8A and CD8B; v) CD4 cells, forkhead box (FOXP3, FOXP6 and CD4; and vi) proliferative cells, marker of proliferation Ki-67.

The fgsea package (version 1.26.0; http://bioconductor.org/fgsea.html) was employed to perform a comprehensive analysis of gene set enrichment based on the outcomes derived from analysis of differential gene expression patterns, following the approach as previously described (24). The msigdbr package was used to retrieve the hallmark gene sets from the extensive Molecular Signatures Database collection (version 7.5.1; http://cran.r-project.org/msigdbr) (19,25).

SCENIC pathway analysis was performed using the motifs database associated with RcisTarget (version 1.20.0; http://www.bioconductor.org/RcisTarget.html) in conjunction with the GRNboost algorithm (SCENIC; version 1.3.1) (26,27). The RcisTarget package was utilized to identify binding sequences for transcription factors within the provided gene list (26). The AUCell package (version 1.22.0; http://bioconductor.org/AUCell.html) was utilized to assess the activity levels of each regulon group across individual cells (26).

Survival of patients with ICC was assessed utilizing the Cancer Genome Atlas dataset through the Gene Expression Profiling Interactive Analysis (GEPIA2) database (http://gepia2.cancer-pku.cn/). The KM-plotter website (https://kmplot.com/analysis/) is the most sophisticated online survival analysis tool, performing all calculations in real time. It contains the gene expression information and clinical case data of 364 patients with liver cancer. The best cutoff value of gene expression was automatically set. The details of this dataset were mentioned by Menyhárt (28). For the Kaplan-Meier survival curves, the follow-up time of patients was ≤5 years, and the log-rank test was performed to observe the statistical difference.

The Metabolic gEne RApid Visualizer (MERAV) website (http://merav.wi.mit.edu/) is designed to analyze human gene expression across a large variety of arrays (4,454 arrays). It contains the gene expression of 7 normal liver tissues and 14 liver primary cancer tissues. The differential expression analysis of SPP1 was conducted directly online.

A total of 2 patients with ICC were thoroughly evaluated to confirm their eligibility (aged >18 years, had no contraindications for surgery or infectious diseases) based on established medical criteria and provided written informed consent prior to participation, ensuring the patients were fully aware of the study's purpose, procedures and potential risks. At last, a 52-year-old female patient and a 63-year-old male patient with ICC from The University of Hong Kong Shenzhen Hospital were enrolled. Detailed information of the patients is presented in Table SI.

Total RNA was isolated from paraffin-embedded tumor tissues through a dewaxing process, followed by the utilization of TransZol reagent (TransGen Biotech Co., Ltd.), according to the manufacturer's protocol. A total of 1 µg RNA was used for cDNA synthesis using the TransScript First-Strand cDNA Synthesis Kit (TransGen Biotech Co., Ltd.). The thermocycling procedure used for template amplification was set as deformation (94°C, 5 sec), annealing (60°C, 15 sec) and extension stage (72°C, 10 sec) for 40 cycles. RT-qPCR reactions were performed in triplicate with the TransStart Tip Green qPCR SuperMix (TransGen Biotech Co., Ltd.) containing SYBR Green dye and detected by the CFX96 real-time PCR system (Bio-Rad Laboratories, Inc.). Primers for PCR were as follows: Human SPP1 forward, 5′-CGAGGTGATAGTGTGGTTTATGG-3′ and reverse, 5′-GCACCATTCAACTCCTCGCTTTC-3′; human beta-actin forward, 5′-CACCATTGGCAATGAGCGGTTC-3′ and reverse, 5′-AGGTCTTTGCGGATGTCCACGT-3′. The relative quantification method was used for quantifying mRNA level of the target gene, beta-actin was the reference gene and used for normalization, and then the normalized values were compared to obtain the fold change. The comparative Cq value (2^−ΔΔCq method) was used to calculate relative quantitative data (29).

The tissues collected were fixed in 10% neutral formalin at room temperature for 24 h. Subsequently, the tissue underwent gradient alcohol dehydration, was cleared with xylene and was embedded in paraffin. Subsequently, 4 µm-thick slices were prepared for IHC staining. In briefly, paraffin sections of excised tumor tissue were dewaxed with xylene and processed using a gradient of ethanol concentrations. Tissue slices were washed with PBS before staining. Endogenous peroxidase activity was quenched with 3% H₂O₂ for 30 min. Subsequently, the tissues were blocked with 3% BSA in PBS for 30 min at room temperature and incubated overnight at 4°C with antibody against human SPP1 (1:100; cat. no. ER1802-16; Hangzhou HuaAn Biotechnology Co., Ltd.). Subsequently, slices were incubated with HRP-conjugated secondary antibody (1:1,000 dilution; cat. no. 7074S; Cell Signaling Technology, Inc.) for 1 h at room temperature after washing with PBS three times. The immune complex was visualized using DAB (Vector Laboratories, Inc.; Maravai LifeSciences) as a chromogen. The reaction was terminated by washing with distilled water, and the sections were counterstained with 0.5% hematoxylin staining solution (cat. no. C0107; Beyotime Institute of Biotechnology) for 5 min at room temperature and mounted with neutral resin. Images were captured using a BX53 microscope (Olympus Corp.).

Statistical analysis and figure generation were performed using R (version 4.3.1; http://cran.rstudio.com/bin/windows/base/old/4.3.1/). Statistical differences in data between paired groups were analyzed using Student's t-tests and the outcomes were presented as the mean ± standard deviation.

To investigate the TME characteristics of ICC, the data of tumor tissues and adjacent tissues collected from 5 patients with ICC were downloaded from the GEO database. The details of the patients with ICC are listed in Table SII. Prior to transcriptomic analysis, histopathological validation had been performed using hematoxylin-eosin staining to confirm the malignant features of tumor regions and normal architecture of adjacent tissues. IHC staining of cytokeratin 19, a cholangiocyte-specific marker, further verified the origin of the tumor cells (Fig. S1A and B). The scRNA-seq data encompassed 4 tumor tissue samples and 3 samples of adjacent normal tissue (1). Finally, after quality control, an scRNA-seq dataset consisting of a total of 20,412 cells from the database was produced.

Using graph-based clustering (30) to analyze the cells from the aforementioned patients with ICC, 10 main cell types were identified, including tumor cells, hepatocytes, B cells, CD4⁺ T cells, CD8⁺ T cells, macrophages, DCs, fibroblast, monocytes and NK cells (Fig. 1A and B) (31). All the cell types were characterized using known markers (Figs. 1C and S1C). A UMAP plot was used to assess the distribution of these cells in tumor tissues (Fig. 1D) and adjacent tissues (Fig. 1F). The distribution of different cell types in the tumor and adjacent tissues from ICC samples indicated that macrophages were significantly increased in tumor tissues, whereas T cells were markedly reduced (Fig. 1E, G and H). The numbers of different cell types in every patient are displayed in Tabl SIII and Table SIV.

To determine the influence of ICC tumor cells on T cells, the presence of infiltrating T lymphocytes (TILs) in ICC tumor tissues was compared with the surrounding non-tumor tissues. A total of 9,654 T cells were shown and 7 subsets of TILs were identified in ICC tissues (Fig. 2A). The expression levels of certain marker genes (Figs. 2B and S2A) and the percentage of each cell type subset in each patient with ICC were analyzed (Fig. 2C). The percentage of subsets of TILs in tumor tissues and the surrounding non-tumor tissues of each patient were assessed (Fig. 2D and E). A UMAP plot was used to measure the distribution of T-cell subsets in tumor tissues (Fig. 2F) and adjacent tissues (Fig. 2G). These data demonstrate that central memory T cells [CD4-interleukin 7 receptor (IL7R) T cells] were decreased, while exhausted T cells [CD4-C-X-C motif chemokine ligand 13 (CXCL13) T cells] and regulatory T cells (CD4-FOXP3 T cells) were increased in the tumor tissue of X23 compared with adjacent tissues (Fig. 2H). Furthermore, the ratio of regulatory T cells to central memory T cells was higher in the TME of ICC tumor tissues compared with adjacent tissues (Fig. 2I). The expression level of CD44 in T cells in tumor tissues was higher compared with that in adjacent tissues (Fig. S2B).

A volcano plot demonstrated that 875 genes were upregulated and 354 genes were downregulated in T cells in tumor tissues compared with adjacent tissues (Fig. 2J). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis demonstrated that differential gene expression mainly occurred in ‘response to interleukin-7, ‘signal transduction by p53’, ‘positive thymic T cell selection’ and ‘DNA damage response’ (Fig. 2K). The expression levels of the top 10 highest differentially expressed genes were measured and their association with the prognosis of patients with ICC was analyzed. It was shown that the basic leucine zipper ATF-like transcription factor gene exhibited increased expression levels in T cells in ICC tumor tissues compared with adjacent tissues and exerted a notable influence on the prognosis of patients with ICC (Fig. S2C and D). This finding was also reported by a previous study (32).

The diffusion pseudotime analysis provided a window for studying cellular dynamic processes, facilitating the prediction of tumor cell development and lineage diversification. The developmental trajectories of the subsets of T cells were analyzed (Figs. 2L and S2E). Pseudotime analysis showed that central memory T cells and regulatory T cells in adjacent tissues had similar gene clusters but had different gene clusters in tumor tissues (Fig. S2E). In addition, most regulatory T cells in tumor tissues were distributed at, or close to, the pseudotime branching point 1 (Fig. S2E), which indicated that these cells were undergoing cell programmed changes such as cell fate differentiation. Thus, both cell proliferation and reprogramming could potentially account for the increased ratio of regulatory T cells to central memory cells in tumor tissues.

To detect the transcriptomic features of ICC tumor cells, scRNA-seq data including tumor cells, hepatocytes and cholangiocytes were analyzed (Fig. 3A). The expression level of each cell group marker gene was assessed (Fig. 3B). CNVs, which comprise amplifications and deletions, could notably accelerate the adaptive evolution and development of cancer (33). The results of the present study demonstrated that ICC tumor cells contained several amplification and deletion variants in the entire chromosome. Specifically, the amplification variants in chromosomes 1, 7, 8, 11 and 21 and the deletion variants in chromosomes 3, 6, 16, 17 and 22 (Fig. 3C). KEGG pathway analysis results demonstrated that ICC tumor cells had different gene expression levels in the ‘PI3K-Akt signaling pathway’, ‘focal adhesion’ cascade, ‘ECM receptor-interaction’ cascade, ‘cell cycle pathway’ and ‘p53 signaling pathway’ (Fig. 3D). Furthermore, a volcano plot was used to investigate the differential gene expression of ICC tumor cells. Finally, 176 genes were identified as upregulated and 228 genes and as downregulated in ICC tumor cells compared with adjacent cells (Fig. 3E).

The SPP1 gene was among the most upregulated genes in ICC tumor cells (Fig. 3F). To further assess the association of SPP1 with ICC tumors, analysis of data from the GEPIA2 disease database was performed, which demonstrated a significant upregulation of SPP1 levels in ICC tumors compared with normal tissue and in stage IV of ICC (Fig. S3A and B). Furthermore, the sequencing results were further verified using biological experiments. SPP1 mRNA levels were significantly increased in tumor tissues compared with adjacent tissues from both patients with ICC (Fig. S3C). The SPP1 protein level was also demonstrated to be upregulated in tumors compared with adjacent tissues using IHC staining (Fig. S3D). The results demonstrate that SPP1 expression is elevated in ICC tumor tissues. Thus, the role of SPP1 in ICC immune microenvironment modulation was evaluated. SPP1 was first identified in osteosarcoma cells to serve a role of mediating osteoblast adhesion (34). The data from normal liver, primary hepatocellular carcinoma (HCC) and primary ICC microarrays showed a significant increase in SPP1 expression levels in liver carcinoma, particularly in ICC, compared with normal liver tissues (Fig. 3G). Furthermore, the expression level of SPP1 had a significant association with adverse prognosis in patients with ICC (Fig. 3H). In addition, S100A2 and stratifin (SFN) were also associated with the prognosis of patients with ICC (Fig. S3E-G) (35,36).

Tumor stage is an important clinical diagnostic index for tumor prognosis. To identify an association between the expression levels of the SPP1 gene and the ICC tumor stage, the expression level of SPP1 in tumor cells at different ICC stages was analyzed. ICC stages T1-2 were defined as early-stage and ICC stages T3-4 were defined as late-stage. Tumor cells from patient X25 were excluded from the datasets due to the absence of tumor tissue data in the public database. Therefore, tumor cells from patients X20, X23 and X24 were used as the early-stage group and tumor cells from patient X18 as late-stage (Fig. 4A). The volcano plot indicated that there were 69 upregulated genes and 101 downregulated genes in late-stage tissues compared with early-stage tissues (Fig. S4A).

The differentially expressed genes and their association with the prognosis of patients with ICC were further investigated. These results demonstrated that the expression level of SPP1 was higher in late-stage ICC compared with earlier stages (Fig. 4B). Furthermore, aldo-keto reductase family 1 member C2 and S100P exhibited increased expression levels, while IL-32 and vitronectin exhibited decreased expression levels in the late-stage ICC tumors compared with early-stage ICC tumors, and all of the aforementioned genes were significantly associated with the prognosis of patients with ICC (Fig. S4B-F). KEGG pathway analysis showed that the early- and late-stage groups exhibited different expression levels of genes associated with ‘cellular senescence’, the ‘p53 signaling pathway’ and ‘ferroptosis’ (Fig. S4G).

Furthermore, a UMAP plot was utilized to analyze the subset of CD4⁺ T cells present in the late and early stages of ICC tumor tissues (Fig. 4C and D). An increase in the ratio of regulatory T cells to central memory T cells in the tumor tissues of late-stage ICC was demonstrated (Fig. 4E). Pseudotime analysis showed distinct gene expression clusters of various T-cell types that were different between the early and late stages of ICC (Figs. 4F and S4H). Compared with early-stage tissues, both central memory T cells and regulatory T cells in late-stage tissues were located at the pseudotime branching points 2 and 3, which suggested that these two subtypes of CD4⁺ T cells underwent cell reprogramming and differentiation during ICC tumor progression.

To identify the interactions between tumor and immune cells, a list of 1,800 reported interactions were compiled from previously published literature. This list included ligand-receptor pairs from several families, such as cytokines, chemokines, TNFs and receptor tyrosine kinases, along with interactions between the ECM and integrins (37). A bubble plot was employed to predict the ligand-receptor interactions that contributed to the signaling between the tumor cells and immune cell clusters in ICC (Fig. 5A). The details of SPP1-CD44, macrophage migration inhibitory factor-CD74 and CD44, CXCL12-C-X-C motif chemokine receptor type 4 and CCL3-C-C motif chemokine receptor type 1 interactions in the microenvironment of ICC were compiled in a chord diagram (Fig. 5B-E).