Mutual exclusivity and co-occurrence patterns of immune checkpoints indicate NKG2A relates to anti-PD-1 resistance in gastric cancer

Abstract

Background

An increasing number of clinical studies have begun to explore combination strategies with immune checkpoint inhibitors, aiming to present new opportunities for overcoming anti-PD-1 treatment resistance in gastric cancer. Unfortunately, the exploration of certain immune checkpoint inhibitor combination strategies has yielded suboptimal results. Therefore, it is necessary to comprehensively analyze the expression patterns of immune checkpoints and identify optimal combination regimens of anti-PD-1 inhibitors with other immune checkpoint inhibitors.

Methods

Leveraging single-cell RNA sequencing (scRNA-seq) and multivariate linear regression interaction models, we dissected the immune checkpoint expression characteristics of CD8⁺ T cells in gastric cancer and the immune checkpoint expression pattern (ICEP) mediating anti-PD-1 treatment resistance. Furthermore, we employed transcription factor analysis and CellOracle to explore the transcriptional regulatory mechanisms governing CD8⁺ T cell differentiation fates. Finally, we utilized Nichenet and spatial transcriptomic analysis to investigate the spatial expression patterns of immune checkpoints.

Results

Interaction analysis indicated that, among the known immune checkpoints, co-expression of NKG2A and PD-1 might exert a more profound inhibitory effect on the proliferative capacity of CD8⁺ T cells. The co-expression analysis revealed differential co-expression pattern of PD-1 and NKG2A, defined as ICEP1 (CD8⁺ T cells co-expressing PD-1, CTLA-4, TIGIT, LAG-3 or CD38) and ICEP2 (CD8⁺ T cells solely expressing NKG2A or co-expressing with other immune checkpoints), reflecting the co-occurrence pattern of PD-1 and the mutual exclusivity of NKG2A. Further, these two ICEP CD8⁺ T cell subsets represented distinct CD8⁺ T cell differentiation fates governed by MSC and RUNX3. Notably, ICEP2 CD8⁺ T cells were associated with anti-PD-1 therapy resistance in gastric cancer. This phenomenon may be attributed to the recruitment of LGMN⁺ macrophages mediated by the CXCL16-CXCR6 signaling pathway.

Conclusion

This study unveiled two distinct ICEPs and the mutually exclusivity and co-occurrence characteristics of CD8⁺ T cells in gastric cancer. The ICEP2 CD8⁺ T cell subset, highly expressed in gastric cancer patients resistant to anti-PD-1 therapy, may be recruited by LGMN⁺ macrophages through CXCL16-CXCR6 axis. These findings provide evidence for NKG2A as a novel immunotherapeutic target in gastric cancer and offer new insights into combination strategies for immune checkpoint inhibitors in gastric cancer.

Graphical Abstract

Introduction

The CheckMate-649 clinical trial solidified the pivotal role of anti-PD-1 therapy in the first-line treatment of advanced gastric cancer and provided compelling evidence for the application of chemotherapy combined with immunotherapy in clinical practice for gastric cancer. However, a subset of gastric cancer patients failed to achieve optimal responses to anti-PD-1 therapy, potentially due to the involvement of immune checkpoints other than PD-1 on T cells, leading to target evasion and resistance to single-agent anti-PD-1/PD-L1 therapy [1]. Consequently, an increasing number of clinical studies have begun exploring “dual immunotherapy” in gastric cancer, which involves co-inhibition of multiple immune checkpoint targets to address target evasion and resistance. Thus, exploring combination regimens of immune checkpoint inhibitors presents new opportunities for overcoming anti-PD-1 treatment resistance in gastric cancer.

Unfortunately, the exploration of certain immune checkpoint inhibitor combination strategies has yielded suboptimal results. The NCT03662659 clinical trial demonstrated that, compared to the combination of anti-PD-1 inhibitor and chemotherapy group, the anti-LAG-3 inhibitor combined with anti-PD-1 inhibitor and chemotherapy showed no significant efficacy advantage in patients with metastatic gastric or gastroesophageal junction adenocarcinoma and had a higher incidence of severe adverse events [2]. This highlights the need for a deeper understanding of the underlying mechanisms linking PD-1 and LAG-3 immune checkpoints in gastric cancer. Another clinical trial, NCT03033576, revealed that in advanced melanoma, the combination of anti-PD-1 and anti-CTLA-4 inhibitors yielded superior efficacy compared to anti-CTLA-4 inhibitor monotherapy. However, there was no statistically significant difference in the levels of tumor-infiltrating CD8⁺ T cells between the two treatment groups, which contradicted findings from basic research studies [3]. This discrepancy may be attributed to the lack of a comprehensive understanding of the intrinsic associations and overall patterns of immune checkpoint expression. Concurrently, current basic research on immune checkpoint inhibitor combinations has primarily focused on immune checkpoints associated with CD8⁺ T cell exhaustion, such as PD-1, LAG-3, CTLA-4, and TIGIT, while the impact of other immune checkpoint expressions on CD8⁺ T cell function has received limited attention [4,5,6]. Therefore, there is an urgent need for a more comprehensive understanding of immune checkpoint expression patterns (ICEPs) in gastric cancer and in-depth exploration of the intrinsic relationships among immune checkpoint expressions.

Single-cell RNA sequencing (scRNA-seq) is a powerful tool for analyzing gene expression at the single-cell level, enabling the study of the functional impact on CD8⁺ T cells when immune checkpoints are expressed individually or simultaneously. In this study, we employed scRNA-seq, multivariate linear regression interaction models, and spatial transcriptomic sequencing to unveil two distinct ICEPs and mutually exclusive and co-occurrence characteristics in gastric cancer CD8⁺ T cells. Notably, the ICEP2 CD8⁺ T cell subset, defined by the core immune checkpoint NKG2A, mediated resistance to anti-PD-1 inhibitor therapy in gastric cancer. Furthermore, through scRNA-seq analysis of intercellular communication, the Cancer Genome Atlas-Stomach Adenocarcinoma (TCGA-STAD) cohort, and spatial transcriptomic sequencing, we identified an immunosuppressive LGMN⁺ macrophage population in the gastric cancer tumor core that recruited ICEP2 CD8⁺ T cells into the tumor-infiltrating region via the CXCL16-CXCR6 chemokine axis. Our study provides new evidence for NKG2A as an immunotherapeutic target in gastric cancer and offers novel insights into combination strategies for immune checkpoint inhibitors in gastric cancer.

Material and methods

Acquisition of single-cell transcriptomic data

The raw sequencing data for the gastric cancer single-cell cohorts GSE183904 and OMIX001073 were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) supported by the National Center for Biotechnology Information (NCBI) and the Open Archive for Miscellaneous Data (OMIX) database (https://ngdc.cncb.ac.cn/omix/releaseList), respectively.

Additionally, to explore the ICEPs mediating anti-PD-1 treatment resistance in gastric cancer patients, we obtained scRNA-seq data from fresh tumor tissues of 17 gastric cancer patients at Nanfang Hospital (Guangzhou, China). All gastric cancer patients signed informed consent forms regarding the use of their clinical information and tissue samples. The detailed clinical information of the single-cell cohorts is provided in Additional file 4: Table S3.

Fresh human tumor tissues were collected and minced into small pieces on ice. Subsequently, the tissues were digested using the Tissue Dissociation Kits (Miltenyi, 130,110,201). The procedure involved mixing the tissues with 5 ml enzyme mix (4.7 mL RPMI1640 + 200μL enzyme H + 100 μl enzyme R + 25 μl enzyme A) and incubating with agitation at 37 ℃ for up to 40 min to obtain a single-cell suspension. The collected cells were then pelleted, resuspended, and filtered through a 30 μm cell strainer to obtain a uniform cell suspension. Samples meeting the following criteria were selected for library construction: viability > 80%, nucleus cell rate > 70%, and clump rate < 20%.

The ChromiumTM Controller and ChromiumTM Single Cell 5′ Reagent Version 2 Kit from 10 × Genomics were employed to conduct library construction experiments. In brief, GemCode Technology was leveraged to encapsulate sorted cells, reagents, and gel beads containing barcoded oligonucleotides into nanoliter-sized GEMs. Polyadenylated mRNA from individual cells was lysed and subjected to barcoded reverse transcription within each GEM. Post-RT GEMs were purified and the cDNA amplified. The cDNA was then fragmented and repaired on the ends, and an A-tail was added to the 5′ end. Adaptors were ligated to the fragments which underwent double-sided size selection with SPRI. Following sample index PCR, another round of double-sided SPRI selection was performed. Real-time quantitative PCR with TaqMan probes was used to assess quality and quantity of the final library. The ultimate products were sequenced on the Xten-PE151 platform from BGIShenzhen, China.

Acquisition of bulk transcriptomic data

The transcriptomic expression values, patient clinical information, and survival data for Stomach Adenocarcinoma (STAD) from TCGA were downloaded from the UCSC Xena database (https://xenabrowser.net/datapages/). The downloaded mRNA expression values were in the format of Transcripts Per Kilobase of exon model per Million mapped reads (TPM). All analyses were performed in R (version 4.2.0).

Acquisition of spatial transcriptomic data

Fresh tumor tissue samples were collected from a gastric cancer patient at Nanfang Hospital, embedded, and cryosectioned to obtain ten frozen sections. These sections were placed in enzyme-free tubes for RNA extraction, with the requirement that the extracted RNA must be of high quality, having an RNA Integrity Number (RIN) greater than 7. To optimize tissue permeabilization conditions, a tissue optimization chip with eight capture areas was utilized. The frozen sections were mounted onto the corresponding areas of the chip, followed by H&E staining and imaging for tissue permeabilization. Under the action of tissue permeabilization enzymes, the released mRNA from cells was captured by probes, forming fluorescently labeled cDNA. The optimal permeabilization time was determined based on the fluorescence imaging results. Once the optimal permeabilization conditions were established, spatial transcriptomic library construction and sequencing could proceed. The required spatial gene expression chips contained either two or four capture areas, each comprising 4992 spots. Each spot contained a unique spatial barcode and millions of nucleotide primers. After permeabilization, the released mRNA from cells hybridized with the barcoded primers, followed by cDNA synthesis and library construction. Depending on the tissue type, each barcode could correspond to 1–10 cells. Prior to spatial transcriptomic sequencing, the frozen sections were mounted onto the capture areas of the gene expression chips and subjected to H&E staining and imaging. Tissue permeabilization was performed according to the previously determined conditions, allowing the released mRNA from cells to be captured by the primers on the chip. By adding reverse transcription reagents (RT Master Mix) and incubating with the sample, full-length barcoded cDNA was generated on the slides. The cDNA was used as a template for PCR amplification. After amplification, the amplified products underwent quality control, including checking the size of the amplified fragments and the yield of the amplified products. Qualified amplified products were further used for sequencing library construction. First, the cDNA was chemically fragmented into approximately 200–300 bp fragments. Subsequently, fragmentation, end-repair, and A-tailing were performed, followed by cDNA fragment selection. Next, the P7 Adapter was ligated, and sample dual-end indexing was introduced through PCR amplification. Finally, fragment selection was carried out to obtain the cDNA library.

Initial processing and quality control of single-cell RNA-seq datasets

CellRanger (version 4.0.0), developed by 10 × Genomics, was used to process raw sequencing data, including alignment, quantification, simple filtering and quality control to generate initial gene expression matrices based on the human reference genome GRCh38. The R package Seurat (version 4.4.0) was then applied for downstream quality control and analysis [7]. High-quality cells were retained by filtering for: (i) cells expressing 500–6000 genes, and (ii) cells with mitochondrial gene UMI counts accounting for less than 20% of total UMI counts. (iii) Genes expressed in at least 3 cells were also kept. After quality control, 65,976 high-quality cells were retained for subsequent analysis.

Downstream analysis of single-cell RNA-seq datasets

Normalization of counts was performed using the NormalizeData function. The FindVariableFeatures function was applied to select the top 2000 most variable genes across cells. Dimensionality reduction was achieved through PCA after gene normalization via ScaleData. Batch effects were corrected across samples using scVI [8]. Clustering was performed at a resolution of 0.6 and visualized through UMAP embedding. Finally, all cells were annotated into 8 cell types (epithelial, B cells, plasma cells, macrophages, fibroblasts, endothelial cells, mast cells, T cells and NK cells) according to cell-specific marker genes.

Multivariate linear regression interaction model

To assess the relationship between two immune checkpoints and proliferation of CD8⁺ T cells, we employed a multivariate linear regression model [9]. First, a proliferation score was computed for each cell in the scRNA-seq data based on a gene signature of T cell proliferation using Ucell. P and O represent expression levels of the two immune checkpoints. The following regression equation was fitted:

$$\text{d }+ {\text{a}}^{*}\text{ P }+ {\text{b}}^{*}\text{ O }+ {\text{c}}^{*} {\text{O}}^{*}\text{ P }=\text{ proliferation}$$

where proliferation is the outcome variable, and d, a, b and c are regression coefficients. P and O are the independent variables representing expression levels of the two checkpoints.

d is the intercept, representing the expected mean proliferation score when both independent variables P and O are zero.

a is the regression coefficient for P, indicating the effect of a single unit change in P on proliferation when O is held constant.

b is the regression coefficient for O, indicating the effect of a single unit change in O on proliferation when P is held constant.

c is the regression coefficient for the interaction term P*O, representing the joint effect of changes in both P and O on proliferation.

This model considers the interaction between the independent variables via the P*O term. The interaction term allows the model to capture non-linear relationships where the joint effect of the variables exceeds their individual effects.

The interaction test t-value was computed as c/StadErr using ordinary least squares estimation. A negative t-value would indicate a synergistic inhibitory effect of co-expression levels of the two checkpoints on CD8⁺ T cell proliferation, whereas a non-significant value suggests no interaction effect.

Cell trajectory analysis

Pseudotime analysis was performed using the R package Monocle3 (version 0.2.0) to infer developmental trajectories and temporal information of macrophages [10]. We use preprocess_cds for data preprocessing, followed by reduce_dimension function to perform UMAP dimensionality reduction. Then we use cluster_cells and plot_cells functions for unsupervised clustering and visualization. Finally, we use the default parameters of the learn_graph function to infer cell differentiation trajectories.

Single-cell transcriptomics NicheNet analysis of cell–cell communication

NicheNet is a computational approach that leverages prior knowledge of ligand-receptor pairs, signaling pathways and gene regulatory networks to prioritize ligands that may regulate the expression of gene sets of interest [11]. Gene sets were defined as upregulated genes in each cell type of interest (non-responders vs responders). Thresholds of P-value ≤ 0.05 and LogFC ≥ 0.50 were applied for each cell type of interest. The analysis included only ligands expressed in ≥ 5% of cells in at least one cell type of interest. This was followed by ligand-receptor activity analysis between CD8⁺ T cell subclusters and macrophages.

Transcription factor analysis

SCENIC, an approach based on co-expression and motifs, was employed to reconstruct gene regulatory networks (GRNs) from single-cell transcriptomics data to predict transcription factors and regulatory elements activated in different cell types [12]. The “hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather” dataset was used for transcription factor annotation. Tools like GRNBoost2 (Gradient Boosting) were applied to predict co-expression modules between transcription factors and target genes, inferring potential transcription factors in GRN. RcisTarget analyzed genes in all co-expression modules to identify enriched transcription factor binding sites, constructing transcription factor-target networks (TF-targets networks) by filtering modules with transcription factor set site enrichment and their targets. Each transcription factor and its potential direct targets were termed a regulon. AUCell scores were calculated for each regulon to determine its activity status in each cell. SCENIC finally identified differences in transcription factor and regulon activities across cell subclusters.

CellOracle gene perturbation modeling

We employed CellOracle[13] to simulate cell identity shifts following transcription factor knockout, computing pseudotime and perturbation scores based on the contextual gene regulatory network constructed from our gastric cancer scRNA-seq data. CellOracle analysis was performed following the official tutorial with default parameters (https://morris-lab.github.io/CellOracle.documentation/).

Prediction of immunotherapy response

The Tumor Immune Dysfunction and Exclusion (TIDE) score is a web-based tool (http://tide.dfci.harvard.edu/login/) [14] used to predict the efficacy of immune checkpoint inhibitors. It quantifies the TIDE score for each patient by calculating the activity of two mechanisms: immune inhibitory factors and T cell functional inactivation.

CIBERSORTx deconvolution of immune cell infiltration from transcriptomic data

To computationally analyze the immune infiltration of TME cell types in the TCGA-STAD cohort corresponding to our scRNA-seq data of gastric cancer immunotherapy response, the online tool CIBERSORTx [15] was used. A reference signature matrix was built from the single-cell RNA-seq dataset and utilized to infer relative cell type proportions in the TCGA-STAD cohort. Pearson correlation analysis was performed to evaluate relationships between cell type infiltrations, considering correlations with P-values < 0.05 to be significant. Through these steps, the immune infiltration of different cell types in the TCGA-STAD cohort was analyzed and described accurately, providing insights into interactions between cell populations in the gastric cancer microenvironment.

Initial processing and quality control of spatial transcriptomics datasets

The raw data was processed using the official 10 × Genomics analysis software Space Ranger (version 2.0.0) (https://support.10xgenomics.com/spatial-gene-expression/software/overview/welcome) for filtering, alignment, quantification and generation of spot gene expression matrices. Seurat R package was used to compute the percentage of mitochondrial genes in each spot and spots with percentages exceeding a threshold of 20% were filtered.

Downstream analysis of spatial transcriptomics datasets

A Seurat spatial object was created from the data using CreateSeuratObject and Read10X_Image functions. Dimensionality reduction was performed using SCTransform to normalize, detect highly variable features and embed the data in SCT analytics. Data was further scaled using ScaleData, followed by PCA dimensionality reduction using RunPCA. Clustering was done at a resolution of 0.6 using FindNeighbors and FindClusters.

Cell type annotation of spatial transcriptomics datasets

The Bayesian method Cell2location [16] leverages single-cell transcriptomics as a reference to perform spatial cell type mapping. It involves two key steps. Firstly, characteristic signatures for reference cell types are computed using negative binomial regression, allowing integration across technologies and batches for improved flexibility and accuracy. Secondly, one or more spatial transcriptomics datasets are parsed to determine abundances for each location and cell type. This step provides detailed spatial distribution insights revealing cell positioning and interactions in tissues. Cell2location first uses negative binomial regression on scRNA-seq annotations to infer characteristic cell type signals. Taking the reference signals and the spatial transcriptomic dataset as input, Cell2location then computes abundances of different cell types across spots in the spatial data.

Spatial transcriptomics Cellchat analysis of cell–cell communication

The R package CellChat (version 1.6.1) was implemented to infer communication between TME subclusters at the spatial transcriptomics level [17]. The gene expression matrix and cell annotations were extracted from the Seurat object using GetAssayData. Subsequent ligand-receptor inference was based on the CellChatDB. Human dataset from the CellChatDB database. Overexpressed genes for each cell type were identified using identifyOverExpressedGenes, and overexpressed interactions using identifyOverExpressedInteractions. The ligand-receptor interactions were projected onto the protein–protein interaction network using projectData. Communication probabilities between cells were computed using computeCommunProb to infer a cell communication interaction network. This revealed the quantity and intensity of interactions between different cell types.

Correlation analysis and survival analysis of bulk transcriptomic data

RNA-seq data and clinical survival information of TCGA-STAD were downloaded from the UCSC Xena database. Log2 transformation was applied for subsequent analysis. Expression values of gene sets of interest were extracted from the dataset. Pearson correlation analysis was performed to determine correlations between gene expressions and a significance of P-value < 0.05 was considered. To calculate the correlation between immune infiltration and patient survival time, the surv_cutpoint function in the R survminer package was used to compute the optimal cut-off of immune infiltration levels in TCGA-STAD patients based on CIBERSORTx estimates. Patients were divided into high and low infiltration groups. Kaplan–Meier survival curves were fitted using survfit and visualized using ggsurvplot to assess the impact of infiltration levels on survival outcomes.

Results

Interaction analysis reveals the differential co-expression effect with PD-1

CD8⁺ T cells play a crucial role in anti-tumor immune responses. As potent effector cells, CD8⁺ T cells possess the ability to directly kill tumor cells as well as activate and regulate other the functions of other immune cells. This makes CD8⁺ T cells an important research subject in tumor immunotherapy. To comprehensively delineate the intrinsic associations between immune checkpoint expressions and the impact of different expression patterns on CD8⁺ T cell function, we obtained gastric cancer single-cell cohorts GSE183904 and OMIX001073 from the GEO and OMIX databases, respectively. After the process described in the Materials and Methods section, including quality control for all samples to exclude low-quality cells, a total of 207,874 cells were retained. Subsequently, we performed batch effect removal and dimensionality reduction clustering analysis, and based on typical marker genes, we classified them into eight major cell types, including epithelial cells (EPITHELIAL, highly expressing EPCAM, KRT19, PROM1, ALDH1A1, CD24), B lymphocytes (B-CELL, highly expressing CD19, MS4A1, CD79A), plasma cells (PLASMA, highly expressing IGHG1, MZB1, CD38, TNFRSF17, SDC1), myeloid cells (MYELOID, highly expressing CD68, CD163, CD14), fibroblasts (FIBROBLAST, highly expressing FGF7, MME), endothelial cells (ENDOTHELIAL, highly expressing PECAM1, VWF), mast cells (MAST, highly expressing KIT), and T lymphocytes and NK cells (LYMPHOID & NK, highly expressing CD3D, CD3E, CD8A, CD8B, CD4 or highly expressing FGFBP2, FCG3RA, CX3CR1) (Additional file 1: Figs. S1A, B).

Next, we extracted a total of 40,381 CD8⁺ T cell subpopulations and further divided them into seven CD8⁺ T cell subgroups, including CD8⁺ naive T cells (CD8⁺ Tn, highly expressing TCF7, LEF1, MAL), CD8⁺ IL7R⁺ memory T cells (CD8⁺ IL7R⁺ Tm, highly expressing IL7R, GPR183, ZFP36L2, CXCR4), CD8⁺ KLRC1⁺ tissue-resident memory T cells (CD8⁺ KLRC1⁺ Trm, highly expressing KLRC1, KLRD1, ZNF683, ITGAE, CXCR6), CD8⁺ GMZK⁺ effector T cells (CD8⁺ GMZK⁺ Teff, highly expressing GZMK, NKG7, GZMA, PDCD1, TIGIT), CD8⁺ GMZH⁺ effector T cells (CD8⁺ GMZH⁺ Teff, highly expressing GZMH, GZMA, IFNG, NKG7), CD8⁺ exhausted T cells (CD8⁺ Tex, highly expressing PDCD1, TIGIT, LAG3, CTLA4, TOX, CD38), and Proliferating CD8⁺ T cells (highly expressing MKI67) (Fig. 1A, B).

Fig. 1

Interaction analysis between immune checkpoints on gastric cancer CD8⁺ T cells. A UMAP plot of CD8⁺ T cell subsets in gastric cancer patients, color-coded by cell type. B Dot plot illustrating the average expression values of known marker genes across different cell types. The size and color of the dots represent the average expression levels of the genes. C Proportion of immune checkpoint-positive cells in CD8⁺ cells gastric cancer patients. D Heatmap of k-means clustering for seven immune checkpoints. E Schematic representation of the multivariate linear regression interaction equation. F Volcano plot displaying the interaction test t-values for the multivariate linear regression interaction equations involving the seven immune checkpoints. G Scatter plot illustrating the relationship between KLRC1, PDCD1, and T cell proliferation

Furthermore, we collected 16 immune checkpoints reported in previous research, including immune checkpoint co-inhibitory molecules PD-1 (PDCD1), CTLA-4 (CTLA4), LAG-3 (LAG3), TIGIT, HAVCR2 (TIM3), NKG2A (KLRC1), and VISTA (VSIR); immune checkpoint co-stimulatory molecules TNFRSF18 (GITR), TNFRSF4 (OX40), TNFRSF9 (4-1BB), TNFRSF14 (HVEM), ICOS (CD278), CD27, CD28; and other immune-related molecules ENTPD-1 (CD39), CD38. By calculating the expression proportions of immune checkpoint molecules in T cells, the results showed that PDCD1, CTLA4, LAG3, TIGIT, KLRC1, VSIR, and CD38 had expression proportions exceeding 10% in T cells. Notably, KLRC1 was primarily expressed on CD8⁺ T cells, with a lower expression proportion on CD4⁺ T cells. In contrast, CTLA4 had a higher expression proportion on CD4⁺ T cells (Fig. 1C, Additional file 1: Fig. S1C). To explore the expression of immune checkpoints in CD8⁺ T cells, we performed k-means clustering based on the expression values of the seven immune checkpoints in CD8⁺ T cells, resulting in 50 subgroups with co-expression of immune checkpoints, suggesting the existence of co-expression phenomena in CD8⁺ T cells (Fig. 1D).

Immune checkpoints exhibit a widespread phenomenon of co-expression, while some checkpoints tend to be mutually exclusive, inducing T cell dysfunction through shared or distinct mechanisms. Whether the impact on T cell function is consistent when different immune checkpoints co-express remains to be investigated. Therefore, we utilized a multivariate linear regression interaction equation to analyze the relationship between the co-expression of two immune checkpoints and CD8⁺ T cell function. First, based on the cell cycle and DNA replication gene set (Additional file 2: Table S1), we quantified the proliferation function of each CD8⁺ T cell using the R package Ucell. Subsequently, we quantified the relationship between CD8⁺ T cell proliferative function and the expression of two immune checkpoints by calculating the interaction test t-value (Coef/StadErr). If the t-value is negative, it indicates that the co-expression levels of the two immune checkpoints have a synergistic inhibitory effect on CD8⁺ T cell proliferative function; otherwise, there is no synergistic inhibitory effect (the interaction equation P-value indicates the statistical significance of the t-value) (Fig. 1E). The results showed that only KLRC1 with PDCD1, and VSIR with KLRC1 had negative t-values, while the t-values between PDCD1, CTLA4, TIGIT, LAG3, and CD38 were all positive (Fig. 1F, G, Additional file 3: Table S2).

Therefore, these results suggest that the co-expression of the immune checkpoints NKG2A and PD-1 on CD8⁺ T cells may maximally inhibit the proliferative function of CD8⁺ T cells. This provides new insights for combination therapy regimens with anti-PD-1 inhibitors in gastric cancer.

Co-expression analysis demonstrates the mutual exclusivity and co-occurrence pattern of immune checkpoints

PD-1 and NKG2A exert reciprocal inhibitory effects on CD8⁺ T cells, suggesting potential inhibition of T cell function through distinct mechanisms. Next, we aim to investigate whether PD-1 and NKG2A co-express or exhibit different expression patterns at the cellular subset level. We first obtained fresh tumor tissue scRNA-seq from 17 gastric cancer patients, including samples from different stages of anti-PD-1 treatment and different best overall responses (BOR). Among these, samples from 8 patients were collected after anti-PD-1 treatment, divided into two efficacy groups: post-treatment non-responders (Post-NR, including post-treatment stable disease (Post-SD)/post-treatment progressive disease (Post-PD), n = 5) and post-treatment responders (Post-R, including post-treatment partial response (Post-PR), n = 3). For the remaining 9 patients, samples were collected through biopsy before anti-PD-1 treatment, and efficacy information was collected after treatment, similarly divided into two efficacy groups: pre-treatment non-responders (Pre-NR, including pre-treatment stable disease (Pre-SD)/ pre-treatment progressive disease (Pre-PD), n = 6) and pre-treatment responders (Pre-R, including pre-treatment partial response (Pre-PR), n = 3) (Fig. 2A, Additional file 4: Table S3). After the process described in the Materials and Methods section, a total of 65,979 cells were retained. Subsequently, batch effect removal and dimensionality reduction clustering analysis were performed, and the cells were classified into 8 major cell types, including 24,273 epithelial cells (EPITHELIAL), 1822 B lymphocytes (B-CELL), 5581 plasma cells (PLASMA), 9894 myeloid cells (MYELOID), 4,816 fibroblasts (FIBROBLAST), 2923 endothelial cells (ENDOTHELIAL), 1340 mast cells (MAST), and 15,327 T lymphocytes and NK cells (LYMPHOID & NK) (Additional file 5: Fig. S2A, B, C, Additional file 6: Table S4).

Fig. 2

Mutual exclusivity and co-occurrence pattern of immune checkpoints in gastric cancer. A Schematic of patient grouping information. B Bar chart illustrating the percentages of gastric cancer patients with individual immune checkpoint expression, co-expression of two immune checkpoints (ICs), co-expression of three immune checkpoints, and co-expression of four immune checkpoints. C Upset plot illustrating patterns of co-occurrence of immune checkpoints. D UMAP plot of CD8⁺ T cell subpopulations in gastric cancer patients, color-coded by cell type. E Heatmap displaying the distribution of ICEPs in cell subpopulations, estimated based on the ratio of observed cell counts to expected cell counts (Ro/e) calculated using a chi-square test

Similarly, after extracting T cells and CD8⁺ T cells and calculating the expression proportions of immune checkpoints, we found that PDCD1, CTLA4, LAG3, TIGIT, KLRC1, VSIR, and CD38 had expression proportions exceeding 10% in T cells. Consistent with the external cohort results, KLRC1 was mainly expressed on CD8⁺ T cells, and CTLA4 had a higher expression proportion on CD4⁺ T cells (Additional file 5: Fig. S2D). Subsequently, when comparing the proportions of each immune checkpoint expressed alone or co-expressed with other immune checkpoints, we found that PDCD1, CTLA4, TIGIT, LAG3, and CD38 exhibited a certain proportion of co-expression with three immune checkpoints, and PDCD1, CTLA4, TIGIT, and LAG3 even showed co-expression with four immune checkpoints. In contrast, compared to other immune checkpoints, KLRC1 had the highest proportion of single expression, reaching 40%. Furthermore, the proportion of KLRC1 co-expressed with three immune checkpoints was the lowest, and there was no co-expression with four immune checkpoints (Fig. 2B-C).

To delineate the co-expression characteristics of immune checkpoints in CD8⁺ T cells, we performed k-means clustering based on the expression values of the seven immune checkpoints in CD8⁺ T cells, resulting in 50 subgroups. Using a threshold of 1, we named the immune checkpoint expression features of the 50 subgroups, ultimately obtaining 34 ICEPs (Additional file 3: Fig. S2E). To further understand the CD8⁺ T cell state represented by the ICEPs, we divided CD8⁺ T cells into 6 subgroups, including CD8⁺ naive T cells (CD8⁺ Tn), CD8⁺ IL7R⁺ memory T cells (CD8⁺ IL7R⁺ Tm), CD8⁺ KLRC1⁺ tissue-resident memory T cells (CD8⁺ KLRC1⁺ Trm), CD8⁺ GMZK⁺ effector T cells (CD8⁺ GMZK⁺ Teff), CD8⁺ GMZH⁺ effector T cells (CD8⁺ GMZH⁺ Teff), and CD8⁺ exhausted T cells (CD8⁺ Tex) (Fig. 2D, Additional file 6: Table S4). By quantifying the subgroup preferences for each ICEP, we found that the subgroups co-expressing PDCD1, CTLA4, TIGIT, LAG3, and CD38 were mainly enriched in CD8⁺ Tex, followed by CD8⁺ GMZK⁺ Teff and CD8⁺ GMZH⁺ Teff. In contrast, the subgroups expressing KLRC1 alone or co-expressing with other immune checkpoints were most highly enriched in CD8⁺ KLRC1⁺ Trm, followed by CD8⁺ IL7R⁺ Tm (Fig. 2E).

In summary, our findings indicated that PD-1, CTLA-4, TIGIT, LAG-3, and CD38 showed a co-occurrence pattern in CD8⁺ T cells, possibly associated with their shared characterization of CD8⁺ T cell exhaustion. Interestingly, NKG2A demonstrated a distinct mutual exclusive expression pattern compared to other immune checkpoints, showing a preference for either singular expression or co-expression with another immune checkpoint. Based on these immune checkpoint expression characteristics, we divided CD8⁺ T cells into two groups: one group co-expressing PD-1, CTLA-4, TIGIT, LAG-3, or CD38 associated with exhausted T cells, which we named immune checkpoint expression pattern (ICEP) 1 CD8⁺ T cells, and the other group expressing NKG2A alone or co-expressing with other immune checkpoints, which we named ICEP2 CD8⁺ T cells.

Immune checkpoint co-expression patterns are associated with distinct differentiation fates of CD8⁺ T cells

We inferred the potential differentiation trajectory of CD8⁺ T cells using Monocle3. The results showed that CD8⁺ Tn cells were at the starting point of the CD8⁺ T cell differentiation trajectory, subsequently branching into two paths. Branch 1 progressed from CD8⁺ Tn to CD8⁺ GMZK⁺ Teff, CD8⁺ GMZH⁺ Teff, and finally to CD8⁺ Tex. Branch 2 followed the path from CD8⁺ Tn to CD8⁺ IL7R⁺ Tm, CD8⁺ KLRC1⁺ Trm, and ultimately to CD8⁺ Tex (Fig. 3A, B). Additionally, we found that the expression of PDCD1, CTLA4, and TIGIT initially decreased during the differentiation from CD8⁺ Tn to CD8⁺ KLRC1⁺ Trm, but reached a peak during the differentiation towards CD8⁺ Tex. In contrast, CD38 and LAG3 exhibited increased expression during the differentiation from CD8⁺ Tn to CD8⁺ GMZK⁺ Teff and CD8⁺ GMZH⁺ Teff, followed by a decrease in expression in CD8⁺ KLRC1⁺ Trm, and a subsequent re-elevation during the differentiation towards CD8⁺ Tex. Conversely, KLRC1 expression gradually increased during the differentiation from CD8⁺ Tn to CD8⁺ KLRC1⁺ Trm but decreased during the differentiation towards CD8⁺ Tex. VSIR exhibited a relatively stable expression pattern throughout the CD8⁺ T cell differentiation process (Fig. 3C). Therefore, we hypothesized that CD8⁺ T cells in gastric cancer patients undergo two distinct differentiation fates, with CD8⁺ T cells on branch 1 overlapping with ICEP1 CD8⁺ T cells, and CD8⁺ T cells on branch 2 overlapping with ICEP2 CD8⁺ T cells, suggesting that the two ICEPs may represent two distinct differentiation fates of CD8⁺ T cells.

Fig. 3

Fate determination and transcriptional regulation mechanisms of gastric cancer CD8⁺ T cells differentiation. A Monocle3 analysis identifying the differentiation trajectory of CD8⁺ T cells, with colors corresponding to pseudotime. B Along the pseudotime axis, density distribution plot of CD8⁺ T cells subpopulations by position. C Dynamic visualization of immune checkpoint expression along the pseudotime analysis axis. D Scatter plot illustrating the transcriptional regulatory network of (left) CD8⁺ Tex subpopulation and (right) CD8⁺ KLRC1⁺ Trm subpopulation. E Box plot showing the expression of transcription factors MSC and RUNX3. F CellOracle simulation depicting the cell fate vector plot after (left) MSC and (right) RUNX3 knockout; *P-value < 0.05; **P-value < 0.01; ***P-value < 0.001; ****P-value < 0.0001

To delineate the upstream regulatory mechanisms of the two ICEPs in CD8⁺ T cells, we analyzed the transcriptional regulatory factors in different CD8⁺ T cell subgroups using the R package SCENIC. The results showed that CD8⁺ Tex cells were highly enriched for transcription factors related to T cell differentiation, such as RORC, NR1D1, PPARG, and MSC [18,19,20] (Fig. 3D). Previous studies have found that the transcription factor MSC can promote the differentiation of regulatory T cells (Tregs), which highly express various immune checkpoints [21]. Therefore, we hypothesized that MSC might also be involved in the regulation of immune checkpoint expression. Further analysis of the TCGA-STAD cohort revealed a positive correlation between MSC expression and the expression of PDCD1, TIGIT, CD38, and LAG3 (Additional file 7: Fig. S3A). Moreover, MSC was relatively highly expressed in the responder group (Fig. 3E), and its expression dynamics gradually increase from CD8⁺ Tn to CD8⁺ Tex during the differentiation process, reaching a peak at CD8⁺ Tex. (Additional file 7: Fig. S3B). In contrast, CD8⁺ KLRC1⁺ Trm cells were highly enriched for RUNX3 (Fig. 3D). Studies have reported that RUNX3 deficiency leads to increased expression of LAG3, TIM3, and CTLA4. Furthermore, RUNX3 has been shown to programmatically control the tissue-resident properties of CD8⁺ T cells in tumors [22], consistent with our findings. Further results showed a positive correlation between RUNX3 and KLRC1 expression (Additional file 7: Fig. S3A), and RUNX3 was relatively highly expressed in the non-responder group (Fig. 3E). Its expression dynamics were also similar to those of KLRC1, gradually increasing during the differentiation from CD8⁺ Tn to CD8⁺ KLRC1⁺ Trm but decreasing during the differentiation towards CD8⁺ Tex (Additional file 7: Fig. S3B).

CellOracle is a machine learning-based transcription factor perturbation model that can simulate cell identity transitions by combining scRNA-seq data with gene regulatory networks (GRNs). To validate the transcriptional regulatory roles of MSC and RUNX3 in the differentiation of CD8⁺ T cells, we utilized the human promoter GRN provided by CellOracle to predict the global gene expression changes in CD8⁺ T cell subgroups upon knockout of MSC and RUNX3, respectively. According to the vector field flow results predicted by CellOracle, the knockout of MSC led to a change in the vector direction of CD8⁺ T cells on branch 1, with most vectors pointing towards CD8⁺ T cells on branch 2, suggesting that the differentiation of ICEP1 CD8⁺ T cells was potentially restricted upon MSC knockout. Conversely, the knockout of RUNX3 resulted in a change in the vector direction of CD8⁺ T cells on branch 2, indicating that the differentiation of ICEP2 CD8⁺ T cells might be restricted upon RUNX3 knockout (Fig. 3F). These results suggest that MSC and RUNX3 may play a role in the differentiation process of CD8⁺ T cells exhibiting the two ICEPs.

In summary, we discovered that the two ICEPs represent two distinct differentiation fates of CD8⁺ T cells. ICEP1 CD8⁺ T cells represent the branch from CD8⁺ Tn to CD8⁺ GMZK⁺ Teff, CD8⁺ GMZH⁺ Teff, and finally to CD8⁺ Tex. ICEP2 CD8⁺ T cells represent the branch from CD8⁺ Tn to CD8⁺ IL7R⁺ Tm, CD8⁺ KLRC1⁺ Trm. Furthermore, MSC and RUNX3 may participate in regulating the differentiation of CD8⁺ T cells exhibiting these two ICEPs, offering new understandings into the upstream regulatory mechanisms of immune checkpoint expression.

NKG2A mutual exclusivity relates to anti-PD-1 resistance in gastric cancer

To investigate whether ICEP mediates resistance to anti-PD-1 therapy in gastric cancer, we assessed the disparity in the proportion of different ICEPs between the non-responder and responder groups. The results showed that the ICEP with the highest ranked difference was the CD8⁺ T cells with high expression of KLRC1 alone, while the pattern with the lowest ranked difference was the CD8⁺ T cells co-expressing high levels of PDCD1, CTLA4, TIGIT, and LAG3 (Fig. 4A). Subsequently, when comparing the expression differences of the seven immune checkpoints and co-inhibitory molecules PDCD1, CTLA4, LAG3, TIGIT, KLRC1, VSIR, and CD38 before and after anti-PD-1 therapy, as well as across different treatment outcomes, we observed that regardless of the time point, the expression of PDCD1, CTLA4, and TIGIT was significantly higher in T cells from the responder group. In contrast, LAG3 expression did not differ significantly between the non-responder and responder groups, consistent with the finding from the NCT03662659 clinical trial that the combination of LAG-3 and PD-1 immune checkpoint inhibitors with chemotherapy did not significantly improve the efficacy in patients with metastatic gastric or gastroesophageal junction adenocarcinoma [2]. Notably, KLRC1 expression was significantly elevated in T cells from the non-responder group at both time points (Additional file 9: Fig. S4A), and this pattern was also observed in CD8⁺ T cells (Fig. 4B). We employed the TIDE score to predict the TCGA-STAD cohort and compared the expression differences of KLRC1 between the immunotherapy-responsive and non-responsive groups. The results indicate that KLRC1 is relatively higher expressed in the non-responder group, highlighting its clinical relevance to the efficacy of gastric cancer immunotherapy (Additional file 9: Fig. S4B). Furthermore, we found that CD8⁺ Tex cells were highly enriched in the responder group, while CD8⁺ KLRC1⁺ Trm cells were highly abundant in the non-responder group at both time points (Additional file 9: Fig. S4C). These results suggest that ICEP2, characterized by NKG2A, in gastric cancer may contribute to the resistance to anti-PD-1 therapy.

Fig. 4

CD8⁺ T cells with mutual exclusivity pattern relates to anti-PD-1 resistance in gastric cancer. A Bar plot showing the intergroup differences in ICEPs in CD8⁺ T cells. B Box plot displaying the intergroup differences in KLRC1 expression in CD8⁺ T cells. C Subpopulations of spatial transcriptomic data. D Spatial localization of tumor cells using Cell2location on the samples, with color intensity representing the infiltration abundance. E Division of samples into intratumoral and peritumoral regions based on the tumor cell localization obtained from Cell2location. F Box plot showing the expression values and statistical differences of PDCD1, CTLA4, TIGIT, and KLRC1 between the peritumoral and intratumoral regions. G Expression values of immune checkpoint markers PDCD1, CTLA4, TIGIT, as well as CD8⁺ T cell, CD4⁺ T cell, ICEP2, and ICEP1 in each clustered subpopulation. H Subdivision of regions into high-density, medium-density, low-density, and peritumoral areas based on tumor cell infiltration levels. I Box plot illustrating the cytotoxicity scores of cells in different regions; *P-value < 0.05; **P-value < 0.01; ***P-value < 0.001; ****P-value < 0.0001

High expression of NKG2A alone does not necessarily indicate stronger immune suppression. Instead, NKG2A-expressing NK cells and γδ2 T cells exhibit immune hyperreactivity, representing enhanced cytotoxic ability. When the NKG2A ligand HLA-E is highly expressed in tumors and binds to NKG2A on NK cells and γδ2 T cells, their immune hyperreactivity is suppressed [23]. To elucidate how ICEP2 CD8⁺ T cells mediate resistance to anti-PD-1 therapy in gastric cancer, we performed spatial transcriptomics on the gastric tumor tissue of a patient with PD after anti-PD-1 treatment. We first performed quality control, normalization, and dimensionality reduction clustering on the spatial transcriptomics sample using the R package Seurat, resulting in the identification of 14 cell subgroups (Fig. 4C). We then conducted a joint analysis of the single-cell and spatial transcriptomics datasets using the Cell2location algorithm to map the identified cell types to their spatial locations within the tissue. To understand the spatial expression heterogeneity of the immune checkpoints associated with anti-PD-1 therapy response, we annotated subgroups 1, 8, 4, 9, 11, 5, 2, and 13 as intratumoral regions based on the spatial enrichment of tumor cells, while subgroups 0, 3, 6, 7, 10, and 12 were annotated as peritumoral regions (Fig. 4D, E). The results showed that the expression of PDCD1 and TIGIT did not differ significantly between the intratumoral and peritumoral regions. Notably, CTLA4 expression was higher in the intratumoral region compared to the peritumoral region, which may be related to the infiltration of CD4⁺ T cells in subgroup 1 that highly express CTLA4 (Fig. 4G). Additionally, we found that KLRC1 expression was higher in the tumor region (Fig. 4F).

To investigate whether the function of ICEP2 CD8⁺ T cells is suppressed by tumor cells, we first defined subgroup 1 as the intratumoral core region based on the degree of tumor cell infiltration. We then extracted the CD8⁺ T cell subgroup based on the marker genes and further divided it into high-density, medium-density, low-density, and peritumoral regions according to the degree of tumor cell infiltration (Fig. 4G, H). Using the Ucell algorithm, we calculated the cytotoxicity scores for these regions (Additional file 8: Table S5). We found that the cytotoxicity score of CD8⁺ T cells in the high-density tumor region was the lowest, although not significantly different from the medium-density region. Interestingly, the cytotoxicity score in the low-density tumor region was significantly higher than the high-density and medium-density regions. This suggests that as ICEP2 CD8⁺ T cells get closer to tumor cells, their functionality may be more strongly suppressed. However, the cytotoxicity score in the peritumoral region was lower than that in the low-density tumor region (Fig. 4I). We speculate that this may be related to the activation of CD8⁺ T cells to a certain extent by tumor cell antigens.

In summary, we observed that CD8⁺ T cells with mutual exclusivity pattern characterized by NKG2A were enriched in the core region of the tumor. Hence, ICEP2 CD8⁺ T cells might undergo inhibition via NKG2A-CD94/HLA-E signaling, a pathway that anti-PD-1 therapy fails to reverse, thereby maintaining the functional suppression of these CD8⁺ T cells.

The immune infiltration of ICEP2 CD8+ T cells is associated with LGMN⁺ macrophages

To explore how the TME contributes to shaping the mutually exclusive expression patterns of immune checkpoints, we used our own scRNA-seq dataset as a training set and employed the deconvolution algorithm CIBERSORTx to calculate the abundance of TME cell types in the TCGA-STAD cohort. To determine the relationships between different cell types, we performed Pearson correlation analysis on the infiltration levels of 13 cell types (Pearson correlation coefficient Rs > 0.3 and P-value < 0.05 were considered significant). We observed that the infiltration of myeloid cells in the TME exhibited the strongest positive correlation with the infiltration of CD8⁺ KLRC1⁺ Trm cells (Fig. 5A), which may be related to the recruitment function of myeloid cells.

Fig. 5

Correlation between ICEP2 CD8 + T cell infiltration and LGMN⁺ macrophage infiltration. A Analysis of correlation between CD8⁺ T cell subpopulations and infiltration of other cell types in the TME in the TCGA-STAD cohort (Pearson correlation coefficient Rs > 0.3 and P-value < 0.05 considered significant, Pie chart showing correlations with a P-value < 0.05). B UMAP plot annotating myeloid cell subsets in gastric cancer patients. C Heatmap displaying the distribution of myeloid cell subpopulations based on the ratio of observed cell counts to expected cell counts (Ro/e) calculated using a chi-square test. The top bar plot shows cell composition, and the right bar plot shows the total number of cells in each subsets. D Violin plot illustrating the functional characteristic scores of monocyte-macrophage subsets. E Analysis of correlation between CD8⁺ T cell subsets and infiltration of other cell types in the TME in the TCGA-STAD cohort. F Kaplan–Meier curve demonstrating the relationship between LGMN⁺ macrophage infiltration and overall survival (OS) of TCGA-STAD patients

We further extracted myeloid cells and reclassified them into eight subtypes through re-clustering and annotation, including two monocyte clusters, three tumor-associated macrophage (TAM) clusters, and three dendritic cell (DC) clusters (Fig. 5B, Additional file 6: Table S4). Among the TAM clusters, the LGMN⁺ macrophages, which had the highest anti-inflammatory gene set score, were most enriched in the Pre-NR group, but their enrichment levels were similar between the Post-R and Post-NR groups. This may be related to the inclusion of two gastric cancer liver metastasis samples in the Post-R group. Previous studies have shown that LGMN can activate pro-MMP-2 to degrade extracellular matrix components and promote tumor invasion and metastasis [24]. Therefore, we speculate that in the gastric cancer liver metastasis samples from the Post-R group, the infiltration of LGMN⁺ macrophages was higher, leading to similar enrichment levels between the Post-R and Post-NR groups (Additional file 11: Fig. S5A). Additionally, the IL1B⁺ macrophages, which highly expressed IL1B, VEGFA, CCL20, and had the highest pro-angiogenic gene set score, were most enriched in the Pre-R group but had lower enrichment levels in both the Post-R and Post-NR groups, suggesting that the infiltration of this subgroup may be suppressed by anti-PD-1 therapy. The TREM2⁺ macrophages were particularly enriched in the Post-R group, while their enrichment levels were lower in the Pre-R and Pre-NR groups, indicating that the infiltration of this subgroup may be influenced by anti-PD-1 therapy. Furthermore, among the identified DC subtypes, the cDC1 cells highly expressing CLEC9A were highly enriched in the Post-NR group but had relatively lower proportions in the Pre-R and Pre-NR groups. Similarly, the cDC2 cells highly expressing CD1C also had lower proportions in the Pre-R and Pre-NR groups. Mature cDCs highly expressing LAMP3, CD274 (PD-L1), and IDO1 were enriched in both the Pre-NR and Post-NR groups (Fig. 5C).

Using the R package Ucell, we scored the functional profiles of macrophage subgroups (Additional file 10: Table S6). The results showed that the LGMN⁺ macrophages had the highest anti-inflammatory score, indicating their immunosuppressive nature, consistent with the findings of Lizhi Pang et al. in glioblastoma [25]. The TREM2⁺ macrophages also had higher anti-inflammatory scores than the IL1B⁺ macrophages. In contrast, the TREM2⁺ macrophages had the highest pro-inflammatory score, followed by the IL1B⁺ macrophages, while the LGMN⁺ macrophages had the lowest score. Additionally, the IL1B⁺ macrophages exhibited the strongest pro-angiogenic ability (Fig. 5D).

Furthermore, we utilized the R package Monocle3 to analyze the developmental trajectories, state transitions, and differentiation processes of myeloid cells. By performing pseudotime analysis on the myeloid cell subgroups, we inferred their potential developmental trajectories. The results showed that monocytes were at the beginning of the trajectory, followed by the divergence of macrophages into two branches. Branch 1 proceeded from monocytes to TREM2⁺ macrophages and ultimately to LGMN⁺ macrophages. Branch 2 went from monocytes to IL1B⁺ macrophages (Additional file 11: Fig. S5B, C). Overall, the LGMN⁺ macrophages exhibited the strongest immunosuppressive capability among all myeloid cell subgroups, while the TREM2⁺ macrophages may represent an intermediate stage in the differentiation from monocytes to LGMN⁺ macrophages.

To comprehensively analyze the infiltration relationships between the two ICEPs of CD8⁺ T cells and other cell types, we re-clustered and annotated the remaining B cells, plasma cells, fibroblasts, and endothelial cell subgroups (Additional file 11: Fig. S5D, Additional file 6: Table S4). We then used CIBERSORTx to evaluate the infiltration levels of all cell clusters identified from our scRNA-seq data in the TCGA-STAD cohort. Subsequently, we performed Pearson correlation analysis on the estimated infiltration levels of the cell subgroups in TCGA-STAD. The results showed that the infiltration of the core subgroup CD8⁺ KLRC1⁺ Trm of ICEP2 CD8⁺ T cells was significantly positively correlated with the infiltration of LGMN⁺ macrophages (Pearson correlation coefficient Rs = 0.59), while the infiltration of the core subgroup CD8⁺ Tex of ICEP1 CD8⁺ T cells did not exhibit a correlation with LGMN⁺ macrophages (Fig. 5E). Additionally, we found that the immune infiltration of LGMN⁺ macrophages was associated with poor prognosis in TCGA-STAD patients (Fig. 5F), while the immune infiltration of CD8⁺ KLRC1⁺ Trm cells was associated with favorable prognosis (Additional file 11: Fig. S5E), further supporting the traits of LGMN⁺ macrophages and the reactive nature of NKG2A⁺ CD8⁺ T cells, which undergo functional suppression upon binding to the HLA-E ligand.

ICEP2 CD8 +T cells may be recruited by LGMN+Macrophages via CXCL16-CXCR6

To further elucidate the reasons behind the correlated immune infiltration of LGMN⁺ macrophages and CD8⁺ KLRC1⁺ Trm cells, we inferred the communication network between these subgroups in our scRNA-seq data using NicheNet. Interestingly, when comparing the communication differences between monocyte, macrophage and CD8⁺ T cell subgroups across treatment outcome groups, we found that in the non-responder group, there was a stronger specific CXCL16-CXCR6 signaling pathway between LGMN⁺ macrophages and CD8⁺ KLRC1⁺ Trm cells compared to the responder group. This may suggest that in anti-PD-1 therapy non-responder patients, LGMN⁺ macrophages recruit ICEP2 CD8⁺ T cells through the enhanced CXCL16-CXCR6 signaling pathway (Fig. 6A).

Fig. 6

LGMN⁺ macrophages recruit ICEP2 CD8⁺ T cells through CXCL16-CXCR6. A NicheNet predicts the ligand-receptor activities within cells and visually displays the differences in interactions between the non-response and response groups to anti-PD-1 treatment. B Spatial localization of LGMN⁺ macrophages using Cell2location on the samples, with color intensity representing the infiltration abundance. C Chord diagram displaying the network of CXCL signaling pathways between subpopulations in the spatial transcriptomic atlas. D Scatter plot showing the expression correlation between LGMN and CXCL16 in the TCGA-STAD cohort. E Scatter plot illustrating the expression correlation between RUNX3 and CXCR6 in the TCGA-STAD cohort

We analyzed spatial transcriptomic sections from PD patients to validate the spatial relationship and interaction between CD8⁺ KLRC1⁺ Trm cells and LGMN⁺ macrophages. We found that subgroups 6 highly expressed PDCD1, CTLA4, TIGIT, and CD8⁺ T cell marker genes, representing a potential enrichment of ICEP1 CD8⁺ T cells. In contrast, subgroups 9 highly expressed KLRC1, KLRD1, and CD8⁺ T cell genes, representing a potential enrichment of ICEP2 CD8⁺ T cells (Fig. 4G). Simultaneously, we observed that LGMN⁺ macrophages were enriched in subgroups 1, which was spatially proximal to the highly infiltrated ICEP2 CD8⁺ T cell subgroups 9 (Fig. 6B). Further analysis using Cellchat revealed intercellular communication between subgroups 1 as the signaling sender and subgroups 9 through the CXCL signaling pathway (Fig. 6C). Additionally, correlation analysis in the TCGA-STAD cohort showed a significant positive correlation between LGMN and CXCL16 expression (Fig. 6D). Therefore, we hypothesize that LGMN⁺ macrophages may recruit ICEP2 CD8⁺ T cells to the tumor core through the CXCL16-CXCR6 signaling pathway. Moreover, correlation analysis in the TCGA-STAD cohort revealed a significant positive correlation between RUNX3 and CXCR6 expression (Fig. 6E), and previous literature has reported that RUNX3 deficiency can significantly downregulate the expression of CCR3 and CCR5 from the CCR family [22]. Consequently, we speculate that RUNX3 not only promotes the differentiation of ICEP2 CD8⁺ T cells but also induces the expression of CCR6 from the CCR family, ultimately enabling the recruitment of ICEP2 CD8⁺ T cells by LGMN⁺ macrophages to the tumor core.

In conclusion, our findings indicated that LGMN⁺ macrophages might recruit ICEP2 CD8⁺ T cells to the tumor core via the CXCL16-CXCR6 axis, contributing to anti-PD-1 resistance in gastric cancer.

Discussion

The combination of immune checkpoint inhibitors provides a new direction for the treatment of advanced gastric cancer patients who are resistant to anti-PD-1 inhibitors. However, the current lack of comprehensive understanding of the expression patterns and intrinsic associations of immune checkpoints poses challenges in selecting combination regimens for immune checkpoint inhibitors. Additionally, previous studies have primarily focused on immune checkpoints related to exhausted T cells, such as PD-1, CTLA-4, TIGIT, and LAG-3, while the impact of other immune checkpoint expressions on CD8⁺ T cell function remains relatively unexplored. In this study, based on a comprehensive analysis utilizing multiple linear regression interaction equations, scRNA-seq, TCGA cohorts, and spatial transcriptomics, we comprehensively analyzed the intrinsic associations and expression characteristics of immune checkpoints in gastric cancer CD8⁺ T cells and revealed the ICEPs associated with anti-PD-1 therapy resistance in gastric cancer.

To elucidate the impact of co-expression of two immune checkpoints on CD8⁺ T cell function, we first discovered through multiple linear regression interaction equations that the immune checkpoints KLRC1 and PDCD1 could synergistically inhibit the proliferative function of CD8⁺ T cells. Subsequently, through co-expression analysis at the single-cell level, we found a high proportion of co-expression among the immune checkpoints PDCD1, CTLA4, TIGIT, LAG3, and CD38, suggesting a co-occurrence expression pattern among these immune checkpoints. Notably, KLRC1 exhibited a low co-expression proportion with these immune checkpoints and was primarily expressed in CD8 +KLRC1⁺ Trm cells, exhibiting tissue-resident expression characteristics, displaying mutual exclusivity with the expression of PDCD1, CTLA4, TIGIT, LAG3, and CD38. Consequently, based on the co-occurrence and mutual exclusivity patterns of immune checkpoint expression, CD8⁺ T cells were classified into two ICEPs: ICEP1 CD8⁺ T cells highly expressing the exhaustion-associated immune checkpoints PDCD1, CTLA4, TIGIT, LAG3, and CD38, and ICEP2 CD8⁺ T cells highly expressing the KLRC1 immune checkpoint. Thus, our study revealed two expression patterns of immune checkpoints in CD8⁺ T cells.

The expression of PDCD1, CTLA4, TIGIT, LAG3, and CD38 represents the exhausted state of CD8⁺ T cells, and our study consistently found that the core subcluster of ICEP1 CD8⁺ T cells was CD8⁺ Tex. In contrast, the core subcluster of ICEP2 CD8⁺ T cells highly expressing KLRC1 was CD8 +KLRC1⁺ Trm. Through pseudotime analysis with monocle3, we discovered that the two ICEPs represent two distinct differentiation fates of CD8⁺ T cells. ICEP1 CD8⁺ T cells represent the branch from CD8⁺ Tn to CD8⁺ GMZK⁺ Teff, CD8⁺ GMZH⁺ Teff, and finally to CD8⁺ Tex (branch 1). ICEP2 CD8⁺ T cells represent the branch from CD8⁺ Tn to CD8⁺ IL7R +Tm, CD8⁺ KLRC1⁺ Trm, and finally to CD8⁺ Tex (branch 2). The regulation of different transcription factors and epigenetics not only influences T cell differentiation, but also affects T cell functional states by modulating immune checkpoint expression. Previous studies have shown that the deficiency of cBAF can promote the differentiation of T cells towards the memory T cell lineage, thereby improving cancer immunotherapy [26]. On the other hand, the transcriptional factor network of IRF4, BATF, NFATc1, NR4A, TOX, and TCF1 may regulate T cell exhaustion by modulating PD-1 expression [27,28,29,30,31]. Through transcription factor analysis and CellOracle gene perturbation analysis, we found that MSC and RUNX3 may play a role in the differentiation of CD8⁺ Tex and CD8⁺ KLRC1⁺ Trm, respectively. Specifically, MSC not only promotes the differentiation of CD8⁺ Tex but also exhibits a positive correlation with PDCD1, TIGIT, CD38, and LAG3. Simulated knockdown of MSC inhibited the differentiation of ICEP1 CD8⁺ T cells. In contrast, RUNX3 promotes the differentiation of CD8⁺ KLRC1⁺ Trm and exhibits a positive correlation with KLRC1 expression. Simulated knockdown of RUNX3 may inhibit the differentiation of ICEP2 CD8⁺ T cells. Therefore, we hypothesize that MSC and RUNX3 may participate in regulating the differentiation of CD8⁺ T cells with two ICEPs, providing new insights into the upstream regulatory mechanisms of immune checkpoint expression.

The core subcluster of ICEP1 CD8⁺ T cells is CD8⁺ Tex, which exhibits high infiltration in the responder group of anti-PD-1 therapy for gastric cancer. Therefore, we speculate that ICEP1 CD8⁺ T cells are more likely to benefit from anti-PD-1 immune checkpoint inhibitor treatment. Additionally, the expression of PD-1, CTLA4, and TIGIT immune checkpoints was relatively high in T cells from both the Pre-R and Post-R groups. Previous studies have shown that the CD28 co-stimulatory signal, which plays a crucial role in T cell activation, is primarily transduced through binding to CD80/86. CTLA-4, a competitive inhibitory receptor for CD28, can also bind to CD80/86, exerting an immunosuppressive effect by attenuating the CD28 signal. Similarly, another T cell co-stimulatory molecule, CD226, transduces activation signals by binding to PVR and PVRL2, but TIGIT, its inhibitory receptor, can also bind to these two ligands, suppressing CD226 signaling through a competitive mechanism. Notably, PD-1 can indirectly inhibit the signal transduction of these two co-stimulatory molecules by recruiting and inactivating the intracellular phosphatase Shp2, thereby impairing its ability to phosphorylate CD28 and CD226 [4]. Therefore, we speculate that anti-PD-1 immune checkpoint inhibitors may enhance the phosphorylation of CD28 and CD226, increasing their competitive ability against CTLA-4 and TIGIT, and ultimately reversing the immunosuppressive effects of CTLA-4 and TIGIT to a certain extent. Additionally, we found that LAG-3, another marker of CD8⁺ T cell exhaustion, exhibited no significant difference between the non-responder and responder groups. Consistently, the NCT03662659 clinical trial found no significant improvement in efficacy when combining LAG-3 and PD-1 immune checkpoint inhibitors with chemotherapy in patients with metastatic gastric or gastroesophageal junction adenocarcinoma [2]. Thus, the application of LAG-3 immune checkpoint inhibitors in advanced gastric cancer warrants further exploration.

CD8⁺ T cells expressing NKG2A possess tissue-resident characteristics and inherent hyperreactivity. When NKG2A forms a heterodimeric complex with CD94 and binds to HLA-E on tumor cells, its hyperreactivity is suppressed. Our study results showed that the infiltration of the core subcluster CD8⁺ KLRC1⁺ Trm of ICEP2 CD8⁺ T cells was associated with favorable prognosis in the TCGA-STAD cohort, which may be related to their inherent hyperreactivity. Moreover, our single-cell analysis of gastric cancer anti-PD-1 treatment cohorts revealed higher infiltration levels of ICEP2 CD8⁺ T cells in the Pre-NR and Post-NR groups of patients. Spatial transcriptomic results from patients with progressive disease after anti-PD-1 therapy further demonstrated that, compared to ICEP1 CD8⁺ T cells, ICEP2 CD8⁺ T cells exhibited better infiltration into the intratumoral regions, and their cytotoxic ability appeared to weaken as their contact with tumor cells increased, reflecting the spatial mutual exclusivity between NKG2A and other immune checkpoints. Therefore, we hypothesize that the function of ICEP2 CD8⁺ T cells may be suppressed upon binding to HLA-E on tumor cells.

Immunotherapy resistance in gastric cancer patients is mediated by the interplay of multiple cells in the TME. Simultaneously, the TME harbors various immunosuppressive cells, such as TAM and Treg cells, which promote tumor evasion and growth by inhibiting the activity of immune cells. Therefore, to explore how the TME influences NKG2A, we analyzed cell–cell communication and found that, compared to other cells, LGMN⁺ macrophages exhibited enhanced chemoattractant CXCL16-CXCR6 signaling with CD8⁺ KLRC1⁺ Trm, suggesting a stronger recruitment of ICEP2 CD8⁺ T cells by LGMN⁺ macrophages. Furthermore, at both the bulk transcriptomic and spatial transcriptomic levels, we observed an infiltration correlation and spatial co-localization between these two cell types, further corroborating that LGMN⁺ macrophages may recruit ICEP2 CD8⁺ T cells into the tumor region via the chemoattractant CXCL16.

Legumain (LGMN) is an enzyme involved in various cellular processes that cleaves peptide bonds at the C-terminal end of asparaginyl residues. This enzyme plays roles in several crucial biological processes, including osteoclast formation, antigen processing, renal function, and brain development [32, 33]. LGMN is widely expressed in various cancers, such as breast, colon, lung, gastric, lymphoma, melanoma, and brain cancers, and its expression correlates with poor prognosis [34,35,36,37,38]. Therefore, LGMN is considered a therapeutic target for these cancers. In our study, the macrophage function gene set scoring results indicate that LGMN⁺ macrophages exhibit strong immunosuppressive properties, consistent with previous reports. Additionally, studies have shown that LGMN, as a highly expressed protease in TAMs, is transcriptionally upregulated by hypoxia-inducible factor 1-α (HIF1α) in the hypoxic TME of glioblastoma. Moreover, LGMN expression is crucial for the immunosuppressive polarization of tumor-associated macrophages [25]. Our study reveals that LGMN⁺ macrophages recruit ICEP2 CD8⁺ T cell infiltration through the chemoattractant CXCL16-CXCR6 signaling axis.

Undeniably, there are several limitations in this study. Firstly, the number of single-cell sequencing samples in the study is limited, particularly in the anti-PD-1 treatment non-responding and responding groups, which are not well-balanced. Thus, a larger cohort is required for validation purposes. Secondly, there is a lack of basic experimental validation regarding the transcriptional regulatory mechanisms underlying the two immune checkpoint expression patterns. Further experimental studies are necessary to expand the understanding of potential biological mechanisms that contribute to the observed mutual exclusivity and co-occurrence patterns of immune checkpoint expression. Thirdly, although the recruitment of ICEP2 CD8⁺ T cells by LGMN⁺ macrophages has been identified, the reasons for the highly activated nature of this recruitment process in the anti-PD-1 non-responding group and the underlying molecular mechanisms remain unclear. To shed light on the specific mechanisms involved, further exploration and molecular biology experiments are necessary.

Conclusion

In conclusion, we identified the mutual exclusivity and co-occurrence patterns of immune checkpoints in gastric cancer CD8⁺ T cells, which classified these cells into two immune checkpoint expression modes. The ICEP2 CD8⁺ T cells, centered around NKG2A, are recruited to the tumor core by LGMN⁺ macrophages via the CXCL16-CXCR6 axis, where they bind to HLA-E on tumor cells, leading to their functional suppression. Notably, anti-PD-1 therapy cannot reverse this state. Our study unveils the two distinct characteristics of immune checkpoints in gastric cancer CD8⁺ T cells and provides strong evidence for NKG2A as a novel immunotherapeutic target in gastric cancer, offering new insights into combination strategies for immune checkpoint inhibitors in gastric cancer treatment.