This study included six patients with gastric adenocarcinoma who were diagnosed through endoscopic biopsy at Seoul National University Hospital from July to August 2021. All patients had advanced gastric cancer, as indicated by preoperative endoscopy and computed tomography scans. Surgical treatment was subsequently undertaken for all patients.

MSI status was assessed at five loci (BAT25, BAT26, D2S123, D5S346, and D17S250) according to the National Institutes of Health guidelines. In brief, tumor and paired normal areas were dissected from formalin-fixed, paraffin-embedded blocks, and genomic DNA was extracted. After polymerase chain reaction (PCR) amplification, the PCR products were run on the ABI 3731 automatic DNA sequencer (Applied Biosystems, Foster City, CA) [14]. GC with instability in two or more, one, or no loci were considered MSI-H, microsatellite instability–low, and MSS, respectively.

Tumor and matched normal tissues were obtained immediately after surgical resection of the stomach to minimize cold ischemic time. Although we collected six tumor samples and paired six normal samples, we selected two normal samples with higher cell viability because of the limitation of budget. The samples were processed for library preparation using Chromium Next GEM Single-Cell 3′ Reagent Kits v3.1 (10× Genomics, Pleasanton, CA). After appropriate cell dissociation and viability assessment (cutoff 80%), cell suspensions were diluted to a targeted count of 10,000. The cells were mixed with master mix and loaded onto a chromium Next GEM chip G with Single Cell 3′ v3.1 Gel Beads and Partitioning Oil. RNA transcripts from single cells were uniquely barcoded and reverse transcribed within droplets. The cDNA molecules were pooled and processed for library construction. The libraries were then quantified using quantitative PCR and qualified using Agilent Technologies 4200 TapeStation (Santa Clara, CA). Finally, the libraries were sequenced using the HiSeq platform (Illumina, San Diego, CA).

Raw gene expression matrix was generated for each sample using CellRanger v2.0.2. Cells that had less than 300 expressed genes or more than 6,000 expressed genes were removed. Cells in which the fraction of mitochondrial genes exceeded 20% also were removed. We removed doublet cells using the DoubletFinder package with PCs 1 to 30, pN set to 0.25, and pK set to 0.092. We identified 3,000 highly variable genes using the FindVariableFeatures function in the Seurat library for initial clustering. The unique molecular identifier (UMI) count per gene were normalized by the total UMI count in each cell and log transformed with the NormalizedData function in Seurat using 10,000 as the scale factor. The effects of the number of detected UMIs, the fraction of mitochondrial genes, and cell cycles on the gene expression values were corrected by regression using the ScaleData function in Seurat.

Before the clustering, we performed batch correction using Harmony across each single-cell sample, and integrate the gene expression matrix of all samples into a whole matrix. After cell-level quality filtering, we used DubStepR packages for the selection of highly variable genes. The same methods were applied for normalization, scaling, and batch correction. To identify contaminated cell clusters, we performed dimensionality reduction with 50 harmony dimensions, n_neighbors set to 50, and min_dist set to 0.1, followed by clustering with a resolution of 0.2.

Representative gastric cancer tumor samples were collected for the spatial transcriptomic sequencing. Samples were cut into 5×5 mm² pieces, embedded in optimal cutting compound media and quickly frozen on dry ice. The frozen tissues were cryosectioned at 10 μm thickness by using the Thermo Scientific CryoStar NX50 cryostat and were placed in the capture area frames on the 10× Visium Spatial slides. Each sample slide was stained with H&E (hematoxylin Dako #S3309, eosin, Dako #CS701, bluing buffer #CS702, Dako, Glostrup, Denmark) and the brightfield images were captured via Leica whole-slide scanner at 10× resolution. Following tissue permeabilization, reverse transcription and cDNA amplification were processed by using Reagent Kit (#PN-1000184, PN-1000193, 10× Genomics). Visium spatial libraries were constructed using Visium spatial library construction kit (#PN-1000184, 10× Genomics) according to the manufacturer’s protocols. Finally, the libraries were sequenced using the Illumina NovaSeq6000 at least 100,000 reads per spot via pair-end 150 bp reading strategy.

Trajectory analysis was conducted using Monocle 3.0, accessible at http://cole-trapnell-lab.github.io/monocle3/. Utilizing the final annotated CD8+ T cell subtype Uniform Manifold Approximation and Projection (UMAP), we employed Monocle 3 to construct learned graphs and determine cell order, without setting the partition_qvalue. To identify genes that change in expression over pseudotime, the graph_test function within Monocle 3 was used, enabling us to pinpoint significant genes false discovery rate < 10^–50.

To identify malignant cell in the epithelial cell clusters, we utilized inferCNV (version 1.14.2), available at https://github.com/broadinstitute/inferCNV. We analyzed the count matrix designating normal sample pit cell and chief cell as the reference. The parameters set included a cutoff of 0.1, data scaling, denoising, a noise filter of 0.12, subcluster analysis mode, an hidden markov model type of i6, and the ward.D2 hierarchical clustering method.

Interaction between CD8+ T cell including effector memory T cell (Tem), exhausted T cell (Tex), proliferative exhausted T cell (pTex) and malignant cell was analyzed using NicheNet [15]. For NicheNet analysis, we set CD8+ T cell as a sender (as the regulator cell) and malignant cell as a receiver (as the target cell). We used the provided ligand-target pairs and ligand-receptor pairs as inputs and we performed differential NicheNet analysis between MSI-H and MSS gastric cancer. We used the gene set which are provided in NicheNet (https://github.com/saeyslab/nichnetr). To identify important interaction signal, we also utilized CellChat between CD8+ T cell and malignant cell [16]. We used the gene set provided by CellChat for reference (https://github.com/sqjin/cellchat).

Due to the low resolution of the Visium technology that couldn’t distinguish CD8+ T cell and malignant cell in the spot level, we calculated CD8+ T cell score using SEAGAL (ver. 2.4) python package. For the Tem we used GZMA, GZMB, IFNG, and PRF1 gene as signature genes, we used PDCD1, LAG3, TIGIT, and CTLA4 for Tex score. Next, we defined each Tem and Tex co-localized spot with the malignant cell when the score was over 0.2.

To quantify interaction between malignant cell and adjacent CD8+ T cell co-localized spot, we calculated the number of adjacent spots identified as CD8+ T cell co-localized spot. This count was the divided by seven which was the number of adjacent spot and core spot of Visium to derive a malignant-CD8+ T cell interaction score.

We utilized CIBERSORTx to quantify the proportions of immune cell types within tumor microenvironment [17]. For signature matrix of 22 immune cell types, we used CIBERSORTx provided signature matrix. Then we compute the fractions of immune cell types in GSE66229 microarray data.

In GSE66229 dataset, we defined MSI status according to MLH1 imunohistochemical result that was provided by original article [18]. Then, we excluded epithelial-mesenchymal transition (EMT) subtype of gastric cancer from the original dataset due to EMT subtype showed distinct characteristics of MSI or MSS samples. Cox proportional hazard analysis in MSI-H and MSS samples was performed using lifeline library of python. For maximal survival differences according to IFNG, GZMB, LAG3, TIGIT, and CTLA4 expression in MSI-H samples, we used MaxStat R package. Overall survival of IFNG, GZMB, LAG3, TIGIT, and CTLA4 expression high and low patients was analyzed by log-rank tests with Kaplan-Meier curves. A p < 0.05 indicated statistical significance.

Droplet-based scRNA-seq was performed on tissue sample obtained from six patients. Clinicopathological features of patents were described in Table 1. After quality control and doublet removal, 45,087 cells from six tumor samples and two normal samples were included in the dataset. After multi-sample batch-effect correction and data normalization, unsupervised cell clustering was performed using the UMAP embedding 45,087 cells. The resulting UMAP plot highlights distinct clusters through color differentiation (Fig. 1A). To gain insights into the identity of each cell type, we superimposed gene expression levels of established markers onto the UMAP plot, providing a view of clusters (Fig. 1B). Further elucidating the composition of the identified clusters, a dot plot demonstrates proportions and average expression levels of marker genes across 10 cell types (Fig. 1C). The single cell atlas comprised of epithelial cell (n=3,188; EPCAM, CDH1), endothelial cell (n=1,270; PLVAP), stromal cell (n=1,326; THY1), CD8+ T cell (n=19,287; CD8A), CD4+ T cell (n=6,764; CD4), B cell (n=5,441; MS4A1), plasma cell (n=2,961; TNFRSF17), mast cell (n=1,354; KIT), neutrophil (n=1,485; CXCL8), and myeloid cell (n=2,011; CD163) (Fig. 1A and C). The UMAP plots demonstrate the differences in cell count across clusters within MSI-H, MSS, and normal samples (Fig. 1D). Examining the broader context, stacked bar plots illustrate the distribution of clusters among samples and demonstrate the fraction of cells within MSI-H, MSS, and normal samples (Fig. 1E and F). In MSI-H samples, a significantly higher proportion of tumor infiltrating CD8+ T cell was observed.

After multi-sample batch-effect correction and data normalization, unsupervised cell clustering was performed, utilizing UMAP to embed 18,150 CD8+ T cell. This analysis identified nine distinct CD8+ T cell clusters, each characterized by unique gene signatures (Fig. 2A). A dot plot elucidates the proportions and average expression levels of marker genes, contributing to an understanding of CD8+ T cell subtypes (Fig. 2B). The single cell atlas comprised of naïve CD8+ T cell, activated naïve CD8+ T cell, effector memory CD8+ T cell (Tem), exhausted CD8+ T cell (Tex), proliferative exhausted CD8+ T cell (pTex), proliferative CD8+ T cell, mucosal-associated invariant T (MAIT) cell, FCER1G+ natural killer T (NKT) cell, and FCGR3A+ NKT cell (Fig. 2A). Violin plots show the expression of TOX, TOX2, CTLA4, and CXCL13 in Tex and pTex (Fig. 2C). The stacked bar plot reveals the distribution of CD8+ T cell subtypes among MSI-H, MSS, and normal samples (Fig. 2D). In MSI-H gastric cancer, Tem, Tex, pTex, and proliferative CD8+ T cell exhibited significantly higher proportions, whereas MAIT cell and NKT cell showed significantly higher proportions in MSS gastric cancer. In Tex and pTex, the exhaustion marker LAG3 was significantly upregulated in MSI-H gastric cancer; however, CXCL13 and TOX were significantly upregulated in MSS gastric cancer (Fig. 2E). Also, the effector function marker IFNG, GZMB, GZMH, and GZMK were significantly upregulated in Tex and pTex of MSI-H gastric cancer (Fig. 2F). The trajectory analysis, originating from naïve CD8+ T cell, demonstrates the differentiation of CD8+ T cell into Tex and pTex as pseudotime progressed (Fig. 2G). Furthermore, it was observed that despite the progression of exhaustion, the expression of IFNG remained elevated in MSI-H gastric cancer. Although the expression of GZMK and GZMH decreased as pseudotime progressed, the expression of these genes was higher in MSI-H gastric cancer (Fig. 2H).

A comprehensive analysis of epithelial cell was conducted to unravel the cellular subtype according to the molecular characteristics (Fig. 3A). We can identify the expression levels of epithelial marker genes across clusters, including TFF1 and MUC5AC for pit cell, LIPF for chief cell, APOB and ALDOB for enterocyte, and CEACAM5 and CCND2 for malignant cell (Fig. 3B). The chromosomal landscape, revealing inferred large-scale copy number variations (CNV) within the epithelial sub-clusters, identified clusters 0, 1, 3, 4, 6, and 8 as malignant cell clusters (Fig. 3C). Based on the gene expression patterns observed in the dot plot and the inferred CNV within the sub-clusters of epithelial cells, we have annotated sub-clusters on the UMAP (Fig. 3D).

To examine the intercellular communication between malignant cell and CD8+ T cell, receptor-ligand interactions were identified. The scaled expression of IFNG and GZMB ligands was observed to be retained in Tex and pTex compared to Tem in MSI-H gastric cancer. However, this distinctive pattern was not observed in MSS gastric cancer (Fig. 4A). Subsequently, the inferred interferon γ (IFN-γ) signaling pathway networks calculated by CellChat demonstrate the enriched network between Tex, pTex and malignant cell in MSI-H gastric cancer and decreased network between Tex, pTex and malignant cell in MSS gastric cancer (Fig. 4B and C).

We further conducted spatial transcriptomic analysis to verify molecular characteristics with spatial distribution. The spatial transcriptomic data shows the expression of EPCAM and CD8A at each spot. Overall, in MSI-H gastric cancer, it was observed that the distribution of CD8A expression spots is diffuse and prominently elevated compared to MSS gastric cancer (Fig. 5A and B). This suggests a significant presence of tumor infiltrating CD8+ T cell in MSI-H gastric cancer.

To assess the interaction between Tem and Tex with malignant cell, co-localization scores for Tem and Tex were evaluated using cor.test and visualized in the spatial transcriptome (Fig. 6A and B). The fraction of co-localized spots for Tem and Tex was found to be significantly higher in MSI-H gastric cancer (Fig. 6C). Based on the spatial distribution of malignant cell and Tem, Tex spots, we calculated interaction scores (Fig. 6D). Significantly higher interaction scores were observed in MSI-H which indicating the abundant presence of Tem and Tex cells in close proximity to malignant cell in MSI-H gastric cancer (Fig. 6E). Finally, we evaluated the IFN-γ pathway between the adjacent Tex and malignant cell, which demonstrates significantly higher interaction in MSI-H gastric cancer. (Fig. 6F).

Further examination was performed by reanalyzing bulk transcriptomic ACRG data (GSE66229) using CIBERSORTx deconvolution, revealing a significantly higher presence of CD8+ T cells in MSI-H compared to MSS gastric cancer (Fig. 7A). In the heat map, it was observed that cytotoxic and exhaustion markers are highly expressed in MSI-H gastric cancer (Fig. 7B). The scatter plot between IFNG and marker genes of exhaustion illustrates a higher coefficient value in MSI-H gastric cancer, indicating elevated expressions of IFNG at the same level of CXCL13 and LAG3 expression (Fig. 7C). Survival analysis demonstrates the significant better survival in MSI-H patients according to the expression of effector function and exhaustion marker genes (Fig. 7D and E).