Crosstalk of Histone and RNA Modifications Identified a Stromal-Activated Subtype with Poor Survival and Resistance to Immunotherapy in Gastric Cancer

Emerging evidence has revealed the pivotal role of epigenetic modifications in shaping the tumor microenvironment (TME). However, crosstalk between different modification types and their clinical relevance in cancers remain largely unexplored. In this study, using ChIP/MeRIP-seq data of seven human gastric cell lines, we systematically characterized the crosstalk of four epigenetic modification types including H3K4me1, H3K4me3, H3K27ac, and N6-methyladenosine (m6A) and identified a recurrent subtype with high FTO expression and low HDAC1 expression across three independent gastric cancer (GC) cohorts, which we named the epigenetic-modification-dysregulated (EMD) subtype. Patients of the EMD subtype were featured with poor survival, stromal activation, and immune suppression. Extensive relevance to clinical characteristics was observed in the EMD subtype, including the Lauren classification, MSI status, histological grade, TNM stage, the Asian Cancer Research Group classification, and the immune/fibrotic classification. An EMD score was then constructed using WGCNA and ssGSEA algorithms, to precisely recognize the EMD subtype and indicate prognosis and response to immunotherapy in multiple independent GC cohorts. Correlations of the EMD score with tumor mutation burden, tumor purity, aneuploidy score, tumorigenic pathways, TME characteristics, and FTO/HDAC1 ratio were measured. In vitro experiments were performed to demonstrate the correlation between FTO and the epithelial–mesenchymal transition pathway, which suggested FTO as a targetable vulnerability for GC patients with a high EMD score. Altogether, by comprehensively analyzing the epigenetic modification patterns of 1518 GC patients, we identified a novel stromal-activated subtype with poor survival and resistance to immunotherapy, which might benefit from the combined immune checkpoint inhibition therapy with FTO inhibition.

1 Introduction

Gastric cancer (GC), the fifth most common cancer and the third most common cause of cancer deaths in the world, remains a non-negligible health problem and social burden globally (Smyth et al., 2020). Patients with early-stage GC featured a 5-year overall survival (OS) rate of more than 60%, for whom surgical resection is the best option (Thrift and El-Serag, 2020). For patients with advanced GC, chemotherapy represented by fluoropyrimidines and platinum significantly improves survival and the quality of life. Unfortunately, due to the high frequency of advanced stage at diagnosis and chemotherapy resistance, the 5-year overall survival rate of advanced GC is still less than 5% (Yuan et al., 2020). As a complex disease with molecular and clinical heterogeneity, GC shows versatile phenotypes in initiation, progression, and even the response to treatment among different patients. Thus, more effective individualized strategies both in diagnosis and therapy for GC are urgently needed to be explored.

Epigenetic modification, mainly including DNA methylation, histone modification, and RNA modification, is a significant regulatory mechanism of diverse physiological or pathological processes (Zhao et al., 2021a). Histone modification mainly refers to the post-translational modifications (PTMs) that occur in the N-terminal tails of histone proteins of nucleosomes, including but not limited to methylation, acetylation, and ubiquitination (Zhao and Shilatifard, 2019). Generally, histone modifications (such as H3K4me1, H3K4me3, and H3K27ac) enriched at the enhancer or promoter region presumably facilitate the transcription process of the targeted genes (Bannister and Kouzarides, 2011). N6-methyladenosine (m6A), defined as methylation of adenosine at the N6 position, is one of the most abundant RNA modification types in eukaryotic species including mammals, plants, insects, yeast, and certain viruses (Gu et al., 2020). Recently, accumulating evidence has suggested the extensive interactions between histone and m6A modifications, which trigger epigenetic remodeling and cause profound impacts on various aspects of cancer progression, including the resistance to medical treatment. (Huang et al., 2019; Li et al., 2020a; Yang et al., 2021a; Li et al., 2021). For example, Li et al. uncovered a SOX4/EZH2/METTL3 axis in TMZ-resistant glioblastoma (GBM), in which EZH2 regulates the METTL3 expression via an H3K27me3 modification-independent manner, and METTL3 leads to nonsense-mediated mRNA decay of EZH2 reversely (Li et al., 2021). Moreover, Li et al. found that the m6A reader YTHDC1 physically interacts with and recruits KDM3B to m6A-associated chromatin regions, promoting H3K9me2 demethylation and gene expression, establishing a direct link between m6A and histone modification (Li et al., 2020a). Similar interactions have also been observed in GC. Yang et al.’s study suggested that HDAC3 regulates the FTO (fat-mass and obesity-associated protein) expression in a FOXA2-dependent manner, thus promoting the proliferation, migration, and invasion of GC cells (Yang et al., 2021a).

The tumor microenvironment (TME), known as the soil of the tumor seed, which plays a pivotal role in tumorigenesis and anti-tumor immunity, has been recently reported to be shaped by various epigenetic modifications (Huang et al., 2012; Gu et al., 2021). For instance, Yin et al. reported that EZH2 depletion increased generation of the IL-15 receptor (IL-15R), CD122 (+) NK precursors, and mature NK progeny from both mouse and human hematopoietic stem and progenitor cells, demonstrating the impact of histone modification H3K27me3 on early NK cell differentiation (Yin et al., 2015). Wang et al. found that METTL3/14-deficient tumors increased cytotoxic tumor-infiltrating CD8⁺ T cells and elevated secretion of IFN-γ, CXCL9, and CXCL10 in the TME in vivo, and inhibition of METTL3/14 enhanced response to anti-PD-1 treatment in pMMR-MSI-L CRC and melanoma (Wang et al., 2020). However, the epigenetic modification patterns and their association with TME, prognosis, and therapeutic response in GC remain poorly investigated.

In this study, we aimed to characterize the extensive crosstalk between histone and m6A modifications in GC, trying to explain the molecular and clinical heterogeneity of GC from the perspective of epigenetic dysregulations.

2 Materials and Methods

2.1 Data Acquisition and Processing

As shown in the flowchart of this study (Figure 1), we obtained ChIP-seq data of three histone modification types (H3K4me1, H3K4me3, and H3K27ac) in four gastric cell lines (GES-1, SNU719, NCC24, and YCC10), and MeRIP-seq data of m6A in three GC cell lines (AGS, BGC823, and SGC7901) from previously published literature reports with the access number of GSE135175 (Okabe et al., 2020), GSE166972 (Chen et al., 2021a), GSE133132 (Yue et al., 2019), and PRJNA595769 in the Gene Expression Omnibus (or the European Nucleotide Archive), respectively (Supplementary Table S1). For ChIP-seq, raw sequencing reads were aligned using Hisat2 (Kim et al., 2019) with default parameters to the hg19. Sambamba (Tarasov et al., 2015) was used to remove PCR duplicates and obtain the uniquely mapped reads. For MeRIP-seq, raw sequencing reads were aligned using Hisat2 with default parameters to the hg38. Samtools (Li et al., 2009) was used to remove the reads with a mapping quality below 30.

FIGURE 1

FIGURE 1. Flow chart of this study. Three types of data from 14 datasets were used in this study: datasets in red refer to the epigenetic data (ChIPseq or MeRIPseq) targeting on histone or m6A modifications, datasets in green refer to the RNAseq or array data with complete clinical information (survival or ICI response), and datasets in yellow refer to the RNAseq data of cell lines with or without FTO depletion. The main method used in each step was also listed.

Totally 1518 GC patients from eight independent cohorts with complete transcriptomic data and clinical information were obtained from The Cancer Genome Atlas TCGA-STAD (stomach adenocarcinoma, https://portal.gdc.cancer.gov), Gene Expression Omnibus [GSE62254 (Cristescu et al., 2015), GSE15459 (Ooi et al., 2009), GSE13861 (Cho et al., 2011), GSE26899 (Oh et al., 2018), GSE26901 (Oh et al., 2018), GSE84437, https://www.ncbi.nlm.nih.gov/gds], and European Nucleotide Archive (PRJEB25780, https://www.ebi.ac.uk) (Supplementary Table S1). Cohorts were divided into exploring set (TCGA-STAD, GSE15459, and GSE62254) and validation set (GSE13861, GSE26899, GSE26901, GSE84437, and PRJEB25780) (Figure 1). Patients with prior or synchronous malignancy diagnoses or patients who survived less than 30 days were excluded from this study. Sva (Jeffrey, 2020) was used to correct the non-biological batch effects among different cohorts.

2.2 Collection of Epigenetic Regulators

Epigenetic regulators, defined as the key genes that play crucial roles in the regulation or function of epigenetic modifications, mainly include the methyltransferases/demethylases of histone methylation (H3K4me1 and K3K4me3), acetylases/deacetylases of histone acetylation (H3K27ac), and writers/readers/erasers of m6A modification. In our study, totally 43 epigenetic regulators (Supplementary Table S1), including histone methyltransferases (SETD1A, SETD1B, KMT2A, KMT2B, KMT2C, and KMT2D), histone demethylases (KDM1A, KDM1B, KDM2A, KDM2B, KDM5A, KDM5B, KDM5C, and KDM5D), histone acetylases (CREBBP and EP300), histone deacetylases (HDAC1, HDAC2, HDAC3, HDAC8, and HDAC11), m6A writers (METTL14, METTL3, METTL16, RBM15, RBM15B, ZC3H13, ZCCHC4, WTAP, CBLL1, and VIRMA), m6A readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, and HNRNPC), and m6A erasers (ALKBH5 and FTO), were collected from the published literature (Zhao and Shilatifard, 2019; Zhao et al., 2021a; Deng et al., 2021) for our study. All the included regulators have been experimentally demonstrated to regulate one or more modification types, most of which were reported to participate in tumorigenesis or progression of GC (Yue et al., 2019; Wu et al., 2021a). Maftools (Mayakonda et al., 2018) was used to identify the mutation (Figure 6F) or co-mutation (Figure 2A) events of epigenetic regulators in the TCGA-STAD cohort. Co-occurrence events with p values less than 0.05 were defined as the co-mutation events (Supplementary Table S2).

FIGURE 2

FIGURE 2. Crosstalk between regulators of histone and m6A modifications. (A) Co-mutation heatmap of epigenetic regulators in the TCGA-STAD cohort. (B) Co-expression heatmap of epigenetic regulators in the TCGA-STAD cohort. (C) Differential H3K4me1 and H3K27ac modifications in the m6A regulator ZCCHC4. Peaks in red and green refer to the H3K4me1 and H3K27ac modifications in GC cell lines, respectively, and gray refers to the corresponding modification in GES-1 as control. The red line at the bottom refers to the promoter region. (D) Differential m6A modification in the histone modification regulator KMT2D. Peaks in blue and orange refer to the immunoprecipitation (IP) signals of AGS and BGC823 cell lines, respectively. Peaks in gray refer to the corresponding input signals. (E) Differential expression of ZCCHC4 between GC cell lines and GES-1 (p < 0.05). (F) Correlations of KMT2D with the m6A regulators including RBM15, RBM15B, and YTHDF1 (Spearman, p < 0.05).

2.3 Identification of the Epigenetic Modification-Related Gene

For histone modifications (H3K4me1, H3K4me3, and H3K27ac), DiffBind (Stark and Brown, 2011) was used to identify the differentially-histone-modfied site (DHMS) between GC cell lines and GES-1 according to the criteria of |fold-change| >1 and p < 0.05. BigWig data were downloaded from the GEO with the access number GSE135175. DESeq2 (Love et al., 2014) was used to identify the differentially expressed gene (DEG) between GC cell lines and GES-1 with the criteria of |log₂FC| >1 and p < 0.05. For m6A modification, exomePeak2 (Wei, 2020) was used to identify the m6A-modified site (MMS) in each GC cell line, according to the criteria of |fold-change|>1 and p < 0.05. The Wilcoxon test was used to identify the correlatively expressed gene (CEG) of the 22 m6A regulators using the criteria of |r|>0.3 and adjusted p < 0.05.

Generally, histone methylations/acetylations (such as H3K4me1, H3K4me3, and H3K27ac) enriched at enhancers or promoters, as well as m6A modifications enriched at coding sequence (CDS) or the 3′ untranslated region (3′ UTR) region, presumably facilitate the expression of target genes (Bannister and Kouzarides, 2011; Shi et al., 2019). Thus, those genes which satisfy one of the following criteria were defined as the epigenetic-modification-related gene (EMRG):

1) At least one DHMS was located in the promoter of the DEG, and the expression of the DEG changed in the same direction as the DHMS in at least one GC cell line (SNU719, NCC24, and YCC10);

2) At least one MMS was located in the CEG in at least one GC cell line (AGS, BGC823, and SGC7901), and the CEG positively correlated (r > 0.3, Wilcoxon test) with at least one of the m6A writers/readers, or negatively correlated (r < −0.3, Wilcoxon test) with at least one of the m6A erasers.

The promoter region was defined as 3 kb upstream and downstream of the transcriptional start site (TSS) in each gene. The correlations between the CEG and m6A regulators were calculated in the TCGA-STAD cohort. Long non-coding RNAs (LncRNAs) of the epigenetic modification-related Gene were defined as the epigenetic-modification-related LncRNA (EMRL), according to the GENCODE annotations (Yang et al., 2021b).

2.4 Unsupervised Clustering in the Exploring Set

By performing univariate Cox hazard analysis in the TCGA-STAD cohort, we selected 34 survival-associated EMRLs (p < 0.01) for further study. Based on the 34 survival-associated EMRLs, ConsensusClusterPlus (Wilkerson and Hayes, 2010) was applied to perform unsupervised clustering in the exploring set with the following parameters: maxK = 7, cluster algorithm = km, and correlation method = Euclidean.

2.5 Construction of the EMD Score

We first identified 1,674 DEGs between the epigenetic-modification-dysregulated (EMD) subtype (C4) and cluster C3 using limma (Ritchie et al., 2015). (Supplementary Table S7) Weighted gene co-expression network analysis (WGCNA) (Langfelder and Horvath, 2008) was then performed based on the DEGs to recognize the most positively correlated (Meturquois, Spearman-r = 0.75) and negatively correlated (MEblue, Spearman-r = −0.61) gene modules relating to cluster C4 (Figures 6A–C). We further chose 147 (set1) and 64 (set2) genes from the turquoise and blue modules, respectively, according to their correlation coefficients with the EMD subtype and corresponding modules (Figure 6D; Supplementary Table S8). The single sample gene set enrichment analysis (ssGSEA) algorithm in R package GSVA (Hänzelmann et al., 2013) was used to estimate the relative abundance (ssGSEA value) of set1 and set2, respectively, for each patient in the GC cohorts. The EMD score was defined as the ratio of the two ssGSEA values. Log transition of the ratio was also performed to make the score less discrete. For each GC patient, the EMD score was calculated as follows:

$EMD score={\log }_2\left\{ ssGSEA score \left(set1\right)/ ssGSEA score \left(set2\right)\right\}.$

2.6 TME Characterization and Function Enrichment Analysis

CIBERSORT (Chen et al., 2018) and ESTIMATE (Yoshihara et al., 2013) were used to evaluate the relative abundance of immune infiltration and stromal components in tumor samples, respectively. Deconvolution results of CIBERSORT were evaluated by a derived p-value (p < 0.05) to filter out the samples with less significant accuracy. The tumor immune dysfunction and exclusion (TIDE) algorithm (Jiang et al., 2018) was used to measure the antitumor immunity features of each patient. ClusterProfiler (Yu et al., 2012) was used to perform gene set enrichment analysis (GSEA), using the hallmark gene sets downloaded from the molecular signatures database (MSigDB).

2.7 Cell Culture and Treatment

Human GC cell lines (SNU719 and SGC7901) were purchased from American Type Culture Collection (Manassas, United States). SNU719 and SGC7901 cells were cultured in RPMI-1640 medium (Gibco) with 10% fetal bovine serum (Gibco). Cells were maintained at 37°C and 5% CO₂. SNU719 and SGC7901 cells were plated onto 6-well plates and reached 70–80% cell confluence on the day of treatment. Cells were divided into the BS group and the control group. The BS group was treated with the FTO inhibitor (Brequinar sodium, V17016, 5 umol/l, InvivoChem), and the DMSO group was treated with an equivalent DMSO concentration as control for 48 h.

2.8 Western Blot

48 h after treatment, the EMT markers in SNU719 and SGC7901 were analyzed by Western blot. In brief, whole protein samples of all groups were extracted, and the concentration was detected by using the Pierce BCA Protein Assay Kit (Thermo 23225). Equal amounts of protein (20 ug per well) were loaded, and the samples were separated by 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis. Then, the proteins were transferred to 0.45 μm polyvinylidene fluoride (PVDF) membranes. The membranes were blocked with 5% nonfat milk and incubated with primary antibodies at 4°C overnight. Later, the membranes were washed with TBST three times and incubated with the secondary HRP antibody at room temperature for 2 h. The primary antibodies used in this research were as follows: E-cadherin (20874-1-AP, 1:2000, Proteintech), N-cadherin (22018-1-AP, 1:2000, Proteintech), Vimentin (10366-1-AP, 1:5,000, Proteintech), and TWIST1 (25465-1-AP, 1:1,000, Proteintech). β-actin (3700T, 1:5,000, CST) was used as a loading control.

2.9 Statistical Analysis

All statistical analyses were performed using R (v4.0.2) (Team, R. C 2020) and its appropriate packages. p-values <0.05 were regarded as statistically significant.

3 Results

3.1 Crosstalk between Regulators of Histone and m6A Modifications

First, we examined the interactions of the four modification types in the TCGA-STAD cohort. Totally, 211 co-mutation events (p < 0.05) and 266 co-expression events (|r|>0.3, p < 0.05) were observed in the 43 epigenetic regulators (Figures 2A, B; Supplementary Table S2). These correlations were further illustrated by the ChIP/MeRIP-seq data from GC cell lines. Specifically, m6A regulators including ZCCHC4, IGF2BP3, and VIRMA were differentially expressed between the GC cell lines and GES-1, which were accompanied by differential histone modifications in their TSS regions (Figure 2; Supplementary Figure S1). Similarly, the histone modification regulators, SETD1A and KMT2D, which showed positive correlations with the m6A regulators (RBM15, RMM15B, and YTHDF1), had differential m6A modification in their exons (Figure 2; Supplementary Figure S1; Supplementary Table S2). These findings revealed the active crosstalk between histone and m6A modifications, both in genome and transcriptome levels.

3.2 Identification of the EMRG and EMRL

To generally characterize the transcriptome landscape altered by the four epigenetic modification types, we then identified the dysregulated genes associated with each of the modification types. Using the criteria described earlier, we finally identified 4,999 EMRGs of H3K4me1, 3,693 EMRGs of H3K4me3, 4,095 EMRGs of H3K27ac, and 6,117 EMRGs of m6A. Previous studies suggested that some histone modifications (such as H3K4me1, H3K4me3, and H3K27ac) enriched at the enhancer or promoter region presumably facilitate the transcription of targeted genes (Bannister and Kouzarides, 2011). Similarly, unbalanced distribution was also observed in the m6A modification, with more than 40% of all modification sites in mRNA being present in 3′ UTRs) (Ke et al., 2015; Kan et al., 2022). Consistent with the previous studies, most of MMSs were located at the 3′UTR, followed by exons, and then the promoter region (Figure 3A). Meanwhile, for the three histone modification types, most of DHMSs were located within the promoter (<1 kb), followed by the promoter (1–2 kb) and promoter (2–3 kb) (Figure 3A). Among the EMRGs, nearly 50% of genes were simultaneously regulated by more than one modification type, which was another evidence for the crosstalk between different modification types. Specifically, about 360 protein-coding genes (PCGs) were co-regulated by the four modification types (Figure 3B).

FIGURE 3

FIGURE 3. Identification of epigenetic modification-related LncRNA (EMRL). (A) Doughnut charts showed the distribution of DHMS or MMS in each modification type. (B)Venn diagram showed the overlaps of EMRG (or EMRL) regulated by different modification types. PCG refers to the protein-coding genes, and LncRNA refers to the long non-coding RNAs. (C) MNX1-AS1 was potentially regulated by H3K4me3 and H3K27ac in GC cell lines. Peaks in yellow and green refer to H3K4me3 and H3K27ac, respectively, while gray ones refer to the corresponding modification in GES-1 as control. The red line at the bottom refers to the promoter region. (D) m6A modification regions of MNX1-AS1 in the AGS cell line. Peaks in blue refer to the immunoprecipitation (IP) signals. Peaks in gray refer to the input signals. (E) Differential expression of MNX1-AS1 in GC cell lines and GES-1 (p < 0.05). (F) Correlations of MNX1-AS1 with the m6A regulators including HNRNPC, HNRNPA2B1, and IGF2BP2 (Spearman, p < 0.05).

As another key member of epigenetic regulation, LncRNA was demonstrated to participate in various tumorigenic processes due to their extensive biological functions in the transcriptome level (Sempere et al., 2021; Kan et al., 2021). To further characterize the epigenetic modification patterns with clinical relevance in GC, we also identified 370 EMRLs (Supplementary Table S3) from the EMRG. Interestingly, three EMRLs (PAXIP1-AS2, MNX1-AS1, and PVT1) were co-regulated by the four modification types (Figure 3B). LINC00511, previously reported as an oncogene in several solid tumors (Wu et al., 2020; Peng et al., 2021; Dong et al., 2021), was significantly modified with H3K4me1 (Supplementary Figure S2) and simultaneously overexpressed in the three GC cell lines (Figure 3F). Similarly, concomitant histone modification and overexpression were observed in LINC01091 (Supplementary Figure S2). Moreover, we found MNX1-AS1, another oncogene in multiple cancers (Wu et al., 2019; Li et al., 2020b; Li et al., 2020c; Wu et al., 2021b), was simultaneously regulated by both histone and m6A modifications, with differential expression across three GC cell lines and positive correlations to multiple m6A regulators (Figures 3C–F). The aforementioned findings shed light on the close connections inside epigenetic regulators including epigenetic modifications and LncRNA, offering a new perspective for exploring the complicated regulation network of the epigenome.

3.3 Unsupervised Clustering Based on EMRL Identified an EMD Subtype of GC

Next, we set out to characterize the modification patterns in the exploring set. Through univariate Cox regression analysis, we selected 34 EMRLs with a significant prognosis value (p < 0.01) for further study (Supplementary Table S4). Unsupervised clustering was then performed based on the 34 survival-associated EMRLs, which divided the exploring set into four distinct clusters (Figures 4A, B; Supplementary Figure S3A; Supplementary Table S5). Surprisingly, extremely similar expression patterns were observed across the three independent GC cohorts (TCGA-STAD, GSE15459, and GSE62254), especially in cluster C4. Specifically, the C4 population of each cohort was characterized by strikingly elevated EMRLs including ZEB1-AS1, NR2F1-AS1, MIR100HG, ZFHX4-AS1, and PART1 and the suppression of EMRLs including HOTTIP, HOXA11-AS, CASC2, and LINC00467. Moreover, a significant difference in prognosis was observed among the four clusters, with cluster C4 having the worst survival and cluster C3 the best (Figure 4C).

FIGURE 4

FIGURE 4. Distinctive epigenetic modification patterns across three GC cohorts recognized by EMRL. (A) Heatmap showed the highly concordant expression patterns of the three GC cohorts in the exploring set. (B) Consensus unsupervised clustering divided the exploring set into four distinct clusters. (C) Kaplan–Meier curves showed the different overall survival of the four clusters. (D) Heatmap showed the distinct expression patterns of the epigenetic regulators in the exploring set. VIRMA and METTL16 were excluded in this section due to lack of expression information in the GSE15459 or GSE62254 cohort.

Consistently, unbalanced modification patterns both of histone and m6A modifications were observed in cluster C4. Specifically, cluster C4 was characterized by the high expression of multiple histone modification writers (SETD1B, KMT2A, and CREBBP) and the low expression of histone modification erasers (KDM1A, KDM1B, HDAC1, and HDAC2). While in the m6A modification, high-expressed erasers (FTO and ALKBH5) and low-expressed readers (YTHDF2, YTHDF3, IGF2BP2, IGF2BP3, and HNRNPC) were observed in cluster C4 (Figure 4D). Then, we further explored the mutation profiling of cluster C4. Interestingly, cluster C4 was also featured with a distinctive mutation pattern, with rarely detected mutations of the epigenetic regulators (Supplementary Figure S3B). Given the distinctive modification patterns in transcriptome and genome levels and the poor prognosis of cluster C4, we named it the epigenetic-modification-dysregulated (EMD) subtype in further study. In this section, using survival-associated EMRLs, we identified a conserved EMD subtype with distinctive modification patterns in GC.

3.4 Molecular and Clinical Characterization of the EMD Subtype

To figure out the mechanism underlying the poor prognosis of the EMD subtype, we comprehensively explored the molecular and clinical characteristics of the EMD subtype by comparing it with cluster C3 (cluster with the best prognosis) in the exploring set. In the GSEA analysis, multiple tumorigenic pathways including the epithelial–mesenchymal transition (EMT, p < 0.001), Hedgehog (p = 0.002), TGF-β (p = 0.01), IL2-STAT5 (p = 0.001), KRAS (p < 0.001), and angiogenesis (p < 0.001) were significantly activated in the EMD subtype (Figure 5A), whereas pathways including DNA repair (p < 0.001) and G2M checkpoint (p < 0.001) were observed to be significantly suppressed in the EMD subtype (Figure 5A).

FIGURE 5

FIGURE 5. Molecular and clinical characteristics of four epigenetic modification patterns. (A) GSEA analysis identified the differentially activated pathways of cluster C4 compared with cluster C3. (B) CIBERSORT analysis estimated the relative ratio of tumor-infiltrating cells in clusters C3 and C4. Each value was defined as the relative ratio compared to the 22 human hematopoietic cell types contained in CIBERSORT. (C) Different TMB, tumor purity, and aneuploidy score among four clusters in the TCGA-STAD cohort. (D) Clinical relevance of the EMD subtype in the TCGA-STAD cohort. (E) Clinical relevance of the EMD subtype in the GSE62254 cohort (*p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001).

Then, we performed CIBERSORT analysis to explore the TME characteristic of the EMD subtype. Compared with cluster C3, the EMD subtype showed a significantly higher ratio of macrophages M2 (p < 0.001) and resting CD4⁺ T cell (p < 0.01) and a significantly lower ratio of macrophages M1 (p < 0.05) and activated CD4⁺ T cell (p < 0.001) (Figure 5B). Moreover, a relatively lower CD4⁺/CD8⁺ T cell ratio was observed in the EMD subtype (p < 0.05, Supplementary Figure S4A), which usually indicates more advanced tumors or immune deficiency in cancer patients. In the TCGA-STAD cohort, we also observed a significant decrease of tumor mutation burden (TMB, p < 0.001), tumor purity (p < 0.001), and aneuploidy score (AS, p < 0.05) in the EMD subtype (Figure 5C).

In addition to the distinctive molecular characteristics, the EMD subtype also showed extensive correlations to various clinical characteristics of GC (Supplementary Table S6). In the TCGA-STAD cohort (Figure 5D), we found the EMD subtype as a good indicator for the IE/F subtype of the newly reported immune/fibrotic classification (Bagaev et al., 2021) (p < 0.001). Besides, more patients with MSS subtype (p < 0.001), poor differentiation (p = 0.002), or diffuse histotype (p < 0.001) were observed in the EMD subtype. While in the GSE62254 cohort (Figure 5E), the EMD subtype mostly overlapped with the MSS/EMT subtype of ACRG classification (Cristescu et al., 2015) (p < 0.001). Similarly, more patients with advanced TNM stage (p < 0.001), perineural invasion (p < 0.001), or diffuse histotype (p < 0.001) were observed in the EMD subtype.

These findings suggested distinctive stromal-activated and immune-suppressed characteristics of the EMD subtype, which may account for the poor prognosis of this subtype.

3.5 Construction of the EMD Score in the Exploring Set

Since clusters C3 and C4 showed marked differences in prognosis and molecular characteristics, we chose clusters C3 and C4 for further study. To expand the applicability of the EMD subtype, an EMD score was then constructed based on the DEGs between the EMD subtype (C4) and cluster C3 using WGCNA (Figures 6A–D) and ssGSEA algorithms. Strikingly, the EMD score was closely related to the prognosis in the exploring set, where the increasing EMD score was accompanied by more death events (Figure 6G) and worse survival (Figure 6E).

FIGURE 6

FIGURE 6. Construction of the EMD score in the exploring set. (A–D) WGCNA analysis based on 1674 DEGs in the exploring set. (E) Kaplan–Meier curves showed the different overall survival between the high- and low-EMD score groups in the exploring set. Patients were separated into the high and low groups according to the median EMD score of each cohort. (F) Mutation landscape of epigenetic regulators in high- and low-EMD score groups in the TCGA-STAD cohort. (G) Risk curves showed more frequent death events with the increasing EMD score in the exploring set. Dots in yellow referred to death events observed in the cohorts. (H) Heatmap showed the close correlations of the EMD score with various molecular characteristics including epigenetic regulators, EMRLs, TME markers, oncogenes, and tumorigenic pathways. All the values were scaled using the Z-score.

Consistent with the EMD subtype, extensive correlations of the EMD score with stromal-activated and immune-suppressed TME characteristics could be observed in the exploring set (Figure 6H; Supplementary Tables S9, S10). Specifically, patients with a high EMD score showed activation of multiple pathways (EMT, TGF-β, and Hedgehog), elevated stromal scores, high abundance of cancer-associated fibroblast (CAF) and mast cells, whereas those with a low EMD score showed low MHC-I expression, low infiltration of activated CD4⁺ T cells, and low M1/M2 ratio. These were further evidenced by upregulated markers of EMT activation (ZEB2, TWIST1, SNAI2, MMP2, CDH2, VIM, etc.), CAF (ACAT2, DES, TNC, PDGFRB, etc.), and immune-suppressed chemokine axis (CXCR4, CXCL12, CCL2, and CCR2) with the increasing EMD score. Moreover, we also observed the overexpression of multiple immune evasion-related genes (CREBBP, HCFC2, TGFBR2, WDR7, ATG14, and PPP2R3C), which were experimentally demonstrated (Lawson et al., 2020), in the patients with a high EMD score. Additionally, in the TIDE analysis, we found the EMD score significantly correlated with the exclusion, dysfunction, CAF, and MSI signature, which implied the potential of the EMD score in predicting the ICI efficacy. Similar to the EMD subtype, a high EMD score indicated a lower frequency of regulator mutation (34%, 51%, Figure 6F). These findings suggested the EMD score as a good indicator for recognizing the EMD subtype in GC.

3.6 Recognizing Patients with EMD Characteristics in the Validation Set Using the EMD Score

We first tested the correlations of the EMD score with the four clusters we identified earlier, ACRG classification, and immune/fibrotic classification. Strikingly, the EMD score showed an excellent performance in recognizing the EMD subtype (AUC = 0.96, Figure 7A) in the exploring set, as well as the MMS/EMT subtype (AUC = 0.95, Figure 7A) of ACRG classification and the IE/F subtype of the immune/fibrotic classification (AUC = 0.86, Figure 7A).

FIGURE 7

FIGURE 7. Recognizing patients with EMD characteristics using the EMD score. (A) Differential EMD score of subtypes in three classification proposals. The ROC below showed the performance of the EMD score in recognizing the EMD, MSS/EMT, and IE/F subtype in the corresponding cohort. (B) Hallmark pathways correlated with the EMD score in the validation set. (C) Kaplan–Meier curves showed the different overall survival between the high and low EMD groups in the validation set. Patients were separated into high and low groups according to the median EMD score of each cohort. (D) Correlations of the EMD score with the ICI outcome in the cohort PRJEB25780. Patients were divided into high or low groups according to the lower quartile (Q1) of the EMD scores (right panel). (E) GSEA analysis suggested significant activated EMT and Hedgehog pathways in non-responders compared with the responder.

To further validate the ability of the EMD score in predicting prognosis and response to ICI treatment, we calculated the EMD score for patients in four independent GC cohorts. As we expected, a high EMD score was closely related to a worse OS (p < 0.05, Figure 7C) in each cohort of the validation set. Strong correlations between the EMD score and multiple pathways (EMT, TGF-β, Hedgehog, etc. Figure 7B; Supplementary Table S11) could also be observed in the validation set. In the cohort PRJEB25780, we also observed a distinct difference in the response rate to immune checkpoint inhibitor (ICI) treatment between high and low EMD score groups (10 vs. 75%, p = 0.0018, Figure 7D). Moreover, non-responders were featured with significant activated EMT (NES = 1.96, p < 0.001) and Hedgehog (NES = 1.85, p < 0.001) pathways compared with the responders (Figure 7E). These findings successfully validated the ability of the EMD score in recognizing the patients with EMD characteristics and the potential of recognizing patients sensitive to ICI treatment in GC.

3.7 FTO as a Potential Target for Patients with High EMD Score

Fat mass- and obesity-associated protein (FTO) was recently reported to participate in multiple tumorigenic processes including immune evasion in several malignant tumors (Li et al., 2017; Zhao et al., 2020; Tao et al., 2021; Zhou et al., 2021), implying the therapeutic potential of FTO inhibition in treating multiple cancers. In this study, high FTO expression and low HDAC1 expression were both observed in patients with high EMD score (Figure 8A). Moreover, FTO and HDAC1 were the most positively correlated (Spearman-r = 0.48, p < 0.001) and negatively correlated (Spearman-r = −0.45, p < 0.001) regulators relating to the EMD score, which also represented the m6A and histone modifications, respectively. To characterize the EMD subtype with both histone and m6A modification features, we measured the FTO/HDAC1 ratio for each GC patient. Interestingly, we found the FTO/HDAC1 ratio closely correlated with the EMD subtype (Supplementary Figure S4B), as well as the EMT (r = 0.52, p < 0.001, Figures 8B, C) and TGF-β (r = 0.34, p < 0.001, Figures 8B, C) pathway activation, which indicated stromal-activated TME characteristics. Also, a high FTO/HDAC1 ratio was accompanied by a low CD4⁺ T cell activation rate and a low M1/M2 ratio (Figure 8D), which indicated immune-suppressed TME characteristics. Moreover, a high FTO/HDAC1 ratio was related to a worse OS in the whole cohort of this study (Figure 8E). These findings revealed the close links of FTO with GC prognosis and multiple tumorigenic pathways.

FIGURE 8

FIGURE 8. FTO may serve as a potential target for patients with a high EMD score. (A) Correlations of the EMD score with FTO and HDAC1 in each cohort of the validation set. (B) Correlations of the FTO/HDAC1 ratio with multiple pathways in the validation set. (C) GSEA analysis confirmed the activation of EMT and Hedgehog pathways in patients with a high FTO/HDAC1 ratio. (D) Correlations of the FTO/HDAC1 ratio with the activated/resting CD4⁺ T cell ratio and M1/M2 ratio in the validation set. Activated/resting CD4⁺ T cell ratio refers to the ratio of activated CD4⁺ T cells and resting CD4⁺ T cells. The M1/M2 ratio refers to the ratio of macrophage M1 and macrophage M2. (E) Kaplan–Meier curves showed the different overall survival between the high and low FTO/HDAC1 ratio group in the whole cohort of this study. Patients were separated into high and low groups according to the median EMD score of the whole cohort. (F) Barplots showed the fold change of EMT and Hedgehog markers under FTO depletion in the AGS cell line (GSE178697). (G) Western blot showed the suppression of the EMT pathway under FTO inhibition in the SNU719 cell line. (H) Western blot showed the suppression of the EMT pathway under FTO inhibition in the SGC7901 cell line (*p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001).

These links were further evidenced by the significant downregulation of activated EMT markers (ZEB2, TWIST1, SNAI2, MMP2, and CDH2) in a GC cell line with FTO knockdown (FTO^−KD) (Figure 8F; Supplementary Table S12). Consistently, EMT (p < 0.001) and TGF-β (p = 0.007) suppressions were also observed in an FTO^−KD colon cancer (CC) cell line (Supplementary Figures S3C–E; Supplementary Table S12). Targeting FTO, as we assumed, might help to improve the immune-suppressed TME caused by the activated EMT and TGF-β pathway.

Brequinar sodium (BS) was reported to bind tightly with FTO protein and inhibit the demethylase activity of FTO (Su et al., 2020). To confirm the FTO’s function of regulating the EMT pathway, we conducted in vitro experiments using two GC cell lines, SNU719 and SGC7901. Strikingly, when treated with BS, markers of EMT activation including N-cadherin, Vimentin, and TWIST1 were downregulated in both cell lines (Figures 8G, H). Although no significant upregulation of E-cadherin (the epithelial marker) was observed in the BS-treated cell lines, we still noticed a remarkable increase of the E-cadherin/N-cadherin ratio after BS treatment, especially in SGC7901 cell line (p < 0.05, not presented). These findings suggested FTO as a promising target for inhibiting the EMT pathway in GC, especially for patients with high EMD score.

4 Discussion

Extensive interactions between epigenetic modification types in shaping the TME of cancers have drawn increasing attention recently. Chen et al. systematically characterized the interactions of four RNA modification types (m6A, N1-methyladenosine, alternative polyadenylation, and adenosine-to-inosine RNA editing) and demonstrated that multilayer alterations of RNA modification “writer” are associated with patient survival and TME cell-infiltrating characteristics (Chen et al., 2021b). Zhao et al. focused on the relationship between LncRNAs and histone/DNA modifications and identified key LncRNA regulators as a prognostic biomarker for breast cancer subtypes (Zhao et al., 2021b). Similarly, regarding GC, Meng et al. (2021) elucidated the interactions between DNA methylation regulators and generated a DMS score for separating GC patients with distinctive prognosis and treatment efficacy (Meng et al., 2021). These studies provided insights into the active interaction networks of epigenetic modifications, suggesting their pivotal roles in shaping the TME across multiple tumors. However, the interactions between histone and RNA modifications and their impact on the TME in GC had not been fully explored yet.

In this study, for the first time, we comprehensively depicted the epigenetic regulation network including histone modification, RNA modification, and LncRNA in GC. Both in genome and transcriptome levels, we characterized the crosstalk between regulators of these four modification types. Zhan et al. found that several important histone regulator genes, including KMT2D, KMT2C, CREBBP, and EP300, are frequently altered in esophageal squamous cell carcinoma (Song et al., 2014). In our study, we further observed frequent co-mutations of these regulators, such as KMT2D-CREBBP, KMT2C-EP300, and CREBBP-EP300 in GC, which implied potential vulnerabilities based on epigenetic-related synthetic lethality in GC (Yang et al., 2019). Moreover, by identifying EMRLs using ChIP/MeRIP-seq data, our study added new evidence to the crosstalk between epigenetic modifications and LncRNA, which is also one of the key members of epigenetic regulation (Kan et al., 2021). Meng et al. defined a DMS score based on specific DNA methylation patterns to recognize GC patients with immune activation status and enhanced efficacy of immunotherapy (Meng et al., 2021). Consistent with Meng’s study, strong links between epigenetic modification patterns and the tumor environment were observed in GC, which in combination suggested epigenetic modification as a promising resource for developing new vulnerabilities for TME-targeting therapy.

With strong complexity and heterogeneity, GC presented disappointing results in most of the clinical trials on novel agents during the last decade (Serra et al., 2019). Thus, increasing molecular classification proposals were developed to characterize the molecular and clinical features of GC so as to optimize individualized diagnosis and treatment. ACRG classification is a globally accepted classification system of GC reported by the Asian Cancer Research Group in 2015, which divides GC into four subtypes: MSI, MSS/EMT, MSS/TM53 (+), and MSS/TP53 (−) (Cristescu et al., 2015). According to the GSE62254 cohort, the MSI subtype had the best prognosis with more than 60% patients of intestinal histotype, while the MSS/EMT subtype had the worst prognosis, over 80% of which were the diffuse histotype. However, the epigenetic modification patterns underlying the MSS/EMT subtype remain largely unexplored. In this study, we surprisingly found the EMD subtype we defined mostly overlapped with the MSS/EMT subtype (91.4%, Supplementary Table S6) in the GSE62254 cohort. Consistently, the EMD subtype had the worst overall survival and mostly consisted of diffuse histotype (80%, Supplementary Table S6). Our findings revealed the distinctive epigenetic modification patterns existing in the MSS/EMT subtype. By identifying the EMD subtype, our study provided not only new strategies for recognizing patients with poor prognosis in GC but also new sights into the epigenetic characterization of the widely used ACRG subtype.

Though with the rapid progress of ICI therapy in recent years, there were still limited GC patients benefiting from the ICI treatment. According to a meta-analysis based on 2003 patients from nine clinical trials, the objective response rate and disease control ratio were 9.9 and 33.3%, respectively, for advanced gastric or gastroesophageal junction (G/GEJ) cancer treated with ICI therapy (Chen et al., 2019). Mechanistically, the immune-suppressed TME may be one of the main causes of the resistance to ICI treatment. Hegde et al. proposed three distinct immunophenotypes (inflamed, immune-excluded, and immune-desert) based on the spatial distribution of CD8⁺ T cells in the TME (Hegde et al., 2016). Immune-excluded tumors were featured with an immune-suppressed TME, represented by T cells clearly embedded in the tumor stromal microenvironment with the presence of low MHC-I expression, TGF-β activation, myeloid inflammation, and angiogenesis (Hegde et al., 2016). Consistent with immune-excluded tumors, the EMD subtype we defined was characterized by the immune-suppressed TME with a low MHC-I expression, TGF-β activation, and angiogenesis activation, which also presented a low response rate to ICI treatment. Moreover, the EMD score we developed showed promising ability in recognizing GC patients with poor prognosis and resistance to ICI therapy. Desbois et al. reported that TGF-β-activated fibroblasts contribute to an immune-suppressed environment by cytokine production in ovarian cancer (Desbois et al., 2020; Desbois and Wang, 2021). In our study, TGF-β activation and high abundance of CAF were both observed in patients with a high EMD score, strongly implying the close links between TGF-β-activated fibroblasts and resistance to ICI treatment in GC.

Recently, combined therapy of ICI with radiotherapy, chemotherapy, or targeted therapy became a prospective strategy for improving the efficacy of cancer treatment. Results from CheckMate 649 showed that nivolumab with chemotherapy improved OS [hazard ratio (HR) 0.71 (98.4% CI 0.59–0.86); p < 0.0001] and PFS [HR 0.68 (98% CI 0.56–0.81); p < 0.0001] in GC patients with a PD-L1 CPS (combined positive score) of five or more when compared with chemotherapy alone (Janjigian et al., 2021). As one of the pivotal roles contributing to the immune-suppressed TME, the EMT pathway provided multiple therapeutic targets for combined therapy of ICI (Dongre et al., 2017; Jiang and Zhan, 2020). Combined ICI therapy with anti-EMT therapy would be a promising strategy, especially for patients with significant EMT activation. In our study, close correlations of the FTO overexpression with a high EMD score and activation of multiple tumorigenic pathways including EMT, TGF-β, and Hedgehog were observed. Moreover, a significant downregulation of the EMT pathway was observed in two FTO^−KD cell lines (Figure 8F, S1C-E), implying the potential links between FTO and the immune-suppressed TME. Su et al. demonstrated that FTO inhibition sensitizes leukemia cells to T cell cytotoxicity and overcomes hypomethylating agent-induced immune evasion (Su et al., 2020). Liu et al. developed a novel FTO inhibitor, Dac51, which can block FTO-mediated immune evasion and synergize with the checkpoint blockade for better tumor control (Liu et al., 2021). In this study, through the in vitro experiments of FTO inhibition conducted in two GC cell lines, we demonstrated that pharmacological inhibition of FTO significantly suppressed the EMT pathway. Our findings provided new clues of FTO’s participation in immune evasion and also suggested FTO as a potential target of combined ICI therapy with anti-EMT therapy in GC.

This study also has some limitations. First, the ability of the EMD score in predicting ICI efficacy needs to be further validated using more GC cohorts in the future. Besides, it remains unclear how FTO regulates the EMT pathway. Further studies are urgently needed to explore the underlying mechanism and verify the efficacy of combined ICI therapy with FTO inhibition in GC.

5 Conclusion

Collectively, we comprehensively characterized the crosstalk between histone and RNA modifications and identified the EMD subtype of GC with poor survival and distinctive TME characteristics. EMD score is a good indicator for prognosis and TME characteristics in GC and also might be a promising tool for recognizing patients suitable for combination ICI therapy.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

CY, YH, and CZ conceived and designed the whole project. CY and JZ collected, processed, and analyzed the data. CD conducted the cell experiments. CY and CD wrote the manuscript. LC and HL reviewed and revised the manuscript. YX, BL, SM, and XJ provided material or technical support for this study. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by Sanming Project of Medicine in Shenzhen (No. SZSM201911010).

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.

Publisher’s Note

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

Acknowledgments

We would like to thank Zhijun Zhou (MD) from the University of Oklahoma Health Sciences Center, and Luyu Xie (PhD) from National Cancer Center Singapore for providing professional suggestions for this study.

Supplementary Material

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

References

Bagaev, A., Kotlov, N., Nomie, K., Svekolkin, V., Gafurov, A., Isaeva, O., et al. (2021). Conserved Pan-Cancer Microenvironment Subtypes Predict Response to Immunotherapy. Cancer Cell 39 (6), 845–e7. doi:10.1016/j.ccell.2021.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Bannister, A. J., and Kouzarides, T. (2011). Regulation of Chromatin by Histone Modifications. Cell Res 21 (3), 381–395. doi:10.1038/cr.2011.22

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B., Khodadoust, M. S., Liu, C. L., Newman, A. M., and Alizadeh, A. A. (2018). Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol. Biol. 1711, 243–259. doi:10.1007/978-1-4939-7493-1_12

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C., Zhang, F., Zhou, N., Gu, Y. M., Zhang, Y. T., He, Y. D., et al. (2019). Efficacy and Safety of Immune Checkpoint Inhibitors in Advanced Gastric or Gastroesophageal junction Cancer: a Systematic Review and Meta-Analysis. Oncoimmunology 8 (5), e1581547. doi:10.1080/2162402X.2019.1581547

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H., Yao, J., Bao, R., Dong, Y., Zhang, T., Du, Y., et al. (2021). Cross-talk of Four Types of RNA Modification Writers Defines Tumor Microenvironment and Pharmacogenomic Landscape in Colorectal Cancer. Mol. Cancer 20 (1), 29. doi:10.1186/s12943-021-01322-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X. Y., Liang, R., Yi, Y. C., Fan, H. N., Chen, M., Zhang, J., et al. (2021). The m6A Reader YTHDF1 Facilitates the Tumorigenesis and Metastasis of Gastric Cancer via USP14 Translation in an m6A-dependent Manner. Front Cel Dev Biol 9, 647702. doi:10.3389/fcell.2021.647702

CrossRef Full Text | Google Scholar

Cho, J. Y., Lim, J. Y., Cheong, J. H., Park, Y. Y., Yoon, S. L., Kim, S. M., et al. (2011). Gene Expression Signature-Based Prognostic Risk Score in Gastric Cancer. Clin. Cancer Res. 17 (7), 1850–1857. doi:10.1158/1078-0432.CCR-10-2180

PubMed Abstract | CrossRef Full Text | Google Scholar

Cristescu, R., Lee, J., Nebozhyn, M., Kim, K. M., Ting, J. C., Wong, S. S., et al. (2015). Molecular Analysis of Gastric Cancer Identifies Subtypes Associated with Distinct Clinical Outcomes. Nat. Med. 21 (5), 449–456. doi:10.1038/nm.3850

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, S., Zhang, H., Zhu, K., Li, X., Ye, Y., Li, R., et al. (2021). M6A2Target: a Comprehensive Database for Targets of m6A Writers, Erasers and Readers. Brief Bioinform 22 (3), bbaa055. doi:10.1093/bib/bbaa055

PubMed Abstract | CrossRef Full Text | Google Scholar

Desbois, M., Udyavar, A. R., Ryner, L., Kozlowski, C., Guan, Y., Durrbaum, M., et al. (2020). Integrated Digital Pathology and Transcriptome Analysis Identifies Molecular Mediators of T-Cell Exclusion in Ovarian Cancer. Nat. Commun. 11 (1), 5583. doi:10.1038/s41467-020-19408-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Desbois, M., and Wang, Y. (2021). Cancer-associated Fibroblasts: Key Players in Shaping the Tumor Immune Microenvironment. Immunol. Rev. 302 (1), 241–258. doi:10.1111/imr.12982

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, L. M., Zhang, X. L., Mao, M. H., Li, Y. P., Zhang, X. Y., Xue, D. W., et al. (2021). LINC00511/miRNA-143-3p Modulates Apoptosis and Malignant Phenotype of Bladder Carcinoma Cells via PCMT1. Front. Cel Dev Biol 9, 650999. doi:10.3389/fcell.2021.650999

CrossRef Full Text | Google Scholar

Dongre, A., Rashidian, M., Reinhardt, F., Bagnato, A., Keckesova, Z., Ploegh, H. L., et al. (2017). Epithelial-to-Mesenchymal Transition Contributes to Immunosuppression in Breast Carcinomas. Cancer Res. 77 (15), 3982–3989. doi:10.1158/0008-5472.CAN-16-3292

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, C., Shi, X., Dai, C., Shen, F., Rocco, G., Chen, J., et al. (2020). RNA m6A Modification in Cancers: Molecular Mechanisms and Potential Clinical Applications. Innov.. 1 (3), 100066. doi:10.1016/j.xinn.2020.100066

CrossRef Full Text | Google Scholar

Gu, Y., Wu, X., Zhang, J., Fang, Y., Pan, Y., Shu, Y., et al. (2021). The Evolving Landscape of N6-Methyladenosine Modification in the Tumor Microenvironment. Mol. Ther. 29 (5), 1703–1715. doi:10.1016/j.ymthe.2021.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hegde, P. S., Karanikas, V., and Evers, S. (2016). The where, the when, and the How of Immune Monitoring for Cancer Immunotherapies in the Era of Checkpoint Inhibition. Clin. Cancer Res. 22 (8), 1865–1874. doi:10.1158/1078-0432.CCR-15-1507

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, H., Weng, H., Zhou, K., Wu, T., Zhao, B. S., Sun, M., et al. (2019). Histone H3 Trimethylation at Lysine 36 Guides m6A RNA Modification Co-transcriptionally. Nature 567, 414–419. doi:10.1038/s41586-019-1016-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Y., Min, S., Lui, Y., Sun, J., Su, X., Liu, Y., et al. (2012). Global Mapping of H3K4me3 and H3K27me3 Reveals Chromatin State-Based Regulation of Human Monocyte-Derived Dendritic Cells in Different Environments. Genes Immun. 13 (4), 311–320. doi:10.1038/gene.2011.87

PubMed Abstract | CrossRef Full Text | Google Scholar

Janjigian, Y. Y., Shitara, K., Moehler, M., Garrido, M., Salman, P., Shen, L., et al. (2021). First-line Nivolumab Plus Chemotherapy versus Chemotherapy Alone for Advanced Gastric, Gastro-Oesophageal junction, and Oesophageal Adenocarcinoma (CheckMate 649): a Randomised, Open-Label, Phase 3 Trial. Lancet 398 (10294), 27–40. doi:10.1016/S0140-6736(21)00797-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeffrey, T. L. (2020). Surrogate Variable Analysis. Biostatistics.

Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T Cell Dysfunction and Exclusion Predict Cancer Immunotherapy Response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Y., and Zhan, H. (2020). Communication between EMT and PD-L1 Signaling: New Insights into Tumor Immune Evasion. Cancer Lett. 468, 72–81. doi:10.1016/j.canlet.2019.10.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Kan, R. L., Chen, J., and Sallam, T. (2022). Crosstalk between Epitranscriptomic and Epigenetic Mechanisms in Gene Regulation. Trends Genet. 38 (2), 182–193. doi:10.1016/j.tig.2021.06.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Kan, R. L., Chen, J., and Sallam, T. (2021). Crosstalk between Epitranscriptomic and Epigenetic Mechanisms in Gene Regulation. Trends Genet. 38, 182–193. doi:10.1016/j.tig.2021.06.01

PubMed Abstract | CrossRef Full Text | Google Scholar

Ke, S., Alemu, E. A., Mertens, C., Gantman, E. C., Fak, J. J., Mele, A., et al. (2015). A Majority of m6A Residues Are in the Last Exons, Allowing the Potential for 3' UTR Regulation. Genes Dev. 29 (19), 2037–2053. doi:10.1101/gad.269415.115

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Paggi, J. M., Park, C., Bennett, C., and Salzberg, S. L. (2019). Graph-based Genome Alignment and Genotyping with HISAT2 and HISAT-Genotype. Nat. Biotechnol. 37 (8), 907–915. doi:10.1038/s41587-019-0201-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R Package for Weighted Correlation Network Analysis. BMC Bioinformatics 9, 559. doi:10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Lawson, K. A., Sousa, C. M., Zhang, X., Kim, E., Akthar, R., Caumanns, J. J., et al. (2020). Functional Genomic Landscape of Cancer-Intrinsic Evasion of Killing by T Cells. Nature 586 (7827), 120–126. doi:10.1038/s41586-020-2746-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Chen, Q., Xue, H., Zhang, L., Wang, K., and Shen, F. (2020). LncRNA MNX1-AS1 Promotes Progression of Intrahepatic Cholangiocarcinoma through the MNX1/Hippo axis. Cell Death Dis 11 (10), 894. doi:10.1038/s41419-020-03029-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Chen, S., Yu, J., Gao, Z., Sun, Z., Yi, Y., et al. (2021). Interplay of M6 A and Histone Modifications Contributes to Temozolomide Resistance in Glioblastoma. Clin. Transl Med. 11 (9), e553. doi:10.1002/ctm2.553

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The Sequence Alignment/Map Format and SAMtools. Bioinformatics 25 (16), 2078–2079. doi:10.1093/bioinformatics/btp352

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Li, Q., Li, D., Shen, Z., Zhang, K., Bi, Z., et al. (2020). Long Non-coding RNA MNX1-AS1 Promotes Progression of Triple Negative Breast Cancer by Enhancing Phosphorylation of Stat3. Front. Oncol. 10, 1108. doi:10.3389/fonc.2020.01108

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Xia, L., Tan, K., Ye, X., Zuo, Z., Li, M., et al. (2020). N6-Methyladenosine Co-transcriptionally Directs the Demethylation of Histone H3K9me2. Nat. Genet. 52 (9), 870–877. doi:10.1038/s41588-020-0677-3

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Liang, G., Xu, H., Dong, W., Dong, Z., Qiu, Z., et al. (2021). Tumors Exploit FTO-Mediated Regulation of Glycolytic Metabolism to Evade Immune Surveillance. Cell Metab 33 (6), 1221–e11. doi:10.1016/j.cmet.2021.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2. Genome Biol. 15 (12), 550. doi:10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayakonda, A., Lin, D. C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: Efficient and Comprehensive Analysis of Somatic Variants in Cancer. Genome Res. 28 (11), 1747–1756. doi:10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, Q., Lu, Y. X., Ruan, D. Y., Yu, K., Chen, Y. X., Xiao, M., et al. (2021). DNA Methylation Regulator-Mediated Modification Patterns and Tumor Microenvironment Characterization in Gastric Cancer. Mol. Ther. Nucleic Acids 24, 695–710. doi:10.1016/j.omtn.2021.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Oh, S. C., Sohn, B. H., Cheong, J. H., Kim, S. B., Lee, J. E., Park, K. C., et al. (2018). Clinical and Genomic Landscape of Gastric Cancer with a Mesenchymal Phenotype. Nat. Commun. 9 (1), 1777. doi:10.1038/s41467-018-04179-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Okabe, A., Huang, K. K., Matsusaka, K., Fukuyo, M., Xing, M., Ong, X., et al. (2020). Cross-species Chromatin Interactions Drive Transcriptional Rewiring in Epstein-Barr Virus-Positive Gastric Adenocarcinoma. Nat. Genet. 52 (9), 919–930. doi:10.1038/s41588-020-0665-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Ooi, C. H., Ivanova, T., Wu, J., Lee, M., Tan, I. B., Tao, J., et al. (2009). Oncogenic Pathway Combinations Predict Clinical Prognosis in Gastric Cancer. Plos Genet. 5 (10), e1000676. doi:10.1371/journal.pgen.1000676

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, X., Li, X., Yang, S., Huang, M., Wei, S., Ma, Y., et al. (2021). LINC00511 Drives Invasive Behavior in Hepatocellular Carcinoma by Regulating Exosome Secretion and Invadopodia Formation. J. Exp. Clin. Cancer Res. 40 (1), 183. doi:10.1186/s13046-021-01990-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Sempere, L. F., Powell, K., Rana, J., Brock, A. A., and Schmittgen, T. D. (2021). Role of Non-coding RNAs in Tumor Progression and Metastasis in Pancreatic Cancer. Cancer Metastasis Rev. 40 (3), 761–776. doi:10.1007/s10555-021-09995-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Serra, O., Galán, M., Ginesta, M. M., Calvo, M., Sala, N., and Salazar, R. (2019). Comparison and Applicability of Molecular Classifications for Gastric Cancer. Cancer Treat. Rev. 77, 29–34. doi:10.1016/j.ctrv.2019.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, H., Wei, J., and He, C. (2019). Where, when, and How: Context-dependent Functions of RNA Methylation Writers, Readers, and Erasers. Mol. Cel 74 (4), 640–650. doi:10.1016/j.molcel.2019.04.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Smyth, E. C., Nilsson, M., Grabsch, H. I., van Grieken, N. C., and Lordick, F. (2020). Gastric Cancer. Lancet 396 (10251), 635–648. doi:10.1016/S0140-6736(20)31288-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Y., Li, L., Ou, Y., Gao, Z., Li, E., Li, X., et al. (2014). Identification of Genomic Alterations in Oesophageal Squamous Cell Cancer. Nature 509 (7498), 91–95. doi:10.1038/nature13176

PubMed Abstract | CrossRef Full Text | Google Scholar

Stark, R., and Brown, G. (2011). DiffBind: Differential Binding Analysis of ChIP-Seq Peak Data. Bioconductor.

Google Scholar

Su, R., Dong, L., Li, Y., Gao, M., Han, L., Wunderlich, M., et al. (2020). Targeting FTO Suppresses Cancer Stem Cell Maintenance and Immune Evasion. Cancer Cell 38 (1), 79–e11. doi:10.1016/j.ccell.2020.04.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, L., Mu, X., Chen, H., Jin, D., Zhang, R., Zhao, Y., et al. (2021). FTO Modifies the m6A Level of MALAT and Promotes Bladder Cancer Progression. Clin. Transl Med. 11 (2), e310. doi:10.1002/ctm2.310

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarasov, A., Vilella, A. J., Cuppen, E., Nijman, I. J., and Prins, P. (2015). Sambamba: Fast Processing of NGS Alignment Formats. Bioinformatics 31 (12), 2032–2034. doi:10.1093/bioinformatics/btv098

PubMed Abstract | CrossRef Full Text | Google Scholar

Team, R.C (2020). A Language and Environment for Statistical Computing. AvaliableAt: https://www.R-project.org

PubMed Abstract | Google Scholar

Thrift, A. P., and El-Serag, H. B. (2020). Burden of Gastric Cancer. Clin. Gastroenterol. Hepatol. 18 (3), 534–542. doi:10.1016/j.cgh.2019.07.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Hui, H., Agrawal, K., Kang, Y., Li, N., Tang, R., et al. (2020). m6 A RNA Methyltransferases METTL3/14 Regulate Immune Responses to Anti-PD-1 therapyA RNA Methyltransferases METTL3/14 Regulate Immune Responses to Anti-PD-1 Therapy. EMBO J. 39 (20), e104514. doi:10.15252/embj.2020104514

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, Z. (2020). exomePeak2: Bias Awared Peak Calling and Quantification for MeRIP-Seq. R package version 1.2.0.

Google Scholar

Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a Class Discovery Tool with Confidence Assessments and Item Tracking. Bioinformatics 26 (12), 1572–1573. doi:10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, F., Zhong, Y., Lang, X. B., Tu, Y. L., and Sun, S. F. (2019). MNX1-AS1 Accelerates the Epithelial-Mesenchymal Transition in Osteosarcoma Cells by Activating MNX1 as a Functional Oncogene. Eur. Rev. Med. Pharmacol. Sci. 23 (19), 8194–8202. doi:10.26355/eurrev_201910_19126

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, J., Chai, H., Shan, H., Pan, C., Xu, X., Dong, W., et al. (2021). Histone Methyltransferase SETD1A Induces Epithelial-Mesenchymal Transition to Promote Invasion and Metastasis through Epigenetic Reprogramming of Snail in Gastric Cancer. Front. Cel Dev Biol 9, 657888. doi:10.3389/fcell.2021.657888

CrossRef Full Text | Google Scholar

Wu, Q. N., Luo, X. J., Liu, J., Lu, Y. X., Wang, Y., Qi, J., et al. (2021). MYC-activated LncRNA MNX1-AS1 Promotes the Progression of Colorectal Cancer by Stabilizing YB1. Cancer Res. 81 (10), 2636–2650. doi:10.1158/0008-5472.CAN-20-3747

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Y., Li, L., Wang, Q., Zhang, L., He, C., Wang, X., et al. (2020). LINC00511 Promotes Lung Squamous Cell Carcinoma Proliferation and Migration via Inhibiting miR-150-5p and Activating TADA1. Transl Lung Cancer Res. 9 (4), 1138–1148. doi:10.21037/tlcr-19-701

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, H., Cui, W., and Wang, L. (2019). Epigenetic Synthetic Lethality Approaches in Cancer Therapy. Clin. Epigenetics 11 (1), 136. doi:10.1186/s13148-019-0734-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, H., Gao, L., Zhang, M., Ning, N., Wang, Y., Wu, D., et al. (2021). Identification and Analysis of an Epigenetically Regulated Five-lncRNA Signature Associated with Outcome and Chemotherapy Response in Ovarian Cancer. Front. Cel Dev Biol 9, 644940. doi:10.3389/fcell.2021.644940

CrossRef Full Text | Google Scholar

Yang, Z., Jiang, X., Zhang, Z., Zhao, Z., Xing, W., Liu, Y., et al. (2021). HDAC3-dependent Transcriptional Repression of FOXA2 Regulates FTO/m6A/MYC Signaling to Contribute to the Development of Gastric Cancer. Cancer Gene Ther. 28 (1-2), 141–155. doi:10.1038/s41417-020-0193-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, J., Leavenworth, J. W., Li, Y., Luo, Q., Xie, H., Liu, X., et al. (2015). Ezh2 Regulates Differentiation and Function of Natural Killer Cells through Histone Methyltransferase Activity. Proc. Natl. Acad. Sci. U S A. 112 (52), 15988–15993. doi:10.1073/pnas.1521740112

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS 16 (5), 284–287. doi:10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, L., Xu, Z. Y., Ruan, S. M., Mo, S., Qin, J. J., and Cheng, X. D. (2020). Long Non-coding RNAs towards Precision Medicine in Gastric Cancer: Early Diagnosis, Treatment, and Drug Resistance. Mol. Cancer 19 (1), 96. doi:10.1186/s12943-020-01219-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Yue, B., Song, C., Yang, L., Cui, R., Cheng, X., Zhang, Z., et al. (2019). METTL3-mediated N6-Methyladenosine Modification Is Critical for Epithelial-Mesenchymal Transition and Metastasis of Gastric Cancer. Mol. Cancer 18 (1), 142. doi:10.1186/s12943-019-1065-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, H., Liu, X., Yu, L., Lin, S., Zhang, C., Xu, H., et al. (2021). Comprehensive Landscape of Epigenetic-Dysregulated lncRNAs Reveals a Profound Role of Enhancers in Carcinogenesis in BC Subtypes. Mol. Ther. Nucleic Acids 23, 667–681. doi:10.1016/j.omtn.2020.12.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, L., Kong, X., Zhong, W., Wang, Y., and Li, P. (2020). FTO Accelerates Ovarian Cancer Cell Growth by Promoting Proliferation, Inhibiting Apoptosis, and Activating Autophagy. Pathol. Res. Pract. 216 (9), 153042. doi:10.1016/j.prp.2020.153042

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Chen, Y., Jin, M., and Wang, J. (2021). The Crosstalk between m6A RNA Methylation and Other Epigenetic Regulators: a Novel Perspective in Epigenetic Remodeling. Theranostics 11 (9), 4549–4566. doi:10.7150/thno.54967

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Z., and Shilatifard, A. (2019). Epigenetic Modifications of Histones in Cancer. Genome Biol. 20 (1), 245. doi:10.1186/s13059-019-1870-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, G., Yan, K., Liu, J., Gao, L., Jiang, X., and Fan, Y. (2021). FTO Promotes Tumour Proliferation in Bladder Cancer via the FTO/miR-576/CDK6 axis in an m6A-dependent Manner. Cell Death Discov 7 (1), 329. doi:10.1038/s41420-021-00724-5

PubMed Abstract | CrossRef Full Text | Google Scholar