Molecular Subgroups of Intrahepatic Cholangiocarcinoma Discovered by Single-Cell RNA Sequencing–Assisted Multiomics Analysis

Intrahepatic cholangiocarcinoma (ICC) is a relatively rare but highly aggressive form of primary liver tumors, following hepatocellular carcinoma, and accounts for 5% to 10% of all primary liver malignancies (1). ICC is a heterogeneous group of malignancies that can be found at any position of the intrahepatic biliary tree, has distinct genetic and genomic features, and therefore causes differences in survival outcomes (1). Surgical resection is an effective approach for ICC. Nevertheless, ICC patients are often diagnosed at a late stage when resection is not possible and only palliative chemotherapy is applicable (2).

The inflammatory environment near the biliary tree may facilitate the development and progression of ICC. Specifically, inflammation in the tumor microenvironment (TME) is promoted by the secretion of various cytokines and chemokines, further aiding in tumor progression and distant metastasis formation. TNFα, IL1β, IL6, TGFβ, and other cytokines have been proven to take part in a multistep process of ICC transformation (3–6). Several studies using bulk tumor transcriptome data have developed different ICC subtype classification systems, which correspond to specific pathway alterations and correlate with different prognoses (7–12). Dysregulated pathways (e.g., erythroblastic leukemia viral oncogene homolog 2 or fibroblast growth factor signaling) were identified as therapeutic targets (9, 13). One study using a cohort of 78 patients characterized four subtypes based on the cellular composition of the TME, which related to distinct prognosis and immune escape mechanisms (14). Nevertheless, the above studies were all based on transcriptome data. One study, thus far, has applied proteomics to identify molecular subtypes in ICC (15). Most of the studies on the ICC TME are based on computational methods to infer the abundance and populations of cellular components in the TME. The detailed comprehensive molecular features of the ICC TME are, however, still not clear, especially at the single-cell level.

In this study, we performed a multiomics analysis on a large series of ICC patient cohorts (n = 110). Three molecular subtypes, “chronic inflammation,” “metabolism,” and “chromatin remodeling,” were identified through the proteomics data. The chronic inflammation subtype was associated with the poorest prognosis in ICC, which was further validated in an independent validation cohort (n = 41). A random forest algorithm identified the key featured gene mutations within each molecular subtype using whole-exome sequencing (WES). Single-cell RNA sequencing (scRNA-seq) was used to generate a detailed characterization of the TME for the “chronic inflammation” ICC subtype. Through this analysis, an APOE⁺C1QB⁺ tumor-associated macrophage (TAM) subtype was identified and shown to affect T-cell inflammation, thus reshaping the TME in ICC into an inflammatory state. In vitro and in vivo experiments further validated the potential of APOE⁺C1QB⁺ TAMs as targets of immunotherapy against ICC.

The sample collection of this study (primary tumor tissue and blood samples) was approved by the Ethics Committee of Nanfang Hospital and the affiliated hospital of Southwest Medical University. The samples were collected during 2014 to 2021. Tumor tissues were formalin-fixed paraffin-embedded, and the blood samples were stored in lipid nitrogen. The patient studies were conducted in accordance with International Ethical Guidelines for Biomedical Research Involving Human Subjects (CIOMS). All participants provided written informed consent to participate in the study. The detailed patient information is available in Supplementary Table S1 (training cohort, n = 110) and Supplementary Table S2 (validation cohort, n = 41). A total of 84.5% of patients in our ICC cohort were HBV infected, which is different from a previous study (16). The two cohorts were constructed from Nanfang Hospital and the affiliated hospital of Southwest Medical University, and the clinicopathologic features are shown in Supplementary Tables S1 and S2. Overall survival (OS) refers to the time that begins at diagnosis and up to the time of death.

All samples were assessed by pathologic examinations to evaluate the proportion of cancer cells and degree of fibrosis/necrosis. The presence or absence of histologic necrosis of each slide was identified according to established criteria (17) and recorded as negative (necrotic area ≤ 5%) or positive (necrotic area > 5%). If tumor tissues were estimated to be greater than 90% on the slides after a thorough pathologic review, the samples were used for proteomics analysis. No perihilar, distal, and combined hepato-cholangiocarcinoma were included in the study. Pan-cancer MSK clinical data and genomic data (patients, n = 7,315; ref. 18) were downloaded from cBioPortal (https://www.cbioportal.org/) for validation.

A total of formalin-fixed paraffin-embedded (FFPE) tumor samples from the ICC training cohort (n = 110) and the ICC validation cohort (n = 41; 0.5–1 mg) were dewaxed and rehydrated and then subjected to acidic hydrolysis with formic acid (FA; Thermo Fisher Scientific). Proteins were denatured with 6 M urea (Sigma–Aldrich) and 2 M thiourea (Sigma–Aldrich) with the assistance of pressure-cycling technology (PCT; Pressure BioSciences Inc; 30 seconds 45,000 psi and 10 seconds ambient pressure, 90 cycles). Proteins were then digested into peptides with trypsin (1:20; Hualishi) and Lys-C (1:80; Hualishi) with the assistance of PCT (50 seconds 20,000 psi and 10 seconds ambient pressure, 120 cycles). The detailed settings have been described previously (19, 20). Peptides were labeled with TMTpro 16 plex (Thermo Fisher Scientific; ref. 21). Each batch contained 15 experimental samples and one pooled sample in the TMT126 channel for normalization.

The fractions (60 per batch) were separated using offline high-pH, reversed-phase fractionation with a Thermo Dionex Ultimate 3000 RSLCnano System (Thermo Fisher Scientific) and then combined to a total of 30 fractions per batch. Subsequently, the fractionated samples were separated with a Thermo Dionex Ultimate 3000 RSLCnano System (Thermo Fisher Scientific) with a reversed C18 column (1.9 μm, 120 Å, 150 mm × 75 μm i.d.; National Center for Protein Sciences). The LC system was operated with buffer A (2% acetonitrile and 0.1% FA in HPLC water) and buffer B (98% acetonitrile and 0.1% FA in HPLC water). The peptides were eluted with a 60-minute LC gradient of 5% to 28% buffer B at a flow rate of 300 nL/minute. The samples were then analyzed with a Q Exactive HF mass spectrometer (Thermo Fisher Scientific) with positive ion mode using the data-dependent acquisition mode. The database search included all reviewed human entries from UniProt (downloaded on April 14, 2020, containing 20,365 proteins) using Proteome Discoverer (version 2.4, Thermo Fisher Scientific). The detailed parameters have been described previously without modification (22, 23). Among the 10,888 proteins, 6,311 proteins were detected in samples with less than 10% missing values and were subjected to further analysis. To evaluate the batch effect, the principal component analysis (PCA) was performed with the protein matrix using R.

The “ConsensusClusterPlus” package from R software was applied to the proteomic data to identify potential molecular subtypes (“chronic inflammation,” “metabolism,” and “chromatin remodeling”) of ICC with the default parameters (24). Differentially expressed protein (DEP) analysis was performed with the “Limma” package using R software to assess the proteomic data using the default parameters (25). An empirical Bayesian method was applied to estimate the fold change between every two molecular subtypes using moderated t tests. The adjusted P values for multiple tests were calculated using the Benjamin–Hochberg correction. Feature protein in each subtype was defined as the highly expressed proteins, with adjusted P < 0.01.

Single-sample gene set enrichment analysis (ssGSEA) was applied on proteomics from the training cohort (n = 110) to evaluate the enrichment scores (“chronic inflammation,” “metabolism,” and “chromatin remodeling”) of each sample using the gene sets downloaded from The Broad Institute (https://www.gsea-msigdb.org/gsea/msigdb/). The “GSVA” (gene set variation analysis) package using R software was used for this analysis with the default paramters (26). The hallmark gene sets were downloaded from The Broad Institute (https://www.gsea-msigdb.org/gsea/msigdb/). All the gene sets of immune cell populations [aDCs, B cells, CD8 T cells, cytotoxic cells, DCs, eosinophils, iDCs, macrophages, mast cells, neutrophils, NK CD56^bright cells, NK CD56^dim cells, NK cells, pDCs, T cells, T helper cells, Tcm (cental memory), Tem (effector memory), TFH cells (T follicular helper), γδ T cells, Th1 cells, Th17 cells, Th2 cells, and Tregs (regulatory T cells)] were obtained from Bindea G and colleagues (27). The score using for survival analysis was defined as high and low using the best-cutoff method. Gene ontology (GO) analysis was performed with the “clusterProfiler” package using R software (28).

Next-generation sequencing analyses were performed according to protocols reviewed and approved by the Ethics Committee of Nanfang Hospital. DNA extraction, library preparation, and target capture enrichment were performed using previously described protocols with minor modifications (29). Briefly, genomic DNA from FFPE tumor samples and matched paired blood samples from patients were isolated using the DNeasy Blood and Tissue Kit (Qiagen) following the manufacturer's protocol. After quantification of DNA with Qubit 3.0 using the dsDNA HS Assay Kit (Life Technologies), a NanoDrop 2000 (Thermo Fisher) was used to evaluate the DNA quality, and samples with A₂₆₀/A₂₈₀ ratio of 1.8–2.1 were used for downstream application.

The KAPA Hyper Prep kit (KAPA Biosystems) was used to prepare libraries according to previously described protocols (30). Briefly, 1 to 2 μg of genomic DNA was sheared by a Covaris M220 instrument into ∼350-bp fragments. End repair, A-tailing, and adaptor ligation of fragmented DNA were performed with the KAPA Hyper DNA Library Prep kit (Roche Diagnostics). Agencourt AMPure XP beads (Beckman Coulter) were used for size selection. The Adaptor-Ligated Library was amplified in a thermal cycler (SimpliAmp Thermal Cycler, Thermo Fisher) using the following program: 98°C for 2 minutes, then 6 cycles of 98°C for 30 seconds, 65°C for 30 seconds, 72°C for 1 minute, followed by one cycle of 72°C for 10 minutes. All primers used during library construction were from IDT (xGen Library Amplification Primer Mix, 96rxns). 500 ng of the amplified library DNA was added into a 1.5-mL LoBind tube. WES was performed with Agilent SureSelect Human All Exon V6 (Agilent Technologies) according to the manufacturer's instruction manual. PCR amplification was performed on captured libraries with KAPA HiFi HotStart ReadyMix (KAPA Biosystems). The KAPA Library Quantification kit (KAPA Biosystems) was used to quantify the purified library. A Bioanalyzer 2100 was used to calculate the fragment size distribution (∼350 bp).

The Illumina NovaSeq 6000 platform (Illumina Inc) was utilized for genomic DNA sequencing in Acornmed Biotechnology Co., Ltd to generate 150-bp paired-end reads. bcl2fastq (v2.19; Illumina) was used to demultiplex the sequencing data. Trimmomatic (31) was used to remove low-quality (quality<15) or N bases. The alignment of the data to the hg19 reference human genome was then performed by the Burrows–Wheeler Aligner (bwa-mem; ref. 32), followed by further processing using the Picard suite (available at https://broadinstitute.github.io/picard/) and the Genome Analysis Toolkit (GATK; ref. 33). GATK-HaplotypeCaller was used to call germline single-nucleotide polymorphisms (SNP), and the GATK-Mutect2 was used to call somatic single-nucleotide variants and indels. Here, Picard software was used to deduplicate the sequence. The mutant allele frequency cutoff was set as 0.5%. Common variants were removed according to dbSNP (https://www.ncbi.nlm.nih.gov/SNP/) and the 1000 Genomes Project (http://www.1000genomes.org). Gene fusions and copy-number variations were identified by FACTERA (34) and ADTEx (35), respectively. For the tumor tissue samples, the log₂ ratio cutoffs for the copy-number gain and the copy-number loss were defined as 2.0 and 0.6, respectively. Tumor mutational burden was defined as the number of somatic, coding, base substitution, and indel mutations per megabase of the examined genome and was calculated as previously described (36). Briefly, all base substitutions were considered, including nonsynonymous and synonymous alterations and indels in the coding region of targeted genes, except for known hotspot mutations in oncogenic driver genes and truncations in tumor suppressors (https://www.oncokb.org).

A random forest algorithm (37) was applied to the WES data to identify the most important mutations associated with the inflammatory response in ICC. Briefly, the in-house WES data set (n = 78) and inflammatory response score from the GSVA method were applied as input. The “ranger” package in R software was used to find the best hyperparameter in the regression process for the random forest model (38). The “e1071” package in R software was then used to generate the random forest model (39).

RPMI-1640 (Gibco; cat. no. 11875-093) with 1 mmol/L protease inhibitor (Solarbio; cat. no. P6730) was used to transport human ICC tumor tissues (n = 3) from the training cohort. Tissues were digested with a dissociation enzyme cocktail prepared by dissolving 2 mg/mL Dispase II (Sigma–Aldrich; cat. 42613-33-2), 1 mg/mL Type VIII collagenase (Sigma–Aldrich; cat. no. C2139), and 1 unit/mL DNase I (NEB, cat. no. M0303S) in phosphate-buffered saline (PBS) with 5% FBS (Gibco; cat. no. 16000-044) for 40 minutes at 37°C. The cells were dissociated by manual vortexing and pipetting and collected every 20 minutes and then filtered using a 40-μm nylon cell strainer (Falcon; cat. no. 352340). Red blood cell lysis buffer (Invitrogen) with 1 unit/mL DNase I was used to remove red blood cells. Finally, the cells were washed in PBS with 0.04% BSA (Sigma–Aldrich; cat. no. B2064). The concentrations of the single-cell suspensions were computed with Countess (Thermo) and adjusted to 1,000 cells/μL. Cells were loaded according to the Chromium single-cell 3′ kit standard protocol to capture 5,000 to 10,000 cells/chip position (V2 chemistry). Library construction, and all the other processes were performed according to the standard manufacturer's protocol.

Processing for RNA isolation and cDNA library preparation from the 3 tumor tissue samples were carried out according to the manufacturer's instructions. A limited dilution approach was applied to capture more than 10,000 cells. The surface of supersaturated beads was filled with oligonucleotide barcodes; thus, a bead could be paired with a cell in a microwell. Once cells were lysed in the cell lysis buffer, the beads were hybridized by polyadenylated RNA molecules. The RNA molecules were then used for reverse transcription to obtain cDNA. Each cDNA was tagged with a special cell label fragment during cDNA synthesis to carry information about the cell of origin. Libraries for scRNA-seq were generated using the BD Rhapsody platform and sequenced on the Illumina Novaseq 6000. STAR software was used to align FASTQs to the human reference genome (GRCh38) using the default parameters (40). For each sample, gene-barcode matrices were generated by counting the unique molecular identifiers and filtering noncell-associated barcodes. Finally, a gene–barcode matrix containing the barcoded cells and the gene-expression counts was generated.

Single-cell transcriptome analysis on the 3 tumor tissue samples was performed with Seurat package from R software with the default parameters (41). Nonlinear dimensional reduction was performed with the t-distributed stochastic neighbor embedding (t-SNE) method using the R software. The featured genes of each cluster were found by “Seurat” (highly expressed genes adjusted P < 0.001). The “GSVA” package from the R software was used to compute the score of each gene set as previously described above. NicheNET analysis was performed with the “nichenetr” package from R software and was used to analyze cell–cell interactions (CCI) by exploring the ligand and target gene pairs (42). Briefly, macrophages and T cells were selected from total cells and CCI analysis was performed using the nichenetr package. The gene-of-interest was selected as the feature gene of CD4⁺ T cells. The regulatory potential score was calculated using network propagation methods (42). Therefore, a total of 19 regulatory pairs were identified with the previously defined network.

Human ICC tumor tissues from the training cohort were fixed with 4% paraformaldehyde (Solarbio, P1110) before embedding them in paraffin (Solarbio, YA0012). Tissue sections (4 μm) were later deparaffinized with xylene (Aladdin, X112054) and rehydrated with graded alcohol (Aladdin, E111991). Antigen epitope retrieval was induced by microwave heating. Tissue sections were then blocked for 1 hour at room temperature in 100 μL per slide goat serum blocking solution (Proteintech, B900780) and PBS (Solarbio, P1020). The sections were incubated with a primary antibody against CD68 (1:200; 76437S, CST) overnight at 4°C. The next day, ICC tissue samples were stained using an anti-mouse/rabbit universal IHC detection kit (pk10006, Proteintech) following the manufacturer's instructions. Mounted sections were examined by light microscopy (Leica), images were analyzed with Image-Pro Plus (version 6.0), and high/low expression was used for staining thresholds.

Resected tumor tissues from ICC patients were fixed with 4% paraformaldehyde before embedding in paraffin. First, the prepared tissue sections (4 μm) were baked in an oven at 65°C for 1 hour, dewaxed with xylene, and rehydrated with graded alcohol. After rehydration, the sections were fixed in 10% neutral buffered formalin (Solarbio, G2161) at room temperature for 20 minutes, placed in an appropriate AR buffer (AKOYA Biosciences, NEL820001KT), and then placed in a microwave for 1 minute at 100% power and then for an additional 15 minutes at 20% power. The slides were blocked with Blocking/Ab Diluent (AKOYA Biosciences, NEL820001KT) after the slides cooled at room temperature and incubated with primary antibody at room temperature for 10 minutes. TBST (Solarbio, T1082) was used to wash the slides three times, and then the slides were incubated with Opal Polymer HRP Ms+Rb (AKOYA Biosciences, NEL820001KT) to cover tissue sections at room temperature for 10 minutes, generally 100 to 300 μL per slide. The slides were rinsed with TBST three times before incubation with Opal Signal Generation (AKOYA Biosciences, NEL820001KT). The microwave treatment, blocking, primary antibody incubation, introduction of Opal Polymer HRP, and signal amplification were repeated. Primary antibodies in this assay including APOE (1:100, 13366S, CST), CD68 (1:100, #76437, CST), and CD4 (1:100, #48274, CST). Once three targets were labeled, DAPI working solution (AKOYA Biosciences, NEL820001KT) was applied in the dark for 5 minutes at room temperature, and the slides were washed with distilled water and TBST before mounting. Finally, a confocal microscope (Nikon) was used to take images of those tissue samples, and the acquired images were analyzed using ImageJ (version 4.0), positive/negative expression was used for staining thresholds.

The following cell lines were purchased from the ATCC recently and maintained in a humidified atmosphere in a 5% CO₂ incubator (Thermo, 3111) at 37°C. HCCC9810, THP1, and HuCCT1 cells were maintained with RPMI-1640 medium (11875093, Gibco), whereas HIBEC, CCLP-1, LICCF, and KMCH1 cells were cultured with Dulbecco's modified Eagle's medium (DMEM; 11965092, Gibco). All media for cell culture were supplemented with both FBS (10%; 10099141, Gibco) and penicillin/streptomycin (1%; 10378016, Gibco). All cell lines used in this study were tested for mycoplasma every two months.

Cells were lysed with RIPA buffer (R0020, Solarbio) supplemented with a proteasome inhibitor (P6731, Solarbio) on ice. Protein concentrations were measured using the BCA protein assay kit (P0012S, Beyotime). Protein samples were separated by sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) and were transferred to nitrocellulose membranes (Merck Millipore, HATF00010). The membranes were then blocked using 4% fat-free milk (LP0031B, Solarbio) in PBS (Solarbio, P1020) for 1 hour at room temperature. The membranes were incubated with the primary antibody (KMT2D, 1:1,000, Proteintech; GAPDH, 1:5,000, Proteintech) overnight at 4°C. PBST (P1031, Solarbio) was used for washing membranes three times. The membranes were then incubated with the secondary antibody (Goat Anti-Rabbit IgG H&L, ab205718, Abcam) for 1 hour at room temperature the following day. The membranes were washed three times with PBST before scanning on a Bio-Rad gel imager according to the manufacturer's instructions. Acquired data were quantitated using Image Lab software, and high/low expression was used as staining thresholds. Three biologically repeated experiments (n = 3) were carried out.

Plasmids (pT3-EF1a-NICD, #46047; PT3-myr-AKT-HA, #31789; and pCMV-CAT-T7-SB100, #34879) were purchased from Addgene. Plasmid (1 μL) was cloned into E. coli DH5α for 30 minutes on ice. Heat shock was performed with water for 1 minute at 42°C. The E. coli were then cultured at 37°C with shaking at 250 RPM for 1 hour and then incubated on ice again for another 5 minutes. A single bacteria colony was obtained by growing the E. coli on the agarose. LB media (200 mL) were used to culture the bacterial colonies in an incubator overnight with shaking at 250 RPM. The medium was centrifugated at 3,000 × g for 5 minutes to obtain the precipitate. 10 mL Buffer S1 (AP-MX-EP-25, Axyprep), 10 mL Buffer S2 (AP-MX-EP-25, Axyprep), and 10 mL Buffer S3K were added sequentially to resuspend the precipitate. The mixed buffer was transferred into the preparation tube followed by centrifugation at 6,000 × g for 5 minutes. Lastly, 500 μL Eluent A was used to dissolve the plasmid.

The animal experiment was approved by the First Affiliated Hospital, Zhejiang University School of Medicine. C57BL/6J mice were purchased from the Model Animal Research Center of Nanjing University (China). All mice were housed in the specific pathogen-free facility of the First Affiliated Hospital, Zhejiang University School of Medicine, with approval from the Institutional Animal Care and Use Committee. To establish a spontaneous ICC mouse model, a plasmid mixture carrying 25 μg pT3-EF1a-NICD, 25 μg pT3-myr-AKT-HA, and 25 μg pCMV-CAT-T7-SB100 per mice was injected at high pressure into the tail vein of the mice, injected with 2 mL plasmid mixture in 7–8 seconds. One week later, the mice were treated with 10 μg/g CSF1R antibody (BE0213, Bio X Cell) via intraperitoneal injection, and the same injections were performed every 3 days thereafter. On day 28, tumors were collected and further analyzed.

The spontaneous ICC tumor tissue sections from mice (4 μm) were deparaffinized with xylene and rehydrated with graded alcohol. Antigen epitope retrieval was induced by microwave heating before blocking in blocking solution (Proteintech, B900780) in PBS for 1 hour at room temperature. Ki67 (1:200; 550609, BD Biosciences) primary antibodies were incubated overnight at 4°C. The following day, tissue sections were washed three times with PBS before staining using an anti-mouse/rabbit universal IHC detection kit (pk10006, Proteintech) according to the manufacturer's instructions. Mounted sections were examined by light microscopy (Leica); images were analyzed with Image-Pro Plus software (version 6.0), and positive/negative expression was used for staining thresholds.

Mouse spontaneous ICC tissue samples were deparaffinized with xylene and rehydrated with graded alcohol as stated above. Antigen epitope retrieval was induced by microwave heating. Tissue sections were permeabilized with 0.2% Triton X-100 (Aladdin, T109027), blocked with 5% BSA for 1 hour at room temperature, and then incubated with the primary antibody (APOE, 1:100, #49285, CST; C1QB, 1:100, abs136795, Absin; F4/80, 1:100, sc-365340, Santa Cruz) overnight at 4°C. The following day, tissue sections were washed with PBS three times. Samples were then stained with fluorescent secondary antibodies (Alexa Fluor 488, 1:400, A-11008, Invitrogen; Alexa Fluor 594, 1:400, A-11032, Invitrogen) in the dark for 1 hour at room temperature and washed with PBS three times before mounting onto slides using Mowoil (81381, Sigma) supplemented with DAPI (Beyotime, C1002) to stain the nucleus. Finally, a confocal microscope (Nikon) was used to take images of ICCA tissue samples. Acquired images were analyzed using ImageJ (version 4.0), and high/low expression was used as staining thresholds.

To induce THP1 cells to differentiate into M0 macrophages, 1×10⁵ THP-1 cells were treated with 80 nmol/L phorbol 12-myristate 13-acetate (PMA; ab120297, Abcam) in 12-well plates. After 48 hours, 2×10⁵ ICCA tumor cells HCCC9810 or CCLP1 were cocultured with M0 macrophages on glass coverslips in 12-well plates. Immunofluorescence was performed. First, the cells were washed with PBS and fixed with 4% paraformaldehyde before permeabilization with 0.2% Triton X-100 at room temperature. Samples were next blocked with 2% BSA for 1 hour at room temperature and then incubated with the primary antibodies (APOE, 1:100, 13366S, CST; CD68, 1:100, #76437, CST; C1QB, 1:100, abs136795, Absin) for 1 hour at room temperature. To examine the colocalization of proteins of interest, cells were stained with Alexa Fluor 488- and 594-labeled fluorescent secondary antibodies (Invitrogen) for 1 hour in the dark at room temperature. Coverslips were washed three times with PBS before mounting onto slides using Mowoil supplemented with DAPI to stain the nucleus. Finally, the cells were examined with a confocal microscope (Nikon). The acquired images were analyzed using ImageJ (version 4.0), and high/low expression was used as staining thresholds.

Mouse spontaneous ICC tumor tissue sections (4 μm) were deparaffinized with xylene and then rehydrated with graded alcohol. Tissue sections were washed three times using PBS before nuclear staining with hematoxylin (Solarbio, G1080) for 30 minutes at room temperature before being washed with PBS three times. The sections were then exposed to 1% ammonia water (Solarbio, G1822) to change the hematoxylin-stained nuclei from a reddish to blue–purple appearance. Subsequently, 75% alcohol was used to rinse the tissue sections for 2 minutes at room temperature, and then the cytoplasm was stained with eosin (Solarbio, G1100) for 1 hour at room temperature. The sections were then directly rinsed with graded alcohol. Finally, xylene was used to displace the anhydrous alcohol before mounting with slides. Sections were examined by light microscopy (Leica), images were analyzed with Image-Pro Plus software (version 6.0), and normal/abnormal structures of the liver were used as staining thresholds.

THP1 cells were treated with 80 nmol/L PMA for 48 hours to induce THP1 to M0 macrophages. Differentiated M0 macrophages were then cocultured with ICC tumor cells, and the cell number was the same as in the previous section. Before the cells were incubated with fluorochrome-labeled antibodies, they were stimulated with 1× Cell Stimulation Cocktail (plus protein transport inhibitors; eBioscience, 00-4975-03) at 37°C for 6 hours. Human peripheral blood mononuclear cells (PBMC) from patients were first isolated by Ficoll density gradient (10771, Sigma). Briefly, isotonic Ficoll solution was added into the blood slowly and centrifuged at 600 × g for 30 minutes. Buffy layers containing cells were collected and washed for PBMC recovery. Human T cells were isolated from PBMCs through negative selection to avoid T-cell activation, with the Pan T-cell isolation kit (130096535, Miltenyi Biotec), LS column (130042401, Miltenyi Biotec), and MidiMACS separator (130042302, Miltenyi Biotec) following the manufacture's protocol, and the purity is more than 90%. Directly before analysis, T cells were first incubated with a Golgi inhibitor cocktail (2 μL of the cocktail was added for every 1 mL of cell culture) for 2 to 4 hours following the manufacturer's protocol (BD, 550583). For surface marker analysis, live cells were resuspended in 1× PBS supplemented with 2% FBS and stained with live/dead staining (BD Bioscience, 564407), BUV395 mouse anti-human CD3 (clone: UCHT1, BD Bioscience, 563546), FITC anti-human CD4 (clone: OKT4, BioLegend, 317408), or APC/Cyanine7 anti-human CD8 (clone: SK1, BioLegend, 344714) at 4°C for 30 minutes. Subsequently, the cells were fixed and permeabilized with Fixation and Permeabilization Solution (554722, BD Bioscience) for 30 minutes and then incubated with PE anti-human TNFα (clone: MAb11, BioLegend, 502909), PE mouse anti-human IFNγ (clone: B27, BD Bioscience, 559327), or PE isotype control (981804, BioLegend) at 4°C for 30 minutes. Data were acquired with an ACEA NovoCyte flow cytometer and analyzed with FlowJo software (version 10).

Mouse ICCA tissue dissociation assays were performed by mechanical and enzymatic digestion in the culture medium containing 0.6 mg/mL collagenase IV (17104019, Gibco), 0.01 mg/mL DNase I (11284932001, Merck), and 2% FBS at 37°C for 1 hour. Single-cell suspensions were resuspended in 36% Percoll (P4937, Sigma) and centrifugated at 300 × g for 5 minutes, followed by red blood cell lysis by using blood lysis buffer (555899, BD) for 10 minutes at room temperature. Cells were resuspended in DMEM supplemented with 2% FBS, and stimulated with Cell Stimulation Cocktail as described above. Cells were then centrifuged at 500 × g for 5 minutes and then incubated with 2.5 μg/mL Fc blocker (BioLegend, 15660) on ice for 20 minutes. For surface marker staining, suspensions were incubated with fluorochrome-labeled antibodies (CD45, CD3, CD4, and CD8) at 4°C for 30 minutes, washed with 1×PBS supplemented with 2% FBS, and centrifuged at 500 × g for 5 minutes. The antibodies included CD45-BV785 (103149; clone: 30-F11; 1:600; BioLegend), CD3-FITC (clone: 17A2; 100203; 1:600; BioLegend), CD4-PerCP-Cy5.5 (clone: RM4-5; 561115; 1:600, BD Biosciences), and CD8-BV605 (clone: 5H10-1; 752632; 1:600; BD Biosciences). For intracellular marker staining, cells were fixed and permeabilized with Fixation and Permeabilization Solution (554722, BD Bioscience) at room temperature for 20 minutes in the dark before incubation with anti-mouse TNFα (clone: MP6-XT22, 506305 BioLegend) or PE isotype control (400408, BioLegend). For the analysis assay, the samples were washed and resuspended in 1×PBS supplemented with 2% FBS. The data were acquired with a five-laser flow cytometer (BD Bioscience) and analyzed with FlowJo software (version 10).

THP1 cells were treated with 80 nmol/L PMA to differentiate the cells into M0 macrophages. The differentiated M0 macrophages were then cocultured with CCLP1 cells, and the cell number was the same as in the previous section. After 48 hours of coculture, cells were isolated and stained with FITC anti-human CD14 (325604, BioLegend) to label macrophages, following the FACS staining methods described previously. After staining, macrophage sorting was performed by Beckman moflo Astrios EQ, and the purity is more than 90%. TRIzol reagent (Invitrogen) was added to extract the total RNA from the cells, according to the manufacturer's instructions. Reverse transcription was performed using PrimeScript reagent (RR047A, TaKaRa) to synthesize cDNA. Quantitative real-time (qRT)-PCR was performed with TB Green Premix Ex Taq (RR820A, TaKaRa) and quantified by a CFX Real-Time PCR Detection System (Bio-Rad), and the template is less than 100 ng. Primers were used as follows: APOE (forward: GCGGCTTGGTAAATGTGCTG, reverse: AATCCCAAAAGCGACCCAGT), C1QB (forward: CAGGGATAAAAGGAGAGAAAGGG, reverse: GGCCGACTTTTCCTGGATTC), and GAPDH (forward: CGGATTTGGTCGTATTGGG, reverse: CTGGAAGATGGTGATGGGATT) was used for normalization, and the 2^–ΔΔCt method was used to determine the relative gene-expression. Three biologically repeated experiments (n = 3) were carried out.

Survival analysis for the training and validation cohorts was performed with the “survival” package from R software (43). The hazard ratio (HR) was determined via univariate Cox regression analysis. Statistical tests were selected based on the specific assumptions relative to the data distribution and its variability. T tests and Kruskal–Wallis tests were performed with the basic function of R software. Sample data were analyzed by a two-tailed Student’ t test to identify statistically significant differences between two groups; the Kruskal–Wallis test was used to identify differences between three groups. The data are represented as the means ± standard error of the mean (SEM). A P value < 0.05 indicated statistical significance.

The data generated in this study have been deposited in iProX (Project number: IPX0003037000).

To comprehensively characterize the molecular features of the ICC tumors, ICC tissues derived from 110 patients (Supplementary Table S1) were procured for proteomics analysis (Fig. 1A). A total of 10,888 proteins were identified and quantified, with low batch variance (Supplementary Fig. S1A). Among the 10,888 proteins, 6,311 proteins were detected in samples with less than 10% missing values and were subjected to further analysis (Supplementary Fig. S1B).

Based on the proteomic matrix of 6,311 proteins, the ICC cohort was divided into three distinct clusters (C1, C2, and C3), with an optimal consensus k-value of 3 (Fig. 1B; Supplementary Fig. S1C and S1D). Each of these clusters was distinguishable in the PCA plot of the patients (Supplementary Fig. S1E). A DEP analysis was performed for each pair of molecular subtypes (Supplementary Fig. S1F–S1H). A Venn diagram was applied to obtain the key DEPs for the three molecular subtypes (Fig. 1C). A total of 85 proteins were identified to be differentially expressed in all comparisons between the three molecular subtypes (Fig. 1C; Supplementary Table S3). Among the 85 proteins, 32 proteins were upregulated in subtype 1, 24 proteins were upregulated in subtype 2, and 29 proteins were upregulated in subtype 3 (Supplementary Table S3). We performed GO analysis for the three protein sets (Supplementary Fig. S1I–S1K). Based on the results from GO analysis, we annotated the three molecular subtypes as “chronic inflammation” (subtype 1), “metabolism” (subtype 2), and “chromatin remodeling” (subtype 3; Fig. 1B and C). Gene set variation analysis (GSVA) was performed on the proteomic data using the related gene sets (ref. 44; Supplementary Table S4), which showed that the inflammatory response score was greatest in the chronic inflammation subtype compared with the other two subtypes (Fig. 1D). Similarly, the metabolism score was greatest in the metabolism subtype, and the chromatin remodeling score was greatest in the chromatin remodeling subtype (Fig. 1E and F). The density distribution of the three ssGSEA scores is shown in Supplementary Fig. S1L. The differentially expressed proteins and pathway alterations are shown in Fig. 1G, indicating the significant different pathway activation in the three molecular subtypes.

When we stratified the ICC patients based on the three subtypes, we observed that ICC patients had significantly different survival rates (Fig. 1H). Specifically, patients belonging to the chronic inflammation subtype exhibited the shortest survival, whereas those with chromatin remodeling subtype survived the longest. No significance was detected between clinicopathologic features and molecular subtypes in the training cohort (Supplementary Table S5). We further validated the molecular subtype in an independent ICC cohort (Supplementary Fig. S2). The proteomics data from the validation cohort were integrated with the training cohort and were subjected to consensus clustering, and the molecular subtypes of samples in the validation cohort were identified. PCA showed distribution of tumor tissues from the validation cohort into the three molecular subtypes (Supplementary Fig. S2A). The clinicopathologic features and molecular subtypes in the validation cohort are shown in Supplementary Table S6. The Kaplan–Meier estimates of OS also suggested the worst prognosis for patients within the chronic inflammation molecular subtype in the independent validation ICC cohort (Supplementary Fig. S2B). Taken together, through proteomic analysis, we were able to stratify ICC into groups with distinct molecular characteristics and patient survival rates.

To further explore the genomic characteristics of the three molecular subtypes, WES was performed using human ICC tumors in the training cohort, and the genomic landscape was depicted for 78 samples (Fig. 2A). ARID1A, KRAS, PBRM1, BRCA2, CREBBP, TP53, KMT2D, ATM, and ATRX were the most frequently mutated genes in ICC tumor tissues (Fig. 2B). A detailed summary of the mutation landscape showed that missense mutations were the most frequent variant classification, and that SNPs were the most frequent variant type (Supplementary Fig. S3A). As expected, RTK-RAS was the most significantly altered pathway in ICC (Supplementary Fig. S3B; ref. 45), and KRAS and BRAF were the most frequently mutated genes in the RTK—RAS pathway (Supplementary Fig. S3C). Nevertheless, the Kaplan–Meier curve and univariate Cox regression indicated an insignificant association between OS and KRAS mutation or RTK—RAS pathway alteration (Supplementary Fig. S3D and S3E).

We then used a random forest–based machine learning method to identify the key gene mutations in each molecular subtype. The GSEA score from Fig. 1D–F was used as the output, and the gene mutations were applied as the input. The feature importance ranking is shown in Fig. 2C–E. KMT2D mutations were found in both chronic inflammation and metabolism subtypes. Among all identified mutations, only mutations in KMT2D were significantly associated with a better prognosis for ICC patients (Fig. 2F). MSK pan-cancer genomic data also verified that patients with KMT2D mutations (KMT2D^mut) had a better prognosis than patients with wild-type KMT2D (KMT2D^wt; Supplementary Fig. S3F). A significantly high proportion of KMT2D mutations was also found in the metabolism subtype (Chi-squared test P < 0.01; Fig. 2G). KMT2D^mut tissue had a greater metabolism ssGSEA score and lower inflammatory ssGSEA score than KMT2D^wt tissues (Fig. 2H and I). We also found that most KMT2D^mut had a missense mutation in exon 39 (Supplementary Table S7). Overall, we identified KMT2D^mut as being associated with the chronic inflammation and metabolism of ICC subtypes.

Next, we dissected the immune landscape of different subtypes of ICC. The abundance of 21 immune cell populations was estimated using the ssGSEA method by assessing the gene signature for each cell type on the proteomic data (Fig. 3A). A significantly greater abundance of macrophages was found in the chronic inflammation subtype (Fig. 3B). IHC staining was performed to validate the results from the proteomic analysis. As expected, more macrophages were found in samples with the chronic inflammation subtype than in the other two subtypes (Fig. 3C and D). We also observed that more macrophages were found in patients in the KMT2D^WT group than in the KMT2D^mut group (Fig. 3E), and IHC staining validated the greater number of macrophages in KMT2D^WT than in KMT2D^mut tissues (Fig. 3F and G).

To further characterize the immune features of the chronic inflammation subtype, we performed scRNA-seq on tumor tissues within the chronic inflammation subtype. A total of 10 cell types, including malignant cells, fibroblasts, endothelial cells and various immune cells, were identified in the ICC tumor microenvironment (TME; Fig. 4A). The featured genes in each cell type were then determined (Supplementary Fig. S4). For instance, SPP1, S100A6, FGG, GDF15, and FGB had high expression in malignant cells. We next focused on TAMs in the chronic inflammation TME because TAMs were found to be a major difference between the three subtypes. Macrophages were further split into two clusters: type 1 and type 2 (Fig. 4B). Several M2 macrophage markers (IL10, IL1B, HLADRA, and CD163) were expressed in both type 1 and type 2 clusters (Supplementary Fig. S5A). The ssGSEA signature scores of M2 and M1 were both greater in type 2 macrophages than in type 1 macrophages (Supplementary Fig. S5B and S5C), suggesting that the traditional M1/M2 classification for macrophages in ICC may not be suitable.

We next characterized the two types of macrophages in several steps. The inflammatory response score was calculated for the two subtypes of TAMs using the GSVA method. Type 2 TAMs showed a significantly greater inflammatory response score than type 1 TAMs (Fig. 4C and D). Differentially expressed genes (DEG) analysis was then performed between the two subgroups of TAMs. APOE, LIPA, CTSD, C1QB, and several other genes were significantly upregulated in type 2 TAMs (Fig. 4E). The upregulated genes in type 2 TAMs were used to perform GO analysis, and we found enrichment for the positive regulation of the inflammatory response and negative regulation of the immune system, suggesting the potential involvement of type 2 macrophages in regulating the inflammatory response in ICC (Fig. 4F). The coexpression of APOE and C1QB in type 2 TAMs was verified by the costaining of APOE/CD68 and C1QB/CD68 (Fig. 4G–I). A higher number of APOE⁺ TAMs was found in the chronic inflammation subtype than in the other two subtypes (Fig. 4J and K). Tissues with wild-type KMT2D exhibited a greater number of APOE⁺ TAMs than tissues with KMT2D mutation (Fig. 4L and M). Kaplan–Meier curve and univariate Cox regression analysis suggested that patients with high expression of APOE and C1QB had significantly worse OS than those in the low expression group (Supplementary Fig. S5D and S5E).

The top 20 upregulated type 2 TAM DEGs from the comparison between type 1 and type 2 TAMs were used to define a type 2 APOE⁺C1QB⁺ TAM signature (Supplementary Table S8), and we performed ssGSEA based on the proteomic data in the training cohort. The chronic inflammation subtype exhibited the highest APOE⁺C1QB⁺ TAM score (Fig. 4N), suggesting that type 2 APOE⁺C1QB⁺ TAMs are the main determinant of the chronic inflammation subtype of ICC. The type 2 APOE⁺C1QB⁺ TAM score was higher in KMT2D^WT patients than in KMT2D^Mut patients (Fig. 4O), and patients with a high type 2 APOE⁺C1QB⁺ TAM signature score exhibited significantly worse OS than those with a low signature score (Supplementary Fig. S5F).

We next asked how KMT2D mutation reshaped TAMs in the ICC TME. First, we screened a series of human ICC cell lines and found that KMT2D was highly expressed by CCLP1 cells, whereas it had low expression by HCCC9810 cells (Fig. 5A; Supplementary Fig. S6A). We thus chose both cell lines and cocultured them with THP1-derived M0 macrophages induced by PMA (46). The expression of APOE and C1QB was determined by IF staining and qRT-PCR (Fig. 5B). IF showed that M0 macrophages cocultured with CCLP1 cells exhibited significantly higher APOE and C1QB expression than M0 macrophages cocultured with HCCC9810 cells (Fig. 5C and D). Macrophages were sorted by FACS, followed by qRT-PCR. M0 macrophages cocultured with ICC cells (HCCC9810 or CCLP1) resulted in a higher expression of APOE and C1QB in cocultivated M0 cells than in M0 cells cultured alone, which confirmed the induction of APOE⁺C1QB⁺ TAMs by ICC malignant cells (Fig. 5E), with M0 macrophages cocultured with CCLP1 cells showing significantly higher APOE and C1QB expression than M0 macrophages cocultured with HCCC9810 cells (Fig. 5D and E).

Next, we investigated how APOE⁺C1QB⁺ TAMs promoted inflammation in the ICC TME. APOE⁺C1QB⁺ TAM and CD4⁺ T-cell ssGSEA scores were calculated using their gene signatures on the proteomic matrix; we found that patients with high APOE⁺C1QB⁺ TAM abundance exhibited higher CD4⁺ T-cell infiltration (Fig. 5F). Multiplex (m)IF on human ICC tissues from the training cohort indicated that ICC tissues with more APOE⁺ TAMs also had greater CD4⁺ T-cell numbers (Fig. 5G and H). We also found many CD4⁺ T cells in the ICC TME-produced TNFα (Fig. 5I).

To investigate whether APOE⁺C1QB⁺ TAMs regulated the inflammatory activity of T cells, we performed a NicheNET analysis on APOE⁺C1QB⁺ TAMs and T cells, modeling their intercellular communication by linking ligands of sender cells to target genes of receiver cells (42). In APOE⁺C1QB⁺ TAMs, IL6, CD40, CXCL16, and several other genes were ranked as the top potential ligands, according to their activity in regulating the target genes of T cells (Fig. 5J). The downstream genes in T cells included AHR, CCL20, CCL5, CSF1, IL7R, IRF1, TNF, TNFRSF1B, and other key important genes involved in the regulation of the inflammatory response of T cells. Because TNF is the target for most ligands from APOE⁺C1QB⁺ TAMs, we thus focused on the effect of APOE⁺C1QB⁺ TAMs on T cells in coculture experiments. KMT2D^high CCLP1 cells and KMT2D^low HCCC9810 cells were cocultured with M0 macrophages and T cells (Fig. 5K). FACS was performed to elucidate the effect of APOE⁺C1QB⁺ macrophages on T cells (gating strategy in Fig. 5L). Our results indicated the upregulation of TNFα in CD4⁺ T cells when cells were cocultured with APOE⁺C1QB⁺ macrophages (induction by the ICC cell lines; Fig. 5M). The coculture of KMT2D^high CCLP1 cells, M0 macrophages, and CD4⁺ T cells resulted in the highest TNFα expression in CD4⁺ T cells (Fig. 5N). Nevertheless, APOE⁺C1QB⁺ macrophages did not improve the percentage of IFNγ⁺CD8⁺ T cells (Supplementary Fig. S6B and S6E). We further checked the effect of APOE⁺C1QB⁺ macrophages on CD8⁺ T cells. APOE⁺C1QB⁺ macrophages did not drive the accumulation of TNFα-secreting CD8⁺ T cells or IFNγ⁺ CD8⁺ T cells (Supplementary Fig. S6C and S6D; S6F and S6G). There was a trend for an increased percentage of TNFα-secreting CD8⁺ T cells and IFNγ⁺ CD8⁺ T cells, but this did not reach significance (Kruskal–Wallis).

Because APOE⁺C1QB⁺ TAMs were the most prominent player driving the chronic inflammation subtype of ICC with the worst prognosis, we tried to target APOE⁺C1QB⁺ TAMs in a mouse ICC tumor model using a blocking CSF1R antibody to deplete TAMs (Fig. 6A). Tumor tissues were harvested, and we found a significantly lower tumor burden in the CSF1R antibody–treated group (Fig. 6B and C). H&E staining was used to assess the properties of the ICC tumor model (Fig. 6D). Ki67 staining suggested lower proliferation in the CSF1R antibody–treated group than in the control group (Fig. 6E). IF staining revealed a lower number of APOE⁺C1QB⁺ TAMs in the CSF1R antibody–treated group than in controls (Fig. 6F and G). The proportion of APOE⁺C1QB⁺ TAMs in all F4/80⁺ cells was quantified. A lower proportion of APOE⁺C1QB⁺ TAMs was found in the CSF1R antibody–treated group than in controls (Fig. 6H). FACS analysis was performed to assess TNFα production by T cells. After CSF1R antibody treatment, APOE⁺C1QB⁺ TAMs were depleted, and less T-cell inflammation was observed, as quantified by the TNF signal (Fig. 6I and J). These results indicated that APOE⁺C1QB⁺ TAMs may be a potential target for immunotherapy against ICC.

ICC is a relatively rare but highly aggressive cancer type that has a low 5-year survival rate. The clinical features and survival outcomes of ICC vary among cohorts (8, 47), thus requiring a novel molecular classification-based prognostic model. However, the relatively rare cases of ICC hinder the construction and validation of a molecular model. In this study, we, therefore, collected a relatively high number of ICC tissues and provided a comprehensive molecular characterization of ICC through system-level multiomics analyses. Three novel molecular subtypes were identified, and the potential underlying molecular mechanisms of the poor prognosis were also explored.

Inflammation has long been related to ICC development and metastasis formation through the modulation of the components of the TME (48, 49). However, the mechanism and effect of inflammation in ICC are still controversial. Daniela Sia and colleagues identified inflammation and proliferation classes as two distinct subtypes in ICC (8). The inflammatory class, defined by the humoral response and imbalance of cytokines, is associated with a more favorable prognosis (8). Nevertheless, several other studies have provided contradictory data and opinions (47, 50, 51). Andersen and colleagues reported two categories of ICC tissues, each containing two subgroups (SGI-SGIV) characterized by distinct gene-expression profiles. Genes differentially expressed between SGI and SGII were mainly involved in the immune response, whereas overexpression of genes involved in the proteasomal activation pathway distinguished SGIII from SGIV (7). Therefore, it is considered that proteasome and anti-inflammatory inhibitors represent a novel therapeutic option for a defined subgroup. Sylvie Job and colleagues constructed another classification system based on the TME of ICC (14). With the help of computational methods, four distinct molecular subtypes, including immune desert, immunogenic (inflamed), myeloid (enriched with M2 macrophages), and mesenchymal, were identified in an ICC cohort (14), which highlighted the importance of macrophages and other immune cell populations in reshaping the TME of ICC. Nevertheless, most of the previous reports focus on bulk transcriptome-based molecular classification systems in ICC, and the underlying mechanisms of the reshaping function of tumor inflammation in ICC warrant further exploration. Consistent with previous studies, we showed that inflammation was the key regulator in ICC. Patients with the chronic inflammation subtype exhibited the worst survival outcome compared with patients in the other two subtypes in our training and validation cohorts. These findings were confirmed by one study using proteomics for molecular subtyping, where four ICC patient subgroups (S1–S4) with subgroup-specific biomarkers were identified (15). The inflammatory and mesenchymal subgroups exhibit a more aggressive phenotype and have a relatively worse prognosis than the other two subgroups (15). Distinct genetic alterations, microenvironment dysregulation, tumor microbiota composition, and potential therapeutics were found in the four subgroups. Taken together, previous studies have mainly focused on molecular subtyping and have opposing opinions on the relationship between inflammatory activity and prognosis in ICC. Therefore, we tried to explore the detailed molecular mechanism of inflammation in ICC at the genomic and single-cell transcription levels.

Through assessment of the genomic landscape, as expected, several key genes, such as ARID1A, KRAS, and PBRM1, were found to be mutated in ICC. The results showed molecular similarities with HCC (52). ARID1A, which only accounts for a 7% mutation rate in the TCGA-HCC cohort, is the most mutated gene in ICC (53). ARID1A performs context-dependent oncogenic and tumor-suppressing functions in liver cancer (54). The detailed mechanism of ARID1A mutation still requires further exploration in ICC. KRAS mutations are associated with poor prognosis in ICC in the Andersen and colleagues cohort (7). However, KRAS mutations did not reach significance in our cohort (P = 0.13). The Chi-squared test also suggested that no significance existed for the association of KRAS mutations and the three molecular subtypes in our cohort. KMT2D is a histone methyltransferase that plays multiple roles in development, cell fate transition, and tumor suppression (55, 56). KMT2D was associated with lower inflammatory activity and mainly occurred within the metabolism subtype by our random forest algorithm. Genes involved in antigen presentation (e.g., PSMA7) and ubiquitination (e.g., PSMA7, PSMA1) were also enriched in the metabolism subtype. Data by Gu and colleagues highlight the correlation between KMT2D and the ubiquitination proteasome system (57). Hence, we considered that KMT2D mutation may be a key event in the metabolism subtype. Our data indicated that KMT2D mutation is associated with a better prognosis than wild-type KMT2D in our ICC cohort, which is supported by an independent study based on Cas9 mutagenic screening (58). One study also indicates that ICC patients harboring KMT2D mutations have higher immunotherapy response rates and prolonged survival outcomes than patients with wild-type KMT2D (59). In vitro experiments showed that the coculture of the KMT2D^high cell line, macrophages, and T cells resulted in the highest TNFα expression in T cells. Hence, we considered that KMT2D mutation is associated with the regulation of inflammatory activities in ICC. Patients harboring KMT2D mutations or bearing ICC tumors that belong to the metabolism subtype may obtain more immunotherapy benefits.

Through ssGSEA of the proteomic data, we found that macrophages were mostly associated with the chronic inflammation subtype of ICC. We thus focused on macrophage populations at the single-cell transcriptome level. scRNA-seq analysis was performed with tumor tissues from the chronic inflammation subtype, and the macrophage populations were further identified. M1 and M2 macrophages have been found in various tumor types with distinct properties (60, 61). M1 macrophages are highly active in inflammatory responses and are commonly considered antitumor macrophages, whereas M2 macrophages can promote tumor cell growth and evasion by dysregulating several key pathways in tumor cells (61). Unlike the traditional M1 and M2 classification, in our study, ICC TAMs were divided into two subtypes with distinct molecular features. Type 2 TAMs showed both higher M1 and M2 scores than type 1 TAMs, suggesting that a new classification method should be used in ICC macrophages. Type 2 TAMs were characterized by high expression of APOE and C1QB, along with a high inflammatory response score, implying the importance of type 2 APOE⁺C1QB⁺ TAMs in reshaping the inflammatory response of the TME. A high abundance of APOE⁺C1QB⁺ TAMs was also found in the chronic inflammation subtype, which once again implied the importance of APOE⁺C1QB⁺ TAMs in the inflammatory response of ICC.

Rheumatoid arthritis (RA) is an autoimmune disease characterized by chronic inflammation in the synovium of the joint tissue, which leads to joint destruction, disability, and shortened life span (62). Interestingly, a study by Zhan and colleagues reveals that 4 subsets of monocytes with distinct molecular signatures exist in RA (63) and that one subset, M3, is enriched with C1QA, C1QB, and other proinflammatory genes. They conclude that the transcriptional profiles of M3 monocytes do not align with known activation states, possibly indicating that M3 represents cell phenotypes tailored to the unique homeostatic needs of the synovium (63). TREM2⁺APOE⁺ macrophages, which are involved in neurodegenerative disorders (64), have been found in tumors, where they display a potentially immunosuppressive role (65). Anti–PD-1–based immunotherapy is more efficient when TREM2⁺APOE⁺ macrophages are targeted by monoclonal antibodies (66). The previous findings and the results of our research confirmed the importance of some subsets of monocyte/macrophage populations (usually characterized by C1QA, C1QB, APOE, and other molecules) in reshaping the inflammatory landscape among different diseases and tumor types.

To further explore the mechanism by which type 2 APOE⁺C1QB⁺ TAMs activated inflammatory activity in the TME of ICC, we explored the interactions of APOE⁺C1QB⁺ TAMs by single-cell transcriptome analysis. We found that APOE⁺C1QB⁺ TAMs promoted TNFα secretion in CD4⁺ T cells through several ligands. To confirm the results of the bioinformatics analysis, we cocultured M0 macrophages, CD4⁺ T cells, and ICC tumor cell lines. FACS revealed that APOE⁺C1QB⁺ TAMs promoted CD4⁺ T-cell inflammation by increasing TNFα secretion. To confirm the results from bioinformatic analysis and in vitro data, we used an in vivo spontaneous ICC murine model. Among target-TAM therapeutics, CSF1R is a promising molecule, and blocking CSF1R leads to significant TAM depletion (67). Pexidartinib, which exhibits selective activity against CSF1R, has been approved for the treatment of adult patients with symptomatic tenosynovial giant cell tumors not amenable to improvement with surgery (68). Nevertheless, the significant toxicity and insufficient effect have limited the application of pexidartinib in several solid tumors. Thus, considering the significant toxicity and the existing compensatory mechanisms, targeting a specific TAM subset may be more advantageous to overcome the drawback of general TAM blockade therapy. We administered anti-CSF1R in the in vivo tumor model, which led to the depletion of APOE⁺C1QB⁺ TAMs and a decrease in T-cell inflammation. Nevertheless, other macrophages/monocytes beyond the APOE⁺C1QB⁺ TAM population were also affected by anti-CSF1R treatment. Developing a novel antibody for targeting the APOE⁺C1QB⁺ TAM population would, therefore, be a promising way to hinder the progression of ICC and improve the immunotherapy response.

To our knowledge, this is the first report to apply scRNA-seq-assisted multiomics analysis to two independent ICC cohorts. Most of the ICC studies are either based on genomic or transcriptomic data in a single center (69, 70). We combined genomic and proteomic information from the in-house ICC cohort to identify key mutations associated with the molecular subtypes of ICC. Our two-center proteomic analysis further confirmed the robustness of the molecular subtyping. Compared with previous studies (14, 71, 72) based on transcriptome data, our scRNA-seq–assisted proteomic analysis might represent substantial progress in the in-depth analysis of the ICC TME at the single-cell level. Nevertheless, there are still some limitations to our study. We illustrated the association between KMT2D mutation and inflammatory activities in the ICC TME. The underlying molecular mechanism is worthy of further exploration. We recruited 151 cases of fresh tumor tissues, but we were not able to perform WES in all of the samples, which may have reduced statistical power in our analyses and, thus, could affect our conclusions. Tissues from two institutions from one country were used to construct the cohort. More centers from different countries could help reduce geographical polymorphisms. Last, the patient number for scRNA-seq was relatively low. It will also be important to show the single-cell landscape for the other two molecular subtypes.

In conclusion, we introduced a novel workflow for the scRNA-seq–assisted molecular subtyping process for ICC. Three molecular subtypes were identified for stratifying ICC patients with different prognoses. In particular, the chronic inflammation subtype associated with the poorest outcome. Mutations in KMT2D were found to be significantly involved in the inflammatory response and metabolic activity of ICC, and an APOE⁺C1QB⁺ TAM subtype showed the capacity to affect T cells, therefore reshaping the inflammatory landscape of the TME in the chronic inflammation subtype of ICC.

T. Guo reports other support from Westlake Omics Inc, nonfinancial support from Pressure Biosciences Inc, Thermo Fisher, Sciex, and Bruker outside the submitted work. No disclosures were reported by the other authors.

X. Bao: Conceptualization, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. Q. Li: Validation, investigation, visualization, methodology. J. Chen: Investigation, visualization, methodology. D. Chen: Data curation. C. Ye: Data curation, methodology. X. Dai: Validation, investigation, methodology. Y. Wang: Validation, investigation, visualization, writing–original draft. X. Li: Data curation, validation, investigation. X. Rong: Resources, data curation, investigation. F. Cheng: Methodology. M. Jiang: Methodology. Z. Zhu: Methodology. Y. Ding: Methodology. R. Sun: Validation, visualization, methodology. C. Liu: Visualization, methodology. L. Huang: Methodology. Y. Jin: Methodology. B. Li: Formal analysis. J. Lu: Formal analysis. W. Wu: Formal analysis. Y. Guo: Formal analysis. W. Fu: Resources, data curation. S. Langley: Methodology. V. Tano: Methodology. W. Fang: Supervision, investigation, writing–review and editing. T. Guo: Data curation, supervision. J. Sheng: Formal analysis, supervision, project administration. P. Zhao: Supervision, project administration. J. Ruan: Supervision, funding acquisition, investigation.

This work was supported in part by the National Natural Science Foundation of China 82101830 (X. Bao), 82102817 (X. Dai), 81874173 (J. Ruan), 81472346 (P. Zhao), 82074208 (P. Zhao), 81972492 (T. Guo), by the National Key R&D Program of China grant 2019YFA0803000 (J. Sheng), 2020YFE0202200 (T. Guo), by the Young Investigator Research Program 588020-D01907 (M. Jiang), and by the Natural Science Foundation of Zhejiang Province LY20H160033 (P. Zhao) and LQ22H160041 (X. Dai). We thank OE Biotech Co., Ltd (Shanghai, China) for providing single-cell RNA-seq.

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.