- 1Liver Cancer Institute, Zhongshan Hospital, Fudan University and Key Laboratory of Carcinogenesis and Cancer Invasion, Ministry of Education, Shanghai, China
- 2Neurosurgery Department of Zhongshan Hospital, Fudan University, Shanghai, China
- 3Department of Hepatobiliary and Pancreatic Surgery, The Second Affiliated Hospital of Nanchang University, Nanchang, China
- 4Department of Thoracic Surgery, Peking University Shenzhen Hospital, Shenzhen, China
Background: Cholangiocarcinoma was a highly malignant liver cancer with poor prognosis, and immune infiltration status was considered an important factor in response to immunotherapy. In this investigation, we tried to locate immune infiltration related genes of cholangiocarcinoma through combination of bulk-sequencing and single-cell sequencing technology.
Methods: Single sample gene set enrichment analysis was used to annotate immune infiltration status in datasets of TCGA CHOL, GSE32225, and GSE26566. Differentially expressed genes between high- and low-infiltrated groups in TCGA dataset were yielded and further compressed in other two datasets through backward stepwise regression in R environment. Single-cell sequencing data of GSE138709 was loaded by Seurat software and was used to examined the expression of infiltration-related gene set. Pathway changes in malignant cell populations were analyzed through scTPA web tool.
Results: There were 43 genes differentially expressed between high- and low-immune infiltrated patients, and after further compression, PNOC and LAIR2 were significantly correlated with high immune infiltration status in cholangiocarcinoma. Through analysis of single-cell sequencing data, PNOC was mainly expressed by infiltrated B cells in tumor microenvironment, while LAIR2 was expressed by Treg cells and partial GZMB+ CD8 T cells, which were survival related and increased in tumor tissues. High B cell infiltration levels were related to better overall survival. Also, malignant cell populations demonstrated functionally different roles in tumor progression.
Conclusion: PNOC and LAIR2 were biomarkers for immune infiltration evaluation in cholangiocarcinoma. PNOC, expressed by B cells, could predict better survival of patients, while LAIR2 was a potential marker for exhaustive T cell populations, correlating with worse survival of patients.
Introduction
Cholangiocarcinoma (CCA) has long been deemed as a malignancy with poor prognosis in liver cancer. Patients conflicted by cholangiocarcinoma often are found in late stages, who were not candidates for surgery and seldom benefit from chemotherapy or comprehensive treatment (1, 2). Though blockade of programmed cell death receptor 1/programmed cell death receptor ligand 1 (PD1/PDL1) axis with mono-antibody, Pembrolizumab and Nivolumab, has shed light on partial patients, who showed high PDL1 expression in tumors, the overall treating efficacy in advanced CCA patients still needs further observation (3–6). Understanding tumor immune microenvironment (TIME) and infiltration status of CCA could better guide the clinical appliance of immunotherapy (7–9).
With the development of single-cell sequencing (scRNA-seq), investigators could further examine gene expression in individual cells and try to locate functional difference between different clinical phenotypes, especially in immune cells that have infiltrated the tumor (10–12). Characterization of CCA immune microenvironment is limited, so in this study to characterize immune cell components in TIME, we combined bulk sequencing data with scRNA-seq data, which could provide a better understanding of functional cell clusters related to disease severity. We found B cell infiltration levels in CCA TIME were related to patients’ overall survival (OS), and propronociceptin (PNOC), which was highly expressed by B cell populations in CCA, could be an independent indicator for better prognosis. Also, in CCA, leukocyte associated immunoglobulin like receptor 2 (LAIR2) was highly expressed by regulatory T cells (Tregs) and part of CD8+/GZMB+ T cells, which could be an indicator of exhaustive immune status in CCA patients. In addition, CCA cell sub-populations demonstrated heterogeneous metabolic and signal transduction activities, in which some CCA cells showed highly activated PD1/PDL1 axis signals, justifying the application of anti-PD1 combining therapy in CCA patients.
Methods
Datasets for Analysis and Derivation of Gene List
In this investigation, dataset of cholangiocarcinoma (CHOL) in database (n = 36) of the Cancer Genome Atlas (TCGA) (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) was used to analyze the differentially expressed genes between high- and low-immune infiltration groups (13). After searching Gene Expression Omnibus (GEO) database, data series with patient count over 100 were located, and the datasets with largest patient counts (GSE32225, n = 149; GSE26566, n = 104) were chosen for further immune infiltration classification (14, 15). Single cell sequencing data of five intrahepatic cholangiocarcinoma patients was procured from dataset of GSE138709 (16). The clinical information for patients’ cohorts were publicly accessible, which does not require additional endorsement from the local ethic committee. The immune meta gene list for 28 immune cell types were downloaded from TISBID database (http://cis.hku.hk/TISIDB/index.php) (17). Workflow of this Investigation was provided in Figure 1.
Calculation of Immune Infiltration Scores in Bulk Sequencing Samples and Analysis of Differentially Expressed Genes Between Groups
In our analysis, single sample gene set enrichment analysis (ssGSEA) for immune infiltration annotation was performed to calculate respective immune infiltration scores of 28 immune cell types, which includes cell types of activated CD4 T cell, activated CD8 T cell, activated dendritic cell, CD56 bright natural killer (NK) cell, central memory CD4 T cell, central memory CD8 T cell, NK cell, NK T cell, type 1 T helper cell, type 17 T helper cell, CD56 dim NK cell, immature dendritic cell, macrophage, myeloid derived suppressive cell (MDSC), neutrophil, plasmacytoid dendritic cell, regulatory T cell (Treg), type 2 T helper cell, activated B cell, eosinophil, gamma delta T cell, immature B cell, mast cell, memory B cell, monocyte and T follicular helper cell (18). Differentially expressed genes between groups were analyzed using edgeR package, contrasting high- with low-immune infiltrated patients (19, 20).
Gene Set Enrichment Analysis
Gene set enrichment analysis (GSEA) was used to demonstrate the altered pathways between patient groups in this study, using software of GSEA v4.1.0 (Broad Institute, Inc., Massachusetts Institute of Technology, and Regents of the University of California) (21). The annotation of changed pathways in this investigation was performed with hallmarks gene set (version: 7.2).
Gene Ontology and Pathway Enrichment
DAVID database was used for gene ontology (GO) analysis, which included biological process, cellular compartment, and molecular function (https://david.ncifcrf.gov/summary.jsp) (22). The protein domains of differentially expressed genes between groups were also analyzed and downloaded from the database for demonstration. REACTOME database was also linked for annotation of significantly changed pathways between groups (www.reactome.org) (23).
Single Cell Data Processing
For single cell sequencing analysis, raw data for GSE138709 were downloaded from portal website, and package of Seurat was used to process data in R (version: 4.0.3) with R studio (version: 1.3.1903) (24–26). The raw data GSE138709 were loaded with Seurat, and cells were filtered with the criteria of >20% mitochondria related genes or more than 6,000 genes expressed. A total of 32,627 cells were included for further analysis, and variable features of each sample were analyzed after normalization. Then we used Seurat function of FindIntegrateionAnchors to merged sample files with common anchors among variables (dims=1:20, k.filter=30) (26). Merged data of cells were clustered into 15 cell populations using function of FindClusters (resolution = 0.3). Respective reduction of cell clustering, including UMAP, TSNE, and PCA, were performed. For cell population annotation, we used the signatures chosen in the original publication (16). For NK and T cell cluster, signatures of CD7, FGFBP2, KLRF1, CD2, CD3D, and CD3E were chosen for annotation. For malignancy and cholangiocyte, signatures of EPCAM, KRT19, KRT7, FXYD2, TM4SF4, and ANXA4 were chosen. For monocytes, CD14 and CD11C were chosen for annotation. For B cell cluster, CD79A, MS4A1 were chosen. For endothelial cells, signatures of ENG and VWF were chosen for annotation. For hepatocytes, APOC3, FABP1, and APOA1 were chosen for annotation. And for fibroblasts, ACTA2 and COL1A2 were chosen for demonstration.
Analysis of Pathway Changes in Malignant Cholangiocarcinoma Cells
To compute and analyze pathway scores in malignant cholangiocarcinoma cells, we used scTPA, which is a web tool for single-cell analysis of activated pathways (http://sctpa.bio-data.cn:8080/index.html) (27, 28). The malignant cell expression matrix was extracted by sample origins in malignancy and cholangiocyte cluster, and then expression matrix was uploaded online. Analyzed results were downloaded for further analysis and demonstration.
Correlation Between Specific Genes and Immune Infiltration Scores
TIMER 2.0 web tool (http://timer.cistrome.org) was used for correlation of gene expression with immune cell infiltration scores, which included scores calculated by CIBERSORT and MCP-counter methods (29–31). Scores of TCGA CHOL sequencing data calculated by other infiltration estimating methods were also downloaded from website for analysis.
Correlation Between Specific Gene Markers
Database GEPIA (http://gepia.cancer-pku.cn) was used for correlation analysis between PNOC, LAIR2, and a series of immune regulators in bulk sequencing data of CHOL and hepatocellular carcinoma (LIHC) in TCGA database (32).
Survival Analysis of Genes in Outside Database
Survival analysis of specific genes was performed in outside database, KMplotter, which is an integrated portal for tumor survival analysis, combining genomics data of microarray with clinical information (33).
Statistics
Survival analysis in this investigation was performed with R packages survival and survminer in R environment, which were used to find the best cutoff values for survival comparison between groups (34). Package of pheatmap was used to construct heat maps (35). Dot plots for correlation analysis and bar plots for GO analysis were generated by packages ggplot2, using spearman correlation test, and GOplot respectively (36, 37). A generalized linear model (GLM) in R was used for prediction of immune infiltration status, using differentially expressed genes, and then stepwise algorithm (backward) was used to choose the appropriate model by an information criterion (AIC) extracted from formerly fitted model (AIC= −2*log L + k* edf; L: likelihood; edf: the equivalent degrees of freedom). Receiver operating characteristics (ROC) were examined using package plotROC. P value under 0.05 was considered significant.
Results
Cholangiocarcinoma Patients in TCGA Dataset Were Clustered Into High- and Low-Immune Infiltrated Groups With Different Prognosis
Using immune gene list for 28 immune infiltrating cell populations, we generated scores for each immune cell type. After clustering cholangiocarcinoma patients according to the calculated scores, we found there was a clearly different immune status between groups (Figure 2A). We also used gene lists for immune stimulators, inhibitors, MHC molecules, chemokine, and chemokine receptors to calculate the corresponding scores, and in high-immune infiltration patients, expression levels for those genes were much higher (Figure 2B). Patients with high immune infiltration showed better prognosis (Figure 2C).
Figure 2 Patients with cholangiocarcinoma were divided into differentially immune infiltrated groups with different prognosis. (A) Clustering of CCA patients according to immune infiltration status calculated by ssGSEA method. (B) Whole scores of chemokine, chemokine receptor, immune stimulator, immune inhibitor, and MHC expression levels between groups. (C) Survival difference between high- and low-immune infiltrated cholangiocarcinoma patients. (D, E) Differentially expressed genes between high- and low-immune infiltrated patients.
Differentially Expressed Genes Between High- and Low-Infiltrated Patients Were Mainly About Immune Functions, and Inflammatory Signals Were Highly Enriched in High-Immune Infiltrated Patients
Differentially expressed genes between high- and low-infiltrated groups were analyzed, and only a set of 43 genes were up-regulated in high-infiltrated patients (Figures 2D, E) Pathway enrichment showed the up-regulated gene set was mainly about inflammatory signals, immune stimulation, and PD1 axis (Figures 3A, B). Gene ontology enrichment for 43 gene set showed those genes were involved in the process of adaptive immune responses and T cell signaling (Figures 3C, D). The protein functional enrichment showed most of the 43 genes were immunoglobulins (Figures 3E, F). Among those genes that were significantly survival-related, all could indicate better overall survival with higher expression (Figures 3G–L). Further gene set enrichment analysis between groups showed hallmarks of complement signaling, IL2/STAT5 signaling, IL6/JAK/STAT3 signaling, inflammatory response signaling, interferon gamma signaling, and TNFA signaling via NFKB were highly enriched (Figures 3M–R).
Figure 3 Functional Enrichment of Differentially Expressed Genes Between High- and Low-Immune Infiltration Groups. (A, B) Pathway enrichment of differentially expressed genes in REACTOME database. (C, D) Gene ontology enrichment of differentially expressed genes. (E, F) Protein function enrichment of differentially expressed genes. (G–L) Among differentially expressed genes, PNOC, TRBC1, TRAV29DV5, IGLV3.16, and AC244205.1 were significantly correlated with CCA patients’ overall survival, while LAIR2 did not achieve significance. (M–R) Signatures of complement pathway, IL2-STAT5 pathway, IL6-Jak-STAT3 pathway, inflammatory response pathway, interferon-gamma response pathway, and TNF via NFKB pathway were highly enriched in high-immune infiltrated patients.
Several Genes Were Associated With Immune Infiltration Status by Stepwise Regression Model
We further calculated immune infiltration scores for datasets of GSE26566 and GSE32225, and after clustering patients into high- and low-infiltration groups, we used backward stepwise regression model to compress the 43 gene set in prediction of immune infiltration status in the two datasets respectively (Table 1). In both models (GSE26566: infiltration score = 6.846 − 0.053*SH2D1A – 0.061*PNOC – 0.021*LAIR2; GSE32225: infiltration score = −1.690 + 0.014*SH2D1A – 0.007*LAIR2 – 0.010*ICOS + 0.019*HEMGN + 0.012*GTSF1L), LAIR2 were related to high-immune infiltration status (Supplementary Figure 4).
Further Demonstration of CCA Tumor Microenvironment Showed PNOC Was Mainly Expressed by B Cell, Which Was Also an Indicator for Better Prognosis
In addition to bulk sequencing analysis, we analyzed the immune microenvironment of intrahepatic cholangiocarcinoma with single cell sequencing dataset GSE138709. We further clustered cell populations into 15 clusters, and using genes CD7, CD3D, KRT19, FXYD2, CD14, CD1C, CD79A, VWF, APOC3, and ACTA2, we classified 15 cell clusters into 7 cell populations, which were fibroblasts, NK and T cells, malignancy and cholangiocytes, endothelial cells, monocytes, hepatocytes, and B cells (Figures 4A–C). Proportions of cells in different tissue types showed most cells in malignancy and cholangiocyte cluster were from tumor samples, while most cells in NK&T cell cluster were form adjacent samples. Also, a high portion of fibroblasts were also seen in tumor samples (Figures 4D–F).
Figure 4 Single Cell Atlas of CCA Patients According to Dataset GSE138709. (A, B) Cell clusters for GSE138709 of five CCA patients. (C) Cell markers for clusters’ annotation. (D) Portions of adjacent and tumor tissues in different cell clusters. (E) Patient portions in different cell clusters. (F) Numbers for cell clusters in dataset after filtration.
We found PNOC was highly expressed in tumors, though the difference between normal and tumor tissues was not significant (Figure 5A). We further examined genes’ expression, selected by stepwise regression models, in single cell populations, and PNOC was mainly expressed by B cell cluster. After sub-clustering, activated B cells, plasma cells, and naive B cells all showed expression of PNOC (Figures 5B–D). We further examined the correlations between PNOC and scores for B cell infiltration in TCGA CHOL samples, calculated by ssGSEA and MCPcounter methods, and results showed PNOC was highly correlated with B cells (Figures 5E–H). Also, we found high B cell infiltration scores in cholangiocarcinoma were related to better prognosis, though survival benefits of high immature B cell and memory B cell scores were not significant due to small sample size (Figures 5I–L). We further used database GEPIA to examine B cell markers’ correlation with PNOC, and results showed PNOC was highly correlated with CD19, CD79A, CD27, and FCRL5 in bulk sequencing data of CHOL (Coefficients >0.95) (Figures 5M–P).
Figure 5 PNOC Was Highly Expressed by B Cell Populations in CCA, and B Cell Infiltration Levels in CCA Indicated Better Overall Survival. (A) PNOC was highly expressed in CCA tumors in TCGA database, though significant difference was not achieved. (B) PNOC was mainly expressed by B cells in single cell levels. (C) Further cluster of B cell populations. (D) Markers for sub B cell populations. (E–H) Correlation between PNOC expression and scores for activated CD8+ T cell, immature B cell, memory B cell, and whole B cell calculated by ssGSEA and MCPcounter methods. (I–L) B cell infiltration levels for activated CD8+ T cell, immature B cell, memory B cell, and whole B cell calculated by ssGSEA and MCPcounter methods were correlated with CCA patients’ overall survival. (M–P) Correlation between CD19, CD79A, FCRL5, CD27, and PNOC in CHOL bulk sequencing samples from GEPIA database.
LAIR2 Was Up-regulated in CCA Samples, Which Was Mainly Expressed by Regulatory T Cells and a Subset of CD8+/GZMB+ T Cells
We further divided NK and T cells populations into sub-clusters, and LAIR2 was found to be expressed by Foxp3+ regulatory T cell and CD8+/GZMB+ T cell clusters (Figures 6A, B). In TCGA CHOL samples, LAIR2 expression was increased in tumor samples (Figure 6C). After further clustering of CD8+/GZMB+ T cells into five sub-clusters, we found four of them showed expression of LAIR2, in which sub-cluster 2 demonstrated higher expression (Figures 6D, E). In addition, in comparison to immune stimulators (CD28, CD40), immune inhibitors (TGFB1, CD96, TIGIT, and LAG3) were highly expressed by all those sub-clusters, especially clusters 1 and 3 (Figure 6F). We further correlated LAIR2 expression with Treg scores and CD8+ T cell scores in TCGA CHOL samples, calculated by CIBERSORT and ssGSEA methods, and results showed LAIR2 was correlated to those cell populations (Figures 6G–J). We correlated LAIR2 with Treg cell markers (CD4, FOXP3, CD25, and CD39) and CD8+ T cell markers (CD8A, GZMB, TIM3, and PD1) in CHOL dataset, which all demonstrated high coefficients. Though the correlation between PD1 and LAIR2 was obvious, the corresponding coefficient did not achieve significance (Figures 6O–R). Considering immune regulation was performed through cooperation of immune regulators, we additionally analyzed the correlation between LAIR2, PNOC, and commonly acknowledged immune regulators, and results showed both LAIR2 and PNOC were significantly highly correlated with a bunch of immune inhibitors and stimulators in TCGA CHOL sequencing samples (Figure 7).
Figure 6 LAIR2 Was Highly Expressed by Regulatory T Cells and CD8+/GZMB+ T Cell Subset. (A) TSNE reduction for demonstration of NK and T cell atlas. (B) Markers for sub T and NK cell populations. (C) LAIR2 expression levels between CCA tumor and normal tissues in TCGA database. (D) Further cluster of CD8+/GZMB+ T cells. (E) LAIR2 expression in further clustered CD8+/GZMB+ T cell sub-populations. (F) Functional markers’ expression levels between further clustered CD8+/GZMB+ T cell sub-populations. (G–J) Correlation between LAIR2 expression and scores for regulatory T cell and CD8+ T cell calculated by CIBERSORT or ssGSEA method. (K–N) Correlation between LAIR2 and Treg markers [CD4, FOXP3, IL2RA (CD25), and ENTPD1 (CD39)] in CHOL bulk sequencing data from GEPIA database. (O–R) Correlation between CD8A, GZMB, HAVCR2 (TIM3), PDCD1 (PD1), and LAIR2 in CHOL bulk sequencing data.
Figure 7 Correlation Between LAIR2, PNOC, and Acknowledged Immune Checkpoints in TCGA CHOL and LIHC Datasets. (A) Correlation between expression of PNOC and immune regulators in CHOL dataset. (B) Correlation between expression of LAIR2 and immune regulators in CHOL dataset. (C) Correlation between expression of PNOC and immune regulators in LIHC dataset. (D) Correlation between expression of LAIR2 and immune regulators in LIHC dataset. (Immune inhibitors were marked with light orange, while immune stimulators were marked with light green. P values over 0.05 were not significant and were marked with white color.)
CCA Cells Demonstrated Heterogeneous Pathway Changes in Single Cell Level, Which Indicated Functional Variance and Malignant Potentials of Different Cancer Cell Clusters
We further extracted malignancy expression matrix from chonlangiocytes’ expression and calculated REACTOME pathway scores for each cell. After clustering of malignant cells according to calculated scores, cells were clustered into 11 populations (Supplementary Figures 2A, B). Of notice, clusters 11, 2, 7, 8, and 9 demonstrated highly malignant traits with high expression of signatures in cell mitotic cycle, IL1 signaling, PD1 signaling, and PI3K signaling (Supplementary Figure 2C).
Discussion
In our analysis, we used bulk sequencing data of cholangiocarcinoma patients in TCGA database to calculate the immune infiltration scores of different immune cell populations, and then we compared expression difference between groups, locating immune infiltration highly associated genes; we found PNOC was mainly expressed by infiltrated B cells, which was survival related, while LAIR2 was mainly expressed by Tregs and partial CD8+/GZMB+ T cells, indicating exhaustive immune status of T cells.
Prepronociceptin (PNOC) was formerly reported to be a pre-protein for a series of products, which act as pain regulators in signal transduction (38). Recent study showed PNOC is involved in long-term opioid response, alcoholic states, and inflammation process (39–44). In cancer investigations, PNOC was related to high inflammatory status and oxidative stress in C6 glioma cells, which also was highly up-regulated in pediatric brainstem ganglioglioma tissues and epithelial ovarian cancer, and in analysis of high-risk gastrointestinal stromal tumor, PNOC was reported as a prognostic biomarker (45–48). In a study of mRNA and microRNA network in colorectal cancer, PNOC and its targeting microRNAs were also prognostic markers for evaluation of patients (49). In our analysis, PNOC was up-regulated in cholangiocarcinoma samples, though the difference didn’t achieve significance. Low expression of PNOC may explain why only one model of GEO datasets included PNOC for prediction of high-immune infiltration; further analysis showed, PNOC was highly expressed by B cell populations in TIME. Expression of PNOC and infiltrating levels of B cell populations in CHOL were both survival-related, and in our analysis, differentially expressed immune genes between high- and low-immune infiltration groups were mainly immunoglobulins, indicating B cell infiltration was crucial in humoral anti-tumor responses. We also examined the prognostic values of PNOC in hepatocellular carcinoma, and we found patients with high PNOC expression also had better overall survival, indicating PNOC could be an independent biomarker for patients’ evaluation (Supplementary Figure 3B). High PNOC expression could predict highly infiltrated TIME. The specific roles of PNOC, expressed by B cells in TIME regulation, still need further experiments to illustrate.
Leukocyte associated immunoglobulin like receptor 2 (LAIR2) was previously reported as a close member to leukocyte associated immunoglobulin like receptor 1 (LAIR1), which is deemed as an immune inhibitor expressed by immune cells (50–52). According to publications, expression of LAIR1 was detected in various immune cell populations, and it has immune receptor tyrosine-based inhibition motif (ITIM), recruiting SHP-1, SHP-2, and Src kinase after phosphorylation (53, 54). Collagen in tumor matrix and damaged tissues is a common ligand for LAIR1 in broad spectrum, inhibiting immune cell functions after ligation, while LAIR2 was found to be a soluble protein with similar extracellular domain, which could block LAIR1 binding by competing ligands (55–58). Former studies also showed in autoimmune diseases, expression of LAIR2 was increased, and genetic single nucleotide polymorphism of LAIR2 was related to susceptibility of autoimmune diseases (59–62). The knowledge of LAIR2 in TIME regulation is limited, however, LAIR2 could interfere platelet activation and adhesion, and secreted LAIR2 could inhibit classical and lectin pathways of complement system in killing pathogens (58, 63). Overexpression of LAIR2 in lung cancer could increase immune infiltration levels and rescue exhaustive CD8+ T cells’ function (58). In our analysis, the mRNA expression of LAIR2 in Tregs and partial CD8+/GZMB+ T cells, in comparison to LAIR1, which was widely expressed by various cell populations, could be an indicator of exhaustive immune status. Former studies showed LAIR1 was expressed by macrophages, dendritic cells, as well as other CD14+ cells in inflammation, and scientists hypothesized LAIR1 could be both a threshold receptor and negative feedback receptor (64, 65). In our analysis, LAIR1 was not related to immune infiltration status in CHOL, and we believed mRNA expression of LAIR2 may be increased to offset LAIR1’s function in feedback loop, highlighting the role of baseline LAIR2 expression. Though LAIR2 in bulk sequencing data of CHOL was not survival-related, high LAIR2 expression in LIHC could indicate worse prognosis (Supplementary Figure 3A). We further examined coefficients for correlations between LAIR2 and other acknowledged immune regulators, such as CD28, LAG3, CD40, CXCR4, and TIGIT, in CHOL and LIHC datasets respectively, most of which achieved significance with high coefficients.
In addition, we analyzed the pathway changes in intrahepatic cholangiocarcinoma cell populations, finding functionally heterogeneous cancer cell clusters. Those malignant cells were clustered into 11 populations, in which several clusters showed high self-replication potentials, while others showed activated PI3K signal cascade through FGFR interaction. Also, cluster 2 and cluster 11 cells showed immune evasion potentials by increasing human lymphocyte associated antigens, and in cluster 7, expression of PDL1 (CD274) was increased. These functionally different cell populations in tumor justify the need of combining immune therapy in cholangiocarcinoma, and PNOC and LAIR2 could be clinical biomarkers for patient evaluation before immune therapy, predicting patients’ survival and tumor immune infiltration accordingly.
There are some limitations of our study. First, though we combined bulk sequencing and single cell sequencing data to characterize TIME of CHOL, protein expression were not examined, and further experiments should be conducted for confirmation. Second, immune regulating roles of PNOC, expressed by B cells, and LAIR2, expressed by Tregs and partial CD8+/GZMB+ T cells, in TIME are still in mist, which shall be further investigated.
Conclusion
High B cell infiltration level could indicate better prognosis in CCA. PNOC was mainly expressed by B cells in TIME and could be an independent indicator for better prognosis. LAIR2 was mainly expressed by Treg and partial CD8+/GZMB T cells, which could be an indicator for exhaustive T cell populations in TME. Both PNOC and LAIR2 were correlated with high immune infiltration levels in CCA patients.
Data Availability Statement
Publicly available datasets were analyzed in this study. This data can be found here: CHOL and LIHC in TCGA database: https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga: GSE32225: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32225, GSE26566: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE26566, GSE138709: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138709.
Ethics Statement
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.
Author Contributions
ZC, MY, JY, LG, BZ, and SL contributed to designing and analyzing process of the investigation, drafting the manuscript afterwards. JL, WZ, BHZ, JG, ZY, and XL helped to collect data and perform analysis correspondingly. JZ, JF, QY, HL, YFX, and YSX reviewed the whole manuscript, adapting the manuscript for final submission. All authors contributed to the article and approved the submitted version.
Funding
This research was funded by National Nature Science Foundation of China (81871924, 81572844, 81572301, 81802893, 81502487, 81972829), Natural Science Foundation of Guangdong Province of China (Grant Nos. 2017A030310641, 2018A030313762), and Foundation of Shanghai Municipal Commission of Health and Family Planning (20174Y0072).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.647209/full#supplementary-material
Supplementary Table 1 | Immune Infiltration Classification of Patients in Datasets of GSE26566 and GSE32225.
Supplementary Table 2 | Coefficients for LAIR2 and PNOC in Association With Treg, CD8+ T cells and B cells in Bulk Sequencing Data of CCA in TCGA Database Calculated by TIMER 2.0.
Supplementary Table 3 | Clinical information for CHOL bulk sequencing samples in TCGA database.
Supplementary Figure 1 | Expression of Immune Inhibitors and Stimulators in High- and Low-Immune Infiltration Patients.
Supplementary Figure 2 | Heterogeneous Pathway Changes in Intrahepatic Cholangiocarcinoma Cells Demonstrated Different Functional Status of Sub-cell Populations. (A) Heatmap for pathway scores of different cholangiocarcinoma cell sub-populations. (B) UMAP reduction for demonstration of cholangiocarcinoma cell sub-populations. (C) Heatmap for genes’ expression in PD1 signaling, Cell Cycle Mitotic signaling, IL-1 signaling, and PI3K-FGFR signaling between groups.
Supplementary Figure 3 | Both of PNOC and LAIR2 Were Related to Overall Survival of HCC Patients. (A) High expression of LAIR2 indicated worse survival in HCC patients. (B) High expression of PNOC indicated better survival in HCC patients.
Supplementary Figure 4 | ROC Plots for Immune Infiltration Models’ Evaluation. (A) ROC curves for regression model of immune infiltration score and each infiltration-related gene in dataset of GSE26566. (B) ROC curves for regression model of immune infiltration score and each infiltration-related gene in dataset of GSE32225 (AUC, area under curve).
References
1. Jain A, Kwong LN, Javle M. Genomic Profiling of Biliary Tract Cancers and Implications for Clinical Practice. Curr Treat options Oncol (2016) 17(11):58. doi: 10.1007/s11864-016-0432-2
2. Chun YS, Javle M. Systemic and Adjuvant Therapies for Intrahepatic Cholangiocarcinoma. Cancer control J Moffitt Cancer Center (2017) 24(3):1073274817729241. doi: 10.1177/1073274817729241
3. Sabbatino F, Villani V, Yearley JH, Deshpande V, Cai L, Konstantinidis IT, et al. PD-L1 and HLA Class I Antigen Expression and Clinical Course of the Disease in Intrahepatic Cholangiocarcinoma. Clin Cancer Res (2016) 22(2):470–8. doi: 10.1158/1078-0432.CCR-15-0715
4. Kriegsmann M, Roessler S, Kriegsmann K, Renner M, Longuespee R, Albrecht T, et al. Programmed cell death ligand 1 (PD-L1, CD274) in cholangiocarcinoma - correlation with clinicopathological data and comparison of antibodies. BMC Cancer (2019) 19(1):72. doi: 10.1186/s12885-018-5254-0
5. Rimassa L, Personeni N, Aghemo A, Lleo A. The immune milieu of cholangiocarcinoma: From molecular pathogenesis to precision medicine. J Autoimmun (2019) 100:17–26. doi: 10.1016/j.jaut.2019.03.007
6. Saeed A, Park R, Al-Jumayli M, Al-Rajabi R, Sun W. Biologics, Immunotherapy, and Future Directions in the Treatment of Advanced Cholangiocarcinoma. Clin colorectal Cancer (2019) 18(2):81–90. doi: 10.1016/j.clcc.2019.02.005
7. Fontugne J, Augustin J, Pujals A, Compagnon P, Rousseau B, Luciani A, et al. PD-L1 expression in perihilar and intrahepatic cholangiocarcinoma. Oncotarget (2017) 8(15):24644–51. doi: 10.18632/oncotarget.15602
8. Ghidini M, Cascione L, Carotenuto P, Lampis A, Trevisani F, Previdi MC, et al. Characterisation of the immune-related transcriptome in resected biliary tract cancers. Eur J Cancer (Oxford Engl 1990) (2017) 86:158–65. doi: 10.1016/j.ejca.2017.09.005
9. Mou H, Yu L, Liao Q, Hou X, Wu Y, Cui Q, et al. Successful response to the combination of immunotherapy and chemotherapy in cholangiocarcinoma with high tumour mutational burden and PD-L1 expression: a case report. BMC Cancer (2018) 18(1):1105. doi: 10.1186/s12885-018-5021-2
10. Ren X, Zhang Z. Understanding tumor-infiltrating lymphocytes by single cell RNA sequencing. Adv Immunol (2019) 144:217–45. doi: 10.1016/bs.ai.2019.08.004
11. Bao X, Shi R, Zhao T, Wang Y, Anastasov N, Rosemann M, et al. Integrated analysis of single-cell RNA-seq and bulk RNA-seq unravels tumour heterogeneity plus M2-like tumour-associated macrophage infiltration and aggressiveness in TNBC. Cancer Immunol Immunother (2020) 70(1):189–202. doi: 10.1007/s00262-020-02669-7
12. Lee HW, Chung W, Lee HO, Jeong DE, Jo A, Lim JE, et al. Single-cell RNA sequencing reveals the tumor microenvironment and facilitates strategic choices to circumvent treatment failure in a chemorefractory bladder cancer patient. Genome Med (2020) 12(1):47. doi: 10.1186/s13073-020-00741-6
13. Tomczak K, Czerwinska P, Wiznerowicz M. The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol (Poznan Poland) (2015) 19(1a):A68–77. doi: 10.5114/wo.2014.47136
14. Andersen JB, Spee B, Blechacz BR, Avital I, Komuta M, Barbour A, et al. Genomic and genetic characterization of cholangiocarcinoma identifies therapeutic targets for tyrosine kinase inhibitors. Gastroenterology (2012) 142(4):1021–31.e15. doi: 10.1053/j.gastro.2011.12.005
15. Sia D, Hoshida Y, Villanueva A, Roayaie S, Ferrer J, Tabak B, et al. Integrative molecular analysis of intrahepatic cholangiocarcinoma reveals 2 classes that have different outcomes. Gastroenterology (2013) 144(4):829–40. doi: 10.1053/j.gastro.2013.01.001
16. Zhang M, Yang H, Wan L, Wang Z, Wang H, Ge C, et al. Single cell transcriptomic architecture and intercellular crosstalk of human intrahepatic cholangiocarcinoma. J hepatology (2020) 73(5):1118–30. doi: 10.1016/j.jhep.2020.05.039
17. Ru B, Wong CN, Tong Y, Zhong JY, Zhong SSW, Wu WC, et al. TISIDB: an integrated repository portal for tumor-immune system interactions. Bioinf (Oxford England) (2019) 35(20):4200–2. doi: 10.1093/bioinformatics/btz210
18. Foroutan M, Bhuva DD, Lyu R, Horan K, Cursons J, Davis MJ. Single sample scoring of molecular phenotypes. BMC Bioinf (2018) 19(1):404. doi: 10.1186/s12859-018-2435-4
19. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinf (Oxford England) (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616
20. McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res (2012) 40:4288–97. doi: 10.1093/nar/gks042
21. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102
22. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc (2009) 4(1):44–57. doi: 10.1038/nprot.2008.211
23. Fabregat A, Jupe S, Matthews L, Sidiropoulos K, Gillespie M, Garapati P, et al. The Reactome Pathway Knowledgebase. Nucleic Acids Res (2018) 46:D649–D55. doi: 10.1093/nar/gkx1132
24. R Core Team. R: A language and environment for statistical computing Vienna, Austria: R Foundation for Statistical Computing (2019). Available at: https://www.R-project.org/.
25. RStudio Team. .RStudio: Integrated Development for R. Boston, MA: RStudio, Inc. (2015). Available at: http://www.rstudio.com/.
26. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, 3rd, et al. Comprehensive Integration of Single-Cell Data. Cell (2019) 177(7):1888–902.e21. doi: 10.1016/j.cell.2019.05.031
27. Zhang Y, Ma Y, Huang Y, Zhang Y, Jiang Q, Zhou M, et al. Benchmarking algorithms for pathway activity transformation of single-cell RNA-seq data. Comput Struct Biotechnol J (2020) 18:2953–61. doi: 10.1016/j.csbj.2020.10.007
28. Zhang Y, Zhang Y, Hu J, Zhang J, Guo F, Zhou M, et al. scTPA: a web tool for single-cell transcriptome analysis of pathway activation signatures. Bioinf (Oxford England) (2020) 36(14):4217–9. doi: 10.1093/bioinformatics/btaa532
29. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337
30. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol (2016) 17(1):218. doi: 10.1186/s13059-016-1113-y
31. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res (2017) 77(21):e108–e10. doi: 10.1158/0008-5472.CAN-17-0307
32. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res (2017) 45(W1):W98–102. doi: 10.1093/nar/gkx247
33. Nagy Á, Lánczky A, Menyhárt O, Győrffy B. Validation of miRNA prognostic power in hepatocellular carcinoma using expression data of independent datasets. Sci Rep (2018) 8(1):9227. doi: 10.1038/s41598-018-27521-y
34. Alboukadel K, Marcin K, Przemyslaw B. survminer: Drawing Survival Curves using ‘ggplot2’. R package version 0.4.8. (2020). https://CRAN.R-project.org/package=survminer.
35. Kolde R. pheatmap: Pretty Heatmaps. R package version 1.0.12. (2019). https://CRAN.R-project.org/package=pheatmap.
37. Walter W, Sánchez-Cabo F, Ricote M. GOplot: an R package for visually combining expression data with functional analysis. Bioinf (Oxford England) (2015) 31(17):2912–4. doi: 10.1093/bioinformatics/btv300
38. Rodriguez-Romaguera J, Ung RL, Nomura H, Otis JM, Basiri ML, Namboodiri VMK, et al. Prepronociceptin-Expressing Neurons in the Extended Amygdala Encode and Promote Rapid Arousal Responses to Motivationally Salient Stimuli. Cell Rep (2020) 33(6):108362. doi: 10.1016/j.celrep.2020.108362
39. Kuzmin A, Bazov I, Sheedy D, Garrick T, Harper C, Bakalkin G. Expression of pronociceptin and its receptor is downregulated in the brain of human alcoholics. Brain Res (2009) 1305 Suppl(Suppl):S80–5. doi: 10.1016/j.brainres.2009.05.067
40. Zhang L, Stuber F, Stamer UM. Inflammatory mediators influence the expression of nociceptin and its receptor in human whole blood cultures. PLoS One (2013) 8(9):e74138. doi: 10.1371/journal.pone.0074138
41. Peppin JF, Raffa RB. Delta opioid agonists: a concise update on potential therapeutic applications. J Clin Pharm Ther (2015) 40(2):155–66. doi: 10.1111/jcpt.12244
42. Carboni L, Romoli B, Romualdi P, Zoli M. Repeated nicotine exposure modulates prodynorphin and pronociceptin levels in the reward pathway. Drug Alcohol Depend (2016) 166:150–8. doi: 10.1016/j.drugalcdep.2016.07.002
43. Guzeldemir-Akcakanat E, Sunnetci-Akkoyunlu D, Orucguney B, Cine N, Kan B, Yılmaz EB, et al. Gene-Expression Profiles in Generalized Aggressive Periodontitis: A Gene Network-Based Microarray Analysis. J periodontology (2016) 87(1):58–65. doi: 10.1902/jop.2015.150175
44. Randesi M, van den Brink W, Levran O, Blanken P, Butelman ER, Yuferov V, et al. Variants of opioid system genes are associated with non-dependent opioid use and heroin dependence. Drug Alcohol Depend (2016) 168:164–9. doi: 10.1016/j.drugalcdep.2016.08.634
45. Chan MH, Kleinschmidt-Demasters BK, Donson AM, Birks DK, Foreman NK, Rush SZ. Pediatric brainstem gangliogliomas show overexpression of neuropeptide prepronociceptin (PNOC) by microarray and immunohistochemistry. Pediatr Blood Cancer (2012) 59(7):1173–9. doi: 10.1002/pbc.24232
46. Jin S, Zhu W, Li J. Identification of key genes related to high-risk gastrointestinal stromal tumors using bioinformatics analysis. J Cancer Res Ther (2018) 14(Supplement):S243–s7. doi: 10.4103/0973-1482.207068
47. Zhang Y, Luan D, Liu Y, Li H, Dong J, Zhang X, et al. Helicid Reverses Lipopolysaccharide-Induced Inflammation and Promotes GDNF Levels in C6 Glioma Cells through Modulation of Prepronociceptin. Chem Biodiversity (2020) 17(7):e2000063. doi: 10.1002/cbdv.202000063
48. Mamoor S. Prepronociceptin is differentially expressed in epithelial ovarian cancer. (2020). doi: 10.31219/osf.io/m382t
49. Jiang D, Xie X, Lu Z, Liu L, Qu Y, Wu S, et al. Establishment of a Colorectal Cancer-Related MicroRNA-mRNA Regulatory Network by Microarray and Bioinformatics. Front Genet (2020) 11:560186. doi: 10.3389/fgene.2020.560186
50. Lebbink RJ, de Ruiter T, Adelmeijer J, Brenkman AB, van Helvoort JM, Koch M, et al. Collagens are functional, high affinity ligands for the inhibitory immune receptor LAIR-1. J Exp Med (2006) 203(6):1419–25. doi: 10.1084/jem.20052554
51. Meyaard L. The inhibitory collagen receptor LAIR-1 (CD305). J Leukocyte Biol (2008) 83(4):799–803. doi: 10.1189/jlb.0907609
52. Zarrin AA, Monteiro RC. Editorial: The Role of Inhibitory Receptors in Inflammation and Cancer. Front Immunol (2020) 11:633686. doi: 10.3389/fimmu.2020.633686
53. Meyaard L, Hurenkamp J, Clevers H, Lanier LL, Phillips JH. Leukocyte-associated Ig-like receptor-1 functions as an inhibitory receptor on cytotoxic T cells. J Immunol (Baltimore Md 1950) (1999) 162(10):5800–4.
54. Meyaard L. LAIR and collagens in immune regulation. Immunol Lett (2010) 128(1):26–8. doi: 10.1016/j.imlet.2009.09.014
55. Lebbink RJ, van den Berg MC, de Ruiter T, Raynal N, van Roon JA, Lenting PJ, et al. The soluble leukocyte-associated Ig-like receptor (LAIR)-2 antagonizes the collagen/LAIR-1 inhibitory immune interaction. J Immunol (Baltimore Md 1950) (2008) 180(3):1662–9. doi: 10.4049/jimmunol.180.3.1662
56. Lebbink RJ, Raynal N, de Ruiter T, Bihan DG, Farndale RW, Meyaard L. Identification of multiple potent binding sites for human leukocyte associated Ig-like receptor LAIR on collagens II and III. Matrix Biol J Int Soc Matrix Biol (2009) 28(4):202–10. doi: 10.1016/j.matbio.2009.03.005
57. Lenting PJ, Westerlaken GH, Denis CV, Akkerman JW, Meyaard L. Efficient inhibition of collagen-induced platelet activation and adhesion by LAIR-2, a soluble Ig-like receptor family member. PLoS One (2010) 5(8):e12174. doi: 10.1371/journal.pone.0012174
58. Peng DH, Rodriguez BL, Diao L, Chen L, Wang J, Byers LA, et al. Collagen promotes anti-PD-1/PD-L1 resistance in cancer through LAIR1-dependent CD8(+) T cell exhaustion. Nat Commun (2020) 11(1):4520. doi: 10.1038/s41467-020-18298-8
59. Olde Nordkamp MJ, van Roon JA, Douwes M, de Ruiter T, Urbanus RT, Meyaard L. Enhanced secretion of leukocyte-associated immunoglobulin-like receptor 2 (LAIR-2) and soluble LAIR-1 in rheumatoid arthritis: LAIR-2 is a more efficient antagonist of the LAIR-1-collagen inhibitory interaction than is soluble LAIR-1. Arthritis Rheum (2011) 63(12):3749–57. doi: 10.1002/art.30612
60. Simone R, Pesce G, Antola P, Merlo DF, Bagnasco M, Saverino D. Serum LAIR-2 is increased in autoimmune thyroid diseases. PLoS One (2013) 8(5):e63282. doi: 10.1371/journal.pone.0063282
61. Camargo CM, Augusto DG, Petzl-Erler ML. Differential gene expression levels might explain association of LAIR2 polymorphisms with pemphigus. Hum Genet (2016) 135(2):233–44. doi: 10.1007/s00439-015-1626-6
62. Farias TDJ, Augusto DG, de Almeida RC, Malheiros D, Petzl-Erler ML. Screening the full leucocyte receptor complex genomic region revealed associations with pemphigus that might be explained by gene regulation. Immunology (2019) 156(1):86–93. doi: 10.1111/imm.13003
63. Olde Nordkamp MJ, Boross P, Yildiz C, Jansen JH, Leusen JH, Wouters D, et al. Inhibition of the classical and lectin pathway of the complement system by recombinant LAIR-2. J Innate Immun (2014) 6(3):284–92. doi: 10.1159/000354976
64. Carvalheiro T, Garcia S, Pascoal Ramos MI, Giovannone B, Radstake T, Marut W, et al. Leukocyte Associated Immunoglobulin Like Receptor 1 Regulation and Function on Monocytes and Dendritic Cells During Inflammation. Front Immunol (2020) 11:1793. doi: 10.3389/fimmu.2020.01793
Keywords: cholangiocarcinoma, immune infiltration, biomarker, single-cell sequencing technology, immunotherapy
Citation: Chen Z, Yu M, Yan J, Guo L, Zhang B, Liu S, Lei J, Zhang W, Zhou B, Gao J, Yang Z, Li X, Zhou J, Fan J, Ye Q, Li H, Xu Y and Xiao Y (2021) PNOC Expressed by B Cells in Cholangiocarcinoma Was Survival Related and LAIR2 Could Be a T Cell Exhaustion Biomarker in Tumor Microenvironment: Characterization of Immune Microenvironment Combining Single-Cell and Bulk Sequencing Technology. Front. Immunol. 12:647209. doi: 10.3389/fimmu.2021.647209
Received: 29 December 2020; Accepted: 03 March 2021;
Published: 24 March 2021.
Edited by:
Wei Wei, Institute for Systems Biology (ISB), United StatesReviewed by:
Chutamas Thepmalee, University of Phayao, ThailandDavid Linnaeus Gibbs, Institute for Systems Biology (ISB), United States
Copyright © 2021 Chen, Yu, Yan, Guo, Zhang, Liu, Lei, Zhang, Zhou, Gao, Yang, Li, Zhou, Fan, Ye, Li, Xu and Xiao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Qinghai Ye, ye.qinghai@zs-hospital.sh.cn; Hui Li, li.hui1@zs-hospital.sh.cn; Yongfeng Xu, xu.yongfeng@zs-hospital.sh.cn; Yongsheng Xiao, xiao.yongsheng@zs-hospital.sh.cn
†These authors have contributed equally to this work and share first authorship