Multiomic Analysis Reveals Comprehensive Tumor Heterogeneity and Distinct Immune Subtypes in Multifocal Intrahepatic Cholangiocarcinoma

Intrahepatic cholangiocarcinoma (ICC) is the second most common primary liver cancer, and its incidence is increasing (1–3). The majority of patients are reported to have multiple lesions or distant metastases when diagnosed, with limited treatment options and dismal prognoses (4–8). Gemcitabine plus cisplatin chemotherapy is currently recommended as the only first-line treatment for these patients, which delivers a median overall survival (OS) time of only 11.7 months and an incidence of severe treatment-related complications of 70% (9).

In recent years, advances in targeted therapy and immunotherapy have transformed the treatment approach for ICC (10). Multiple FGFR inhibitors have demonstrated improved survival rates for patients with ICC in phase I/II trials (11–15). The FIGHT 202 trial revealed an overall survival time of 21.1 months for patients with locally advanced or metastatic cholangiocarcinoma with FGFR2 fusions or rearrangements (16). On the basis of these promising results, the FDA approved Pemazyre (pemigatinib) as the first targeted drug for patients with cholangiocarcinoma in April 2020. Other potential druggable targets in ICC are isocitrate dehydrogenase (IDH)-1 and -2 alterations. The ClarIDHy trial showed that ivosidenib could significantly prolong the progression-free survival (PFS) time of patients with pretreated IDH1-mutant cholangiocarcinoma compared with that achieved with a placebo (17). In addition, several studies have revealed a robust and durable treatment response to anti-PD-1 immunotherapy in patients with advanced ICC and known microsatellite instability (MSI) or mismatch repair (MMR) deficiency (18, 19). However, the overall response rates to targeted drugs and anti-PD-1 therapies for patients with ICC were less than 20% and 10%, respectively (10). The high heterogeneity of ICC tumors may be one of the reasons for the unsatisfactory efficacy of the current therapies.

A recent study evaluated the intratumoral genomic heterogeneity of ICC based on whole-exome sequencing (WES) data (20); however, little is known about the tumor origin and intertumoral genomic heterogeneity in multifocal ICC. Tumor heterogeneity in patients with multiple lesions may be more complex than that in patients with solitary tumors (21). In addition, the transcriptional, epigenomic, and immune heterogeneity of multifocal ICC also remains unclear, leaving the mechanisms of intertumoral heterogeneity in response to treatment poorly understood. Therefore, understanding the intratumoral and intertumoral heterogeneity of multifocal ICC at the genomic, transcriptomic, epigenomic, and local tumor microenvironment levels is critical for guiding personalized treatment decision making.

To profile the heterogeneity of multifocal ICC at the abovementioned levels and to identify the level suitable for ICC treatment, we performed a multiomic analysis, including WES, bulk RNA sequencing (RNA-seq), single-cell RNA sequencing (scRNA-seq), Illumina 850K Human Methylation microarray, multiple immunofluorescence staining, and multispectral imaging, on patients with multifocal ICC. In addition, we suggest a clinical practical ICC immune classification system with prognostic and potential decision-making value.

This study was approved by the Institutional Ethics Review Board of the First Affiliated Hospital, Sun Yat-sen University, and was conducted according to the Declaration of Helsinki. Written informed consent was obtained from each patient in the study. Sixteen patients with multifocal ICC who underwent resection before any adjuvant therapy from January 2015 to January 2021 were included, including 15 with bulk sequencing analysis and one with scRNA-seq analysis. The clinicopathologic information of these 16 patients is listed in Supplementary Table S1. Two experienced pathologists confirmed the diagnosis through histological examination. Tumor tissues (2–6 tumors per patient, 61 tumors in total) and adjacent nontumor liver tissues of 15 patients were snap-frozen in liquid nitrogen within 15 minutes after resection for subsequent WES, RNA-seq, and Illumina 850K Human Methylation microarray. For 6 of these 15 patients, at least two regions were sampled from the largest tumor of each patient for evaluating intratumor heterogeneity. Besides, we obtained another five normal bile duct epithelial tissue samples for methylation microarray analysis as the tissue control. A detailed sample and experiment illustration is shown in Fig. 1A and B.

To evaluate the heterogeneity of different tumor nodules within a patient with ICC in response to immunotherapy, 6 patients with advanced ICC who received anti-PD-1 therapy from September 2019 to January 2021 were included. The inclusion and exclusion criteria were described in Supplementary Materials and Methods and Supplementary Fig. S1. RECIST 1.1 was used to evaluate the treatment response of all the tumors, in which complete response (CR), partial response (PR), stable disease (SD), or progressive disease (PD) were classified. In addition, to assess whether the developed immune classification was associated with patients’ prognosis, we included an ICC cohort (n = 111, validation cohort 1) with prognosis and RNA-seq data from our center (detailed clinical characteristics in Supplementary Table S2), and another cohort (n = 118, validation cohort 2) with prognosis and gene expression microarray data from the public dataset (GSE89749) released by Jusakul and colleagues (22). To evaluate the association between immune classification and immunotherapy response, we included 15 patients who were treated with at least two cycles of anti-PD-1 monotherapy and had pretreatment qualified tumor tissues for RNA-seq and analysis. Patients who achieved CR and PR were defined as responders, and those with SD and PD were defined as nonresponders.

Total DNA and RNA was extracted from the above snap-frozen tissues using the QIAGEN DNeasy Blood & Tissue Kit (Qiagen), TRizol and QIAQuick PCR Extraction Kit (Qiagen), respectively. DNA and RNA libraries were later constructed and loaded on NovaSeq sequencing platform (Illumina), and 150 bp paired-end reads were generated.

Raw sequencing reads were processed by Cutadapt (RRID: SCR_011841) to remove adapter sequences. FASTQC (RRID:SCR_014583) was then used to evaluate the quality of each sequencing read, and low-quality reads were discarded in our analysis. Qualified reads were mapped to the human reference genome version GRCh37 (hg19) using BWA (RRID:SCR_010910; ref. 23) by default parameters. On the basis of BWA alignment, somatic mutations were called by the MuTect2 (24), and were further filtered by the criteria listed in Supplementary Materials and Methods. The selected mutations were regarded as high confident somatic mutations and were then annotated with ANNOVAR (RRID:SCR_012821; ref. 25) and subjected to subsequent analyses. Next, copy-number alterations (CNA) were determined using ADTEx (RRID:SCR_012059; ref. 26). The potential driver genes in our analysis were selected from COSMIC cancer gene census (27) and other published works (22, 28, 29). The druggable genetic alterations include 136 molecular targets, which are extracted from the TARGET database (30). Phylogenetic trees were constructed based on the mutated sequences using Maximum Likelihood, Neighbor-Joining, and Maximum Parsimony algorithm in MEGA7. The inferred trees were manually recolored, with branch/trunk lengths proportional to the number of mutations. The mutated driver genes were also marked in phylogenetic tree trunk or branch. To analyze the clonal evolution in our sample, clonal analysis was performed by PyClone (RRID: SCR_016873). Then, phylogenetic trees of each clonal cluster were constructed using CITUP (31). Finally, population structures were visualized by MapScape (32).

Raw sequencing reads were processed with Cutadapt (RRID: SCR_011841) to remove adapter sequences. FASTQC (RRID: SCR_014583) was then applied to evaluate the sequencing quality of each read. Qualified reads were mapped to the reference human genome (hg19) by HISAT (RRID: SCR_015530; ref. 33) with default parameters. RSeQC (RRID: SCR_005275; ref. 34) was used to measure gene expression abundance as reads per kilobase per million mapped reads (RPKM). Differentially expressed genes between normal and tumor samples were identified by Student t test with a significance level of P < 0.01. A curated list of 66 immune markers that encompass cell-surface markers of different immune cell populations were collected from the genomic sequencing project of HCC by The Cancer Genome Atlas (TCGA; ref. 35). On the basis of the expression level of these collected immune markers, patient samples were clustered into high- and low-immune group using hclust (RRID: SCR_009154) function in R. The immune cell composition was analyzed by CIBERSORT (RRID: SCR_016955; ref. 36) based on gene expression data, and differences of immune cell infiltration were identified between high- and low-immune group using Wilcoxon rank-sum test.

Methylation data were processed by ChAMP package (RRID:SCR_012891; ref. 37). The somatic methylation changes in an individual was defined as probes with methylation level changes larger than 0.2 compared to those in adjacent nontumor tissues. We further defined probes with methylation level upregulated in tumor tissue as hypermethylation events. Similarly, downregulated probes in tumor tissue were regarded as hypomethylation events.

Then, the somatic methylation changes were mapped to the promoter region of each gene; an average methylation level was calculated for quantifying the promoter methylation. To reveal the correlation between promoter methylation and RNA expression, we calculated their spearman correlation among different tumors in each patient. Genes with negative correlations to their promoter methylations were selected to plot a heatmap.

To further infer the regulatory role of DNA methylation, the MIRA score (38) was calculated for each tumor using the annotated pathway from MSigDB C2 set. According to the cluster results from RNA-seq analysis, the pathway activity regulated by DNA methylation was compared between different clusters.

Finally, to test whether DNA methylation can alter neoantigen depletion, the methylation level in the corresponding gene promoter region was compared between all the neoantigens and their match controls using a chi-square test.

To evaluate the intertumor heterogeneity (ITH) of each patient, the ITH index proposed by Zhang and colleagues was estimated (39). For a patient with multiple tumors, the ITH index was calculated as the mean Jaccard distance between variant set of two different tumors (Supplementary Materials and Methods). As defined by Zhang and colleagues, an ITH index ranges from 0 to 1. An ITH index approaching 0 indicates low ITH. Conversely, ITH index approaching 1 suggests high ITH.

According to the median ITH index, we divided the patients with multifocal ICC into low-ITH group and high-ITH group. The mutation status and CNAs of the known driver genes were compared between these two groups. Also, the related immune pathways and immune escape status (immunoediting and HLA-LOH) were compared. To uncover the potential regulation mechanism in the two ITH groups, we further identified differentially activated pathways at both transcriptomic and epigenomic level using gene set variation analysis (GSVA) and MIRA, respectively. Specifically, the two-sample t test and Wilcoxon test was applied to infer statistical significance.

To infer the HLA genotype, the BAM file of normal tissue was inputted to POLYSOLVER (40) for analysis. The LOH events for class I HLA genes were predicted by LOHHLA (41) based on the genotyping result from POLYSOLVER. An LOH event for a given HLA allele was determined if the estimated copy number was less than 0.5 and the P value was smaller than 0.05.

Neoantigens were characterized from nonsynonymous missense mutations detected by WES. Mutated peptides of 9–11 amino acids in length were inputted to NetMHCpan 4.0 (RRID:SCR_018182; ref. 42) to predict their binding affinity to HLA alleles. When the binding affinity of a given candidate antigenic peptide was predicted lower than 500 nmol/L, we considered it as a neoantigen. If at least five RNA-seq reads were mapped to the mutation position, and at least three contained the mutated base, we considered the corresponding neoantigen was expressed in the tumor sample.

The immunoediting score for each sample was computed on the basis of the method developed by Rooney and colleagues (43). The background model of immunoediting was constructed using tumor samples from patients with liver cancer (LIHC) in TCGA project. If the immunoediting score of a given sample was calculated as smaller than 1, we defined it as an immunoediting event (Supplementary Materials and Methods).

The 4-μm-thick slides cut from the formalin-fixed paraffin embedded blocks were dewaxed, rehydrated through, and fixed in neutral buffered formalin (10% NBF). Slides were stained to enable the simultaneous visualization of six marker Abs anti-CK19 (Catalog # ab52625, RRID:AB_2281020), anti-CD20 (Catalog # ab9475, RRID:AB_307267), anti-CD8 (Catalog # ab93278, RRID:AB_10563532), anti-CD68 (Catalog # 76437, RRID:AB_2799882), anti-PD-1 (Catalog # 84651, RRID:AB_2800041), and anti-CD4 (Catalog # ab133616, RRID:AB_2750883), on the same slide using PANO 7-plex IHC Kit (Catalog # 0004100100; Panovue). After performing antigen retrieval and blocking, we sequentially incubated six primary antibodies for 30 to 60 minutes at 37°C. Then the incubation of HRP-conjugated secondary antibody and tyramide signal amplification (TSA) with Opal was followed. Six staining cycles were performed for the following antibodies/fluorescent dye combinations: anti-CK19/Opal-540, anti-CD20/Opal-650, anti-CD8/Opal-690, anti-CD68/Opal-560, anti-PD-1/Opal-620, and anti-CD4/Opal-520. Finally, all slides were stained with 4′-6′-diamidino-2-phenylindole (DAPI; Sigma-Aldrich) for 5 minutes and enclosed with Mounting Medium 0022001010 (Panovue).

The stained slides were scanned using the TissueFAXS platform (TissueGnostics) at 10× magnification, and the scans were combined to build a single stack image. Spectral libraries were established and single-stained slides were applied to extract the spectrum of autofluorescence of tissues and each fluorescein, respectively. The library was then used to unmix the multispectral images with the StrataQuest software (TissueGnostics).

Spatial analysis of multispectral imaging data was conducted on reconstructed whole slides. To analyze the spatial colocalization of different cell phenotypes, we estimated a wide range of metrics with the StrataQuest software (TissueGnostics), including the distance and number of mutual neighbors for each pair of cell phenotypes, as well as the distances of the interested cells to the nearest neighbors. Besides, we output the detailed coordinates of all the cells of interest from the StrataQuest software, and performed advanced spatial statistics analysis using the statistics environment R and the package spatstat to evaluate the spatial organization of each pair of cell phenotypes. The function G(r), which is the cumulative distribution of the distance from a typical random cell X to its nearest cell Y, is estimated by the statistics environment R. The argument r is the radius of the area in which G(r) is assessed. The function G(r) and its derived curves can be used to distinguish clustered or dispersed spatial patterns when compared with a random point process model. When the observed G(r) curve was deviated from the theoretical G(r) curve, the spatial organization of the cell phenotype pair may be segregated.

For scRNA-seq, tumor tissues were digested to single-cell suspensions and CD45⁺ cells were enriched. Single cells were processed with the GemCode Single Cell Platform using the Gemcode Gel Bead, Chip and Libraray Kits (10× Genomics) following the manufacturer's protocol with modifications. Libraries were constructed according to the standard 10× Genomics protocol (Single Cell 5′ Reagent Kits v5.2 User Guide) and performed the paired-end 2 × 150 bp (PE150) sequencing on Novaseq 6000 (llumina).

Raw data of scRNA-seq were processed by Cell Ranger (version 6.0.0, 10× Genomics) pipeline (44) to obtain the feature-barcode unique molecular identifier (UMI) matrices. The human genome GRCh38 was used as reference in our study. Quality control, expression normalization, dimensionality reduction, and clustering were performed using the Seurat package (45). For each cell cluster, we manually annotated the cell type based on the 10 most expressed genes in this cluster. The CellMarker database (46) was used to assist the annotation process.

To identify potential markers for immune classification, we first clustered the multifocal ICC samples into two groups based on 66 known immune cell–related genes using hierarchical clustering. Differential expression analysis was subsequently performed between these two immune clusters, and candidate genes were ranked according to their expression fold change between the high- and low-immune groups. The 15 most differentially expressed genes were selected for the next-step analysis. To allow comparing expression measurement from different experiments and further remove technical biases, we first performed rank normalization on our gene expression dataset. On the basis of the rank-normalized data, we evaluated the discriminative capability of different combinations of selected genes on two immune groups using the area under ROC curves (AUC). By preserving the largest AUC value, we can select the optimal combination of immune marker genes for classification. Finally, we selected the combination of CD8B and ICOS as potential markers to divide the high-immune and low-immune group.

Next, to find the optimal cutoff value for classification, we calculated the Youden index in the ROC curve, and determined the rank-normalized cutoff value of CD8B as 35.2531 and ICOS as 39.5525. Using the optimal thresholds of CD8B and ICOS, we next divided samples from another ICC cohort into high-immune and low-immune groups. Specifically, for a given sample, if the expression level of these two genes both exceeded the cutoff values, we classified such samples into high-immune group, otherwise into low-immune group. After classification, the prognosis of patients from these two groups were compared using the Kaplan–Meier estimator, and the significance was evaluated using log-rank test.

To extend our classification model to microarray platform, another classifier was constructed using the expression level of CD8B and ICOS. Firstly, to remove biases between RNA-seq and microarray data, a quantile normalization method (47) was performed. In the normalization step, the expression levels of the multifocal ICC samples were used as reference, and the expression levels in microarray platform was transformed into the same distribution. Logistic regression was then conducted on the multifocal ICC samples to distinguish the high-immune and low-immune group. Again, using Youden index, an optimal cutoff value of 0.98 was selected for final classification of immune status. Another validation cohort released by Jusakul and colleagues (22) was downloaded from Gene Expression Omnibus (GEO) database (GSE89749). The constructed immune classifier was applied in this microarray dataset to predict high-immune and low-immune tumors. Survival analysis was also performed to compare prognosis of the two predicted groups using the Kaplan–Meier estimator. The association between the immune classification and treatment response to immunotherapy was also evaluated (detailed in Supplementary Materials and Methods).

The sequencing data reported in this paper have been deposited in the Genome Sequence Archive (48) in National Genomics Data Center (49), Beijing Institute of Genomics (China National Center for Bioinformation), Chinese Academy of Sciences, under accession number HRA000984, that are publicly accessible at http://bigd.big.ac.cn/gsa-human. Other detailed methods are described in Supplementary Materials and Methods.

To investigate the evolutionary trajectory of multifocal ICC, WES was performed on all the included tumor samples and adjacent nontumor liver tissue samples. The anatomical location of each tumor sample is described in Supplementary Fig. S2. WES achieved a median depth of 262× for tumors (Supplementary Table S3), identifying 2,676 nonsynonymous and 1,145 synonymous somatic mutations (Supplementary Table S4; Supplementary Fig. S3). Approximately 28.4% of the somatic mutations detected in the tumor genome were expressed at the transcriptome level (Supplementary Fig. S4).

We first assessed the nonsilent somatic mutations of all tumors and found that all tumors within individual patients shared some common mutations (Fig. 1C; Supplementary Fig. S5). To further examine the evolutionary trajectory of these multifocal ICC tumors, we constructed phylogenetic trees for each patient based on all the somatic mutations (Fig. 1D; Supplementary Fig. S6). Although a dominant pattern of branched evolution was observed, all tumors had common trunks, indicating that all the tumors within an individual patient were of the same origin. The robustness of phylogenetic trees was validated by the other three methods, including maximum likelihood, neighbor-joining, and maximum parsimony (Supplementary Fig. S7). In addition, CNAs were shared across multiple tumors within a patient in all cases (Fig. 1E; Supplementary Fig. S8). Collectively, these results indicated that the multiple tumors in each of our included patients with ICC might be intrahepatic metastatic tumors of the same origin.

Next, to further elucidate the clonal evolution of these ICC tumors with intrahepatic metastasis, we identified 74 clonal genotypes for 15 patients based on variation in their allele fractions (Supplementary Fig. S9). Detailed results are in the Supplementary Results.

Then, we depicted intratumoral and intertumoral heterogeneity at the genomic level for all patients with multifocal ICC with intrahepatic metastasis. The ITH index was estimated (39), and the median ITH index was 26% (range, 7%–57%). According to the median ITH index, we divided the patients with multifocal ICC into a low-ITH group and a high-ITH group (Fig. 2). A higher percentage of truncal mutations was observed in the low-ITH group than in the high-ITH group (67.7% vs. 24.0%, P < 0.001, two-sample t test), whereas the tumor mutational burden (TMB) was not significantly different between the two ITH groups (P = 0.438, Wilcoxon test; Fig. 2A). For the high-ITH group, the median proportion of unique somatic single-nucleotide variants (SNV) and indels within a tumor was 59.6% and that among multiple tumors in a patient was 74.4% (Supplementary Fig. S10). For the low-ITH group, the median corresponding proportions were 28.7% and 30.0%, respectively (Supplementary Fig. S10). Then, we evaluated mutations of the known ICC driver genes in these two ITH groups (Fig. 2A; Supplementary Table S5; refs. 22, 27–29). Quantificationally, 66.7% of the driver genes were truncal mutations for multiple tumors within a patient in the high-ITH group, while 94.4% of them were in the low-ITH group (Fig. 2B; Supplementary Fig. S6). When comparing the difference between these two groups, we found that mutations associated with epigenetic modifiers (P = 0.004, Pearson chi-square test), mutations in TGFβ signaling genes (P = 0.002, Pearson chi-square test), and mutations in WNT signaling pathway genes (P = 0.032, continuity correction of chi-square test) were more common in the high-ITH group (Fig. 2A). In terms of druggable target alterations (i.e., IDH1/2), 71.4% (25/35) were truncal events (Fig. 2B; Supplementary Fig. S6), which is much higher than that of HCC (20%) observed in previous study (50). However, when comparing the two ITH groups, the proportion of truncal druggable target alterations was lower in the high-ITH group than in the low-ITH group (average, 57.1% vs. 92.3%; Fig. 2B), suggesting that patients with ICC with different ITH degrees may theoretically respond differently to the corresponding targeted drugs. We also examined the mutational signatures associated with potential environmental exposure and found divergence within a tumor and among multiple tumors in both the low-ITH group and in the high-ITH group (Supplementary Fig. S11). Subsequently, we assessed CNA across tumors within a patient and between ITH groups. For the high-ITH group, the median percentage of unique CNAs was 67.3% within a tumor and 66.4% among multiple tumors in a patient (Fig. 2A; Supplementary Fig. S12). For the low-ITH group, the median corresponding percentages were 21.7% and 20.0%, respectively (Fig. 2A; Supplementary Fig. S12). When comparing the reported common CNAs for ICC in the two ITH groups (51), tumors in the high-ITH group harbored a higher percentage of CNAs than tumors in the low-ITH group, for example, 9q34.2 (P = 0.003, Pearson chi-square test), 18q12.2 (P = 0.006, Pearson chi-square test), and 1p36.13 (P = 0.029, Pearson chi-square test; Fig. 2A).

At the transcriptomic level, the RNA expression data analysis corroborated the tumor heterogeneity identified at the DNA level (Supplementary Table S6; Supplementary Fig. S13). We analyzed the differentially expressed genes between all tumors and normal samples and found more upregulated genes than downregulated genes (1,816 vs. 338). Similar to the genomic analysis, we next evaluated the differentially expressed RNAs in the two ITH groups (Supplementary Fig. S14). For the high-ITH group, the median proportion of unique differentially expressed RNAs was 42.0% within a tumor and 48.2% among different tumors within a patient. For the low-ITH group, the corresponding proportions were 26.5% and 47.6%, respectively. Then, we compared the biological pathway differences between the two ITH groups by performing GSVA. We found that oncogenic pathways, including ERBB signaling (P = 0.020, Wilcoxon test), mTORC1 signaling (P = 0.036, Wilcoxon test), apoptosis (P = 0.038, two-sample t test), and PI3K/AKT/mTOR signaling pathways (P = 0.040, two-sample t test), were more active in the high-ITH group (Fig. 2C).

Similarly, at the epigenome level, intratumor and intertumor heterogeneity of the hypermethylated and hypomethylated loci in the promoter region was detected (Supplementary Fig. S15). For the high-ITH group, the median proportions of the unique hypermethylated and hypomethylated loci within a tumor were 27.3% and 42.1%, and those among multiple tumors of a patient were 80.4% and 87.8%, respectively (Supplementary Fig. S16). For the low-ITH group, the corresponding proportions of the unique hypermethylated and hypomethylated loci within a tumor were 32.1% and 27.0%, and those among multiple tumors of a patient were 39.6% and 46.1% (Supplementary Fig. S16). To further evaluate the difference in DNA methylation changes in biological pathways in the two ITH groups, we compared the degree of promoter hypomethylation of genes in the biological pathways by MIRA score. We found that MIRA scores in several oncogenic pathways (i.e., apoptosis, WNT/β-catenin signaling, cell adhesion pathways) and metabolic pathways (i.e., glycosaminoglycan biosynthesis) were higher in the high-ITH group (all P < 0.001, two-sample t test), suggesting that these pathways were more active in the high-ITH group (Fig. 2D).

Immunotherapy represents a promising option for ICC, and the tumor immune microenvironment has been reported to be associated with immunotherapy (52). To investigate the immune microenvironment heterogeneity of multifocal ICC, bulk transcriptome analysis, scRNA-seq analysis, and multispectral imaging were performed. First, we estimated the cellular composition of immune cell infiltration in sequenced samples with the CIBERSORT method and found that the proportion of immune cells was similar across different tumors within a patient in those in both the high-ITH group and the low-ITH group, except for patient 1149 (all P > 0.05, Friedman test; Fig. 3A; Supplementary Table S7). The proportion of immune cells was also not significantly different between the high- and low-ITH groups (P = 0.074, Friedman test; Supplementary Table S7). Deconvolution of immune cells showed that the abundant immune cells were CD8⁺ T cells, CD4⁺ T cells, macrophages, and B cells in our patients with ICC. Then, multiplex immunostaining of paraffin sections for CD4, CD8, CD68, CD20, PD-1, and CK19 was performed to characterize the density of these cell phenotypes. Multispectral analysis showed that the densities of CD8⁺ T cells (P = 0.060, one-way ANOVA), CD68⁺ macrophages (P = 0.062, one-way ANOVA), and CD4⁺ T cells (P = 0.181, one-way ANOVA), except for CD20⁺ B cells (P < 0.01, one-way ANOVA), were similar among multiple tumors within the same patient (Fig. 3B and C), which was consistent with the results of the transcriptomic analysis. To further confirm the findings from the bulk RNA-seq data, we performed scRNA-seq on five tumor samples from one patient, and a total of 16,917 cells passed the quality control stage (2,471–4,143 cells; median: 3,605 cells/sample), with 5,677 total mapped reads and 1,375 detected genes on average. We identified and visualized 13 clusters using Umap (Uniform Manifold Approximation and Projection for Dimension Reduction; Fig. 3D). ICC malignant cells were confirmed by the detection of CNVs as inferred by scRNA-seq data (Supplementary Fig. S17). The identified immune cells mainly consisted of CD4⁺ T cells (IL7R and CCR7), CD8⁺ T cells (CD8A), and myeloid cells (CD68 and CD163), including macrophages and B cells (MS4A1 and CD79A; Fig. 3D; Supplementary Figs. S18 and S19). All these cell subtypes were shared among the tumors within the patient (Fig. 3E). Moreover, the proportions of these major immune cell subtypes were similar across the intrahepatic multiple tumors in this patient (P > 0.05, Friedman test; Fig. 3F). Similarly, in the bulk RNA-seq data, we found that multiple tumors shared the same major immune cell subtypes and similar proportions in the immune fractions (P > 0.05, Friedman test; Supplementary Fig. S20), consistent with the scRNA-seq data.

We further compared the immune pathway enrichments between the two ITH groups and found no significant difference between them (all P > 0.05, two-sample t test; Supplementary Fig. S21). We subsequently investigated immune-evasion mechanisms in the sequenced tumor samples, including immunoediting, HLA-LOH, and mutations of genes in the antigen-presenting pathway. First, we did not find any mutations in the core genes (i.e., ß2m) of the antigen-presenting pathway in our samples, whereas the prevalence of immunoediting and HLA-LOH events were 3.3% and 6.6%, respectively. Moreover, we found that immunoediting and HLA-LOH events were similar across the different tumors within individual patients in the majority of the patients, implying less heterogeneity in multifocal ICC in terms of tumor immune evasion capabilities (Supplementary Fig. S22).

Given the relatively low heterogeneity of the immune microenvironment across multifocal ICC tumors, we next evaluated whether different tumors in individual patients with multifocal ICC exhibited similar responses to anti-PD-1 immunotherapy. Six patients with multifocal ICC who received anti-PD-1 immunotherapy were enrolled and closely monitored (Supplementary Fig. S1). Among them, 5 patients presented consistent responses of different tumor nodules of a patient to immunotherapy (Fig. 3G; Supplementary Fig. S23). More detailed clinical information is provided in Supplementary Table S8.

Although immunotherapy is a promising treatment for ICC, only approximately 10% of patients respond to immunotherapy and identifying patients who are likely to benefit from immunotherapy is important. According to the above findings that the immune microenvironment is relatively similar among multiple tumors within an individual patient in both ITH groups, we speculated that multiple tumors in a patient might have the same immune classification. Thus, developing an immunophenotypic classification to identify potential candidates is critical for guiding precise immunotherapy decision making. In light of this, we conducted unsupervised hierarchical clustering of gene expression of the reported immune markers for all tumor samples (35). All tumor samples were divided into two clusters, one high-immune cluster and one low-immune cluster (Fig. 4A). Notably, these two clusters were applicable to nearly all ICC lesions from the same patient and from multiple patients, suggesting that intertumoral heterogeneity did not change the classification. Multiplex immunofluorescence staining revealed that higher densities of CD8⁺ T cells (P < 0.001), CD4⁺ T cells (P = 0.031), and CD20⁺ B cells (P = 0.017) were infiltrated in tumor samples from patients in the high-immune cluster than in those from patients in the low-immune cluster (Fig. 4B). To further elucidate the potential interactions between tumor cells and different immune cell phenotypes in these two immune clusters, we used the topographic information derived from multispectral imaging and estimated the mutual neighbor distances between each cell pair (i.e., CK19⁺ and CD8⁺, CK19⁺ and CD8⁺PD1⁺, CK19⁺ and CD68⁺, CK19⁺ and CD68⁺PD1⁺, CK19⁺ and CD20⁺, CD19⁺ and CD20⁺PD1⁺, CK19⁺ and PD1⁺). The mutual neighbor distances between CK19⁺ tumor cells and immune cells (i.e., CD8⁺ T cells and CD20⁺ B cells) were shorter for tumors with more immune cell infiltrates than for tumors with fewer immune cell infiltrates (mean, 19.93 μm vs. 48.60 μm, P = 0.045; 20.13 μm vs. 50.74 μm, P = 0.038; two-sample t test; Fig. 4C). Then, we discerned differential spatial organization of CD8⁺ T cells and tumor cells in these two distinct immune clusters, accompanied by their representative multifluorescent images (Fig. 4D) and cumulative distributions of cross-type distances, G(r) (Fig. 4E). The density of CD8⁺ T cells surrounding the tumor cells within a radius of 25 μm was much higher in the high-immune cluster (mean, 160.1 vs. 55.81, P < 0.001, two-sample t test; Fig. 4F). The spatial organization of CD8⁺ T cells and tumor cells in the high-immune cluster represented a clustered pattern, whereas that in the low-immune cluster represented a dispersed pattern. In addition, we found that the tumors in the high-immune cluster had upregulated immune pathways, such as the T-cell receptor signaling pathway and cytokine–cytokine receptor interaction pathway (Supplementary Fig. S24). In addition, these tumors presented a higher expression of the IFN signature (P < 0.001, two-sample t test; Fig. 4G), which has been demonstrated to be predictive of anti-PD-1 immunotherapy response in multiple solid malignancies (53). Taken together, these results imply that the tumors in the high-immune cluster have higher immune cell infiltration, closer tumor–immune cell interactions, and upregulated immune pathway enrichment than those in the low-immune cluster.

Next, we wondered whether this immune classification was clinically valuable and could be applied to other cohorts. Thus, we first compared the expression fold change of immune markers between these two immune clusters and selected the 15 most differentially expressed genes as candidates. To allow for comparisons of gene expression measurements from different experiments and to remove technical noise, we performed rank normalization on the training dataset by replacing each observation with its fractional rank. After evaluating their discriminative capabilities on the two immune clusters using the area under the ROC curve (AUC), the combination of CD8B and ICOS was selected as the final marker. We found that the levels of CD8B and ICOS were significantly different between the high-immune cluster and low-immune cluster in the training cohort (P < 0.001 vs. P = 0.011, two-sample t test; Fig. 5A), suggesting their ability to distinguish these two immune subtypes. The optimal rank-normalized cutoff values of CD8B and ICOS were determined to be 35.25 and 39.56, respectively, using Youden index. For a given sample, if the rank-normalized expression level of these two genes both exceeded the cutoff value, we classified the sample into the high-immune cluster; otherwise, we classified it into the low-immune cluster. Then, we used CD8B and ICOS levels to classify ICCs into two immune clusters in another independent cohort (n = 111; Fig. 5B). With the optimal cutoff values, these two markers could successfully classify the tumor samples into one high-immune cluster and one low-immune cluster (Fig. 5C), which was consistent with our above clustering. Moreover, the levels of CD8B and ICOS were also significantly different between the two immune clusters in this validation cohort 1 (both P < 0.001, two-sample t test; Fig. 5A). Subsequently, to eliminate the batch effect and platform differences, we performed batch correction on our validation cohort 1 and performed immune cluster classification using these two markers. As expected, the cutoff values remained effective (Supplementary Fig. S25), further suggesting the robustness of our model. Survival analysis further revealed that the two subtypes clustered by these markers displayed distinct overall survival (OS; P = 0.016, log-rank test) and recurrence-free survival (RFS; P = 0.008, log-rank test) rates in our validation cohort 1 (Fig. 5D). Multivariable Cox regression analysis showed that immune classification was a significant independent factor of the RFS rate (HR = 0.24; 95% CI, 0.08–0.70; P = 0.009) and a marginally significant factor of the OS rate (HR = 0.16; 95% CI, 0.02–1.19; P = 0.074; Fig. 5E). To broaden the applicability of our immune classifier, a model dedicated to handle microarray data was constructed on the basis of the two previously selected markers. A quantile normalization method (47) was used to first normalize the expression value in microarray platform so that it can distribute the same as in RNA-seq. Then, a logistic regression model was trained using levels of CD8B and ICOS. Again, a cutoff value of 0.98 is selected to classify the high- and low-immune clusters. We tested the classifier and cutoff values in another microarray dataset (validation cohort 2), and found that this new model can distinguish the high-immune tumors from the low-immune tumors (Supplementary Fig. S26), indicating a potential application of our selected markers in the microarray platform. Further survival analysis suggested that the high-immune group had better prognosis than the low-immune group, although did not reach statistical significance (Supplementary Fig. S26). In addition, the immune classification was also tested to determine whether it was associated with immunotherapy responses in another cohort treated with anti-PD-1 monotherapy (n = 15). One patient was categorized into the high-immune cluster and achieved PR to immunotherapy, 14 patients were categorized into the low-immune cluster, and 13 of them did not respond to treatment (SD or PD; Fig. 5F). Collectively, these results suggest that this immune classification approach is reliable, and can be used for patient stratification. Moreover, the immune classification may be associated with ICC immunotherapy efficacy, but is required to be further confirmed in future larger cohorts.

Considering the distinct immune status of the two immune clusters, we then explored the role of genomic alteration and epigenetic modification in the immune profiles of these two distinct ICC clusters. First, we analyzed mutations in known immune-related genes in all the tumors within the two immune clusters (54, 55). We found that 88.2% (15/17) of immune-related gene mutations were shared among different tumors in the same patient (Supplementary Fig. S27). Then, we found that there was no significant difference in the TMB between these two clusters (P = 0.34, Wilcoxon test; Supplementary Fig. S28). Moreover, the tumor samples in our study were all microsatellite stable.

Next, we analyzed the correlation between the DNA methylation level and RNA expression in the tumor samples. The average promoter methylation level was quantified and compared with the corresponding mRNA expression level. The results showed that the expression level of mRNA was negatively correlated with promoter DNA methylation in 1,080 to 6,845 genes in each patient, implying that hypermethylation of the promoter region can regulate the downregulation of the corresponding genes (Fig. 6A; Supplementary Fig. S29). To further investigate the impact of DNA methylation changes on immune signaling, we subsequently evaluated the degree of promoter hypomethylation of genes in the immune pathways by the MIRA score for tumors in different immune clusters. The MIRA scores in several immune-active pathways (i.e., activation of T and B cells, etc.) were higher in tumors of the high-immune cluster (all P < 0.01; Fig. 6B). However, the MIRA scores in suppressive immune pathways, including TGFβ signaling and IL2–STAT5 signaling, were significantly lower in tumors of the high-immune cluster (both P < 0.01, two-sample t test; Fig. 6B). These findings were validated in the TCGA dataset (n = 36) and a public dataset (n = 138, GSE89749; ref. 22), which concordantly showed that the MIRA scores in the immune-active pathways (i.e., activation of T and B cells) in the high-immune cluster were higher (all P < 0.01, two-sample t test; Fig. 6C). The above results imply that the epigenetic modification of the expression of genes in the immune pathways might be involved in the regulation of immune activity in ICC and thus deserves further investigation.

Finally, we explored the relationship between the expression of neoantigens and the methylation level of neoantigens. For genes that carried neoantigenic mutations, genes that were not expressed harbored a 5.3-fold increase in promoter hypermethylation when compared with those genes that were expressed (P < 0.001, chi-square test; Fig. 6D). To determine whether the observed hypermethylation was neoantigen specific, the level of promoter hypermethylation between all neoantigens and the same genes without neoantigens was compared. Both the nonexpressed neoantigens (OR = 0.63; P < 0.001, chi-square test; Fig. 6E) and expressed neoantigens (OR = 0.64, P = 0.16, chi-square test; Fig. 6F) were shown to display no enriched hypermethylation in the promoter regions when compared with the same genes in nonmutated samples. These results indicate that immune pressure might select for promoter hypermethylation in reducing neoantigen expression in ICC.

In this study, we demonstrated high intertumoral heterogeneity of the genomic profile, RNA expression signature, and DNA methylation profile in patients with multifocal ICC in the high-ITH group. The immune microenvironment across multiple tumors within a patient was relatively less heterogeneous in both the high- and low-ITH groups, and consistent treatment responses of multiple tumors within a patient to anti-PD-1 immunotherapy were observed. In light of this, we developed and validated an immune classification that was predictive of patient prognosis and associated with immunotherapy efficacy for ICC patients. Finally, although genomic alterations were similar between different immune clusters, promoter DNA methylation might contribute to different immune profiles of two immune subtypes by regulating immune gene expression and neoantigen depletion.

According to the mutations and structural aberrations, the different lesions in our patients with multifocal ICC were all intrahepatic metastases. Although multiple tumors in an individual were of the same origin, a dominant branched tumor evolution was observed with evident genomic heterogeneity. On the basis of the median ITH index, we divided patients into high- and low-ITH groups. The heterogeneities among different tumors within the same patient in genomic profiles, RNA expression, and DNA epigenetic profiles were relatively high in the high-ITH group. Moreover, we observed that patients in the high-ITH group had a higher prevalence of driver gene mutations associated with epigenetic modifiers (i.e., ARID1A) and mutations in genes in the WNT and TGFβ pathways, suggesting a strong epigenetic component and a role for WNT signaling and TGFβ signaling in ICC carcinogenesis and tumor heterogeneity (56–58). Regarding druggable alterations, 71.4% were truncal events, which was significantly higher than that of sorafenib-targeted alterations in HCC in our previous research (50). This finding is supported by the promising results of the Fight 202 trial and the ClarIDHy trial (16, 17). Pemigatinib, an FGFR inhibitor, can deliver a disease control rate of 82% and a median OS time of 21.1 months for patients with locally advanced or metastatic cholangiocarcinoma who have FGFR2 fusions or rearrangements. However, the percentage of truncal druggable target alterations in patients in the high-ITH group was only approximately 50%, much lower than that in patients in the low-ITH group (92.3%). These findings suggest that targeted therapy might not be effective for all tumor lesions in nearly half of the patients with multifocal ICC in the high-ITH group, emphasizing the importance of multiple-tumor sampling of all foci after surgery when treating patients with multifocal HCC. Moreover, tumors in the high-ITH group were found to have upregulated oncogenic pathways, such as the ERBB signaling pathway. Consistently, the promoter hypomethylation level of genes in the oncogenic pathways was also higher in the high-ITH group. Previous studies have revealed that upregulation of these pathways is associated with tumor progression and metastasis (59, 60), implying that tumors in the high-ITH group might be more aggressive than those in the low-ITH group.

In this study, we found a relatively high degree of intertumoral heterogeneity within a patient at the genomic, transcriptomic, and epigenomic levels in the high-ITH group, whereas the immune microenvironment heterogeneity of multifocal ICC within a patient appears to be lower in both the low- and high-ITH groups. The similarity of the local immune microenvironment among tumors within a patient suggested that manipulating the ICC microenvironment might be a suitable strategy because these treatments might theoretically affect all tumors within an individual. However, the proportion of patients in the high-immune cluster was only approximately 13% in our study, which was consistent with the low treatment response rate to immunotherapy among patients with ICC (10). Therefore, developing an immunophenotypic classification to identify patients in the high-immune cluster is important for guiding immunotherapy decision making. All the ICC samples in our study could be clustered into two distinct immune subtypes, and multiple tumors within a patient were almost classified into the same immune subtype; only 1 patient had tumors of different subtypes. Thus, patients with ICC could be correctly classified into one of the immune subtypes according to our immune clustering method, providing evidence to determine whether the patients are potential candidates for immunotherapy. Our established immune classification was well validated in two cohorts with different platforms although different cutoffs were needed in RNA-seq and microarray datasets due to different features of these two platforms. Also, the two immune subtypes were still clearly identified before and after batch effect correction, implying the reliability and generalization of the immune classification. However, future studies with larger cohorts are still required to confirm this classification system.

For the nodules in the high-immune cluster, more infiltration of CD8⁺ T cells, macrophages, and B cells; closer tumor–immune cell interaction; and upregulated IFN signature were observed. All of these factors were reported to be predictive of the anti-PD-1 therapeutic response across multiple solid malignancies (53, 61, 62). For tumor nodules with high-immune infiltration levels and active immune pathways, we hypothesized that these tumors might be more likely to respond to immunotherapy than other nodules. In the cohort treated with anti-PD-1 immunotherapy, the immune classification was associated with the treatment response to immunotherapy, with patients in the high-immune subtype achieving immunotherapy response, whereas patients in the low-immune subtype had no response. However, this finding needs to be further confirmed in future studies with larger cohorts. In addition, for resectable tumors, emerging studies have shown that neoadjuvant immunotherapy might have superior clinical efficacy to adjuvant application (63), and it is important to identify patients who are most likely to respond to neoadjuvant immunotherapy. Whether the immune classification that partially reflects the tumor immunity in our study plays a role in predicting which patients with ICC will benefit from neoadjuvant immunotherapy remains unknown and is worth further investigation in future clinical trials with multifocal ICC patients in the different immune subtypes. In future trials evaluating neoadjuvant immunotherapy for ICC, pretreatment biopsy for in-depth scientific studies to explore drug mechanism-of-action and efficacy biomarkers should be considered. Finally, these two immune clusters were associated with different OS and RFS times for patients with ICC, indicating that our immune classification could improve the patients’ prognostic stratification for better management of ICC.

Considering the different immune microenvironments among these two immune clusters, we explored the potential factors related to different immune profiles of ICCs. Our study showed that mutations in immune-related genes were mostly similar among different nodules within individual patients despite there being very few different mutations across the nodules. Epigenetic modification plays an important role in the occurrence and development of ICC (61), but little is known about the role of DNA methylation in the immune profiles of ICC. Therefore, we subsequently evaluated the impact of promoter DNA methylation on different immune profiles of ICC. We found that the promoter hypomethylation level of genes in the immune-active pathways (i.e., 18-gene IFNγ signature) for tumors in the high-immune cluster was significantly higher, while the promoter hypomethylation level of genes in the immunosuppressive pathways was lower for tumors in the high-immune cluster. CCL5, one of the genes in the 18-gene IFNγ signature, is an important chemokine recruiting CD8⁺ T cells into the local tumor microenvironment (64). Promoter hypomethylation of the CCL5 gene could lead to upregulation of the expression level of the CCL5 chemokine, which might further attract CD8⁺ T cells to infiltrate the ICC microenvironment and thus contribute to the immune “high” state of ICC. This hypothesis is supported by recent research showing that increased expression of CCL5 by epigenetic treatment could increase T-cell infiltration and promote the memory and effector T-cell phenotypes (65). Our findings suggest that although the immune-related gene mutations are similar across different nodules, the DNA methylation modification of immune genes might contribute to different immune activities among multiple nodules. In addition, we revealed a 5.3-fold increase in promoter hypermethylation of genes carrying neoantigenic mutations that were not expressed compared with those in genes that were expressed, indicating that promoter hypermethylation of neoantigens is an important mechanism of neoantigen repression in ICC. This finding is consistent with observations in lung cancer (62). Overall, epigenetic modification appears to play a critical role in the regulation of the immunogenomic profile of ICC. Our data indicate that epigenetic reprogramming therapy might deserve exploration as an alternative therapy for ICC with low immune cell infiltration. Further studies are warranted to further investigate the molecular mechanisms of DNA methylation in regulating the immune microenvironment of ICC.

In conclusion, we comprehensively characterized the landscape of genomic, transcriptional, epigenomic, and immune heterogeneity of multifocal ICC. The characteristics of druggable target alterations and the immune similarity among different lesions within a patient might help explain the clinical responses to targeted drugs and immunotherapy in patients with multifocal ICC. Moreover, we proposed an ICC immunophenotypic classification approach with prognostic and potential decision-making value. Future studies with larger cohorts are required to confirm this classification system.

No disclosures were reported.

S. Chen: Conceptualization, data curation, formal analysis, supervision, funding acquisition, investigation, methodology, writing–original draft, project administration, writing–review and editing. Y. Xie: Data curation, formal analysis, investigation, visualization, methodology, writing–original draft. Y. Cai: Data curation, formal analysis, investigation, visualization, methodology, writing–original draft. H. Hu: Data curation, formal analysis, investigation, visualization, methodology. M. He: Conceptualization, supervision, investigation, visualization, writing–review and editing. L. Liu: Data curation, investigation. C. Liao: Investigation, visualization. Y. Wang: Investigation, visualization. J. Wang: Formal analysis, investigation, visualization. X. Ren: Investigation, visualization. Q. Zeng: Investigation. H. Peng: Resources, investigation. S. Shen: Resources, investigation. S. Li: Resources, investigation. D. Li: Resources, investigation. J. Lai: Resources. B. Peng: Resources. J. Ren: Investigation, methodology. M. Kuang: Conceptualization, resources, supervision, funding acquisition, project administration, writing–review and editing. S. Peng: Conceptualization, resources, supervision, writing–review and editing.

The work is supported by the National Science Fund for Distinguished Young Scholars (grant to M. Kuang, No. 81825013), the National Natural Science Foundation of China (grant to S. Chen, No. 81801703; grant to M. Kuang, No. 81771958), the Guangdong Natural Science Foundation (grant to S. Chen, Nos. 2018A030310328 and 2021A1515010450), and the Kelin Outstanding Young Scientist of the First Affiliated Hospital, Sun Yat-sen University (grant to S. Chen, No. R08030).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.