Novel Molecular Subgroups of Intrahepatic Cholangiocarcinoma Discovered by Single-Cell RNA Sequencing-Assisted Multi-Omics Analysis

doi:10.21203/rs.3.rs-729739/v1

Download PDF

Article

Novel Molecular Subgroups of Intrahepatic Cholangiocarcinoma Discovered by Single-Cell RNA Sequencing-Assisted Multi-Omics Analysis

https://doi.org/10.21203/rs.3.rs-729739/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Intrahepatic cholangiocarcinoma (ICC) is a relatively rare but highly aggressive tumour type that responds poorly to chemotherapy and immunotherapy. Comprehensive molecular characterization of ICC is essential for the development of novel therapeutics. We performed a comprehensive multi-omics analysis of ICC via proteomic, whole-exon sequencing (WES) and single-cell RNA sequencing (scRNA-seq). We identified three molecular subtypes with deteriorating prognosis in ICC: chromatin remodelling, metabolism, and inflammation. The inflammation subtype was associated with a poor prognosis. Our random forest algorithm revealed that the mutation of KMT2D frequently occurred in the metabolism subtype and were associated with lower inflammatory activity. scRNA-seq further identified a novel APOE⁺C1QB⁺ macrophage subtype, which showed the capacity of promoting the inflammation subtype and contributing to a poor prognosis in ICC. Taking together, with a single-cell transcriptome-assisted multi-omics analysis, we identified novel molecular subtypes of ICC and valided APOE⁺C1QB⁺ TAMs as potential novel immunotherapy targets against ICC.

Cancer Biology

Comprehensive Molecular Characterization

Whole-exon Sequencing

Inflammation Molecular Subtype

Intrahepatic cholangiocarcinoma (ICC) is a relatively rare but highly aggressive form of primary liver tumour following hepatocellular carcinoma (HCC) and accounts for 5–10% of all primary liver malignancies ¹. 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 ¹. Surgical resection is the effective approach for ICC. Nevertheless, ICC patients are oftentimes diagnosed at a late stage when resection is not possible and only palliative chemotherapy is applicable ².

The inflammatory environment near the biliary tree may facilitate the development and progression of ICC. Specifically, inflammation in the tumour microenvironment (TME) is promoted by the secretion of various cytokines and chemokines, further aiding in tumour progression and distant metastasis formation. TNFα, IL-1β, IL-6, TGFβ and other cytokines have been proven to take part in a multi-step process of ICC transformation ³^,⁴^,⁵^,⁶. The detailed comprehensive molecular features of the ICC TME are, however, still not clear, especially at the single-cell level. Furthermore, the heterogeneous prognoses of ICC patients require updated subtypes classification for effective personalized treatment.

In this study, we performed a multi-omics analysis on a large series of ICC patient cohorts. Three molecular subtypes, “inflammation”, “metabolism”, and “chromatin remodelling”, were identified through the proteomics data. The “inflammation” subtype was found to be associated with the poorest prognosis in ICC. A random forest algorithm identified the key featured gene mutations within each molecular subtype using whole-exon sequencing (WES). Single-cell RNA sequencing (scRNA-seq) was used to generate a detailed characterization of the TME of the three molecular subtypes of ICC. Through this analysis, a novel APOE⁺C1QB⁺ tumour-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 novel targets of immunotherapy against ICC.

Three molecular subtypes of ICC discovered by proteomic analysis

To comprehensively characterize the molecular features of the ICC tumours, ICC tissues derived from 110 patients (Supplementary Table 1) were procured for proteomics analysis (Figure 1A). A total of 10,888 proteins were identified and quantified with low batch variance (Figure. S1A-S1B).

Three molecular subtypes of ICC were identified by consensus clustering based on the proteomic data (Figure 1B, Figure S1C-S1D): each of these clusters was distinguishable in the principal component analysis (PCA) plot of the patients, including all the detected proteins (Figure S1E). A Differentially expressed gene (DEG) analysis was performed for each pair of molecular subtypes (Figure S1F-S1H). A Venn plot was applied to obtain the key DEGs for the three molecular subtypes (Figure 1C). A total of Eighty-five genes were identified to be differentially expressed in all comparisons between the three molecular subtypes (Figure 1C and Table S2). Based on a GO analysis of these 85 genes (Figure S1I-S1K), we annotated the three molecular subtypes as “inflammation”, “metabolism”, and “chromatin remodelling” subtypes (Figure 1B-1C). A Gene Set Enrichment Analysis (GSEA) using the hallmark gene set ⁷ showed that the inflammatory response score was greatest in the inflammation subtype compared with the other two (Figure 1D). Similarly, the metabolism score was greatest in the metabolism subtype, and the chromatin remodelling score was greatest in the chromatin remodelling subtype (Figure 1E-1F). More details of the DEGs and activated pathways are shown in Figure. 1G. The three subtypes of ICC patients had significantly different survival rates (Figure 1H). Specifically, patients belonging to the metabolism subtype exhibited the shortest survival, while the chromatin remodelling ones survived the longest. Through proteomic analysis, we thus stratified ICC into groups with different molecular characteristics and survival rates.

KMT2D^mutis the key regulator in the “Inflammation” and “Metabolism” ICC subtypes

To further explore the genomic characteristics of the three molecular subtypes, WES was performed, and the genomic landscape was depicted on 78 samples (Figure 2A). ARID1A, KRAS, PBRM1, BRCA2, CREBBP, TP53, KMT2D, ATM, and ATRX were the most frequently mutated genes in the ICC tumour tissues (Figure 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 (Figure S2A). As expected, RTK-RAS was the pathway that was the most significantly altered pathway in ICC ⁸ (Figure S2B). KRAS and BRAF were the most frequently mutated genes in the RTK-RAS pathway (Figure S2C). Nevertheless, the Kaplan-Meier curve and univariate Cox regression indicated an insignificant association between overall survival (OS) and KRAS mutation or RTK-RAS pathway alteration (Figure S2D-S2E).

We used a random forest-based machine learning method to identify the key gene mutations in each molecular subtype. The feature importance ranking is shown in Figure. 2C-2E. KMT2D mutation was found in both inflammation and metabolism subtypes. Among all identified mutations, KMT2D’s was the only one to be significantly associated with a better prognosis for ICC patients (Figure 2F). A significantly high proportion of KMT2D mutation was also found in the metabolism subtype (Figure 2G). KMT2D^mut tissue had a greater metabolism score and lower inflammatory score than KMT2D^wt tissues (Figure 2H-2I). In addition, we found that most KMT2D^mut had a missense mutation in exon 39 (Table S3). Overall, we identified KMT2D^mutas the key regulator in the inflammation and metabolism ICC subtypes.

KMT2D^mutis associated with less TAMs in the ICC TME

Next, we dissected the immune landscape of different subtypes of ICC. The abundance of 21 immune cell populations was estimated by the ssGSEA method by assessing the gene signature for each cell type (Figure 3A). A significantly greater proportion of macrophages was found in the inflammation subtype (Figure 3B). We also observed that more macrophages were found in the KMT2D^WTpatients than in the KMT2D^mutpatients (Figure 3C).

To further characterize the difference in the TME between KMT2D^WTand KMT2D^mut patients, we performed mutation panel analysis for 10 patients using punctured specimens (Figure 3D and 3E), and then patients were selected based on the mutation profile (Figure 3E). Specimens from 2 KMT2D^WTand 1 KMT2D^mut patients undergoing surgery were used for scRNA-seq. A total of ten cell types, including malignant cells, fibroblasts, endothelial cells and various immune cells, were identified in the ICC TME (Figure 3F and Figure S3). The featured genes in each cell type were then determined (Figure S3). For instance, SPP1, S100A6, FGG, GDF15, and FGB were highly expressed in malignant cells (Figure S3). Again, fewer TAMs were found in the KMT2D^mutICC TME (Figure 3G). Thus, based on both the proteomic and the scRNA-seq analysis of the ICC TME, we found that the KMT2D^mutICC TME was associated with fewer TAMs in the ICC TME.

APOE⁺C1QB⁺TAMs are correlated with a poor prognosis and contribute to the inflammation subtype of ICC

To further understand the molecular mechanism underlying the ICC molecular subtypes, we next focused on the TAMs in the ICC TME as TAM was found to be a major difference between the three subtypes.

The macrophages were further split into two clusters: type 1 and type 2 (Figure 4A). The scores for each hallmark gene set (http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp) was calculated, and unsupervised clustering was performed. The results indicated that the two types of TAMs can be clearly separated by an unsupervised clustering on the score of hallmark gene sets, highlighting the distinct molecular mechanisms of the two types of TAMs in ICC (Figure 4B).

Several M2 macrophage markers (IL10, IL1B, HLA-DRA, and CD163) were expressed in both type 1 and type 2 clusters (Figure S4A). The signature scores of M2 and M1 were both greater in type 2 rather than type 1 macrophages (Figure S4B-S4C). The traditional M1/M2 classification for macrophages in ICC may thus not be suitable.

We next further characterized the two types of macrophages by several steps. The inflammatory response score was calculated for the two subtypes of TAMs. Type 2 TAMs showed a significantly greater inflammatory response score than type 1 TAMs (Figure 4C-4D). DEG analysis was then performed between the two subgroups of TAMs. The results indicated that APOE, LIPA, CTSD, C1QB and several other genes were significantly upregulated in type 2 TAMs (Figure 4E). We found a positive regulation of the inflammatory response, and a negative regulator of the immune system in type 2 TAMs, suggesting the potentially important role of type 2 macrophages in regulating the inflammatory response in ICC (Figure 4F). The high expression of APOE and C1QB in type 2 TAMs was verified by the co-staining of APOE and CD68, and C1QB and CD68 (Figure 4G-4I). Furthermore, a higher proportion of APOE⁺C1QB⁺ TAMs was found in the inflammation subtype than in the other two (Figure. 4J). The Kaplan-Meier curve and the univariate Cox regression analysis suggested that patients with high expression of APOE and C1QB had significantly worse OS than those in the low-risk group (Figure S3D-S3E).

With the top upregulated type 2 TAM DEGs from the comparison between type 1 and type 2 TAM as a type 2 APOE⁺C1QB⁺ TAM signature (Table S4), we performed ssGSEA based on the proteomic data. Our results indicated that the APOE⁺ macrophage score was significantly correlated with the inflammatory response score, suggesting the potential involvement of APOE⁺C1QB⁺ TAMs in the inflammatory activity of ICC tissues (Figure 4K). The inflammation subtype exhibited the highest APOE⁺ macrophage score (Figure 4L), suggesting that type 2 APOE⁺C1QB⁺ TAMs are the main determinant of the inflammation subtype of ICC. Consistently with previous findings, the type 2 APOE⁺C1QB⁺ TAM score was higher in KMT2D^WT patients than in KMT2D^Mut patients (Figure 4M). Furthermore, patients with a high type 2 APOE⁺C1QB⁺ TAM signature score exhibited significantly worse OS than those with a low one (Figure S4G).

KMT2D^mut reduces the number of TNF-α-producing CD4⁺ T cells induced by APOE⁺C1QB⁺ TAMs

We next ask how KMT2D mutation shaped TAMs in the ICC TME. First, we screened a series of human ICC cell lines and found that KMT2D was highly expressed by CCLP1 while it had a low expression by HCCC9810 (Figure 5A). We thus chose both two cell lines and co-cultured them with M0 macrophages induced by PMA (Figure 5B) ⁹. M0 macrophages co-cultured with ICC cells (HCCC9810 or CCLP1) resulted in a higher expression of APOE and C1QB in co-cultivated M0 cells than M0 cells cultured alone, which confirmed the induction of APOE⁺C1QB⁺ TAMs by ICC malignant cells (Figure 5C-5E). Furthermore, M0 macrophages co-cultured with CCLP1 showed significantly higher APOE and C1QB expression than M0 macrophages co-cultured with HCCC9810 (Figure 5C-E).

Next, we investigated how APOE⁺C1QB⁺ TAMs promoted inflammation in the ICC TME. With the APOE⁺C1QB⁺ TAM and CD4⁺ T cell abundance scores calculated by ssGSEA using their gene signatures, patients with high APOE⁺C1QB⁺ TAM abundance were found to exhibit higher CD4⁺ T cell infiltration (Figure 5F) and many CD4⁺ T cells in ICC TME produced TNF-α (Figure 5G).

To investigate whether APOE⁺C1QB⁺ TAMs regulate the inflammatory activity of T cells, we performed a NicheNET analysis on APOE⁺C1QB⁺ TAMs and T cells, modelling their intercellular communication by linking ligands of sender cells to target genes of receiver cells ¹⁰. 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 (Figure 5H). 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. As 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 terms of TNF by a co-culture experiment (Figure 5I). KMT2D^high CCLP1 cells and KMT2D^low HCCC9810 cells were co-cultured with M0 macrophages and T cells (Figure 5I). FACS was performed to elucidate the effect of APOE⁺C1QB⁺ macrophages on T cells. Our results showed the upregulation of TNF-α in T cells by co-culturing them with APOE⁺ C1QB⁺ macrophages (induction by the ICC cell lines) (Figure 5J). The co-culture of CCLP1, M0 macrophages, and T cells resulted in the highest TNF-α expression in T cells (Figure 5K).

Targeting APOE+C1QB⁺TAMs as a novel immunotherapy against ICC

Since APOE⁺ C1QB⁺TAM was the most important player driving the “Inflammation” subtype of ICC with the worst prognosis, we tried to target the APOE⁺ C1QB⁺TAM. A mouse ICC tumour model was constructed, followed by CSF1R antibody injection to perform TAM depletion (Figure 6A). Significantly less tumour burden was observed in the CSF1R antibody-treated group (Figure 6B-6C). HE staining was used to assess the properties of the ICC tumour model (Figure6D). Ki67 IHC staining suggested lower proliferation properties was found in CSF1R antibody-treated group than in the control group (Figure 6E). IF staining revealed a lower number of APOE⁺C1QB⁺ TAM in the CSF1R antibody-treated group than in the control (Figure 6F). The proportion of APOE⁺C1QB⁺ TAMs in all F4/80 positive cells was quantified. A lower proportion of APOE⁺C1QB⁺ TAM was found in the CSF1R antibody-treated group than that in the control (Figure 6G-6H). FACS analysis was performed to check TNF-α production by T cells. After the CSF1R antibody treatment, APOE⁺C1QB⁺ TAM was depleted and a lower T cell inflammation was observed as quantified through the TNF signal (Figure 6I-6J). These results indicate the potential role of APOE⁺C1QB⁺ TAM in driving the ICC “Inflammation” subtype.

In this study, we provided a comprehensive molecular characterization of ICC through system-level multi-omics analyses. ICC is a relatively rare but highly aggressive cancer type that has a low 5-year survival rate. The clinical feature and survival outcome of ICC varies among cohorts ¹¹^,¹², , thus requiring a novel molecular classification-based prognosis model. 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 identified three novel molecular subtypes. The potential molecular mechanism was also explored.

Inflammation has been long related to ICC development and metastasis formation through the modulation of the components of the TME ¹³^,¹⁴. Nevertheless, the mechanism and effect of inflammation in ICC are still controversial. Daniela Sia et al identified inflammation and proliferation classes as two distinct subtypes in ICC ¹². The inflammatory class, defined by the humoural response and imbalance of cytokines, was associated with a more favourable prognosis ¹². Nevertheless, several other studies provided contradictory data and opinions ¹¹^,¹⁵^,¹⁶. In our study, we confirmed that inflammation was the key regulator. Patients in the “inflammation” subtype exhibited the worst survival outcome compared with patients with the other two subtypes. These contradictory findings regarding the inflammatory activity in ICC lead us to further explore the detailed molecular mechanism of inflammation in ICC from 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 ¹⁷. Nevertheless, ARID1A, which only accounts for a 7% mutation rate in the TCGA-HCC cohort, was the most mutated gene in ICC. ARID1A performs context-dependent oncogenic and tumour-suppressing functions in liver cancer ¹⁸. The detailed mechanism of ARID1A mutation still requires further exploration in ICC. The detailed mechanism of ARID1A mutation still requires further exploration in ICC. KMT2D was one of the most common mutations in pan-cancer, and always correlated with a poor prognosis ¹⁹. Nevertheless, the KMT2D mutation indicated a better prognosis than wild-type KMT2D in our ICC cohort. Our random forest algorithm found that KMT2D was associated with a lower inflammatory activity and was associated with the “metabolism” subtype. Our conclusion is supported by an independent study based on Cas9 mutagenic screening ²⁰.

Through the ssGSEA analysis of the proteomic data, we found that macrophages were mostly associated with the inflammation subtype of ICC. We thus focused on macrophage cell populations in the single-cell transcriptome level. Interestingly, the mutation of KMT2D was found to be associated with low macrophage abundance at the bulk and single-cell levels. The detailed molecular mechanism of the relationship of KMT2D mutation with macrophages is worthy of further exploration.

M1 and M2 macrophages have been identified in various tumour types with distinct properties ²¹^,²². M1 macrophages are highly active in inflammatory responses and are commonly considered anti-tumour macrophages, while M2 macrophages can promote tumour cell growth and evasion by dysregulating several key pathways in the tumour cells ²². Unlike the traditional M1 and M2 classification, in our study, the ICC tumour macrophages were divided into two subtypes with distinct molecular features. Type 2 TAMs showed both higher M1 and M2 scores than the type 1 TAMs, suggesting a new classification method should be used in ICC macrophages. Type 2 TAM was 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 TME. Importantly, a high abundance of APOE⁺C1QB⁺ TAMs was found in the “Inflammation” subtype, which once again implied the importance of APOE⁺C1QB⁺ TAMs in the inflammatory response of ICC.

To further explore the mechanism of type 2 APOE⁺C1QB⁺ TAM in activating the inflammatory activity in the TME of HCC, we explored the cell-cell interaction (CCI) between APOE⁺C1QB⁺ TAMs and T cells by single-cell transcriptome. CCI analysis revealed that APOE⁺C1QB⁺ TAMs promoted TNF expression in T cells through several ligands. To confirm the results of the bioinformatics analysis, we co-cultured M0 macrophages and ICC tumour cell lines. The expression of APOEandC1QB was found to be higher in co-cultured macrophages than in macrophages alone, which confirmed the induction of APOE⁺C1QB⁺ TAMs by ICC malignant cells. FACS revealed APOE⁺C1QB⁺ TAMs promoted T cells inflammation by increasing the TNF signal. Our in vivo tumour model revealed the depletion of APOE⁺C1QB⁺ TAMs that decreased T cells inflammation. Reducing the abundance of type 2 APOE⁺C1QB⁺ TAM would, therefore, be a promising way to hinder the progression of HCC and improve the immunotherapy response.

There are still some limitations to our study. Although we illustrated the immune cell atlas of ICC by scRNA-seq, the chosen markers may limit the identification of some unknown immune cell phenotypes. Further, we recruited 110 cases of fresh tumor tissues, but we were not able to perform WES in all of the samples, which may have reduced the reliability of our conclusion. Additionally, we proposed our molecular classification of ICC contained information for clinical decision making, and provided our clinical practice to further prove this system. However, the cohort size was still relatively small.

In conclusion, we introduced novel molecular subtypes for stratifying ICC patients with different prognoses. In particular, the “Inflammation” subtype was associated with the poorest outcome. The mutation of KMT2D was found to be significantly involved in the inflammatory response and metabolism activity of ICC. A novel APOE⁺C1QB⁺ TAMs subtype showed the capacity to affect T cells, therefore reshaping the inflammation of TME in the “inflammatory response” subtype of ICC.

Sample collection

This study was approved by the Ethics Committee of Nanfang Hospital (NFEC-2019-057) and the First Affiliated Hospital (2021-Ani-58), Zhejiang University School of Medicine. All participants provided written informed consent to take part in the study. The detailed patient information is available in Table S1. Sample tissues were obtained from the biopsies of 110 patients diagnosed with ICC. All samples were assessed by pathological examinations.

Tandem mass tag (TMT) -based proteomic analysis

A total of Formalin-fixed paraffin-embedded (FFPE) samples (0.5-1 mg) were dewaxed and rehydrated and then subjected to acidic hydrolysis with formic acid (FA). Proteins were denatured with 6 M urea (Sigma-Aldrich, Germany) and 2 M thiourea (Sigma-Aldrich, Germany), and then digested into peptides with trypsin (1:20, Hualishi, Beijing, China) and Lys-C (1:80, Hualishi, Beijing, China) with the assistance of pressure-cycling technology (PCT) ²³^,²⁴. Peptides were labelled with TMTpro^TM 16 plex (Thermo Fisher Scientific™, San Jose, USA) ²⁵. 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 with a Thermo Dionex Ultimate 3000 RSLCnano System 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 and then analysed with a Q Exactive HF mass spectrometer using the data-dependent acquisition (DDA) mode (Thermo Fisher Scientific™, San Jose, USA). The database searching included all reviewed human entries from UniProt (downloaded on 14 April 2020, containing 20,365 proteins) using Proteome Discoverer (version 2.4, Thermo Fisher Scientific™, Waltham, MA). The detailed parameters were described previously ²⁶^,²⁷.

Molecular subtype identification

The “ConsensusClusterPlus” package was applied to the proteomic data to identify potential molecular subtypes of ICC ²⁸. Differentially expressed protein (DEP) analysis was performed with the “Limma” package to assess the proteomic data ²⁹. 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.

Single-sample GSEA (ssGSEA) and gene ontology (GO) analysis

ssGSEA was applied to evaluate the enrichment scores of each sample. The “GSVA” package was used for this analysis ³⁰. The hallmark gene sets were downloaded from The Broad Institute. Gene sets of immune cells populations were came from Bindea G et al. ³¹. The correlation between the risk and the enrichment scores was performed using Spearman’s coefficient. GO analysis was performed with the “clusterProfiler” package ³².

DNA extraction, library preparation, and targeted enrichment

Next-generation sequencing (NGS) analyses were performed according to protocols reviewed and approved by the Ethics Committee of the First Affiliated Hospital of Zhejiang University and Nanfang Hospital. DNA extraction, library preparation, and target capture enrichment were performed using previously described protocols with minor modifications ³³. Briefly, genomic DNA from white blood cells was extracted using the DNeasy Blood & Tissue Kit (Qiagen) and was applied as the normal control to remove germline variations. After deparaffinizing the FFPE samples with xylene, genomic DNA was extracted using the QIAamp DNA FFPE Tissue Kit (Qiagen). 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.

The KAPA Hyper Prep kit (KAPA Biosystems) was used to prepare libraries according to previously described protocols ³⁴. Briefly, 1-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. DNA libraries were generated by polymerase chain reaction (PCR) followed by purification using Agencourt AMPure XP beads.

WES was performed with a customized xGen lockdown probe panel (Integrated DNA Technologies). The blocking reagents were human cot-1 DNA (Life Technologies) and xGen Universal Blocking Oligos (Integrated DNA Technologies). Dynabeads M-270 (Life Technologies) and the xGen Lockdown Hybridization and Wash kit (Integrated DNA Technologies) were used to perform the capture reaction. 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.

Sequencing and bioinformatics analysis

The HiSeq4000 platform (Illumina) was used to sequence the target enriched libraries with 2×150-bp paired-end reads. bcl2fastq (v2.19) was employed to demultiplex the sequencing data. Trimmomatic ³⁵ 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) ³⁶, followed by further processing using the Picard suite (available at: https://broadinstitute.github.io/picard/) and the Genome Analysis Toolkit (GATK) ³⁷. Single nucleotide polymorphism (SNP) and insertion/deletion (indel) calling for GATK was performed with VarScan2 ³⁸ and HaplotypeCaller/UnifiedGenotyper, and the mutant allele frequency (MAF) cutoff was set as 0.5%. Common variants were removed according to dbSNP and the 1000 Genome project. The patients’ whole blood controls were used to filter out the germline mutations.

Gene fusions and copy number variations (CNVs) were identified by FACTERA ³⁹ and ADTEx ⁴⁰, respectively. For the 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. TMB 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 ⁴¹. 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 tumour suppressors.

Random forest algorithm for mutation importance ranking

A random forest algorithm was applied to the WES data to identify the most important mutations associated with the inflammatory response in ICC. Briefly, the gene mutation dataset and inflammatory response were applied as input. The “ranger” package was used to find the best hyperparameter in the regression process for the random forest model ⁴². The “e1071” package was then used to generate the random forest model ⁴³.

Tissue dissociation

RPMI 1640 (Gibco, Cat. no. 11875–093, US), with 1 mM protease inhibitor (Solarbio, Cat. no. P6730, China), was used to transport ICC tissues. Tissues were digested with a dissociation enzyme cocktail prepared by dissolving 2 mg/mL Dispase II (Sigma-Aldrich, Cat.42613-33-2 US), 1 mg/mL Type VIII Collagenase (Sigma-Aldrich, Cat. no. C2139, US), and 1 unit/mL DNase I (NEB, Cat. no. M0303S, US) in PBS with 5% foetal bovine serum (FBS; Gibco, Cat. no. 16000–044, US) for 40 min at 37°C. The cells were dissociated and collected every 20 min, and then filtered using a 40 μm nylon cell strainer (Falcon, Cat. no. 352340, US). Red blood cell lysis buffer (Invitrogen, Cat. US) with 1 unit/mL DNase I was used for removing red blood cells. Finally, cells were washed in PBS with 0.04% Bovine Serum Albumin (BSA; Sigma-Aldrich, Cat. no. B2064, US).

Library preparation and sequencing

The concentration of the single-cell suspension was computed with countess (Thermo) and adjusted to 1000 cells/μL. Cells were loaded according to the Chromium single-cell 3′ kit standard protocol to capture 5,000-10,000 cells/chip position (V2 chemistry). The library construction and all the other processes were performed according to the standard manufacturer’s protocol.

scRNA-seq data processing

Illumina HiSeq X Ten was used to obtain single-cell libraries using 150 nt paired-end sequencing. Reads were processed with the Cell Ranger 2.1.0 pipeline using default parameters. The STAR software was used to align FASTQs to the human reference genome (hg19). For each sample, gene-barcode matrices were generated by counting the unique molecular identifiers (UMIs) and filtering non-cell-associated barcodes. Finally, a gene-barcode matrix containing the barcoded cells and the gene expression counts was generated.

scRNA-seq data analysis

Single-cell transcriptome analysis was performed with “Seurat” with the default parameters ⁴⁴. A nonlinear dimensional reduction was performed with the t-distributed stochastic neighbour embedding (t-SNE) method. The gene features of each cluster were found by “Seurat”. The “GSVA” package was used to compute the score of each gene set as previously described. The “nichenetr” package was adopted to analyze the cell-cell interaction (CCI) by exploring the ligand and target gene pairs ¹⁰. The regulatory potential score was calculated using network propagation methods.

Cell lines and culture

The following cell lines were purchased from American Type Culture Collection (ATCC) and maintained in a humidified atmosphere in the CO2 incubator (Thermo, 3111) at 37 °C: ICC HCCC9810, CCLP-1, HuCCT1, RBE, KMCH1, LICCF, intrahepatic bile duct epithelial HIBEC, promonocytic THP1, and Jurkat. HCCC9810, THP1, HuCCT1, and Jurkat were maintained with RPMI-1640 media, while HIBEC, CCLP-1, LICCF, and KMCH1 were cultured with DMEM (Dulbecco’s Modified Eagle’s Medium) media, supplemented with both fetal bovine serum (10%, Gibco, USA) and penicillin/streptomycin (1%, Gibco, USA).

IHC staining for tissue sections

CCA tumor tissues were fixed with 4% paraformaldehyde before embedding them in paraffin. Tissue sections (4 μm) were later deparaffinized with xylene and rehydrated with graded alcohol. Antigen epitope retrieval was induced by microwave heating. Tissue sections were then blocked for 1h, at room temperature, in goat serum blocking solution (Proteintech, B900780, China) and phosphate buffered saline (PBS). Before primary antibody overnight incubation at 4°C, the next day ICCA tissue samples were stained using an anti-mouse/rabbit universal immunohistochemical detection kit (pk10006, proteintech) following the manufacturer’s instructions. Mounted sections were examined by light microscopy (Leica), and images were analyzed with Image-Pro Plus (version 6.0).

Immunofluorescence (IF) staining for tissue sections

Tissue samples were deparaffinized with xylene and rehydrated with graded alcohol. Antigen epitope retrieval was induced by microwave heating. Tissue sections were permeabilized with 0.2% Triton X-100 and then blocked with 5% BSA, for 1h, at room temperature, and then incubated with the primary antibody, overnight, at 4°C. The following day, tissue sections were washed with PBS three times. Samples were then stained with fluorescent secondary antibodies (Invitrogen, USA) in the dark, for 1h, at room temperature, and then washed with PBS three times before mounting them onto slides using Mowoil supplemented with DPAI to stain the nucleus. Finally, a confocal microscope (Nikon, Japan) was used to take images of ICCA tissue samples. Acquired images were analyzed using ImageJ (version 4.0).

IF staining for cell culture

To induce THP1 cells to differentiate into M0 macrophages,1 X 105 THP-1 cells were treated with 80 nM Phorbol 12-myristate 13-acetate (PMA) once plated into 12 well plates, after 48h, 2 x 105 ICCA tumor cells were plated for co-cultured with M0 macrophages on glass coverslips in 12 well plates. Immunofluorescence assays were performed as follows. First, cells were washed with PBS and fixed with 4% paraformaldehyde, before permeabilizing them with 0.2% Triton X-100 at room temperature. Samples were next blocked with 2% BSA for 1h, at room temperature, and then incubated with the primary antibody for 1h at room temperature. To examine the co-localization of proteins of interest (CD68, APOE, and C1QB), cells were stained with Alexa Fluor ® 488 and 594-labled fluorescent secondary antibody (Invitrogen, USA) for 1h in the dark at room temperature. Coverslips were washed three times with PBS before mounting them onto slides using Mowoil supplemented with DPAI to stain the nucleus. Finally, cells were examined with a confocal microscope (Nikon, Japan). The acquired images were analyzed using ImageJ (version 4.0).

In vivo tumour model

C57BL/6J mice were purchased from the Model Animal Research Center of Nanjing University (China). All mice were housed in the specific pathogen free (SPF) facility of the First Affiliated Hospital, Zhejiang University School of Medicine, with approval from the Institutional Animal Care & Use Committee (IACUC). To establish a spontaneous ICC mouse model, a plasmid mixture carrying T3-EF1a-NICD, PT3-myr-AKT-HA, and SB100 was injected at high pressure into the tail vein of the mice. One week later, the mice were treated with 10 μg/g CSF1R antibody via intraperitoneal injection, and the same injections were performed every 3 days thereafter. On day 28, tumors were collected and further analyzed.

Hematoxylin and eosin (HE) staining

Mouse tumor tissue sections (4μm) were deparaffinized with xylene, and then rehydrated with graded alcohol. Tissue sections were washed for three times using PBS before stained nucleus with hematoxylin for 30 min at room temperature before washed with PBS for three times. Then the sections were exposed to the ammonia water 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 min at room temperature, then stained cytoplasm with eosin for 1h at room temperature. And then used graded alcohol directly rinsed after the eosin staining. Finally, the xylene was used to displace the anhydrous alcohol before mounted with slides. sections were examined by light microscopy (Leica), images were analyzed with Image-Pro Plus software (version 6.0).

Ki67 IHC staining

Prepared tumor tissue sections(4μm) deparaffinized with xylene and rehydrated with graded alcohol. Antigen epitope retrieval were induced by microwave heating before blocked in blocking solution (Proteintech, B900780, China) in phosphate buffered saline (PBS) for 1h at room temperature. Then Ki67 primary antibody (1:200, 550609, BD Biosciences) were incubated for overnight at 4°C, the following day, tissue sections were washed three times with PBS before stained using anti-mouse/rabbit universal immunohistochemical 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).

Fluorescence-activated cell scanning (FACS) analysis

THP1 cells were treated with 80 nM PMA for 48 h to induce THP1 to M0 macrophages and then co-cultured with ICC tumour cells. Before the cells were incubated with fluorochrome-labelled antibody, they were stimulated with Cell Stimulation Cocktail (plus protein transport inhibitors) (eBioscience, 00-4975-03, USA) at 37 °C for 6 h. For surface marker analysis, live cells were re-suspended in 1×PBS supplemented with 2% FBS and stained with anti-human CD3 (UCHT1, BioLegend, 300406) at 4 °C for 30 min. Subsequently, the cells were fixed and permeabilized with Fixation and Permeabilization Solution (BD Bioscience, America) for 30 min and then incubated with anti-human TNF-α (MAb11, Biolegend, 502909) at 4 °C for 30 min. Data were acquired with an ACEA NovoCyte flow cytometer and analysed with FlowJo software (version 10).

Mouse ICCA tissue dissociation assays were performed as previously described. Cells were resuspended with DMEM (Dulbecco’s Modified Eagle’s Medium) media supplemented with 2% FBS, stimulated with Cell Stimulation Cocktail (plus protein transport inhibitors) (eBioscience, 00-4975-03, USA) for 6 h at 37 °C, centrifugated at 500 xg for 5 min, and then incubated with 2.5 μg/ml Fc blocker (Biolegend, 15660) on ice for 15 min. Next, suspensions were incubated with fluorochrome-labeled antibodies (CD45, CD3, CD4, and CD8) at 4 °C for 30 min, then washed with 1×PBS supplemented with 2% FBS and centrifuged at 500 xg for 5 min. Cells were fixed and permeabilized with Fixation and Permeabilization Solution (BD Bioscience, America) at room temperature for 20 min in the dark, before incubating them with anti-mouse TNFα (MP6-XT22, Biolegend, 506305). For 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, Fortessa) and analysed with FlowJo software (version 10).

Statistical analysis

Survival analysis was performed with the “survival” package ⁴⁵. 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. Sample data were analyzed by a two-tailed Student's t-test to identify statistically significant differences between two groups; the Kruskal-Walis test was used to identify differences between three groups. The data were represented as means ± standard error of mean (SEM). A p-value < 0.05 indicated statistical significance.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China 81874173 (J.R.), 81472346 (P.Z.), 82074208 (P.Z.), 81972492 (T.G.), by the National Key R&D Program of China grant 2019YFA0803000 (J.S.), 2020YFE0202200 (T.G.), by the Young Investigator Research Program 588020-D01907 (M.J.), by the Natural Science Foundation of Zhejiang Province LY20H160033 (P.Z.).

Author contributions

X.B., Q.L., J.C., T.G., J.S., P.Z., J.R. designed the research, X.B., Q.L., J.C., D.C., C.Y. developed experimental methods, performed most experiments and analyzed data. X.M., Y.W., X.L., X.R., .F.C., .M.J., Z.Z., R.S., L.H., W.W., Y.G. and J.W. provided technical supports. W.F., V.T. and W.F. supervised this work. X.B., Q.L., J.C., D.C. and C.Y. wrote the manuscript. T.G., J.S., P.Z. and J.R. edited the manuscript.

Competing interests

TNG is the shareholder of Westlake Omics Inc. LLH is the employee of Westlake Omics Inc. All other authors declare no competing interests.

Materials and correspondence

Correspondence and requests for should be addressed to T.G., J.S., P.Z. and J.R.

1. Gupta A, Dixon E. Epidemiology and risk factors: intrahepatic cholangiocarcinoma. Hepatobiliary surgery and nutrition 6, 101 (2017).

2. Kelley RK, Bridgewater J, Gores GJ, Zhu AX. Systemic therapies for intrahepatic cholangiocarcinoma. Journal of hepatology 72, 353-363 (2020).

3. Balkwill F. TNF-α in promotion and progression of cancer. Cancer and Metastasis Reviews 25, 409-416 (2006).

4. Okabe H, et al. Identification of CXCL5/ENA‐78 as a factor involved in the interaction between cholangiocarcinoma cells and cancer‐associated fibroblasts. International journal of cancer 131, 2234-2241 (2012).

5. Yadav A, Kumar B, Datta J, Teknos TN, Kumar P. IL-6 promotes head and neck tumor metastasis by inducing epithelial–mesenchymal transition via the JAK-STAT3-SNAIL signaling pathway. Molecular Cancer Research 9, 1658-1667 (2011).

6. Sato Y, et al. Epithelial-mesenchymal transition induced by transforming growth factor-β1/Snail activation aggravates invasive growth of cholangiocarcinoma. The American journal of pathology 177, 141-152 (2010).

7. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database hallmark gene set collection. Cell systems 1, 417-425 (2015).

8. Boerner T, et al. Genetic Determinants of Outcome in Intrahepatic Cholangiocarcinoma. Hepatology, (2021).

9. Daigneault M, Preston JA, Marriott HM, Whyte MK, Dockrell DH. The identification of markers of macrophage differentiation in PMA-stimulated THP-1 cells and monocyte-derived macrophages. PloS one 5, e8668 (2010).

10. Browaeys R, Saelens W, Saeys Y. NicheNet: modeling intercellular communication by linking ligands to target genes. Nature methods 17, 159-162 (2020).

11. Sellers CM, Uhlig J, Ludwig JM, Stein SM, Kim HS. Inflammatory markers in intrahepatic cholangiocarcinoma: Effects of advanced liver disease. Cancer medicine 8, 5916-5929 (2019).

12. Sia D, et al. Integrative molecular analysis of intrahepatic cholangiocarcinoma reveals 2 classes that have different outcomes. Gastroenterology 144, 829-840 (2013).

13. Kumar M, Zhao X, Wang XW. Molecular carcinogenesis of hepatocellular carcinoma and intrahepatic cholangiocarcinoma: one step closer to personalized medicine? Cell & bioscience 1, 1-13 (2011).

14. Hasita H, et al. Significance of alternatively activated macrophages in patients with intrahepatic cholangiocarcinoma. Cancer science 101, 1913-1919 (2010).

15. Gomez D, Morris‐Stiff G, Toogood GJ, Lodge JPA, Prasad KR. Impact of systemic inflammation on outcome following resection for intrahepatic cholangiocarcinoma. Journal of surgical oncology 97, 513-518 (2008).

16. Tsilimigras DI, et al. The systemic immune-inflammation index predicts prognosis in intrahepatic cholangiocarcinoma: an international multi-institutional analysis. Hpb 22, 1667-1674 (2020).

17. Eichenmüller M, et al. The genomic landscape of hepatoblastoma and their progenies with HCC-like features. Journal of hepatology 61, 1312-1320 (2014).

18. Sun X, et al. Arid1a has context-dependent oncogenic and tumor suppressor functions in liver cancer. Cancer cell 32, 574-589. e576 (2017).

19. Martínez-Jiménez F, et al. A compendium of mutational cancer driver genes. Nature Reviews Cancer 20, 555-572 (2020).

20. Wang G, et al. CRISPR-GEMM Pooled Mutagenic Screening Identifies KMT2D as a Major Modulator of Immune Checkpoint Blockade. Cancer Discov 10, 1912-1933 (2020).

21. Bronte V, Murray PJ. Understanding local macrophage phenotypes in disease: modulating macrophage function to treat cancer. Nature medicine 21, 117-119 (2015).

22. Wang N, Liang H, Zen K. Molecular mechanisms that influence the macrophage M1–M2 polarization balance. Frontiers in immunology 5, 614 (2014).

23. Zhu Y, et al. High-throughput proteomic analysis of FFPE tissue samples facilitates tumor stratification. Mol Oncol 13, 2305-2328 (2019).

24. Gao H, et al. Accelerated Lysis and Proteolytic Digestion of Biopsy-Level Fresh-Frozen and FFPE Tissue Samples Using Pressure Cycling Technology. J Proteome Res 19, 1982-1990 (2020).

25. Li J, et al. TMTpro reagents: a set of isobaric labeling mass tags enables simultaneous proteome-wide measurements across 16 samples. Nat Methods 17, 399-404 (2020).

26. Shen B, et al. Proteomic and Metabolomic Characterization of COVID-19 Patient Sera. Cell 182, 59-72 e15 (2020).

27. Nie X, et al. Multi-organ proteomic landscape of COVID-19 autopsies. Cell 184, 775-791 e714 (2021).

28. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572-1573 (2010).

29. Smyth GK. Limma: linear models for microarray data. In: Bioinformatics and computational biology solutions using R and Bioconductor). Springer (2005).

30. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC bioinformatics 14, 7 (2013).

31. Bindea G, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 39, 782-795 (2013).

32. Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics: a journal of integrative biology 16, 284-287 (2012).

33. Tong L, et al. Tumor-derived DNA from pleural effusion supernatant as a promising alternative to tumor tissue in genomic profiling of advanced lung cancer. Theranostics 9, 5532 (2019).

34. Yang Z, et al. Investigating novel resistance mechanisms to third-generation EGFR tyrosine kinase inhibitor osimertinib in non–small cell lung cancer patients. Clinical Cancer Research 24, 3097-3107 (2018).

35. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114-2120 (2014).

36. Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:13033997, (2013).

37. DePristo MA, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature genetics 43, 491 (2011).

38. Koboldt DC, et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome research 22, 568-576 (2012).

39. Newman AM, et al. FACTERA: a practical method for the discovery of genomic rearrangements at breakpoint resolution. Bioinformatics 30, 3390-3393 (2014).

40. Amarasinghe KC, Li J, Halgamuge SK. CoNVEX: copy number variation estimation in exome sequencing data using HMM. In: BMC bioinformatics). Springer (2013).

41. Fang W, et al. Comprehensive genomic profiling identifies novel genetic predictors of response to anti–PD-(L) 1 therapies in non–small cell lung cancer. Clinical Cancer Research 25, 5015-5026 (2019).

42. Wright MN, Ziegler A. ranger: A fast implementation of random forests for high dimensional data in C++ and R. arXiv preprint arXiv:150804409, (2015).

43. Meyer D, et al. Package ‘e1071’. The R Journal, (2019).

44. Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome biology 19, 15 (2018).

45. Therneau TM, Lumley T. Package ‘survival’. Survival analysis Published on CRAN 2, 3 (2014).

There is NO Competing Interest.

SupplementaryMaterials.pdf
figure s1, figure s2, figure s3, figure s4, table s1, table s2, table s3, table s4

Download PDF

Version 1

posted

You are reading this latest preprint version

Novel Molecular Subgroups of Intrahepatic Cholangiocarcinoma Discovered by Single-Cell RNA Sequencing-Assisted Multi-Omics Analysis

Status:

Version 1

Abstract

Figures

Introduction

Results

Discussion

Methods

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1