Intercellular Molecular Crosstalk Networks within Invasive and Immunosuppressive Tumor Microenvironment Subtypes Associated with Clinical Outcomes in Four Cancer Types

Abstract

Heterogeneity is a critical basis for understanding how the tumor microenvironment (TME) contributes to tumor progression. However, an understanding of the specific characteristics and functions of TME subtypes (subTMEs) in the progression of cancer is required for further investigations into single-cell resolutions. Here, we analyzed single-cell RNA sequencing data of 250 clinical samples with more than 200,000 cells analyzed in each cancer datum. Based on the construction of an intercellular infiltration model and unsupervised clustering analysis, four, three, three, and four subTMEs were revealed in breast, colorectal, esophageal, and pancreatic cancer, respectively. Among the subTMEs, the immune-suppressive subTME (subTME-IS) and matrix remodeling with malignant cells subTME (subTME-MRM) were highly enriched in tumors, whereas the immune cell infiltration subTME (subTME-ICI) and precancerous state of epithelial cells subTME (subTME-PSE) were less in tumors, compared with paracancerous tissues. We detected and compared genes encoding cytokines, chemokines, cytotoxic mediators, PD1, and PD-L1. The results showed that these genes were specifically overexpressed in different cell types, and, compared with normal tissues, they were upregulated in tumor-derived cells. In addition, compared with other subTMEs, the expression levels of PDCD1 and TGFB1 were higher in subTME-IS. The Cox proportional risk regression model was further constructed to identify possible prognostic markers in each subTME across four cancer types. Cell-cell interaction analysis revealed the distinguishing features in molecular pairs among different subTMEs. Notably, ligand-receptor gene pairs, including COL1A1-SDC1, COL6A2-SDC1, COL6A3-SDC1, and COL4A1-ITGA2 between stromal and tumor cells, associated with tumor invasion phenotypes, poor patient prognoses, and tumor advanced progression, were revealed in subTME-MRM. C5AR1-RPS19, LGALS9-HAVCR2, and SPP1-PTGER4 between macrophages and CD8+ T cells, associated with CD8+ T-cell dysfunction, immunosuppressive status, and tumor advanced progression, were revealed in subTME-IS. The spatial co-location information of cellular and molecular interactions was further verified by spatial transcriptome data from colorectal cancer clinical samples. Overall, our study revealed the heterogeneity within the TME, highlighting the potential pro-invasion and pro-immunosuppressive functions and cellular infiltration characteristics of specific subTMEs, and also identified the key cellular and molecular interactions that might be associated with the survival, invasion, immune escape, and classification of cancer patients across four cancer types.

Keywords: C5AR1-RPS19; COL1A1-SDC1; COL6A3-SDC1; LGALS9-HAVCR2; SPP1-PTGER4; cell-cell interaction; immunosuppressive; single-cell RNA-seq; tumor invasion; tumor microenvironment subtypes.

Figure 1

TME subtypes across cancers. (a) UMAP plots of cells from normal and tumor tissue of CRC patients, showing 6 clusters indicating major cell types and 2 clusters indicating cells derived from tumor and normal tissues. Each cluster is shown in a different color. (b) The same as shown in a but in BC patients. (c) The three cellular modules on the basis of correlations of cell subclusters from tumors with positive (Spearman correlation; correlation coefficient r > 0.3 and FDR < 0.05, in red), negative (r < −0.3 and FDR < 0.05, in blue), or non-significant (white) pairwise correlation for infiltration in CRC samples. (d) The four cellular modules on the basis of correlations of cell subclusters from tumors in BC samples. (e,f) subTME score in normal and tumor sample groups in CRC and BC. *** p < 0.001, * p < 0.05, Wilcoxon tests.

Figure 2

Definition of TME subtypes. (a) Boxplot showing the anti-inflammatory score among myeloid cells in BC and CRC samples. *** p < 0.001, Kruskal–Wallis test. (b) Boxplot showing the checkpoint molecules and T-cell exhaustion score among subTMEs in BC. ** p < 0.01, Kruskal–Wallis test. (c) Boxplot showing the myCAF score among fibroblasts in PDAC and ESCA samples. *** p < 0.001, Kruskal–Wallis test. (d) Boxplot showing the pEMT score among epithelial cells in PDAC and ESCA samples. *** p < 0.001, Kruskal–Wallis test. (e) Boxplot showing the matrix remodeling and EMT score among subTMEs in PDAC. ** p < 0.01, Kruskal–Wallis test. (f) Definitions of the TME subtypes across cancers. (g) Scatterplot showing the Spearman correlation between immune escape score and subTME-IS score in BC, CRC, and PDAC, the error band indicates 95% confidence interval.

Figure 3

KM survival analysis, risk score assessment using the subTME gene markers, and time-dependent ROC curves in CRC and PDAC datasets. (a) KM survival analysis of high- and low-risk samples. (b) Relationship between the survival status/risk score rank and survival time/risk score rank. (c) Time-dependent ROC curve for overall survival of the CRC datasets. (d–f) The same as in (a–c) but in PDAC datasets.

Figure 4

Cell-cell interaction across cell subclusters within TME subtypes. (a) Social graph depicting the number of interactions between cell types within subTME-IS in BC. (b) Social graph depicting the number of interactions between cell types within subTME-IS in CRC. (c) Enrichment of selected ligand-receptor interactions in CRC for the cell-type pairs within subTME-IS. (d) Enrichment of selected ligand–receptor interactions in ESCA for the cell-type pairs within subTME-MRM. The bubbles shown in the figure indicate p < 0.05; the color of the bubbles indicates the level of mean gene expression of ligand–receptor.

Figure 5

EMT phenotype of tumor cells regulated by stromal cells. (a) Heatmap showing the Spearman correlation between gene expression and pEMT score in tumor cells, p < 0.05. (b,c) Scatterplot showing the Spearman correlation of the gene expression in fibroblasts and cell proportion of tumor cells in BC and CRC, the error band indicates 95% confidence interval. (d) Bubble chart showing the gene expression of ligand–receptor pairs between normal and tumor samples derived from fibroblasts and tumor cells, respectively, in single-cell data; the color of the bubbles indicates the level of gene expression of a ligand or receptor in sender or receiver cells; the bubbles shown in the figure indicate p < 0.05, Wilcoxon tests. (e) Boxplot showing the gene expression of ligand–receptor pairs between normal and tumor samples in TCGA data, *** p < 0.001, ** p < 0.01, Wilcoxon tests. (f) Scatterplot showing the Spearman correlation of the gene expression in TCGA-BC tumor data.

Figure 6

Clinical prognosis of ligand–receptor interactions between fibroblasts and tumor cells. (a) The Kaplan–Meier curve shows overall survival of COL1A1-SDC1, COL1A2-SDC1, and COL6A3-SDC1 in BC and PDAC patients. (b) Scatterplot showing the Spearman correlation of the gene expression in TCGA-PDAC tumor data, the error band indicates 95% confidence interval.

Figure 7

Immunosuppressive microenvironment regulated by TAMs. (a) Heatmap showing the Spearman correlation between gene expression and M2-like macrophage polarization score as well as CD163 expression in macrophages; color in red indicating positive (Spearman correlation; FDR < 0.05), color in blue indicating negative (FDR < 0.05), color in white indicating non-significant (FDR > 0.05). (b) Scatterplot showing the Spearman correlation of the gene expression in tumor cells and cell proportion of TAMs in CRC and PDAC, the error band indicates 95% confidence interval (c) Scatterplot showing the Spearman correlation of the gene expression in TAMs and cell proportion of CD8+ T cells in CRC, the error band indicates 95% confidence interval. (d) Bubble chart showing the gene expression of ligand–receptor pairs derived from tumor cells and TAMs or TAMs and CD8+ T cells, respectively, between normal and tumor samples in single-cell data in PDAC; the color of the bubbles indicates the level of gene expression of a ligand or receptor in sender or receiver cells; the bubbles shown in the figure indicate p < 0.05, Wilcoxon tests. (e) Boxplot showing the gene expression of ligand–receptor pairs between normal and tumor samples in TCGA-ESCA data, *** p < 0.001, ** p < 0.01, Wilcoxon tests. (f) Heatmap showing the Spearman correlation between gene expression levels in TAMs in CRC and PDAC, p < 0.05.

Figure 8

Co-localization of cell types in spatial transcriptomics data. (a) Spatial feature plots of signature score of ecm_myCAF (left) and Cancer_Malig (middle) in tissue sections and Spearman correlation of signature score of ecm_myCAF and Cancer_Malig (right) in CRC patient #19. (b) Spatial feature plots of signature score of Macro_APOE (left) and CD8Teff_2 (middle) in tissue sections and Spearman correlation of signature score of Macro_APOE and CD8Teff_2 (right) in patient #19. (c) The same as a but in patient #36. (d) The same as b but in patient #36.

Figure 9

Schema of cellular interactions between cell types of interest within both subTME-MRM and subTME-IS.