Miaomiao Huo1
Min Zhang1
Jingyao Zhang1
Yong Wang2
Ting Hu1Tianyu Ma1Yinuo Wang1
Baowen Yuan1Hao Qin1
Xu Teng3
Hefen Yu3
Wei Huang3*
Yan Wang1*- 1Key Laboratory of Cancer and Microbiome, State Key Laboratory of Molecular Oncology, National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
- 2Department of Ultrasound, National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
- 3Beijing Key Laboratory of Cancer Invasion and Metastasis Research, Department of Biochemistry and Molecular Biology, School of Basic Medical Sciences, Capital Medical University, Beijing, China
Background: Cancers arise from genetic and epigenetic abnormalities that affect oncogenes and tumor suppressor genes, compounded by gene mutations. The N6-methyladenosine (m6A) RNA modification, regulated by methylation regulators, has been implicated in tumor proliferation, differentiation, tumorigenesis, invasion, and metastasis. However, the role of m6A modification patterns in the tumor microenvironment of gastric cancer (GC) remains poorly understood.
Materials and methods: In this study, we analyzed m6A modification patterns in 267 GC samples utilizing 31 m6A regulators. Using consensus clustering, we identified two unique subgroups of GC. Patients with GC were segregated into high- and low-infiltration cohorts to evaluate the infiltration proportions of the five prognostically significant immune cell types. Leveraging the differential genes in GC, we identified a “green” module via Weighted Gene Co-expression Network Analysis. A risk prediction model was established using the LASSO regression method.
Results: The “green” module was connected to both the m6A RNA methylation cluster and immune infiltration patterns. Based on “Module Membership” and “Gene Significance”, 37 hub genes were identified, and a risk prediction model incorporating nine hub genes was established. Furthermore, methylated RNA immunoprecipitation and RNA Immunoprecipitation assays revealed that YTHDF1 elevated the expression of DNMT3B, which synergistically promoted the initiation and development of GC. We elucidated the molecular mechanism underlying the regulation of DNMT3B by YTHDF1 and explored the crosstalk between m6A and 5mC modification.
Conclusion: m6A RNA methylation regulators are instrumental in malignant progression and the dynamics of tumor microenvironment infiltration of GC. Assessing m6A modification patterns and tumor microenvironment infiltration characteristics in patients with GC holds promise as a valuable prognostic biomarker.
1 Introduction
RNA methylation plays a significantly role in normal cellular homeostasis and pathological conditions (Han et al., 2020), in which N6-methyladenosine (m6A) is the most prevalent in eukaryotic cells and has gained increasing attention because of its presence in mRNA, lncRNAs, and miRNA (Dominissini et al., 2012; Meyer et al., 2012; Ontiveros et al., 2019; Huang et al., 2020). Recently, m6A modification, a dynamic and reversible epigenetic process, has attracted considerable attention. This modification is orchestrated by regulators commonly referred to as “writers”, “erasers”, and “readers”. The methylation process is specifically catalyzed by methyltransferases, or “writers”, which encompass enzymes such as METTL3, METTL14, WTAP, and METTL16. Conversely, the demethylation process is executed by demethylases such as FTO and ALKBH5, termed “erasers”. Additionally, there is a set of RNA-binding proteins encompassing, but not restricted to, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1 (Zhang et al., 2020). Numerous studies have reported that m6A modifications are common in cancer. Such modifications profoundly influence tumorigenesis and tumor progression by disrupting cellular pathways; promoting cell proliferation, self-renewal, and tumor metastasis; and leading to aberrations in immunomodulation (Pinello et al., 2018; Tong et al., 2018; Chen et al., 2020; Dong and Cui, 2020; Haruehanroengra et al., 2020). Furthermore, programmed cell death pathways, which are closely associated with the cancer initiation, progression, and resistance, have highly complex links to m6A modification (Liu et al., 2022).
Gastric cancer (GC) is one of the most prevalent digestive tract cancers globally, with an incidence rate of 5.6% and a mortality rate of 7.7%, ranks the top five in both categories (Sung et al., 2021). With progress in biological information technology and medical means, genome analysis has become the main method for identifying new biological targets in GC (Cancer Genome Atlas Research, 2014; Oh et al., 2018). The overall abundance of m6A mRNA in human GC tissues is significantly higher compared to normal tissue (Wang Q. et al., 2020). Studies have found that tumor progression is not solely attributed to genetic and epigenetic modifications in tumor cells. The tumor microenvironment (TME), on which cancer cells rely for growth and survival, plays a pivotal role. Cancer cells modify various biological behaviors through direct and indirect interactions. As our comprehension of the complexities and variations within TME deepening, increasing evidence highlights the pivotal role in tumor progression and immune evasion, as well as its influence on responses to immunotherapy (Quail and Joyce, 2013; Ali et al., 2016). Consequently, there is a growing emphasis on studying biomarkers that can predict responses to immune checkpoint blockade therapies, aiming to enhance precision immunotherapy strategies. Virtually, m6A modification alterations intrinsically affect immune cells and extrinsically affect immune cell responses in the TME (Li et al., 2022). m6A is closely related to macrophage phenotype and dysfunction (Zhu X. et al., 2023). In GC, m6A modification plays a non-negligible role in characterizing TME infiltration, both in terms of diversity and complexity (Zhang et al., 2020).
In this study patient samples are divided into two subgroups by conducting a consistent cluster analysis of the expression profiles of 31 m6A RNA methylation factors from patients with GC found in the TCGA database. Based on the proportion of infiltration from the five immune cell types related to prognosis, patients with GC were categorized into high-infiltration and low-infiltration groups, where there was a significant difference in prognosis. Using WGCNA analysis of the differential genes in GC, we identified a module associated with both the m6A methylation cluster and the immune infiltration classification. Hub genes were isolated from this module using module membership (MM) and gene significance (GS). By leveraging these hub genes, we employed LASSO regression to develop a risk prediction model. Based on the risk score of this model, we further categorized the samples into high-risk and low-risk groups. Furthermore, we validated the isolated hub genes and found that the m6A reader, YTHDF1, elevated the expression of DNMT3B, synergistically promoting the initiation and development of GC.
2 Materials and methods
2.1 GC datasets source and preprocessing
We extracted gene expression datasets and associated clinical annotations from two renowned repositories: UCSC Xena (accessible at https://xenabrowser.net/datapages/) and cBioportal (available at http://www.cbioportal.org/). We included a cohort of 267 gastric cancer patient samples along with 32 samples derived from healthy individuals. These samples were repleted with pertinent survival data and detailed tumor staging, all sourced from The Cancer Genome Atlas-Stomach Adenocarcinoma (TCGA-STAD) database download from UCSC Xena. A comprehensive list of these details is provided in Table 1. We also integrated the GC dataset GSE62254 from the Gene Expression Omnibus, constituting an additional 300 GC patient samples; patients without survival information were excluded from further evaluation. The comprehensive sample data are presented in Table 2. The workflow of data analysis is shown in Figure 1.
Table 1. Clinical characteristics of the patients from TCGA-STAD.
Table 2. Clinical characteristics of the patients from GSE62254.
Figure 1. Flow chart of analysis.
2.2 Consistent cluster analysis
A list of the 31 m6A regulators is provided in the Supplementary Material (Table S1). Utilizing the R package “ConsensusClusterPlus” (Wilkerson and Hayes, 2010), we performed a consistent clustering analysis with a robust setting of 100 iterations and an 80% resampling rate, leveraging the Pearson correlation as the chosen distance metric. This rigorous analysis stratified the samples into two distinct clusters: Cluster 1 and Cluster 2.
A Principal Component Analysis (PCA) was performed to visually discern and compare the underlying variations between the two clusters. This dimensionality reduction technique allowed us to capture the essence of variance in the data and offered a clearer perspective on the differences in expression patterns. The Wilcoxon test was used to ascertain the specific m6A RNA methylation regulators that exhibited significant differential expression between Cluster 1 and Cluster 2. This non-parametric statistical test, tailored for datasets that did not necessarily follow a normal distribution, enabled us to rigorously compare the expression levels of each regulator between the two defined clusters.
2.3 Survival analysis of m6A cluster and subgroup functional pathway analysis
Survival analysis encompasses statistical methods dedicated to exploring the expected time until one or more events occur. We used the Kaplan-Meier method to generate survival curves, and the log-rank test was used to identify prognostic factors correlated with survival, with significance determined at a threshold of p < 0.05. Furthermore, for enrichment analysis, we leveraged the “GSVA” package in R, which employs an unsupervised, non-parametric approach (Hanzelmann et al., 2013; Shen et al., 2022). Significance level of p < 0.05 was considered statistically significance.
2.4 Estimation of the proportion of TME cell infiltration in GC
CIBERSORT (accessible at https://cibersort.stanford.edu/) was synergistically paired with the LM22 signature matrix, facilitating the estimation of the proportion of human hematopoietic cell phenotypes within the 22 samples categorized from both the high-risk and low-risk patient cohorts. Notably, the cumulative proportion of all the inferred immune cell types within each sample was 1. Subsequently, to identify the immune cells with significant prognostic implications, the proportions of these various immune cells were subjected to univariate Cox regression analyses.
2.5 Identification of differentially expressed genes (DEGs)
Differentially expressed genes between 267 tumor samples and 32 normal samples were identified using the Limma R package (Smyth, 2004). Genes were considered differentially expressed based on thresholds of |log2FC| > 0.585 and FDR < 0.01.
2.6 Weighted gene co-expression network analysis (WGCNA)
WGCNA is an advanced bioinformatics approach aimed at deciphering the complex patterns of gene expression data (Yin et al., 2020). By calculating the correlation between the module eigengenes and the clinical traits of interest, biologically relevant modules were identified. Within the modules significantly associated with clinical traits, hub genes were identified based on their connectivity, highlighting those with potential key roles in the module’s biological function.
2.7 Construction and validation of the prognosis model
For the identified hub genes, the LASSO method was employed to select prognostically relevant genes and construct a prognostic model. The systematically derived risk-scoring formula was as follows:
Based on the median risk score derived from the model or the optimal cutoff value calculated using the surv_cutpoint function, patient samples were stratified into high-risk and low-risk groups. Kaplan-Meier survival analysis was used to assess the predictive capability of the model.
2.8 Cell culture and transfection
Cell lines AGS, BGC-823 and HGC-27 were obtained from Shanghai Institute of Biochemical Cell Science, Chinese Academy of Sciences. BGC-823 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) and AGS cells were cultured in DMEM/F12. All the culture mediums were supplemented with 10% fetal bovine serum (FBS), 100 units/mL penicillin, and 100 mg/mL streptomycin (Gibco, 15140-122, United States). The cells were cultured in a constant temperature incubator equilibrated with 5% CO2 at 37°C. The sources and culture conditions of other cell lines are detailed in Supplementary Table S2. All experiments were performed with mycoplasma-free cells.
All plasmids were verified by DNA sequencing and transfected using TurboFect™ reagent (ThermoFisher, R0531, USA). SiRNAs for gene knockdown were transfected using Lipofectamine RNAiMAX reagent (Invitrogen, 13778-150, United States). The siRNA sequences used are listed in the Supplementary Table S3.
2.9 Real-time quantitative polymerase chain reaction (real-time qPCR) analysis
Total RNA was extracted using TRIzol reagent. Subsequently, reverse transcription was performed using the PrimeScript™ RT Master Mix (Takara, RR036A, Japan), according to the manufacturer’s instructions. qPCR was conducted using FastStart Universal SYBR Green Master (Roche, 4913914001, United States), where β-actin was analyzed as the loading control. The relative expression of target genes was calculated using the 2−ΔΔCT method. The primers used are listed in Supplementary Table S4.
2.10 Western blot
Total proteins were separated via SDS-PAGE. The proteins were transferred onto the PVDF membrane and blocked using 5% skim milk at room temperature for 1 h, followed by immunoblotting with the indicated antibodies overnight, including anti-YTHDF1 (Proteintech, Cat No.17479-1-AP, China) and anti-DNMT3B (Cell Signaling Technology, 57868S, United States). After incubation with secondary antibodies for 1 h at room temperature, the membranes were washed and transferred onto an X-ray radiographic cassette and treated with ECL Super Signal™ West Pico PLUS (ThermoFisher, 34580, USA). Subsequently, the membranes were blotted onto X-ray films for visualization.
2.11 Methylated RNA immunoprecipitation (MeRIP) and RIP
The MeRIP-qPCR assay was conducted using a MeRIP assay kit, according to the manufacturer’s instructions (Bersinbio, China). The RIP assay was performed using approximately 2 × 107 cells per sample, and the specific experimental steps were based on previously methods reported (Gagliardi and Matarazzo, 2016). The m6A sites of DNMT3B were predicted using SRAMP (http://www.cuilab.cn/sramp) (Zhou et al., 2016), and primers containing m6A sites were subsequently designed. The primers used are listed in Supplementary Table S5.
2.12 EdU assay
GC cells with overexpression or depletion of the indicated genes and control were seeded into 96-well plates at a density of 6 × 104 cells per well. DNA proliferation was detected after the cells were cultured overnight using an EdU assay kit (RiboBio, Bes5203-1, China). Images were acquired using an inverted fluorescence microscope and statistically analyzed using ImageJ software.
2.13 Cell migration and invasion assay
For the cell migration assay, Transwell chamber filters were placed in 24-well plates. Cells transfected with indicated siRNAs were suspended in serum-free medium, and 8 × 104 cells were seeded into the upper chamber of the wall, whereas the lower chamber was cultured in medium containing 15% FBS. Following incubation for 24 h, the cells were fixed with 4% paraformaldehyde for 20 min, then stained with 0.1% crystal violet after which the cells in the upper chamber were removed. Images were acquired under an inverted microscope and statistically analyzed using ImageJ software. For the cell invasion assay, Transwell chamber filters wrapped in 10% Matrigel (Corning, 354234, United States) were used.
2.14 Statistical analysis
Data analyses were performed using GraphPad Prism (version 9.1.1) and results are displayed as the mean ± SD. Student’s t-test was used to compare differences between the two groups. Potential m6A modification sites were predicted using SRAMP (http://www.cuilab.cn/sramp). Survival curves were plotted using the Kaplan-Meier “survival” package in R (version 3.4.3), where the log-rank test was used to assess statistical significance. Statistical significance was set at p < 0.05.
3 Results
3.1 Landscape of TME in GC and infiltration characteristics in distinct m6A modification patterns
m6A methylation modification patterns mediated by 31 regulators were analyzed in patient-derived GC samples. The m6A methylation mediated by regulators classified as “writers”, “erasers” and “readers” is a dynamic reversible process (Supplementary Figure S1A). We employed the R package “ConsensusClusterPlus”, a tool specifically designed for the robust class discovery and visualization of gene expression datasets, and all GC samples into two distinct groups (Figure 2A). PCA was used to evaluate these groups, revealing significant differences between them (Figure 2B). We utilized the PAM clustering method, with the sample correlation coefficient calculated using Pearson correlation. As depicted in the cumulative consistency distribution map (Supplementary Figures S1B, S1C), there was a noticeable increase in the broken line beyond K = 2, prompting us to categorize all GC samples into two distinct groups. We conducted a Wilcoxon test on the expression levels of 31 m6A regulators across both groups. The findings highlighted that 17 methylation factors, including YTHDC2, IGF2BP2, YTHDC1, and HNRNPC, exhibited significant variance in expression between two sample groups (Figure 2C).
Figure 2. Landscape of TME in gastric cancer and infiltration characteristics in distinct m6A modification patterns. (A) All gastric cancer samples were divided into two groups (K = 2); (B) Under PCA algorithm, the two subgroups cluster 1 and cluster 2 showed significant difference; (C) The expression of 31 m6A regulators between 2 m6A regulators genes cluster: cluster 1 and cluster 2. Cluster 1, Orange; Cluster 2, blue. The upper and lower ends of the boxes represented inter quartile range of values. The lines in the boxes represented median value, and black dots showed outliers; (D) Kaplan-Meier curves for disease progress probability gastric cancer patients from two clusters; (E) Subgroup analysis estimating clinical prognostic value in different types of immune cell infiltration. The length of the horizontal line represents the 95% confidence interval for each group; (F) Kaplan-Meier curves for patients with high and low cohort. Log-rank test shows an overall p = 0.023; (G) Heat map of distribution of five kinds of immune cells in different immune invasion groups; (H) Gene expression level of m6A methylation regulator in different immune infiltration groups; (I) Correlation between the expression of m6A regulator and the proportion of immune cells. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001; two-tailed unpaired t-test.
Survival outcomes were assessed using the Kaplan-Meier log-rank test. Factors defined as p < 0.05 were determined to be prognostic determinants pertinent to survival rates. Incorporating the survival data of the patients, a pronounced disparity in prognosis was observed between the two sample clusters: cluster 1 (N = 46) and cluster 2 (N = 221) (Figure 2D). Interestingly, no significant variance in mutational count was observed between these clusters (Supplementary Figure S1D). Pathway enrichment analysis was performed on distinct subgroups using the GSVA package in R, with significance adjusted to p < 0.05. The findings revealed associations of distinct subgroups with several pathways and functional modalities, including the p53 pathway, interferon γ response, and reactive oxygen species pathways (Table 3; Supplementary Figure S1E).
Table 3. Pathway enrichment analysis of distinct subgroups.
Using the CIBERSORT algorithm in tandem with the LM22 signature matrix, we ascertained the proportional abundance of 22 immune cell subtypes within GC specimens. Univariate Cox regression analysis was performed on 19 immune cell categories. Three immune cell types, namely CD4 naïve (detected in only three samples), T cell gamma delta (present in 18 samples), and eosinophils (identified in 21 samples), were excluded from the analysis because of their infrequent occurrence, rendering them unsuitable for a statistically robust Cox regression analysis. Our investigation revealed that naïve B cells, plasma cells, activated CD4+ memory T cells, and five additional immune cell types exhibited a pronounced correlation with prognosis (Figure 2E). By incorporating patient survival information, we observed a pronounced difference in prognosis between the two sample types characterized by distinct immune infiltration (Figure 2F). Subsequently, based on the infiltration metrics of these five salient immune cells, patient specimens were stratified into two discrete clusters, termed high-infiltration and low-infiltration, using the K-means unsupervised clustering technique. Within these classifications, the sample distribution was 196 with high infiltration and 71 with low infiltration (Figure 2G).
According to the immune infiltration grouping, Wilcoxon test was performed on the expression of 31 m6A regulators in the two groups of samples. Our analysis revealed that the expression of m6A methylation factors, notably YTHDC2, METTL16, and YTHDF1, was significantly different between the two groups (Figure 2H). Subsequently, the Pearson correlation coefficient between the 31 m6A RNA methylation factors and the infiltration proportion of the 22 immune cell types was calculated using the R package psych. Remarkably, most m6A methylation factors were correlated with the infiltration proportions of certain immune cells, particularly CD4 memory-activated T cells and follicular helper T cells (Figure 2I).
3.2 Correlation of hub genes with m6A regulators and immune infiltration
The gene expression profile data of 267 patient samples and 32 control samples from the TCGA database were used for differential gene screening and differential gene expression analysis. According to the multiple differences (|log2FC| > 0.585) and significance threshold (FDR < 0.01), 4,391 DEGs were screened using the R package Lima, including 3,232 upregulated and 1,159 downregulated genes (Figures 3A, B).
Figure 3. Correlation of hub genes with m6A regulators and immune infiltration. (A) Display of differential gene volcano map; (B) The heat map shows the top 100 |log2FC| of differential genes; (C) Gene dendrogram and module colors; (D) Module-trait relationship; (E) Correlations between hub genes. Negative correlation was marked with blue and positive correlation with red; (F) The correlation between immune infiltration types and hub genes DNMT3B, ZNF64, CSTF1; (G) The correlation between m6A regulators and hub genes DNMT3B, ZNF64, CSTF1; (H) Analysis of the differences in immune infiltration and SCNAs among hub genes. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001; two-tailed unpaired t-test.
Leveraging 4,391 DEGs, we constructed a weighted gene co-expression network using the R package WGCNA. Cluster analysis indicated the presence of an outlier sample; therefore, subsequent analyses focused on the remaining 266 samples. Our analysis confirmed that the constructed co-expression network adhered to a scale-free topology. Specifically, the logarithm log (k) of nodes with connectivity k exhibited a negative correlation with the logarithm log [P (k)] of the node occurrence probability, achieving a correlation coefficient exceeding 0.8. To ensure the scale-free nature of the network, the optimal soft-thresholding power was determined to be β = 4 (Supplementary Figures S2A, S2B). Subsequently, the expression matrix was transformed into an adjacency matrix, which was then converted into a TOM.
Genes were clustered using the average link hierarchical clustering method. By adhering to the hybrid dynamic tree-cutting standard, a minimum module size of 30 genes was established. After determining the gene modules using the dynamic tree-cutting approach, the eigenvectors for each module were computed. Module clustering was then performed, amalgamating closely related modules into unified modules with a set threshold of height = 0.25. This process resulted in the identification of 13 distinct modules (Figure 3C). The statistics for the number of genes in each module are shown in Supplementary Figure S2C. We assessed the Pearson correlation coefficient between the ME of each module and the phenotypic traits of the samples. The significance of the module increased with increasing correlation coefficients. Supplementary Figure S2D shows the eigengenes of each module with the accompanying list denoting the phenotypic characteristics of each sample. Subsequently, the significance value for each gene module was calculated (Figure 3D). A heightened GS value underscores a module’s increased relevance to cluster1 samples. The “Green” module was the most significant module, and the related genes were displayed in pathway enrichment analysis (Supplementary Figure S2E). From the amalgamated insights derived from the module-phenotypic correlation analyses, the green module emerged as the most pertinent module in relation to m6A clusters and immune infiltration. Using the criteria of MM >0.7 and GS > 0.1, 37 hub genes were identified in the green module, which exhibited strong interrelations (Figure 3E).
Subsequently, we evaluated the correlation between the identified 37 hub genes and the proportion of immune-infiltrating cells and m6A methylation regulators. Most hub genes demonstrated significant correlations with the proportions of Macrophages M0, as exemplified by DNMT3B, ZNF64, and CSTF1 (Figure 3F; Supplementary Figures S2F, S3A). Furthermore, m6A regulators were also correlated with the hub genes DNMT3B, ZNF64, and CSTF1 (Figure 3G; Supplementary Figures S2G, S3B). Most of the hub genes were significantly correlated with YTHDF1 and IGF2BP1. Notably, DNMT3B was positively correlated with IGF2BP1, YTHDF1, and VIRMA, and negatively correlated with YTHDC2. Furthermore, we investigated the effects of the hub genes on immune cell infiltration using the TIMER database. Different types of somatic copy number alterations regulate immune cell infiltration into the GC microenvironment. Hub genes and YTHDF1, which are significantly associated with these genes, markedly affected various types of immune-infiltrating cells Figure 3H; Supplementary Figure S2H).
3.3 Construction of a prognostic model based on hub genes
LASSO regression analysis was used to analyze the trajectory of the independent variables, where the x-axis represents the logarithm of the variable λ and the y-axis denotes the coefficients of the independent variables (Figure 4A). Figure 4B shows the confidence intervals corresponding to each lambda value within the LASSO regression framework. In our subsequent analyses, the survival information of the nine hub genes identified in GC were analyzed, among which the expression of DNMT3B is significantly associated with poor prognosis (Figure 4C).
Figure 4. Construction and validation of a prognostic model based on hub genes. (A) The trajectory of independent variables in Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis of 9 hub genes in gastric cancer; (B) Line plot LASSO regression analysis of 9 hub genes in gastric cancer; (C) Kaplan-Meier analysis of DNMT 3B high- and low-expression groups for gastric cancer from TCGA; (D) ROC curve of 1-year, 3-year and 5-year survival probability in gastric cancer from TCGA; (E) Kaplan-Meier analysis of high-risk and low-risk groups for gastric cancer from TCGA; (F) Classify gastric cancer patients into different risk groups according to the median risk score; (G) Distribution of risk score in gastric cancer; (H) Heatmap of the 9 prognosis-related hub genes expression profiles combined with clinical traits in the high-risk and low-risk groups in the prognostic model; (I) Kaplan-Meier analysis of high-risk and low-risk groups for gastric cancer from GSE62254 dataset; (J) Classify gastric cancer patients from GSE62254 dataset into different risk groups according to the median risk score; (K) Distribution of risk score in gastric cancer; (L) Heatmap of the 9 prognosis-related hub genes expression profiles GSE62254 dataset; (M) ROC curve of RiskScore, Age, Gender and Stage in gastric cancer from TCGA.
The risk prediction model was assessed using a ROC curve analysis, and the area under the ROC curve values for GC patient samples at 1, 3, and 5 years were 0.645, 0.678, and 0.748, respectively (Figure 4D). We computed the risk score for each sample, and the samples were divided into high- and low-risk group based on the median score. The two groups showed significant difference in prognosis. Kaplan-Meier curve analysis of disease progression probability was conducted to evaluate the effects of the low- and high-risk score groups on prognosis (Figure 4E). Results indicated that the prognosis of high-risk group was worse than that of low-risk one (Figures 4F, G). The heatmap generated from the expression profiles of the nine hub genes in the prognostic model illustrate the expression trends of the differential hub genes as the risk score of the sample increased (Figure 4H).
To validate the efficacy of the risk-scoring model, we conducted a validation using the GSE6254 dataset. Risk scores were computed, based on the optimal cutoff value, samples were categorized into high- and low-risk groups. Similarly, a significant prognostic difference was observed between groups. Kaplan-Meier curve analysis of disease progression probability depicted the low- and high-risk score groups for prognosis in GSE62254 (Figure 4I). The prognosis of high-risk patients was worse than that of low-risk ones (Figures 4J, K). The expression profiles from the GSE62254 dataset of the nine hub genes in the prognostic model are presented in a heatmap (Figure 4L). Compared to Figure 4H, the differential hub genes displayed similar trends as the risk score of the sample increased. Using ROC curves to display different clinical factors and the “RiskScore” to distinguish the survival probability of samples, the results showed that “RiskScore” had the best classification effect (Figure 4M).
3.4 Validation of hub genes
To verify the outcome prediction value of these nine hub genes in GC, we obtained their expression patterns in GC samples from the HPA database (https://www.proteinatlas.org/). Consistent with the above results, high expression of DNMT3B, TM9SF4, and TTI1 was observed in GC, whereas the expression of ZFP64, TOMM34, CSTF1, ACTR5, STK35, and SS18L1 was lower in GC (Figure 5A). Among these, DNMT3B, ACTR5 and TM9SF4 demonstrated significant prognostic differences when expressed at high or low levels (Supplementary Figure S4). To further investigate the expression of these nine hub genes in GC, we measured their expression in the normal human gastric epithelial cell GES-1 and in the cancer cell SNU-5, AGS, BGC-823, SGC-7901, MGC-803, and HGC-27 by qPCR. DNMT3B expression was significantly higher in tumor cells than normal epithelial cell GES-1 (Figures 5B, C).
Figure 5. Validation of hub genes. (A) Immunohistochemistry (IHC) of 9 hub genes in gastric cancer and normal samples from the HPA database (https://www.proteinatlas.org/); (B) Expression profile of 9 hub genes (ZFP64, TOMM34, DNMT3B, CSTF1, TM9SF4, TTI1, ACTR5, STK35, SS18L); (C) The mRNA expression of DNMT3B in the normal human gastric epithelial cell line GES-1 and gastric cancer cell lines SNU-5, AGS, BGC-823, SGC-7901, MGC-803, HGC-27; (D) m6A dot blot of total RNA in the normal human gastric epithelial cell line GES-1 and gastric cancer cell lines SNU-5, AGS, BGC-823, SGC-7901, MGC-803, HGC-27; (E) m6A dot blot of total RNA in the in adjacent and tumor tissues of gastric cancer samples; (F) Analysis of the mRNA expression of YTHDF1 and IGF2BP1 in the normal human gastric epithelial cell line GES-1 and gastric cancer cell lines SNU-5, AGS, BGC-823, SGC-7901, MGC-803, HGC-27 via RT-qPCR; (G) Analysis of expression of YTHDF1 and IGF2BP1 in the normal human gastric epithelial cell line GES-1 and gastric cancer cell lines SNU-5, AGS, BGC-823, SGC-7901, MGC-803, HGC-27, β-actin served as a loading control; (H) Analysis of the proteins from AGS and BGC-823 cells transfected with vector, YTHDF1, or NC, two different siRNA against YTHDF1 were analyzed by Western blot; (I) Analysis of the mRNA expression of YTHDF1 in AGS and BGC-823 cells transfected with vector, YTHDF1, or NC, two different siRNA against YTHDF1 were analyzed by RT-qPCR; (J) Analysis of the mRNA expression of DNMT3B in AGS and BGC-823 cells transfected with vector, YTHDF1, or NC, two different siRNA against YTHDF1 were analyzed by RT-qPCR; (K) Me-RIP assays; (L) RIP assays; (M) 5mC DNA were analyzed in AGS and BGC-823 cells transfected with vector, YTHDF1, or NC, two different siRNA against YTHDF1 by 5mC dot blot. Error bars represent the mean ± SD of three independent experiments. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001; two-tailed unpaired t-test.
We analyzed the modification levels of m6A in normal epithelial and GC cells using m6A dot blot experiments. Interestingly, the modification level of m6A was higher than that in normal epithelial cells (Figure 5D). Similarly, in RNA extracted from clinical tissue samples, m6A modification in tumor tissue was observed to be higher than in adjacent tissues (Figure 5E). Based on these findings, we hypothesize that in the pathogenesis of gastric cancer, there is an upregulation of m6A modifications, which in turn regulates key genes critical to the oncogenic processes, thereby facilitating the progression of cancer.
Previous results reported a positive correlation between the expression of DNMT3B and that of m6A “reader” YTHDF1 and IGF2BP1. We speculate that the modification of m6A may lead to an increase in the expression of DNMT3B, during which YTHDF1 or IGF2BP1 recognizes the m6A modification site of DNMT3B mRNA, thus stabilizing DNMT3B mRNA and promoting the occurrence and development of GC. Therefore, the relationship between m6A modification and DNMT3B expression should be further investigated in future research. We detected the expression of YTHDF1 or IGF2BP1 in normal epithelial and tumor cell lines. YTHDF1 was highly expressed in cancer cells, while IGF2BP1 did not exhibit statistically significant differences (Figures 5F, G). To further investigate the molecular mechanism by which YTHDF1 regulates DNMT3B, YTHDF1 was overexpressed and knocked down in AGS and BGC-823 cells. As expected, DNMT3B exhibited changes consistent with those of YTHDF1 (Figures 5H–J). Similarly, we also examined the expression changes of DNMT3B after overexpression or knockdown of IGF2BP1. The results showed that IGF2BP1 did not significantly affect the mRNA expression levels of DNMT3B (Supplementary Figure S5A). Furthermore, by overexpressing or knocking down the m6A methyltransferase METTL3 to alter the overall m6A modification level in cells, we examined nine hub genes, including DNMT3B. The results showed that an increase in the overall m6A modification level upregulated ZFP64, TOMM34, DNMT3B, TM9SF4, whereas a decrease downregulated DNMT3B (Supplementary Figure S5B, S5C).
To validate the association between YTHDF1 and DNMT3B mRNA, RIP and MeRIP assays were conducted. Two pairs of primers were designed respectively, according to the two highest confidence m6A methylation sites predicted by SRAMP. Subsequently, primer #1 and primer #2 were identified by qPCR (Figure 5K). Then, a direct interaction between YTHDF1 and DNMT3B mRNA was validated by RIP-qPCR with anti-YTHDF1 antibody (Figure 5L). Then we detected the DNA methylation levels in cells overexpressed or knocked down YTHDF1 using 5mC dot blot. We also performed methylene blue staining as a nucleic acid loading control. Accompanied by the overexpression of YTHDF, DNMT3B protein also increases, leading to DNA methylation levels increasing. Similarly, knocking down YTHDF1 also led to a decrease in DNA methylation levels in AGS and BGC-823 cell lines (Figure 5M). Collectively, YTHDF1 promotes DNMT3B protein expression by a direct interaction with DNMT3B mRNA, which resulting a crosstalk between RNA methylation and DNA methylation in GC.
3.5 YTHDF1 promotes the progression of GC by regulating DNMT3B
Studies have shown that both YTHDF1 and DNMT3B facilitate the tumorigenesis in GC (Wong et al., 2019; Chen et al., 2021; Bai et al., 2022). Consistent with previous reports, the overexpression of YTHDF1 and DNMT3B, respectively, promotes the proliferation of GC cells AGS and BGC-823. Simultaneously, knocking down YTHDF1 results in an inhibition of GC cells (Figures 6A, B). The efficiency of gene overexpression or knockdown is depicted on the right side of the figure. EdU incorporation assays were conducted in AGS and BGC-823 cells transfected with the vector, YTHDF1, DNMT3B, as well as two different siRNAs targeting YTHDF1 or DNMT3B, to assess DNA proliferation. Both YTHDF1 and DNMT3B promoted GC cell proliferation, and the cell proliferation promoted by YTHDF1 overexpression can be reversed by knocking down DNMT3B. On the contrary, knocking down YTHDF1 also attenuated the cell proliferation promotion ability of DNMT3B (Figure 6C).
Figure 6. YTHDF1 promotes the occurrence and development of gastric cancer by regulating DNMT3B. (A,B) Growth curve analysis was performed in AGS cells transfected with indicated genes; (C) EdU incorporation assays were performed in AGS and BGC-823 cells transfected with indicated genes; (D) Migration assays assay; (E) Invasion assays; (F) One step TUNEL apoptosis assays. Error bars represent the mean ± SD of three independent experiments. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001; two-tailed unpaired t-test.
To further validate the functions of YTHDF1 and DNMT3B in GC, we conducted migration and invasion experiments in AGS and BGC-823 cells. The findings demonstrate that overexpression of YTHDF1 and DNMT3B enhances the migration and invasion of GC cells. Moreover, the migration and invasion abilities promoted by YTHDF1 overexpression can be reversed by knocking down DNMT3B, and knocking down YTHDF1 also reduced the migration and invasion promotion by DNMT3B (Figures 6D, E). To investigate the effects of YTHDF1 and DNMT3B on cell apoptosis, we stained cells transfected with the vector, YTHDF1, DNMT3B, as well as two different siRNAs targeting YTHDF1 or DNMT3B with dUTP labeled with a red fluorescent probe Cyanine 3. The results showed that both YTHDF1 and DNMT3B can inhibit cell apoptosis, while knocking down YTHDF1 or DNMT3B can accelerate cell apoptosis caused by YTHDF1 and DNMT3B (Figure 6F). Thus, YTHDF1 promotes the progression of GC by upregulating DNMT3B, and inhibiting DNMT3B contributes to reduce the tumor promotion ability of YTHDF1.
4 Discussion
Cancer progression, driven by genetic and epigenetic aberrations, has received considerable attention. Our exploration of the role of m6A modification in GC has further underscored the functions of m6A methylation factors, both in the realm of genetic modulation and their interplay with the tumor microenvironment. Controlled by methylation regulators, m6A RNA modification patterns influence many key oncological processes, such as tumor proliferation, differentiation, and metastasis, with profound implications (Wang T. et al., 2020; Jiang et al., 2021). Analyzing m6A modification patterns in GC samples revealed distinct subgroups for diagnosis, guiding therapeutic strategies and ensuring personalized patient care.
Our findings highlight the critical role of immune cell infiltration in GC prognosis. The differences in prognosis between the high- and low-infiltration cohorts highlight the crucial influence of the TME on disease outcomes. This is congruent with current oncology paradigms that emphasize the significance of immune-tumor interactions. The identification of the “green” module via WGCNA showcases a direct association between the m6A RNA methylation cluster and immune infiltration patterns. This intersection of epigenetic regulation and immune dynamics underscores the holistic nature of cancer progression.
Establishing a risk prediction model based on a selected group of hub genes has potential clinical utility. Upon validation in larger and more diverse patient cohorts, this model could emerge as a reliable tool for risk stratification in patients with GC. The nine hub genes used to establish the risk prediction model are vital for cancer progression. ZFP64 enhances the activation of the p65 subunit, thereby promoting the production of pro-inflammatory and type-I interferons by Toll-like receptor-activated macrophages (Wang et al., 2013). Overexpression of ZFP64 promotes the proliferation of lung adenocarcinoma cells through activating the Notch pathway and is associated with poor prognosis (Jiang et al., 2020). ZFP64 functions as a transcription factor that promotes the expression of Galectin-1 (GAL-1), contributing to stem-cell-like properties and an immunosuppressive tumor environment. This activity enhances resistance to the chemotherapy drug nab-paclitaxel, playing a key role in the progression and chemoresistance of gastric cancer (Zhu et al., 2022). CD4+ T cells specifically induce the expression of mitochondrial TOMM34 (Gerner et al., 2019), and the role of TOMM34 in cancer cell growth suggests its potential in anti-cancer drug development or colorectal cancer diagnosis (Shimokawa et al., 2006). TOMM34 was identified as differentially expressed between intestinal-type and diffuse-type gastric cancer, suggesting it plays a role in the distinct molecular pathways of these cancer subtypes. This involvement may relate to processes important in cancer progression, such as adaptation to stress and resistance to therapy (Tanabe et al., 2020). DNMT3B is widely overexpressed in non-small cell lung cancer (NSCLC) and may be a potential molecular biomarker for personalized therapy (Samakoglu et al., 2012). DNMT3B influences tumor development through its enzymatic activity. Specifically, S-nitrosylation of DNMT3B reduces its enzymatic activity, leading to an abnormal upregulation of the Cyclin D2 gene (CCND2), which is necessary for the proliferation of certain tumor cells (Okuda et al., 2023). In gastric cancer, DNMT3B promotes tumor progression by methylating the MYH11 gene, thereby decreasing its expression and allowing the increase of TNFRSF14, which supports cancer development. This highlights DNMT3B as a potential target for cancer therapy (Wang et al., 2021). CSTF1, pivotal in DNA damage repair, is linked to increased breast cancer risk in BRCA2 mutation carriers due to CSTF1 mutations (Paolillo et al., 2015). In the study on gastrointestinal stromal tumors (GISTs), CSTF1 was involved in a fusion with Aurora kinase A (AURKA). This suggests that CSTF1, through this fusion, could play a role in the progression or behavior of GISTs, although specific mechanisms were not detailed (Denu et al., 2024). TM9SF4, primarily involved in cell adhesion and innate immunity, is overexpressed in a small subset of patients with metastatic melanoma, acute myeloid leukemia, and myelodysplastic syndromes (Paolillo et al., 2015). TM9SF4 was identified as a key gene in the regulatory network affecting response to cisplatin and fluorouracil treatment. Its specific role isn’t detailed, but its prominence in the network suggests it may influence mechanisms underlying chemoresistance (Sun et al., 2021). The interaction of Tti1 with mTOR in both mTORC1 and mTORC2 complexes regulates autophagy suppression (Kaizuka et al., 2010). ACTR5 has a pro-tumorigenic effect in neuroblastoma, and the knockdown of ACTR5 reduces cell proliferation and differentiation abilities (Veschi et al., 2017). STK35 regulates apoptosis and proliferation in osteosarcoma cells in osteosarcoma, exhibiting oncogenic properties (Wu et al., 2018). STK35 has been linked to immune signatures in gastric cancer, suggesting it may impact the immune response and effectiveness of immunotherapy. This connection highlights STK35 as a potential target for improving treatment outcomes (He and Wang, 2020). SS18L1 is associated with the occurrence and development of endometrial serous carcinoma (Saglam et al., 2020). SS18L1 has been identified as having copy number variations significantly linked to tumor metastasis. This association suggests that SS18L1 may influence the spread of gastric cancer, making it a potential marker or target for therapeutic strategies (Zhu et al., 2020). While these genes are implicated in tumor development and progression, their specific mechanisms in GC remain incompletely understood. Therefore, the nine hub genes identified through integration of m6A modification characteristics and the TME may play important roles in the prognostic assessment of patients with GC.
Furthermore, we found that DNMT3B was positively correlated with IGF2BP1 and YTHDF1, upon evaluating the correlation of these 37 hub genes with m6A methylation regulators and the proportion of immune-infiltrating cells. In a series of molecular experiments and cellular phenotypic validations, we demonstrated that DNMT3B and YTHDF1 cooperate to promote the proliferation, invasion, and metastasis of GC cells. YTHDF1 plays an important role in GC progression, and its functions and molecular mechanisms have been extensively investigated. High expression of YTHDF1 is associated with more aggressive tumor progression and poor prognosis in GC. Engineered small extracellular vesicles targeting YTHDF1 efficiently suppress GC progression and metastasis through epigenetic and immune modulation (You et al., 2023). The loss of YTHDF1 in gastric tumors potentiates the antitumor immune response by promoting the infiltration of mature dendritic cells (Bai et al., 2022). Elevated YTHDF1 expression also acts as a shield against the antitumor effects of chemotherapy and immunotherapy (Chen et al., 2022). Moreover, YTHDF1 overexpression holds clinical diagnostic significance across various cancers, including NSCLC, breast cancer, cervical cancer, GC, and colorectal cancer (Zhu Y. et al., 2023). YTHDF1 is significantly associated with metastatic gene signatures through ARHGEF2 translation and RhoA signaling activation in colorectal cancer (Wang et al., 2022). YTHDF1 directly targets p65 mRNA, promoting p65 protein overexpression without altering mRNA levels in Ythdf1-KO cells (Bao et al., 2023). YTHDF1 promotes cancer stem cell renewal and resistance to tyrosine kinase inhibitors in hepatocellular carcinoma (HCC), which enhances the stability and translation of m6A-modified NOTCH1 mRNA, leading to increased expression of NOTCH1 target genes. YTHDF1 drives HCC stemness and drug resistance, making it a potential therapeutic target for HCC treatment (Zhang et al., 2024). YTHDF1 promotes migration, invasion, and osteoblast adhesion and induces osteoclast differentiation of cancer cells in vitro and in vivo by inducing EZH2 and CDH11 translation (Wang et al., 2024).
In conclusion, our study highlights the intricate ties between m6A RNA methylation and TME dynamics in GC. As we move toward precision medicine, such insights will be pivotal in driving therapeutic innovations and improving patient outcomes. In addition, utilizing MeRIP and RIP experiments, we elucidated the molecular mechanism underlying the regulation of DNMT3B expression by m6A “reader”.
YTHDF1 and explored the crosstalk between m6A modification and 5mC modification in GC cells. Although our findings are promising, further investigations are essential to fully understand the mechanistic underpinnings and translate these insights into applicable clinical strategies.
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
MH: Data curation, Formal Analysis, Resources, Writing–original draft. MZ: Supervision, Writing–review and editing. JZ: Methodology, Writing–review and editing. YoW: Supervision, Writing–review and editing. TH: Validation, Writing–original draft. TM: Visualization, Writing–original draft. YiW: Software, Writing–original draft. BY: Visualization, Writing–review and editing. HQ: Supervision, Writing–review and editing. XT: Methodology, Writing–review and editing. HY: Resources, Writing–review and editing. WH: Funding acquisition, Investigation, Methodology, Project administration, Writing–review and editing. YaW: Funding acquisition, Investigation, Methodology, Project administration, Writing–review and editing.
Funding
The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by the National Natural Science Foundation of China [grant numbers 41931291, 42125707, 81974268, 82273403]; Major State Basic Research Development Program of China [grant number 2022YFA1103402]; Non-profit Central Research Institute Fund of Chinese Academy of Medical Sciences [grant number 2021-RC310-006]; and Chinese Academy of Medical Sciences Innovation Fund for Medical Sciences [grant number 2021-I2M-1-018].
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.
The reviewer LS declared a shared affiliation with the author YW to the handling editor at the time of review.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2024.1445321/full#supplementary-material
Abbreviations
DMEM, Dulbecco’s modified Eagle’s medium; GC, Gastric cancer; GS, Gene significance; ME, Module eigengene; MM, Module membership; NC, Negative control; NSCLC, Non-small cell lung cancer; PCA, Principal Component Analysis; TME, Tumor microenvironment; TCGA, The Cancer Genome Atlas-Stomach Adenocarcinoma; TOM, Topological Overlap Measure; WGCNA, Weighted Gene Co-expression Network Analysis.
References
Ali, H. R., Chlon, L., Pharoah, P. D., Markowetz, F., and Caldas, C. (2016). Patterns of immune infiltration in breast cancer and their clinical implications: a gene-expression-based retrospective study. PLoS Med. 13 (12), e1002194. doi:10.1371/journal.pmed.1002194
Bai, X., Wong, C. C., Pan, Y., Chen, H., Liu, W., Zhai, J., et al. (2022). Loss of YTHDF1 in gastric tumors restores sensitivity to antitumor immunity by recruiting mature dendritic cells. J. Immunother. Cancer 10 (2), e003663. doi:10.1136/jitc-2021-003663
Bao, Y., Zhai, J., Chen, H., Wong, C. C., Liang, C., Ding, Y., et al. (2023). Targeting m(6)A reader YTHDF1 augments antitumour immunity and boosts anti-PD-1 efficacy in colorectal cancer. Gut 72 (8), 1497–1509. doi:10.1136/gutjnl-2022-328845
Cancer Genome Atlas Research, N. (2014). Comprehensive molecular characterization of gastric adenocarcinoma. Nature 513 (7517), 202–209. doi:10.1038/nature13480
Chen, D., Cheung, H., Lau, H. C., Yu, J., and Wong, C. C. (2022). N(6)-Methyladenosine RNA-binding protein YTHDF1 in gastrointestinal cancers: function, molecular mechanism and clinical implication. Cancers (Basel) 14 (14), 3489. doi:10.3390/cancers14143489
Chen, X., Xu, M., Xu, X., Zeng, K., Liu, X., Pan, B., et al. (2020). METTL14-mediated N6-methyladenosine modification of SOX4 mRNA inhibits tumor metastasis in colorectal cancer. Mol. Cancer 19 (1), 106. doi:10.1186/s12943-020-01220-7
Chen, X. Y., Liang, R., Yi, Y. C., Fan, H. N., Chen, M., Zhang, J., et al. (2021). The m(6)A reader YTHDF1 facilitates the tumorigenesis and metastasis of gastric cancer via USP14 translation in an m(6)a-dependent manner. Front. Cell Dev. Biol. 9, 647702. doi:10.3389/fcell.2021.647702
Denu, R. A., Joseph, C. P., Urquiola, E. S., Byrd, P. S., Yang, R. K., Ratan, R., et al. (2024). Utility of clinical next generation sequencing tests in KIT/PDGFRA/SDH wild-type gastrointestinal stromal tumors. Cancers (Basel) 16 (9), 1707. doi:10.3390/cancers16091707
Dominissini, D., Moshitch-Moshkovitz, S., Schwartz, S., Salmon-Divon, M., Ungar, L., Osenberg, S., et al. (2012). Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature 485 (7397), 201–206. doi:10.1038/nature11112
Dong, Z., and Cui, H. (2020). The emerging roles of RNA modifications in glioblastoma. Cancers (Basel) 12 (3), 736. doi:10.3390/cancers12030736
Gagliardi, M., and Matarazzo, M. R. (2016). RIP: RNA immunoprecipitation. Methods Mol. Biol. 1480, 73–86. doi:10.1007/978-1-4939-6380-5_7
Gerner, M. C., Niederstaetter, L., Ziegler, L., Bileck, A., Slany, A., Janker, L., et al. (2019). Proteome analysis reveals distinct mitochondrial functions linked to interferon response patterns in activated CD4+ and CD8+ T cells. Front. Pharmacol. 10, 727. doi:10.3389/fphar.2019.00727
Han, X., Wang, M., Zhao, Y. L., Yang, Y., and Yang, Y. G. (2020). RNA methylations in human cancers. Semin. Cancer Biol. 75, 97–115. doi:10.1016/j.semcancer.2020.11.007
Hanzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7
Haruehanroengra, P., Zheng, Y. Y., Zhou, Y., Huang, Y., and Sheng, J. (2020). RNA modifications and cancer. RNA Biol. 17 (11), 1560–1575. doi:10.1080/15476286.2020.1722449
He, Y., and Wang, X. (2020). Identification of molecular features correlating with tumor immunity in gastric cancer by multi-omics data analysis. Ann. Transl. Med. 8 (17), 1050. doi:10.21037/atm-20-922
Huang, H., Weng, H., and Chen, J. (2020). m(6)A modification in coding and non-coding RNAs: roles and therapeutic implications in cancer. Cancer Cell 37 (3), 270–288. doi:10.1016/j.ccell.2020.02.004
Jiang, J., Zhang, J., Fu, K., and Zhang, T. (2020). Function and mechanism exploration of zinc finger protein 64 in lung adenocarcinoma cell growth and metastasis. J. Recept Signal Transduct. Res. 41, 457–465. doi:10.1080/10799893.2020.1825490
Jiang, X., Liu, B., Nie, Z., Duan, L., Xiong, Q., Jin, Z., et al. (2021). The role of m6A modification in the biological functions and diseases. Signal Transduct. Target Ther. 6 (1), 74. doi:10.1038/s41392-020-00450-x
Kaizuka, T., Hara, T., Oshiro, N., Kikkawa, U., Yonezawa, K., Takehana, K., et al. (2010). Tti1 and Tel2 are critical factors in mammalian target of rapamycin complex assembly. J. Biol. Chem. 285 (26), 20109–20116. doi:10.1074/jbc.M110.121699
Li, X., Ma, S., Deng, Y., Yi, P., and Yu, J. (2022). Targeting the RNA m(6)A modification for cancer immunotherapy. Mol. Cancer 21 (1), 76. doi:10.1186/s12943-022-01558-0
Liu, L., Li, H., Hu, D., Wang, Y., Shao, W., Zhong, J., et al. (2022). Insights into N6-methyladenosine and programmed cell death in cancer. Mol. Cancer 21 (1), 32. doi:10.1186/s12943-022-01508-w
Meyer, K. D., Saletore, Y., Zumbo, P., Elemento, O., Mason, C. E., and Jaffrey, S. R. (2012). Comprehensive analysis of mRNA methylation reveals enrichment in 3’ UTRs and near stop codons. Cell 149 (7), 1635–1646. doi:10.1016/j.cell.2012.05.003
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
Okuda, K., Nakahara, K., Ito, A., Iijima, Y., Nomura, R., Kumar, A., et al. (2023). Pivotal role for S-nitrosylation of DNA methyltransferase 3B in epigenetic regulation of tumorigenesis. Nat. Commun. 14 (1), 621. doi:10.1038/s41467-023-36232-6
Ontiveros, R. J., Stoute, J., and Liu, K. F. (2019). The chemical diversity of RNA modifications. Biochem. J. 476 (8), 1227–1245. doi:10.1042/BCJ20180445
Paolillo, R., Spinello, I., Quaranta, M. T., Pasquini, L., Pelosi, E., Lo Coco, F., et al. (2015). Human TM9SF4 is a new gene down-regulated by hypoxia and involved in cell adhesion of leukemic cells. PLoS One 10 (5), e0126968. doi:10.1371/journal.pone.0126968
Pinello, N., Sun, S., and Wong, J. J. (2018). Aberrant expression of enzymes regulating m(6)A mRNA methylation: implication in cancer. Cancer Biol. Med. 15 (4), 323–334. doi:10.20892/j.issn.2095-3941.2018.0365
Quail, D. F., and Joyce, J. A. (2013). Microenvironmental regulation of tumor progression and metastasis. Nat. Med. 19 (11), 1423–1437. doi:10.1038/nm.3394
Saglam, O., Tang, Z., Tang, G., Medeiros, L. J., and Toruner, G. A. (2020). KAT6A amplifications are associated with shorter progression-free survival and overall survival in patients with endometrial serous carcinoma. PLoS One 15 (9), e0238477. doi:10.1371/journal.pone.0238477
Samakoglu, S., Deevi, D. S., Li, H., Wang, S., Murphy, M., Bao, C., et al. (2012). Preclinical rationale for combining an EGFR antibody with cisplatin/gemcitabine for the treatment of NSCLC. Cancer Genomics Proteomics 9 (2), 77–92.
Shen, X. T., Xie, S. Z., Xu, J., Yang, L. Y., and Qin, L. X. (2022). Pan-cancer analysis reveals a distinct neutrophil extracellular trap-associated regulatory pattern. Front. Immunol. 13, 798022. doi:10.3389/fimmu.2022.798022
Shimokawa, T., Matsushima, S., Tsunoda, T., Tahara, H., Nakamura, Y., and Furukawa, Y. (2006). Identification of TOMM34, which shows elevated expression in the majority of human colon cancers, as a novel drug target. Int. J. Oncol. 29 (2), 381–386. doi:10.3892/ijo.29.2.381
Smyth, G. K. (2004). Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 3, Article3. doi:10.2202/1544-6115.1027
Sun, J., Zhao, J., Yang, Z., Zhou, Z., and Lu, P. (2021). Identification of gene signatures and potential therapeutic targets for acquired chemotherapy resistance in gastric cancer patients. J. Gastrointest. Oncol. 12 (2), 407–422. doi:10.21037/jgo-21-81
Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660
Tanabe, S., Quader, S., Ono, R., Cabral, H., Aoyagi, K., Hirose, A., et al. (2020). Molecular network profiling in intestinal- and diffuse-type gastric cancer. Cancers (Basel) 12 (12), 3833. doi:10.3390/cancers12123833
Tong, J., Cao, G., Zhang, T., Sefik, E., Amezcua Vesely, M. C., Broughton, J. P., et al. (2018). m(6)A mRNA methylation sustains Treg suppressive functions. Cell Res. 28 (2), 253–256. doi:10.1038/cr.2018.7
Veschi, V., Liu, Z., Voss, T. C., Ozbun, L., Gryder, B., Yan, C., et al. (2017). Epigenetic siRNA and chemical screens identify SETD8 inhibition as a therapeutic strategy for p53 activation in high-risk neuroblastoma. Cancer Cell 31 (1), 50–63. doi:10.1016/j.ccell.2016.12.002
Wang, C., Liu, X., Liu, Y., Zhang, Q., Yao, Z., Huang, B., et al. (2013). Zinc finger protein 64 promotes Toll-like receptor-triggered proinflammatory and type I interferon production in macrophages by enhancing p65 subunit activation. J. Biol. Chem. 288 (34), 24600–24608. doi:10.1074/jbc.M113.473397
Wang, J., Xu, P., Hao, Y., Yu, T., Liu, L., Song, Y., et al. (2021). Interaction between DNMT3B and MYH11 via hypermethylation regulates gastric cancer progression. BMC Cancer 21 (1), 914. doi:10.1186/s12885-021-08653-3
Wang, Q., Chen, C., Ding, Q., Zhao, Y., Wang, Z., Chen, J., et al. (2020a). METTL3-mediated m(6)A modification of HDGF mRNA promotes gastric cancer progression and has prognostic significance. Gut 69 (7), 1193–1205. doi:10.1136/gutjnl-2019-319639
Wang, S., Gao, S., Zeng, Y., Zhu, L., Mo, Y., Wong, C. C., et al. (2022). N6-Methyladenosine reader YTHDF1 promotes ARHGEF2 translation and RhoA signaling in colorectal cancer. Gastroenterology 162 (4), 1183–1196. doi:10.1053/j.gastro.2021.12.269
Wang, S., Xu, L., Wang, D., Zhao, S., Li, K., Ma, F., et al. (2024). YTHDF1 promotes the osteolytic bone metastasis of breast cancer via inducing EZH2 and CDH11 translation. Cancer Lett. 597, 217047. doi:10.1016/j.canlet.2024.217047
Wang, T., Kong, S., Tao, M., and Ju, S. (2020b). The potential role of RNA N6-methyladenosine in Cancer progression. Mol. Cancer 19 (1), 88. doi:10.1186/s12943-020-01204-7
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
Wong, C. C., Kang, W., Xu, J., Qian, Y., Luk, S. T. Y., Chen, H., et al. (2019). Prostaglandin E(2) induces DNA hypermethylation in gastric cancer in vitro and in vivo. Theranostics 9 (21), 6256–6268. doi:10.7150/thno.35766
Wu, Z., Liu, J., Hu, S., Zhu, Y., and Li, S. (2018). Serine/threonine kinase 35, a target gene of STAT3, regulates the proliferation and apoptosis of osteosarcoma cells. Cell Physiol. Biochem. 45 (2), 808–818. doi:10.1159/000487172
Yin, X., Wang, P., Yang, T., Li, G., Teng, X., Huang, W., et al. (2020). Identification of key modules and genes associated with breast cancer prognosis using WGCNA and ceRNA network analysis. Aging (Albany NY) 13 (2), 2519–2538. doi:10.18632/aging.202285
You, Q., Wang, F., Du, R., Pi, J., Wang, H., Huo, Y., et al. (2023). m(6) A reader YTHDF1-targeting engineered small extracellular vesicles for gastric cancer therapy via epigenetic and immune regulation. Adv. Mater 35 (8), e2204910. doi:10.1002/adma.202204910
Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020). m(6)A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol. Cancer 19 (1), 53. doi:10.1186/s12943-020-01170-0
Zhang, X., Su, T., Wu, Y., Cai, Y., Wang, L., Liang, C., et al. (2024). N6-Methyladenosine reader YTHDF1 promotes stemness and therapeutic resistance in hepatocellular carcinoma by enhancing NOTCH1 expression. Cancer Res. 84 (6), 827–840. doi:10.1158/0008-5472.CAN-23-1916
Zhou, Y., Zeng, P., Li, Y. H., Zhang, Z., and Cui, Q. (2016). SRAMP: prediction of mammalian N6-methyladenosine (m6A) sites based on sequence-derived features. Nucleic Acids Res. 44 (10), e91. doi:10.1093/nar/gkw104
Zhu, M., Zhang, P., Yu, S., Tang, C., Wang, Y., Shen, Z., et al. (2022). Targeting ZFP64/GAL-1 axis promotes therapeutic effect of nab-paclitaxel and reverses immunosuppressive microenvironment in gastric cancer. J. Exp. Clin. Cancer Res. 41 (1), 14. doi:10.1186/s13046-021-02224-x
Zhu, X., Tang, H., Yang, M., and Yin, K. (2023a). N6-methyladenosine in macrophage function: a novel target for metabolic diseases. Trends Endocrinol. Metab. 34 (2), 66–84. doi:10.1016/j.tem.2022.12.006
Zhu, Y., Li, J., Yang, H., Yang, X., Zhang, Y., Yu, X., et al. (2023b). The potential role of m6A reader YTHDF1 as diagnostic biomarker and the signaling pathways in tumorigenesis and metastasis in pan-cancer. Cell Death Discov. 9 (1), 34. doi:10.1038/s41420-023-01321-4
Keywords: STAD, m6A modification, TME infiltration, WGCNA, LASSO
Citation: Huo M, Zhang M, Zhang J, Wang Y, Hu T, Ma T, Wang Y, Yuan B, Qin H, Teng X, Yu H, Huang W and Wang Y (2024) Prognostic analysis of patients with gastric cancer based on N6-methyladenosine modification patterns and tumor microenvironment characterization. Front. Pharmacol. 15:1445321. doi: 10.3389/fphar.2024.1445321
Received: 07 June 2024; Accepted: 24 July 2024;
Published: 09 August 2024.
Edited by:
Yichao Hou, Shanghai Jiaotong University School of Medicine, ChinaReviewed by:
Jing Wu, Anhui University of Science and Technology, ChinaLei Shi, Tianjin Medical University, China
Jing Nie, People’s Liberation Army General Hospital, China
Copyright © 2024 Huo, Zhang, Zhang, Wang, Hu, Ma, Wang, Yuan, Qin, Teng, Yu, Huang and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Wei Huang, weihuang@ccmu.edu.cn; Yan Wang, yanwang@cicams.ac.cn