Landscape of four different stages of human gastric cancer revealed by single-cell sequencing

Abstract

BACKGROUND

Gastric cancer (GC) poses a substantial risk to human health due to its high prevalence and mortality rates. Nevertheless, current therapeutic strategies remain insufficient. Single-cell RNA sequencing (scRNA-seq) offers the potential to provide comprehensive insights into GC pathogenesis.

AIM

To explore the distribution and dynamic changes of cell populations in the GC tumor microenvironment using scRNA-seq techniques.

METHODS

Cancerous tissues and paracancerous tissues were obtained from patients diagnosed with GC at various stages (I, II, III, and IV). Single-cell suspensions were prepared and analyzed using scRNA-seq to examine transcriptome profiles and cell-cell interactions. Additionally, quantitative real-time polymerase chain reaction (qRT-PCR) and flow cytometry were applied for measuring the expression of cluster of differentiation (CD) 2, CD3D, CD3E, cytokeratin 19, cytokeratin 8, and epithelial cell adhesion molecules.

RESULTS

Transcriptome data from 73645 single cells across eight tissues of four patients were categorized into 25 distinct cell clusters, representing 10 different cell types. Variations were observed in these cell type distribution. The adjacent epithelial cells in stages II and III exhibited a degenerative trend. Additionally, the quantity of CD4 T cells and CD8 T cells were evidently elevated in cancerous tissues. Interaction analysis displayed a remarkable increase in interaction between B cells and other mast cells in stages II, III, and IV of GC. These findings were further validated through qRT-PCR and flow cytometry, demonstrating elevated T cells and declined epithelial cells within the cancerous tissues.

CONCLUSION

This study provides a comprehensive analysis of cell dynamics across GC stages, highlighting key interactions within the tumor microenvironment. These findings offer valuable insights for developing novel therapeutic strategies.

Key Words: Gastric cancer; Single-cell RNA sequencing; Tumor microenvironment; Cancer progression; Cell-cell interactions

Core Tip: Our single-cell atlas of gastric cancer (GC) is capable of providing a clearer view of the connection between the tumor microenvironment and paracancerous tissues, revealing the various cell progression in GC development and providing potential solutions for inhibiting cancer progression.

INTRODUCTION

Gastric cancer (GC) is one of the most common cancers worldwide, with high mortality rates and profound impacts on human health. Recent statistics indicate that in 2020, GC ranked fourth in global cancer diagnosis rates, and second in cancer-related mortality, with 1.03 million new cases and 720000 deaths reported annually. Notably, Asia exhibits the highest prevalence of GC, accounting for 70% of global cases[1]. Unfortunately, early GC is characterized by cancerous cells that do not penetrate the submucosa or spread to lymph nodes, resulting in delayed detection until the disease advances to a more severe stages and manifests noticeable symptoms[2]. The majority of individuals with GC are typically diagnosed at advanced and metastatic stages, frequently exhibiting resistance to pharmacological interventions. Despite advancements in surgical approaches, outcomes remain suboptimal, with a grim prognosis indicated by data showing a median survival of less than one year post-treatment for GC patients[3]. Metastatic spread presents a significant obstacle in GC management, as it can disseminate to distant regions via the bloodstream, lymphatic system, and peritoneum, subsequently leading to further metastatic spread to additional organs[4]. In clinical practice, individuals with GC are typically categorized into I (well-differentiated), II (moderately differentiated), III (poorly differentiated), and IV (undifferentiated) stages[5]. Patients diagnosed with stage T3 and T4 tumors exhibit the worst prognosis relative to those at other stages[6]. Early detection and treatment of GC are essential in enhancing patient prognosis. Advancements in research have revealed promising prospects for biologically targeted therapies in ameliorating the morbidity and mortality among GC patients[3]. GC is a gastrointestinal disorder with a multifaceted and intricate development mechanism, and its pathogenesis remains incompletely elucidated. Therefore, it is imperative to uncover the molecular mechanisms and associations between various stages of GC to identify novel targets for enhancing diagnostic and therapeutic strategies.

The progress in high-throughput sequencing technology has led to the discovery of an increasing number of GC-associated proto-oncogenes. This has yielded valuable insights into both GC treatment and the identification of genes linked to metastasis or prognosis[7,8]. However, this relatively macroscopic sequencing approach masks the variations in gene expression among different cell populations and overlooks the changes in cell states during GC progression[9]. The emergence of single-cell RNA sequencing (scRNA-seq), however, has paved the way for investigating individual cell changes during cancer progression[10]. Wu et al[11] presented a comprehensive single-cell and spatial characterization of colorectal liver metastases through scRNA-seq. They identified a group of metabolically activated MRC1 + CCL18 + M2-like macrophages located at the metastatic site. Yeo et al[12] gathered glioblastoma tissues at varying stages to elucidate the changes in the immune microenvironment during glioblastoma progression. Sun et al[13] revealed the single-cell landscape of precancerous lesions in oral squamous cell carcinoma initiation through this method. These findings are essential in comprehending the entire process of cancer development and advancement, thereby offering valuable insights for devising effective therapeutic strategies for cancer treatment.

In recent years, there has been a surge in scRNA-seq research targeting GC. Notably, Li et al[14] identified a subset of cancer-associated fibroblasts with pro-invasive characteristics, which are linked to poor clinical prognoses in GC patients. Additionally, Jiang et al[15] investigated the single-cell transcriptomics of organ-specific metastasis in GC cases. The above investigations provide a biological basis for potential treatment approaches for GC. However, the specific distribution and dynamic changes of different cell populations within the tumor microenvironment during GC progression remain unclear. Therefore, this study characterized the cellular composition of human GC at four distinct stages through scRNA-seq, offering novel perspectives on GC pathogenesis.

MATERIALS AND METHODS

Tissue sample collection

This study recruited four GC patients undergoing therapeutic surgery at our hospital during a specific period. They were diagnosed with stages I, II, III, and IV GC, respectively. The samples were obtained from the primary tumor and paracancerous normal gastric tissues, at a minimum distance of 10 cm from the tumor boundary. Prior to scRNA-seq sample collection, all patients were histologically confirmed and had not received chemotherapy. The samples were classified as follows: Stage I cancer tissue (YYF-1) and paracancerous tissue (YYF-3), stage II cancer tissue (GWD-1) and paracancerous tissue (GWD-3), stage III cancer tissue (WXZ-2) and paracancerous tissue (WXZ-3), and stage IV cancer tissue (XZJ-1) and paracancerous tissue (XZJ-3). All individuals enrolled in the study provided informed consent. Approval for the research protocols was granted by the ethics committee of the Tumor Hospital Affiliated to Xinjiang Medical University (2024-96).

Single-cell suspension preparation

Fresh GC and adjacent tissues were dissected into fragments of roughly 1 mm³, followed by placing in Dulbecco's Modified Eagle Medium with 10% fetal bovine serum (FSD500, Excell Bio, Suzhou, China). Tissues were enzymatically digested at 37 °C for 30 minutes with rotation using the MACS^® Tissue Dissociation Kit (130-095-929, Miltenyi Biotec, Beijing, China), as per the instructions. Upon filtrating through a 70 μm cell strainer, cell suspension underwent 5-minute centrifugation at 400 × g. A red blood cell lysis buffer (R1010, Solarbio, Beijing, China) was employed for lysing the red blood cells. Suspension of the cells was conducted on sorting buffer consisting of phosphate-buffered saline (PBS) and 2% fetal bovine serum following two washes with PBS. Subsequently, chilled fetal bovine serum (10 mL) was introduced into the tube, followed by 10-minute centrifugation at 250 g. After discarding the supernatant, the resulting pellet was reconstituted in magnesium- and calcium-free PBS (5 mL) containing 0.04% (w/v) bovine serum albumin. The preparation of a single-cell suspension was achieved at a concentration of 700–1200 cells per microliter with viability exceeding 85%, as determined using the Countess II FL Automated Cell Counter (AMQAF1000, Thermo Fisher Scientific, Shanghai, China).

Library construction

Single-cell suspensions underwent 7-AAD Viability Staining Solution (00-6993-50, eBiocience) for flow cytometry analysis using the BD FACSAria III instrument. Subsequently, the isolation of 1 × 10⁵ viable cells was carried out, followed by transfer to 1.5 mL tubes containing sorting buffer. Then, automated counting was performed using the Countess II Automated Cell Counter. Following the manufacturer's instructions, the isolated single cells were subjected to processing via the GemCode Single Cell Platform with the GemCode Gel Bead, Chip, and Library Kits (10 × Genomics). A cell load of 1 × 10⁵ was utilized for each sample. Within the GemCode instrument, cells were encapsulated in Gel Bead-in-Emulsion for cell lysis and barcod attachment through reverse transcription. Subsequent steps involved RNA amplification, shearing, and ligation of 3' adaptors and sample indexes. The resulting libraries were sequenced on the Illumina HiSeq 4000 sequencing platform utilizing 150 bp paired-end sequencing.

scRNA-seq data processing

The gene expression matrix for individual samples was generated using CellRanger (v2.0.2) software. To ensure high-quality data, the excluded cells involved those with fewer than 200 or more than 7000 detected genes, or those with over 10% of their total reads derived from mitochondrial genes. This filtering contributes to eliminating low-quality or dying cells. Scrublet (v0.1) was employed for identifying and removing potential doublets attributed to the simultaneous capture of multiple cells. Normalization of unique molecular identifier (UMI) counts for each gene was achieved by dividing the UMI counts by the total UMI counts within each cell. Log-transformation was conducted using the NormalizeData function in Seurat, with a scale factor of 10000. Subsequently, the ScaleData function in Seurat was applied for addressing the impacts of the detected UMI counts, mitochondrial gene proportion, and cell cycle on gene expression values.

Single-cell data clustering

The Seurat v3.2.2 R package was adopted for single-cell data analysis. A canonical correlation assay was initially selected to address batch effects across experiments, followed by integration of the gene expression matrix from all samples into a unified matrix. The commonly expressed genes were identified using the FindVariableGenes function, with an average expression of 0.05-5 and a dispersion of at least 0.5, to determine the prevalent cell types in the dataset. Subsequently, visualization of dimensionality reduction was conducted using the Uniform Manifold Approximation and Projection Plotting function. The FindClusters function was applied for grouping all cells in both cancer and normal samples, with a resolution set at 0.4. The RunTSNE function was utilized for visualizing t-distributed stochastic neighbor embedding data. Additionally, cell clustering was conducted using Harmony (v1.0), and a comparison of the results with those from Seurat demonstrated consistent clustering outcomes between the two methods.

Reclustering of different cell types

The gene expression variability among samples was assessed through a two-way analysis of variance to distinguish various cell subtypes or states within homogenous cell populations. Following removing these genes from the comprehensive list of variants, a secondary clustering analysis was conducted. Clustering of the cells was performed at varying resolutions (0.4, 0.6, and 0.8), with the optimal resolution for discriminating cell subtypes or states manually determined.

Quasi-temporal analysis and screening of differentially expressed genes

Based on the Seurat version 3.2.2 analysis, a single-cell trajectory was established through Monocle 2, an R package. The process involved generating a cell data set object using the new cell data set function with parameters from the negbinomial.size function. Gene expression levels were assessed using the dispersion Table function, excluding genes expressed in fewer than 10 cells at levels below 0.001. The Find Variable Features function with the 'vst' method was adopted for identifying highly variable genes, which were provided to Monocle for inferring cell trajectories. The Wilcoxon rank-sum test was employed for identifying differentially expressed genes (DEGs) within specific clusters by comparing the cluster cells to all other cells. The Bonferroni correction was applied for adjusting the P value, with statistical significance for DEG identification set at P < 0.05. Marker genes were defined as those with over a 2-fold higher expression frequency relative to other clusters. Cell types within each cluster were determined by annotating classical marker genes and top DEGs.

Cell-cell interactions

CellPhoneDB v. 2.1.4 was utilized for analyzing cell-cell interactions. Genes expressed in fewer than 10 cells were excluded. Based on specific receptors and ligands, CellPhoneDB v2.1.4 was utilized to import the remaining data for categorization, illustrating the cell-cell interaction network. CellChat, an R package, was employed for predicting interactions between different cell types. Such a process involved creating a CellChat object, and functions like create CellChat, identify Over Expressed Genes, and identify Over Expressed Interactions were adopted for data preprocessing. Subsequently, the calculation of cell communication was conducted using the aggregate Net, with visualization achieved through methods like netVisual_circle, netVisual_aggregate, and netVisual_diffInteraction. Significance was assessed using permutation tests. The interaction groups were refined by excluding the average value of all cell subset interaction groups following initial assessment with default parameters.

Quantitative real-time polymerase chain reaction

Quantitative real-time polymerase chain reaction (qRT-PCR) was adopted for assessing the mRNA levels of cluster of differentiation (CD) 2, CD3D, CD3E, cytokeratin (KRT) 19, KRT8, and epithelial cell adhesion molecule (EPCAM) in cancer and paracancerous tissue. TRIzol Reagent (QIAGEN Invitrogen, United States) was employed for extracting total RNA, followed by reverse transcription with the SuperScript III First-Strand Synthesis System (Invitrogen, Carlsbad, CA, United States). Polymerase chain reaction was performed with SYBRR Green PCR Master Mix (Takara Bio, Japan), with GAPDH as the internal control. The 2^-△△CT method was applied for calculating the relative mRNA expression. The primer sequences employed were provided below: CD2, Forward: 5'-AGGGTCATCACACACAAGTGG-3'; Reverse: 5'-CAGGCAGAGAGCCCATTACTC-3'; CD3D, Forward: 5'-ACTCTGCTCCTTGCTTTGGG-3'; Reverse: 5'-AGGTTCACTTGTTCCGAGCC-3'; CD3E, Forward: 5'-TCTCCGTTCAGTTCCCTCCT-3'; Reverse: 5'-TACCCAGTCCATCCCCAGAG-3'; KRT19, Forward: 5'-AAATCAGTACGCTGAGGGGC-3'; Reverse: 5'-CGACCTCCCGGTTCAATTCT-3'; KRT8, Forward: 5'-CTTACCTCCCTCCATGCTGTC-3'; Reverse: 5'-TAGCCACAGGGTCCCTAAGAG-3'; EPCAM, Forward: 5′-AGCTGGCCGTAAACTGCTTT-3'; Reverse: 5'-GCCCCTTCAGGTTTTGCTCT-3'; GAPDH, Forward: 5'-TGTTGCCATCAATGACCCCTT-3', Reverse: 5'-CTCCACGACGTACTCAGCG-3'.

Flow cytometry

Following mince, cancer and adjacent tissues underwent DNase I (0.2 mg/mL) and collagenase (1 mg/mL). Upon passing through a 70 μm cell strainer, a cell concentration of 2 × 10⁷ cells/mL was achieved for the suspension. Subsequently, the cell suspension (100 μL) was dispensed into the 96-well plates, after which the respective fluorescently labeled antibody was added for incubation at a reduced temperature. A flow cytometer (Attune Flow Cytometers, Thermo Fisher Scientific, United States) was adopted for measuring the cell counts. The antibodies utilized in this process involved Alexa Fluor^® 647 Mouse monoclonal [OX34] to CD2 (ab256295, Abcam, United States), Alexa Fluor^® 488 Rabbit monoclonal [EP776(2)Y] to CD3 zeta (phospho Y83) (ab237451, Abcam, United States), Alexa Fluor^® 555 Rabbit monoclonal [EPR20532-222] to EpCAM (ab313668, Abcam, United States), and Alexa Fluor^® 488 Rabbit monoclonal [EP1580Y] to KRT19 (ab192643, Abcam, United States).

Statistical analysis

GraphPad Prism version 9.0 software (GraphPad Software, San Diego, CA, United States) was utilized for statistical analysis and graphical representation. Data were presented as mean ± SE of the mean. Data were first assessed for normality using the Shapiro-Wilk test. An unpaired t-test was employed for analyzing the normally distributed data, while the Mann-Whitney test was applied for non-parametric data. P < 0.05 represented the statistical significance threshold.

RESULTS

Global single-cell atlas of four different stages of GC

scRNA-seq was conducted on viable cells from cancer and paracancerous tissues from four cancer patients. A total of 73645 single cells were identified across the eight tissue samples, leading to the classification of these cells into 25 distinct clusters following dimensionality reduction and unsupervised cell clustering (Figure 1A). The distribution of cell cluster proportions varied across tissue samples, as revealed by subsequent analysis. Both cancerous tissues at different stages and paracancerous tissues exhibited distinct cell cluster percentages (Figure 1B), suggesting highly intricate alterations in cell cluster proportions during GC development[16]. Further annotation allocated the 25 cell clusters into 10 distinct types, including T/NK cells, B cells, plasma cells, myeloid cells, neutrophils, mast cells, epithelial cells, fibroblasts, myofibroblasts, and endothelial cells. Analysis of cell percentages across different stages revealed variations among tissue samples, particularly in the proportions of epithelial cells and T/NK cells (Figure 1C). Furthermore, the percentages of these 10 cell types in cancer and paracancerous tissues were examined for exploring the specific cellular composition and heterogeneity in both tissues. Overall, T/NK cells were evidently raised in cancerous tissues, along with a notable reduction in epithelial cells (Figure 1D). The above outcomes indicated varying distribution of cell types at various developmental stages within the tissues. Notably, significant changes were observed in the populations of T/NK cells and epithelial cells (Figure 1E). Such a phenomenon was also evident when comparing cancerous tissues with paracancerous tissues (Figure 1F). The marker proteins identified by the 10 distinct cell types are presented in Figure 1G and H. Specifically, CD2, CD3, and CD3D were indicative of T/NK cells, while KRT19, KRT8, and EPCAM represented epithelial cells.

Figure 1 Global single-cell atlas of four different stages of gastric cancer. A: Single-cell atlas of 73645 cells captured from 25 cell clusters. Tissue samples are shown on the left, and cell clusters on the right, with colors indicating cell types; B and C: Composition of cell clusters and cell types within various tissue samples (n = 8); D: Proportion of different cell types in cancer (n = 4) and para cancerous tissues (n = 4); E: Tissue sample distribution in different cell types; F: Distribution of cancer and para cancerous tissues across varying cell types; G and H: Dot Plot and Heatmap of the corresponding marker genes identifying the cell types.

Tumor cell identification and epithelial cell re-clustering

Based on the aforementioned findings, significant differences in epithelial cell proportion were observed in cancerous tissues relative to paracancerous tissues, as well as within tumor tissues at varying stages. To further characterize the distribution of epithelial cells during GC progression, malignant tumor cells were initially classified based on their malignancy scores[17]. The outcomes displayed that the vast majority of malignant tumor cells were concentrated in epithelial cells (Figure 2A). Upon identifying the tumor cells, the different cell type distributions were examined in the cancerous tissues of patients at various stages. The findings revealed a notable elevation in tumor cell proportion in stages I, II, and IV, accompanied by a remarkable decline in stage III (Figure 2B). Given these, re-clustering was conducted on epithelial cells (Figure 2C), malignant cells (Figure 2D), and non-malignant cells (Figure 2E), respectively, each exhibiting 15 cell clusters. Distinct gene expression patterns were fully exhibited by these clusters (Figure 2F). The data presented indicated the presence of 15 distinct cellular states within epithelial cells throughout GC progression. In cluster 1, a lower abundance of tumor cells exhibiting malignant characteristics was observed.

Figure 2 Tumor cell identification and epithelial cell re-clustering. A: Uniform manifold approximation and projection (UMAP) plot after cancer cell identification, with tissue samples on the left and cell types on the right; B: Proportion of different cell types in tissue samples from patients at various stages after cancer cell identification; C-E: UMAP of epithelial cell re-clustering (C, left: Tissue samples, right: Cell cluster), malignant cell re-clustering (D, left: Malignant cells, right: Cell types), and non-malignant cell re-clustering (E: left: Non-malignant cells, right: Cell types); F: Differential analysis of re-clustering between malignant tumor cells and non-malignant tumor cells.

CD4 T cell re-clustering and quasi-temporal analysis

GC progression is significantly influenced by T cells[18]. Prior results indicated the presence of multiple T cell subtypes within cancerous tissues, paracancerous tissues, and various stages of cancer. To further explore these observations, an initial examination was conducted on T cell distribution across 25 distinct cell clusters. The findings demonstrated a noticeable presence of T cells across multiple cell clusters (Supplementary Figure 1A). Among these, CD4 and CD8 T cells were found to be essential in combating tumor cells, with their distribution across different cell clusters showing clear separation without overlap (Supplementary Figure 1B and C).

CD4 T cell re-clustering (Figure 3A) revealed a markable elevation in cell clusters 2, 3, 4, and 6 in cancerous tissues as opposed to para cancerous tissues (Figure 3B). Notably, cluster 6 presented high expression of proliferation markers proliferating cell nuclear antigen and marker of proliferation Ki-67, while there was an elevation in levels of immune checkpoint markers PDCD1 and C-C motif chemokine receptor in cluster 4 and CTL4A and T-cell immunoreceptor with immunoglobulin and immunoreceptor tyrosine-based inhibitory motif domains in cluster 3 (Figure 3C). These outcomes suggested a shift of CD4 T cells toward more proliferative and exhausted states in the tumor microenvironment, potentially contributing to immune evasion by the tumor.

Figure 3 CD4 T cell re-clustering and quasi-temporal analysis. A: Uniform manifold approximation and projection of CD4 T cell re-clustering, with tissues on the left and cell types on the right; B: Proportion of different CD4 T cell clusters in cancerous tissues and para cancerous tissues; C: Dot Plot illustrating the mean expression levels of marker genes in six main identified cell clusters; D: Quasi-temporal analysis displaying quasi-temporal distribution of cell trajectories (top left), cell trajectories in differentiated states (top right), sample cell trajectories (bottom left), subset cell trajectories (bottom right); E: Scatter plot of differential gene expression in differentiation states; F: Heatmap of hierarchical clustering of branched differential genes.

The differentiation trajectories of CD4 T cells were further examined using quasi-temporal analysis. The outcomes illustrated five distinct stages of differentiation, with a progressive transition of CD4 T cells from para cancerous tissues to cancerous tissues. During stage III of differentiation, clusters 3, 4, and 6, which predominated in cancerous tissues, exhibited a pronounced clustering tendency (Figure 3D). A marked up-regulation was observed in genes, such as B-cell translocation gene 2, C-C motif chemokine ligand 20, dual specificity phosphatase (DUSP) 2, heat shock protein family D member 1, heat shock protein family E member 1, Jun proto-oncogene, AP-1 transcription factor subunit, and regulator of G protein signaling 1, with an elevation in regenerating family member 1-alpha during the initial stage of differentiation (Figure 3E). These findings suggested that CD4 T cells underwent a dynamic differentiation process as they transitioned from the para cancerous environment to the cancerous environment, potentially influencing tumor progression and immune response. Based on the hierarchical clustering heatmap of branching differential genes, these cells showed different expression states (Figure 3F). Subsequently, the top 10 branches with the most pronounced changes in expression levels were identified. Additionally, the expression trends at stages I and V differed from those observed at the other stages upon examining the expression patterns of these genes along the differentiation trajectories (Supplementary Figure 2A). Trajectory diagrams were generated to illustrate the gene expression patterns from the top 10 differentially expressed branches across various stages of cancerous tissues, para cancerous tissues, and distinct CD4 T cell subsets. These diagrams provided a quasi-temporal line of gene expression changes. Notably, a consistent decline in gene expression over time was observed in cancerous tissues at advanced stages and in specific cell clusters, as displayed by the analysis (Supplementary Figure 2B-D).

CD8 T cell re-clustering and quasi-temporal analysis

Next, re-clustering and quasi-temporal analysis were conducted on CD8 T cells, dividing them into seven distinct cell clusters (Figure 4A). As opposed to para cancerous tissues, clusters 0, 1, 2, and 6 were markedly raised in cancerous tissues, indicating their potential roles in anti-tumor immunity (Figure 4B). An elevation in lymphocyte activation gene 3 Levels was observed in cluster 6, which was highly expressed in cancerous tissues, as well as an increase in Methyl-CpG-binding domain protein 2 expression in cluster 2. Clusters 1 and 0 were marked by raised levels of CD27, granzyme B, CTL4A, and PDCD1 (Figure 4C). These findings suggested that CD8 T cells in cancerous tissues might exhibit both cytotoxic and exhausted phenotypes, potentially limiting their effectiveness in tumor clearance.

Figure 4 CD8 T cell re-clustering and quasi-timing analysis. A: Uniform manifold approximation and projection of CD8 T cell re-clustering, with tissues on the left and cell types on the right; B: Proportion of different CD8 T cell clusters in cancerous tissues and para cancerous tissues; C: Dot Plot illustrating the mean expression levels of marker genes across six main identified cell clusters; D: Quasi-temporal analysis depicting quasi-temporal distribution of cell trajectories (top left), cell trajectories in differentiated state (top right), sample cell trajectories (bottom left), and subset cell trajectories (bottom right); E: Scatter plot of differential gene expression in differentiation states; F: Heatmap of hierarchical clustering of branching differential genes.

Quasi-temporal analysis of CD8 T cells revealed three distinct stages of differentiation, transitioning from para cancerous to cancerous tissues (Figure 4D). During the initial phase of differentiation, there was a notable up-regulation in the expression levels of activating transcription factor 3, DUSP1, DUSP2, heat shock protein 90 alpha family class B member 1, heat shock protein family A (Hsp70) member (HSPA) 1A, HSPA1B, HSPA6, heat shock protein family E member 1, and Jun proto-oncogene, AP-1 transcription factor subunit, along with elevated progastricsin expression during the stage III of differentiation (Figure 4E). This progression suggested that CD8 T cells underwent functional changes when transitioning them from paracancerous to cancerous environments. Such a finding could reflect their role in the tumor immune response. The hierarchical clustering of differential genes presented different expression states in these cells (Figure 4F). Subsequently, the top 10 branches exhibiting the most significant variations in expression levels were identified. The expression patterns of these genes along the differentiation trajectories varied between stage I and the subsequent stages II and III (Supplementary Figure 3A). Trajectory diagrams were generated to demonstrate the expression patterns of various DEGs within the top 10 branches across various stages of cancerous tissues, para cancerous tissues, and different CD8 T cell subsets. These diagrams were utilized for depicting a pseudo-time line, highlighting a downward trend in expression levels within advanced-stage cancerous tissues, particularly evident in stages II and III (Supplementary Figure 3B-D).

Re-clustering of myeloid cells, endothelial cells, fibroblasts, and myofibroblasts

Given the aforementioned findings, myeloid cells, endothelial cells, fibroblasts, and myofibroblasts exhibited varying proportions in patients across various stages. Subsequently, reclassification was conducted. Notably, myeloid cells were assigned to 11 distinct cell clusters (Figure 5A), endothelial cells to nine cell clusters (Figure 5B), fibroblasts to 13 cell clusters (Figure 5C), and myofibroblasts to nine cell clusters (Figure 5D). These observations underscored the intricate involvement of these cell types in GC progression.

Figure 5 Reclassification of myeloid cells, endothelial cells, fibroblasts, and myofibroblasts. Uniform manifold approximation and projection. A: Myeloid cell; B: Endothelial cell; C: Fibroblast; D: Myofibroblast; re-clustering, with tissues on the left and cell types on the right.

Cell-cell interactions

Cellular interactions among various cell populations were analyzed at distinct time points during the specified period. B cells directly interacted with mast cells, myeloid cells, and neutrophils across different stages (Figure 6). Notably, the interactions between cells varied over time. Specifically, during the initial stage of para cancerous tissue, limited interactions were found between B cells and mast cells, as well as other cell types. However, as cancer progressed to stages II, III, and IV, the interactions between B cells and other cell types were notably increased. Furthermore, interactions between myeloid cells and neutrophils with other cell types were more frequent in para cancerous tissues, while these interactions were evidently dropped in cancerous tissues.

Figure 6 Cell-cell interactions. The rows and columns in the heatmap representing different cell types across various tissue samples, respectively. The colors indicating the degree of intercellular interaction, ranging from blue (0) to red (4), with higher values denoting stronger intercellular activity.

Increased T cells and decreased epithelial cells in GC tissues

In this study, qRT-PCR was employed for validating the findings of the scRNA-seq analysis. Specifically, assessment was conducted on the expression profiles of T cell marker genes (CD2, CD3D, and CD3E) and epithelial cell marker genes (KRT19, KRT8, and EPCAM) in cancerous tissues and para cancerous tissues. The findings demonstrated that relative to para cancerous tissues, the CD2, CD3D, and CD3E expression was evidently u-regulated in cancerous tissues, along with a notable decline in the KRT19, KRT8, and EPCAM expression (Figure 7A). Additionally, flow cytometry was adopted for quantifying the number of T cells and epithelial cells. The trends in T cells and epithelial cells in cancerous tissues were consistent with the qRT-PCR findings. The outcomes displayed that as opposed to para cancerous tissues, the CD2 and CD3 cells were markedly raised in cancerous tissues, accompanied by reduced EPCAM and KRT19 cells (Figure 7B and C).

Figure 7 Trends in T cells and epithelial cells in gastric cancer tissues. A: Quantitative real-time polymerase chain reaction for assessing the marker gene expression levels of T cells and epithelial cells; B and C: Flow cytometry for analyzing the quantities of T cell and epithelial cell marker genes. ^aP < 0.001, ^bP < 0.0001 vs para cancerous tissues.

DISCUSSION

The precise therapeutic strategy development has been hindered due to the limited understanding of GC pathogenesis, thereby impacting treatment efficacy. Therefore, there is an urgent need for a more comprehensive and detailed understanding of the biological changes within the tumor microenvironment throughout GC progression. In this study, scRNA-seq was employed for unveiling the diversity and gene expression profiles of various cell populations in GC at varying stages, delineating the trajectories and dynamics of different cell types during cancer progression.

A total of 73645 individual cells were detected across the eight samples. They were categorized into 25 distinct cell clusters, representing a higher number of cell clusters relative to prior scRNA-seq studies on GC[19,20]. This study revealed significant variations in the proportions of cell clusters between cancerous and para cancerous tissues at various stages, suggesting that interactions between cancer cells and other cellular components within the tumor microenvironment are essential in driving the malignant transformation of non-cancerous cells[21,22]. Such a finding was further supported by identifying pronounced differences in epithelial cell distribution. Additionally, malignant tumor cells were predominantly localized within the epithelial population, highlighting the potential for epithelial cells to undergo carcinogenesis during cancer progression[23,24].

The increased CD4 and CD8 T cells in cancerous tissues was a notable finding. CD4 T cells, particularly in clusters 3 and 4, exhibited up-regulation of key immune modulators, such as CTL4A and PDCD-1, which are known to inhibit T cell activation and facilitate immune evasion by tumor cells[25-28]. This finding suggests that while CD4 T cells infiltrate the tumor microenvironment, their functional exhaustion may limit their ability to effectively counteract tumor growth. Interestingly, the quasi-temporal analysis highlighted the transition of CD4 T cells from para cancerous tissues to a cancerous tissue-like state, further emphasizing the shift towards an immunosuppressive environment. These results indicate that immune checkpoint inhibitors (ICIs) targeting PD-1/PD-L1 and CTLA-4 may be valuable in reactivating these exhausted T cells and improving the immune response against GC[29,30].

CD8 T cells are the primary effector cells in the anti-tumor immune response[31,32], which were remarkably increased in cancerous tissues, particularly in clusters 0, 1, 2, and 6. While the expression of cytotoxic markers, such as granzyme B and CD27, suggested that these cells retain their cytotoxic potential, the up-regulated PDCD-1 in cluster 0 indicated an exhausted phenotype. Such outcomes were similar to those observed in CD4 T cells. Our findings highlighted the intricate role of CD8 T cells in the tumor microenvironment. In this context, they initially act as cytotoxic cells but gradually enter an exhausted state, potentially limiting their anti-tumor effectiveness. Such a finding emphasized the fact that combining ICIs with therapies aimed at enhancing CD8 T cell function could improve anti-tumor efficacy by reversing immune exhaustion and increasing the cytotoxic activity of these cells.

Apart from ICIs, this study also suggested that targeting other components of the tumor microenvironment might be beneficial. The interactions among immune cells, myeloid cells, endothelial cells, fibroblasts, and myofibroblasts were demonstrated to be essential in promoting tumor progression. These interactions might facilitate tumor survival and invasion by promoting immune suppression, angiogenesis, and malignant transformation of epithelial cells. Targeting the tumor-stroma interface could disrupt these supportive interactions, potentially preventing immunosuppression and inhibiting tumor growth. Anti-stromal therapies or inhibitors that block tumor-stroma interactions could be integrated with ICIs to provide a more comprehensive treatment approach, addressing both the immune evasion mechanisms and the supportive role of the tumor microenvironment.

Despite these findings providing significant insights into potential therapeutic strategies, several challenges must be addressed for clinical translation. Firstly, GC heterogeneity presents a challenge for broad application, as significant variations in the tumor microenvironment may exist between patients. Identifying biomarkers that predict response to ICIs or stroma-targeting therapies is critical for developing personalized treatment strategies. Determining the optimal timing and therapy sequencing presents another challenge; for instance, administering stroma-targeting agents prior to ICIs can enhance immune cell infiltration, and is crucial for maximizing treatment efficacy. Overcoming resistance to ICIs may require combining them with other immune-modulating agents or tumor microenvironment-targeting therapies. Future clinical trials should focus on these combination approaches to improve patient outcomes.

However, several limitations must be acknowledged, which may influence the generalizability and depth of our findings.

Firstly, the relatively small sample size of this study, with only four patients representing the four GC stages, may reduce statistical power and may not fully capture GC heterogeneity. Future studies should include larger patient cohorts with multiple time points to offer a more comprehensive view of the dynamic changes in the tumor microenvironment.

Secondly, all patients in this study came from a single geographical region, which may introduce potential bias related to genetic or environmental factors specific to this population. Future research should involve patients from diverse regions and backgrounds to better understand the broader applicability of our findings.

Thirdly, while qRT-PCR and flow cytometry were applied for verifying key scRNA-seq findings, additional validation methods, including immunohistochemistry and functional assays, can enhance the robustness of our conclusions.

Furthermore, despite the successful tracking of various cell types involved in GC development, the challenge of differentiating distinct cell subtypes and states based on spatial context persists. Spatial transcriptomics or other spatial mapping technologies can complement our current single-cell analysis and provide additional insights into how the tumor microenvironment interacts with specific cell types in situ.

In summary, although this study offers a foundational understanding of GC cell dynamics, it is critical for translating these findings into clinical practice and developing targeted therapies by overcoming these limitations through larger and more diverse cohorts as well as complementary methodologies.

CONCLUSION

Overall, our single-cell atlas of GC provides a clearer depiction of the connection between the tumor microenvironment and paracancerous tissues, revealing the progression of various cells in GC development and providing a potential solution for inhibiting cancer development.