Comprehensive Analysis of PD-L1 Expression, Immune Infiltrates, and m6A RNA Methylation Regulators in Esophageal Squamous Cell Carcinoma

Background: Esophageal squamous cell carcinoma (ESCC) is one of the most common cancer types and represents a threat to global public health. N6-Methyladenosine (m6A) methylation plays a key role in the occurrence and development of many tumors, but there are still few studies investigating ESCC. This study attempts to construct a prognostic signature of ESCC based on m6A RNA methylation regulators and to explore the potential association of these regulators with the tumor immune microenvironment (TIME).

Methods: The transcriptome sequencing data and clinical information of 20 m6A RNA methylation regulators in 453 patients with ESCC (The Cancer Genome Atlas [TCGA] cohort, n = 95; Gene Expression Omnibus [GEO] cohort, n = 358) were obtained. The differing expression levels of m6A regulators between ESCC and normal tissue were evaluated. Based on the expression of these regulators, consensus clustering was performed to investigate different ESCC clusters. PD-L1 expression, immune score, immune cell infiltration and potential mechanisms among different clusters were examined. LASSO Cox regression analysis was utilized to obtain a prognostic signature based on m6A RNA methylation modulators. The relationship between the risk score based on the prognostic signature and the TIME of ESCC patients was studied in detail.

Results: Six m6A regulators (METTL3, WTAP, IGF2BP3, YTHDF1, HNRNPA2B1 and HNRNPC) were observed to be significantly highly expressed in ESCC tissues. Two molecular subtypes (clusters 1/2) were determined by consensus clustering of 20 m6A modulators. The expression level of PD-L1 in ESCC tissues increased significantly and was significantly negatively correlated with the expression levels of YTHDF2, METL14 and KIAA1429. The immune score, CD8 T cells, resting mast cells, and regulatory T cells (Tregs) in cluster 2 were significantly increased. Gene set enrichment analysis (GSEA) shows that this cluster involves multiple hallmark pathways. We constructed a five-gene prognostic signature based on m6A RNA methylation, and the risk score based on the prognostic signature was determined to be an independent prognostic indicator of ESCC. More importantly, the prognostic value of the prognostic signature was verified using another independent cohort. m6A regulators are related to TIME, and their copy-number alterations will dynamically affect the number of tumor-infiltrating immune cells.

Conclusion: Our study established a strong prognostic signature based on m6A RNA methylation regulators; this signature was able to accurately predict the prognosis of ESCC patients. The m6A methylation regulator may be a key mediator of PD-L1 expression and immune cell infiltration and may strongly affect the TIME of ESCC.

Introduction

Esophageal cancer (EC) is a highly aggressive malignancy and is the eighth most common malignant tumor in the world; furthermore, the mortality rate of EC ranks sixth worldwide and is still rising (1–3). As with the National Central Cancer Registry of China (NCCR) statistics, Chinese EC patients compose up to 70% of all EC cases worldwide (4). Among the two main histopathological subtypes of EC, esophageal squamous cell carcinoma (ESCC) is more common than esophageal adenocarcinoma (5). Neoadjuvant chemotherapy or chemoradiotherapy have been utilized as the standard treatments for ESCC, improving the surgical effect (6). Although related immunotherapies involving esophageal cancer are still in the preliminary stages of research, some related immunosuppressants have entered clinical trials and have exhibited long-lasting antitumor activity and controllable adverse reactions. Studies have found that PD-L1 positive patients with advanced and metastatic esophageal cancer are very likely to be sensitive to immunotherapy (7). Abnormal levels of N6-methyladenosine (m6A) methylation also play vital roles in the progression of various cancers (8). Accurately predicting the prognosis of ESCC is the key to the success of clinical diagnosis and treatment and individualized medication. Therefore, the identification of novel and reliable prognostic molecular signatures from multiple dimensions is very important for selecting the most appropriate treatment strategy and improving the poor prognosis of ESCC patients.

Numerous studies in recent years have shown that N6-methyladenosine (m6A) methylation is a commonly seen modification in eukaryotic messenger RNA (mRNA) and strongly affects many basic biological processes, such as cell differentiation, tissue development and tumorigenesis (9–12). The level of m6A methylation is regulated by methyltransferases, demethylases, and binding proteins. m6A RNA modification can be catalysed enzymatically by various methyltransferases, known as m6A “writers” (METTL3, METTL14, METTL16, KIAA1429, WTAP, RBM15, RBM15B and ZC3H13). N6-methyladenosine in RNA can be removed by demethylases, known as m6A “erasers” (FTO and ALKBH5). Proteins that selectively bind m6A can be defined as m6A “readers” (YTHDC1, YTHDC2, IGF2BP2, IGF2BP3, YTHDF1, YTHDF2, YTHDF3, HNRNPA2B1, HNRNPC, and RBMX) that exert regulatory functions by selective recognition of methylated RNA (13). In addition, emerging evidence suggests that m6A modulators have cancer promoter or inhibitory effects in the development of various malignant tumors. Zhang et al. demonstrated the m6A modulator-mediated methylation modification pattern and the tumor microenvironment infiltration characteristics of gastric cancer (14–16). Han et al. (17) confirmed that m6A methylation can prolong the neoantigen-specific immunity mediated by YTHDF1. YTHDF1 may be a potential therapeutic target and an important mediator of tumor immune evasion. These findings suggest that both m6A methylation regulators and the tumor immune microenvironment (TIME) may affect the prognosis of cancer patients. The TIME may affect the patient’s response to immune checkpoint inhibitors, so PD-L1 expression in ESCC should be further considered to evaluate tumor immunity. However, to the best of our knowledge, there are few studies on the use of m6A methylation regulators to predict the prognosis of ESCC. In addition, the correlation between m6A methylation regulators and PD-L1, the expression of PD-L1 in ESCC, and the abundance of immune infiltrating cells need to be fully studied.

In this study, RNA sequencing data, clinical information and immune cell data of ESCC patients were obtained from The Cancer Genome Atlas (TCGA). According to the expression of m6A RNA methylation regulators, consensus cluster analysis was performed, and a gene signature and risk model were constructed to predict the prognosis of ESCC patients more accurately. Although there have been reports of consensus cluster analysis on patients with esophageal cancer (18), the study did not separate the subgroups of esophageal squamous cell carcinoma and further verification, and the number of m6A RNA methylation regulators included in the study was not large enough. To enhance the prediction performance of the m6A methylation-related gene signature, we verified it in another ESCC cohort in the Gene Expression Omnibus (GEO) database. In addition, we also studied the infiltration of immune cells in the TIME, the expression of PD-L1 in ESCC, and the correlation with m6A RNA methylation regulators.

Materials and Methods

Data Acquisition

Clinical information, including sex, T stage, N stage, M stage, TNM stage, survival information and RNA-seq expression profiles from ESCC (95 patients) cohorts of TCGA database (https://portal.gdc.cancer.gov), was used in this study. We also selected gene expression profiles of patient-derived ESCC tissue (GSE53625, GPL18109, n = 358) from the GEO database to validate the candidate prognostic gene signature identified from TCGA data. This data set contains gene sequencing information and clinical information of cancer tissues and adjacent normal tissues from 358 ESCC patients in China.

m6A RNA Methylation Regulator Collection

According to previously published research reports, we collected 20 m6A RNA methylation regulators (METTL3, METTL14, METTL16, KIAA1429, WTAP, RBM15, RBM15B, ZC3H13, YTHDC1, YTHDC2, IGF2BP2, IGF2BP3, YTHDF1, YTHDF2, YTHDF3, HNRNPA2B1, HNRNPC, RBMX, FTO and ALKBH5) for further research.

Bioinformatic Analysis

The Search Tool for the Retrieval of Interacting Genes/Proteins database (STRING, version 11.0, http://string-db.org/) has powerful functions in studying gene interactions and visualization (19). We entered 20 m6A RNA methylation regulators into STRING to understand their interactions and performed functional annotation analysis to initially explore the biological process (BP), cell composition (CC), molecular function (MF) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways involved in these regulators. Pearson correlation analysis was used to elucidate the correlation between different m6A RNA methylation regulators.

To functionally explore the biological properties of m6A regulators in ESCC, we used the “ConsensusClusterPlus” package (http://www.bioconductor.org/, 1000 iterations and resampling rate of 80%) to divide ESCC patients from the TCGA database into different groups. Gene set enrichment analysis (GSEA) was utilized to understand the biological processes involved in different subgroups. Hallmarks in GSEA were used to identify predefined gene sets; 5000 permutations were performed according to the gene set to determine p-values. A pathway with a p-value < 0.01 and a false discovery rate (FDR) < 0.25 was considered to be significant, as described in the Results section.

CIBERSORT (http://cibersort.stanford.edu/), a deconvolution algorithm based on gene expression, can calculate the composition of immune cells from the gene expression profile of complex tissues (20). We used CIBERSORT software to calculate the infiltration level of 22 immune cells based on the expression profile data of ESCC in the TCGA database. Subsequently, the ESTIMATE algorithm (“estimate” package in R) was used to calculate the immune score of each patient, and the difference in the immune score between the two cluster subgroups was evaluated (21).

The least absolute shrinkage and selection operator (LASSO) Cox regression model includes all m6A methylation regulators to construct a strong prognostic signature and calculate the coefficient of each gene. The coefficients obtained from the LASSO regression algorithm were used to yield the following risk score equation: risk score = sum of coefficients × m6A regulator expression level. According to the score, the gene signature with the strongest ability to predict the prognosis of ESCC patients was obtained. The ESCC patients were divided into high- and low-risk groups according to the median risk score, and the difference in overall survival (OS) between the risk score groups was evaluated. Through the receiver operating characteristic (ROC) curve and the area under the ROC curve (AUC value), the accuracy of the genetic signature for predicting prognosis was evaluated. Univariate and multivariate Cox regression analyses were used to identify independent prognostic factors for ESCC patients. Importantly, we verified this prognostic signature using a cohort of 358 ESCC samples from the GEO database.

Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/) was used to analyze the six main immune cells in tumors to evaluate the effect of copy number alternations (CNAs) of m6A regulators on the level of immune cell infiltration.

Statistical Analysis

R (version 3.6.3; The R Foundation for Statistical Computing) and SPSS 23.0 (IBM Corp., New York, USA) were used for data analysis and statistics. m6A regulator expression and clinical data were analysed by χ² test or Fisher’s exact test. Student’s t-test and one-way ANOVA were used to separately perform the group comparisons of two subgroups and more than two subgroups. Kaplan-Meier survival analysis and log-rank tests were used to analyse the differences in OS between different risk score groups. The correlation of gene expression was evaluated by Spearman’s R and statistical significance. An absolute value of R greater than 0.25 was considered relevant, and a p-value < 0.05 was considered to be significant.

Results

Expression of m6A RNA Methylation Regulators in ESCC

To understand the expression pattern of m6A RNA methylation regulators between ESCC tumors and normal tissues, we drew a heatmap based on 20 m6A methylation genes in the TCGA database. The red or green in the figure indicates relatively high or low expression, respectively (Figure 1A). The expression levels of writers (METTL3 and WTAP) and readers (i.e., IGF2BP3, YTHDF1, HNRNPA2B1 and HNRNPC) were notably higher in ESCC tissues than in normal adjacent tissues (Figure 1B) (P < 0.05). The abnormal expression of these m6A RNA methylation regulators may indicate an important biological role in the occurrence and development of ESCC.

FIGURE 1

Figure 1 The expression levels of m⁶A RNA methylation regulators between tumor and normal samples in TCGA ESCC cohort. (A) Heatmap of m⁶A RNA methylation regulator expression level in each sample; (B) the expression difference of m6A RNA methylation regulator between tumor and normal samples. * means p<0.05; ** means p<0.01; **** means p<0.0001.

Correlation and Functional Enrichment Between m6A RNA Methylation Regulators

With the STRING database used to further understand the interaction between the 20 m6A RNA methylation regulators, a PPI network was obtained. After deleting the isolated genes without interaction, we determined that the PPI network contained 20 nodes and 112 edges, as shown in Figure 2A. KIAA1429 and METTL3 appeared to be the hub genes of the interaction network. As shown in Table 1, we performed Gene Ontology (GO) analysis on these 20 RNA methylation regulators to obtain a preliminary understanding of their biological functions. Regulation of mRNA metabolic process, N6-methyladenosine-containing RNA binding and RNA N6-methyladenosine methyltransferase complex were the most significantly enriched GO items. In addition, we observed that all m6A RNA methylation regulators were generally positively correlated, and KIAA1429 had the highest correlation with YTHDF3 (Figure 2B) (r = 0.74).

FIGURE 2

Figure 2 Interaction and correlation between m⁶A RNA methylation regulators in ESCC. (A) A PPI network was constructed to evaluate the interaction between m⁶A RNA methylation regulators; (B) the correlations among m⁶A RNA methylation regulators were analyzed by Pearson correlation.

TABLE 1

Table 1 Functional annotation of 20 m6A methylation regulators.

Consensus Clustering Identified Two Clusters of Patients With ESCC

A consensus cluster consisting of 20 m6A RNA methylation regulators was constructed by using the “ConsensusClusterPlus” package. Figures 3A, B show the relative change of the cumulative distribution function (CDF) of the consensus cluster from k = 2 to 9 and the area under the CDF curve from k = 2 to 9, and k = 2 is proven to be the most suitable choice to divide the ESCC patient cohort into two clusters (Figure 3C). The tracking plot for k = 2 to k = 10 is demonstrated in Figure 3D. In addition, we studied the relationship between clustering subgroups and clinicopathological parameters of ESCC patients. The results showed that there was no significance between cluster 1 and cluster 2 and patient age, sex, T stage, N stage, M stage, or TNM stage (Figure 3E).

FIGURE 3

Figure 3 Consensus clustering identified two ESCC patient clusters and their relationship with clinicopathological parameters. (A) Consensus clustering cumulative distribution function (CDF) for k=2 to 9; (B) relative change in area under CDF curve for k=2 to 9; (C) The ESCC cohort from TCGA was divided into two distinct clusters when k=2; (D) distribution of each sample when k ranges from 2 to 9; (E) comparison of the relationship between the clinicopathological characteristics of two clusters.

Association of PD-L1 With m6A RNA Methylation

We assessed the difference in PD-L1 expression between tumor and normal tissues in ESCC patients. Compared with normal adjacent tissues, the expression of PD-L1 in ESCC tissues was significantly increased (Figure 4A) (P < 0.001). In the cluster subtypes of cluster 1 and cluster 2 that we constructed, the expression difference of PD-L1 was not significant (Figure 4B). In addition, the expression of PD-L1 in ESCC patients was significantly negatively correlated with the expression levels of YTHDF2, METL14 and KIAA1429 (Figure 4C). The ratio of 22 immune cell types between the two subgroups was analysed (Figure 5).

FIGURE 4

Figure 4 Association of PD-L1 with m6A RNA methylation and the landscape of immune cell infiltration in ESCC. (A) PD-L1 expression was significantly higher in ESCC cohort from TCGA; (B) the expression level of PD-L1 in cluster1/2 subtypes; (C) the correlation of PD-L1 with m⁶A methylation regulators in ESCC cohort from TCGA. **** means p<0.0001.

FIGURE 5

Figure 5 The infiltrating levels of 22 immune cell types in cluster1/2 in ESCC. * means p<0.05; ** means p<0.01.

Immune Cell Infiltration in Consensus Cluster Subgroups of Esophageal Squamous Cell Carcinoma

Subsequently, we evaluated the immune scores of the ESCC immune microenvironment based on the ESTIMATE algorithm and found that the immune scores of the two m6A RNA methylation-related clusters were significantly different (Figure 6A). Our research shows that cluster 2 exhibits high levels of CD8 T cells, resting mast cells and regulatory T cell (Treg) infiltration (Figures 6B–D) (P < 0.05). To explore the underlying regulatory mechanism that led to the TIME difference between the two subgroups, we performed GSEA. The results showed that the first 5 hallmark pathways that were significantly enriched in cluster 1 included bile acid metabolism, fatty acid metabolism, myogenesis, oxidative phosphorylation, and reactive oxygen species pathway, while the hallmark pathways involved in cluster 2 included E2F targets, the epithelial-mesenchymal transition, the G2/M checkpoint, the mitotic spindle and TNFA signalling via NFKB (Figures 6E, F).

FIGURE 6

Figure 6 Differences in the level of immune cell infiltration between the two clusters in ESCC cohort and the biological pathways involved. (A) immune score in the cluster 1/2 subtypes; (B–D) the infiltrating levels of the CD8 T cells (B), resting mast cells (C), and regulatory T cell (D) in two clusters; (E, F) GSEA shows the first 5 signaling pathways involved in cluster 1 (E) and cluster 2 (F). * means p<0.05; ** means p<0.01.

Construction of Prognostic Signatures Based on m6A RNA Methylation Genes

To explore the prognostic value of these 20 m6A RNA methylation regulators in ESCC, we performed univariate Cox regression. The results demonstrated that METTL16 (P = 0.005), KIAA1429 (P = 0.002), RBM15 (P = 0.006), IGF2BP3 (P = 0.02), YTHDF1 (P = 0.041), YTHDF3 (P = 0.005) and ALKBH (P = 0.043) were significantly correlated with OS (Figure 7A). The hazard ratio (HR) of these genes with prognostic value was less than 1. Subsequently, the LASSO algorithm was used to obtain the coefficient of each prognostic gene (Figures 7B, C). According to the minimum standard, 5 m6A regulators (HNRNPC, RBM15, IGF2BP3, METTL16 and KIAA1429) were selected to construct a prognostic signature, and the risk score of each ESCC patient was calculated. The formula was as follows: risk score = (0.0933 * HNRNPC expression) − (0.2370 * RBM15 expression) − (0.2949 * IGF2BP3 expression) − (0.5176 * METTL16 expression) − (0.6240 * KIAA1429 expression). According to the median risk score, ESCC patients were divided into low-risk and high-risk groups.

FIGURE 7

Figure 7 Construction of the prognostic signature based on TCGA ESCC cohort. (A) Univariate analysis of 20 m6A RNA methylation regulators to identify the genes that significantly correlated with OS; (B, C) The prognostic signature constructed by the minimum criterion of LASSO Cox regression algorithm; (D) The Kaplan-Meier curve shows that the risk score based on the prognostic signature of m6A RNA methylation is significantly correlated with OS in ESCC patients; (E) time-dependent ROC curves was applied to assess the predictive efficiency of the signature in TCGA.

Risk Score Based on the Prognostic Signature Is an Independent Prognostic Factor in the ESCC Cohort of TCGA and GEO

To verify the prognostic value of risk grouping in ESCC patients, Kaplan-Meier survival analysis was performed. The results showed that the OS of patients in the high-risk group was significantly lower than that in the low-risk group (Figure 7D) (P < 0.0001). Subsequently, ROC curves were drawn to evaluate the specificity and sensitivity of the prognostic signatures associated with m6A RNA methylation regulators. The results showed that the areas under the curve (AUCs) at 1, 3 and 5 years were 0.73, 0.73 and 0.78, respectively, indicating good prediction performance (Figure 7E). Importantly, we used the GSE53625 dataset from the GEO database to verify the prognostic value of the prognostic signature in ESCC. The results also showed that the high-risk score group was significantly related to worse prognosis in patients (Figure 8A). The 1-, 3-, and 4-year area under the curve (AUC) values of the ROC curve were 0.6, 0.6, and 0.69, respectively, exhibiting a good distinguishing performance for ESCC patients (Figure 8B).

FIGURE 8

Figure 8 Verification of the prognostic value of m6A methylation regulator-related signatures in the GEO database. (A) The high risk score of the signature in the GEO database indicates poor OS; (B) time-dependent ROC curve in the GEO database confirms the predictive efficiency of the prognostic signature.

Moreover, we performed univariate and multivariate Cox regression analyses in the two cohorts of TCGA and GEO to determine whether the risk score based on prognostic markers is an independent prognostic indicator for ESCC patients. In the TCGA cohort, after univariate analysis obtained factors related to OS in ESCC patients, multivariate Cox regression analysis showed that risk score (P<0.001, HR = 6.665), N stage (P = 0.018, HR = 2.123), TNM stage (P = 0.035, HR = 1.797), and sex (P = 0.029, HR = 9.628) were identified as independent prognostic factors (Figure 9A). The risk scores, OS and OS status distributions of 95 ESCC patients from the TCGA database are shown in Figure 9C. The prognostic value of the risk score calculated with the prognostic signatures of 5 m6A RNA methylation modulators was verified in 358 ESCC patients from the GEO database. The univariate analysis showed that the risk score (P < 0.001, HR = 2.718), N stage (P < 0.001, HR = 1.438) and TNM stage (P < 0.001, HR = 1.994) were significantly correlated with OS, and subsequent multivariate Cox regression analysis demonstrated that the risk score (P < 0.001, HR = 2.769) and TNM stage (P = 0.021, HR = 2.013) were independent prognostic factors for ESCC patients (Figure 9B). The risk scores, OS and OS status distributions of ESCC patients verified in the GEO database are shown in Figure 9D.

FIGURE 9

Figure 9 Evaluation of prognostic values of the risk scores. Univariate and multivariate Cox regression analysis of the risk scores in TCGA. (A) and GEO (B) database. The distributions of risk scores and OS status in TCGA (C) and GEO (D) database.

Prognostic Risk Scores Correlated With Stage, Immune Score, and Clinicopathological Parameters in ESCC

The relationship between risk score and clinical characteristics and cluster subgroups was further evaluated. It can be seen from the heatmap that the high immune score generally corresponds to the high expression of the 5 regulators. HNRNPC and IGF2BP3 are highly expressed in cluster 1. There were significant differences between the high-risk group and the low-risk group in clustering subtypes (P <0.001) and age (P <0.05) of ESCC (Figure 10A). We also further examined the relationship between risk score and subtype, immune score and TNM staging. The results show that the risk score of cluster 2 is significantly higher than that of cluster 1 (P <0.001, Figure 10B). Although the high and low immune scores and risk scores are not significant, it can be seen that the median risk score of the high immune score group is higher than that of the low immune score group (Figure 10C). In addition, TNM staging and risk score were not found to have significant statistical significance (Figure 10D). These findings indicate that the risk score of ESCC patients may profoundly influence clinical outcomes.

FIGURE 10

Figure 10 Prognostic risk scores correlated with stage, immune score, and clinicopathological parameters in ESCC. (A) Heatmap and clinicopathologic parameters of high and low risk groups; The relationship between risk score and cluster 1/2 (B), immune score (C) and TNM staging (D).

Effect of Genetic Alterations of the m6A Regulator Signatures on Immune Cell Infiltration

The relationship between the risk scores of nine immune cell types and the level of infiltration was analyzed to evaluate the impact of five m6A regulator-based signatures on the ESCC immune microenvironment. Perhaps due to the limitation of sample size, we only found a significant negative correlation between the risk score and the infiltration level of macrophages (P =0.039) and neutrophils (P =0.019, Figure 11). Risk signatures based on m6A regulators are likely to have a potential role in regulating the immune microenvironment of ESCC.

FIGURE 11

Figure 11 Relationships between the risk score and infiltration abundances of nine immune cell types. (A) B cell; (B) CD4+ T cell; (C) activated memory CD4+ T cell; (D) CD8+ T cell; (E) follicular helper T cell; (F) macrophage; (G) myeloid dendritic cell; (H) neutrophil; (I) regulatory T cell.

In order to preliminarily clarify the potential mechanism of risk score and different immune cell infiltration, the influence of somatic cell CNAs based on m6A regulator on immune cell infiltration was further analyzed. The identified CNAs of m6A regulator signatures, including arm-level deletion, high amplication, and arm-level gain, significantly affected the infiltration level of B cells, CD4 + T cells, CD8 + T cells, neutrophils and dendritic cells in ESCC (Figure 12). This series of studies further show that m6A methylation regulators have a key regulatory effect on TIME in ESCC patients.

FIGURE 12

Figure 12 Effect of the genetic alterations of m6A regulator-relevant signature on the immune cell infiltration. (A) RBM15; (B) KIAA1429; (C) HNRNPC; (D) IGF2BP3. * means p<0.05; ** means p<0.01; *** means p<0.001.

Discussion

ESCC is a highly malignant tumor. Surgery combined with radiotherapy, chemotherapy, and targeted therapy are currently the main treatments for ESCC. However, due to the high rate of recurrence and metastasis, the survival results for this cancer are far from satisfactory (22–24). Although related immunotherapies involving esophageal cancer are still in the preliminary research stage, some related inhibitors have entered clinical trials and have shown long-lasting antitumor activity and controllable adverse reactions, which indicates that the immune microenvironment of ESCC warrants further exploration (7). At the same time, m6A methylation, which is the most common form of mRNA modification, has been proven to promote or suppress cancer in many tumor types (25), but there are very few related studies in ESCC. Therefore, it is necessary to further explore the role played by m6A methylation in ESCC and the infiltration of the ESCC TIME. In addition, the effect of m6A methylation on the ESCC TIME has not been fully elucidated.

In this study, the expression patterns, prognostic values and effects on the TIME of m6A RNA methylation regulators in ESCC were explored. The expression levels of m6A “writers” (METTL3 and WTAP) and “readers” (i.e., IGF2BP3, YTHDF1, HNRNPA2B1 and HNRNPC) in ESCC were significantly higher than those in normal tissues, and elevated expression levels of 2 genes (METTL3 and IGF2BP3) have been reported to be an independent prognostic factor for ESCC patients (26, 27); other genes have not been studied. In STRING, a protein-protein interaction (PPI) network composed of 20 m6A RNA methylation regulators was generated, and the biological functions involved in the regulators were preliminarily analysed through GO functional annotation. In addition, KIAA1429 and YTHDF3, which have the highest correlation coefficients, have not been reported in tumors, and they are likely to play important roles in the occurrence and development of ESCC.

Next, two molecular subtypes (clusters 1/2) were determined by consensus clustering of 20 m6A methylation regulators. Perhaps due to insufficient sample size, there was no significant difference between cluster 1/2 subtypes and clinicopathological parameters of ESCC patients. We also determined that the expression of PD-L1 in ESCC tumor tissues was significantly higher than that in normal tissues. This phenomenon is not obvious in the cluster 1/2 subtype. In addition, it has been observed that PD-L1 is significantly negatively correlated with YTHDF2, METL14, and KIAA1429. These regulators are likely to predict the efficacy of immunotherapy in ESCC patients, which requires further research. We also studied the immune cell infiltration of the TIME in the cluster 1/2 subtype based on m6A regulators in ESCC patients. The results showed that the infiltration level of CD8 T cells, resting mast cells and regulatory T cells (Tregs) in cluster 2 was higher than that in cluster 1. Studies have shown that high mast cell density is related to the progression of ESCC and reduced postoperative survival. The high concentration of Tregs in ESCC can lead to immune escape and promote tumor progression (28, 29). The immune score calculated based on the ESTIMATE algorithm is also significantly higher in cluster 2, which indicates a significant difference in the TIME of ESCC patients. GSEA of cluster 1 was performed to explore the biological processes involved, and the results showed that it involved multiple carcinogenic pathways of digestive tract tumors (30–33). Similarly, we determined the signaling pathway involved in the regulator in cluster 2, and the relationship between it and the level of immune cell infiltration in the tumor microenvironment warrants further exploration.

Furthermore, we also constructed a five-gene prognostic signature consisting of HNRNPC, RBM15, IGF2BP3, METTL16, and KIAA14297 from m6A methylation regulators, and the calculated risk score showed good performance in predicting the prognostic outcome of ESCC patients. In addition to the risk score, N stage, TNM stage and sex were also independent prognostic factors. Importantly, the potential of m6A methylation modulator-related prognostic signatures was verified in another ESCC cohort from the GEO database. The GSE53625 verification cohort of 358 patients also shows the strong prognostic potential of this signature. A high risk score is significantly related to poor OS in ESCC patients and is an independent prognostic factor for ESCC patients. Among these risk markers, HNRNPC can interact with LBX2-AS1 and enhance the stability of ZEB1 and ZEB2 mRNA, thereby promoting the migration and development of ESCC (34). Another study showed that the high expression of HNRNPC is significantly related to the poor OS of patients with lung adenocarcinoma, and it is likely to be an oncogene in patients with lung adenocarcinoma and breast cancer (35, 36). IGF2BP3 is highly expressed in many tumors, including ESCC, lung adenocarcinoma, colon cancer, and gastric cancer, and leads to poor prognosis (26, 37–39). There is no research report on KIAA1429 in ESCC, but some scholars pointed out that KIAA1429 can promote the progression of liver cancer, breast cancer and osteosarcoma, leading to poor prognosis (40–42). For RBM15 and METL16, there are currently few reports in tumors, but they may be potential prognostic biomarkers.

At present, the effect of m6A methylation regulators on immune cell infiltration in TIME is still unclear. In this study, the risk score based on the prognostic signatures of five m6A regulators was significantly negatively correlated with the infiltration level of macrophages and neutrophils. High immune scores often correspond to high expression of HNRNPC and IGF2BP3. Studies have shown that HNRNPC and HNRNPK in the subfamily of heterogeneous ribonucleoproteins (hnRNPs) can regulate the recruitment and activation of neutrophils and macrophages respectively, thereby affecting the systemic immune response (43, 44). The research on the other three risk signatures and immune cell infiltration remains unclear. For other m6A regulators, Han et al. showed that in mouse tumors lacking YTHDF1, the level of CD8 + T and NK cell infiltration increased, thereby enhancing the cross-expression of tumor antigens in vivo and the cross-priming of CD8 + T cells (17). Li et al. reported that the loss of METT3 or METT14 triggers the disorder of T cell proliferation and differentiation, thereby reducing the sensitivity of interleukin 7 (IL-7) in vivo (45). Our study also showed that CNAs of m6A methylation regulators, including arm level deletion, high amplication and arm level gain, significantly affect the level of immune cell infiltration in ESCC. The m6A methylation regulator is likely to have a crucial regulatory effect on TIME in ESCC patients.

This study has several limitations. First, the results may be affected by the small sample size. It is necessary to further improve the sample size, sequencing data and clinical information of ESCC patients in future research. In addition, our conclusions are based on the results of bioinformatic analysis of datasets containing genetic and other molecular information from patient tissues, which need to be further verified in clinical studies.

In summary, this study systematically evaluated the expression of m6A RNA regulators in ESCC, their correlation with PD-L1, and potential regulatory mechanisms. Two ESCC subtypes (clusters 1/2) were obtained through consensus clustering of m6A regulators, and the difference in the level of immune cell infiltration in the TIME was determined. m6A RNA regulators may improve the responsiveness of ESCC patients to immunotherapy by regulating the TIME and expression of PD-L1. More importantly, we constructed a prognostic signature containing 5 genes based on m6A RNA methylation, and the risk score was determined to be an independent prognostic factor in two different ESCC cohorts, indicating that the prognostic signature is a promising tool for predicting the survival outcome of ESCC patients.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding authors.

Author Contributions

JH, SG, WG, and FT conceived, designed, and supervised the study. WG, FT, and QH drafted the manuscript. WG, FT, QH, ZW, FS, GZ, ZY, RL, and QX collected the data. QH, ZW, FS, GZ, ZY, RL, and QX performed all data analysis. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by National Key R&D Program (2018YFC1312105), the Institutional Fundamental Research Funds (2018PT32033), the Ministry of Education Innovation Team Development Project (IRT-17R10), and the Beijing Hope Run Special Fund of Cancer Foundation of China (LC2019B15).

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.

Acknowledgments

We sincerely thank the data provided by TCGA and GEO databases, and thank all the staff of thoracic surgery for their support during the study.

References

1. Gupta B, Kumar N. Worldwide Incidence, Mortality and Time Trends for Cancer of the Oesophagus. Eur J Cancer Prev (2017) 26(2):107–18. doi: 10.1097/CEJ.0000000000000249

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Siegel RL, Miller KD, Jemal A. Cancer Statistics, 2018. CA Cancer J Clin (2018) 68(1):7–30. doi: 10.3322/caac.21442

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Smyth EC, Lagergren J, Fitzgerald RC, Lordick F, Shah MA, Lagergren P, et al. Oesophageal Cancer. Nat Rev Dis Primers (2017) 3:17048. doi: 10.1038/nrdp.2017.48

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Chen W, Zheng R, Baade PD, Zhang S, Zeng H, Bray F, et al. Cancer Statistics in China, 2015. CA Cancer J Clin (2016) 66(2):115–32. doi: 10.3322/caac.21338

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Abnet CC, Arnold M, Wei WQ. Epidemiology of Esophageal Squamous Cell Carcinoma. Gastroenterology (2018) 154(2):360–73. doi: 10.1053/j.gastro.2017.08.023

PubMed Abstract | CrossRef Full Text | Google Scholar

6. de Gouw D, Klarenbeek BR, Driessen M, Bouwense SAW, van Workum F, Fütterer JJ, et al. Detecting Pathological Complete Response in Esophageal Cancer After Neoadjuvant Therapy Based on Imaging Techniques: A Diagnostic Systematic Review and Meta-Analysis. J Thoracic Oncol Off Publ Int Assoc Study Lung Cancer (2019) 14(7):1156–71. doi: 10.1016/j.jtho.2019.04.004

CrossRef Full Text | Google Scholar

7. Shah MA, Kojima T, Hochhauser D, Enzinger P, Raimbourg J, Hollebecque A, et al. Efficacy and Safety of Pembrolizumab for Heavily Pretreated Patients With Advanced, Metastatic Adenocarcinoma or Squamous Cell Carcinoma of the Esophagus: The Phase 2 KEYNOTE-180 Study. JAMA Oncol (2019) 5(4):546–50. doi: 10.1001/jamaoncol.2018.5441

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Chen M, Wong CM. The Emerging Roles of N6-Methyladenosine (m6A) Deregulation in Liver Carcinogenesis. Mol Cancer (2020) 19(1):44. doi: 10.1186/s12943-020-01172-y

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the Human and Mouse M6a RNA Methylomes Revealed by M6a-Seq. Nature (2012) 485(7397):201–6. doi: 10.1038/nature11112

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Li Z, Weng H, Su R, Weng X, Zuo Z, Li C, et al. FTO Plays an Oncogenic Role in Acute Myeloid Leukemia as a N(6)-Methyladenosine RNA Demethylase. Cancer Cell (2017) 31(1):127–41. doi: 10.1016/j.ccell.2016.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Su R, Dong L, Li C, Nachtergaele S, Wunderlich M, Qing Y, et al. R-2hg Exhibits Anti-Tumor Activity by Targeting FTO/m(6)A/MYC/CEBPA Signaling. Cell (2018) 172(1-2):90–105.e23. doi: 10.1016/j.cell.2017.11.031

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Yang L, Liu X, Song L, Su G, Di A, Bai C, et al. Melatonin Restores the Pluripotency of Long-Term-Cultured Embryonic Stem Cells Through Melatonin Receptor-Dependent m6A RNA Regulation. J Pineal Res (2020) 69(2):e12669. doi: 10.1111/jpi.12669

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Luo J, Liu H, Luan S, He C, Li Z. Aberrant Regulation of mRNA m⁶A Modification in Cancer Development. Int J Mol Sci (2018) 19(9):2515. doi: 10.3390/ijms19092515

CrossRef Full Text | Google Scholar

14. Li H, Su Q, Li B, Lan L, Wang C, Li W, et al. High Expression of WTAP Leads to Poor Prognosis of Gastric Cancer by Influencing Tumour-Associated T Lymphocyte Infiltration. J Cell Mol Med (2020) 24(8):4452–65. doi: 10.1111/jcmm.15104

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Yang DD, Chen ZH, Yu K, Lu JH, Wu QN, Wang Y, et al. METTL3 Promotes the Progression of Gastric Cancer Via Targeting the MYC Pathway. Front Oncol (2020) 10:115. doi: 10.3389/fonc.2020.00115

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m(6)A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer. Mol Cancer (2020) 19(1):53. doi: 10.1186/s12943-020-01170-0

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Han D, Liu J, Chen C, Dong L, Liu Y, Chang R, et al. Anti-Tumour Immunity Controlled Through mRNA m(6)A Methylation and YTHDF1 in Dendritic Cells. Nature (2019) 566(7743):270–4. doi: 10.1038/s41586-019-0916-x

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Xu LC, Pan JX, Pan HD. Construction and Validation of an M6a RNA Methylation Regulators-Based Prognostic Signature for Esophageal Cancer. Cancer Manage Res (2020) 12:5385–94. doi: 10.2147/CMAR.S254870

CrossRef Full Text | Google Scholar

19. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. String v11: Protein-Protein Association Networks With Increased Coverage, Supporting Functional Discovery in Genome-Wide Experimental Datasets. Nucleic Acids Res (2019) 47(D1):D607–d13. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

21. 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. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Cheng YF, Chen HS, Wu SC, Chen HC, Hung WH, Lin CH, et al. Esophageal Squamous Cell Carcinoma and Prognosis in Taiwan. Cancer Med (2018) 7(9):4193–201. doi: 10.1002/cam4.1499

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Hirano H, Kato K. Systemic Treatment of Advanced Esophageal Squamous Cell Carcinoma: Chemotherapy, Molecular-Targeting Therapy and Immunotherapy. Japanese J Clin Oncol (2019) 49(5):412–20. doi: 10.1093/jjco/hyz034

CrossRef Full Text | Google Scholar

24. Hu X, Wu D, He X, Zhao H, He Z, Lin J, et al. circGSK3β Promotes Metastasis in Esophageal Squamous Cell Carcinoma by Augmenting β-Catenin Signaling. Mol Cancer (2019) 18(1):160. doi: 10.1186/s12943-019-1095-y

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Su Y, Huang J, Hu J. m(6)A RNA Methylation Regulators Contribute to Malignant Progression and Have Clinical Prognostic Impact in Gastric Cancer. Front Oncol (2019) 9:1038. doi: 10.3389/fonc.2019.01038

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Wakita A, Motoyama S, Sato Y, Nagaki Y, Fujita H, Terata K, et al. IGF2BP3 Expression Correlates With Poor Prognosis in Esophageal Squamous Cell Carcinoma. J Surg Res (2020) 259:137–44. doi: 10.1016/j.jss.2020.10.024

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Xia TL, Yan SM, Yuan L, Zeng MS. Upregulation of METTL3 Expression Predicts Poor Prognosis in Patients With Esophageal Squamous Cell Carcinoma. Cancer Manage Res (2020) 12:5729–37. doi: 10.2147/CMAR.S245019

CrossRef Full Text | Google Scholar

28. Fakhrjou A, Niroumand-Oscoei SM, Somi MH, Ghojazadeh M, Naghashi S, Samankan S. Prognostic Value of Tumor-Infiltrating Mast Cells in Outcome of Patients With Esophagus Squamous Cell Carcinoma. J Gastrointest Cancer (2014) 45(1):48–53. doi: 10.1007/s12029-013-9550-2

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zheng Y, Chen Z, Han Y, Han L, Zou X, Zhou B, et al. Immune Suppressive Landscape in the Human Esophageal Squamous Cell Carcinoma Microenvironment. Nat Commun (2020) 11(1):6268. doi: 10.1038/s41467-020-20019-0

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Ashton TM, McKenna WG, Kunz-Schughart LA, Higgins GS. Oxidative Phosphorylation as an Emerging Target in Cancer Therapy. Clin Cancer Res an Off J Am Assoc Cancer Res (2018) 24(11):2482–90. doi: 10.1158/1078-0432.CCR-17-3070

CrossRef Full Text | Google Scholar

31. Currie E, Schulze A, Zechner R, Walther TC, Farese RV Jr. Cellular Fatty Acid Metabolism and Cancer. Cell Metab (2013) 18(2):153–61. doi: 10.1016/j.cmet.2013.05.017

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Li T, Apte U. Bile Acid Metabolism and Signaling in Cholestasis, Inflammation, and Cancer. Adv Pharmacol (San Diego Calif) (2015) 74:263–302. doi: 10.1016/bs.apha.2015.04.003

CrossRef Full Text | Google Scholar

33. Prasad S, Gupta SC, Tyagi AK. Reactive Oxygen Species (ROS) and Cancer: Role of Antioxidative Nutraceuticals. Cancer Lett (2017) 387:95–105. doi: 10.1016/j.canlet.2016.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhang Y, Chen W, Pan T, Wang H, Zhang Y, Li C. LBX2-AS1 is Activated by ZEB1 and Promotes the Development of Esophageal Squamous Cell Carcinoma by Interacting With HNRNPC To Enhance the Stability of ZEB1 and ZEB2 mRNAs. Biochem Biophys Res Commun (2019) 511(3):566–72. doi: 10.1016/j.bbrc.2019.02.079

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Guo W, Huai Q, Zhang G, Guo L, Song P, Xue X, et al. Elevated Heterogeneous Nuclear Ribonucleoprotein C Expression Correlates With Poor Prognosis in Patients With Surgically Resected Lung Adenocarcinoma. Front Oncol (2021) 10:3145. doi: 10.3389/fonc.2020.598437

CrossRef Full Text | Google Scholar

36. Wu Y, Zhao W, Liu Y, Tan X, Li X, Zou Q, et al. Function of HNRNPC in Breast Cancer Cells by Controlling the dsRNA-induced Interferon Response. EMBO J (2018) 37(23):e99017. doi: 10.15252/embj.201899017

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Guo W, Huai Q, Wan H, Guo L, Song P, Gao S, et al. Prognostic Impact of IGF2BP3 Expression in Patients With Surgically Resected Lung Adenocarcinoma. DNA Cell Biol (2021) 40(2):316–31. doi: 10.1089/dna.2020.6136

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Xu W, Sheng Y, Guo Y, Huang Z, Huang Y, Wen D, et al. Increased IGF2BP3 Expression Promotes the Aggressive Phenotypes of Colorectal Cancer Cells In Vitro and Vivo. J Cell Physiol (2019) 234(10):18466–79. doi: 10.1002/jcp.28483

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zhou Y, Huang T, Siu HL, Wong CC, Dong Y, Wu F, et al. IGF2BP3 Functions as a Potential Oncogene and is a Crucial Target of miR-34a in Gastric Carcinogenesis. Mol Cancer (2017) 16(1):77. doi: 10.1186/s12943-017-0647-2

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Han Q, Yang J, Yang H, Li C, Li J, Cao Y. KIAA1429 Promotes Osteosarcoma Progression by Promoting Stem Cell Properties and is Regulated by miR-143-3p. Cell Cycle (Georgetown Tex) (2020) 19(10):1172–85. doi: 10.1080/15384101.2020.1749465

CrossRef Full Text | Google Scholar

41. Lan T, Li H, Zhang D, Xu L, Liu H, Hao X, et al. KIAA1429 Contributes to Liver Cancer Progression Through N6-methyladenosine-dependent Post-Transcriptional Modification of GATA3. Mol Cancer (2019) 18(1):186. doi: 10.1186/s12943-019-1106-z

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Qian JY, Gao J, Sun X, Cao MD, Shi L, Xia TS, et al. KIAA1429 Acts as an Oncogenic Factor in Breast Cancer by Regulating CDK1 in an N6-methyladenosine-independent Manner. Oncogene (2019) 38(33):6123–41. doi: 10.1038/s41388-019-0861-z

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Bhandary YP, Velusamy T, Shetty P, Shetty RS, Idell S, Cines DB, et al. Post-Transcriptional Regulation of Urokinase-Type Plasminogen Activator Receptor Expression in Lipopolysaccharide-Induced Acute Lung Injury. Am J Respir Crit Care Med (2009) 179(4):288–98. doi: 10.1164/rccm.200712-1787OC

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Ostareck DH, Ostareck-Lederer A. RNA-Binding Proteins in the Control of LPS-Induced Macrophage Response. Front Genet (2019) 10:31. doi: 10.3389/fgene.2019.00031

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Li HB, Tong J, Zhu S, Batista PJ, Duffy EE, Zhao J, et al. m(6)A mRNA Methylation Controls T Cell Homeostasis by Targeting the IL-7/STAT5/SOCS Pathways. Nature (2017) 548(7667):338–42. doi: 10.1038/nature23450

PubMed Abstract | CrossRef Full Text | Google Scholar