The Crosstalk Between Malignant Cells and Tumor-Promoting Immune Cells Relevant to Immunotherapy in Pancreatic Ductal Adenocarcinoma

Abstract

Background: Pancreatic ductal adenocarcinoma (PDAC) is dominated by an immunosuppressive microenvironment, which makes immune checkpoint blockade (ICB) often non-responsive. Understanding the mechanisms by which PDAC forms an immunosuppressive microenvironment is important for the development of new effective immunotherapy strategies. Methods: This study comprehensively evaluated the cell-cell communications between malignant cells and immune cells by integrative analyses of single-cell RNA sequencing data and bulk RNA sequencing data of PDAC. A Malignant-Immune cell crosstalk (MIT) score was constructed to predict survival and therapy response in PDAC patients. Immunological characteristics, enriched pathways, and mutations were evaluated in high- and low MIT groups. Results: We found that PDAC had high level of immune cell infiltrations, mainly were tumor-promoting immune cells. Frequent communication between malignant cells and tumor-promoting immune cells were observed. 15 ligand-receptor pairs between malignant cells and tumor-promoting immune cells were identified. We selected genes highly expressed on malignant cells to construct a Malignant-Immune Crosstalk (MIT) score. MIT score was positively correlated with tumor-promoting immune infiltrations. PDAC patients with high MIT score usually had a worse response to immune checkpoint blockade (ICB) immunotherapy. Conclusion: The ligand-receptor pairs identified in this study may provide potential targets for the development of new immunotherapy strategy. MIT score was established to measure tumor-promoting immunocyte infiltration. It can serve as a prognostic indicator for long-term survival of PDAC, and a predictor to ICB immunotherapy response.

Keywords: cell-cell communication; immunocyte infiltration; immunotherapy; pancreatic ductal adenocarcinoma; single cell RNA-seq.

FIGURE 1

ScRNA-seq and bulk RNA-seq reveals the immune microenvironment in PDAC. (A), heatmap showing the expression of 22 functional gene expression signatures in 33 cancer types from TCGA. (B), pie chart showing the proportion of each cluster for the normal samples (top panel) and tumor samples (bottom panel). (C), Uniform Manifold Approximation and Projection (UMAP) plot of CD4⁺ T cells (top panel), CD8⁺ T cells (middle panel) and myeloid cells (bottom panel), color-coded for subclusters (left panel) and tissue type (right panel), respectively. (D), UMAP plot of CD4⁺ T cells (top panel), CD8⁺ T cells (lower left panel), myeloid cells (lower right panel), color-coded for the average expression of functional gene expression signature. (E), violin plots showing the normalized expression of immune related genes in CD4⁺ T cell (left panel), CD8⁺ T cell (middle panel) and myeloid (right panel) subclusters.

FIGURE 2

Characteristics and evolution of epithelial cells during cancer progression and their ligand-receptor communications with immune cells in PDAC. (A), UMAP plot of epithelial cells, color-coded for the subclusters of epithelial cells. (B), Dot plot showing the expression of marker genes in each subcluster. (C), UMAP plot of epithelial cells, color-coded for the CNV score calculated by infercnv. (D), UMAP of epithelial cells, color-coded for the four major subclusters: acinar cells, malignant cells, ductal-normal cells and ductal-tumor cells. (E), dot plot showing the enriched pathways in the four major subclusters. The color of each dot represents the normalized enrichment score (NES), while the size of the dot represents p-value. (F), trajectory analysis showing the pseudotime of the four major subclusters. (G), line plot showing the expression of cell adhesion molecules (ITGA3, CD47, MDK, ITGB1, CD55, and CDH1) along the pseudotime. Each line with different color represents a gene. (H), circos plot showing the ligand-receptor interactions between immune cells and ductal-normal cells (left panel), ductal-tumor cells (middle panel), and malignant cells (right panel). (I), UMAP of epithelial cells showing the expression of known checkpoint genes (CD274, PVR, CD276 and ADORA2A). (J), the significantly enriched ligand-receptor pairs between in CD4⁺ T cells (left panel)/CD8⁺ T cells (middle panel)/myeloid cells (right panel) and ductal-normal cells/ductal-tumor cells/malignant cells.

FIGURE 3

Construction of Malignant-Immune-Talk (MIT) score to measure tumor-promoting immune microenvironment. (A), boxplot showing the differential expression of the 14 ligands from the above 15 ligand-receptor pairs between TCGA tumor and normal samples. *, p -value< 0.05; **, p -value< 0.01; ***, p-value < 0.001 from two-sided paired Student’s t-test. (B), immunohistochemical (IHC) staining of FAM3C, LGASL9, ANXA1, SPP1, CD99 and LAMP1 in normal tissue (top panel) or PDAC tumor tissue (bottom panel). Scale bar, 50 μm. (C), barplot showing the pathways enriched for the 14 ligands from the above 15 ligand-receptor pairs. (D), barplot showing the coefficients of the 7 genes selected from the lasso cox model. (E), Kaplan-Meier curve showing the overall survival difference between patients with high MIT socre and patients with low MIT score in the TCGA PAAD cohort. p-value was calculated by the log-rank test. (F), Kaplan-Meier curve showing the overall survival difference between patients with high MIT socre and patients with low MIT score in four independent cohorts of PDAC (GSE62452, GSE28735, GSE57495 and GSE85916). p-value was calculated by the log-rank test. (G), box plot showing the MIT scores among different clinical grades. (H), Heat map showing the differential expression of therapy-related genes between high and low MIT groups. (I), box plot showing the different MIT scores between chemotherapy responder group (pD/SD) and chemotherapy non-responder group (pR/CR). (J), box plot showing the difference in estimated of IC50 for four EGFR inhibitors between high and low MIT group.

FIGURE 4

MIT score is associated with tumor-promoting immune microenvironment. (A), enriched pathways for the upregulated genes in MIT high group compared to MIT low group. (B), the two enriched immune hallmark pathways in MIT high group compared to MIT low group derived from GSEA analysis. (C), heatmap showing the differential expression of immune-associated genes between MIT high and low group. The color bar represents the level of expression (Standardized by TPM). (D), boxplot showing the difference in proportions of immune cell types (derived from EPIC algorithm) between MIT high and low groups. *, p-value < 0.05; **, p-value < 0.01; ***, p-value < 0.001 from two-sided paired Student’s t-test. (E), boxplot showing the difference in proportions of immune cell types (derived from xCell algorithm) between MIT high and low groups. *, p-value < 0.05; **, p-value < 0.01; ***, p-value < 0.001 from two-sided paired Student’s t-test. (F), barplot showing the correlation between the MIT score and functional gene expression signatures in PDAC. (G), The correlation between the proportion of tumor-promoting immune clusters (CD4_CXCL13, CD4_FOXP3, CD8_ISG15, CD8_CXCL13, Mac_C1QA, Mac_SPP1, Mac_ISG15 and cDC3_LAMP3) and MIT score.

FIGURE 5

MIT-associated oncogenic events and their relevance to ICB response. (A), recurrent copy number events (GISTIC2 Q < 0.1) in MIT-high group (top panel) and MIT-low group (bottom panel). (B), comparison of mutation rates of individual genes between the MIT high and MIT low groups. *, p-value < 0.05; **, p-value < 0.01; ***, p-value < 0.001 (C), frequency of the mutation of KRAS, TP53, PDZRN3, ADAMTSL4 and CDKN2A in MSK-IMPACT dataset with ICB clinical information. (D), Kaplan-Meier curves showing the survival difference between patients with no mutation and patients having mutation in MSK-IMPACT dataset with ICB treatment. The p-value is computed via a two-sided log-rank test. Have mutation indicates patients with any of KRAS, TP53, PDZRN3, ADAMTSL4 and CDKN2A mutations. (E), boxplot showing the TMB difference between MIT-high group and MIT-low group. *, p -value< 0.05 from two-sided paired Student’s t-test. (F), The Kaplan-Meier curve showing the survival difference among four subgroups of patients. Have mutation representing patients with any of KRAS, TP53, PDZRN3, ADAMTSL4 and CDKN2A mutations. p-values were calculated using stratified log-rank test.

FIGURE 6

MIT score can be a prognostic factor related to the long-term efficacy of immunotherapy. (A), Heatmap shows the correlation between the MIT score and functional gene expression signatures in 33 cancer types from TCGA. (B–F), Kaplan Meier curves showing the survival difference between MIT-high and MIT-low group (left panel) and bar chart showing the difference in proportion of responders and non-responders between MIT-high and MIT-low group (right panel) in BLCA, CCRCC, SKCM patients, respectively. The p-values for survival difference are calculated via a two-sided log-rank test, and the p-values for proportion difference are calculated from fisher’s exact test.