Cross-talk of m6A methylation modification and the tumor microenvironment composition in esophageal cancer

Background: Esophageal cancer (EC) remains a significant clinical challenge, characterized by its aggressive nature and poor prognosis. Current therapeutic strategies, including targeted therapies, have limitations due to the complex interplay between tumor heterogeneity and the tumor microenvironment (TME). However, the specific contributions of N6-methyladenosine (m⁶A) methylation to the TME in EC are yet to be fully elucidated.

Methods: Through comprehensive bioinformatics analyses, a detailed examination of m⁶A regulators were conducted in EC using datasets from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO). Single-cell RNA sequencing (scRNA-seq) and a consensus clustering algorithm was employed to classify m⁶A modification patterns and analyze their relationships with immune cell infiltration and clinical outcomes. Additionally, an m⁶A scoring system was developed based on principal component analysis to assess the prognostic value of identified m⁶A modification patterns.

Results: The findings revealed two distinct m⁶A modification clusters associated with divergent TME characteristics and immune infiltration profiles. Patients exhibiting the immune-inflamed phenotype (m⁶A cluster B) demonstrated significantly improved survival compared to those with the immune-excluded phenotype (m⁶A cluster A). Notably, m⁶A scores correlated positively with immune cell presence and related with adverse prognostic outcomes, indicating their potential as predictive biomarkers for immunotherapy responses. A low m⁶A score indicated a better response to immunotherapy.

Conclusion: This study highlights the critical role of m⁶A methylation in shaping the TME and influencing immune dynamics in EC. The m⁶A score developed herein provides a novel quantitative tool for predicting tumor behavior and treatment efficacy, paving the way for more personalized immunotherapeutic strategies in clinical practice. This scoring system illustrates a strong correlation of EC with TME immune cell composition, suggesting potential as a biomarker for targeted therapeutic strategies for EC.

1 Introduction

N6-methyladenosine (m⁶A) has been recognized as a crucial RNA modification. m⁶A modifications are dynamically regulated by various regulators, including methyltransferase complex writers such as METTL3, METTL14, METTL16, RBM15, RBM15B, VIRMA, WTAP and ZC3H13, and several binding proteins, such as FMR1, HNRNPC, HNRNPA2B1, IGFBP1/2/3, LRPPRC, RBMX, YTHDC1/2, and YTHDF1/2/3, have been identified as readers, as well as demethylases erasers such as FTO and ALKBH5. Numerous studies have demonstrated that aberrant expression of m⁶A core modification and reading proteins is implicated in diverse physiological and pathological processes, including biological growth and development, DNA damage repair, biological rhythms, angiogenesis, and various types of tumors (1). In recent years, substantial progress has been achieved in m⁶A epitranscriptomics, underscoring its pivotal roles in cancer initiation and progression by modulating RNA stability, mRNA splicing, microRNA processing, and mRNA translation (2). Unlike genetic events, m⁶A modifications are reversible, making epigenetic regulation particularly interesting for the development of new therapeutic technologies for cancer treatment.

Esophageal cancer has the sixth highest cancer-related mortality rate, but research data on this disease are limited compared to other cancers (3–5). Esophageal cancer is characterized by its aggressive nature and dismal 5-year survival rate, which stands at only 18% (6). Recent advances in the identification of molecular markers specific to esophageal cancer have led to the development of novel targeted therapy approaches by targeting these markers (7–12). However, inhibitors have the potential to cause primary or acquired resistance in patients who receive these treatments (13–16). Furthermore, in esophageal cancer, a diverse array soluble immunosuppressive factors and cells with immunosuppressive properties can interfere with immune effector cells, thereby creating a distinct immunosuppressive microenvironment.

Multiple factors influence the outcome of multi-modality treatments. An individual tumor’s intrinsic features are crucial to shaping its immune microenvironment and affecting the effectiveness of immunotherapy in esophageal cancer (17). As our understanding of the tumor microenvironment deepens, we increasingly recognize the significance of immune cell subsets in tumor development and the identification of potential therapeutic targets. The microenvironment in esophageal cancer is complex, comprising of various components such as NK cells, tumor-associated macrophages, dendritic cells (DCs), myeloid-derived suppressor cells (MDSCs), neutrophils, mast cells (MCs), eosinophils, endothelial cells, tumor angiogenesis, and cancer-associated fibroblasts (CAFs) (18, 19). Extensive exploration has been conducted on the utilization of clinical immunotherapy approaches that target innate immune cells as adjuvant therapies in conjunction with surgical resection and chemoradiotherapy for the treatment of diverse cancers. The strategies encompass the utilization of immune checkpoint inhibitors and Chimeric Antigen Receptor T-Cell Immunotherapy (20, 21). Analyzing the characteristics of cells within the tumor microenvironment to predict immune infiltration is crucial for exploring new immunization strategies and studying responses to existing immune checkpoint inhibitors (22, 23). Recent research has categorized the microenvironments of tumors in cancer patients into three fundamental immune profiles: tumors that are immune inflamed (“hot”), immune excluded, and immune desert (“cold”). These profiles suggest different treatment options, excluding esophageal cancer, and provide valuable insights for effective therapeutic interventions (24, 25). To conclude, a meticulous and all-encompassing examination of the esophageal cancer tumor microenvironment, coupled with the determination of the corresponding tumor immunophenotype, can prove to be a valuable approach in directing immunotherapy and forecasting its effectiveness (20, 21, 23).

Several studies have substantiated the significant involvement of m⁶A modification in the development of tumor microenvironment (TME) diversity and complexity, a phenomenon that cannot be entirely elucidated by the RNA degradation mechanism (26). m⁶A modulators affected inflammation infiltrates and neovascularization of tumor tissues in human abdominal aortic aneurysm samples, where the markers METT14, FTO, and YTHDF3 are strongly colocalized with CD45+ leukocytes and CD3+ T cells, as well as CD68+ macrophages (27). Similarly, a METTL3/FTO-m⁶A methylation-mediated generation of M1/M2 macrophages from murine bone marrow-derived macrophages (BMDMs) has been described (28, 29). A new class of drugs targeting DNA methyltransferases (DNMTs) has been shown to successfully restore coordinated immune responses in solid tumors by triggering MHC 1 and interferon (IFN)-triggered immune-related signaling (30, 31). However, most studies, which are constrained by the state of technology, focus on just one or two m⁶A regulators, which is insufficient to describe the intricate functions of regulators in tumors. These research were made feasible by the ongoing development and collection of transcriptomics and genomic data, which offer a wealth of tools and resources for the study of m⁶A regulators and immune modulation (32).

In current study, we conducted a comprehensive analysis of m⁶A modifications and identified two distinct patterns of modifications, termed m⁶A clusters. These patterns were associated with different survival advantages and exhibited characteristics relevant to the TME, immune cell infiltration, and transcriptome (33). The observation that the TME characteristics linked with each m⁶A modification pattern closely corresponded with the manifestations and features of immune exclusion and immune inflammation, respectively, was of significant interest. This indicates a significant influence of m⁶A modification on individual tumor microenvironments (34, 35). Furthermore, a scoring system was devised to evaluate singular m⁶A modifications, facilitating the prediction of prognosis and the efficacy of immunosuppressive therapy. The strong correlation between m⁶A modifications and TME immune cell infiltration suggests that these modifications could serve as important prognostic markers and guide immunotherapy decisions in esophageal cancer.

2 Materials and methods

2.1 Esophageal cancer data acquisition and preprocessing

The Supplementary Figure S1 depicts the workflow employed in this study. The esophageal cancer samples’ transcriptional and clinical feature data were obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) databases and The Cancer Genome Atlas (TCGA, 2022.12.01, https://portal.gdc.cancer.gov/). Two distinct cohorts of esophageal cancer (ESCA), namely TCGA-ESCA and GSE13898, were used for further analysis. The RNA-Seq data obtained from the TCGA cohort underwent additional processing, resulting in the conversion of the data into transcripts per kilobase (TPM). Retrieve the normalization matrix file from the GEO database and employ R’s SVA package to address batch effects across distinct datasets. Obtain the survival duration and survival outcome data of two cohorts, with the exclusion of samples with survival periods less than 31 days and incomplete survival information. The Cancer Genome Atlas database was utilized to obtain somatic mutations, and copy number variation data for esophageal cancer were obtained from the UCSC Xena database (http://xena.ucsc.edu/) (36).

2.2 Classification according to 23 m⁶A regulators

These regulators include eight writers (METTL3, METTL14, METTL16, RBM15, RBM15B, VIRMA, WTAP and ZC3H13), 13 readers (FMR1, HNRNPC, HNRNPA2B1, IGFBP1/2/3, LRPPRC, RBMX, YTHDC1/2 and YTHDF1/2/3), and two erasers (FTO and ALKBH5). These modulators have been reported to affect or modulate the performance of RNA (Supplementary Figure S2). The expression levels of these 23 m⁶A regulators were utilized for unsupervised clustering analysis to identify distinct subtypes of m⁶A methylation modifications and classify patients for further analysis. The consensus clustering technique, implemented with the R package ConsensusClusterPlus, was utilized to calculate the number of clusters and assess their stability (37, 38).

2.3 scRNA-seq data processing

We analyzed the dataset GSE196756 about Esophageal Squamous Cell Carcinoma (ESCC) cells from the GEO repository (39). The data were sourced from Homo sapiens, with the data platform being GPL24676. We picked specific ESCC samples (GSM5900215,GSM5900216,GSM5900217,GSM5900218,GSM5900219,GSM5900220) for analysis. The R package: “Seurat” was used to analyze the transcript count matrix for quality control and preliminary data exploration (40). The filtering threshold was set as follows: Excluding genes detected in less than 3 cells, excluding cells with < 200 genes detected, Excluding cells with > 20% mitochondrial gene expression. We addressed batch differences using LogNormalize, Harmony and Principal Component analysis (PCA) helped us cluster cells based on variable genes via Seurat’s “FindNeighbors” and “FindClusters” functions. Uniform t-distributed Stochastic Neighbor Embedding (t-SNE) helped visualize this.

2.4 Estimation of immune infiltrating cells in TME

The R software package GSVA was utilized to conduct GSVA enrichment analysis to look at variations in m⁶A modification patterns in biological processes. The GSVA technique is a nonparametric and unsupervised approach that is predominantly employed to assess alterations in the activity of biological processes and pathways within samples (41). The gene sets from “c2.cp.kegg.v7.2.symbols” were downloaded from the MSigDB database for performing GSVA analysis. The present study utilized Single Sample Gene Set Enrichment Analysis (ssGSEA) in the R software package GSVA to assess the infiltration rates of 24 immune cells across various m⁶A regulator clusters. The differences between different m⁶A regulator clusters were assessed using the Wilcox test, and survival analysis was conducted to examine their association with patient outcomes.

2.5 Gene expression differences among phenotypes modified with m⁶A

Using a consensus clustering algorithm, we were able to divide esophageal cancer patients into two distinct subtypes according to m⁶A regulator expression. We revealed that the relationship between the two m⁶A clusters and immune landscape by CIBERSORT, EPIC, MCPCOUNTER, QUANTISEQ, TIMER, and XCELL algorithms (42). The differentially expressed genes (DEGs) between these two m⁶A-modified clusters were subsequently identified using the Limma package. The significance criterion for determining differential genes was set at a p-value < 0.05.

2.6 Differentially expressed genes enriched in functional pathways and functions

An important bioinformatics tool for gene annotation and analysis is the Gene Ontology (GO). It encompasses three categories: cellular component (CC), biological process (BP), and molecular function (MF). The Kyoto Encyclopedia of Genes and Genomes (KEGG) database serves as an integrative platform for genomic, chemical, and system function data, enabling the correlation of gene catalogs with higher-level system functions across various levels, including the cell, species, and ecosystem. To annotate the DEGs and gain insights into their biological functions, we utilized the clusterProfiler package, a widely used R package for functional enrichment analysis. The clusterProfiler package offers convenient functions to perform GO and KEGG enrichment analyses. For the study to be meaningful, the p-value must be less than 0.05 and the q-value must be less than 0.05.

2.7 The construction of the m⁶A score

We created a scoring system for m⁶A based on PCA to measure the patterns of m⁶A change in specific esophageal cancer patients. Genes demonstrating significant prognostic effects were selected from the different m⁶A modification clusters, based on which clustering of samples and construction of m⁶A scores were performed using a univariate Cox regression model. We determined the number of gene clusters and ensured stability using the consensus clustering algorithm. Marker scores for m⁶A-related genes were generated using PCA, and the first and second principal components were extracted as the marker scores. The method emphasizes the scores based on the collective behavior of highly correlated or inversely correlated genes within significant gene clusters, while minimizing the impact of genes that do not align with other members of the cluster. PCA is a dimensionality reduction method typically used to reduce the dimensionality of a dataset by transforming a large number of variables into fewer variables that still contain most of the information in the set (43, 44). We used the following method to define m⁶A scores: m⁶Ascore = Σ(PCA1i + PCA2i), the variable “i” denotes the final gene expression linked to the m⁶A phenotype (45, 46).

2.8 Evaluate the m⁶A scoring model

To evaluate the clinical applicability and reliability of the m⁶A score, receiver operating characteristic (ROC) curves were utilized to predict the outcomes at 1 year, 3 years, and 5 years. Initially, the ROC curve was constructed using all samples, followed by a separate analysis focusing on the TCGA-ESCA cohort to compare the prognostic predictive performance of the m⁶A score against other clinical variables. Correlations between the m⁶A score, clinical variables, and prognosis were examined using both univariate and multivariate Cox regression analyses. The purpose of the study was to examine the potential of the m⁶A score as a standalone predictive marker for esophageal cancer. Significance at the p < 0.05 level is usually used to determine statistical significance in a forest plot diagram. Furthermore, a nomogram was constructed using eight indicators (age, gender, tumor grade, stage T, N, M, pathologic stage, and m⁶A score) to anticipate the patient’s 1-, 3-, and 5-year survival rates. The predictive performance of the nomogram was evaluated using ROC curves. The R packages timeROC, rms, survival and survminer were employed for the necessary calculations and graphical representation.

2.9 Data research on genome mutations

The frequency of copy number variation (CNV) for the 23 m⁶A regulators in the TCGA-ESCA cohort was computed to assess the occurrence of copy number increases or losses. Copy number variation plots were generated using the R package Circos to visualize CNV patterns of m⁶A regulators of human chromosomes. The Tumor Mutation Burden (TMB) was computed by quantifying the aggregate count of nonsynonymous mutations present in the TCGA-ESCA cohort. The R package maftools was employed to create an oncoprint, which visually represents the gene mutation landscape. Using these approaches, the copy number variation map and oncoprint provide insights into the copy number alterations and mutation profiles of the m⁶A regulators in esophageal cancer based on the TCGA-ESCA cohort.

2.10 Tumor immune dysfunction and exclusion prediction and IC50 estimation

The Tumor Immune Dysfunction and Exclusion (TIDE) model, developed by researchers (47, 48), is used to evaluate the clinical efficacy of immune checkpoint inhibition therapy. The TIDE model provides prediction scores that reflect the likelihood of a patient’s response to immune checkpoint inhibition. Higher TIDE prediction scores are associated with a poorer response to immune checkpoint inhibition therapy. This model helps clinicians and researchers assess the potential effectiveness of immune checkpoint inhibition in individual patients.

2.11 Collect critical information for ICI-based cohorts

A systematic search was conducted to identify publicly available gene expression profiles of patients undergoing immune checkpoint inhibitor (ICI) therapy. The search aimed to identify datasets that included detailed clinical and pathological information. Ultimately, three immunotherapeutic cohorts were included in our study: metastatic melanoma patients treated with nivolumab (anti-PD-1 monoclonal antibody) (49) or ipilimumab (anti-CTLA-4 monoclonal antibody) (50), and patients who have been diagnosed with metastatic urothelial carcinoma (mUC) and have received treatment with the anti-PD-L1-targeting drug Atezolizumab (51). The ESCA-specific immunotherapy-treated cohort GSE165252 was found (n=45 ESCA patients treated with anti-PD-1 monoclonal antibody Atezolizumab), which contains the binary information on immune therapy response (response and non-response groups). We curated gene expression profiles from pre-therapy biopsy samples and transformed them into TPM (Transcripts Per Million) format. These datasets provide valuable information for our study on the response to ICI therapy and associated gene expression patterns.

2.12 Sensitivity analysis of anticancer drugs

For the study of molecular therapies for cancer and gene mutations, relevant data from the Genomics of Cancer Drug Sensitivity (GDSC) database were downloaded (52). This database offers a valuable resource for studying drug sensitivity in various cancer types. We utilized the pRRophetic package to obtain cell line gene mutation data and IC50 values associated with various anticancer drugs from GDSC, allowing us to analyze the correlation between patients with high- and low-risk scores and their sensitivity to different anticancer drugs. Through this analysis, we were able to examine the correlation between patients exhibiting high- and low-risk m⁶A scores and their responsiveness to a diverse array of anticancer medications (53). By leveraging these resources, we aimed to gain insights into the association between m⁶A modification patterns and the response to specific anticancer therapies, further enhancing our understanding of personalized cancer treatment approaches.

2.13 Cell transfection and cell line establishment

Esophageal carcinoma cell lines KYSE510 and TE-1 were cultured in RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS) at a temperature of 37 degrees Celsius and under a 5% CO₂ atmosphere. In preparation for cell transfection, these cell lines were seeded into 6-well culture plates and incubated overnight to allow for attachment and initial growth. On the subsequent day, once the cells had reached a confluence of 20%-30%, transfection was performed with siRNAs at a final concentration of 50 nmol/L using Lipofectamine 2000 (Invitrogen), following the protocol provided by the manufacturer.

2.14 Explore and validate potential oncogenic functions of RBMX In ESCA

Supplemental Experimental Procedures include the following information: Western blot for protein expression, Plate clone formation assay, EdU assay for cell proliferation detection, Wound healing assay for assessing cellular migration, Transwell migration/invasion assay (Supplementary Data Sheet S1: Supplemental Experimental Procedures).

2.15 Statistical analyses

Statistical analyses and graphical representations were conducted using R version 4.3.1. The Wilcoxon rank sum test, a statistical method, is useful for assessing and contrasting dissimilarities between two groups. The correlation between m⁶A regulators and prognosis was assessed with univariate Cox regression models and Kaplan-Meier survival analysis. The selection of cutoff points for the m⁶A score was performed by repeatedly testing all possible cutoffs using the survminer package in R, aiming to identify the maximum rank statistic. Partition the samples into groups based on their m⁶A scores, with one group consisting of high scores and the other of low scores. Prognosis was assessed using the Kaplan-Meier method, and the log-rank test was assessed for both cohorts. At the same time, there are other statistical methods for targeted analysis. Heatmaps were generated using the pheatmap package in R. All statistical tests were two-tailed, and a p-value less than 0.05 was considered statistically significant.

3 Results

3.1 Mutation of m⁶A regulators, immune infiltration, and construction of the prognostic landscape

Our study included 23 m⁶A regulators. Firstly, we calculated the frequency of mutations in the 23 regulators in ESCA. The 23 m⁶A regulators exhibited low mutation frequencies, with only 23 (12.5%) out of 184 ESCA samples from the TCGA database showing genetic alterations. Mutation information for each gene in each sample was presented in the waterfall plot, with different colors and specific annotations at the bottom representing the different mutation types. Interestingly, the Oncoplot analysis revealed that ZC3H13 displayed the highest mutation rate, predominantly characterized by nonsense mutations, while YTHDC2 had a mutation frequency of 2% (Figure 1A).

Figure 1

Figure 1. Landscape of genetic and expression variation of m⁶A regulators in esophagus cancer (A) Mutations of 23 m⁶A regulators in the TCGA-ESCA cohort. Each column represented individual patients. The upper barplot showed TMB, the number on the right indicated the mutation frequency in each regulator. The right barplot showed the proportion of each variant type. The stacked barplot below showed fraction of conversions in each sample. (B) In an TCGA-ESCA cohort, we looked at the CNV mutation rates of 23 m⁶A regulators. The findings were represented by red and green dots to represent increased and absent frequencies, respectively. (C) The precise chromosomal locations of CNVs in m⁶A regulators across all 23 chromosomes. (D) Interactions and prognostic implications of 23 m⁶A regulators in ESCA. The three types of m⁶A regulatory modifiers are represented by different colors: eraser in red, reader in orange, and writer in gray. The size of the circles corresponds to the prognostic relevance of each m⁶A modulator. The lines connecting the regulators indicate their interactions, with positive correlations in pink and negative correlations in blue. Prognostic risk factors are highlighted in purple, while prognostic protective factors are shown in green. (E) The expression of 23 m⁶A regulators between normal tissues and gastric tissues. Tumor, red; Normal, blue. The upper and lower ends of the boxes represented interquartile range of values. The lines in the boxes represented median value, and red and blue dots showed outliers. The asterisks represented the statistical p value (*P < 0.05; **P < 0.01; ***P < 0.001).

Moreover, our analysis of copy number variations (CNVs) in the 23 m⁶A regulators highlighted the prevalence of CNV mutations in ESCA. Notably, YTHDC1, FMR1, VIRMA, YTHDF1, METTL3, WTAP, HNRNPA2B1, RBMX, HNRNPC, and IGFBP1 exhibited a high frequency of CNV amplification, while HNRNPC, RBM15, YTHDF2, IGFBP2, and RBM15B displayed extensive CNV deletions (Figure 1B). The chromosomal alterations associated with these CNVs are visually depicted in Figure 1C. In order to evaluate the influence of m⁶A regulators on patient prognosis, a Kaplan-Meier survival analysis was conducted. The findings demonstrated significant correlations between the prognosis of ESCA patients and 8 m⁶A regulators (Supplementary Figure S3). Additionally, seventeen modifiers with prognostic value in ESCA patients were identified using univariate Cox regression analysis (Supplementary Table S1). Furthermore, our analysis of the m⁶A network revealed the intricate interactions, connectivity, and prognostic significance of m⁶A regulators in ESCA (Figure 1D). Our findings indicate noteworthy correlations between the expression levels of m⁶A regulatory factors within the same functional class, as well as significant associations among the three distinct types of regulatory factors. This interplay is likely to contribute to the generation of distinct m⁶A modification patterns, which play a crucial role in the initiation and progression of cancer.

As well, we identified CNV alterations as a potential underlying cause of disrupted expression of m⁶A regulatory factors. To further investigate this, we compared the gene expression levels of the 23 m⁶A regulators between normal and tumor tissues (54). In ESCA tissues, m⁶A regulatory factors’ expression exhibiting CNV amplification (such as METTL3, WTAP, VIRMA, YTHDC1, YTHDF1, HNRNPC, FMR1, and HNRNPA2B1) was significantly higher than in normal esophageal tissues, while the expression of IGFBP2 was lower (Figure 1E). Collectively, these analyses underscore the noteworthy diversity in the genetic and expression profiles of m⁶A regulators detected between normal and ESCA specimens. These findings emphasize the critical role of dysregulated expression of m⁶A regulators in the development and progression of ESCA.

3.2 scRNA-seq analysis

Single-cell RNA sequencing (scRNA-seq) of 25,796 immune and 8,197 non - immune cells from three primary tumor and paired normal samples in ESCA patients was generated by using 10x Genomics platform. Before filtration, there were 33993 cells in the 6 ESCA samples. For GSE196756, counts were normalized and technical covariates (mitochondrial percentage) were regressed out using the LogNormalize method (default settings), and batch effects across samples (6 ESCA patients) were corrected for using Harmony with theta = 2 to preserve biological variance (55). We then performed data normalization and quality control, and finally selected the top 2000 highly expressed and variable genes for further analysis. PCA used to reduce the dimensionality of the data showed no clear tendency for cells to separate. Nonlinear dimensionality reduction was performed using the t-SNE algorithm, which successfully clustered the cells into 13 clusters (Figure 2A). We then annotated all clusters and identified 9 cell types (Figure 2B). Furthermore, the expression levels of 23 m⁶A modulators were most abundant in B cells and T cells (Figure 2C). WTAP, ZC3H13, YTHDC1, HNRNPC, HNRNPA2B1 and RBMX are expressed in most cell types.

Figure 2

Figure 2. scRNA-seq data analysis. (A) The t-SNE algorithm divided the cells into 13 clusters by principal components. (B) The tSNE plot revealing 13 clusters was annotated into 9 different cell types. (C) The expression of 23 m⁶A regulators in 9 cell types.

3.3 Twenty-three regulator-mediated isoforms of m⁶A methylation

Using the ConsensusClusterPlus R software package, we performed patient classification based on the expression of the 23 m⁶A regulators to delineate distinct m⁶A-modified subtypes. Our analysis revealed two subtypes: subtype A consisting of 111 cases and subtype B consisting of 75 cases (Figures 3A–D and Supplementary Table S2). Notably, patients belonging to m⁶A regulator group B exhibited significantly longer survival compared to those in m⁶A regulator group A (P = 0.019, Figure 3E). We generated a heatmap to visualize the expression patterns of the 23 m⁶A regulators in the two m⁶A-modified subtypes (Figure 3F).

Figure 3

Figure 3. Identification of m⁶A methylation modification subtypes. (A) Consensus clustering and generated a CDF with the number of subtypes ranging from k = 2 to 9. (B) The heat map of sample clustering under k = 2 in 2 independent ESCA cohorts. (C) The relative change in the area under the CDF curve for values of k ranging from 2 to 9. (D) Principal component analysis of transcriptome profiles of two m⁶A modification patterns (E) Survival analyses for the two m⁶A modification patterns based on 186 patients with esophagus cancer from TCGA-ESCA and GEO cohorts (GSE13898) including 111 cases in m⁶Acluster-A, 75 cases in m⁶Acluster-B, Kaplan-Meier curves with Log-rank p value 0.019 showed a significant survival difference among two m⁶A modification patterns. The overall survival rate of cluster B in m⁶A cluster A and B subclusters is better. (F) Unsupervised clustering of 23 m⁶A regulators in two cohorts with heatmap analysis of m⁶A cluster, tumor stage, survival status, and age. Red is high expression, blue is low expression.

3.4 TME cell infiltration characteristics in distinct m⁶A modification patterns

To gain insights into the underlying biomolecular signatures associated with the different m⁶A-modified phenotypes, we integrated the expression profiling data of both TCGA-ESCA and GSE13898 cohorts and performed differential expression analysis using the Limma R software package. This analysis identified 2599 DEGs, which were subsequently annotated using the clusterProfiler R package. The DEGs were found to be enriched in several important biological processes, including T cell activation, regulation of immune effector process, neutrophil-mediated immunity, mesenchyme development, mesenchymal cell differentiation, leukocyte transendothelial migration, chemokine signaling pathway, and VEGF signaling pathway (Figures 4A, B, Supplementary Tables S3, S4).

Figure 4

Figure 4. Comparison of the enrichment analysis for immune cells and immune pathways between two m⁶A clusters. (A) Functional annotation of the genes with different expressions between cluster A and cluster B using GO terms. (B) Pathway of KEGG enrichment of DEGs between two m⁶A clusters. (C) The abundance of each TME infiltrating cell in two m⁶A modification patterns. The upper and lower ends of the boxes represented interquartile range of values. The lines in the boxes represented median value, and the dots showed outliers. The asterisks represented the statistical p value (*P < 0.05; **P < 0.01; ***P < 0.001). (D) The heatmap was used to visualize these KEGG enrichment pathways, and blue represented activated pathways and yellow represented inhibited pathways. (E) The heatmap was used to visualize these immune cells. (F) The heatmap was used to visualize these immune-related pathways. (G) The box plot figure demonstrates the enrichment scores for clusters A (red) and B (yellow) across several biological processes, highlighting statistically significant differences (*P < 0.05; **P < 0.01; ***P < 0.001, ****P < 0.0001).

In order to examine the biological alterations linked to diverse m⁶A modification patterns, a comparative analysis of immune cell composition in the TME was performed. The findings indicate that the A subcluster exhibited a higher degree of infiltration by memory B cells, immature B cells, T helper 1 (Th1) cells, activated memory CD4+ T cells, and regulatory T cells (Treg). On the other hand, m⁶A cluster B exhibited significantly increased infiltration of natural killer cells and neutrophils (Figure 4C).

We further employed GSVA enrichment analysis to gain insights into the biological activity associated with these distinct m⁶A modification patterns. The results of our observations indicate that m⁶A cluster A exhibits noteworthy enrichment in pathways linked to stroma and cancer metastasis, including ECM-receptor interaction, focal adhesion, and others. On the other hand, m⁶A cluster B showed enrichment in metabolic pathways such as histidine metabolism, fatty acid metabolism, propanoate metabolism, glycolysis, fructose, and mannose metabolism (Figure 4D).

Interestingly, GSVA enrichment analysis revealed that m⁶A cluster A exhibited significant enrichment in adaptive immune cell infiltration, encompassing memory B cells, activated memory CD4+ T cells, immature B cells, Th1 cells, regulatory T cells (Treg), and stromal activation (Figures 4C–G, 5E). Surprisingly, despite the higher immune cell infiltration, patients with this m⁶A modification pattern did not demonstrate a survival advantage (Figure 4F). Prior research has detected an immune-excluded phenotype within tumors, wherein immune cells exist in the stroma encircling nests of tumor cells, yet are unable to penetrate the tumor parenchyma. T-cell suppression is known to occur when the stroma in the TME is activated. Hence, our speculation is that the stromal activation observed in cluster A suppresses the antitumor effect of immune cells in patients with ESCA. The aforementioned conjecture was subsequently substantiated through analyses that demonstrated a marked increase in stromal activity within cluster A, which encompassed the activation of epithelial-mesenchymal transition (EMT), transforming growth factor beta (TGF-β), and WNT pathways, all of which were found to be statistically significant (Figure 4G).

Figure 5

Figure 5. Construction of m⁶A score. (A) Differences in m⁶A scores between two m⁶A subclusters. (B) Differences in m⁶A scores between the three gene clusters. (C) The Sankey diagram illustrates the association between m⁶A score, m⁶A clusters, gene clusters, and survival outcomes. (D) Correlations between the m⁶A score and tumor-infiltrating immune cells using Spearman’s analysis. The positive and negative correlations are marked with red and blue, respectively. (E) Correlations between m⁶A score and the known biological gene signatures using Spearman analysis. The negative correlation was marked with blue and positive correlation with red. (F) Kaplan-Meier curve showing overall survival probability between high and low m⁶A score groups.

We have integrated immune deconvolution tools such as CIBERSORT, EPIC, MCP_COUNTER, QUANTISEQ, TIMER and XCELL to distinct immune microenvironments characterize two m⁶A clusters (Supplementary Figure S4). Comparative analysis of immune infiltration patterns between the two m⁶A clusters revealed significant heterogeneity. Cluster A exhibited higher infiltration of immunosuppressive regulatory T cells (Tregs, P < 0.01 by CIBERSORT/QUANTISEO) and exhausted CD8+ T cells (PD-1+Tim-3+, P < 0.05), whereas cluster B showed elevated cytotoxic CD8+ T cells (Granzyme B+, P < 0.001 by TIMER). Pro-tumor M2 macrophages were enriched in cluster A (P < 0.001 across CIBERSORT/QUANTISEO/XCELL), while cluster B had higher M1 macrophages (P <0.05), suggesting divergent macrophage polarization states. CAFs were markedly increased in cluster A (P < 0.001 by EPIC/MCP-counter), correlating with elevated ECM remodeling scores (e.g., collagen cross-linking, P= 0.002). These findings were robust across multiple deconvolution algorithms (CIBERSORT, EPIC, XCELL).

Based on the comprehensive analyses conducted, it is intriguing to note that the two m⁶A modification patterns exhibit distinct characteristics in terms of TME cell infiltration. Cluster A is associated with an immune-excluded phenotype, characterized by the infiltration of adaptive immune cells and stromal activation. On the other hand, cluster B corresponds to an immune-inflamed phenotype, characterized by the infiltration of innate immune cells and metabolic reprogramming (Figures 4C–G). These findings suggest that m⁶A methylation modifications may be involved in tumor metabolism, EMT, immune regulation, and have close associations with tumor initiation and progression.

3.5 Characteristics of clinical and transcriptome traits in m⁶A-related phenotypes

Despite the successful categorization of ESCA patients into two subtypes through a consistent clustering algorithm utilizing m⁶A regulator expression, the genetic alterations responsible for these phenotypes and their prognostic implications remain inadequately comprehended. To gain deeper insights, we conducted univariate Cox regression analysis on the 2599 DEGs identified between the previously established m⁶A clusters. A total of 80 survival-related genes were identified through this analysis, which we referred to as the m⁶A-related signature genes (Supplementary Table S5). Through unsupervised clustering analysis using representative m⁶A-associated marker genes, we identified three stable transcriptomic phenotypes, denoted as gene clusters A, B, and C (Figures 6A–C; Supplementary Table S6). The predictive importance of these gene subgroups was then investigated by fusing transcriptome data with survival data. Based on Kaplan-Meier analysis and log-rank test, it was observed that patients assigned to gene cluster B displayed a favorable prognosis (Figure 6D). A heat map was generated to visually depict the clinical characteristics of 80 m⁶A -related signature genes and the expression of m⁶A subgroups in the three gene clusters (Figure 6E). Notably, the three m⁶A gene clusters exhibited significant differential expression of m⁶A regulatory factors, which aligns with the methylation modification process and supports the predicted effects of m⁶A (Figure 6F).

Figure 6

Figure 6. Consensus clustering of m⁶A -related gene subtypes (A) Consensus matrices of 80 m⁶A phenotype-related genes according to TCGA and GEO cohort for k = 3. (B) Consensus clustering CDF with the number of subtypes k = 2 to 9. (C) Consensus clustering cumulative distribution function k = 2 to 9. (D) Kaplan-Meier overall survival curves for patients in three m⁶A-related gene clusters. (E) Heatmap showing the correlation between the expression levels of the DEGs derived from 3 m⁶A clusters and sex, age, m⁶A clusters, tumor stage, survival status and gene clusters. Red is high expression, blue is low expression. (F) The expression of 23 m⁶A regulators in three gene cluster. The upper and lower ends of the boxes represented interquartile range of values (**P < 0.01; ***P < 0.001).

3.6 Establishment of m⁶A score and its association with tumor microenvironmental features

While previous analyzes have yielded valuable insights into the impact of m⁶A methylation on immune cell infiltration status and tumor prognosis, accurate prediction of m⁶A methylation patterns in individual patients remains a challenging task. To address this challenge, the PCA score was employed to compute the m⁶A score, which also provides a quantitative assessment of the modified m⁶A landscape in patients with ESCA. Figure 5A illustrates that patients in m⁶A cluster B exhibit lower m⁶A scores compared to those in m⁶A cluster A, and Figure 5B demonstrates that patients in gene cluster B have lower m⁶A scores than those in gene clusters A and C. We have depicted the process of m⁶A score construction in a Sankey diagram (Figure 5C). In order to evaluate the association between the m⁶A score and tumor-infiltrating immune cells, a Spearman’s analysis was conducted and the outcomes were presented in a heatmap (Figure 5D), revealing a positive correlation between the m⁶A score and the presence of immune cells within the tumor microenvironment. Additionally, we examined the correlation between the m⁶A score and known signal pathway signatures. The resulting correlation matrix heatmap demonstrated that the m⁶A score exhibited significant positive associations with signatures related to EMT, stromal activity, DNA repair, antigen processing machinery, and the TGF-β pathway (Figure 5E). In addition, we conducted an evaluation of the prognostic relevance of the m⁶A score. Through implementation of the Kaplan-Meier survival analysis, it was determined that patients exhibiting low m⁶A scores experienced a more favorable prognosis in contrast to those with high m⁶A scores (Figure 5F). This indicates that the implementation of the m⁶A-score-based computation proficiently delineates the prognosis of patients.

In addition, we also found that in T0, T1-2, or T3–4 stage, N0, N1-2, or N3 stage, M0 or M1 stage, male or female, young or old patients, and patients, lower m⁶A score showed more significant survival advantage, which means that m⁶A score can also be used to access various clinical features of patients, such as age, gender, or clinical stage subgroup (Supplementary Figure S5).

3.7 Verification and clinical evaluation of m⁶A score

To validate the m⁶A score, we conducted ROC curve analysis for 1-year, 3-year, and 5-year intervals and calculated the corresponding area under the curve (AUC) values. The results showed that all three ROC curves in the total sample cohort (Figure 7A) and the separate TCGA-ESCA cohort (Figure 7B) showed AUC values exceeding 0.67. Furthermore, when comparing the m⁶A score with other clinical features, the AUC value of the m⁶A score was found to be the highest (Figure 7C).

Figure 7

Figure 7. Verification and clinical evaluation of m⁶A score. (A) For all samples, the area under the curve (AUC) values of 1-, 3-, and 5-year ROC curves for m⁶A scores have been calculated. (B) For TCGA-ESCA cohort, the area under the curve (AUC) values of 1-, 3-, and 5-year ROC curves for m⁶A scores have been calculated. (C) Comparison of 1-year ROC curves of the m⁶A score model with other clinical features. (D) Univariate COX regression analysis of clinicopathological parameter and m⁶A score. (E) Multivariate Cox regression analysis of clinicopathological parameter and m⁶A score. (F) The nomogram is used to forecast the probabilities of 1-year, 3-year, and 5-year survival rates. (G) The AUC values of the nomogram’s 1-year, 3-year, and 5-year ROC curves.

The findings of the univariate Cox regression analysis indicate that the stage, stage M, stage N, and m⁶A score possess prognostic potential (Figure 7D). Additionally, the multivariate Cox regression analysis reveals that both the stage and m⁶A score exhibit independent prognostic value (Figure 7E). To quantitatively assess individual risks in the clinical setting, the integration of multiple clinical indicators can be achieved through a nomogram. In this study, we developed a nomogram for predicting patients’ overall survival (OS) at 1-year, 3-year, and 5-year intervals (Figure 7F). The predictive performance of the nomogram was evaluated using ROC curve analysis. The present study determined the AUC values of the ROC curves for 1-year, 3-year, and 5-year intervals to be 0.784, 0.831, and 0.801, respectively (Figure 7G). These results indicate that the m⁶A score may serve as a promising clinical predictor and, when integrated with other clinical factors, could potentially improve the prognostic precision and clinical outcomes for patients diagnosed with ESCA.

3.8 Somatic variation correlates with m⁶A score

The potential of TMB as a tumor marker for immune checkpoint therapy in patients has been demonstrated. Given the clinical significance of TMB, an analysis was conducted to investigate the genetic characteristics within each subgroup, as defined by the m⁶A score, and their association with TMB. Patients were divided into two subgroups based on TMB.

Based on the results depicted in Figures 8A, B, it was observed that both TP53 (86% vs. 71%) and TTN (44% vs. 32%) exhibited a higher rate of somatic mutation in the group with a high m⁶A score, suggesting a potential association with the poorer prognosis observed in this group (Figure 8C). Subsequently, we assessed the combined prognostic value of these scores in stratifying ESCA patients. Survival analysis revealed that the TMB status did not influence predictions based on the m⁶A score, consistently demonstrating a survival advantage in the low m⁶A score group (Figure 8D). The results of this study contribute to a more thorough comprehension of the impact of the m⁶A score on genomic variability, presenting innovative perspectives for investigating potential associations between m⁶A methylation modification and somatic mutations. These findings demonstrated that distinct m⁶A modification patterns significantly influenced tumor immune phenotypes and may serve as predictive biomarkers for anti-PD-1/PD-L1 immunotherapy response efficacy. It has also been revealed that the m⁶A score is indirectly used to predict the success of immunotherapy.

Figure 8

Figure 8. Correlation between m⁶A score and TMB and TIDE in the ESCA cohort. (A, B) OncoPrint for gene mutations in high and low m⁶A score groups. in the high m⁶A score group. (C) Kaplan-Meier curve showing overall survival probability between high TMB group and low TMB group (P < 0.05). (D) TMB and m⁶A scores were used in a stratified survival analysis. (E) Violin plot of differential expression of other immune checkpoints between groups with high and low m⁶A scores. (*P < 0.05, **P < 0.01,***P < 0.001) (F) TIDE differences in the TCGA cohort between high and low m⁶A score groups. (G) TIDE differences in the GSE13898 cohort between high and low m⁶A score groups.

3.9 M⁶A score predicts the possibility of benefit from immunotherapy

Subsequently, the differences in the levels of other immune checkpoints between the high and low m⁶A score groups were compared. The high m⁶A score group had higher expression of CTLA4, CD70, TNFSF14, ICOS, CD80, TNFRSF9, HAVCR2, CD200, NRP1, TNFSF15, TNFSF4, CD40, TNFRSF14, LGALS9, CD86, ADORA2A, and CD28, while the low m⁶A score group had higher expression of BTLA (Figure 8E).

The use of ICI therapy, specifically CTLA-4/PD-1 inhibitors, has resulted in a significant advancement in antitumor treatment. Alongside established predictors such as TMB, PD-L1, and MSI (56, 57), newly discovered indicators such as TIDE are extensively utilized and highly recommended for assessing immune response. Our analysis further demonstrated a noteworthy reduction in TIDE within the low m⁶A score group, as evidenced by the TIDE distribution in TCGA-ESCA and GSE13898 (both P < 0.01) (Figures 8F, G). As a result of these findings, it is inferred that tumor m⁶A modification patterns play an important role in mediating immune responses in tumors.

Based on the significant correlation between m⁶A scores and immune responses, our subsequent investigation aimed to assess whether m⁶A modification signatures could serve as predictive markers for patient response to ICI therapy in three separate immunotherapy cohorts. Firstly, a high m⁶A score exhibited significantly shorter survival time (HR, 1.845 [95% CI, 1.254 to 2.714], P = 0.013, Figure 9A) and a markedly clinical response in an anti-PD-L1 therapy in a cohort of metastatic urothelial carcinoma (response rate, low vs. high m⁶A score, 53% vs. 19%, Figure 9B) (51). This result was also identified in both the anti-PD-1 cohort (49) and anti-CTLA-4 cohort (50). Patients belonging to the high m⁶A score group demonstrated noteworthy clinical drawbacks and a considerably reduced lifespan (anti-PD-1, HR, 2.886 [95% CI, 1.002 to 8.314], P = 0.018. (Figure 9C) anti-CTLA-4, HR, 2.141 [95% CI, 1.018 to 4.503], P = 0.035, Figure 9E). The significant therapeutic benefits and immune response to ICI treatment were confirmed in patients with a low m⁶A score compared to those with a high m⁶A score (anti-PD-1, response rate, low vs. high m⁶A score, 33% vs. 18%, Figure 9D; anti-CTLA-4, response rate, low vs. high m⁶A score, 40% vs. 32%, Figure 9F). The m⁶A score of the GSE165252 cohort was further validated, the significant therapeutic benefits and immune response to anti-PD-1 treatment were confirmed in patients with a low m⁶A score compared to those with a high m⁶A score (anti-PD-1, response rate, low vs. high m⁶A score, 39% vs. 18%) (Supplementary Figure S6). The m⁶A score is also associated with patient response to immunotherapy and can be used to predict patients’ prognoses. In conclusion, the m⁶A score serves as a promising prognostic indicator in ESCA and may also provide guidance for ICI treatment in clinical practice.

Figure 9

Figure 9. The m⁶A risk score predicts immunotherapeutic benefits. (A) Survival difference analysis of patients with high and low m⁶A risk score in the IMvigor210 cohort. P = 0.013. (B) Rate of clinical response to anti-PD-L1 immunotherapy in high or low m⁶A risk score groups in the IMvigor210 cohort. (C) Kaplan-Meier curves for high and low m⁶A risk score patient groups in the Riaz et al. cohort. Log-rank test, P = 0.018. (D) The fraction of patients with clinical response to anti-PD-1 immunotherapy (Riaz et al. cohort) in low or high m⁶A risk score groups. CR/PR vs. SD/PD: 33% vs. 67% in the low m⁶A risk score groups, 18% vs. 82% in the high m⁶A risk score groups. (E) Kaplan-Meier curves for high and low m⁶A risk score patient groups in the Vanallen et al. cohort. Log-rank test, P = 0.035. (F) The fraction of patients with clinical response to anti-CTLA-4 immunotherapy in low or high m⁶A risk score groups of Vanallen et al. cohort. CR/SD vs. PD: 40% vs. 60% in the low m⁶A risk score groups and 32% vs. 68% in the high m⁶A risk score groups. CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease.

3.10 Sensitivity analysis of patients with ESCA to different small molecule drugs based on m⁶A risk score

We performed an estimation of IC50 values and assessed the drug sensitivities of chemotherapeutic drugs for a cohort of 186 ESCA patients, utilizing data from the TCGA and GEO databases. The estimation process employed the “pRRophetic” R package, which utilized the expression profiles of the patients. Then, IC50 values were compared between the groups with high and low m⁶A scores. The IC50 values are utilized to assess the cellular response of various cell lines to a total of 138 distinct chemotherapeutic and small molecule anticancer drugs. The research found statistically significant differences (P < 0.05) between patients with high and low m⁶A risk scores in the IC50 values of several chemotherapeutic drugs and small molecule anticancer medicines. Notably, Bortezomib, Camptothecin, Cytarabine, Erlotinib, Gefitinib, Gemcitabine, Metformin, Methotrexate, and Paclitaxel exhibited particularly noteworthy differences (Figures 10A–I; Supplementary Figure S7).

Figure 10

Figure 10. Sensitivity of the m⁶A risk score to different chemotherapy drugs and small molecule anticancer drugs was analyzed based on the GDSC database. (A) Bortezomib. (B) Camptothecin. (C) Cytarabine. (D) Erlotinib. (E) Gefitinib. (F) Gemcitabine. (G) Metformin. (H) Methotrexate. (I) Paclitaxel.

3.11 RBMX’s impact on ESCC cell proliferation and migration

To establish the mechanistic link between m⁶A modification and malignant progression in ESCC, we prioritized RBMX for functional interrogation based on its central position in the m⁶A regulatory network. Bioinformatics analysis identified RBMX as a hub gene in protein-protein-interaction network and co-expressed with key m⁶A regulators (METTL3, FTO, YTHDF2). RBMX expression levels were quantified in the ESCC cell lines KYSE510 and TE-1, revealing a notable reduction in protein expression following RBMX knockdown (Figure 11A). The clone formation assay demonstrated that the knockdown of RBMX significantly impeded the proliferative capacity of ESCC cells (Figure 11B). EDU staining corroborated these findings, indicating a significant decrease in the proliferative activity of si-RBMX-transfected KYSE510 and TE-1 cells (Figures 11C, D). The wound healing assay further illustrated that, after 48 hours, the wound closure ability of si-RBMX-transfected KYSE510 and TE-1 cells was markedly diminished compared to the Si-NC control group (Figure 11E). Additionally, migration and invasion assays were conducted to evaluate the impact of RBMX on ESCC cell motility. The knockdown of RBMX in KYSE510 and TE-1 cells led to a significant reduction in both the invasive and migratory capabilities of the cells (Figures 11F, G).

Figure 11

Figure 11. Impact of altered RBMX expression on proliferation and invasion of esophageal cancer cells. (A) RBMX protein expression levels in Kyse510 and TE-1 cells. (B) plate clone formation experiments assays in transfected Kyse510 cells and TE-1 cells. (C) EDU assays in transfected Kyse510 cells (D) and TE-1 cells. (E) Wound healing assays in transfected Kyse510 cells and TE-1 cells. (F) Invasion assay in transfected Kyse510 cells and TE-1 cells. (G) migration assay in transfected Kyse510 cells and TE-1 cells (*P < 0.05; **P < 0.01; ***P < 0.001, ****P < 0.0001).

4 Discussion

Accumulating evidence from various studies emphasizes the important role of m⁶A methylation modification in the immune process of organisms. Further investigation is necessary to achieve a thorough comprehension of the immune cell infiltration within the TME in ESCA that is mediated by multiple m⁶A regulators. Therefore, it is crucial to clarify the characteristics of immune cell infiltration in relation to diverse m⁶A modification patterns. This will enhance our understanding of the TME and antitumor immune responses within it, and offer approaches for risk stratification and clinical management of patients with esophageal cancer. This study identified two distinct modification patterns with the assistance of 23 m⁶A regulators. The mRNA transcriptome differences observed between these patterns were found to be significantly associated with T cell activation, regulation of immune effector processes, neutrophil-mediated immunity, mesenchyme development, mesenchymal cell differentiation, leukocyte transendothelial migration, the Chemokine signaling pathway, and the VEGF signaling pathway (Figures 4A, B). The two patterns exhibited markedly different TME cell-infiltrating characteristics. Cluster A was categorized as an immune-excluded phenotype, marked by the infiltration of adaptive immune cells and stromal activation. On the other hand, cluster B was classified as an immune-inflamed phenotype, characterized by the infiltration of innate immune cells and metabolic reprogramming. The immune-inflamed phenotype, also known as “hot tumors,” is distinguished by substantial immune cell infiltration within the TME (24, 58, 59). Despite the presence of a significant number of immune cells in the immune-excluded phenotype, their distribution is limited to the stromal compartment surrounding the tumor cell nests, rather than infiltrating the tumor parenchyma. The stromal compartment may be localized to the tumor periphery or may extend into the tumor, potentially leading to the misinterpretation that immune cells are present within the tumor (60, 61). Consistent with the established definitions, our findings revealed that cluster A exhibited a pronounced stromal activation status, including elevated expression of EMT and TGF-β pathways (Figures 4D, E, G), which are associated with T-cell suppression. The observed TME cell-infiltrating characteristics in each cluster reinforce the validity of our immune phenotype classification based on distinct m⁶A modification patterns. Consequently, after comprehensively exploring the TME cell-infiltrating characteristics induced by distinct m⁶A modification patterns, it is not surprising that cluster A, despite having activated innate immunity, exhibited poorer prognosis.

The stromal activation in cluster A (e.g., TGF-β, EMT) and its link to immune exclusion are supported by recent studies on CAF subtypes (60–62). Cluster A exhibits an immunosuppressive stromal microenvironment owing to the enrichment of TGF-β and multiple EMT-related pathways and (Figures 4D, G). The tumor stroma, particularly CAFs and their remodeled extracellular matrix (ECM), plays a pivotal role in shaping the immunosuppressive TME by regulating T cell infiltration and function through both physical barrier and molecular mechanisms (60). The physical barrier prevents cytotoxic T cells from contacting cancer cells, creating an “immune-excluded” microenvironment (60). In addition to structural constraints, CAF heterogeneity further exacerbates immune evasion. Distinct CAF subpopulations may drive divergent stromal remodeling patterns: certain subsets promote the formation of rigid, cross-linked stroma that impedes T cell migration, while others secrete immunosuppressive factors (60, 61). Three functional subtypes of CAFs have been identified in non-small cell lung cancer (NSCLC) based on their heterogeneity (63). These functional disparities among CAFs are driven by their intrinsic TGF-β signaling. This CAF functional classification correlates with patients’ clinical responses to targeted therapies and is also associated with the tumor immune microenvironment (63). Notably, RNA modification “writers” (e.g., m⁶A/m¹A regulators) appear to influence CAF activation states, as evidenced by the association between high “Writers-Score”, poor prognosis, and suppressive immune infiltration (e.g., M2 macrophages, EMT) (62). These findings suggest that epigenetic reprogramming of CAFs may reinforce immune exclusion by coupling matrix stiffness with broader immunosuppressive signals, such as PD-L1 upregulation. Thus, stromal activation drives CAF heterogeneity and immune exclusion via coordinated ECM remodeling (e.g., collagen cross-linking, fibronectin deposition) and epigenetic reprogramming (e.g., m⁶A-mediated RNA stabilization of TGF-β signaling components). Therapeutically targeting these matrix-driven immunosuppressive mechanisms—such as through ECM-degrading enzymes (e.g., collagenase) or epigenetic inhibitors—could dismantle the stromal-T cell barrier, thereby enhancing the efficacy of T cell–mediated antitumor immunity. The study by Du et al. elucidates that RBMX stabilizes IL-33 mRNA through a liquid-liquid phase separation mechanism, thereby activating the TGF-β signaling pathway. This process orchestrates the bidirectional regulation of tumor plasticity and the immunosuppressive microenvironment, providing a theoretical foundation for developing precision therapeutic strategies targeting the RBMX/TGF-β axis (64). We hypothesize that RBMX, as an m⁶A reader, regulates CAF crosstalk by enhancing RNA stability of the IL-33/TGF-β axis and increasing stromal stiffness, while simultaneously suppressing immune-activating signals (e.g., CXCL10-STAT1) to impair T cell function (61). This dual mechanism aligns with clinical observations in cluster A patients, where despite high immune cell infiltration, T cells are predominantly confined to stromal regions and exhibit significantly reduced survival rates. Such an “immune-excluded” phenotype closely mirrors the stroma-mediated immune privilege phenomenon proposed by Joyce et al. (60).

The immune landscape analysis underscores how m⁶A modification patterns shape tumor-immune interactions (Supplementary Figure S4). M⁶A cluster A (immune-excluded): Dominated by Tregs, M2 macrophages, and CAFs, this phenotype aligns with TGF-β-driven stromal activation. The concomitant suppression of cytotoxic lymphocytes (evidenced by low CD8+/NK cell ratios, P < 0.001) may explain poorer immunotherapy responses observed in this subgroup. M⁶A cluster B (immune-inflamed): Enriched for cytotoxic T/NK cells and immunostimulatory dendritic cells, this cluster demonstrates the potential of m⁶A modulation to overcome immune desertification. Notably, the M1/M2 macrophage balance (P < 0.001) mirrors metabolic reprogramming linked to m⁶A-regulated pathways.

The expression and function of m⁶A modulator genes in these cells may play an important role in regulating the tumor microenvironment. Especially in immune cells, the expression of m⁶A modulator genes may affect the function and activity of immune cells, thereby regulating tumor immune responses. Most of the 23 m⁶A modulator genes are distributed in epithelial cells, B cells, and T cells. As important components of the immune system, B cells and T cells play an important role in the tumor microenvironment. The expression level and functional status of m⁶A modulator genes may affect the activity, proliferation, migration and cell fate decisions of these immune cells, thereby affecting the efficacy of tumor immunotherapy. In addition, the expression of m⁶A modulator genes in tumor cells and epithelial cells may also directly affect tumor development and treatment response. Epithelial cells are often the cells of origin of tumors, and the expression of m⁶A modulator genes in these cells may regulate the proliferation, invasion, and metastasis capabilities of tumor cells. Recent advancements such as spatial transcriptomics and proteomics, exemplified by works utilizing techniques like spatial CITE-seq (65), multimodal tri-omics (66), and spatially resolved CRISPR screens (67), offer powerful methodologies for dissecting complex interactions within the tumor microenvironment. These technologies could provide novel insights into the spatial and functional dynamics of m⁶A methylation modifications and their impact on immune infiltration and cancer progression, potentially unveiling new therapeutic avenues.

M⁶A related characteristic genes were identified as DEGs associated with the prognosis of ESCA. By employing the m⁶A signature genes, we have successfully categorized the samples into three distinct subtypes of m⁶A -related genes, which exhibit significant associations with stromal and immune activations. Therefore, it is imperative to conduct a comprehensive evaluation of m⁶A modification patterns to enhance our understanding of TME cell infiltration characterization. To mitigate inter-individual variations, quantification of the m⁶A modification pattern among m⁶A -modified tumors is necessary. To this end, we have devised a set of scoring systems, referred to as the m⁶A gene signature, to evaluate the m⁶A modification pattern. The m⁶A modification pattern associated with the immune-excluded phenotype demonstrated a higher m⁶A score, whereas the immune-inflamed phenotype exhibited a lower m⁶A score.

Our findings align with previous studies on the TME, supporting the notion that m⁶A methylation modifications play a vital role in influencing distinct immune properties within the TME. Scoring models constructed using specific biomarkers modified by m⁶A have been successfully used in gastric cancer and colorectal cancer, providing improved clinical treatment selection and prognosis assessment for cancer patients (62, 68, 69). The findings suggest that the m⁶A score possesses the capacity to serve as a comprehensive metric for assessing the m⁶A modification pattern of individual tumors, and may be employed in the investigation of tumor immunophenotype and TME immune cell infiltration. Additionally, the validation of the m⁶A score through the TCGA-ESCA cohort highlights its considerable potential as a prognostic indicator for patients afflicted with ESCA. The nomogram, incorporating the m⁶A score along with other clinical variables, demonstrated effective predictive capabilities for patient prognosis.

Furthermore, our m⁶A score demonstrated a superior predictive capability in the context of immunotherapy for esophageal cancer. These findings were robustly corroborated in the IMvigor210 cohort, as well as in cohorts receiving anti-PD-1 and anti-CTLA-4 treatments, where the immune phenotype had been established (49, 50, 56). We could also predict the efficacy of adjuvant chemotherapy and the patients’ clinical response to anti-PD-1/PD-L1 immunotherapy through the m⁶A score.

The evaluation of genes that may drive mutations in tumors is an essential method for exploring the fundamental mechanisms of tumorigenesis and progression. Furthermore, it contributes to the rational selection of cancer diagnosis and treatment strategies. In our study, we observed a significant increase in the mutation rates of TP53 and TTN in the high m⁶A score group. TP53 mutations are prevalent in various cancer types and play a critical role in cancer progression. Loss or mutation of TP53 in cancer cells can disrupt T cell recruitment and impair T cell activity, aiding immune evasion and accelerating cancer growth in the process. Research on esophageal cancer has revealed that the absence of TP53, which encodes the P53 protein, Consequently, there is an augmentation of regulatory T cell (Treg) infiltration in both paracancerous and intratumoral tissues (70). On the other hand, TTN mutations have been associated with poor immune infiltration and worse prognosis in liver hepatocellular carcinoma, colorectal cancer, and ovarian serous cystadenocarcinoma (71–73). Notably, TTN mutations are frequently detected in solid tumors and have been correlated with increased TMB. Moreover, TTN mutations have been found to be associated with the objective response to immune checkpoint blockade (ICB) therapy (74).These findings highlight the potential impact of TP53 and TTN mutations in modulating the immune response within the tumor microenvironment and their relevance to clinical outcomes. Understanding the role of these mutations in tumor biology can provide valuable insights for the development of targeted therapies and immunotherapeutic strategies in cancer treatment.

The study elucidates the role of RBMX in ESCC, focusing specifically on its impact on cell proliferation and migration. The findings suggest that the expression levels of RBMX are critical for the malignant behavior of ESCC cells. In the KYSE510 and TE-1 ESCC cell lines, significant reductions in RBMX protein expression were observed following knockdown. This indicates that RBMX may play a crucial role in maintaining the cancerous state of these cells. The plate colony formation assay revealed that RBMX knockdown significantly impaired the proliferative capacity of ESCC cells, highlighting its potential as a therapeutic target. EDU staining, which assesses DNA synthesis during the S phase, confirmed the reduced proliferative activity in cells with lower RBMX expression. These findings support the hypothesis that RBMX is a key regulator of cell cycle progression in ESCC. The wound healing assay demonstrated that RBMX knockdown significantly diminished the wound closure ability of ESCC cells, underscoring its role in cell migration, which is crucial for cancer invasion and metastasis. Migration and invasion assays further indicated significant reductions in both the invasive and migratory capabilities of ESCC cells following RBMX knockdown. These observations suggest that RBMX is central to ESCC cell motility, a key factor in the metastatic spread of cancer. In summary, this study provides evidence that RBMX has multiple influences on ESCC, impacting both cell proliferation and migration. These findings indicate that RBMX may serve as a promising target for therapeutic intervention in ESCC. Additional research is needed to elucidate the molecular mechanisms through which RBMX exerts its effects and to investigate the potential of RBMX-targeted therapies for treating ESCC.

The role of RBMX in tumors is highly tissue-specific. In hepatocellular carcinoma (HCC) and T-cell lymphoma, elevated RBMX expression enhances tumor progression and chemoresistance by stabilizing oncogenic long non-coding RNAs (lncRNAs), such as BLACAT1, or modulating RNA metabolism (75, 76). In contrast, in bladder cancer, RBMX exhibits an oncogenic effect by inhibiting hnRNP A1-mediated PKM splicing (77). This paradox indicates that the function of RBMX may rely on the tissue-specific expression of its interacting partners, such as hnRNP A1 and specific lncRNAs. RBMX has been linked to chemoresistance in both T-cell lymphoma and small-cell lung cancer (75, 78), suggesting that it may affect treatment responses in esophageal cancer, particularly in platinum-resistant ESCC, by modulating DNA damage repair and apoptotic pathways, such as those involving the BCL2 family.

Research by Tuersun and Bei has emphasized that RBMX is a significant prognostic biomarker in various cancers, including esophageal cancer, where its expression correlates with tumor progression and poor clinical outcomes (79, 80). Investigating how RBMX influences alternative splicing and m6A methylation, particularly in relation to other RNA-binding proteins such as TRA2A, may reveal new insights into the biology of esophageal cancer and resistance to therapies like sorafenib (80). The interaction of RBMX with splicing factors such as TRA2A and hnRNP A1 offers deeper insights into the regulatory networks governing esophageal cancer progression. RBMX’s role in m⁶A methylation may contribute to the dynamic regulation of oncogenic lncRNAs, thereby influencing tumor biology. Future investigations should examine the mechanistic pathways by which RBMX influences alternative splicing and m⁶A modification across a broader range of cancers. Longitudinal studies are needed to assess its prognostic value over extended periods.

Our research has several limitations that should be acknowledged. Firstly, although we included 23 well-known m⁶A regulators reported in the literature, the significance of incorporating recently identified regulators to enhance the precision of m⁶A methylation pattern identification is incontrovertible. Incorporating additional regulators into the model can potentially improve the comprehensive understanding of m⁶A modifications. Secondly, while immunotherapy has shown significant benefits for some patients with low m⁶A scores, it is important to recognize that not all patients with low scores derive equal benefit. To enhance the predictive accuracy, it would be valuable to integrate additional clinicopathological features into the analysis. By incorporating these features, we can better identify patients who are more likely to respond favorably to immunotherapy. Thirdly, although we obtained a relatively large sample size of 186 ESCA patients from various cohorts, it is important to acknowledge that a larger and independent prospective cohort of ESCA patients who have undergone immunotherapy is required to validate our findings. Prospective trials with a substantial patient cohort are required to provide a more definitive demonstration of the prognostic value of the m⁶A score in relation to the response to immunotherapy. Furthermore, our study focused on a holistic analysis of the tumor microenvironment without further distinguishing between tumor, immune, and stromal components. This lack of component-specific analysis may mask certain subtype-specific information, which is a limitation of our study. Future investigations should consider dissecting the tumor microenvironment into its individual components to gain deeper insights into the interactions and contributions of different cell types. Lastly, we primarily aimed to propose molecular subtypes associated with m⁶A methylation across the tumor microenvironment and subsequently develop a scoring system. Furthermore, clinical analysis revealed that the m⁶A score, when combined with other clinical indicators, can serve as a valuable adjunct to existing variables and effectively predict patient prognosis. Addressing these limitations through further research and validation will enhance the scientific significance and clinical applicability of our findings.

This study offers novel insights into the clinical application of immunotherapy, presenting potential implications for its use in the field. One potential avenue for the development of novel immunotherapy drugs or treatment strategies involves the modulation of m⁶A modification patterns through the targeting of m⁶A regulators or m⁶A -related marker genes. This approach may serve to reverse unfavorable immune cell infiltration in the tumor microenvironment, thereby converting immune cold tumors into hot tumors (81). These findings aid in the identification of distinct immune phenotypes, thereby enhancing our understanding of patient response to immunotherapy. This information may help with the clinical use of customized immunotherapy for the treatment of cancer (82). We also demonstrated that patients with high m⁶A scores had increased resistance to immunotherapy, which may lead to different treatment effects of classical chemotherapeutics in different patients.

5 Conclusions

We assessed the landscape of m⁶A methylation modifications mediated by 23 regulators based on 186 ESCA samples. The variety and complexity of immune infiltration in the TME are closely connected to m⁶A methylation modifications. An m⁶A score has been developed to offer a comprehensive evaluation of the m⁶A modification pattern and immune infiltration features within a singular tumor. This score also helps determine the tumor’s immune phenotype, providing new insights and directions for identifying potential therapeutic targets.

Data availability statement

All code and data related to this study are available via GitHub: https://github.com/mintsun0035/m6aTMEESCARBMX.

Ethics statement

The studies involving humans were approved by the Research Ethics Committee of Hubei University of Medicine (Shiyan, Hubei, China). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study. Ethical approval was not required for the studies on animals in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions

MS: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. PS: Conceptualization, Data curation, Formal analysis, Writing – original draft. JY: Project administration, Resources, Software, Supervision, Writing – review & editing. HZ: Conceptualization, Data curation, Formal analysis, Funding acquisition, Software, Writing – review & editing. YL: Investigation, Software, Validation, Visualization, Writing – review & editing. RC: Conceptualization, Resources, Validation, Visualization, Writing – review & editing. YF: Software, Supervision, Validation, Visualization, Writing – review & editing. LZ: Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This research was supported by the National Natural Science Foundation of China (81902498), the Natural Science Foundation of Hubei Province of China (2019CFB177, 2023AFB1008), Natural Science Foundation of Hubei Provincial Department of Education (Q20182105), Chen Xiao-Ping Foundation for the Development of Science and Technology of Hubei Provincial (CXPJJH11800001-2018333), The Foundation of Health and Family Planning Commission of Hubei Province (WJ2021Q007), Innovation and Entrepreneurship Training Program (201810929005, 201810929009, 20181092 9068, 201813249010, S201910929009, S201910929045, S202013249005, S202013249008 and 202010929009) and The Scientific and Technological Project of Taihe Hospital (2021JJXM009), and joint supported by Hubei Provincial Natural Science Foundation and Shiyan of China (2025AFD178).

Acknowledgments

We thank the Cancer Genome Atlas and the Gene Expression Omnibus for providing the data for this study.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1572810/full#supplementary-material

Supplementary Figure 1 | Integrated analysis and study design flowchart.

Supplementary Figure 2 | The Metascape enrichment network is visually represented through a visualization that highlights similarities both within and between clusters of terms. The clustering of terms is indicated by consistent color coding.

Supplementary Figure 3 | Kaplan-Meier survival curves of overall survival in ESCA cohort according to the expression value of YTHDF2, YTHDF1, RBMX, LRPPRC, IGFBP3, IGFBP1, FMR1 or ALKBH5 mRNA level in each tumor sample, the optimal value in each cohort was chosen as the cut-off point.

Supplementary Figure 4 | Comparative immune landscape analysis of m⁶A modification clusters in esophageal cancer. Heatmap depicting immune cell infiltration patterns between m⁶A cluster A (left) and cluster B (right) as quantified by four deconvolution algorithms (CIBERSORT, EPIC, MCP_COUNTER, QUANTISEQ, TIMER and XCELL). Rows represent immune cell subsets grouped by lineage (T cells, B cells, myeloid cells, stromal cells), while columns represent individual samples. Color scale indicates relative abundance (z-score normalized).

Supplementary Figure 5 | Relationship between the m⁶A score and different clinical characteristics. Kaplan-Meier curves showing the differences in survival depending on the m⁶A score and different clinical characteristics. (A) T0; (B) T1–2; (C) T 3–4; (D) N0; (E) N1–2; (F) N3; (G) M0; (H) M1; (I) male; (J) female; (K) age less than or equal to 65 years; (L) age above 65 years.

Supplementary Figure 6 | The m⁶A risk score predicts immunotherapeutic benefits in the GSE165252 cohort. The fraction of patients with clinical response to anti-PD-1 immunotherapy in low or high m⁶A risk score groups.

Supplementary Figure 7 | Sensitivity of the m⁶A risk score to different chemotherapy drugs and small molecule anticancer drugs was analyzed based on the GDSC database. (A) Axitinib, (B) Bexarotene, (C) Bicalutamide, (D) Vinblastine, (E) Bosutinib, (F) Bryostatin.1, (G) Lapatinib, (H) Imatinib, (I) Elesclomol.

Abbreviations

m⁶A, N6-methyladenosine; TME, tumor microenvironment; DCs, dendritic cells; MDSCs, myeloid-derived suppressor cells; MCs, mast cells; CAR, chimeric antigen receptor; BMDMs, bone marrow-derived macrophages; DNMTs, drugs targeting DNA methyltransferases; TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; ESCA, esophageal cancer; TPM, transcripts per kilobase million; CNV, copy number variation; ssGSEA, Single-sample gene set enrichment analysis; DEGs, differentially expressed genes; MF, molecular function; BP, biological process; CC, cellular component; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; PCA, principal component analysis; ROC, receiver operating characteristic; TMB, tumor mutation burden; TIDE, Tumor Immune Dysfunction and Exclusion; ICI, immune checkpoint inhibitor; mUC, metastatic urothelial cancer; EMT, Epithelial–mesenchymal transition.

References

1. Miao YQ, Chen W, Zhou J, Shen Q, Sun Y, Li T, et al. N(6)-adenosine-methyltransferase-14 promotes glioma tumorigenesis by repressing argininosuccinate synthase 1 expression in an m6A-dependent manner. Bioengineered. (2022) 13:1858–71. doi: 10.1080/21655979.2021.2018386

PubMed Abstract | Crossref Full Text | Google Scholar

2. Taketo K, Konno M, Asai A, Koseki J, Toratani M, Satoh T, et al. The epitranscriptome m6A writer METTL3 promotes chemo- and radioresistance in pancreatic cancer cells. Int J Oncol. (2018) 52:621–9. doi: 10.3892/ijo.2017.4219

PubMed Abstract | Crossref Full Text | Google Scholar

3. Enzinger PC and Mayer RJ. Esophageal cancer. New Engl J Med. (2003) 349:2241–52. doi: 10.1056/NEJMra035010

PubMed Abstract | Crossref Full Text | Google Scholar

4. Zhang Y. Epidemiology of esophageal cancer. World J Gastroenterol. (2013) 19:5598–606. doi: 10.3748/wjg.v19.i34.5598

PubMed Abstract | Crossref Full Text | Google Scholar

5. Oh TK, Kim K, Jheon SH, Do SH, Hwang JW, Jeon YT, et al. Long-term oncologic outcomes, opioid use, and complications after esophageal cancer surgery. J Clin Med. (2018) 7:33. doi: 10.3390/jcm7020033

PubMed Abstract | Crossref Full Text | Google Scholar

6. Xi Y, Shen Y, Wu D, Zhang J, Lin C, Wang L, et al. CircBCAR3 accelerates esophageal cancer tumorigenesis and metastasis via sponging miR-27a-3p. Mol Cancer. (2022) 21:145. doi: 10.1186/s12943-022-01615-8

PubMed Abstract | Crossref Full Text | Google Scholar

7. Altiok S, Mezzadra H, Jagannath S, Tsottles N, Rudek MA, Abdallah N, et al. A novel pharmacodynamic approach to assess and predict tumor response to the epidermal growth factor receptor inhibitor gefitinib in patients with esophageal cancer. Int J Oncol. (2010) 36:19–27.

PubMed Abstract | Google Scholar

8. Kang X, Chen K, Li Y, Li J, D’Amico TA, and Chen X. Personalized targeted therapy for esophageal squamous cell carcinoma. World J Gastroenterol. (2015) 21:7648–58. doi: 10.3748/wjg.v21.i25.7648

PubMed Abstract | Crossref Full Text | Google Scholar

9. Shah MA. Update on metastatic gastric and esophageal cancers. J Clin Oncol. (2015) 33:1760–9. doi: 10.1200/jco.2014.60.1799

PubMed Abstract | Crossref Full Text | Google Scholar

10. Doi T, Shitara K, Naito Y, Shimomura A, Fujiwara Y, Yonemori K, et al. Safety, pharmacokinetics, and antitumour activity of trastuzumab deruxtecan (DS-8201), a HER2-targeting antibody-drug conjugate, in patients with advanced breast and gastric or gastro-oesophageal tumours: a phase 1 dose-escalation study. Lancet Oncol. (2017) 18:1512–22. doi: 10.1016/s1470-2045(17)30604-6

PubMed Abstract | Crossref Full Text | Google Scholar

11. Flynn M, Young K, Cunningham D, and Starling N. The evolving immunotherapeutic landscape in advanced oesophagogastric cancer. Ther Adv Med Oncol. (2018) 10:1758835918786228. doi: 10.1177/1758835918786228

PubMed Abstract | Crossref Full Text | Google Scholar

12. Babar L, Kosovec JE, Jahangiri V, Chowdhury N, Zheng P, Omstead AN, et al. Prognostic immune markers for recurrence and survival in locally advanced esophageal adenocarcinoma. Oncotarget. (2019) 10:4546–55. doi: 10.18632/oncotarget.27052

PubMed Abstract | Crossref Full Text | Google Scholar

13. Zhu ZJ, Pang Y, Jin G, Zhang HY, Wang WH, Liu JW, et al. Hypoxia induces chemoresistance of esophageal cancer cells to cisplatin through regulating the lncRNA-EMS/miR-758-3p/WTAP axis. Aging. (2021) 13:17155–76. doi: 10.18632/aging.203062

PubMed Abstract | Crossref Full Text | Google Scholar

14. Goltsov AA, Fang B, Pandita TK, Maru DM, Swisher SG, and Hofstetter WL. HER2 confers resistance to foretinib inhibition of MET-amplified esophageal adenocarcinoma cells. Ann Thorac surgery. (2018) 105:363–70. doi: 10.1016/j.athoracsur.2017.09.003

PubMed Abstract | Crossref Full Text | Google Scholar

15. Marcone S, Buckley A, Ryan CJ, McCabe M, Lynam-Lennon N, Matallanas D, et al. Proteomic signatures of radioresistance: Alteration of inflammation, angiogenesis and metabolism-related factors in radioresistant oesophageal adenocarcinoma. Cancer Treat Res Commun. (2021) 27:100376. doi: 10.1016/j.ctarc.2021.100376

PubMed Abstract | Crossref Full Text | Google Scholar

16. Macedo-Silva C, Miranda-Gonçalves V, Lameirinhas A, Lencart J, Pereira A, Lobo J, et al. JmjC-KDMs KDM3A and KDM6B modulate radioresistance under hypoxic conditions in esophageal squamous cell carcinoma. Cell Death disease. (2020) 11:1068. doi: 10.1038/s41419-020-03279-y

PubMed Abstract | Crossref Full Text | Google Scholar

17. Fang P, Zhou J, Liang Z, Yang Y, Luan S, Xiao X, et al. Immunotherapy resistance in esophageal cancer: Possible mechanisms and clinical implications. Front Immunol. (2022) 13:975986. doi: 10.3389/fimmu.2022.975986

PubMed Abstract | Crossref Full Text | Google Scholar

18. Schiffmann LM, Plum PS, Fuchs HF, Babic B, Bruns CJ, and Schmidt T. Tumor microenvironment of esophageal cancer. Cancers (Basel). (2021) 13:4678. doi: 10.3390/cancers13184678

PubMed Abstract | Crossref Full Text | Google Scholar

19. Cui K, Hu S, Mei X, and Cheng M. Innate immune cells in the esophageal tumor microenvironment. Front Immunol. (2021) 12:654731. doi: 10.3389/fimmu.2021.654731

PubMed Abstract | Crossref Full Text | Google Scholar

20. Sun H and Sun C. The rise of NK cell checkpoints as promising therapeutic targets in cancer immunotherapy. Front Immunol. (2019) 10:2354. doi: 10.3389/fimmu.2019.02354

PubMed Abstract | Crossref Full Text | Google Scholar

21. Demaria O, Cornen S, Daëron M, Morel Y, Medzhitov R, and Vivier E. Harnessing innate immunity in cancer therapy. Nature. (2019) 574:45–56. doi: 10.1038/s41586-019-1593-5

PubMed Abstract | Crossref Full Text | Google Scholar

22. Binnewies M, Roberts EW, Kersten K, Chan V, Fearon DF, Merad M, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med. (2018) 24:541–50. doi: 10.1038/s41591-018-0014-x

PubMed Abstract | Crossref Full Text | Google Scholar

23. Galon J and Bruni D. Approaches to treat immune hot, altered and cold tumours with combination immunotherapies. Nat Rev Drug discovery. (2019) 18:197–218. doi: 10.1038/s41573-018-0007-y

PubMed Abstract | Crossref Full Text | Google Scholar

24. Chen DS and Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. (2017) 541:321–30. doi: 10.1038/nature21349

PubMed Abstract | Crossref Full Text | Google Scholar

25. Abou Khouzam R, Goutham HV, Zaarour RF, Chamseddine AN, Francis A, Buart S, et al. Integrating tumor hypoxic stress in novel and more adaptable strategies for cancer immunotherapy. Semin Cancer Biol. (2020) 65:140–54. doi: 10.1016/j.semcancer.2020.01.003

PubMed Abstract | Crossref Full Text | Google Scholar

26. Li M, Zha X, and Wang S. The role of N6-methyladenosine mRNA in the tumor microenvironment. Biochim Biophys Acta Rev cancer. (2021) 1875:188522. doi: 10.1016/j.bbcan.2021.188522

PubMed Abstract | Crossref Full Text | Google Scholar

27. He Y, Xing J, Wang S, Xin S, Han Y, and Zhang J. Increased m6A methylation level is associated with the progression of human abdominal aortic aneurysm. Ann Trans Med. (2019) 7:797. doi: 10.21037/atm.2019.12.65

PubMed Abstract | Crossref Full Text | Google Scholar

28. Gu X, Zhang Y, Li D, Cai H, Cai L, and Xu Q. N6-methyladenosine demethylase FTO promotes M1 and M2 macrophage activation. Cell signalling. (2020) 69:109553. doi: 10.1016/j.cellsig.2020.109553

PubMed Abstract | Crossref Full Text | Google Scholar

29. Liu Y, Liu Z, Tang H, Shen Y, Gong Z, Xie N, et al. The N(6)-methyladenosine (m(6)A)-forming enzyme METTL3 facilitates M1 macrophage polarization through the methylation of STAT1 mRNA. Am J Physiol Cell Physiol. (2019) 317:C762–c775. doi: 10.1152/ajpcell.00212.2019

PubMed Abstract | Crossref Full Text | Google Scholar

30. Li H, Chiappinelli KB, Guzzetta AA, Easwaran H, Yen RW, Vatapalli R, et al. Immune regulation by low doses of the DNA methyltransferase inhibitor 5-azacitidine in common human epithelial cancers. Oncotarget. (2014) 5:587–98. doi: 10.18632/oncotarget.1782

PubMed Abstract | Crossref Full Text | Google Scholar

31. Wrangle J, Wang W, Koch A, Easwaran H, Mohammad HP, Vendetti F, et al. Alterations of immune response of Non-Small Cell Lung Cancer with Azacytidine. Oncotarget. (2013) 4:2067–79. doi: 10.18632/oncotarget.1542

PubMed Abstract | Crossref Full Text | Google Scholar

32. Chong W, Shang L, Liu J, Fang Z, Du F, Wu H, et al. m(6)A regulator-based methylation modification patterns characterized by distinct tumor microenvironment immune profiles in colon cancer. Theranostics. (2021) 11:2201–17. doi: 10.7150/thno.52717

PubMed Abstract | Crossref Full Text | Google Scholar

33. Zhang B, Wu Q, Li B, Wang D, Wang L, and Zhou Y. mA regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. J Mol Cancer. (2020) 19:53. doi: 10.1186/s12943-020-01170-0

PubMed Abstract | Crossref Full Text | Google Scholar

34. Li X, Ma S, Deng Y, Yi P, and Yu J. Targeting the RNA mA modification for cancer immunotherapy. J Mol Cancer. (2022) 21:76. doi: 10.1186/s12943-022-01558-0

PubMed Abstract | Crossref Full Text | Google Scholar

35. Chong W, Shang L, Liu J, Fang Z, Du F, Wu H, et al. mA regulator-based methylation modification patterns characterized by distinct tumor microenvironment immune profiles in colon cancer. J Theranostics. (2021) 11:2201–17. doi: 10.7150/thno.52717

PubMed Abstract | Crossref Full Text | Google Scholar

36. Şenbabaoğlu Y, Michailidis G, and Li J. Critical limitations of consensus clustering in class discovery. J Sci Rep. (2014) 4:6207. doi: 10.1038/srep06207

PubMed Abstract | Crossref Full Text | Google Scholar

37. Wilkerson M and Hayes D. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. J Bioinf. (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | Crossref Full Text | Google Scholar

38. Yu Z, Wong H, and Wang H. Graph-based consensus clustering for class discovery from gene expression data. J Bioinf. (2007) 23:2888–96. doi: 10.1093/bioinformatics/btm463

PubMed Abstract | Crossref Full Text | Google Scholar

39. Shi K, Li Y, Yang L, Zhang Z, Guo D, Zhang J, et al. Profiling transcriptional heterogeneity of epithelium, fibroblasts, and immune cells in esophageal squamous cell carcinoma by single-cell RNA sequencing. FASEB J. (2022) 36:e22620. doi: 10.1096/fj.202200898R

PubMed Abstract | Crossref Full Text | Google Scholar

40. Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. (2021) 184:3573–87. doi: 10.1016/j.cell.2021.04.048

PubMed Abstract | Crossref Full Text | Google Scholar

41. Wang X, Park J, Susztak K, Zhang N, and Li M. Bulk tissue cell type deconvolution with multi-subject single-cell expression reference. J Nat Commun. (2019) 10:380. doi: 10.1038/s41467-018-08023-x

PubMed Abstract | Crossref Full Text | Google Scholar

42. Sturm G, Finotello F, Petitprez F, Zhang JD, Baumbach J, Fridman WH, et al. Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinf (Oxford England). (2019) 35:i436–45. doi: 10.1093/bioinformatics/btz363

PubMed Abstract | Crossref Full Text | Google Scholar

43. Huang L, Zhu J, Kong W, Li P, and Zhu S. Expression and prognostic characteristics of m6A RNA methylation regulators in colon cancer. Int J Mol Sci. (2021) 22:2134. doi: 10.3390/ijms22042134

PubMed Abstract | Crossref Full Text | Google Scholar

44. Zhu M, Cui Y, Mo Q, Zhang J, Zhao T, Xu Y, et al. Characterization of m(6)A RNA methylation regulators predicts survival and immunotherapy in lung adenocarcinoma. Front Immunol. (2021) 12:782551. doi: 10.3389/fimmu.2021.782551

PubMed Abstract | Crossref Full Text | Google Scholar

45. Zeng D, Li M, Zhou R, Zhang J, Sun H, Shi M, et al. Tumor microenvironment characterization in gastric cancer identifies prognostic and immunotherapeutically relevant gene signatures. Cancer Immunol Res. (2019) 7:737–50. doi: 10.1158/2326-6066.Cir-18-0436

PubMed Abstract | Crossref Full Text | Google Scholar

46. Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J, et al. Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst. (2006) 98:262–72. doi: 10.1093/jnci/djj052

PubMed Abstract | Crossref Full Text | Google Scholar

47. Wu Z, Li S, and Zhu X. The mechanism of stimulating and mobilizing the immune system enhancing the anti-tumor immunity. Front Immunol. (2021) 12:682435. doi: 10.3389/fimmu.2021.682435

PubMed Abstract | Crossref Full Text | Google Scholar

48. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. J Nat Med. (2018) 24:1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | Crossref Full Text | Google Scholar

49. Riaz N, Havel JJ, Makarov V, Desrichard A, Urba WJ, Sims JS, et al. Tumor and microenvironment evolution during immunotherapy with nivolumab. Cell. (2017) 171:934–949 e16. doi: 10.1016/j.cell.2017.09.028

PubMed Abstract | Crossref Full Text | Google Scholar

50. Van Allen EM, Miao D, Schilling B, Shukla SA, Blank C, Zimmer L, et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science. (2015) 350:207–11. doi: 10.1126/science.aad0095

PubMed Abstract | Crossref Full Text | Google Scholar

51. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. (2018) 554:544–8. doi: 10.1038/nature25501

PubMed Abstract | Crossref Full Text | Google Scholar

52. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. (2013) 41:D955–61. doi: 10.1093/nar/gks1111

PubMed Abstract | Crossref Full Text | Google Scholar

53. Geeleher P, Cox N, and Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One. (2014) 9:e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | Crossref Full Text | Google Scholar

54. Lu X, Li R, Ying Y, Zhang W, and Wang W. Gene signatures, immune infiltration, and drug sensitivity based on a comprehensive analysis of m6a RNA methylation regulators in cervical cancer. J Trans Med. (2022) 20:385. doi: 10.1186/s12967-022-03600-7

PubMed Abstract | Crossref Full Text | Google Scholar

55. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. (2019) 16:1289–96. doi: 10.1038/s41592-019-0619-0

PubMed Abstract | Crossref Full Text | Google Scholar

56. Balar AV, Galsky MD, Rosenberg JE, Powles T, Petrylak DP, Bellmunt J, et al. Atezolizumab as first-line treatment in cisplatin-ineligible patients with locally advanced and metastatic urothelial carcinoma: a single-arm, multicentre, phase 2 trial. Lancet. (2017) 389:67–76. doi: 10.1016/s0140-6736(16)32455-2

PubMed Abstract | Crossref Full Text | Google Scholar

57. Chen H, Yang M, Wang Q, Song F, Li X, and Chen K. The new identified biomarkers determine sensitivity to immune check-point blockade therapies in melanoma. Oncoimmunology. (2019) 8:1608132. doi: 10.1080/2162402X.2019.1608132

PubMed Abstract | Crossref Full Text | Google Scholar

58. Turley SJ, Cremasco V, and Astarita JL. Immunological hallmarks of stromal cells in the tumour microenvironment. Nat Rev Immunol. (2015) 15:669–82. doi: 10.1038/nri3902

PubMed Abstract | Crossref Full Text | Google Scholar

59. Gajewski TF, Woo SR, Zha Y, Spaapen R, Zheng Y, Corrales L, et al. Cancer immunotherapy strategies based on overcoming barriers within the tumor microenvironment. Curr Opin Immunol. (2013) 25:268–76. doi: 10.1016/j.coi.2013.02.009

PubMed Abstract | Crossref Full Text | Google Scholar

60. Joyce JA and Fearon DT. T cell exclusion, immune privilege, and the tumor microenvironment. Science. (2015) 348:74–80. doi: 10.1126/science.aaa6204

PubMed Abstract | Crossref Full Text | Google Scholar

61. Salmon H, Franciszkiewicz K, Damotte D, Dieu-Nosjean MC, Validire P, Trautmann A, et al. Matrix architecture defines the preferential localization and migration of T cells into the stroma of human lung tumors. J Clin Invest. (2012) 122:899–910. doi: 10.1172/jci45817

PubMed Abstract | Crossref Full Text | Google Scholar

62. Chen H, Yao J, Bao R, Dong Y, Zhang T, Du Y, et al. Cross-talk of four types of RNA modification writers defines tumor microenvironment and pharmacogenomic landscape in colorectal cancer. Mol Cancer. (2021) 20:29. doi: 10.1186/s12943-021-01322-w

PubMed Abstract | Crossref Full Text | Google Scholar

63. Hu H, Piotrowska Z, Hare PJ, Chen H, Mulvey HE, Mayfield A, et al. Three subtypes of lung cancer fibroblasts define distinct therapeutic paradigms. Cancer Cell. (2021) 39:1531–1547 e10. doi: 10.1016/j.ccell.2021.09.003

PubMed Abstract | Crossref Full Text | Google Scholar

64. Du Y, Zhan J, Chen S, Liu Q, Song J, Li J, et al. IDDF2024-ABS-0311 RBMX promotes tumor plasticity and immune suppression by activating IL-33/TGF-β signaling niche in colorectal cancer. Gut. (2024) 73:A193–8. doi: 10.1136/gutjnl-2024-IDDF.123

Crossref Full Text | Google Scholar

65. Liu Y, DiStasio M, Su G, Asashima H, Enninful A, Qin X, et al. High-plex protein and whole transcriptome co-mapping at cellular resolution with spatial CITE-seq. Nat Biotechnol. (2023) 41:1405–9. doi: 10.1038/s41587-023-01676-0

PubMed Abstract | Crossref Full Text | Google Scholar

66. Fan R, Zhang D, Rodriguez-Kirby L, Lin Y, Song M, Wang L, et al. Spatial dynamics of mammalian brain development and neuroinflammation by multimodal tri-omics mapping. Res square. (2024). doi: 10.21203/rs.3.rs-4814866/v1

PubMed Abstract | Crossref Full Text | Google Scholar

67. Baysoy A, Tian X, Zhang F, Renauer P, Bai Z, Shi H, et al. Spatially Resolved in vivo CRISPR Screen Sequencing via Perturb-DBiT. bioRxiv: preprint server Biol. (2024). doi: 10.1101/2024.11.18.624106

PubMed Abstract | Crossref Full Text | Google Scholar

68. Zhang B, Wu Q, Li B, Wang D, Wang L, and Zhou YL. m(6)A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. (2020) 19:53. doi: 10.1186/s12943-020-01170-0

PubMed Abstract | Crossref Full Text | Google Scholar

69. Sun M, Xie M, Zhang T, Wang Y, Huang W, and Xia L. m(6)A methylation modification patterns and tumor microenvironment infiltration characterization in pancreatic cancer. Front Immunol. (2021) 12:739768. doi: 10.3389/fimmu.2021.739768

PubMed Abstract | Crossref Full Text | Google Scholar

70. Mariathasan S, Turley S, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. J Nat. (2018) 554:544–8. doi: 10.1038/nature25501

PubMed Abstract | Crossref Full Text | Google Scholar

71. Wang R, Hu X, Liu X, Bai L, Gu J, and Li Q. Construction of liver hepatocellular carcinoma-specific lncRNA-miRNA-mRNA network based on bioinformatics analysis. PLoS One. (2021) 16:e0249881. doi: 10.1371/journal.pone.0249881

PubMed Abstract | Crossref Full Text | Google Scholar

72. Gomes FC and Figueiredo ERL. Social, genetics and histopathological factors related to titin (TTN) gene mutation and survival in women with ovarian serous cystadenocarcinoma: bioinformatics analysis. Genes (Basel). (2023) 14:1092. doi: 10.3390/genes14051092

PubMed Abstract | Crossref Full Text | Google Scholar

73. Zhao L and Fan W. Construction of a TTN mutation-based prognostic model for evaluating immune microenvironment. Cancer Stemness Outcomes Colorectal Cancer Patients. (2023) 2023:6079957. doi: 10.1155/2023/6079957

PubMed Abstract | Crossref Full Text | Google Scholar

74. Jia Q, Wang J, He N, He J, and Zhu B. Titin mutation associated with responsiveness to checkpoint blockades in solid tumors. JCI Insight. (2019) 4:e127901. doi: 10.1172/jci.insight.127901

PubMed Abstract | Crossref Full Text | Google Scholar

75. Schumann FL, Bauer M, Gross E, Terziev D, Wienke A, Wickenhauser C, et al. RBMX protein expression in T-cell lymphomas predicts chemotherapy response and prognosis. Cancers (Basel). (2021) 13:4788. doi: 10.3390/cancers13194788

PubMed Abstract | Crossref Full Text | Google Scholar

76. Song Y, He S, Ma X, Zhang M, Zhuang J, Wang G, et al. RBMX contributes to hepatocellular carcinoma progression and sorafenib resistance by specifically binding and stabilizing BLACAT1. Am J Cancer Res. (2020) 10:3644–65.

PubMed Abstract | Google Scholar

77. Yan Q, Zeng P, Zhou X, Zhao X, Chen R, Qiao J, et al. RBMX suppresses tumorigenicity and progression of bladder cancer by interacting with the hnRNP A1 protein to regulate PKM alternative splicing. Oncogene. (2021) 40:2635–50. doi: 10.1038/s41388-021-01666-z

PubMed Abstract | Crossref Full Text | Google Scholar

78. Zhang Z, Zhang C, Yang Z, Zhang G, Wu P, Luo Y, et al. m(6)A regulators as predictive biomarkers for chemotherapy benefit and potential therapeutic targets for overcoming chemotherapy resistance in small-cell lung cancer. J Hematol Oncol. (2021) 14:190. doi: 10.1186/s13045-021-01173-4

PubMed Abstract | Crossref Full Text | Google Scholar

79. Tuersun H, Liu L, Zhang J, Maimaitizunong R, Tang X, and Li H. m6A reading protein RBMX as a biomarker for prognosis and tumor progression in esophageal cancer. Transl Cancer Res. (2023) 12:2319–35. doi: 10.21037/tcr-23-84

PubMed Abstract | Crossref Full Text | Google Scholar

80. Bei M, Hao S, Lin K, Chen Q, Cai Y, Zhao X, et al. Splicing factor TRA2A contributes to esophageal cancer progression via a noncanonical role in lncRNA m(6) A methylation. Cancer Sci. (2023) 114:3216–29. doi: 10.1111/cas.15870

PubMed Abstract | Crossref Full Text | Google Scholar

81. Zhao B, Roundtree I, and He C. Post-transcriptional gene regulation by mRNA modifications. J Nat Rev Mol Cell Biol. (2017) 18:31–42. doi: 10.1038/nrm.2016.132

PubMed Abstract | Crossref Full Text | Google Scholar

82. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. J Nat Commun. (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | Crossref Full Text | Google Scholar