m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer
Molecular Cancer volume 19, Article number: 53 (2020)
The epigenetic regulation of immune response has been demonstrated in recent studies. Nonetheless, potential roles of RNA N6-methyladenosine (m6A) modification in tumor microenvironment (TME) cell infiltration remain unknown.
We comprehensively evaluated the m6A modification patterns of 1938 gastric cancer samples based on 21 m6A regulators, and systematically correlated these modification patterns with TME cell-infiltrating characteristics. The m6Ascore was constructed to quantify m6A modification patterns of individual tumors using principal component analysis algorithms.
Three distinct m6A modification patterns were determined. The TME cell-infiltrating characteristics under these three patterns were highly consistent with the three immune phenotypes of tumors including immune-excluded, immune-inflamed and immune-desert phenotypes. We demonstrated the evaluation of m6A modification patterns within individual tumors could predict stages of tumor inflammation, subtypes, TME stromal activity, genetic variation, and patient prognosis. Low m6Ascore, characterized by increased mutation burden and activation of immunity, indicated an inflamed TME phenotype, with 69.4% 5-year survival. Activation of stroma and lack of effective immune infiltration were observed in the high m6Ascore subtype, indicating a non-inflamed and immune-exclusion TME phenotype, with poorer survival. Low m6Ascore was also linked to increased neoantigen load and enhanced response to anti-PD-1/L1 immunotherapy. Two immunotherapy cohorts confirmed patients with lower m6Ascore demonstrated significant therapeutic advantages and clinical benefits.
This work revealed the m6A modification played a nonnegligible role in formation of TME diversity and complexity. Evaluating the m6A modification pattern of individual tumor will contribute to enhancing our cognition of TME infiltration characterization and guiding more effective immunotherapy strategies.
In all living organisms, as the third layer of epigenetics, more than 150 RNA modifications including 5-methylcytosine (m5C), N6-methyladenosine (m6A) and N1-methyladenosine (m1A) have been identified [1, 2]. Among these modifications, m6A RNA methylation, which are widely found in the mRNA, lncRNA as well as miRNA, is recognized as the most prominent and abundant form of internal modifications in eukaryotic cells, of whose abundance account for 0.1–0.4% total adenosine residues [3,4,5]. Similar to the modification of DNA and protein, m6A modification is a kind of dynamic reversible process in mammalian cells, which is regulated by methyltransferases, demethylases and binding proteins, also known as “writers”, “erasers” and “readers” . The formation process of m6A methylation is catalyzed by methyltransferases consisting of RBM15, ZC3H13, METTL3, METTL14, WTAP and KIAA1429, while the removal process is mediated by demethylases including FTO and ALKBH5. In addition, a group of specific RNA-binding proteins composed of YTHDF1/2/3, YTHDC1/2, HNRNPA2B1, LRPPRC, FMR1 and so on can recognize m6A motif, thus affecting m6A functions [7, 8]. The in-depth understanding of these regulators would help reveal the role and mechanism of m6A modification in post-transcriptional regulation. It has been reported that the m6A regulators play a crucial role in a variety of biological functions in vivo [9,10,11]. Increasing evidence demonstrated that dysregulated expression and genetic changes of m6A regulators were correlated with the disorders of multiple biological process including dysregulate cell death and proliferation, developmental defects, tumor malignant progression, impaired self-renewal capacity, and immunomodulatory abnormality [12,13,14].
Immunotherapy represented by immunological checkpoint blockade (ICB, PD-1/L1 and CTLA-4) has demonstrated astounding clinical efficacy in a small percentage of patients with durable responses. Unfortunately, the majority of patients experienced minimal or no clinical benefit, far from a met clinical need . Traditionally, the tumor progression has been considered as a multi-step process that only involves the genetic and epigenetic variation in tumor cells. However, numerous studies have shown that the microenvironment in which tumor cells depend for growth and survival also play a crucial role in the tumor progression. The tumor part was composed of a complex tumor microenvironment (TME) that not only contained cancer cells but also stromal cells such as resident fibroblasts (cancer associated fibroblast; CAF) and macrophages, and distant recruited cells such as infiltrating immune cells (myeloid cells and lymphocytes), bone marrow-derived cells (BMDCs) such as endothelial progenitor and hematopoietic progenitor cells, secreted factors such as cytokines, chemokines, growth factors, and new blood vessels. Of these, five distinct myeloid populations including tumor-associated macrophages (TAM), tumor-associated neutrophils (TANs), dendritic cells, myeloid-derived suppressor cells (MDSCs) and Tie2-expressing monocytes comprised the tumor-associated myeloid cells (TAMCs) . Cancers cells elicited multiple biological behavior changes through direct and indirect interactions with other TME components such as inducing proliferation and angiogenesis, inhibiting apoptosis, avoiding hypoxia as well as inducing immune tolerance. As the understanding of the diversity and complexity of tumor microenvironment has deepened, emerging evidence reveals its critical role in the tumor progression, immune escape, and its effect on response to immunotherapy. Predicting the response to ICB based on the characterization of TME cell infiltration is a key procedure on increasing the success of existing ICBs and exploiting novel immunotherapeutic strategies [17, 18]. Therefore, by comprehensively parsing the TME landscape heterogeneity and complexity, different tumor immune phenotypes are likely to be identified, and the abilities of guiding and predicting immunotherapeutic responsiveness would also improve. Additionally, the promising biomarkers could be revealed, which will prove highly effective in recognizing patients’ response to immunotherapy and will benefit the search for new therapeutic targets [19, 20].
Recently, several studies have revealed the special correlation between TME infiltrating immune cells and m6A modification, which can’t be explained via RNA degradation mechanism. Dali et al. reported that binding of YTHDF1 to the transcripts encoding lysosomal proteases modified by m6A methylation improved the translational efficiency of lysosomal cathepsins in dendritic cells (DCs), while suppression of cathepsins in DC significantly strengthened its ability to cross-present tumor antigens, which in turn enhanced tumor infiltrating CD8+ T cell antitumor response. And YTHDF1 inhibition also improved the therapeutic efficacy of anti-PD-L1 blockade . The study of Huamin et al. revealed that METTL3-mediated m6A modification promoted the activation and maturation of DCs. Declining expression of co-stimulatory molecules CD80 and CD40 resulted by METTL3 specific depletion reduced the ability of stimulating T cell activation. And down-regulation of Tirap inhibited the transmission of the TLR4/NF-κB signaling pathway and decreased the secretion of pro-inflammatory cytokines . In addition, some studies have focused on the intrinsic oncogenic pathways induced by dysregulated expression and genomic variation of m6A regulators. For example, Qiang et al. found that METTL3 overexpression promote gastric cancer (GC) malignant progression and liver metastasis through angiogenesis and glycolysis pathway .
However, the above studies have necessarily been confined to only one or two m6A regulators and cell types owing to technical limitations, while the antitumor effect is characterized by numerous tumor suppressor factors that interact in a highly coordinated manner. Therefore, comprehensive recognizing the TME cell infiltration characterizations mediated by multiple m6A regulators will contribute to enhancing our understanding of TME immune regulation. In this study, we integrated the genomic information of 1938 gastric cancer samples to comprehensively evaluate the m6A modification patterns, and correlated the m6A modification pattern with the TME cell-infiltrating characteristics. We revealed three distinct m6A modification patterns, and surprisingly found that the TME characteristics under these three patterns were highly consistent with the immune-excluded phenotype, immune-inflamed phenotype and immune-desert phenotype, respectively, suggesting the m6A modification played a nonnegligible role in shaping individual tumor microenvironment characterizations. For that, we established a set of scoring system to quantify the m6A modification pattern in individual patients.
Gastric cancer dataset source and preprocessing
The workflow of our study was shown in Figure S1A. Public gene-expression data and full clinical annotation were searched in Gene-Expression Omnibus (GEO) and the Cancer Genome Atlas (TCGA) database. Patients without survival information were removed from further evaluation. In total, 7 eligible GC cohorts (GSE15459, GSE34942, GSE57303, GSE62254/ACRG, GSE84437, GSE26253 and TCGA-STAD (The Cancer Genome Atlas-Stomach Adenocarcinoma)) were gathered in this study for further analysis. For microarray data from Affymetrix®, we downloaded the raw “CEL” files and adopted a robust multiarray averaging method with the affy and simpleaffy packages to perform background adjustment and quantile normalization. For microarray data from other platforms, the normalized matrix files were directly downloaded. As to datasets in TCGA, RNA sequencing data (FPKM value) of gene expression were downloaded from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/) using the R package TCGAbiolinks , which was specifically developed for integrative analysis with GDC data . Then FPKM values were transformed into transcripts per kilobase million (TPM) values. Batch effects from non-biological technical biases were corrected using the “ComBat” algorithm of sva package. The baseline information of all eligible GC datasets was summarized in Table S1. The somatic mutation data was acquired from TCGA database. The GSE62717 dataset from ACRG cohort was downloaded for Copy Number Variation (CNV) analysis. Data were analyzed with the R (version 3.6.1) and R Bioconductor packages.
Unsupervised clustering for 21 m6A regulators
Owing to the few m6A regulators detected by Illumina HumanRef-8 WG-DASL v3.0 platform, we did not include GSE26253 cohort for clustering analysis. A total of 21 regulators were extracted from five integrated GEO datasets for identifying different m6A modification patterns mediated by m6A regulators. These 21 m6A regulators included 8 writers (METTL3, METTL14, RBM15, RBM15B, WTAP, KIAA1429, CBLL1, ZC3H13), 2 erasers (ALKBH5, FTO) and 11 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, HNRNPA2B1, HNRNPC, FMR1, LRPPRC, ELAVL1). Unsupervised clustering analysis was applied to identify distinct m6A modification patterns based on the expression of 21 m6A regulators and classify patients for further analysis. The number of clusters and their stability were determined by the consensus clustering algorithm . We used the ConsensuClusterPlus package to perform the above steps and 1000 times repetitions were conducted for guaranteeing the stability of classification .
Gene set variation analysis (GSVA) and functional annotation
To investigate the difference on biological process between m6A modification patterns, we performed GSVA enrichment analysis using “GSVA” R packages. GSVA, in a non-parametric and unsupervised method, is commonly employed for estimating the variation in pathway and biological process activity in the samples of an expression dataset . The gene sets of “c2.cp.kegg.v6.2.symbols” were downloaded from MSigDB database for running GSVA analysis. Adjusted P with value less than 0.05 was considered as statistically significance. The clusterProfiler R package was used to perform functional annotation for m6A-related genes, with the cutoff value of FDR < 0.05.
Estimation of TME cell infiltration
We used the ssGSEA (single-sample gene-set enrichment analysis) algorithm to quantify the relative abundance of each cell infiltration in the GC TME. The gene set for marking each TME infiltration immune cell type was obtained from the study of Charoentong, which stored various human immune cell subtypes including activated CD8 T cell, activated dendritic cell, macrophage, natural killer T cell, regulatory T cell and so on (Table S2) [28, 29]. The enrichment scores calculated by ssGSEA analysis were utilized to represent the relative abundance of each TME infiltrating cell in each sample.
Identification of differentially expressed genes (DEGs) between m6A distinct phenotypes
To identify m6A-related genes, we classified patients into three distinct m6A modification patterns based on the expression of 21 m6A regulators. The empirical Bayesian approach of limma R package was applied to determine DEGs between different modification patterns . The significance criteria for determining DEGs was set as adjusted P value < 0.001.
Generation of m6A gene signature
To quantify the m6A modification patterns of individual tumor, we constructed a set of scoring system to evaluate the m6A modification pattern of individual patients with gastric cancer—the m6A gene signature, and we termed as m6Ascore. The procedures for establishment of m6A gene signature were as follows:
The DEGs identified from different m6Aclusters were firstly normalized among all ACRG samples and the overlap genes were extracted. The patients were classified into several groups for deeper analysis by adopting unsupervised clustering method for analyzing overlap DEGs. The consensus clustering algorithm was utilized for defining the number of gene clusters as well as their stability. Then, we performed the prognostic analysis for each gene in the signature using univariate Cox regression model. The genes with the significant prognosis were extracted for further analysis. We then conducted principal component analysis (PCA) to construct m6A relevant gene signature. Both principal component 1 and 2 were selected to act as signature scores. This method had advantage of focusing the score on the set with the largest block of well correlated (or anticorrelated) genes in the set, while down-weighting contributions from genes that do not track with other set members. We then define the m6Ascore using a method similar to GGI [31, 32]:
where i is the expression of m6A phenotype-related genes.
Correlation between m6A gene signature and other related biological processes
Mariathasan et al. constructed a set of gene sets that stored genes associated with some biological processes, including (1) immune-checkpoint; (2) antigen processing machinery; (3) CD8 T-effector signature; (4) epithelial-mesenchymal transition (EMT) markers including EMT1, EMT2 and EMT3; (5) Angiogenesis signature; (7) pan-fibroblast TGFb response signature (Pan-F-TBRS); (8) WNT targets; (9) DNA damage repair; (10) mismatch repair; (11) Nucleotide excision repair; (12) DNA replication; (13) Antigen processing and presentation [33,34,35]. We them performed a correlation analysis to further reveal the association between m6A gene signature and some related biological pathways.
Collection of immune-checkpoint blockade genomic and clinical information
We performed a systematical search for the immune checkpoint blockade gene expression profiles, which could be publicly obtained and reported with complete clinical information. Two immunotherapeutic cohorts were finally included in our study: advanced urothelial cancer with intervention of atezolizumab, an anti-PD-L1 antibody (IMvigor210 cohort) , and metastatic melanoma treated with pembrolizumab, an anti-PD-1 antibody (GSE78220 cohort downloaded from GEO) . For IMvigor210 cohort, based on the Creative Commons 3.0 License, the complete expression data and detailed clinical annotations could be obtained from http://research-pub.Gene.com/imvigor210corebiologies. The raw count data were normalized by the DEseq2 R package and then the count value was transformed into the TPM value. For GSE78220 cohort, after standardization using limma package, the FPKM data of gene expression profiles was also converted to the more comparable TPM value among samples.
Correlations coefficients between the TME infiltrating immune cells and expression of m6A regulators were computed by Spearman and distance correlation analyses. One-way ANOVA and Kruskal-Wallis tests were used to conduct difference comparisons of three or more groups . On the basis of the correlation between m6Ascore and patients’ survival, the cut-off point of each dataset subgroup was determined using the survminer R package. The “surv-cutpoint” function, which repeatedly tested all potential cut points in order for finding the maximum rank statistic, was applied to dichotomize m6Ascore, and then patients were divided into high and low m6Ascore groups based on the maximally selected log-rank statistics to decrease the batch effect of calculation. The survival curves for the prognostic analysis were generated via the Kaplan-Meier method and log-rank tests were utilized to identify significance of differences. We adopted a univariate Cox regression model to calculate the hazard ratios (HR) for m6A regulators and m6A phenotype-related genes. The independent prognostic factors were ascertained through a multivariable Cox regression model. Patients with detailed clinical data were eligible for final multivariate prognostic analysis. The forestplot R package was employed to visualize the results of multivariate prognostic analysis for m6Ascore in ACRG cohort and TCGA-STAD cohort. The specificity and sensitivity of m6Ascore were assessed through receiver operating characteristic (ROC) curve, and the area under the curve (AUC) were quantified using pROC R package. The waterfall function of maftools package was used to present the mutation landscape in patients with high and low m6Ascore subtype in TCGA-STAD cohort. The R package of RCircos was adopted to plot the copy number variation landscape of 21 m6A regulators in 23 pairs of chromosomes . All statistical P value were two-side, with p < 0.05 as statistically significance. All data processing was done in R 3.6.1 software.
Landscape of genetic variation of m6A regulators in gastric cancer
A total of 21 m6A regulators including 8 writers, 2 erasers and 11 readers were finally identified in this study. Figure 1a summarized the dynamic reversible process of m6A RNA methylation mediated by regulators as well as their potential biological functions for RNA. We first summarized the incidence of copy number variations and somatic mutations of 21 m6A regulators in GC. Among 433 samples, 101 experienced mutations of m6A regulators, with frequency 23.33%. It was found that the ZC3H13 exhibited the highest mutation frequency followed by KIAA1429, while both demethylases (FTO and ALKBH5) as well as METTL3 did not show any mutations in GC samples (Fig. 1b). Further analyses revealed a significant mutation co-occurrence relationship between ELAVL1 and KIAA1429, YTHDF1 and ZC3H13, along with RBM15 and YTHDC1 (Figure S1B). The investigation of CNV alteration frequency showed a prevalent CNV alteration in 21 regulators and most were focused on the amplification in copy number, while ELAVL1, YTHDF2 and FMR1 had a widespread frequency of CNV deletion (Fig. 1c). The location of CNV alteration of m6A regulators on chromosomes was shown in Fig. 1d. Based on the expression of these 21 m6A regulators, we could completely distinguished GC samples from normal samples (Fig. 1e). To ascertain whether the above genetic variations influenced the expression of m6A regulators in GC patients, we investigated the mRNA expression levels of regulators between normal and GC samples, and found that the alterations of CNV could be the prominent factors resulting in perturbations on the m6A regulators expression. Compared to normal gastric tissues, m6A regulators with amplificated CNV demonstrated markedly higher expression in GC tissues (e.g. CBLL1 and FTO), and vice versa (e.g. ELAVL1 and YTHDF2) (Fig. 1c and f). The above analyses presented the highly heterogeneity of genetic and expressional alteration landscape in m6A regulators between normal and GC samples, indicating that the expression imbalance of m6A regulators played a crucial role in the GC occurrence and progression.
m6A methylation modification patterns mediated by 21 regulators
Five GEO datasets with available OS data and clinical information (GSE15459, GSE34942, GSE57303, GSE62254/ACRG and GSE84437, Table S1) were enrolled into one meta-cohort. A univariate Cox regression model revealed the prognostic values of 21 m6A regulators in patients with gastric cancer (Figure S1C). The comprehensive landscape of m6A regulator interactions, regulator connection and their prognostic significance for GC patients was depicted with the m6A regulator network (Fig. 2a and Table S3). We found that not only the m6A regulators in the same functional category presented a remarkably correlation in expression, but also a significant correlation was shown among writers, erasers, and readers. We also demonstrated that whether tumors with a high writer gene expression exhibits a low eraser gene expression actually depended on the different writer and eraser genes (Figure S2A-S2H). It was found that tumors with a high expression of writer genes (WATP and RBM15) showed a low expression of eraser gene FTO, while the high expression of WATP and RBM15 did not affect the expression of another eraser gene ALKBH5 (Figure S2A-S2B). Tumors with a high expression of writer gene METTL14, METTL3, KIAA1429 and ZC3H13 showed a high expression of eraser gene FTO, and METTL14, METTL3 and ZC3H13 also did not interfere with ALKBH5 expression, while KIAA1429 shared a common trend in gene expression with ALKBH5. In addition, the change of RBM15B expression did not affect the expression of these two eraser genes (Figure S2C-S2H). Considering the relatively higher mutation frequency of writer gene ZC3H13, we analyzed the difference in expression of eraser genes between ZC3H13-mutant and wild types. Of these, ALKBH5 was significantly up-regulated in ZC3H13-mutant tumors compared to wild-type tumors, while FTO was significantly down-regulated (Figure S2I).
The above results indicated that cross-talk among the regulators of writers, readers, and erasers may play critical roles in the formation of different m6A modification patterns and TME cell-infiltrating characterization between individual tumors.
The R package of ConsensusClusterPlus was used to classify patients with qualitatively different m6A modification patterns based on the expression of 21 m6A regulators, and three distinct modification patterns were eventually identified using unsupervised clustering, including 389 cases in pattern A, 348 cases in pattern B and 322 cases in pattern C. We termed these patterns as m6Acluster A-C, respectively (Figure S2J and Table S4). Prognostic analysis for the three main m6A modification subtypes revealed the particularly prominent survival advantage in m6Acluster-B modification pattern (Fig. 2b).
TME cell infiltration characteristics in distinct m6A modification patterns
To explore the biological behaviors among these distinct m6A modification patterns, we performed GSVA enrichment analysis. As shown in Fig. 2c and Table S5, m6Acluster-A was markedly enriched in stromal and carcinogenic activation pathways such as ECM receptor interaction, TGF beta signaling pathway, cell adhesion and MAPK signaling pathways. m6Acluster-B presented enrichment pathways associated with immune fully activation including the activation of chemokine signaling pathway, cytokine-cytokine receptor interaction, T cell receptor signaling pathway and Toll like receptor signaling pathways (Fig. 2c). While m6Acluster-C was prominently related to immune suppression biological process (Fig. 2d). To our surprise, subsequent analyses of TME cell infiltration indicated m6Acluster-A was remarkably rich in innate immune cell infiltration including natural killer cell, macrophage, eosinophil, mast cell, MDSC, plasmacytoid dendritic cell (Fig. 3a and Table S4). However, patients with this m6A modification pattern did not show a matching survival advantage (Fig. 2b). Previous studies demonstrated that tumors with immune-excluded phenotype also showed the presence of abundant immune cells, while these immune cells were retained in the stroma surrounding tumor cell nests rather than penetrate their parenchyma. The activation of stroma in TME were considered T-cell suppressive . The results from GSVA analyses have revealed cluster A modification pattern was significantly associated with stromal activation. Therefore, we speculated that stromal activation in cluster A inhibited the antitumor effect of immune cells. Subsequent analyses showed that stroma activity was significantly enhanced in cluster A such as the activation of epithelial-mesenchymal transition (EMT), transforming growth factor beta (TGFb) and angiogenesis pathways, which confirmed our speculation (Fig. 3b) Based on the above analyses, we were surprised to find three m6A modification patterns had significantly distinct TME cell infiltration characterization. Cluster A was classified as immune-excluded phenotype, characterized by innate immune cell infiltration and stromal activation; cluster B was classified as immune-inflamed phenotype, characterized by adaptive immune cell infiltration and immune activation; cluster C was classified as immune-desert phenotype, characterized by the suppression of immunity (Figs. 2c-d and 3a-b). We then used the CIBERSORT method, a deconvolution algorithm using support vector regression for determining the immune cell type in tumors, to compare the component differences of immune cells among the three m6A modification patterns. We found that there were no significant differences on the compositions of TME cell types between the three m6A modification patterns, which suggested that m6A methylation modification did not change TME infiltrating-cell types of tumors (Figure S2K).
We then examined the specific correlation between each TME infiltration cell type and each m6A regulator using spearman’s correlation analyses (Figure S3A). We focused on the regulator KIAA1429, a m6A methyltransferases, and revealed its significantly negative correlation with numerous TME infiltrating immune cells. We used ESTIMATE algorithm to quantify the overall infiltration of immune cells between high and low KIAA1429 expression patients. The results showed that low expression of KIAA1429 exhibited high immune scores, which meant that the TME with low expression of KIAA1429 existed a significantly increased immune cell infiltration, thus confirming the above findings (Figure S3B). We then explored the specific difference of 23 TME infiltrating immune cells between high and low KIAA1429 expression patients. We found tumors with low expression of KIAA1429 presented significantly increased infiltration in 23 TME immune cells compared to patients with high expression (Figure S3C). Recent studies paid special attention to the mechanism of m6A modification regulating the activation of dendritic cells (DCs). DCs, which are responsible for antigen presentation and the activation of naive T cells, are a bridge connecting innate and adaptive immunity, and their activation depending on the high expression level of MHC molecules, costimulatory factors and adhesion factors . Our study indicated that tumors with low expression of KIAA1429 showed significant more enrichment of TME DCs infiltration including activated DCs, immature DCs, and plasmacytoid DCs. We also noted that the decreased expression of KIAA1429 resulted in the comprehensively elevated expression of MHC molecules, costimulatory molecules, and adhesion molecules (Figure S3D). Subsequent pathway enrichment analyses, as expected, tumors with low KIAA1429 expression exhibited an obvious enhancement in immune activation pathways including the pathway of antigen processing and presentation, C-type lectin receptor, NOD-like receptor, T cell receptor, Toll-like receptor and NF-κB signaling pathway (Figure S3E). It was interesting that the immune-related pathway enhancements were accompanied by the increased expression of immunological checkpoint molecules PD1/L1 (Figure S3D-S3E). So we investigated whether the expression of KIAA1429 regulator affected the therapeutic efficacy of immune checkpoint blockade. In anti-PD-L1 immunotherapy cohort (IMvigor210), a survival benefit trend was observed in patients with low expression of KIAA1429 (Figure S3F). In anti-PD-1 immunotherapy cohort (GSE78220), we did not observe a significantly prolonged survival owing to the few samples (Figure S3G). From above, we could speculate that KIAA1429-mediated m6A methylation modification may promote the activation of TME DCs, thus enhancing the intratumoral antitumor immune response.
m6A methylation modification patterns in ACRG cohort
To further explore the characteristics of these m6A modification phenotypes in the different clinical traits and biological behaviors, we fixed attention on the ACRG cohort, which comprised 300 gastric cancer patients and offered the most comprehensive clinical annotation. Similar to all GC datasets clustering, unsupervised clustering also discovered three fully distinct patterns of m6A modification in ACRG cohort (Figure S4A-S4D and Fig. 3c-d). There was significant distinction existed on the m6A transcriptional profile among three different m6A modification patterns (Fig. 3d). m6Acluster A was characterized by the increased expression of FTO and HNRNPA2B1, and presented variable decreases in other m6A regulators; m6Acluster B showed high expression of ELAVL1, HNRNPC, LRPPRC, METTL14, RBM15, RBM15B, YTHDC2 and YTHDF2; and m6Acluster C exhibited significant increases in the expression of FMR1 IGF2BP1, WTAP, ZC3H13 and YTHDF1. Patients with EMT molecular subtypes were characterized by the m6Acluster-A methylation modification patterns, while MSI subtypes were characterized by the m6Acluster-B modification patterns. We also noted that tumors with m6Acluster-A patterns presented poorer differentiation and were enriched in the diffuse histological subtype. A better tumor differentiation was observed in the m6Acluster-B and m6Acluster-C patterns, which were enriched in the intestinal histological subtype. In gastric cancer, the EMT molecular subtype and diffuse histological type was markedly linked to a poorer survival, while MSI linked to a better clinical outcome. Therefore, the tumors characterized by m6Acluster-A modification patterns were significantly correlated with stromal activation, high malignancy and rapid progression (Fig. 3c). One-way ANOVA test also confirmed the remarkable differences on m6A regulator expression between three key m6A modification patterns. Prognostic analysis also revealed m6Acluster B to be markedly related to prolonged survival, while m6Acluster A and m6Acluste C were characterized by poorer survival (Figure S4E-S4F). Consistent with the above findings, most patients with EMT subtypes were clustered into m6Acluster A and almost no EMT subtypes were in m6Acluster B, which confirmed again that m6Acluster A was significantly relevant to the stromal activation and m6Aclustre B relevant to the immune activation (Fig. 3e and Table S6).
Generation of m6A gene signatures and functional annotation
To further investigate the potential biological behavior of each m6A modification pattern, we determined 718 m6A phenotype-related DEGs using limma package (Figure S4G and Table S7). The clusterProfiler package was used to perform GO enrichment analysis for the DEGs. The biological processes with significant enrichment were summarized in Table S8. Surprisingly, these genes showed enrichment of biological processes remarkably related to m6A modification and immunity, which confirmed again that m6A modification played a non-negligible role in the immune regulation in tumor microenvironment (Fig. 3f). To further validate this regulation mechanism, we then performed unsupervised clustering analyses based on the obtained 718 m6A phenotype-related genes in order to classify patients into different genomic subtypes. Consistent with the clustering grouping of m6A modification patterns, the unsupervised clustering algorithm also revealed three distinct m6A modification genomic phenotypes and we named these three clusters as m6A gene cluster A-C, respectively (Figure S5A-S5D, Fig. 4a and Table S6). This demonstrated that three distinct m6A methylation modification patterns did exist in gastric cancer. We observed that tumors in m6A gene cluster C patterns also exhibited poorer differentiation and were enriched in the diffuse histological subtype. The opposite patterns were observed in m6A gene cluster A and cluster B. Patients with alive status or MSI subtypes were mainly concentrated in the m6A gene cluster A, while patients with clinical stage IV or EMT molecular subtypes were characterized by the m6A gene cluster C patterns (Fig. 4a). Analysis also indicated three distinct gene clusters were characterized by different signature genes (Fig. 4a). Eighty-eight of three hundred patients with gastric cancer were clustered in gene cluster A, which were proved to be related to better prognosis. While patients in gene cluster C (105 patients) experienced the outcome of poorer prognosis. An intermediate prognosis was observed in gene cluster C, with 107 patients clustered (Fig. 4b). In the three m6A gene clusters, the prominent differences in the expression of m6A regulators were observed, which was in accordance with the expected results of m6A methylation modification patterns (Fig. 4c).
Characteristics of clinical and transcriptome traits in m6A-related phenotypes
To reveal the role of m6A-related phenotypes in the TME immune regulation, we studied the expression of chemokine and cytokine characterizing three gene clusters. The selected cytokine and chemokine were extracted from published literature, of which, TGRB1, SMAD9, TWIST1, CLDN3, TGFBR2, ACTA2, COL4A1, ZEB1 and VIM were considered to be associated with the transcripts of transforming growth factor (TGF)b/EMT pathway. PD-L1, CTLA-4, IDO1, LAG3, HAVCR2, PD-1, PD-L2, CD80, CD86, TIGIT and TNFRSF9 were considered to be related to the transcripts of immune checkpoints. TNF, IFNG, TBX2, GZMB, CD8A, PRF1, GZMA, CXCL9 and CXCL10 were to be correlated with the transcripts of immune activation [29, 32]. We found the mRNAs relevant to TGFb/EMT pathway were significantly upregulated in gene cluster C, which demonstrated that this cluster was deemed as stromal-activated group. While gene cluster A showed high expression of mRNAs related to immune activation transcripts. This suggested that gene cluster A could be classified as the immune-activation group (Figure S5F-S5H). To better depict the function of m6A signature genes, we examined the known signatures in patients with gastric cancer (Figure S5E). The results also confirmed that gene cluster C was characterized by the status of stromal activation and cancer promotion, and gene cluster A was significantly related to immune activation status (Figure S5E-S5H). Consistent with the above findings, as shown in Fig. 4d and Table S6, almost all (41 of 45, 91%) patients with EMT subtype (molecular subtypes in ACRG cohort) were classified into gene cluster C, which was relevant to the worse survival outcome.
The above results showed again that m6A methylation modification played a non-negligible regulation role in shaping different TME landscapes. However, these analyses were only based on the patient population and could not accurately predict the pattern of m6A methylation modification in individual patients. Considering the individual heterogeneity and complexity of m6A modification, based on these phenotype-related genes, we constructed a set of scoring system to quantify the m6A modification pattern of individual patients with gastric cancer, we termed as m6Ascore. The alluvial diagram was used to visualize the attribute changes of individual patients (Fig. 4d). To better illustrate the characteristics of m6A signature, we also tested the correlation between the known signatures and the m6Ascore (Fig. 4e and Table S9). Kruskal-Wallis test revealed significant difference on m6Ascore between m6A gene clusters. Gene cluster A showed the lowest median score while gene cluster C had the highest median score, which indicated that low m6Ascore could be closely linked to immune activation-related signatures, whereas high m6Ascore could be linked to stromal activation-related signatures (Fig. 4f). More importantly, m6Acluster A showed the significantly increased m6Ascore compared to the other clusters and m6Acluster B presented the lowest median score (Fig. 4g). The analyses for the activity of stroma-related pathways indicated high scores were significantly associated with enhanced activation of stromal pathways (Fig. 4h). In addition, patients with EMT subtypes also showed the lowest m6Ascore compared to other three ACRG molecular subtypes (Fig. 5a). The above results strongly suggested that low m6Ascore was significantly correlated with immune-activation and high m6Ascore was correlated with stromal-activation. The m6Ascore could better evaluate the m6A modification patterns of individual tumor, and further evaluate tumors’ TME cell-infiltration characterization, in order to distinguish the true and false nature of TME immune infiltration.
Next, we sought to further identify the value of m6Ascore in predicting patients’ outcome. With the cutoff value 0.0291 determined by survminer package, patients were divided into low or high m6Ascore group. Patients with low m6Ascore demonstrated a prominent survival benefit (HR 3.0 (2.12–4.21); Fig. 5b), with 5-year survival rate twice than patients with high m6Ascore (69.4% vs 33.5%). We tested whether the m6Ascore could serve as an independent prognostic biomarker for gastric cancer. Multivariate Cox regression model analysis, which included the factors of patients’ age, gender, TNM status, histological type, MSI status, TP53 status and ACRG molecular subtypes, confirmed m6Ascore as a robust and independent prognostic biomarker for evaluating patient outcomes (HR 2.54(1.71–3.8); Figure S6A). We specifically examined the ability of m6Ascore signature to predict the efficacy of adjuvant chemotherapy in patients with gastric cancer. We found that patients with low m6Ascore showed significant therapeutic advantages among patients who also received adjuvant chemotherapy, with 5-year survival rate 77.5% vs 59.2% (Fig. 5c). Another results obtained indicated that the prediction power of m6Ascore was not interfered by adjuvant chemotherapy, and both in patients receiving chemotherapy or not, low m6Ascore group always showed the obvious survival advantage (Fig. 5c). In addition, we revealed that younger patients, diffuse histological subtype and advanced patients were significantly associated with a higher m6Ascore, which meant that these patients were characterized by the m6Acluster-A modification patterns and immune-excluded phenotype, with a poorer clinical outcome. These results demonstrated m6Ascore could be also used to evaluate certain clinical characteristics of patients such as MSI status, molecular subtypes, histological subtypes as well as clinical stage, etc. (Figure S6B).
Characteristics of m6A modification in TCGA molecular subtypes and tumor somatic mutation
A comprehensive molecular landscape has been constructed for gastric cancer by TCGA project, which classified gastric cancer into four molecular subtypes including genome stable (GS), microsatellite instability (MSI), EBV infection, and chromosomal instability (CIN). We evaluated the difference of m6Ascore between these molecular subtypes. The higher m6Ascore was obviously concentrated on GS subtype and showed a worse survival in patients, while the lower m6Ascore was concentrated on the subtypes of MSI and EBV infection, which was related to better survival (5-year survival rate, 25.9% vs 43.3%; HR 1.81(1.26–2.62); Fig. 5d-e). The highly microsatellite instability subtype, characterized by better prognosis, was significantly correlated with lower m6Ascore, whereas MSI-Low and MSS had a higher m6Ascore (Fig. 5f). Multivariate analysis for TCGA-STAD cohort also confirmed that m6Ascore could act as an independent prognostic biomarker in gastric cancer (Figure S6C). Previous studies indicated that patients with EBV-positive gastric cancer have been shown to respond to anti-PD-1/L1 antibodies in several studies in spite of the lower MSI or tumor mutation burden (TMB) [41, 42]. In our study, EBV infected patients were markedly associated with lower m6Ascore than CIN and GS subtypes as well as EBV non-infected patients, which implied m6Ascore signature could be a more effective biomarker for the prediction of immunotherapeutic efficacy than MSI and TMB in patients with gastric cancer (Fig. 5e-f). Further research showed that tumors with MSI subtype were mainly characterized by the m6Acluster-B methylation modification patterns, while tumors with MSS subtype were characterized by the m6Acluster-C modification patterns (Figure S7A). The m6A regulators ALKBH5, CBLL1, ELAVL1, FMR1, HNRNPC, KIAA1429, METTL14, RBM15, RBM15B, WTAP, YTHDC1, YTHDC2, YTHDF2 and YTHDF3 were significantly up-regulated in MSI subtypes compared to MSS subtypes, while IGF2BP1, YTHDF1 and ZC3H13 were markedly down-regulated (Figure S7B). For EB virus infection, patients with EBV-positive were mainly characterized by the m6Acluster-A methylation modification patterns, while EBV-negative patients did not show a characteristic pattern of m6A methylation modification. There were no significant difference on m6A modification patterns between EBV-negative patients (Figure S7C). In addition, we found the m6A regulators IGF2BP1, KIAA1429, LRPPRC, YTHDF3 and ZC3H13 were remarkably down-regulated in EBV-negative patients than EBV positive patients, while FTO was significantly down-regulated (Figure S7D). The above results suggested that the potential mechanisms on the change of m6A modification patterns mediated by EBV and MSI etc. may be that these factors changed the status of m6A regulators. These findings could contribute to enhancing our understanding of the mechanisms of the formation of m6A modification pattern differences in tumors.
Then, we analyzed the distribution differences of somatic mutation between low and high m6Ascore in TCGA-STAD cohort using maftools package. As shown in Fig. 5g-h, low m6Ascore group presented more extensive tumor mutation burden than the high m6Ascore group, with the rate of the 10th most significant mutated gene 25% versus 10%. The TMB quantification analyses confirmed the low m6Ascore tumors was markedly correlated with a higher TMB (Figure S7E). The m6Ascore and TMB also exhibited a significant negative correlation (Figure S7F). Accumulated evidence demonstrated patients with high TMB status presented a durable clinical response to anti-PD-1/PD-L1 immunotherapy. Therefore, the above results indirectly demonstrated that the difference in tumor m6A modification patterns could a crucial factor that mediated the clinical response to anti-PD-1/PD-L1 immunotherapy. And the values of m6Ascore in predicting immunotherapeutic outcomes were also indirectly confirmed.
The clinical trials as well as preclinical researches have revealed patients with higher somatic TMB were correlated with enhanced response, long-term survival and durable clinical benefit when treated with immune checkpoint blockade therapy. The individual altered genes could mediate resistance or sensitivity to immunotherapy. For specific altered genes in TCGA-TSAD such as ARID1A and PIK3CA, mutant type had significantly lower m6Ascore compared to wild type, whereas there was no significant difference in m6Ascore between wild and mutant types in TP53 and RHOA (Fig. 5f). These results would provide novel perspective for exploring the mechanisms of m6A methylation modification in the tumor somatic mutations, shaping of TME landing, and roles in immune checkpoint blockade therapy.
m6A modification patterns in the role of anti-PD-1/L1 immunotherapy
In order to further test the stability of m6Ascore model, we applied m6Ascore signature established in ACRG cohort to other independent gastric cancer cohorts to verify its prognostic value (GSE84437, HR 1.89(1.44–2.49); GSE15459, HR 2.05(1.30–3.22); GSE34942, HR 1.52(0.64–3.63); Figure S8A-S8C). The combined set of all GEO cohorts was validated (HR 1.94(1.62–2.31); Figure S8D). The ability of m6Ascore to predict relapse-free survival was also evaluated (GSE26253, HR 1.33(0.98–1.80); GSE62254, HR 2.53(1.75–3.65); Figure S8E-S8F). Next, we continued to extend the m6Ascore signature to all digestive system tumors including cholangiocarcinoma, colon adenocarcinoma, pancreatic adenocarcinoma, esophageal carcinoma and liver hepatocellular carcinoma (HR 1.4(1.17–1.68); Figure S8G). These data indicated m6A modification patterns correlated with better clinical benefit. The predictive advantage evaluated with ROC curves was especially reflected in elderly patients (Figure S8H-S8I).
Immunotherapies represented by PD-L1 and PD-1 blockade has undoubtedly emerged a major breakthrough in cancer therapy. We investigated whether the m6A modification signature could predict patients’ response to immune checkpoint blockade therapy based on two immunotherapy cohorts. In both anti-PD-L1 cohort (IMvigor210) and anti-PD-1 cohort (GSE78220), patients with low m6Ascore exhibited significantly clinical benefits and a markedly prolonged survival (Fig. 6a-g; IMvigor210, HR 1.73(1.20–2.48), Fig. 6a; GSE78220, HR 4.58(1.23–17.10), Fig. 6d). The significant therapeutic advantages and clinical response to anti-PD-1/L1 immunotherapy in patients with low m6Ascore compared to those with high m6Ascore were confirmed (Fig. 6b-c and e-g). In addition, patients with low m6Ascore showed a obviously high expression of PD-L1, which indicated a potential response to anti-PD-1/L1 immunotherapy (Fig. 6h). Further research revealed that regulatory T-cells and TME stroma were significantly activated in tumors with high m6Ascore, which mediated immune tolerance of tumors (Fig. 6i). Tumor neoantigen burden, closely linked to immunotherapeutic efficacy, was also assessed. We found patients with combination of low m6Ascore and high neoantigen burden showed a great survival advantage (Fig. 6j). The above implied that the quantification of m6A modification patterns was a potential and robust biomarker for prognosis and clinical response assessment of immunotherapy (Fig. 6k). The immune phenotypes of tumors in the IMvigor210 cohort has been detected, so we investigated the difference of m6Ascore among different phenotypes. We found that higher m6Ascore was remarkably associated with exclusion and desert immune phenotypes, and checkpoint inhibitors were difficult to exert antitumor effect in these phenotype (Fig. 6l). In summary, our work strongly indicated that m6A methylation modification patterns was significantly correlated with tumor immune phenotypes and response to anti-PD-1/L1 immunotherapy, and the established m6A modification signature would contribute to predicting the response to anti-PD-1/L1 immunotherapy.
Increasing evidence demonstrated that m6A modification took on an indispensable role in inflammation, innate immunity as well as antitumor effect through interaction with various m6A regulators. As most studies focus on single TME cell type or single regulator, the overall TME infiltration characterizations mediated by integrated roles of multiple m6A regulators are not comprehensively recognized. Identifying the role of distinct m6A modification patterns in the TME cell infiltration will contribute to enhancing our understanding of TME antitumor immune response, and guiding more effective immunotherapy strategies.
Here, based on 21 m6A regulators, we revealed three distinct m6A methylation modification patterns. These three patterns had significantly distinct TME cell infiltration characterization. Cluster A was characterized by the activation of innate immunity and stroma, corresponding to immune-excluded phenotype; cluster B was characterized by the activation of adaptive immunity, corresponding to immune-inflamed phenotype; cluster C was characterized by the suppression of immunity, corresponding to immune-desert phenotype. The immune-excluded and immune-desert phenotypes could be regarded as non-inflamed tumors. The immune-inflamed phenotype, known as hot tumor, show by a large number of immune cell infiltration in TME [39, 43, 44]. Although the immune-excluded phenotype also showed the presence of abundant immune cells, the immune cells were retained in the stroma surrounding tumor cell nests rather than penetrate their parenchyma. The stroma could be confined to the tumor envelope or may penetrate the tumor itself, making the immune cells appear to be really inside the tumor [45,46,47]. The immune-desert phenotypes were associated with immune tolerance and ignorance, and lack of activated and priming T-cell . Consistent with the above definitions, we found cluster A exhibited a significant stroma activation status, including the highly expressed angiogenesis, EMT and TGF-β pathways, which were considered T-cell suppressive. Combined with the TME cell-infiltrating characteristics in each cluster, it confirmed the reliability of our classification of immune phenotypes for different m6A modification patterns. Therefore, after fully exploring the TME cell–infiltrating characterization induced by distinct m6A modification patterns, it was not surprising that cluster A had the activated innate immunity but poorer prognosis.
Further, in this study, the mRNA transcriptome differences between distinct m6A modification patterns have been proved to be significantly associated with m6A and immune related biological pathways. These differentially expressed genes were considered as m6A-related signature genes. Similar to the clustering results of the m6A modification phenotypes, three genomic subtypes were identified based on m6A signature genes, which were also significantly correlated with stromal and immune activation. This demonstrated again that the m6A modification was of great significance in shaping different TME landscapes. Therefore, a comprehensive assessment of the m6A modification patterns will enhance our understanding of TME cell-infiltrating characterization. Considering the individual heterogeneity of m6A modification, it was urgently demanded to quantify the m6A modification patterns of individual tumor. For that, we established a set of scoring system to evaluate the m6A modification pattern of individual patients with gastric cancer—the m6A gene signature. The m6A modification pattern characterized by immune-excluded phenotype exhibited a higher m6Ascore, while the pattern characterized by immune-inflamed phenotype showed a lower m6Ascore. In addition, In IMvigor210 cohort with the determined immune phenotype, these results were well validated . This suggested m6Ascore was a reliable and robust tool for comprehensive assessment of individual tumor m6A modification patterns, which could be used to further determine the TME infiltration patterns, that was, tumor immune phenotypes. Integrated analyses also demonstrated that m6Ascore was an independent prognostic biomarker in gastric cancer. Patients with EBV and MSI subtypes, sensitive to checkpoint immunotherapy , was significantly related to lower m6Ascore. Considering the low mutation burden but high immune infiltration in EBVþ tumors , our m6Ascore showed a predictive advantage in precision immunotherapy for gastric cancer.
Our data also revealed a markedly negative correlation between m6Ascore and tumor mutation burden. Consistent with previous studies, EMT and GS molecular subtypes demonstrated the lowest m6Ascore, underlining the core role of stromal activation in resistance to checkpoint immunotherapy [49, 50]. This indicated that response to checkpoint immunotherapy was not only associated with antigen processing, and improved cytolytic activity, but also related to suppression of angiogenesis, fibroblast activation, TGF beta pathway components and the EMT. Previous studies confirmed that the EMT- and TGFbeta-related pathway activation resulted in decreased trafficking of T-cell into tumors as well as their weakened tumor killing effects [33, 49]. The above suggested that the activated stromal TME in the activated immune TME could mediate therapeutic resistance to immune-checkpoint blockade, as well as influence the individual precise immunotherapy of gastric cancer. In this work, we showed m6A methylation modification patterns played a nonnegligible role in shaping different stromal and immune TME landscape, implying m6A modification could affect the therapeutic efficacy of immune checkpoint blockade. The m6A gene signature with integrated various biomarkers including mutation load, neoantigen load, PD-L1 expression, stromal and immune TME and MSI status, could be the more effective predictive strategy for immunotherapy. We also confirmed the predictive value of the m6Ascore in two cohort with anti-PD-1 and anti-PD-L1 immunotherapy. A significantly difference on m6Ascores existed between non-responders and responders.
In short, in clinical practice, the m6Ascore could be used to comprehensively evaluate the m6A methylation modification patterns as well as their corresponding TME cell infiltration characterization within individual patient, further to determine the immune phenotypes of tumors and guide the more effective clinical practice. We also demonstrated the m6Ascore could be utilized for assessing patients’ clinicopathological features including stages of tumor inflammation, tumor differentiation levels, clinical stages, histological subtypes, molecular subtypes, genetic variation, MSI status, EBV infection and tumor mutation burden etc. The detailed relationships between m6Ascore and clinicopathological features could be found in our study. Similarly, m6Ascore could act as an independent prognostic biomarker for predicting patients’ survival. We could also predict the efficacy of adjuvant chemotherapy and the patients’ clinical response to anti-PD-1/PD-L1 immunotherapy through m6Ascore. More importantly, this study has yielded several novel insights for cancer immunotherapy that targeting m6A regulators or m6A phenotype-related genes for changing the m6A modification patterns, and further reversing the adverse TME cell infiltration characterization, that was the transformation of “cold tumors” into “hot tumors”, may contribute to exploiting the development of novel drug combination strategies or novel immunotherapeutic agents in the future. Our findings provided novel ideas for improving the patients’ clinical response to immunotherapy, identifying different tumor immune phenotypes and promoting personalized cancer immunotherapy in the future.
In conclusion, this work demonstrated the extensive regulation mechanisms of m6A methylation modification on tumor microenvironment. The difference of m6A modification patterns was a factor that could not be ignored to cause the heterogeneity and complexity of individual tumor microenvironment. The comprehensive evaluation of individual tumor m6A modification pattern will contribute to enhancing our understanding of tumor microenvironment cell-infiltrating characterization and guiding more effective immunotherapy strategies.
Copy number variation
Differentially expressed genes
Gene set variation analysis
Immunological checkpoint blockade
Pan-fibroblast TGFb response signature
Principal component analysis
Single-sample gene-set enrichment analysis
The Cancer Genome Atlas
Transforming growth factor beta
Boccaletto P, Machnicka MA, Purta E, Piatkowski P, Baginski B, Wirecki TK, de Crecy-Lagard V, Ross R, Limbach PA, Kotter A, et al. MODOMICS: a database of RNA modification pathways. 2017 update. Nucleic Acids Res. 2018;46:D303–7.
Roundtree IA, Evans ME, Pan T, He C. Dynamic RNA modifications in gene expression regulation. Cell. 2017;169:1187–200.
Alarcon CR, Lee H, Goodarzi H, Halberg N, Tavazoie SF. N6-methyladenosine marks primary microRNAs for processing. Nature. 2015;519:482–5.
Patil DP, Chen CK, Pickering BF, Chow A, Jackson C, Guttman M, Jaffrey SR. m(6)A RNA methylation promotes XIST-mediated transcriptional repression. Nature. 2016;537:369–73.
Zhao BS, Roundtree IA, He C. Post-transcriptional gene regulation by mRNA modifications. Nat Rev Mol Cell Biol. 2017;18:31–42.
Yang Y, Hsu PJ, Chen YS, Yang YG. Dynamic transcriptomic m(6)A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res. 2018;28:616–24.
Chen XY, Zhang J, Zhu JS. The role of m(6)A RNA methylation in human cancer. Mol Cancer. 2019;18:103.
Tong J, Flavell RA, Li HB. RNA m(6)A modification and its function in diseases. Front Med. 2018;12:481–9.
Liu F, Clark W, Luo G, Wang X, Fu Y, Wei J, Wang X, Hao Z, Dai Q, Zheng G, et al. ALKBH1-mediated tRNA demethylation regulates translation. Cell. 2016;167:816–28 e816.
Wang S, Sun C, Li J, Zhang E, Ma Z, Xu W, Li H, Qiu M, Xu Y, Xia W, et al. Roles of RNA methylation by means of N(6)-methyladenosine (m(6)A) in human cancers. Cancer Lett. 2017;408:112–20.
Wang CX, Cui GS, Liu X, Xu K, Wang M, Zhang XX, Jiang LY, Li A, Yang Y, Lai WY, et al. METTL3-mediated m6A modification is required for cerebellar development. PLoS Biol. 2018;16:e2004880.
Fu Y, Dominissini D, Rechavi G, He C. Gene expression regulation mediated through reversible m(6)A RNA methylation. Nat Rev Genet. 2014;15:293–306.
Pinello N, Sun S, Wong JJ. Aberrant expression of enzymes regulating m(6)A mRNA methylation: implication in cancer. Cancer Biol Med. 2018;15:323–34.
Tong J, Cao G, Zhang T, Sefik E, Amezcua Vesely MC, Broughton JP, Zhu S, Li H, Li B, Chen L, et al. m(6)A mRNA methylation sustains Treg suppressive functions. Cell Res. 2018;28:253–6.
Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, Powderly JD, Carvajal RD, Sosman JA, Atkins MB, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med. 2012;366:2443–54.
Pitt JM, Marabelle A, Eggermont A, Soria JC, Kroemer G, Zitvogel L. Targeting the tumor microenvironment: removing obstruction to anticancer immune responses and immunotherapy. Ann Oncol. 2016;27:1482–92.
Quail DF, Joyce JA. Microenvironmental regulation of tumor progression and metastasis. Nat Med. 2013;19:1423–37.
Ali HR, Chlon L, Pharoah PD, Markowetz F, Caldas C. Patterns of immune infiltration in breast cancer and their clinical implications: a gene-expression-based retrospective study. PLoS Med. 2016;13:e1002194.
Binnewies M, Roberts EW, Kersten K, Chan V, Fearon DF, Merad M, Coussens LM, Gabrilovich DI, Ostrand-Rosenberg S, Hedrick CC, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med. 2018;24:541–50.
Fang H, Declerck YA. Targeting the tumor microenvironment: from understanding pathways to effective clinical trials. Cancer Res. 2013;73:4965–77.
Han D, Liu J, Chen C, Dong L, Liu Y, Chang R, Huang X, Liu Y, Wang J, Dougherty U, et al. Anti-tumour immunity controlled through mRNA m(6)A methylation and YTHDF1 in dendritic cells. Nature. 2019;566:270–4.
Wang H, Hu X, Huang M, Liu J, Gu Y, Ma L, Zhou Q, Cao X. Mettl3-mediated mRNA m(6)A methylation promotes dendritic cell activation. Nat Commun. 2019;10:1898.
Wang Q, Chen C, Ding Q, Zhao Y, Wang Z, Chen J, Jiang Z, Zhang Y, Xu G, Zhang J, et al. METTL3-mediated m(6)A modification of HDGF mRNA promotes gastric cancer progression and has prognostic significance. Gut. 2019. [Epub ahead of print].
Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I, et al. TCGAbiolinks: an R/bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016;44:e71.
Hartigan JAWM. Algorithm AS 136: a K-means clustering algorithm. Appl Stat. 1979;28:9.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26:1572–3.
Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.
Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 2017;18:248–62.
Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, Schinzel AC, Sandy P, Meylan E, Scholl C, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009;462:108–12.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.
Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J, Nordgren H, Farmer P, Praz V, Haibe-Kains B, 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.
Zeng D, Li M, Zhou R, Zhang J, Sun H, Shi M, Bin J, Liao Y, Rao J, Liao W. Tumor microenvironment characterization in gastric cancer identifies prognostic and immunotherapeutically relevant gene signatures. Cancer Immunol Res. 2019;7:737–50.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE III, Koeppen H, Astarita JL, Cubas R, et al. TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554:544–8.
Rosenberg JE, Hoffman-Censits J, Powles T, van der Heijden MS, Balar AV, Necchi A, Dawson N, O'Donnell PH, Balmanoukian A, Loriot Y, et al. Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: a single-arm, multicentre, phase 2 trial. Lancet. 2016;387:1909–20.
Senbabaoglu Y, Gejman RS, Winer AG, Liu M, Van Allen EM, de Velasco G, Miao D, Ostrovnaya I, Drill E, Luna A, et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. 2016;17:231.
Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, Berent-Maoz B, Pang J, Chmielowski B, Cherry G, et al. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell. 2016;165:35–44.
Hazra A, Gogtay N. Biostatistics series module 3: comparing groups: numerical variables. Indian J Dermatol. 2016;61:251–60.
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28:1747–56.
Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. 2017;541:321–30.
Qian C, Cao X. Dendritic cells in the regulation of immunity and inflammation. Semin Immunol. 2018;35:3–11.
Kim ST, Cristescu R, Bass AJ, Kim KM, Odegaard JI, Kim K, Liu XQ, Sher X, Jung H, Lee M, et al. Comprehensive molecular characterization of clinical responses to PD-1 inhibition in metastatic gastric cancer. Nat Med. 2018;24:1449–58.
Panda A, Mehnert JM, Hirshfield KM, Riedlinger G, Damare S, Saunders T, Kane M, Sokol L, Stein MN, Poplin E, et al. Immune activation and benefit from avelumab in EBV-positive gastric cancer. J Natl Cancer Inst. 2018;110:316–20.
Turley SJ, Cremasco V, Astarita JL. Immunological hallmarks of stromal cells in the tumour microenvironment. Nat Rev Immunol. 2015;15:669–82.
Gajewski TF, Woo SR, Zha Y, Spaapen R, Zheng Y, Corrales L, Spranger S. Cancer immunotherapy strategies based on overcoming barriers within the tumor microenvironment. Curr Opin Immunol. 2013;25:268–76.
Gajewski TF. The next hurdle in cancer immunotherapy: overcoming the non-T-cell-inflamed tumor microenvironment. Semin Oncol. 2015;42:663–71.
Joyce JA, Fearon DT. T cell exclusion, immune privilege, and the tumor microenvironment. Science. 2015;348:74–80.
Salmon H, Franciszkiewicz K, Damotte D, Dieu-Nosjean MC, Validire P, Trautmann A, Mami-Chouaib F, Donnadieu E. 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.
Kim JM, Chen DS. Immune escape to PD-L1/PD-1 blockade: seven steps to success (or failure). Ann Oncol. 2016;27:1492–504.
Tauriello DVF, Palomo-Ponce S, Stork D, Berenguer-Llergo A, Badia-Ramentol J, Iglesias M, Sevillano M, Ibiza S, Canellas A, Hernando-Momblona X, et al. TGFbeta drives immune evasion in genetically reconstituted colon cancer metastasis. Nature. 2018;554:538–43.
Voron T, Marcheteau E, Pernot S, Colussi O, Tartour E, Taieb J, Terme M. Control of the immune response by pro-angiogenic factors. Front Oncol. 2014;4:70.
This study was funded by the Natural Science Foundation of China (81871763), KeyResearch and Development Projects of Jiangsu Province (BE2017681), Jiangsu Youth Medical Talent Project (QNRC2016702) and the National Undergraduate Training Programs for Innovation (201910304030Z and 201910304028Z).
Ethics approval and consent to participate
The patient data in this work were acquired from the publicly available datasets whose informed consent of patients were complete.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Overview of study design and prognostic analysis of 21 m6A regulators. Figure S2. Correlation between writer gene expression/mutation and eraser gene expression. Figure S3. Correlation between TME infiltration cells and m6A regulators and the roles of KIAA1429 in activation of dendritic cells. Figure S4. Unsupervised clustering of 21 m6A regulators in the ACRG gastric cancer cohort. Figure S5. Characteristics of cytokine transcriptome, chemokine transcriptome and known signatures in distinct gene clusters. Figure S6. The prognostic value of m6Ascore and correlation between the clinicopathological features and m6Ascore. Figure S7. The effect of microsatellite status and EB virus infection on the three m6A modification patterns and m6A regulators. Figure S8. Prognostic value of m6Ascore in gastric cancer cohorts and digestive cancer cohorts.
Table S1. Basic information of datasets included in this study for identifying distinct m6A methylation modification patterns. Table S2. The gene sets used in this work for marking each TME infiltration cell type. Table S3. Spearman correlation analysis of the 21 m6A modification regulators. Table S4. Estimating relative abundance of tumor microenvironment cells in 1059 gastric cancer patients by the Single-Sample Gene-Set. Table S5. The activation states of biological pathways in distinct m6A modification patterns by GSVA enrichment analysis. Table S6. The changes of m6Aclusters, ACRG molecular subtypes, gene clusters and m6Ascore. Table S7. Prognostic analysis of 718 m6A phenotype-related genes using a univariate Cox regression model. Table S8. Functional annotation for m6A phenotype -related genes (Gene Ontology-Biological process). Table S9. Spearman correlation between m6Ascore and other known signatures within the gastric cancer.
About this article
Cite this article
Zhang, B., Wu, Q., Li, B. et al. m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer 19, 53 (2020). https://doi.org/10.1186/s12943-020-01170-0