Single-Cell Genomic Characterization Reveals the Cellular Reprogramming of the Gastric Tumor Microenvironment

Gastric cancer is the fifth most common cancer and the third leading cause of cancer-related deaths worldwide (1). The current histopathologic classification scheme designates gastric cancers as either intestinal or diffuse according to the morphology, differentiation, and cohesiveness of glandular cells. Intestinal gastric cancer is preceded by changes in the gastric mucosa called the Correa cascade that progresses through inflammation, metaplasia, dysplasia, and adenocarcinoma (2). Diffuse gastric cancers lack intercellular adhesion and exhibit a diffuse invasive growth pattern. Recent integrated genomic and proteomic analyses including by The Cancer Genome Atlas (TCGA) and the Asian Cancer Research Group (ACRG) have refined the classification of gastric cancer into distinct molecular subtypes that include the intestinal and diffuse classifications (3, 4). Regardless of the histopathologic or molecular subtype, gastric cancers are not isolated masses of cancer epithelial cells. Rather, these tumors have a complex morphology where cancer cells are surrounded by the tumor microenvironment (TME), a cellular milieu containing diverse cell types such as fibroblasts, endothelial cells, and immune cells.

Increasingly, it is recognized that the cellular features of the TME play an important role in enabling tumors to proliferate and metastasize. A major component of the TME that influences tumor cell survival as well as response to treatments such as immune checkpoint blockade is the diverse and deregulated cellular states of the immune cells (5). Thus, the cellular characterization of the TME provides a more sophisticated picture of the context of tumor cell growth within its tissue of origin, characteristics of immune infiltrate and intercellular interactions.

The major objective of this study was to determine the specific cellular and transcriptional features that distinguish the gastric cancer TME from normal gastric tissue. We sought to define these differences at the resolution of single cells with single-cell RNA-seq (scRNA-seq). We delineated cell-specific features that are otherwise lost when using “bulk” methods in which molecular analytes cannot be attributed to their cell-of-origin. We accomplished this by using an extensive analytic framework (Fig. 1A; refs. 6–9) that revealed changes in transcriptional states, regulatory networks, and intercellular communication between matched gastric tumor and normal tissue from the same patients, together with peripheral blood mononuclear cells (PBMC) from a subset of patients. Our study identified cellular and biological features that are specific to the TME and thus offer insights that may help infer new therapeutic targets.

All samples were acquired with informed consent under an approved Institutional Review Board protocol from Stanford University (Stanford, CA) as surgical resections or endoscopic biopsies. Matched normal tissue was obtained from sites displaced at least several centimeters from the tumor and was confirmed to lack tumor cells on histopathology review.

Tissues were collected immediately after resection and dissociated into a single-cell suspension (Supplementary Methods). PBMCs were isolated using density gradient centrifugation. Single-cell libraries were generated from cell suspensions using Chromium Single Cell 3′ Library & Gel Bead Kit v2 (10x Genomics) as per the manufacturer's protocol and sequenced on Illumina sequencers (Illumina; Supplementary Table S1). Cell Ranger v 3.0 (10x Genomics) “mkfastq” and “count” commands were used with default parameters and alignment to GRCh38 to generate matrix of unique molecular identifier (UMI) counts per gene and associated cell barcode. Datasets are available under dbGAP identifier phs001818.v1.p1. Cell Ranger outputs are available at https://dna-discovery.stanford.edu.

We used Seurat (v2.3.4; ref. 9) to create data objects from the matrix outputs. We removed cells that expressed fewer than 200 genes, had greater than 20% mitochondrial genes, or had number of UMI in an outlier range indicative of potential doublets (Supplementary Table S1). We excluded genes detected in fewer than three cells. Data were normalized to log scale using the “NormalizeData” function with a default scale parameter of 10,000. Highly variable genes were identified using the “FindVariableGenes” function with parameters for x.low.cutoff = 0.0125, x.high.cutoff = 6 and y.cutoff = 0.5. The effects of variation in sequencing depth were regressed out by including “nUMI” as a parameter in the “ScaleData” function. These variable genes were used as input for PCA using the “RunPCA” function. The first 20 principal components (PC) and a resolution of 0.8 were used for clustering using “FindClusters.” UMAP was used for two-dimensional representation of first 20 PCs with “RunUMAP.”

Differential gene expression for identifying markers of a cluster relative to all other clusters or compared with a specific cluster was determined using the “FindAllMarkers” or “FindMarkers” functions, respectively. Parameters provided for these functions were: genes detected in at least 25% cells, and differential expression threshold of 0.25 log fold change using Wilcoxon rank sum test with P < 0.05 following Bonferroni correction. We compared the marker genes for each cluster to literature-based markers of cell lineages to assign a cell lineage per cluster (Supplementary Table S2).

Individual Seurat data objects were merged iteratively using the “MergeSeurat” function after filtering doublets identified by DoubletFinder, an R package that enables computational identification of doublets (Supplementary Methods; ref. 10). The merged object was processed as described above with library preparation batch and number of UMIs (Supplementary Table S1) included as parameters for regression in the “ScaleData” function to regress batch effects and variation in sequencing depth, respectively. The “DoHeatmap,” “FeaturePlot,” “DimPlot,” “DotPlot,” and “VlnPlot” functions were used for visualization.

For a secondary cluster analysis of each cell lineage from this aggregated dataset, clusters of interest were identified and subset using “SubsetData” with parameter “do.clean” set to true. Detection of variable genes, scaling with UMI regression, PCA, clustering, and UMAP were repeated as described above. Following this step, we removed clusters with coexpression of cell lineage markers as multiplets (Supplementary Table S3). Proportions of each cell type relative to the total number of cells in the sample were compared for tumor and normal sites using two proportions Z-test and represented as box plots after reclustering analysis for each lineage.

A complete description of methods is available under Supplementary Methods.

We obtained tissues from surgical resections or endoscopic biopsies from seven patients with gastric cancer and one patient with gastrointestinal metaplasia (GIM). These tissue samples represented paired gastric tumor and normal tissue from the same patient derived from the same anatomic region. Four tumors were located at the gastroesophageal junction (GEJ) and four in the body and antrum. We analyzed matched PBMCs from two patients. On the basis of histopathology review, the gastric cancer tumors had intestinal, diffuse, or mixed features (Supplementary Table S4; Fig. 1B; Supplementary Fig. S1A). We classified tumors into microsatellite stability (MSS) or microsatellite instability (MSI) molecular subtypes based on expression of the DNA mismatch repair proteins MLH1, MSH2, MSH6, and PMS2 according to IHC (Supplementary Table S4; Supplementary Methods). Histopathology showed that three patients (P5931, P6207, P6342) had active gastritis or intestinal metaplasia in paired nonmalignant tissue.

With scRNA-seq, we obtained transcriptional profiles of 32,407 single cells from tumors or metaplasia, 18,657 single cells from paired normal tissue, and 5,103 PBMCs (Supplementary Table S1). To determine gene expression changes, we employed a series of analytic steps for each individual dataset including quality filtering, principal component analysis (PCA), and graph-based clustering (Supplementary Methods; refs. 9, 11). We employed uniform manifold approximation and projection (UMAP) to reduce the dimensionality of this data and allow the visualization of cell-type clusters defined by their transcriptional profiles.

Differentially expressed (DE) genes were identified as genes expressed in greater than 25% of cells in a cluster and having a log fold change greater than 0.25, using a cutoff of P < 0.05 following Bonferroni correction. We compared the DE genes from each cluster to known marker genes of various cell types (Supplementary Table S2). This information enabled us to link clusters to specific cell types including epithelial cells (expressing PGC, TFF1, MUC5AC, EPCAM, GIF, CHGA), fibroblasts (THY1, DCN, COL4A1, FAP), endothelial cells (PECAM, ENG, VWF, SELE), immune cells (PTPRC) such as CD4 T (CD3D, CD4, IL7R), cytotoxic T (CD3D, CD8A, CD8B), regulatory T (FOXP3, IL2RA), NK (NKG7, GNLY), B (MS4A1), plasma cells (immunoglobulin genes), mast cells (TPSAB1), and macrophages (CD68, CD14, FCGR3A). Examination of cells expressing markers of disparate cell types facilitated the computational detection of doublets (Fig. 1C and D; Supplementary Methods).

To determine the similarities and differences among all of the samples and patients, we aggregated tumor, metaplastic, normal sites, and PBMCs scRNA-Seq data into a single data matrix. To reduce experimental variance, we regressed out batch effects from library preparation and variation in sequencing depth. Also, we removed computationally detected doublets (Supplementary Methods).

Clustering analysis identified 40 distinct clusters that were specific for a variety of cell types including epithelial, stromal (fibroblasts, endothelial cells), lymphocytes, macrophages, and mast cells; they were found among all of the patient samples (Supplementary Figs. S1B and S1C; Fig. 1E). On closer examination, each cell type had multiple distinct transcriptional states. For additional characterization, we aggregated data for each cell type across samples and conducted a clustering analysis (Supplementary Methods).

We detected differences between tumor, metaplastic, and normal epithelial cells; differences in tumor epithelium derived from different patients; as well as subclonal heterogeneity within an individual tumor. Reclustering analysis of epithelial cells from our integrated dataset revealed three subclasses. The first subclass consisted of normal gastric epithelial cells; over 80% were derived from normal gastric samples (Fig. 2A–C). Normal epithelial cells were detected in all samples regardless of their origin from tumor, normal, or metaplastic tissue (Supplementary Fig. S2C; Supplementary Table S5).

The second subclass consisted of tumor-specific epithelial cells, which we label as the “GC type 1.” Approximately 98% of these cells originated from tumor samples. Interestingly, each cluster in this subclass was dominated by a single patient (Supplementary Fig. S2C; Supplementary Table S5); this indicated the extent of intertumor heterogeneity among all of the gastric cancers, meaning that each individual tumor had distinct transcriptional properties.

The third subclass involved epithelial cells derived from gastric cancer as well as normal tissue, which we label as “GC type 2.” It also contained around 58% of cells from patient P6649 who had gastric metaplasia but did not have gastric cancer (Supplementary Fig. S2C). We detected a range covering 11% to 68% of cells from normal tissues of patients with background metaplasia or gastritis on histology (P5931, P6207, P6342; Supplementary Table S4; Supplementary Fig. S2C; Supplementary Table S5) in this subclass. Given their origin, these cells are likely to represent metaplastic or dysplastic epithelium, which had transcriptional features that overlapped with tumor epithelium.

As an independent determination of gastric tumor versus normal epithelial cells, we employed a different method called scPred that uses a supervised machine-learning algorithm and thus provides highly accurate cell-type assignment (12). We randomly subset tumor and normal site–derived cells into two. One subset from each of these classes was included in the training dataset to build the scPred prediction model. We tested this model on all of the remaining cells for all samples. When we compared the scPred results to our Seurat analysis, the results were concordant. 89% and 77% of cells in the “GC type 1” and “GC type 2” subclasses were not classified as normal gastric epithelial cells, indicating that machine learning confirms the distinction between normal, gastric cancer type I, and gastric cancer type II cells (Fig. 2D; Supplementary Fig. S2D). Thus, scPred confirmed the transcriptional differences among the epithelial subclasses identified by our clustering analysis.

Differential expression analysis distinguished the three subclasses of epithelial cells. In normal gastric epithelium, we identified distinct mucosal populations such as pit, mucous neck, zymogen secreting chief, intrinsic factor producing parietal, and neuroendocrine cells (Supplementary Fig. S2B). In contrast, gastric cancer type I and gastric cancer type II subclasses downregulated some of these gastric mucosa marker genes such as MUC6, TFF2, TFF1, and MUC5AC (Supplementary Table S6). The gastric cancer type I subclass had increased expression of intestinal mucosa markers TFF3, FABP1, SPINK4, MUC13, and REG4. The gastric cancer type II subclass had significantly increased expression of previously identified gastric cancer marker genes KRT7 and KRT17 (3), SOX4, and HES1 that have been implicated in metaplasia pathogenesis (13). Compared with normal cells, both gastric cancer type I and type II had upregulation of gene sets that included pathways for Myc, DNA repair, and Notch signaling (Supplementary Fig. S2E). However, only gastric cancer type I cells had higher upregulation of EMT and KRAS signaling.

Using scRNA-seq, one can detect amplifications or deletions at chromosome arm level by analyzing concomitant increase or decrease in gene expression, respectively (Supplementary Methods; ref. 14). To identify them, we analyzed each patient's normal epithelial subclass cells against those from matched gastric cancer type I or gastric cancer type II subclass (Fig. 2E; Supplementary Fig. S3A). Copy number changes were inferred according to the posterior probability for each cell to belong to one of the two components with lower or higher gene expression indicative of deletion or amplification, respectively.

A diverse spectrum of chromosome arm imbalances was present in these tumors. For patient P6207, we detected amplifications in 7p, 7q, 8q, and 9q with a deletion in 10p. Patient P6342′s tumor had amplification of 20q with deletion of 4q. Patient P5846′s tumor had an amplification of 19q. Patient P5866′s tumor had deletion of 16q. Patient P5931′s tumor also had amplifications of 7p and 7q, which we have successfully identified using single-cell DNA sequencing (15). For patient P6649 with metaplasia and P6592, samples contained only a small number of cells with significant copy number changes. We excluded the sample from patient P6709 from this analysis because we only detected 21 epithelial cells from the tumor site, possibly indicative of response to neoadjuvant chemotherapy (Supplementary Figs. S2B and S1A; Supplementary Table S4). These CNA results provide additional orthogonal confirmation of the distinction between normal and tumor epithelial cells in our gastric cancer dataset.

We analyzed gastric cancer type I and gastric cancer type II cells across patients for the transcriptional activation of various oncogenic pathways. We observed significant differences in activation levels reaffirming the transcriptional heterogeneity between patients (Fig. 2F). Interestingly, these activation profiles did not cluster according to differences between the molecular subtypes of MSI and MSS.

The individual patient tumors contained multiple clusters of gastric cancer type I and gastric cancer type II cells, an indication of intratumor subclonal heterogeneity. Pathway activation analysis grouped these clusters into three to five subpopulations for each patient's tumor confirming a subclonal structure (Fig. 2G; Supplementary Fig. S4). For example, heterogeneity in P6207 was characterized by differences in cell cycle, KRAS pathway, and Wnt activation (Supplementary Fig. S2G). We postulate that these subpopulations may provide a distinct growth advantage to the tumor.

TME macrophages had distinct cellular and gene expression changes compared with normal gastric tissue consistent with the observation that macrophages acquire heterogenous phenotypes depending on their activating stimulus (16). Macrophage phenotypes are called M1 or M2 with anti- and protumorigenic functions, respectively. However, the gene expression signatures we observed did not fall in line with either the canonical M1 or M2 classifications.

Specifically, we detected various subclasses of myeloid lineage cells as seen with 11 distinct clusters across all patients (Fig. 3A–D; Supplementary Fig. S5A; Supplementary Table S7). The nine monocyte–macrophage clusters were defined by marker gene expression of CD14, FCGR3A, and CD68. The two dendritic cell (DC) clusters were defined by marker gene expression of CLEC4C, ID2, IRF4, and CD83 (Fig. 3C). Notably, macrophages were significantly enriched in tumor compared with normal gastric tissue (Fig. 3E).

We examined the expression of marker genes for M1 (e.g., CCL19, TNF, CCL5) and M2 (e.g., MRC1, CCL18, CCL13, CD163) states across the macrophage clusters (Fig. 3F; ref. 16). Expression of M1/M2 genes did not distinguish the clusters. Moreover, these genes were coexpressed in the same cluster. This suggests that the transcriptional heterogeneity was independent of the M1/M2 classification. We identified the significantly DE genes across all clusters to assess heterogeneous phenotypes (Supplementary Table S8). This revealed that heterogeneity was related to differences in the expression of HSP family genes; THBS1; chemokines including CCL20, CCL18, and CCL3; matrix metalloproteinase genes; complement family; and cell-cycle regulation genes (Fig. 3G). Clusters also showed differential enrichment of hallmark gene set activity confirming their distinct transcriptional programs (Supplementary Fig. S5B).

PBMC monocytes clustered distinctly from tumor or normal macrophages indicating transcriptional differences (Fig. 3B and D). As it stands, conventional single-cell clustering methods do not delineate the dynamic process of cell differentiation. A new type of method called trajectory analysis uses single-cell gene expression to determine the transition among specific cell-type lineages and states. Our trajectory analysis of tissue macrophages and PBMC monocytes yielded three different trajectory states (Fig. 3H). Monocytes were present in a single trajectory state with few tissue macrophages. The majority of tissue macrophages differentiated along two distinct states. This suggests that tumor-infiltrating macrophages differentiate from monocytes but retain some fundamental similarities to macrophages within normal tissue.

DCs had two subclasses. One had genes that define plasmacytoid DCs including IL3RA (CD123) and CLEC4C (CD303); this subclass was detected predominantly in PBMCs (Fig. 3C and D). A second was enriched in the TME and showed significant DE of activated DC gene markers CD83, CCR7, IL7R, and ID2 (Supplementary Table S9; Supplementary Fig. 5C; refs. 17, 18). This subclass expressed the chemokines CCL22, CCL17, CCL19, and IL32, which are associated with recruitment of naïve T cells. Moreover, it highly expressed IDO1, characteristic of an immunosuppressive phenotype (19, 20). This result represented a novel gene expression program in TME infiltrating DCs not previously described.

We compared activity levels of 1,558 experimentally derived immunologic gene signatures containing the term “macrophage,” “DC,” or “monocyte” (Supplementary Methods; Supplementary Table S10) to these gene expression profiles. The identity of monocyte and DCs were confirmed using this approach. Also, these results provided gene set information indicating the activation phenotype of tumor-specific DCs. Each macrophage cluster was enriched for gene sets derived from a variety of experimental conditions. Hence, macrophage heterogeneity likely reflects stimulus-based context within the TME.

Regulatory genes controlling expression of a group of genes are referred to as regulons. We identified regulons for these different transcriptional cell states (6). This analysis identified transcriptional regulators such as IRF4 in DCs and also revealed a distinct set of regulatory genes defining the various macrophage populations including NFKB1, ETS2, CREM, REL, STAT1, FOXO3, etc. (Supplementary Table S11). Our data provide direct in vivo evidence that tumor-specific macrophages exist in a continuum of stimulus-dependent functional states regulated by a specific set of genes rather than the M1/M2 paradigm.

Exhausted T cells (TEx) were a prominent feature of the gastric TME compared with normal tissue. The TME was also significantly enriched for Treg cells compared with normal. Initially, we identified CD4 helper T, CD8 T, NK, Treg, plasma, and B cells using classic immunophenotyping markers. We filtered clusters with high expression of HSP family genes and lacking lineage markers (Supplementary Fig. S6A–S6D; Supplementary Methods).

Cytotoxic CD8 T lymphocytes (CTL) were distributed across five different clusters indicative of transcriptionally distinct cell states (Supplementary Fig. S7A). They were observed among all patients with the majority detected in P5866, P5931, and P6207 (Supplementary Table S12). We examined each cluster for the dominant sample origin (normal, tumor, or PBMCs), expression of naïve markers (CCR7, SELL), tissue effector memory markers (CD69, ITGAE, ITGA1), and cytotoxic genes (GZMB, GZMA, PRF1, IFNG, NKG7; Fig. 4A and B). CTL subclasses included naïve tumor CTLs (cluster 0), effector normal CTLs (cluster 1), effector PBMC CTLs (cluster 19), and two subclasses of tumor effector CTLs (clusters 6, 22). Differential expression analysis identified distinct signatures among these subclasses (Supplementary Table S13).

Within cancer, one observes TEx cells with reduced cytotoxic activity and expression of inhibitory receptors (21). In the TME, both CTL subclasses had low expression of PRF1 and IFNG while expressing granzymes (Fig. 4B), indicative of a TEx phenotype (22). Furthermore, they expressed multiple immune checkpoints [PDCD1, TIGIT, CD96, HAVCR2 (TIM3), LAG3, CTLA4, VSIR]; the expressions levels of these checkpoints is an indicator of the degree of immunosuppression. These cells also expressed multiple co-stimulatory molecules [ICOS, TNFRSF18 (GITR), CD27, TNFRSF9, CD40, CD226, TNFRSF4; Fig. 4B].

We compared the transcriptional profiles of tumor, normal, and PBMC CTLs to previously published gene signatures to further understand the tumor TEx phenotype. First, we used gene signatures from T-cell exposure to viruses. This experimental setup distinguishes the transcriptional program of naïve cells unexposed to virus, effector cells responding to the lymphocytic choriomeningitis virus (LMCV), and TEx cells responding to a specific clone 13 of LCMV (22). This result confirmed our assignment of naïve, effector, and TEx classes across the five clusters. The two tumor-specific TEx clusters had greater enrichment for the exhaustion signatures generated with LCMV clone 13 (Supplementary Fig. S7B). Also, we compared the CTL transcriptional profiles to three independently derived TEx profiles from single-cell analysis of human tumors, one from mouse tumors, together with signatures for cytotoxicity and proliferation (Fig. 4C; refs. 23–25). This result confirmed the low-cytotoxicity, high-exhaustion phenotype in both subclasses of tumor-TEx cells. Interestingly, these TEx subclasses were derived mainly from tumor samples from P5931, P6207, P6592, and P6709 that received neoadjuvant therapy (Supplementary Fig. S7A).

The two subclasses of tumor-TEx cells differed in the extent of their exhaustion and proliferation gene expression program (Fig. 4C). The subclass with higher proliferation (cluster 22) was associated with greater exhaustion indicative of an active immune response, which has previously been associated with terminal exhaustion (26). Immune checkpoint or co-stimulatory molecule gene expression was not significantly different between the two. The second subclass (cluster 6) had lower exhaustion and proliferation plus significantly higher expression of CXCL13 as previously identified in tumor but not in viral TEx (Supplementary Fig. S7C; ref. 27). These cells expressed high RBPJ, NR3C1, and BATF that are regulators of CD8 T-cell fate (28). Our results thus demonstrated that effector CTLs in the TME are exhausted unlike normal tissue or PBMCs with two distinct subclasses characterized by high CXCL13 expression or proliferation.

To verify our findings, we conducted multiplex immunofluorescence staining for pan-cytokeratin/SOX-10 (epithelial cells), CD45RO (memory T cells), CD3 (T cells), and PD-1 (TEx) in tumors from four patients where adequate tissue samples were available (Fig. 4D; Supplementary Fig. S8A). We detected these cell lineages in all samples including cellular subpopulations of effector T cells (CD45RO, CD3 positive) and exhausted effector T cells (PD-1, CD45RO, CD3 positive) based on coexpression analysis (Supplementary Fig. S8B). The stromal cells and macrophages were also apparent as they lacked expression of these markers.

Tregs were significantly enriched in the gastric TME compared with normal gastric tissue, thus indicating an important mode of immunosuppression (Supplementary Fig. S9A and S9B). Majority of cells were detected in P5866, P5931, P6207, and P6592 (Supplementary Table S12). The Treg cells had two distinct subclasses distinguished by higher expression of markers of proliferation (eg. MKI67, TYMS) in one subclass (cluster 22). In addition, Tregs expressed several immune checkpoint and co-stimulatory molecules representing potential targets to modulate their function (Supplementary Fig. S9C).

CD4 T cells were represented by four subclasses, mainly represented in P5866, P5931, and P6207 (Supplementary Table S12). Three of these were characterized by expression of naïve markers (CCR7, SELL) originating from the PBMCs or normal gastric tissue (Supplementary Fig. S9D and S9E). Effector CD4 T cells (cluster 6) lacking naïve marker gene expression were found in both normal and tumor tissue (Supplementary Fig. S9D and S9E). They had significantly higher expression of GZMA, GZMB, CXCL13, BATF, and HLA genes (Supplementary Fig. S9E). These cells are likely to represent follicular helper-like CD4 cells that express CXCL13 and are associated with tertiary lymphoid structures (29, 30).

Three subclasses of NK cells were detected in PBMCs (cluster 2) and both tumor and normal tissue (clusters 13, 14) (Supplementary Fig. S10A). They were derived prominently from P5866, P5931, and P6207 (Supplementary Table S12) and also contained a mix of rare populations of invariant NK and innate lymphoid cells (Supplementary Fig. S10B; ref. 31). Cells in tumor and normal sites clustered together indicating transcriptional similarity. Tumor-site cells expressed cytotoxic molecules such as GZMA, XCL2, CCL5, PRF1, CCL3, and CCL4, indicating potential for an antitumor response (Supplementary Fig. S10C). Cells also expressed several inhibitory and co-stimulatory molecules including TNFRSF18 (GITR), CD96, and KIR2DL4 expression representing targets for modulating their function.

B cells from gastric TME and normal sites clustered together indicating transcriptional similarity (Supplementary Fig. S10D). However, plasma cell clusters showed significant differences in the expression of genes encoding immunoglobulin isotypes with increased IgA encoding genes in normal tissue and IgG in TME (Supplementary Fig. S10E and S10F).

Next, we analyzed the regulatory genes or regulons that control these lymphocyte subpopulations (Supplementary Table S14; ref. 6). The CXCL13-high tumor-TEx cells had significant enrichment of FOXO1 activity that is required for post-antigen expansion of CD8 T cells (32). In CXCL13-high CD4 T cells, the BATF gene network activity was prominently high validating their similarity to recently described CXCL13-producing helper T cells (33). We successfully identified FOXP3 and BATF gene network enrichment in Tregs confirming the accuracy of this approach. We discovered additional transcriptional regulators of TME Treg fate including KDM5B, MAF, IKZF2, SOX4, BCL3, etc. These regulons are of potential translational value given the interest in targeting epigenetics for immunotherapy (34).

We discovered transcriptional reprogramming of stromal cells in the tumor compared to normal that allows the generation of a tumor-specific extracellular matrix (ECM). Our analysis of the stromal cells across all samples identified fibroblasts (THY1, DCN, COL4A1, FAP), endothelial cells (PECAM, ENG, VWF, SELE), and pericytes (RSG5, PDGFRB; Fig. 5A and C). The observed heterogeneity among fibroblasts is likely to be driven by patient-specific factors (Supplementary Fig. S11A; Supplementary Table S15). Endothelial clusters additionally differed in the expression of genes encoding secretory factors (ESM1, ANGPT2), tip cell markers (COL4A1, COL4A2, DLL4, MARCKSL1), and stalk cell markers (ACKR1, CD36, SELP, VWF; ref. 35) and were mainly represented in P5866, P6342, P6592, and P6709 (Supplementary Fig. S11B and S11C; Supplementary Table S15).

All three cell types were enriched in tumor tissue compared with normal (Fig. 5B and D). Stromal cells are responsible for the production and maintenance of ECM that provides mechanical support to cells and also influences their growth. Genes encoding for components or regulators of the ECM have previously been identified as the “matrisome” (36, 37). This gene group consists of core factors (collagens, proteoglycans, and ECM glycoproteins) that make up the ECM and an associated program (ECM regulators, secretory factors, and ECM-affiliated proteins). To determine the phenotypical differences of stromal cells in normal or tumor tissue, we compared their DE genes to the matrisome gene expression program. Tumor-specific fibroblasts, pericytes, and endothelial cells expressed diverse ECM core and associated components (Fig. 5E; Supplementary Table S16). In addition, fibroblasts in tumors had significant overexpression of ACTA2 compared with normal tissue, indicative of their contractile ability (Supplementary Fig. S11D).

Stromal cells at tumor or normal sites had significantly different regulatory genes or regulons (Supplementary Table S17). For example, tumor-specific endothelial cells had greater activity of SOX18 and SOX7, which are known regulators of a variety of endothelial cell processes (38). Tumor-specific fibroblasts had high activity of the EGR2 gene that can influence fibrosis (39). Tumor-specific pericytes were enriched for FOXF2 activity that is known to regulate pericyte differentiation (40). Hence, our approach distinguished differences in both the gene expression program and its regulators between tumor and normal stromal cells.

We discovered a TME-specific intercellular communications network that can potentially affect cellular behavior. First, we identified significant receptor-ligand interactions between different cell types using CellPhoneDB (8). Then, we compared these networks between tumor and normal tissue to generate a TME-specific interactome (Fig. 6A; Supplementary Table S18).

Stromal cells were among the most prolific interactors. Prominent communication between them and epithelial cells occurred through various integrin receptor interactions with collagen, fibronectin, and THBS1 ligands (Fig. 6B). Bidirectional interactions between ephrin receptor family and ligands were detected in epithelial and endothelial cells. Among growth factor signaling that can promote cancer cell proliferation, we detected EGFR and MET receptors on epithelial cells together with respective ligands on stromal cells. We also detected significant EGFR interactions in fibroblasts.

Fibroblasts were a prominent source of Wnt ligands with expression of corresponding receptors on tumor epithelium, endothelial cells, fibroblasts, and pericytes. This included a LGR4–RSPO3 interaction that has the potential to regulate stemness (Fig. 6C), validating our previous discovery of fibroblast-derived RSPO3 in an organoid model of gastric cancer (41).

Autocrine and paracrine Notch signaling, a known regulator of angiogenesis (42), was evident in endothelial cells. Interactions promoting Notch signaling were also significant in fibroblasts (Fig. 6D). Angiogenic receptors KDR, FLT1, FLT4, PDGFB, and TEK on endothelial cells and pericytes had significant autocrine and paracrine interactions with their respective ligands (Fig. 6E). Among the interactome were 19 cytokines including chemokines, interleukins, TNFs, and their corresponding receptors that can influence immune cell fates.

For this study, we leveraged paired distal normal tissue and PBMCs to analyze the cellular dysregulation and biological changes in the gastric cancer TME. With single-cell gene expression analysis, we demonstrated that gastric cancer TME leads to a series of dramatic cellular changes compared with matched normal stomach mucosa. Specifically, we noted increases in cell numbers of stromal cells and Tregs in the TME. We also identified transcriptional cell states unique to the TME including in DCs and two subclasses of exhausted CTLs. Among our novel findings, we discovered that gene expression profiles for gastric cancer TME macrophages are heterogenous and not confined to a binary M1/M2 designation. We demonstrated that TME stromal cells encode for a specific ECM composition not found in normal tissue. We identified novel gene regulatory networks and TME-specific intercellular communication.

We validated previously described changes in normal, metaplastic, and tumor epithelial cells (13, 43) and were additionally able to elucidate intratumor heterogeneity by examining the activity of various cancer-promoting mechanisms within the tumor cells. The diversity in activation profiles suggests that targeting multiple clones with different combination strategies may be necessary to eradicate them completely.

Immune cell lineages were detected across all patients. However, a limitation of our analysis is that the majority of cells were derived from two patients (P5866 and P5931) with two sequencing replicates (Fig. 1F; Supplementary Table S1). Immunosuppression in gastric cancer TME was evident by the increased proportion of Tregs compared with normal tissue. We identified several checkpoint and co-stimulatory molecules on these cells, and their transcriptional regulators. These regulators can be investigated to understand Treg biology and to derive therapeutic targets. Indeed, recent evidence indicates that anti-CTLA4 activity is a consequence of Treg depletion in the TME (44). We detected expression of multiple immune checkpoints on cytotoxic T cells similar to other studies (23, 45). These checkpoints were also detected on helper T and Treg subsets. Thus, it is important to understand the effects of immune checkpoint blockade on distinct T-cell subpopulations that express the same target. Our analysis revealed transcriptional regulators responsible for these states. Plasma cells in tumor tissue expressed IgGs rather than IgAs that were detected in paired normal tissue that have been associated with a procancer role by influencing myeloid cell Fc receptors (46).

Increased immune cell signaling has previously been associated with EBV molecular subtype of gastric cancer (3). Meanwhile, MSI subtype cancers have been demonstrated to have higher responses to immune checkpoint blockade owing to the favorable immune milieu generated by higher neoantigen burden (47). In our study, immune cell subtypes did not cluster clearly according to the MSI/MSS status of samples. This resembles findings from recent studies in gastric cancer where molecular subtypes could be differentiated by gene signatures based on a combination of immunosuppression, immune activation, and stromal activation processes rather than the leucocyte fractions within the TME (48, 49). This raises the possibility that a combination of transcriptional states in stromal and immune cells, genomic alterations in tumor cells, and their resultant interactions determines the final composition and transcriptional cell states of the TME in each molecular subtype. Utilizing scRNA-seq on a larger cohort of gastric cancer has the potential to address the mechanisms by which tumor genomics influences the TME composition and possibly responses to immunotherapy.

It has been demonstrated that neoadjuvant chemotherapy leads to increased expression of CD4, CD8, PD1, PD-L1, and TIM-3 proteins in gastric cancer TME (50). We observed presence of TEx cells in all patients that received neoadjuvant chemotherapy. Increasing the cohort size will also help to address the influence of neoadjuvant treatment on the TME transcriptional cell states that we have observed in this study.

Interactome analysis demonstrated protumor effects of TME components and also the influence of cancer cells on the TME. Tumor-specific interactome represents potential treatment targets to inhibit cancer proliferation, overcome the immunosuppressive microenvironment, and restore the cancer immunity cycle. In addition, although some targets such as Wnt inhibition have previously been regarded only in the context of tumor epithelial cells, our analysis demonstrates that this might have implications for the TME.

Our study did not consider spatial context and might be affected by the dissociation process. The use of dual single-cell proteomics and transcriptomics is likely to provide a more refined analysis of immune cell subtypes (51).

No potential conflicts of interest were disclosed.

Conception and design: A. Sathe, G. Poultsides, H.P. Ji

Development of methodology: A. Sathe, B.T. Lau, H.P. Ji

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A. Sathe, B.T. Lau, J. Chen, G. Poultsides, H.P. Ji

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): A. Sathe, S.M. Grimes, C. Suarez, H.P. Ji

Writing, review, and/or revision of the manuscript: A. Sathe, S.M. Grimes, R.J. Huang, G. Poultsides, H.P. Ji

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): H.P. Ji

Study supervision: H.P. Ji

We are grateful to all patients who participated in the study. We thank Christine Handy, Alison Almeda, and Christina Wood-Bouwens for assistance in sample collection and documentation and Ann Renschler for assistance in sequencing. We thank Dr. Katir Patel from Ultivue, Inc. for assistance in multiplex immunofluorescence staining. This work was supported by NIH grants P01HG000205 (HPJ, CWB, and SMG), R01HG006137 (HPJ), and U01CA217875 (HPJ and AS). HPJ also received support from the American Cancer Society (124571-RSG-13-297-01), the Clayville Foundation, and the Gastric Cancer Foundation.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.