Skip to main content
Download PDF

Single-cell N6-methyladenosine-related genes function within the tumor microenvironment to affect the prognosis and treatment sensitivity in patients with gastric cancer

Cancer Cell International volume 24, Article number: 44 (2024) Cite this article

Abstract

Background

Gastric cancer (GC) ranks fifth for morbidity and third for mortality worldwide. The N6-methyladenosine (m6A) mRNA methylation is crucial in cancer biology and progression. However, the relationship between m6A methylation and gastric tumor microenvironment (TME) remains to be elucidated.

Methods

We combined single-cell and bulk transcriptome analyses to explore the roles of m6A-related genes (MRG) in gastric TME.

Results

Nine TME cell subtypes were identified from 23 samples. Fibroblasts were further grouped into four subclusters according to different cell markers. M6A-mediated fibroblasts may guide extensive intracellular communications in the gastric TME. The m6A-related genes score (MRGs) was output based on six differentially expressed single-cell m6A-related genes (SCMRDEGs), including GHRL, COL4A1, CAV1, GJA1, TIMP1, and IGFBP3. The protein expression level was assessed by immunohistochemistry. We identified the prognostic value of MRGs and constructed a nomogram model to predict GC patients’ overall survival. MRGs may affect treatment sensitivity in GC patients.

Conclusion

Our study visualized the cellular heterogeneity of TME at the single-cell level, revealed the association between m6A mRNA modification and intracellular communication, clarified MRGs as an independent risk factor of prognosis, and provided a reference for follow-up treatment.

Background

Gastric cancer (GC) is the fifth most common malignancy and the third-leading cause of cancer-related mortality globally [1]. Recently, immune checkpoint inhibitors (ICIs) have significantly improved the long-term prognosis in advanced GC patients, especially those with EBV-positive and microsatellite instable (MSI) [2]. However, a large population received limited benefits or develop drug resistance with a median overall survival (OS) of merely 9–10 months [3]. The in-depth understanding of the molecular mechanisms will contribute to bringing new insights for GC treatment.

There has been an increasing interest in post-translational aberration. Among over 150 chemical modifications, N6-methyladenosine (m6A) RNA methylation, a new level of epigenetic regulation, is the most abundant modification of protein-coding and non-coding RNAs in eukaryotes [4, 5]. Dynamic m6A modification relies on readers, writers, and erasers, respectively responsible for m6A’s function, methylation, and demethylation [6]. Roles of m6A mRNA regulation on tumor cells have been gradually revealed, modulating tumor progression, stemness, and invasive ability [4]. Recent studies provide evidence of m6A mRNA in mediating ICI resistance [7, 8]. It is thereafter that m6A modification-related genes may be potential therapeutic targets to optimize immunotherapy.

Tumorigenesis cannot occur without the favor of tumor microenvironment (TME), comprising not just tumor cells but also cancer-associated fibroblasts (CAFs), infiltrating immune cells, vascular cells, mesenchymal cells [9]. Characterization of the cellular composition of the gastric TME is the fundamental to understand immunosuppression, angiogenesis, and distant metastasis [10, 11]. For example, fibroblasts in the TME have been proven to affect the efficacy of ICIs [12]. The high degree of TME heterogeneity represent important barriers to popularizing ICIs. But the cellular milieu contributing to the malignant nature of the gastric TME remains to be elucidated. To date, how m6A modification-related genes function within the gastric TME was poorly understood.

Here, we combined single-cell and bulk RNA sequencing to stratify the cellular milieu of the gastric TME, explore how m6A-related genes affect intracellular communication in the TME, and screen which m6A-related genes are associated with prognosis and guide follow-up treatment in GC.

Materials and methods

Data acquisition

52 publicly available human gastric cancer single-cell mRNA sequence (scRNA-seq) datasets (GSE150290) were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo) [13, 14], including 23 tumor tissues and 29 normal adjacent tissues. The bulk mRNA sequence datasets were acquired from The Cancer Genome Atlas program-stomach adenocarcinoma (TCGA-STAD) (https://www.cancer.gov/ccg/research/genome-sequencing/tcga) as a training cohort [15], including 343 tumor samples and 31 normal tissues. The GSE66229 was used as a validation cohort. All data analyzed in our study are freely available in public domain or published literatures.

Visualization of cell types in the TME

The Seurat R package 4.0.3 was used for processing scRNA-seq data. The quality standards are as follow: (1) each gene expressed in at least three cells; (2) cells with < 200 detected genes were excluded; (3) cells with ≥ 10% mitochondria-expressed genes were excluded. The top 2000 variable genes were used to normalize RNA counts. Principle component analysis (PCA) was performed to identify significant principle components (PCs). Moreover, we used t-distributed stochastic neighbor embedding (t-SNE) algorithm for dimensionality reduction to further obtain the top PCs. Finally, we annotated and visualized the cell types and subtypes constituting the TME of gastric carcinoma.

Screening single-cell m6A-related genes

We used Seurat’s FindAllMarkers to analyze differential expression to get marker genes in each cell type [16]. We downloaded 701 experimentally verified m6A-related genes from RM2Target [17]. Then, we obtained single-cell m6A-related genes (SCMRGs) by intersecting cell markers and m6A-related genes. The expression of SCMRGs was scored using the AUCell package [18].

Functional enrichment analysis

We performed The Gene Set Variation Analysis (GSVA) to calculate the enrichment score in each cell type. The gene set with a P-value < 0.05 was considered significantly enriched. The clusterProfile package was used to conduct Gene Oncology (GO) [19] and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses [20].

Pseudotime trajectory analysis

Monocle 2 R package revealed pseudo-time trajectories of fibroblasts [21]. We used the DDRTree method for dimensionality reduction, projecting single cells onto space and ordering into a trajectory with branch points. The dynamic expression heatmap was performed with “plot-pseudo-time-heatmap”.

Cell-to-cell communication

According to the ligand-receptor interaction database, CellChat can analyze intracellular communication networks of scRNA-seq data annotated as different cell clusters [22]. The number and strength of ligand-receptor interactions were calculated. The “netVisual-circle function” is used to visualize the intensity and number of interactions between cell types, and the “netVisual-bubble function” is used to visualize receptor-ligand pairs between cell clusters. Ligand-receptor pairs with a P-value < 0.05 were filtered to explore the interaction between different cell types.

Identification of SCMRDEGs

DESeq2 package [23] was used to analyze differentially expressed genes (DEGs) between gastric tumor tissues and para-cancer tissues. Significant DEGs were identified by the ‘limma’ package with | logFC | > 1 and adjusted P-value < 0.05. The intersection of SCMRG and bulk RNA-seq DEGs was used to obtain differentially expressed single-cell m6A-related genes (SCMRDEGs). In the TCGA-STAD cohort, “ConsensusClusterPlus” R package was conducted to cluster GC patients into three clusters based on the expression of SCMRDEGs [24]. The optimal cluster number k = 3 was selected based on cumulative distribution function (CDF). The expression of SCMRDEGs and immune cell infiltration in 3 clusters were compared.

Establishment and validation of the nomogram model

We used univariate Cox and LASSO analyses to screen SCMRDEGs associated with prognosis. An optimal model was constructed with 6 SCMRDEGs, output as m6A-related gene score (MRGs).

$$MRGs=\text{coefficient}+{\sum }_{i}^{1}\beta i\text{*expG}i$$

We used the “survivalROC” R package to verify the predictive efficacy of MRGs in training and validation cohorts, displayed with Kaplan-Meier (K-M) and receiver operating characteristics (ROC) curves [25]. Subsequently, based on univariate and multivariate Cox analyses, we integrated MRGs and classical parameters to establish a nomogram model. Calibration curves were employed to validate the stability of the model.

Immunohistochemistry

A total of 37 pairs of GC and para-cancerous tissue samples were collected from patients who underwent surgical excision at the First Affiliated Hospital of Zhengzhou University and were used for immunohistochemistry (IHC) between September 2016 and April 2017 [26]. Briefly, IHC assays were performed to assess the expression level of selected SCMRDEGs (GHRL, COL4A1, CAV1, GJA1, TIMP1, and IGFBP3). 4 μm thick paraffin-embedded tissues sections were stained with the GHRL (YT1900, 1:200, Immunoway, USA), COL4A1 (CY1657, 1:150, abways, China), CAV1 (CY5021, 1:150, abways, China), GJA1 (26980-1-AP, 1:200, proteintech, China), TIMP1 (CY6712, 1:150, abways, China), and IGFBP3 (10189-2-AP, 1:200, proteintech, China) antibody per manufacturer’s instructions. The tissue sample collection was approved by the Ethics Committee of Scientific Research of the First Affiliated Hospital of Zhengzhou University. All patients were informed and have written the informed consents.

Immune infiltration analysis

The MRG-score was divided into high MRGs type and low MRGs type according to the median MRG-score in the training cohort. To assess the relationship between MRGs and immunological state in the TME, we estimated the “StromalScore, ImmuneScore, EstimateScore, and TumorPurity” using the “ESTIMATE” functions of “IOB” R package. The proportional infiltration levels of 22 types of immune cells in two MRGs types were quantified and compared using “CIBERSORT” [27].

Chemotherapy sensitivity and immune checkpoint analysis

We used the “oncoPredict” package to analyze different drug sensitivity between high- and low- MRGs groups [28]. Six chemotherapeutic drugs with the most apparent difference in half-maximal concentration (IC50) were selected to display. We also analyzed the correlation between MRG score and immune checkpoints, aiming to select the potential immunotherapy-sensitive populations.

Statistical analyses

Statistical analyses were conducted using R software (version 4.1.2). K-M curves and log-rank tests were applied for survival analyses. A two-sided P-value < 0.05 was considered statistically significant.

Results

Visualization of TME cell types in GC

The overall flow diagram was displayed in Fig. 1. To characterize the TME in gastric cancer, we performed single-cell RNA-seq analysis on 23 tumor samples from the GEO dataset. All cells were divided into 21 clusters (0–20) (Fig. 2A). All clusters could be further divided into 9 cell subtypes, including epithelial cells, endotheliocytes, fibroblasts, myeloid cells, T/NK cells, B cells, plasma cells, pericytes, and mast cells (Fig. 2B). The proportion of each cell lineage in 23 tumor samples was visualized in Fig. 2C. We observed that there are significant differences in the proportion of cell types in different tumor samples, reflecting the high heterogeneity of gastric cancer at the cellular level. The underlying mechanism of heterogeneity might contribute to optimizing individualized treatment. We also plotted different cell types and their corresponding high-expression genes in the heatmap (Fig. 2D). Canonical cell markers identified by previous studies are displayed in Fig. 2E. The cell subtypes we have identified corresponded to their signature cell markers (Fig. 2F).

Fig. 1
figure 1

Schematic diagram of the overall study design. SCMRG, single cell m6A-related genes; DEGs, differentially expressed genes; SCMRDEGs, single cell m6A-related differentially expressed genes; MRGs, m6A-related gene score; TME, tumor microenvironment

Full size image
Fig. 2
figure 2

The association of cellular heterogeneity and m6A-related genes in gastric TME. (A) t-SNE plot of 21 cell clusters. (B) t-SNE plot of 9 cell types. (C) Distribution of 9 cell types in each included sample. (D) Heatmap of highly expressed genes in each cell type. (E) Marker genes of each cell type verified by published studies. (F) Expression levels of marker genes in the 9 cell types. (G) Differentially expressed cell markers in GSE150290. (H) Identification of single-cell m6A-related genes (SCMRGs); (I) t-SNE plot of expression level of SCMRDEGs in each cell type. (J) Ridge map of expression level of SCMRDEGs in each cell type. K. Gene ontology (GO) analysis of 114 SCMRGs. L. Kyoto Encylopedia of Genes and Genomes (KEGG) analysis of 114 SCMRGs.

Full size image

Identification of single-cell m6A-related genes (SCMRG)

We used the Seurat package to list differentially expressed genes (DEGs) of 9 cell types relative to others, indicating that each cell type has its specific expression profile (Fig. 2G). We identified 114 single-cell m6A-related genes (SCMRGs) by intersecting 701 experimentally verified m6A-related genes with DEGs (Fig. 2H).

The overall expression level of 114 SCMRGs in each cell type was displayed in a UMAP map (Fig. 2I). The higher SCMRG expression was mainly in the interstitial cells. As shown in the ridge map (Fig. 2J), the expression of SCMRGs was highest in fibrocytes, followed by myeloid cells and pericytes. M6A-related genes might play a vital role in the interstitial cells.

To analyze the biological function and associated pathways, we performed GO (Fig. 2K) and KEGG (Fig. 2L) analyses. Regarding biological process (BP), 114 SCMRGs are related to extracellular matrix structural constituent, integrin binding, and glycosaminoglycan binding. For cellular components (CC), these genes are mostly enriched in collagen-containing extracellular matrix and endoplasmic reticulum lumen. For molecular function (MF), these SCMRGs are associated with external encapsulating structure organization, extracellular structure organization, and extracellular matrix organization. KEGG analysis showed that 114 SCMRGs are related to focal adhesion, intestinal immune network for IgA production, and PI3K-AKT signaling pathway.

M6A-mediated fibroblasts guide intracellular communication of the TME

Considering the highest SCMRG expression in fibroblasts, we extracted them for dimensionality reduction clustering. Fibroblasts were divided into 11 sub-clusters (Fig. 3A) and annotated into 4 types: iCAFs, myCAFs, apCAFs, and lipo-fibroblasts (Fig. 3B). There was no difference in SCMRG expression among the four fibrocyte subsets (Fig. 3C). The cell markers of fibroblasts identified by published literature were visualized as a bubble plot (Fig. 3D). The Fibrocyte subsets we identified corresponded to their cell markers (Fig. 3E).

Fig. 3
figure 3

M6A-related mRNA modification modifies the features of fibroblasts. (A) t-SNE plot of 11 fibroblast subclusters in the TME. (B) t-SNE plot of 4 fibroblast subtypes in the TME. (C) t-SNE plot of MRG score in 4 fibroblast subtypes. (D) Marker genes of each fibroblast subtype. (E) Expression level of marker genes in 4 fibroblast subtypes. F-H. Differentiation trajectory of fibroblasts, colored for cell types (F), states (G), and pseudotime (H). I. Heatmap of dynamic changes of 114 SCMRGs with pseudotime. J. The number of interactions among different TME cell subtypes. K. The strength of interactions among different TME cell subtypes. L. Receptor-ligand pairs that regulate fibroblasts. M. Receptor-ligand pairs regulated by fibroblasts

Full size image

We performed a pseudo-time analysis (Fig. 3F-H) to clarify the differentiation transition of fibrocyte subtypes. These cells formed a continuum with 9 cell states. During pseudo-time, iCAFs and apCAFs gradually decreased while lipo-fibroblast stepwise increased. It is noted that mCAFs did not change significantly. The heatmap (Fig. 3I) showed the 114 SCMRG expression dynamic changes along the pseudo-time. M6A-related genes might mediate the differentiation transition among different fibroblast subtypes.

Cell-chat analysis was undertaken to explore diverse interactions between fibroblasts and other cells in the TME. The number and strength of intracellular communication were presented in Fig. 3J and K, respectively. Fibroblasts had a relatively higher intensity and frequency interplay with other cell types. There are the most interactions between fibrocytes and pericytes and the most potent interactions between fibrocytes and myeloid cells. By analyzing receptor-ligand pairs (Fig. 3L-M), fibroblasts mainly serve as signal senders than receivers. Specifically, fibrocytes mainly received regulatory signals from peripheral cells, endothelial cells, and tumor cells. In contrast, fibrocytes could send regulatory signals to regulate all other cell types, consistent with the results in Fig. 3J-K. Therefore, m6A-mediated fibroblasts may guide intracellular communications among TME cells, supporting its potential as a novel target.

Other cell types in the TME

Myeloid cells are crucial components of the TME, exerting both tumor-stimulating and suppressing roles. Subgroup analysis successfully divided 13 distinct clusters (Figure S1A) and revealed 7 myeloid cell subtypes: macro-THBS1, monocytes, macro-APOE, cDC1-CD1C, DC-LAMP3, pDC-LILRA4, cDC1-XCR1 (Figure S1B). We collected specific markers of myeloid cells in Figure S1C. Also, we used a violin plot to validate the expression of typical cell markers in 7 cell subsets (Figure S1D). These markers could well distinguish each myeloid cell subset. We observed that there was no significant difference among the scores of 114 SCMRGs in different myeloid subgroups (Figure S1E).

Although MRG scores of immune cells were not high, T/NK cells were extracted to investigate whether there were differences in MRG scores among different subgroups. Similarly, we further classified 12 cell clusters (Figure S2A) into 5 cell types, including CD4+ T cell, CD8+ T cell, cycling T cell, regulatory T cell, and NK cell (Figure S2B). T/NK cell markers were collected through single-cell-related studies (Figure S2C). The violin plot was performed to verify the expression level of cell markers in each cell type (Figure S2D). It is noted that there was no significant difference in 114 SCMRG scores among 5 types of T/NK cells from the UMAP plot (Figure S2E).

Identification of SCMRDEGs

We used the TCGA-STAD transcriptome data to assess differentially expressed genes (DEGs) with the cut-off value of |log2-fold change| ≥2 and adjusted P-value < 0.05. A total of 4661 DEGs (2371 upregulated and 2290 downregulated) were identified between normal and tumor samples (Fig. 4A). By intersecting 4661 DEGs and 114 SCMRGs, we finally obtained 32 single-cell m6A-related differentially expressed genes (SCMRDEGs), with 16 upregulated and 16 downregulated genes (Fig. 4B). A heatmap (Fig. 4C) showed expression differences of 32 SCMRDEGs between tumor and normal tissues. GHRL, KLF4, and DUSP5 expression levels were relatively lower in tumor tissues, whereas DSP, PERP, and MMP9 were highly expressed in tumor tissues. Furthermore, GO analysis indicated that 32 SCMRDEGs were enriched in metalloendopeptidase activity, serine hydrolase activity, cell-substrate junction, focal adhesion, external encapsulating structure organization, and extracellular structure organization (Fig. 4D).

Fig. 4
figure 4

Identification of single cell m6A-related differentially expressed genes (SCMRDEGs). (A) Volcanic map for differentially expressed genes (DEGs) between normal and tumor tissues. (B) Venn diagram indicates overlapped SCMRDEGs between DEGs and single-cell m6A-related genes. (C) Heatmap of expression levels of DEGs in TCGA-STAD dataset. (D) GO analysis of 32 SCMRDEGs. (E) Delta area under the cumulative distribution function (CDF) curve. (F) Heatmap of clustering at consensus k = 3. (G) CDF curves of different consensus k-value. (H) Principle component analysis (PCA) of three clusters. (I) Immune infiltration of 22 immune cells in the three clusters. (J) Differences of immune checkpoints among three clusters

Full size image

Consensus clustering analysis

To clarify the association between SCMRDEGs and GC subtypes, we performed consensus clustering analysis with 343 TCGA-STAD samples. When K > 3, we observed that the area under the cumulative distribution function (CDF) curve did not increase significantly (Fig. 4E). The optimal division was achieved when k = 3 based on the clustering heatmap (Fig. 4F) and consensus CDF (Fig. 4G). Therefore, k = 3 was selected to divide TCGA-STAD samples into three subtypes: MC1, MC2 and MC3. Principle component analysis (PCA) results demonstrated that these 3 samples could not be completely separated by SCMRDEGs but with a clear degree of distinguish (Fig. 4H). The CIBERSORT algorithm compared the immune infiltration score among 3 clusters. As shown in Fig. 4I, Naïve B cells, resting mast cells, resting memory CD4+ T cells, and regulatory T cells were enriched in the cluster-MC1; M0-macrophages, activated mast cells, and neutrophils were enriched in the cluster-MC2; activated dendritic cells and activated memory CD4+ T cells were enriched in the cluster-MC3. Except for JUNB, SOX4, C12orf75, and KLF4, the expression levels of immune checkpoints were significantly different in all three clusters (Fig. 4J). Notably, the expression level of all immune checkpoints in cluster MC2 were relatively high. It is estimated that patients in cluster MC2 are more likely to respond to immunotherapy.

Establishment and validation of a nomogram model

Univariate Cox regression analysis screened 10 genes that were associated with the prognosis of GC patients, including COL1A1, TIMP1, COL3A1, ACTA2, TIMP3, IGFBP3, COL4A1, CAV1, GJA1, and GHRL. The selected genes were subsequently included in the Lasso analysis to select the optimal model, which output m6A-related gene score (MRGs) based on 6 SCMRDEGs (GHRL, COL4A1, CAV1, GJA1, TIMP1, and IGFBP3) (Fig. 5A-B). The formula is shown as follows:

Fig. 5
figure 5

Construction and validation of a MRGs-based nomogram model. A-B. The selection of prognostic SCMRDEGs based on the optimal parameter 𝜆 that was obtained in LASSO analysis. C-D. K-M survival curves of high-MRGs and low-MRGs groups in training and validation cohorts. E-F. ROC curves of MRGs in training and validation cohorts. G. Multivariate cox regression analysis in the training cohort. H. A nomogram model was constructed with MRGs, age, and AJCC staging. I-K. Calibration curves were used to validate the predictive accuracy of our nomogram for predicting 1-, 3-, and 5-year OS. L. IHC assessed the protein expression level of GHRL, COL4A1, CAV1, GJA1, TIMP1, and IGFBP3.

Full size image
$$MRGs=\text{coefficient}+{\sum }_{i}^{1}\beta i\text{*expG}i$$

Patients were divided into high and low MRGs groups based on the median value. Kaplan-Meier curves were plotted with THE TCGA-STAD training and GEO validation cohorts. Results showed there were significant prognostic differences (P < 0.05) between high- and low-MRGs groups (Fig. 5C-D). The AUCs for predicting 1-, 3-, and 5-year OS in the training cohort were 0.619, 0.632, 0.754 (Fig. 5E), and those in the validation cohort were 0.558, 0.587, 0.583 (Fig. 5F). This finding indicated that MRGs had a favorable predictive ability in the prognosis of GC patients.

Furthermore, the univariate and multivariate Cox regression analyses identified that MRGs and age were independent factors affecting the OS of GC patients (Fig. 5G). AJCC staging also had a significant impact on prognosis. Therefore, a prognostic nomogram model was established based on MRGs, age, and AJCC staging (Fig. 5H). The calibration curves (Fig. 5I-K) of the 1-, 3- 5-year OS demonstrated that the observed results were consistent with the predicted values, verifying the stability of our nomogram model.

The IHC assays (Fig. 5L) verified the differential expressions of six SCMRDEGs in the tumor and adjacent tissues. Of the 24 tissue sections, proteins encoded by COL4A1, IGFBP3, and TIMP1 were slightly up-regulated in tumor tissues, while GHRL and GJA1 proteins were relatively down-regulated in tumor tissues than para-tumor tissues. No significant difference was observed in the expression of CAV1 protein between tumor and normal tissues, partly because it is mainly expressed in the stroma.

Association of MRGs and immune infiltration

To infer the relationship between MRGs and immune infiltration, we used the “estimate” R package to calculate Stromal Score (Fig. 6A), Immune Score (Fig. 6B), ESTIMATE Score (Fig. 6C), and Tumor Purity (Fig. 6D). The high-MRGs group had a relatively higher immune score and stromal score compared to low-MRGs group, but the corresponding tumor purity was lower. Moreover, naïve B cells, resting dendritic cells, resting mast cells, monocytes, and resting memory CD4+ T cells were enriched in the high-MRGs group, whereas M0-macrophages, activated memory CD4+ T cells, follicular helper T cells, and γδ T cells had a higher infiltration in the low-MRGs group (Fig. 6E). Also, we found correlations between 6 SCMRDEGs and immune infiltrating cells (Fig. 6F).

Fig. 6
figure 6

Immune infiltration and drug sensitivity analyses. A. Stromal score. B. Immune score. C. ESTIMATE score. D. Tumor purity. (E) Immune infiltration of 22 immune cells between high- and low-MRGs groups. (F) The correlation between six SCMRDEGs and immune cell infiltration. G-L. Violin diagrams of the top six drugs with the most significant difference in drug sensitivity between high-MRGs and low-MRGs groups. M. A heatmap displayed the correlations between 47 immune checkpoint molecules and 6 SCMRDEGs.

Full size image

Association of MRGs and drug sensitivity

We used the oncoPredict package to analyze drug sensitivity and obtain the half-maximal inhibitory concentration (IC50) of all chemotherapeutic drugs corresponding to each sample. The violin diagrams (Fig. 6G-L) showed the top 6 drugs with the most significant difference in drug sensitivity between the high-MRGs and the low-MRGs groups. Noteworthily, the high-MRGs group was more sensitive to these 6 chemotherapeutic agents than the low-MRGs group. The heatmap (Fig. 6M) displayed the correlations between 47 immune checkpoints and 6 SCMRDEGs that comprise the MRG score. As the darker color indicates a stronger correlation, TIMP1 was significantly negatively correlated with TNFRSF14 (P < 0.05), and GHRL was significantly negatively correlated with TNFSF9 (P < 0.05). A significantly positive correlation was observed between IGFBP3 and CD276 (P < 0.05). COL4A1, CAV1, and GJA1 were positively correlated with PDCD1LG2 (P < 0.05), while both CAV1 and GJA1 were negatively correlated with LGALS9 (P < 0.05).

Discussion

There have been extensive studies on the correlation between m6A mRNA modification and the pathogenesis of colorectal cancer [29,30,31,32,33]. Still, fewer studies provided an m6A mRNA modification landscape of the GC microenvironment at the single-cell level. In the present study, we leverage the advantage of combined single-cell and bulk RNA sequencing to clarify the heterogeneity of TME in GC, screened significant m6A-related genes as independent risk factors, and predicted treatment response in patients with distinct MRGs. This unique perspective may provide a deeper understanding of how m6A mRNA modification functions within the TME to impact the prognosis of GC patients.

Recent studies have revealed that the high cellular heterogeneity of the TME supports uncontrolled growth and facilitates immune evasion of solid tumors [34, 35]. In our study, nine cell types were identified preliminarily: epithelial cells, endotheliocytes, fibroblasts, myeloid cells, T/NK cells, B cells, plasma cells, pericytes, and mast cells. We found that cell types with high MRGs are mainly stromal cells. Fibroblasts ranked the highest MRG score, followed by pericytes and myeloid cells. This finding is consistent with subpopulations harboring the highest intensity and frequency of intracellular communications in the TME. It is inferred that m6A-related genes may guide extensive and enhanced communications between stromal cells in gastric TME.

Cancer-associated fibroblasts (CAF), a crucial stromal cell component, were further classified as iCAFs, myCAFs, apCAFs, and lipo-fibroblasts, based on specific molecular characteristics [36]. By analyzing receptor-ligand pairs, we found that fibroblasts are more likely to be a sender than a receiver in the cell-chat analysis, sending upstream signals to all other cell types. Also, a recent study manifested that m6A-mediated fibroblasts communicate with epithelial cells more extensively than non-m6A-mediated fibroblasts in GC [34]. Growing evidence suggested that CAFs may drive an immunosuppressive TME via secreting CXCL2, IL6, and CCL2 [37,38,39,40]. Activated fibroblasts are associated with non-response to immunotherapy in pancreatic and breast cancers [41, 42]. It is reasonable to speculate that m6A-mediated fibroblasts may form immunosuppressive interplays with other TME cells. Furthermore, m6A-related modification might mediate the trajectory process of fibroblast subtypes, regulating various biological processes in GC tumorigenesis and progression. Hence, targeting m6A-related fibroblasts is a promising strategy to modulate gastric TME and overcome drug resistance.

We obtained 32 single-cell m6A-related differentially expressed genes by integrating bulk RNA sequencing data. Consensus clustering analysis was performed to divide GC samples into three clusters (MC1, MC2, and MC3) with distinct m6A expression patterns. We observed that the expression of all immune checkpoints in cluster-MC2 was relatively higher, indicating better immunotherapy responses than other clusters. To further evaluate the predictive value of these SCMRDEGs, we combined univariate and LASSO regression analyses to determine six candidate SCMRDEGs (GHRL, COL4A1, CAV1, GJA1, TIMP1, and IGFBP3) for calculating MRGs. We used the IHC to verify the protein expression level encoded by six SCMRDEGs. Moreover, patients were divided into high- and low-MRGs groups according to the median value of MRGs. By K-M survival and ROC curve analyses, MRGs showed a good predictive ability in the prognosis of GC patients. Patients with high-MRGs have a shorter OS than those in the low-MRGs group. Ultimately, we established and validated a nomogram model based on MRGs, age, and AJCC staging to predict GC patients’ 1-, 3-, and 5-year OS.

Additionally, we assessed the relationship between MRGs and immunological state. Compared with the low-MRGs group, the high-MRGs group tends to have a higher immune and stromal score but lower tumor purity. Multiple resting immune cells showed enrichment in the high-MRGs group, whereas activated immune cells displayed a relatively higher infiltration in the low-MRGs group. Patients in the low-MRGs group may drive a quicker anti-tumor immune response and receive more benefits from immunotherapy. We also comprehensively evaluated the correlations between MRGs and drug sensitivity. Patients in the high-MRGs group were more sensitive to chemotherapy than those in the low-MRGs group. These chemotherapeutic agents include BMS.754807-2172, NU7441-1038, AZD8186-1918, JQ1-2172, AZD8055-1059, and Dasatinib-1079. Therefore, we could use MRGs to predict different responses to immunotherapy and chemotherapy among distinct GC patients.

Some limitations in our study should be highlighted. Firstly, detailed mechanisms of m6A mRNA modifications in the multiple TME cells need experimental validation. Secondly, single-cell analysis needs to be more in-depth, and more clinical samples are required. We will verify our findings with a larger sample size and longer follow-up periods in future studies. Thirdly, given that these analyses were based on published databases, it is necessary to validate our findings in real-world cohorts, ensuring the robustness of MRGs as a predictive marker for GC prognosis. Nonetheless, the single-cell and bulk RNA sequencing analyses provide a novel perspective on how m6A mRNA modifications function within the heterogenous TME.

Conclusion

This study integrated scRNA-seq and bulk RNA-seq data to identify m6A-modified cellular heterogeneity of TME, reveal m6A-mediated fibroblasts guiding intercellular communication of gastric TME, determine the prognostic value of MRGs, and evaluate the effects of MRGs on treatment sensitivity. This study provides a new perspective for an in-depth understanding of the characteristics of m6A mRNA modification in the TME cell subtype, which is a critical step for clinical practice and individualized therapy.

Data availability

The datasets of this article were generated from the TCGA database and GEO database.

References

  1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer statistics 2020: GLOBOCAN estimates of incidence and Mortality Worldwide for 36 cancers in 185 countries. Cancer J Clin. 2021;71(3):209–49.

    Article  Google Scholar 

  2. Kim ST, Cristescu R, Bass AJ, Kim KM, Odegaard JI, Kim K, et al. Comprehensive molecular characterization of clinical responses to PD-1 inhibition in metastatic gastric cancer. Nat Med. 2018;24(9):1449–58.

    Article  CAS  PubMed  Google Scholar 

  3. Olino K, Park T, Ahuja N. Exposing hidden targets: combining epigenetic and immunotherapy to overcome cancer resistance. Sem Cancer Biol. 2020;65:114–22.

    Article  CAS  Google Scholar 

  4. Li M, Zha X, Wang S. The role of N6-methyladenosine mRNA in the tumor microenvironment. Biochim et Biophys acta Reviews cancer. 2021;1875(2):188522.

    Article  CAS  Google Scholar 

  5. Roignant JY, Soller M. M(6)A in mRNA: an ancient mechanism for fine-tuning gene expression. Trends Genet. 2017;33(6):380–90.

    Article  CAS  PubMed  Google Scholar 

  6. Jiang X, Liu B, Nie Z, Duan L, Xiong Q, Jin Z, et al. The role of m6A modification in the biological functions and diseases. Signal Transduct Target Therapy. 2021;6(1):74.

    Article  CAS  Google Scholar 

  7. Wang L, Hui H, Agrawal K, Kang Y, Li N, Tang R, et al. M(6) a RNA methyltransferases METTL3/14 regulate immune responses to anti-PD-1 therapy. EMBO J. 2020;39(20):e104514.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Yin H, Zhang X, Yang P, Zhang X, Peng Y, Li D, et al. RNA m6A methylation orchestrates cancer growth and metastasis via macrophage reprogramming. Nat Commun. 2021;12(1):1394.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Papalexi E, Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity. Nat Rev Immunol. 2018;18(1):35–45.

    Article  CAS  PubMed  Google Scholar 

  10. Gonzalez H, Hagerling C, Werb Z. Roles of the immune system in cancer: from tumor initiation to metastatic progression. Genes Dev. 2018;32(19–20):1267–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Sahai E, Astsaturov I, Cukierman E, DeNardo DG, Egeblad M, Evans RM, et al. A framework for advancing our understanding of cancer-associated fibroblasts. Nat Rev Cancer. 2020;20(3):174–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Bruni D, Angell HK, Galon J. The immune contexture and immunoscore in cancer prognosis and therapeutic efficacy. Nat Rev Cancer. 2020;20(11):662–80.

    Article  CAS  PubMed  Google Scholar 

  13. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets–update. Nucleic Acids Res. 2013;41(Database issue):D991–5.

    CAS  PubMed  Google Scholar 

  14. Kim J, Park C, Kim KH, Kim EH, Kim H, Woo JK, et al. Single-cell analysis of gastric pre-cancerous and cancer lesions reveals cell lineage diversity and intratumoral heterogeneity. NPJ Precision Oncology. 2022;6(1):9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Mounir M, Lucchetta M, Silva TC, Olsen C, Bontempi G, Chen X, et al. New functionalities in the TCGAbiolinks package for the study and integration of cancer data from GDC and GTEx. PLoS Comput Biol. 2019;15(3):e1006701.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Zhang X, Lan Y, Xu J, Quan F, Zhao E, Deng C, et al. CellMarker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res. 2019;47(D1):D721–d8.

    Article  CAS  PubMed  Google Scholar 

  17. Bao X, Zhang Y, Li H, Teng Y, Ma L, Chen Z, et al. RM2Target: a comprehensive database for targets of writers, erasers and readers of RNA modifications. Nucleic Acids Res. 2023;51(D1):D269–d79.

    Article  CAS  PubMed  Google Scholar 

  18. Aibar S, González-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods. 2017;14(11):1083–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Gene Ontology Consortium. Going forward. Nucleic Acids Res. 2015;43(Database issue):D1049–56.

    Article  Google Scholar 

  20. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. 2017;14(10):979–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. 2021;12(1):1088.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Love MI, Huber W, Anders S. Moderated estimation of Fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinf (Oxford England). 2010;26(12):1572–3.

    CAS  Google Scholar 

  25. Heagerty PJ, Lumley T, Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2000;56(2):337–44.

    Article  CAS  PubMed  Google Scholar 

  26. Wang Z, Chen C, Ai J, Shu J, Ding Y, Wang W, et al. Identifying mitophagy-related genes as prognostic biomarkers and therapeutic targets of gastric carcinoma by integrated analysis of single-cell and bulk-RNA sequencing data. Comput Biol Med. 2023;163:107227.

    Article  CAS  PubMed  Google Scholar 

  27. Zeng D, Ye Z, Shen R, Yu G, Wu J, Xiong Y, et al. IOBR: Multi-omics Immuno-Oncology Biological Research to Decode Tumor Microenvironment and signatures. Front Immunol. 2021;12:687975.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. 2021;22(6).

  29. Chen RX, Chen X, Xia LP, Zhang JX, Pan ZZ, Ma XD, et al. N(6)-methyladenosine modification of circNSUN2 facilitates cytoplasmic export and stabilizes HMGA2 to promote colorectal liver metastasis. Nat Commun. 2019;10(1):4695.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Li T, Hu PS, Zuo Z, Lin JF, Li X, Wu QN, et al. METTL3 facilitates tumor progression via an m(6)A-IGF2BP2-dependent mechanism in colorectal carcinoma. Mol Cancer. 2019;18(1):112.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Wu Y, Yang X, Chen Z, Tian L, Jiang G, Chen F, et al. M(6)A-induced lncRNA RP11 triggers the dissemination of colorectal cancer cells via upregulation of Zeb1. Mol Cancer. 2019;18(1):87.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Chen H, Gao S, Liu W, Wong CC, Wu J, Wu J, et al. RNA N(6)-Methyladenosine methyltransferase METTL3 facilitates colorectal Cancer by activating the m(6)A-GLUT1-mTORC1 Axis and is a therapeutic target. Gastroenterology. 2021;160(4):1284–300e16.

    Article  CAS  PubMed  Google Scholar 

  33. Sun L, Wan A, Zhou Z, Chen D, Liang H, Liu C, et al. RNA-binding protein RALY reprogrammes mitochondrial metabolism via mediating miRNA processing in colorectal cancer. Gut. 2021;70(9):1698–712.

    Article  CAS  PubMed  Google Scholar 

  34. Gao Y, Wang H, Chen S, An R, Chu Y, Li G, et al. Single-cell N(6)-methyladenosine regulator patterns guide intercellular communication of tumor microenvironment that contribute to colorectal cancer progression and immunotherapy. J Translational Med. 2022;20(1):197.

    Article  CAS  Google Scholar 

  35. Xie S, Cai Y, Chen D, Xiang Y, Cai W, Mao J, et al. Single-cell transcriptome analysis reveals heterogeneity and convergence of the tumor microenvironment in colorectal cancer. Front Immunol. 2022;13:1003419.

    Article  CAS  PubMed  Google Scholar 

  36. Galbo PM Jr., Zang X, Zheng D. Molecular features of Cancer-associated fibroblast subtypes and their implication on Cancer Pathogenesis, Prognosis, and Immunotherapy Resistance. Clin cancer Research: Official J Am Association Cancer Res. 2021;27(9):2636–47.

    Article  CAS  Google Scholar 

  37. Kalluri R. The biology and function of fibroblasts in cancer. Nat Rev Cancer. 2016;16(9):582–98.

    Article  CAS  PubMed  Google Scholar 

  38. Miyake M, Hori S, Morizawa Y, Tatsumi Y, Nakai Y, Anai S, et al. CXCL1-Mediated Interaction of Cancer cells with Tumor-Associated macrophages and Cancer-Associated fibroblasts promotes Tumor Progression in human bladder Cancer. New York, NY: Neoplasia; 2016;18(10):636–46.

    Google Scholar 

  39. Kobayashi H, Enomoto A, Woods SL, Burt AD, Takahashi M, Worthley DL. Cancer-associated fibroblasts in gastrointestinal cancer. Nat Reviews Gastroenterol Hepatol. 2019;16(5):282–95.

    Article  Google Scholar 

  40. Zhao Q, Huang L, Qin G, Qiao Y, Ren F, Shen C, et al. Cancer-associated fibroblasts induce monocytic myeloid-derived suppressor cell generation via IL-6/exosomal miR-21-activated STAT3 signaling to promote cisplatin resistance in esophageal squamous cell carcinoma. Cancer Lett. 2021;518:35–48.

    Article  CAS  PubMed  Google Scholar 

  41. Dominguez CX, Müller S, Keerthivasan S, Koeppen H, Hung J, Gierke S, et al. Single-cell RNA sequencing reveals stromal evolution into LRRC15(+) myofibroblasts as a determinant of patient response to Cancer Immunotherapy. Cancer Discov. 2020;10(2):232–53.

    Article  CAS  PubMed  Google Scholar 

  42. Kieffer Y, Hocine HR, Gentric G, Pelon F, Bernard C, Bourachot B, et al. Single-cell analysis reveals fibroblast clusters linked to Immunotherapy Resistance in Cancer. Cancer Discov. 2020;10(9):1330–51.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by National Natural Science Foundation of China (No. 82273381) and National Natural Science Foundation of China (No. 82203852).

Author information

Author notes

  1. Zehua Wang and Chen Chen contributed equally to this work.

Authors and Affiliations

  1. Department of Oncology, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, 450052, China

    Zehua Wang, Chen Chen, Yihan Liu, Haoyue Cao, Yongxu Jia & Yanru Qin

  2. Department of Pathology, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, 450052, China

    Jiao Shu

  3. The First Affiliated Hospital of Nanchang University, Nanchang, 330006, China

    Jiaoyu Ai

Authors
  1. Zehua Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Chen Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Jiao Shu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Jiaoyu Ai
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Yihan Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Haoyue Cao
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Yongxu Jia
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Yanru Qin
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

YQ, YJ, ZW, and CC conceived and designed the study. ZW, CC, and YL acquired, analyzed, and interpreted of data. ZW and CC drafted the manuscript. YL, JS, HC and JA critically revised of the manuscript for important intellectual content. YQ obtained funding. YQ and YJ provided administrative, technical, and material support. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Yongxu Jia or Yanru Qin.

Ethics declarations

Ethics approval and consent to participate

Written informed consent for the use of tissue samples was obtained from all patients or their legal guardians. The study was approved by the Ethics Committee of Scientific Research of the First Affiliated Hospital of Zhengzhou University.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

: Fig. S1. Myeloid cells in the TME. (A) t-SNE plot of 13 myeloid cell subclusters in the TME. (B) t-SNE plot of 7 myeloid cell subtypes in the TME. (C) Marker genes of each myeloid cell type. (D) Expression level of marker genes in 7 myeloid cell subtypes. (E) t-SNE plot of MRG score in 7 myeloid cell subtypes

Supplementary Material 2

: Fig. S2. T/NK cells in the TME. (A) t-SNE plot of 12 T/NK cell subclusters in the TME. (B) t-SNE plot of 5 T/NK cell subtypes in the TME. (C) Marker genes of each T/NK cell type. (D) Expression level of marker genes in 5 T/NK cell subtypes. (E) t-SNE plot of MRG score in 5 T/NK cell subtype

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, Z., Chen, C., Shu, J. et al. Single-cell N6-methyladenosine-related genes function within the tumor microenvironment to affect the prognosis and treatment sensitivity in patients with gastric cancer. Cancer Cell Int 24, 44 (2024). https://doi.org/10.1186/s12935-024-03227-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12935-024-03227-2

Keywords

Cancer Cell International

ISSN: 1475-2867

Contact us