Revealing the transcriptional heterogeneity of organ-specific metastasis in human gastric cancer using single-cell RNA Sequencing

Abstract

Background: Deciphering intra- and inter-tumoural heterogeneity is essential for understanding the biology of gastric cancer (GC) and its metastasis and identifying effective therapeutic targets. However, the characteristics of different organ-tropism metastases of GC are largely unknown.

Methods: Ten fresh human tissue samples from six patients, including primary tumour and adjacent non-tumoural samples and six metastases from different organs or tissues (liver, peritoneum, ovary, lymph node) were evaluated using single-cell RNA sequencing. Validation experiments were performed using histological assays and bulk transcriptomic datasets.

Results: Malignant epithelial subclusters associated with invasion features, intraperitoneal metastasis propensity, epithelial-mesenchymal transition-induced tumour stem cell phenotypes, or dormancy-like characteristics were discovered. High expression of the first three subcluster-associated genes displayed worse overall survival than those with low expression in a GC cohort containing 407 samples. Immune and stromal cells exhibited cellular heterogeneity and created a pro-tumoural and immunosuppressive microenvironment. Furthermore, a 20-gene signature of lymph node-derived exhausted CD8⁺ T cells was acquired to forecast lymph node metastasis and validated in GC cohorts. Additionally, although anti-NKG2A (KLRC1) antibody have not been used to treat GC patients even in clinical trials, we uncovered not only malignant tumour cells but one endothelial subcluster, mucosal-associated invariant T cells, T cell-like B cells, plasmacytoid dendritic cells, macrophages, monocytes, and neutrophils may contribute to HLA-E-KLRC1/KLRC2 interaction with cytotoxic/exhausted CD8⁺ T cells and/or natural killer (NK) cells, suggesting novel clinical therapeutic opportunities in GC. Additionally, our findings suggested that PD-1 expression in CD8⁺ T cells might predict clinical responses to PD-1 blockade therapy in GC.

Conclusions: This study provided insights into heterogeneous microenvironment of GC primary tumours and organ-specific metastases and provide support for precise diagnosis and treatment.

Keywords: HLA-E-KLRC1/KLRC2; gastric cancer; metastasis; single-cell RNA sequencing; tumoural heterogeneity.

FIGURE 1

Single‐cell transcriptome profiles of GC primary tumour, metastasis and adjacent non‐tumour samples. (A) scRNA‐seq and data analyses. Ten fresh human tissue samples from six patients were collected, including three primary GC samples (PT; i.e., PT1, PT2, and PT3), one adjacent non‐tumour sample (NT; i.e., NT1) and six metastasis samples (M). For the metastases, we obtained two liver tumour samples (Li; i.e., Li1, Li2), two metastatic lymph nodes samples (LN; i.e., LN1, LN2), one peritoneal tumour sample (P; i.e., P1), and one ovary tumour sample (O; i.e., O1) during gastroscopy, biopsy, or surgical resection. (B) t stochastic neighbour embedding (tSNE) projection within each patient and sample origin. (C) tSNE showing seven cell types for the 42 968 cells. (D) Heatmap of highly variable genes for seven major lineages. (E) Heatmap of marker genes for seven major lineages. (F) Proportion of each cell type in NT, PT, M, Li, LN, P, and O

FIGURE 2

Four subclusters of malignant epithelial cells in GC and their characterisation associated with metastasis. (A) tSNE of 1615 malignant and 128 non‐malignant epithelial cells. (B) tSNE of the four malignant epithelial clusters G0–G3. (C) Relative proportion of G0–G3 cells in PT (PT1, PT2, PT3), M, Li (Li1, Li2), LN (LN1, LN2), P (P1) and O (O1). (D) IPA results showing the top enriched canonical pathways in G0, G2 or G3 cells based on the DEGs. z‐score > 0 indicates that the pathway was activated; z‐score < 0 indicates that the pathway was inhibited; z‐score = 0 indicates that the pathway was neither activated nor inhibited. (E) Dot plot showing the diseases and bio functions of G0–G3 cells based on DEGs using IPA analysis. (F) Violin plot showing the expression of CD44 in G0–G3 cells. ^**** p < .0001. (G) Unsupervised transcriptional trajectory of G0–G3 subsets predicted by Monocle 2. (H–J) High level of G0‐, G1‐ and G3‐associated genes predicted poor prognosis in the TCGA‐STAD.htseq_counts.tsv dataset (n = 407 patients). Log‐rank p < .05 was considered statistically significant

FIGURE 3

Various immune responses are mediated by T and B cells during GC progression. (A) tSNE of T cells. Tregs: regulatory T cells; CD4⁺ T_EM: effector memory CD4⁺ T cells; GADD45B⁺ Th1‐like CD4⁺ T cells: GADD45B⁺ T helper type 1‐like CD4⁺ T cells; CD8⁺ T_EM: effector memory CD8⁺ T cells; MAIT: mucosal‐associated invariant T. (B) Relative proportion of each T cluster in NT, PT and M. NT: NT1; PT: PT1, PT2, PT3; M: Li1, Li2, LN1, LN2, O1, P1. (C) Violin plot showing expression of IFNG, GZMB, and KLRG1 in MAIT cells in different samples. ^* p < .05; ^*** p < .001; ^**** p < .0001; ns means p > .05. (D) Immunofluorescence staining indicates the co‐expression of CD3, TCR Vα7.2, and DAPI (nuclei) on MAIT cells in M and PT. (E) Bubble plots exhibiting significant interactions between cancer cells and T/NK cells as well as MAIT cells in M/PT and T/NK cells by the ligand–receptor pair HLE‐A‐KLRC1/KLRC2. (F) Violin plot showing the expression of PDCD1 in CD8⁺ T cells across different patients. ^**** p < .0001. (G) tSNE of B cells. (H) Relative proportion of each B cluster in NT, PT, O, P, LN, and Li. (I) IPA results showing the diseases and bio functions in T cell‐like B cells based on the DEGs. (J) Bubble plots exhibiting significant interactions between T cell‐like B cells and other cell groups by the ligand–receptor pairs HLA‐E‐KLRC1/KLRC2

FIGURE 4

Endothelial cells promote angiogenesis and create immune resistance in GC. (A) tSNE of the four endothelial cell clusters E0–E3. (B) Relative proportion of E0–E3 from NT, PT, and M. (C) Dot plot showing the diseases and bio functions of E0–E3 based on DEGs by IPA. (D) Violin plots showing the pathway scores of top enriched canonical pathways (IPA) in E0–E3 based on the DEGs. ^* p < .05; ^*** p < .001; ^**** p < .0001. (E) Bubble plots exhibiting significant interactions between endothelial cells and other cell groups by ligand–receptor pairs. (F–G) High levels of E0‐ and E3‐associated genes predicted poor prognosis in the TCGA‐STAD.htseq_counts.tsv dataset (n = 407 patients). Log‐rank p < .05 was considered statistically significant

FIGURE 5

iCAFs could be identified in the TME of GC and are associated with tumour invasion. (A) tSNE of iCAFs (F0) and mCAFs (F1). (B) Relative proportion of iCAFs and mCAFs in NT, PT, O, P, LN, and Li. (C) Heatmap of DEGs between iCAFs and mCAFs. (D) Dot plot showing iCAFs expressing PDGFRA and CXCL12 and mCAFs expressing RGS5 and ACTA2. (E) GO enrichment analysis showing the top enriched functions in iCAFs and mCAFs based on the DEGs. The functions in the red frame mean the common functions between iCAFs and mCAFs; functions in the blue frame represent unique functions of iCAFs or mCAFs. (F) Violin plot showing the expression of MMP2 and MMP14 in iCAFs and mCAFs. ^**** p < .0001. (G) Dot plot showing the expression level of growth factors across iCAFs and mCAFs. (H) High levels of iCAF‐associated genes, PDGFRA and MMP2, predicted poor prognosis in the TCGA‐STAD.htseq_counts.tsv dataset (n = 407 patients). Log‐rank p < .05 was considered statistically significant. (I) Bubble plots exhibiting significant interactions between iCAFs and other cell groups by ligand–receptor pairs of cytokines

FIGURE 6

Myeloid cells are abundant during GC progression. (A) tSNE of dendritic cells (DCs) in NT, PT, and M. (B) Counts and relative proportions of each DC cell cluster in NT, PT, O, P, LN, and Li. (C) Violin plots showing the expression of GZMB, LILRA4, LILRB4, LAMP3, CD68, CD80, and CD83 in each DC cluster. ^* p < .05; ^** p < .01; ^*** p < .001; ^**** p < .0001; ns means p > .05. (D) tSNE of macrophage cells (Ma) in NT, PT, and M samples. (E) Counts and relative proportions of each macrophage cell cluster in NT, PT, O, P, LN, and Li. (F) Violin plots showing the expression of VEGFA, SPARC, HLA‐DPB1, HLA‐DPA1, FABP5, TNF, and IL1B in macrophage clusters (Ma0‐Ma3). ^** p < .01; ^*** p < .001; ^**** p < .0001; ns means p > .05. (G) tSNE of monocyte clusters (Mo0‐Mo3) in NT, PT, and M. (H) Counts and relative proportions of each monocyte cell cluster in NT, PT, O, P, LN, and Li. (I) Violin plots showing the expression of CD14 and FCGR3A (CD16) in monocyte clusters. ^* p < .05; ^*** p < .001; ^**** p < .0001; ns means p > .05. (J) Counts and relative proportions of each neutrophil cluster in NT, PT, O, P, LN, and Li. (K) tSNE of neutrophils clusters (N0‐N5) in NT, PT, and M. (L) Violin plots showing the expression of LYZ in each neutrophil cluster. ^**** p < .0001. (M) RNAscope staining indicates the co‐expression of LYZ, NAMPT, and DAPI (nuclei) on Paneth‐like cells in PT

FIGURE 7

T and B lymphocyte subclusters vary among different organ metastases. (A) tSNE of CD8⁺ T_EM, exhausted CD8⁺ T, CD8⁺ T_EM, and Tregs coloured with metastasis type and unsupervised cluster number. (B) Heatmap of top 20 DEGs in exhausted CD8⁺ T cells in the different samples. (C) The violin plot shows a significant difference in the mean score of the 20‐gene signature from lymph node‐derived exhausted CD8⁺ T cells between the non‐lymph node metastasis (n = 17 intestinal‐histology and mixed‐histology patients) and lymph node metastasis (n = 146 intestinal‐histology and mixed‐histology patients) groups in the primary GC dataset (GSE62254). The box plots represent the median, bottom and top quantiles, whiskers correspond to 1.5× the interquartile range. Log‐rank p < .05 was considered statistically significant. (D) Prognostic significance of the 20‐gene signature in exhausted CD8⁺ T cells derived from the lymph node metastasis samples was validated in the GC dataset (GSE84437, n = 433 patients). Log‐rank p < .05 was considered statistically significant