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.

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.

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).

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):

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).

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.

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: E M D   s c o r e = log 2 {   s s G S E A   s c o r e   ( s e t 1 ) /   s s G S E A   s c o r e   ( s e t 2 ) } .

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).

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.

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.

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.

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.

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).

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.

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).

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.

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).

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.

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).

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.

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).

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.

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.

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.