Global cellular microenvironment landscape in gastric cancer
To characterize the TME in gastric cancer, we performed bulk tissue RNA sequencing (RNA-seq) and whole exome-sequencing together with single-cell RNA-seq (scRNA-seq) of tumor and matched non-malignant gastric tissue samples from 24 treatment-naive gastric cancer patients (Fig. 1a, Additional file 1: S1a and Additional file 2: Table S1). For identifying the cell types constituting the gastric TME, samples were rapidly digested into single cells, depleted for EPCAM-positive epithelial cells by fluorescent-activated cell sorting (FACS) to enrich for all cells of non-epithelial origin and analyzed using scRNA-seq. The resulting 96,623 cells were clustered into eleven major cell types and further into 81 subtypes (Fig. 1b, Additional file 1: Fig. S1f and “Methods”). The major clusters were annotated according to the expression of defining marker genes as either B cells (CD19, MS4A1), plasma cells (IGHG1, CD79A), CD4+ T cells (CD3D, CD4), CD8+ T cells (CD3D, CD8A), natural killer (NK) cells (NCR1, FGFBP2), myeloid cells (CD14, CD68), mast cells (TPSAB1, TPSB2), endothelial cells (RAMP2, PECAM1), fibroblasts (DCN, LUM), mural cells (PDGFRB, ACTA2), glial cells (PLP1, SOX10), or epithelial cells (PGC, PGA3) (Fig. 1c). To characterize the differences between the TME and normal gastric stroma in the following, we quantified the change in abundance and cellular activity associated with malignant transformation for each of the stromal cell types.
Tumor stroma is characterized by the presence of activated fibroblasts
To determine if any fibroblast and related mural cell subtypes are associated with malignant transformation, we performed re-clustering of the corresponding 16,492 cells in our dataset which resulted in 15 cell clusters (Fig. 2a). By evaluating the expression of specific marker genes across these clusters, we found mural cells to be clearly separated from fibroblasts and to be comprised of pericytes characterized for instance by the expression of RGS5, and smooth muscle cells expressing among other genes ACTA2 and DES (Additional file 1: Fig. S2a). Fibroblasts on the other hand were distinctly divided into resting and activated cells as could be inferred from the expression of markers for resting fibroblasts such as S100A4, CFD, and DPT and markers for active cells like POSTN, CXCL14, and COL3A1 [9, 16,17,18] (Fig. 2b, S2a). A closer look at the individual cell subtypes revealed multiple unique genes for each cluster, underlining the heterogeneity of both gastric fibroblast and mural cells (Fig. 2b). Applying consensus non-negative matrix factorization (cNMF) [19], we identified distinct gene expression programs for the different fibroblast subtypes confirming that the different clusters indeed represent either different cell subtypes or cell activation states (Additional file 1: Fig. S3a-b). The fibroblast subtypes further displayed different characteristics in extracellular matrix-related gene expression, pointing to possible distinct functions in shaping the extracellular matrix (Additional file 1: Fig. S3c-h). Cluster fraction comparison between malignant and non-malignant samples highlighted subtype F13-CTHRC1 as highly cancer-associated fibroblasts (CAFs) while revealing a drastic reduction of several resting fibroblast clusters in the gastric TME (Fig. 2b and S2b). In line with this finding, expression of CTHRC1 is upregulated in tumors from many different indications in the cancer genome atlas (TCGA) cohort (Additional file 1: Fig. S2c). Comparing subtype F13-CTHRC1 to a recently published fibroblast compendium revealed a strong similarity of these cells to activated myofibroblasts specifically expressing for instance COL3A1, ACTA2, and CTHRC1 [18]. Diffusion map analysis between resting and activated fibroblasts positioned resting fibroblasts F01-SLPI and F13-CTHRC1 at both ends of a trajectory (Additional file 1: Fig. S2d), pointing towards activated fibroblasts emerging from resting fibroblasts.
Expression of CTHRC1-positive fibroblast-specific genes is linked to poor survival
In addition to changes in cell type frequencies, many molecular pathways have been described to be deregulated in gastric tumors compared to healthy normal tissue [20]. To better characterize transcriptional changes associated with malignant transformation, we integrated our bulk and single-cell RNA-seq data and identified upregulated (designated by capital U) and downregulated (designated by a capital D) gene clusters that harbor different expression patterns across the single-cell populations (“Methods” and Additional file 1: Fig. S4a). Multiple upregulated and downregulated gene clusters specifically expressed in individual fibroblast subsets were identified (Additional file 1: Fig. S4b-c). Among upregulated clusters, U1 was expressed in resting and activated fibroblasts and strongly enriched with genes associated with gene ontology (GO) terms linked to remodeling of the extracellular matrix. In addition, clusters U6 and U7 were highly expressed in L12-ANGPT2 pericytes and F13-CTHRC1 CAFs, respectively (Fig. 2c, S4c). Cluster U6 thereby contained genes like ACAN and COL5A3 which are involved in extracellular matrix organization as well as genes like ANGPT2, EGFL6, PDGFRB, and PGF (Additional file 1: Fig. S4b, S4d) which are connected to growth factor responses and angiogenesis in the GO ontology. These genes have also been reported as highly expressed in multiple cancer types and as being associated with angiogenesis and tumor invasion [21]. Cluster U7 contained genes related to fibroblast proliferation (WNT5A and WNT2), matrix remodeling (MMP1 and MMP3), and cell migration (PDPN and TWIST1) (Fig. 2d, S4b,d). All genes from cluster U7 were specifically expressed in F13-CTHRC1 activated fibroblasts and were associated with poor survival in an independent gastric cancer cohort (Fig. 2e) as well as in several cancer indications from TCGA (Additional file 1: Fig. S2e-f). This is in alignment with previous reports that implicated these genes in cancer progression and bad outcome [22,23,24,25,26].
Tumor endothelial cells swap immune attraction for angiogenic pathways
Endothelial cells are key components of the TME since supply of energy and nutrients through active blood circulation, rather than passive diffusion, is required for continuous tumor growth [27]. Upon analyzing the 3684 endothelial cells in our data, we found nine clusters that corresponded to six different endothelial cell subtypes (Fig. 3a). These endothelial clusters were associated with arteries (GJA5, TSPAN2), capillaries (CA4, BTNL9), immature endothelial cells (HSPG2, VWA1), tip cells (PGF), lymphatic cells (PROX1, LYVE1), and vein cells (ACKR1) which could be further delineated into postcapillary cells (CPE), activated (POSTN), and IL6-expressing cells (Fig. 3b). This observation is in line with a recently reported taxonomy of endothelial cells obtained from healthy human lung and lung carcinoma samples [28], demonstrating the consistent character of these vasculature-associated cell types across tissues. Cluster fraction comparison between normal and malignant gastric tissue identified clusters EN10-SERPINE1 and EN03-ESM1, hereafter referred to as activated endothelial cells, as almost exclusively tumor-specific cell types (Fig. 3b, S2g).
To define the functionality of tumor endothelial cells, we first performed differential gene expression analysis between activated endothelial cells and the remaining endothelial cells (Additional file 4: Table S3). Subsequent gene set enrichment analysis performed on the significantly differentially expressed genes revealed an increase of angiogenic and cell motility response pathways together with a decrease of antigen-presenting and immune defense genes in both EN10-SERPINE1 and EN03-ESM1, indicating that immune cell attraction through endothelial cells is reduced in gastric cancer (Fig. 3c). This confirms observations from other cancer indications that point to a decrease of antigen-presenting pathways in endothelial cells [29, 30] and could be a reason for the synergistic effect of anti-angiogenic therapy together with ICI [31].
In line with angiogenesis being critical for sustained tumor growth, we speculated that these cell types might be linked to poor survival and found that SERPINE1 expression was associated with a drastically worsened survival in the independent gastric cancer cohort from TCGA (Fig. 3d). SERPINE1 is a serine protease inhibitor that mainly functions as a regulator of cell adhesion and spreading and has been reported as upregulated in multiple other cancer indications such as head and neck squamous cell carcinoma, leading to a poor prognosis [32]. SERPINs have been described as important genes in cancer and vascular co-option and therefore play an important role in cancer progression leading to tumor metastasis [33]. We found a comparable difference in survival using the signature of the 20 most specific genes of the EN10-SERPINE1 cluster (Fig. 3e). Together, these findings indicate that EN10-SERPINE1 endothelial cells could be targeted to block the metastatic spread of gastric tumors.
Macrophages display an inflammation dichotomy
To shed light on myeloid populations present in gastric tumors, we re-clustered the 10,546 myeloid cells in our dataset and identified fifteen subclusters which were subsequently determined to be either macrophages (CD68, CD14), dendritic cells (FLT3, FCER1A), or neutrophils (CSF3R) (Fig. 4a). The four dendritic cell clusters in our data matched dendritic cells found in hepatocellular and colorectal carcinoma [34, 35] and were annotated as conventional cDC1 M16-CLEC9A, cDC2 M03-CD1C, plasmacytoid M18-CLEC4C, and M17-LAMP3 dendritic cells (Fig. 4b). Through differential gene expression analysis, macrophage clusters were identified as either proinflammatory (S100A8, S100A9 [36], IL1B [37], and CXCL8 [38]), anti-inflammatory (APOE [39], MAF [40], C1QB, and SEPP1 [41]), or tissue-resident macrophages (F13A1 [42] and CCL2 [43]) (Fig. 4b). Tissue-resident macrophages were found to be distinct from anti-inflammatory macrophages through the expression of various genes including FOLR2, CCL2, LYVE1, SEPP1, and F13A1, all of which have been reported as markers of tissue-resident macrophages [44]. Myeloid cells were drastically increased in the TME with only tissue-resident macrophages higher in healthy stomach tissue. In particular, the gastric TME was characterized by a significant increase in M07-APOE, M11-SPP1, and M04-C3 anti-inflammatory macrophages (Fig. 4c, S5).
Surprisingly, we found the expression of the above pro- and anti-inflammatory gene signature to be anti-correlated across patient samples in both our data as well as in the independent stomach adenocarcinoma cohort from TCGA (Fig. 4d). In line with a recent publication that linked poor survival to the presence of anti-inflammatory myeloid cells [45], the expression of our anti-inflammatory signature was significantly associated with reduced survival in the TCGA gastric cancer cohort (Fig. 4e). However, the abundance of the corresponding macrophage populations as determined by CIBERSORTx deconvolution of the bulk RNA-seq samples was not predictive for outcome (data not shown).
Lymphocytes divert towards immunosuppressive and differentiation phenotypes
In contrast to CD8+ T cells, whose fraction stayed nearly unchanged between tumor and normal samples, we observed an over fivefold increase in CD4+ T cells in the gastric tumor samples. To understand which T cell subsets most differentially invaded the gastric TME, we reclustered the 12,905 CD4+ T cells into eight subclusters which we then assigned to six T helper cell subtypes including regulatory (FOXP3, IL2RA), helper 17 (IL17A), naïve (CCR7), central memory (CM) (ANXA1, CCR7), effector memory (EM) (ANXA1, CCL5), and follicular T helper cells (CXCL13) (Fig. 5a and S6a). In addition, the 31,705 CD8+ lymphocytes and 2,658 NK cells reclustered into thirdteen subclusters, 12 of which corresponded to eight different CD8 T cell subtypes and one corresponded to NK cells. The CD8+ T cells were designated as naïve (IL7R), effector (GLNY), effector memory (ANXA1, GZMK), effector memory expressing CD45RA (LMNA, CREM), resident memory (ANXA1, GZMB), exhausted (HAVCR2), intraepithelial (CD160) [46], and CD8+ T cells with a tertiary lymphoid signature (CXCL13) (Fig. 5b and Additional file 1: Fig. S6b). The expression profile associated with these T cell clusters closely resembled those of T cell subtypes reported in other cancer indications [47, 48]. Immunosuppressive Th17 and Treg cells were among the most increased CD4+ T cells while naive CD8+ T cells were reduced in tumors (Fig. 5c), likely reflecting the activation and expansion of T cells in the tumor. Besides gastric cancer, disproportionate CD4+ over CD8+ T cell ratios have also been observed in CRC [30, 49].
Our bulk and single-cell integration found one gene cluster U4 upregulated in gastric cancer that was preferentially expressed in T01-ICOS (Fig. 5d). Upon observing this cluster in more detail, we found many genes associated with activated regulatory T cell like FOXP3 and CCR8, inhibition of apoptosis such as PMAIP1, and immune suppression such as CD274, ALOX15, and SERPINB9 [50], emphasizing a strong activation of regulatory T cells and immune suppression in the gastric TME.
We identified 16,883 plasma cells and annotated them through gene signatures for heavy and light chain as well as kappa and lambda immunoglobulins, as IgAλ, IgAκ, IgGλ, IgGκ, IgMλ, IgMκ (Additional file 1: Fig. S7a). Comparing the ratio of plasma cell clusters between distal non-malignant gastric tissue and gastric tumor showed a high increase of IgG isotypes in gastric cancer whereas IgA isotypes decreased, suggesting a systemic change in the immune microenvironment (Additional file 1: Fig. S7b). The increase of IgG-positive plasma cells in gastric cancer is intriguing and has been documented in previous studies on other cancer indications [51].
Cell communication network uncovers potential drivers of gastric cancer development
To investigate cellular interactions in the gastric cancer tumor microenvironment, we constructed a cell communication network based on ligand and matching receptor expression information from our single-cell atlas. The microenvironment was found to display a wide range of communications across different cell subtypes. Communication models derived from the cells of the TME and from normal gastric tissue thereby differed significantly in their structure, with different cell subtypes ranked most central in the interaction networks for TME and normal tissue, and F13-CTHRC1 as a central network hub in the TME but not in normal tissue (Additional file 1: Fig. S8c-d). F13-CTHRC1 thereby exhibited the highest centrality index in the TME network, with significant connections to many cell subtypes (Fig. 6a, S8a) including other fibroblasts, vasculature-associated endothelial cells, and immune cells (Fig. 6b, 6e). Many genes upregulated in gastric cancer participated in the F13-CTHRC1 as well as in endothelial cell-specific interactions, while none of the downregulated genes were interacting, further pointing towards the existence of a tumor-specific cell communication program (Fig. 6d). Notably, F13-CTHRC1 was predicted to communicate with the tumor-enriched EN10-SERPINE1 and EN03-VWA1 endothelial cells through WNT5A-MCAM and LRP1-SERPINE1 interactions, both of which are highly expressed by these cell types and have been characterized as molecular switches for cell motility and angiogenesis [52, 53]. Moreover, by engaging endothelial cells and macrophages via integrin signaling, F13-CTHRC1 might not only foster tumor remodeling through angiogenic stimulation but might also enhance the development of anti-inflammatory macrophages [54] (Fig. 6c). Finally, F13-CTHRC1 harbors immunosuppressive interactions with T cells and DCs through TSLP-IL7R signaling, respectively (Additional file 1: Fig. S9c, S9d). TSLP is responsible for the expansion of regulatory T cells [55] and initiates T helper 2 responses in the tumor which are associated with worse survival prognosis [56]. Together, these findings indicate a potential role for F13-CTHRC1 as driver of gastric cancer progression through angiogenesis-stimulating and immune suppressive interactions with endothelial cells and immune cells.
Top ligand-receptor interactions between gastric tumor cells and the TME point to mechanisms of immune evasion
Although we applied EPCAM+ cell depletion before single-cell sequencing to enrich for cells of the TME, our dataset also contains 7210 epithelial cells. These cells clustered into fifteen groups (Additional file 1: Fig. S10a) which were annotated as gland mucus cells (GMCs) (MUC6 and TFF2), pit mucus cells (PMCs) (MUC5AC and TFF1), parietal cells (ATP4A and ATP4B), chief cells (PGA3 and PGA4), enterocytes (FABP1 and APOA1), and malignant cells (Additional file 1: Fig. S10a-b). The latter were identified via the copyKAT algorithm (Additional file 1: Fig. S10d) that quantifies the level of aneuploid chromosomal regions in cells. Malignant cells (designated EP11-MSC-CLDN7) were found only in the tumor samples (Additional file 1: Fig. S10a-b) and expressed the highest levels of EPCAM (Additional file 1: Fig. S10e) although more sporadic expression of this marker was also seen in the other epithelial cell clusters including enterocytes and GMCs, indicating that our EPCAM cell depletion has been incomplete. High EPCAM levels in EP11-CLDN7 correlated positively with the presence of cell cycle markers in this cluster further confirming the malignant nature of these cells. (Additional file 1: Fig. S10e).
Malignant cell cluster EP11-CLDN7 was connected to 33 cell types in the cell-interaction network (Fig. 6A). The top interacting cell types were thereby CD8 + T cells as well as neutrophils and F13-CTHRC1 fibroblasts (Additional file 1: Fig. S10f). The interaction between EP11-CLDN7 and CD8+ T cells was dominated by the interactions of LGALS9—HAVCR2 and HVEM—CD160 both of which are well-known checkpoints that inhibit CD8 + T cell functionality (Additional file 1: Fig. S10g), suggesting an immunosuppressive function of the gastric tumor cells. EP11-CLDN7 and F13-CTHRC1 cells are predicted to interact through FZD5–WNT5A, reinforcing the hypothesis of F13-CTHRC1 potentially facilitating tumor cell migration [57]. Lastly, neutrophils were connected to EP11-CLDN7 cells via chemokine-receptor pairs CXCL1/2/3/5/8–CXCR1/2 which could be an important axis for the attraction of immune suppressive neutrophils [58, 59].
Cell type signatures predict patient outcome in an independent gastric cancer cohort
Our precise mapping of diverse non-malignant cells in the gastric cancer microenvironment provided a comprehensive reference for interrogating the impact of different cell populations on patient prognosis from other types of data, especially the bulk transcriptome data from the TCGA gastric cancer cohort. Using transcriptomic information from our single-cell data as reference, we used CibersortX [57] to infer the fraction of each cell subtype in the TCGA samples and subsequently assessed the survival differences between patients with higher and lower fractions of each cell subtype. Notably, patients with higher fraction of tumor-associated cell subtypes F13-CTHRC1 and EN10-SERPINE1 had significantly reduced survival time (Fig. 7a, b). In addition, several other endothelial and epithelial subtypes also displayed significant negative impact on gastric cancer patient survival (Additional file 1: Fig. S8), pointing towards being potentially relevant in gastric cancer progression. In contrast, the fraction of cDC1 dendritic cells (M16-cDC-CLEC9A) had significant positive impact on patient survival (Fig. 7c). Although being a rather rare population in the tumor microenvironment, intratumoral cDC1 cells are considered critical for antitumor immunity as they attract, stimulate, and support tumor-infiltrating cytotoxic T cells through different molecular signals [60,61,62]. Consistent with this functional role [63], in our network analysis we found M16-CLEC9A to interact with a wide collection of CD8+ T cell and NK cell populations through a number of communication signals including most prominently the XCL1/XCL2-XCR1 receptor-ligand axis (Fig. 7d). Accordingly, combining these inversely correlated signatures resulted in an even more pronounced survival difference between F13-enriched/M16-depleted and F13-depleted/M16-enriched patients (Fig. 7e).
Integration of gastric cancer immunotherapy bulk RNA-seq reveals cellular and molecular predispositions to response
Immunotherapy via blocking the PD1/PDL1 axis or via blocking CTLA4 signaling has become standard of care in several cancer indications and has yielded promising results in gastric cancer as well [64, 65]. However, many gastric cancer patients do not profit from ICI treatment [66]. To identify what cell types and TME signaling pathways might contribute to resistance to ICI treatment, we obtained bulk RNA-seq data from a recent gastric cancer study [6] that evaluated gene expression in tumor samples from 45 patients treated with anti-PD-1 (pembrolizumab) therapy. Differential gene expression analysis yielded 696 and 1382 genes significantly higher expressed in pretreated samples from the 12 responding and 33 non-responding patients, respectively (Fig. 8a).
Bi-clustering the set of 696 genes associated with response across our single-cell expression data yielded 8 gene groups, designated R1 to R8 (Fig. 8b, S12a). Cluster R1 was exclusively expressed in T cells and NK cells and was strongly enriched for genes associated with IFN-γ signaling [67]. Genes associated with peptide presentation via MHC I complexes (R7) [68] and MHC II complexes (R8) [69] were also found to be overexpressed in responders and specifically associated with lymphoid cells and myeloid cells, respectively (Additional file 1: Fig. S12b). Interestingly, cluster R8 contained prominent marker genes such as C1QA and C1QB and was most strongly expressed in M11-SPP1 anti-inflammatory macrophages, suggesting a positive influence of these cells on immunotherapy response (Fig. 8c).
Bi-clustering of the 1382 non-responder genes in our single-cell data yielded 13 different clusters, designated NR1 to NR13, with high specificity to endothelial cells and fibroblast as well as epithelial and plasma cells (Fig. 8b, S12c). Gene set enrichment analysis identified as main associated pathways angiogenesis in endothelial cells, and IL6 signaling in myeloid cells (Additional file 1: Fig. S12d). Fibroblast have been linked to immunotherapy resistance in several cancer types [70, 71]. Correspondingly, here we found four different gene clusters, NR10, NR11, NR12, and NR13 highly specific for pericytes, resting fibroblasts, F13-CTHRC1, and activated fibroblast, respectively (Fig. 8d). Cluster NR12 was specifically enriched for pathways involving extracellular matrix degradation and expressed several metalloproteinases including MMP9 which has been shown to mediate anti-PD1 resistance in melanoma.
To understand the distribution of these gene clusters in gastric cancer, we calculated a signature score for all eight responder and 13 non-responder clusters for each of the 295 gastric cancer RNA-seq samples from the gastric cancer TCGA cohort. Hierarchical bi-clustering of the scores clearly separated response and non-response gene signatures and revealed differential expression of each signature across different patient groups (Additional file 1: Fig. S12e). Non-responder clusters NR3, NR5, and NR8 were thereby found to be clearly distinct from the other non-responder gene clusters likely due to their association with plasma and epithelial cells. Classifying the TCGA tumors based on recently identified molecular subtypes revealed that responder signatures were enriched in ICI sensitive EBV and MSI patients while non-response signatures were enriched in CIN and GS patients (Additional file 1: Fig. S12f).
To investigate if any of the 11 main cell types or 81 subtypes identified in the gastric TME could be used to predict response to immunotherapy, we analyzed the power of a binary classification model consisting of the top marker genes from each cell type (see “Methods”) to predict treatment outcome in the above gastric cancer immunotherapy cohort [6]. With an area under the ROC curve of 0.85, the presence of T27-cycling CD8+ T cells followed by the presence of pan CD8 T cells (AUROC of 0.81) was found to best predict treatment outcome in this cohort (Fig. 8e). Presence of F14-ADAM28 activated fibroblasts followed by presence of pan-fibroblasts was best predictive for non-response with AUROCs of 0.82 and 0.7, respectively (Fig. 8e) and up to 100% specificity (Additional file 1: Fig. S12h-i and Additional file 7: Table S6).