Multimodule characterization of immune subgroups in intrahepatic cholangiocarcinoma reveals distinct therapeutic vulnerabilities.

BACKGROUND: Immune microenvironment is well recognized as a critical regulator across cancer types, despite its complex roles in different disease conditions. Intrahepatic cholangiocarcinoma (iCCA) is characterized by a tumor-reactive milieu, emphasizing a deep insight into its immunogenomic profile to provide prognostic and therapeutic implications.
METHODS: We performed genomic, transcriptomic, and proteomic characterization of 255 paired iCCA and adjacent liver tissues. We validated our findings through H&E staining (n=177), multiplex immunostaining (n=188), single-cell RNA sequencing (scRNA-seq) (n=10), in vitro functional studies, and in vivo transposon-based mouse models.
RESULTS: Integrated multimodule data identified three immune subgroups with distinct clinical, genetic, and molecular features, designated as IG1 (immune-suppressive, 25.1%), IG2 (immune-exclusion, 42.7%), and IG3 (immune-activated, 32.2%). IG1 was characterized by excessive infiltration of neutrophils and immature dendritic cells (DCs). The hallmark of IG2 was the relatively higher tumor-proliferative activity and tumor purity. IG3 exhibited an enrichment of adaptive immune cells, natural killer cells, and activated DCs. These immune subgroups were significantly associated with prognosis and validated in two independent cohorts. Tumors with KRAS mutations were enriched in IG1 and associated with myeloid inflammation-dominated immunosuppression. Although tumor mutation burden was relatively higher in IG2, loss of heterozygosity in human leucocyte antigen and defects in antigen presentation undermined the recognition of neoantigens, contributing to immune-exclusion behavior. Pathological analysis confirmed that tumor-infiltrating lymphocytes and tertiary lymphoid structures were both predominant in IG3. Hepatitis B virus (HBV)-related samples tended to be under-represented in IG1, and scRNA-seq analyses implied that HBV infection indeed alleviated myeloid inflammation and reinvigorated antitumor immunity.
CONCLUSIONS: Our study elucidates that the immunogenomic traits of iCCA are intrinsically heterogeneous among patients, posing great challenge and opportunity for the application of personalized immunotherapy. © Author(s) (or their employer(s)) 2022. Re-use permitted under CC BY-NC. No commercial re-use. See rights and permissions. Published by BMJ.

Intrahepatic cholangiocarcinoma (iCCA) represents the second most common type of primary liver cancer, which is generally featured by a highly desmoplastic tumor microenvironment (TME). Sequencing efforts have improved the understanding of the molecular pathogenesis of iCCA; however, the comprehensive landscape of iCCA microenvironment phenotype remains largely unexplored.

Our multimodule study reveals three distinct immune subgroups (IG1, IG2, and IG3) with specific immunological and non-immune profiles associated with different patient outcomes. Myeloid inflammation/KRAS mutation, antigen presentation/proliferation, and immune checkpoint could serve as potential therapeutic targets for the treatments of IG1, IG2, and IG3-related iCCA, respectively. Hepatitis B virus (HBV) infection could probably alleviate myeloid inflammation in iCCA, thus activates antitumor immunity and improves the efficacy of immunotherapies.

These findings extend our understanding of the immunological profile of iCCA. Defined immune subgroups provide novel insights into the heterogeneous TME and underlying therapeutic strategies.

Introduction

Intrahepatic cholangiocarcinoma (iCCA) accounts for 10%–15% of all primary liver cancer, with a gradual increase over the past decades.1 Despite the high postoperative recurrence rate, surgery is prior for patients diagnosed at an early stage. Unfortunately, most patients are diagnosed at advanced stages, when only palliative chemotherapy is available, with a median survival shorter than 12 months. Emerging molecular targeted therapies against FGFR2 fusion and IDH1/2 mutations show modest survival benefits in a limited subgroup of patients with iCCA.2 Hence, identification of novel and effective therapeutics remains an urgent need.

In general, iCCA is featured by a highly desmoplastic tumor microenvironment (TME) with excessively infiltrating immune and stromal cells. Available evidence indicates that the inflammatory microenvironment is associated with iCCA progression and poor prognosis,3 highlighting TME as a potential target for the treatment of iCCA. Indeed, recent phase II/III clinical trials confirm that immune checkpoint inhibitors (ICIs) plus chemotherapy is promising in patients with iCCA.4 5 However, the responses and survival benefits are only limited to a small subset of patients, probably owing to the heterogeneous and complex TME. Therefore, a better understanding of the immunogenomic profile is urgently needed for designing novel immunotherapies to improve responses and patient outcomes.

Patient stratification based on immune profile can help delineate immunologically heterogeneous tumors and identify suitable patients for immunotherapy.6 7 We and others have applied bulk-tumor transcriptomic or proteomic profiles to define clinically relevant iCCA subtypes,8–13 collectively converging on intertumor immune heterogeneity. Nevertheless, the underlying cancer-cell-intrinsic properties that dictate the diversity in immune landscape remain less characterized. Integrated multidimensional and multimodule analyses of preclinical models and patient samples are needed and expected to deeply characterize the immunological portraits and identify novel personalized immunotherapeutic opportunities in each molecular subtype of iCCA.

In this study, we analyzed the transcriptomic profiles of TME using the largest single-centered iCCA cohort and classified 255 patients into three immune subgroups (IG1, IG2, IG3). By integrating multiple complementary approaches, including whole-exome sequencing (WES), proteomics, single-cell RNA sequencing (scRNA-seq), immunostaining, and mouse models, we delineated the distinct genetic and molecular features among these subgroups. IG1 was an immune-suppressed subgroup with the worst prognosis, harboring the highest myeloid infiltration and KRAS mutations; IG2 was characterized by tumor proliferation and antigen presentation defects; IG3 was featured by the activation of antitumor immunity and enrichment of tertiary lymphoid structures (TLSs), with a good prognosis. The substantial differences in immunological profile suggested that each immune subgroup need specific therapeutic strategies.

Methods

The iCCA samples from 255 patients used for this study were collected from Zhongshan Hospital of Fudan University. Patients treated between December 2014 to October 2018 were randomly selected from our prospectively collected tissue bank. Detailed clinicopathological features are summarized in online supplemental table S1.

Full details of multi-omics analysis, multiplexed immunostaining of tissue microarray (TMA), pathological examination, mouse model construction, scRNA-seq analysis, functional experiments, and statistical analysis are provided in the online supplemental materials and methods.

Results

Immune-centered classification of iCCA

To characterize the immune profile of iCCA with prognostic significance, multi-omics data derived from the Fudan University (FU-iCCA) cohort (n=255) were obtained (online supplemental figure S1A and online supplemental table S1).11 Consistent with iCCA’s inflammatory and desmoplastic features, our cohort had a relatively higher immune infiltration among The Cancer Genome Atlas (TCGA) Pan-Cancer datasets (online supplemental figure S1B). Unsupervised hierarchical clustering based on the mRNA levels of 170 prognostic genes from TME signature14 revealed three Immune subGroups, namely IG1, IG2, and IG3 (figure 1A; online supplemental figure S1C–E and online supplemental table S1). Pairing the transcriptomic and proteomic data from this signature created 103 mRNA–protein pairs, which showed an overall positive correlation (median r=0.6, multiple-test adjusted p<0.01; online supplemental figure S1F). Of note, the abundance of corresponding 103 proteins also clustered these patients into similar immune subgroups with a comparable overlap (online supplemental figure S1G).

Figure 1

Immune-centered classification of intrahepatic cholangiocarcinoma (iCCA). (A) Unsupervised hierarchical clustering of the Fudan University (FU)-iCCA (n=255) cohort based on the 170 genes. The enriched pathways in IG1 and IG3 were shown on the right (χ2 test and adjusted analysis of variance (ANOVA)). (B) Kaplan-Meier curves for overall survival (upper panel) and recurrence-free survival (bottom panel) among the three immune subgroups (log-rank test). (C) HR and p values of immune subgroups and covariates in multivariate analysis for overall survival. (D) The relative abundance of tumor microenvironment (TME) subsets identified from xCell (Fisher’s exact test or ANOVA). (E) Univariate analysis for overall survival of each cell subset. (F) Pathological tumor-infiltrating lymphocyte (TIL) estimates (n=177) plotted for each tumor sample against indicated signatures. (G) Representative immunostaining images showing the indicated immune subsets. (H) Quantification of staining intensities for the indicated immune markers (n=188, ANOVA).

As expected, the immune subgroups significantly differed in overall survival (p=1e-09, log-rank test) and recurrence-free survival (p=4e-07, log-rank test), where IG1 had the worst and IG3 showed the best prognosis (figure 1B). Multivariable Cox analysis further authenticated the immune subgroups as a significantly independent prognosticator (HR, 2.4; 95% CI, 1.33 to 4.5; p=4e-03; figure 1C). Carbohydrate antigen 19-9 (CA19-9, p=6.8e-11) and carcinoembryonic antigen (CEA, p=1.2e-06) were significantly higher in IG1, whereas albumin (ALB, p=9.9e-03) levels increased gradually from IG1 to IG3 (online supplemental figure S1H). IG1 had more tumors with necrosis (p=0.012) and advanced tumor node metastases (TNM) stage (p=2.7e-04) (figure 1D and online supplemental figure S1H). Of note, lymph node metastasis (p=0.019) and vascular invasion (p=0.038) were predominated in IG1, whereas endothelial cells were enriched in IG3. Increased lymphangiogenesis and reduced angiogenesis could make iCCA more invasive and impair treatment efficacy.15 Thus, immunostaining confirmed that compared with IG3, where extensive CD31+ blood vessels were observed, IG1 showed a marked increase in podoplanin positive (PDPN+) lymphatic vessels (online supplemental figure S1I).

Within this prognostic TME signature, IG1 was enriched with pathways of neutrophil degranulation and epithelial–mesenchymal transition (EMT), suggesting an inflammatory and invasive phenotype. In contrast, leucocyte-mediated cytotoxicity and lymphocyte-mediated immunity pathways were predominant in IG3, indicating an immune-active milieu (figure 1A). As expected, xCell14 deconvolution of the RNA-seq data showed that IG1 was dominant by myeloid cells such as neutrophils and immature DCs (iDC), whereas IG3 was abundant in adaptive immune cells, natural killer (NK) cells, and activated DCs (aDC). Immune infiltration in IG2 was relatively low and heterogeneous, accompanied by the scarcity of stromal cells, which might account for the highest tumor purity (p=2.2e-04; figure 1D). Generally, the abundance of neutrophils (p=0.011) and iDC (p=0.045) predicted poor survival, whereas the accumulations of lymphocytes, such as CD4+ T cells, CD8+ T effector/memory (Tem) cells, and B cells, were associated with good prognosis (log-rank test; figure 1E).

Further pathological examination confirmed that tumors in IG3 contained more tumor-infiltrating lymphocyte (TIL) estimates (figure 1D; online supplemental figure S1J, K and online supplemental table S1). Indeed, TIL estimates correlated negatively with tumor purity (p=5.3e-03) and positively with immune deconvolution (CD8+ T, p=1.6e-06; CD4+ T, p=5.6e-05; B cell, p=3.8e-04; figure 1F). Similarly, multiplex immunostaining on TMAs confirmed that the intensities of CD15, CD66b, and CD1a were the strongest in IG1, whereas CD3, CD8, and CD20 were primarily expressed in IG3 (figure 1G and H; online supplemental figure 1L and online supplemental table S1), remarkably consistent with the xCell estimation.

We then performed unsupervised clustering based on the overall mRNAs with the top 25% median absolute deviation and identified three subgroups that were somehow concordant with the immune subgroups (figure 1A). Specifically, pairing immune subgroups and overall mRNA subgroups (RGs) showed that IG1 and IG3 highly overlapped with RG1 and RG2, respectively. Leveraging our TME gene signature, we also classified Jusakul et al’s and Job et al’s iCCA cohorts and observed similar patient subgrouping and survival stratification (online supplemental figure 2A–D). Likewise, clustering our RNA-seq data with the transcriptomic signatures from four previous studies8–10 12 resulted in a significant concordant patient allocation (online supplemental figure 2E–H).

Distinct molecular features among the three immune subgroups

To further illustrate the key biological processes in patient subgroups, we performed single sample Gene Set Enrichment Analysis (GSEA) using the 50 hallmark gene sets16 (online supplemental figure S3A and online supplemental table S2). Analysis of principal-component feature showed that PC1 could be explained by tumor purity with immune and stromal gene set vectors opposite to oncogenic and other vectors, whereas PC2 discriminated immune from stromal gene sets (online supplemental figure S3A–C). Unsupervised hierarchical clustering further verified the association between signaling pathways and immune subgroups, as well as separated clustering of immune and other gene sets (online supplemental figure S3D).

Pathway enrichment analysis incorporating both transcriptomics and proteomics data demonstrated distinct molecular features among the immune subgroups (figure 2A). IG1 was high in inflammatory response, interleukin 6 (IL-6)/janus kinase (JAK)/signal transducer and activator of transcription 3 (STAT3) signaling, and tumor necrosis factor-α (TNF-α)/nuclear factor kappa-B (NF-κB) signaling, indicating that these tumors may be driven by smoldering inflammation. Hypoxia and glycolysis pathways, which have been implicated in myeloid cell-mediated tumor immune evasions,17 18 were also predominant in IG1. IG2 were characterized by enrichment of cell-cycle transcriptional programs (G2M checkpoint, E2F targets, Myc targets) and DNA repair pathway, demonstrating the high proliferative activity. Considering the highest tumor burden in IG2, we adjusted the tumor purity19 and still found enriched proliferation-associated pathways in IG2 (online supplemental figure S3E and online supplemental table S2). IG3 was associated with downregulation of glycolysis and hypoxia pathways, as well as upregulation of the oxidative phosphorylation (OXPHOS) pathway, which has been highlighted as the critical aspects of immune activation20 (figure 2A).

Figure 2

Distinct molecular features among the three immune subgroups. (A) Altered pathways at transcriptome and proteome levels among the three immune subgroups. (B) Heatmaps depicting mRNA and protein levels of various markers involved in cancer-promoting or cancer-inhibitory inflammation (adjusted analysis of variance (ANOVA)). (C) Heatmaps depicting mRNA and protein levels of S100A family members (adjusted ANOVA). (D) Signature scores among the three immune subgroups (adjusted ANOVA). (E) Comparisons of COX-2 mRNA levels among the three immune subgroups (Wilcoxon rank-sum test). (F) Kaplan-Meier curves for overall survival based on mRNA abundance of COX-2 (log-rank test). (G and H) Pearson correlation coefficient of COX-2 mRNA level with COX-2-associated inflammatory signature (COX-IS) genes (G) and the xCell enrichment scores of indicated cell subtypes (H).

We further focused on cancer-promoting and cancer-inhibitory inflammation. IG1 was enriched in protumor inflammatory cytokines and chemokines, including IL-1A, IL-1B, IL-6, CXCL3, CXCL5, and CCL18 (figure 2B). The mRNA and protein levels of most S100A family members, including S100A8/A9, were upregulated in IG1, which may be involved in cancer-promoting inflammation21 (figure 2C). Concurrently, increased mRNA and protein levels of CD8A, GZMA, and PRF1 confirmed the upregulation of cancer-inhibitory inflammation in IG3 (figure 2B). Moreover, a panel of immune signatures18 22–24 collectively showed that IG3 displayed elevated T-cell survival and T-effector signatures relative to other subgroups, whereas IG1 exhibited elevated myeloid-derived suppressor cells (MDSC), M2/M1 macrophage, myeloid inflammation, and proinflammation signatures (figure 2D and online supplemental table S2). Overall, our stratification identified two immune subgroups with antagonistic inflammatory phenotypes.

Recently, the cyclooxygenase 2 (COX-2)/prostaglandin E2 (PGE2) axis has been reported to be involved in regulating the antagonistic inflammatory phenotypes in a Pan-Cancer study.25 Indeed, the transcript levels of PTGS2, encoding for COX-2, were significantly upregulated in IG1 (figure 2B, E), consistent with the excessive infiltration of PGE2-associated neutrophils and iDCs.26 27 As expected, patients with higher COX-2 mRNA levels had a worse prognosis (p=1e-04, log-rank test; figure 2F). Strong positive correlations between COX-2 mRNA levels and cancer-promoting factors28 including IL-1A, IL-1B, and IL-6 were detected (figure 2G and online supplemental table S2). Infiltration of neutrophils and iDCs positively correlated with COX-2 mRNA level, whereas plasma cells, CD4+ T cells, and CD8+ Tem cells showed negative correlations (figure 2H and online supplemental table S2). The COX-2-associated inflammatory signature (COX-IS) was a predictor for response to ICIs in multiple cancers,28 and COX-IS decreased sequentially from IG1 to IG3 (online supplemental figure S3F). Additionally, IG3 had the highest cytolytic score29 and immunophenoscore30 (online supplemental figure S3G, H), further supporting the hypothesis that IG3 could be more likely to respond to ICIs.

KRAS mutation-associated myeloid inflammation was dominant in IG1

We further investigated genomic alterations that could be associated with immunological features. IG2 had more frequent somatic copy number variations (SCNAs) (online supplemental figure S4A), in agreement with previous reports that SCNA burden negatively correlated with immune cell infiltration31 (online supplemental figure S4B). When selecting subtype-specific mutated genes, we found that nearly half of tumors in IG1 were KRAS mutants (figure 3A, B and online supplemental table S3). KRAS mutation was negatively associated with most xCell-derived adaptive immune and NK cell signatures, but positively correlated with iDC and neutrophil signatures (figure 3C and online supplemental table S3). GSEA at mRNA and protein levels both identified neutrophil degranulation to be the pathway most strongly associated with KRAS mutation (figure 3D and online supplemental table S3). Consistent with previous report,32 KRAS mutation significantly and positively correlated with several inflammatory molecules, such as COX-2, IL-1A, IL-1B, IL-6, and S100A8/A9, resembling the immunological feature of IG1 (figure 3E and online supplemental figure S4C). Tumors with KRAS mutation also displayed elevated myeloid inflammatory, COX-IS, and MDSC signatures, as well as decreased cytolytic signature (figure 3F, online supplemental figure 4D–F and online supplemental table S3). In Jusakul’s cohort,33 KRAS mutation also led to elevated neutrophil infiltration, COX-2/S100A8 expression, and neutrophil degranulation pathway (online supplemental figure S4G–J).

Figure 3

The association between KRAS mutation and myeloid inflammation. (A) Genes with significant differences in mutational frequency among the three immune subgroups (Fisher’s exact test). (B) Somatic variants along KRAS protein in the Fudan University-intrahepatic cholangiocarcinoma (FU-iCCA) cohort. (C) Association between mutation profiles and immune/stromal signatures from xCell. Only significantly associations were shown (Wilcoxon rank-sum test). (D) Pathway enrichment analysis based on differentially expressed genes (left panel) and proteins (right panel) that associated with KRAS mutations. (E, F) COX-2 mRNA level (E) and myeloid inflammation score (F) in tumor samples with and without KRAS mutation (Wilcoxon rank-sum test). (G, H) Immunohistochemistry of COX-2, Ly6G, S100A8, and S100A9 (G) and quantification of their staining intensities (H). Representative data of triplicate experiments (mean±SEM). Scale bar, 100 µm. (I) Boxplot showing the fractions of indicated tumor microenvironment (TME) composition in AY (n=6) and AYK (n=6) tumors (analysis of variance (ANOVA)). (J) Violin plot showing the exhausted score in T/natural killer (NK) cells (left panel) and myeloid inflammation and s100a8/s100a9 expression in neutrophils (right panel) in AY and AYK samples (ANOVA).

Considering that PTGS2 was the most significantly upregulated gene associated with myeloid inflammation both in IG1 and KRAS mutant tumors, we hypothesized that KRAS mutation may upregulate COX-2 to produce PGE2, which promoted myeloid inflammation. Indeed, ectopic expression of KRAS in iCCA cells significantly upregulated the mRNA and protein levels of COX-2 and COX-1, as well as increased PGE2 levels in the culture supernatant (online supplemental figure 4K–O). Transwell assay revealed that, compared with KRAS, the culture supernatants from KRAS significantly enhanced neutrophil migration (online supplemental figure S4P), which could be hampered by COX-2 inhibitor, celecoxib (online supplemental figure S4Q). KRAS culture supernatant-treated neutrophils showed stronger immunosuppressive activity, as assessed by IL-2 production and apoptosis of preactivated Jurkat T cells (online supplemental figure S4R, S).

We next generated iCCA mouse models using oncogenic drivers myr-AKT and YAPS127A in combination with KRASWT (AY) or KRASG12D (AYK) via the sleeping beauty system and hydrodynamic tail vein injection34 (online supplemental figure S5A). Immunohistochemistry confirmed that AYK, but not AY tumors, were conspicuously infiltrated with neutrophils (positive for Ly6G, S100A8, and S100A9), while devoid of CD8+ T cells and B cells (figure 3G, H and online supplemental figure S5B, C). scRNA-seq analyses revealed that the proportions of neutrophil and macrophage were significantly higher in AYK tumors, whereas T/NK cells, B cell, and stromal cells were more abundant in AY tumors, consistent with multi-omics’ results (figure 3I and online supplemental figure 5D–F). Neutrophils in AYK tumors showed higher myeloid inflammation and elevated expression of inflammatory markers like S100a8, S100a9, Ptgs2, and Il1b, accompanied by stronger exhausted phenotypes of T/NK cells (figure 3J and online supplemental figure S5G). Of note, KRAS mutation upregulated the expression of COX-2 both in tumor cells and the microenvironment (figure 3G, H and online supplemental figure S5G). Thus, excessive infiltration of neutrophils may account for KRAS mutation-associated myeloid inflammation. Overall, KRAS mutation had a critical role in subverting antitumor immunity to myeloid inflammation-associated immunosuppression (online supplemental figure S5H).

Antigen presentation defects are associated with immune exclusion in IG2

Defects in antigen presentation could interrupt neoantigen recognition and lead to tumor immune evasion, irrespective of tumor mutation burden (TMB). We thus explored the underlying roles of these avenues of immune evasion among the three immune subgroups (figure 4A). Despite relatively higher TMB and tumor neoantigen burden (TNB) in IG2, defects in antigen presentation through loss of heterozygosity in human leucocyte antigen (HLA LOH, p=2e-03) or alterations of the other components of antigen presentation machinery (APM, p=0.065) were frequently observed in IG2 (figure 4A, B; online supplemental figure S6A and online supplemental table S4). Such defects may hamper the cross-priming of naive T cells and subsequent tumor recognition by primed T cells, which may explain the reduction in T cell receptor/B cell receptor (TCR/BCR) diversity and RNA-seq reads mapping to VDJ loci in IG2 (online supplemental figure 6B–D).

Figure 4

Antigen presentation defects occur frequently in IG2. (A) Antigen presentation machinery (APM) defects in each immune subgroup (Pearson’s χ2 test or Fisher’s exact test). (B, C) Comparison of tumor mutation burden (TMB) (B) and immunoediting score (C) among the three immune subgroups (Wilcoxon rank-sum test). (D) Immunoediting score were plotted by HLA LOH status and immune classification (Wilcoxon rank-sum test). (E, F) Kaplan-Meier curves for overall survival (E) and recurrence-free survival (F) based on loss of heterozygosity in human leukocyte antigen (HLA LOH) status (log-rank test). (G) Model illustrating how HLA LOH may lead to immune escape in IG2.

We further calculated the immunoediting score for each tumor by quantifying the ratio of observed to expected numbers of predicted HLA class-I binding neoantigens, with a score of <1 indicating the presence of immunoediting.35 Immunoediting scores decreased continuously from IG1 to IG3 (figure 4C and online supplemental table S4). Within each immune subgroup, only IG2 tumors with HLA LOH had a significantly higher immunoediting score, suggesting that tumor subclones in IG2 may be prone to HLA LOH-associated immune evasion (figure 4D). Overall, 15% (n=38) of patients exhibited HLA LOH, which had decreasing trends of overall survival (p=0.069, log-rank test) and recurrence-free survival (p=0.024, log-rank test) compared with the HLA intact patients (figure 4E, F).

The expression of immune checkpoints may reflect a cancer-adaptive immune response to an active immune system. We found a subgroup-specific pattern of these molecules, where IG2 had relatively few exclusively upregulated immune checkpoints (online supplemental figure S6E). HLA LOH correlated with some immune checkpoint molecules, such as programmed death ligand 1, and refined TMB as a biomarker of ICIs.36 Comparing the multi-omics profiles between tumors with and without HLA LOH revealed differentially expressed mRNA and proteins that regulated metabolic reprogramming, cell-cycle transcriptional programs, and immune processes. HLA LOH was enriched with higher TMB/TNB, higher tumor purity, decreased interferon gamma pathway, and decreased major histocompatibility complex class I signatures. Likewise, T-cell survival and T-effector signatures were significantly reduced in tumors with HLA LOH (online supplemental figure S6F). Therefore, APM defects, especially HLA LOH, might involve in immune anergy and afterwards immune-exclusion despite the relatively high TMB in IG2 (figure 4G).

Tertiary lymphoid structures are related to antitumor immunity in IG3

TLSs provide a pivotal microenvironment for generating antitumor immune responses. We used a previous TLS signature37 to estimate TLS distribution and found that the mRNA levels of those genes, together with estimated TLS scores, were relatively upregulated in IG3 (figure 5A, B and online supplemental table S5). Pathological examination also confirmed that intratumoral TLSs were significantly accumulated in IG3 (figure 5C, D). Patients with both primary follicles (FL-I) and secondary follicles (FL-II) had higher estimated TLS scores compared with those with only Aggregates (Agg) and no TLSs (TLSs-), confirming the accuracy of TLS score (figure 5E and online supplemental figure S7A). Consistent with our latest report,38 intratumoral TLSs predicted a better prognosis (figure 5F–H), and TLSs+ (Agg, FL-I, and FL-II) patients were characterized by early TNM stage (p=0.018), lack of microvascular invasion (p=0.033), absence of intrahepatic metastasis (p=0.011), and small tumor size (p=8.1e-04; online supplemental figure S7A and online supplemental table S5).

Figure 5

Tertiary lymphoid structures (TLSs) are enriched in IG3. (A) TLS-associated features in each immune subgroup (adjusted analysis of variance (ANOVA)). (B) Comparison of TLS scores among the three immune subgroups (Wilcoxon rank-sum test). (C) Histological appearance of intratumoral TLSs. (D) Comparison of H&E based intratumoral TLS maturation degrees among the three immune subgroups (Fisher’s exact test). (E) Comparison of TLS scores among H&E-based intratumoral TLS maturation degrees (Wilcoxon rank-sum test). (F–H) Kaplan-Meier curves for overall survival based on TLS score (F), H&E-based intratumoral TLSs (G), and H&E based intratumoral TLS maturation degrees (H) (log-rank test). (I) Diagram showing the multi-omics profiles of regulators of TLS formation and function.

Multi-omics data indicated that TLSs+ tumors showed specific downregulation of pathways of G2M checkpoint, inflammatory response, and EMT, whereas upregulated ones in metabolism and peroxisome (online supplemental figure S7A). As expected, T cells and B cells, as well as coinhibitors, were enriched in TLSs+ tumors, suggesting an immune-activated phenotype (online supplemental figure S7B, C). Moreover, both cytolytic scores and COX-IS predicted a better response to immunotherapy for TLSs+ iCCA patients (online supplemental figure S7D, E). Thus, immunosurveillance in IG3 might be partially attributed to intratumoral TLSs (figure 5I), which predicted a better response to immunotherapy in this subgroup.

HBV infection is negatively associated with myeloid inflammation in iCCA

Hepatitis B virus (HBV) infection is the well-known risk factor for iCCA; however, little is known about the immunogenomic profile of HBV-related iCCA. We found that patients with HBV infection (n=68) were under-represented in IG1 (p=0.022), mutually exclusive with KRAS mutation (p=3.3e-03) and coexisted with TP53 mutation (p=0.036) (figure 6A; online supplemental figure S8A and online supplemental table S6). HBV-positive tumors displayed an elevated cytolytic score, decreased myeloid inflammation signature and decreased COX-2 mRNA levels (figure 6B). Recent study proposed that ICIs showed efficacy for patients with HBV-positive hepatocellular carcinoma (HCC), rather than non-viral HCC.39 Considering the negative impact of myeloid inflammation on immunotherapy,40 the COX-IS score predicted worse responses to ICIs in HBV-negative subgroup (figure 6B). The transcriptome of patients with HBV-positive iCCA resembled the expression profile of other solid tumors that responded to ICIs41 (figure 6C). Thus, the decreased myeloid inflammation observed in HBV-positive iCCA may define a population that could benefit from ICIs.

Figure 6

Influences of hepatitis B virus (HBV) infection on the immune microenvironment of intrahepatic cholangiocarcinoma (iCCA). (A) Associations of HBV infection with indicated features (analysis of variance (ANOVA), Pearson’s χ2 test, or Fisher’s exact test). (B) Comparisons of representative signatures in patients with and without HBV infection (Wilcoxon rank-sum test). (C) Subclass mapping with other solid tumors treated with anti-PD-1 mAb (anti-programmed cell death protein 1 monoclonal antibody). (D) The proportion of indicated immune cells in tumor samples. (E) Uniform Manifold Approximation and Projection (UMAP) plot showing the subtypes of T/natural killer (NK) cells. (F) Boxplot showing the fractions of T/NK cells (ANOVA). (G) UMAP plot showing the subtypes of myeloid derived cells. (H) Boxplot showing the fractions of myeloid-derived cells (ANOVA). (I) Violin plot indicating the mRNA levels of COX-2 (upper panel) and proinflammatory scores (bottom panel) in myeloid subgroups. (J, K) Pathway enrichment analysis based on differentially expressed genes that associated with HBV infections in myeloid cells (J), CD4+ T cells, NK cells, and macrophages (K). (L) Kaplan-Meier curves of overall survival based on HBV infection status (log-rank test).

We additionally collected tumors and adjacent nontumor tissues from five HBV-positive and five HBV-negative iCCAs for scRNA-seq analyses (online supplemental figures S1A and S8B). Increased proportions of T/B/NK cells and reduced proportions of myeloid cells were observed in HBV-positive tumors (figure 6D), but HBV infection did not impact the proportions of each cluster in adjacent nontumor tissues (data not shown). We then classified the 10 patients into the closest immune subgroups according to the enriched gene sets (online supplemental figure S8C). At the CI of 0.7, four patients were excluded, with one (P06T), three (P03T/P09T/P10T), and two (P13T/P15T) samples grouping into IG1, IG2, and IG3, respectively. P06T, the only one remaining in IG1, was also the only patient with KRAS mutation among all the 10 iCCAs. The mRNA levels of inflammatory molecules, including COX-2, IL1A, S100A8, and S100A9, were higher in tumor or myeloid subsets in P06T (online supplemental figure S8D). The mRNA levels of cytolytic molecules, such as CD8A, GZMA, and PRF1, were downregulated in T cells in patient P06T, confirming the immunosuppressive TME in IG1 (online supplemental figure S8D).

Next, we performed unsupervised clustering of T/NK cells and myeloid cells. A total of 13 clusters emerged within T/NK cells (figure 6E and online supplemental figure S8E). The proportions of inflammatory CD4 T_CD15442 were significantly higher in IG1, whereas the proportion of cytolytic NK_GNLY cluster slightly increased from IG1 to IG3 (online supplemental figure S8G). Both CD4_SOCS343 and NK_CD16044 were enriched in HBV-positive tumors, implying that HBV infection alleviated inflammatory responses and enhanced cytolytic activities (figure 6F). The reclustering of myeloid lineage revealed 14 populations (figure 6G and online supplemental figure S8F). Notably, Macro_IL1B was the myeloid cluster with the most significant decrease in proportion after HBV infection, characterized by the exclusively high expression of IL1A, IL1B, IL6, and TNF (figure 6H and online supplemental figure S8F). COX-2 mRNA levels and the proinflammatory signatures were significantly higher in Macro_IL1B than the other myeloid clusters, demonstrating that Macro_IL1B was the primary source of myeloid inflammation (figure 6I). Indeed, myeloid cells from HBV-related iCCA showed downregulation of inflammatory pathways, such as complement and coagulation cascades and interferon alpha response (figure 6J). A smaller fraction of cycling cells was observed in HBV-positive than in HBV-negative samples, reflecting the relatively decreased proliferative capacities (figure 6H). Altogether, HBV-positive patients harbored a unique TME, with increased proportions of CD4_SOCS3 and NK_CD160 clusters, and reduced Macro_IL1B cluster and myeloid cell proliferation, which contributed to antitumor immunity.

Differentially expressed gene (DEG) analysis revealed that the top downregulated genes in macrophages of HBV-positive tumors included CCL13, CTSC, FCGR2A, C1QA, and C1QB (online supplemental figure S8H, I). NK cells in HBV-positive tumors upregulated several cytolytic-related genes, such as KLRC2, GZMK, and CXCR6 (online supplemental figure S8J). In HBV-positive tumors, macrophages were characterized by downregulated inflammatory responses, whereas NK cells upregulated cytotoxicity and allograft rejection responses (figure 6K and online supplemental figure 8H–J). Collectively, HBV infection in iCCA may alleviate myeloid inflammation and promote antitumor immunity. Hence, HBV infection was slightly associated with longer overall survival in our cohort (p=0.091, log-rank test; figure 6L), consistent with a previous study.45

Discussion

Tremendous sequencing studies have improved our understanding of the molecular pathogenesis of iCCA.8–11 33 However, little is known about the immunogenomic profile of iCCA and how to leverage this information to maximize responses to immune therapies. Herein, we provided an immunogenomic characterization of iCCAs, involving WES, RNA-seq, proteomics, scRNA-seq and mouse model analyses that could enable the discovery of novel biological insights and identification of promising targets for precise oncological practice.

We identified three distinct immune subgroups with diverse clinical, immune, genetic, and molecular features. IG1 was defined by the enrichment of neutrophils and iDCs, as well as KRAS mutations. KRAS mutation is associated with the recruitment of myeloid cells and MDSCs into the TME, which then drives immune suppression.46 Consistently, KRAS mutation may probably be the primary cause of myeloid inflammatory and immunosuppression in IG1. IG2 was an immune-exclusion subgroup with high tumor purity, strong proliferative activity, and defects in APM. Some tumors in IG2 could potentially have an ancestral immune-active microenvironment, and then gradually became “exclusion” during the enrichment of HLA LOH (figure 4G). IG3 was dominated by adaptive immune cells, aDCs, and NK cells, demonstrating an immune-activated phenotype with the enrichment of TLSs. The immunogenomic features of the three immune subgroups were reproducible in multiple independent cohorts,12 33 supporting that our classification was suitable and reliable.

Chronic inflammation induced by HBV infection increases the risk of HCC and iCCA. Using high-dimensional CyTOF analysis, Lim et al, revealed that HBV-related HCC microenvironment is more immunosuppressive and exhausted than that of non-viral-related HCC.47 In contrast, we found that HBV infections alleviated myeloid inflammation and enhanced the cytolytic activity of NK cells in iCCA, indicating that this common etiology may lead to distinct microenvironments between HCC and iCCA. Recently, three randomized phase III clinical trials in more than 1600 patients with advanced HCC revealed that immunotherapy was superior in patients with viral-related HCC.39 However, it is still unknown whether underlying viral hepatitis impacts on ICIs’ response in iCCA. We here proposed that HBV infection in iCCA could counteract the excessive myeloid inflammation and potentially reinvigorate antitumor immunity, probably favoring the efficacy of ICIs. The reduced myeloid inflammation in patients with virus-related iCCA might be due to the selection for KRAS wild type tumors during persistent HBV infection.48 Another rational explanation is the negative correlation between chronic HBV infection and innate immune response, especially the function of neutrophils.49

Therapeutically, this study aimed to characterize iCCA immune subgroups amenable to specific individualized therapies (figure 7 and online supplemental table S7). IG1 might be responsive to myeloid targeted therapies, such as C-X-C motif chemokine receptor 2 (CXCR2) and colony stimulating factor receptor (CSFR) inhibitors to deplete or reprogram tumor-associated neutrophils and M2 macrophages, respectively. Recent success in selective inhibitors of KRAS G12C50 may spark a renewed enthusiasm in targeting KRAS mutations in this immune-suppressive subgroup. IG2 may acquire better outcomes from specific CDK inhibitors, as well as strategies that regulate HLA and APM transcription. The enrichment of TLSs in IG3 and its potential role in the response to ICIs treatment may help guide clinical application of ICIs. Notably, the three immune subgroups showed distinct expression patterns of immune checkpoints, underscoring the need for personalized ICIs.

Figure 7

Summary of immunogenomic features of intrahepatic cholangiocarcinoma. Radar charts represented mean Z scores for the molecular features and clinical characteristics describing immune-suppressive, immune-exclusion, and immune-activated subgroups. Therapeutic targets of each subgroup were proposed.

Our study had some limitations. The use of a retrospective cohort with varied clinicopathological features and treatments may require further prospective validation. In addition, immunostaining on TMAs for estimating immune cell infiltration was likely to over-represent or under-represent the density or distribution of each immune subset. In conclusion, our integrative multimodule analyses provide a comprehensive understanding of iCCA’s immune landscape. The three immune subgroups with specific immune and non-immune profiles may have implications for the design of clinical trials using novel therapeutic strategies in iCCA.