The human GC cell line BGC-823 was purchased from the Institute of Cell Research of the Chinese Academy of Sciences in Shanghai. The BGC-823 cells were cultured in DMEM (Procell Life Science&Technology Co.,Ltd) supplemented with with 10% fetal bovine serum (FBS; Gibco) and 1% penicillin/streptomycin (HuaXiao gene Technology) in a 37°C, 5% CO₂ incubator with a humidified atmosphere. Lentivirus-mediated METTL3-overexpressing and control were purchased from Jikai Gene Chemical Technology Co.,Ltd (Shanghai, China). Establishment and screening of METTL3-overexpressing stable cell models according to the instructions of Jikai Gene Technology Co., Ltd. Briefly, when the cells in the logarithmic growth phase, the control (NC group) and METTL3 overexpressing lentivirus were added in BGC-823 cells, puromycin (4 μg/mL) was used to screen the successfully transfected cells for the subsequent experiments. The western blotting analyses were performed as described previously (15), and the primary antibodies used in western blotting experiments are as follows: β-Actin (66009-1-Ig, Proteintech); AVEN (25846-1-AP, Proteintech); DNAJB1 (13174-1-AP, Proteintech); METTL3 (15073-1-AP, Proteintech).

Proteins were extracted from control and overexpressed METTL3-overexpressed BGC-823 cells, respectively. The tandem analysis of the high-precision mass spectrometer was performed by Novogene Technology Co., Ltd. (Beijing, China). The MS proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD036773. The results of protein quantitation were statistically analyzed by T-test. The proteins whose expression levels significantly differed between experimental and control groups (p<0.05 and fold-change ≥ 1.2) were defined as differentiallyexpressed proteins (DEPs) (16).

The total RNA of the control and METTL3 overexpression group were extracted by Trizol for m6A-tagged mRNA quantification. Then, m6A epitranscriptomic microarrays were performed by Aksomics, Inc (Shanghai, China), and data were acquired using the Agilent feature extraction software (version 11.0.1.1). We have uploaded the RAW data to the public repository Gene Expression Omnibus (GSE213499). Screening criteria for multiple variations and statistically significant (P-value) thresholds were used to identify differences between the two m6A-methylated RNA groups. Stratified clustering was performed to determine differences in m6A methylation patterns. For the analysis of the results of m6A methylation microarray, the fold-change ≥ 1.2, P< 0.05 was defined as DMGs (17).

The gene expression profiles of 804 GC patients were downloaded from the GSE84437 dataset of the TCGA and GEO (http://www.ncbi.nlm.nih.gov/geo/) databases, and the “limma” package in R was used to correct the data. Mutation data from GC patients is processed by a Perl script. The gene expression profile of a total of 804 patients in TCGA and GEO were both in FPKM format and log2 normalized. To further explore the biological significance of proteins and genes, we used the hypergeometric distribution of “clusterprofiler” package to analyze and visualize data for enrichment analysis of functions and pathways. In GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, the adjust P< 0.05 and FDR q< 0.25 are considered significant. The gene sets of oxidative phosphorylation, carbohydrate metabolism, amino acid metabolism and lipid metabolism were obtained from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/). Z-score is calculated using the “GOplot” package.

Protein-protein interactionPPI analysis of METTL3 was performed using STRING (https://string-db.org/) database. In addition, we use Cytoscape was applied (version3.8.2) to visualize the relationship between DEPs.

To explore the function of molecules regulated by METTL3, GC patients were c-ategorized into various clusters using the “ConsensusClusterPlus” package (50 iteratio-ns, 80% resampling rate Pearson correlation, https://www.bioconductor.org/packages/release/bioc/html/ConsensusClusterPlus.html).

Univariate and multivariate cox proportional risk analyses were conducted using the “survival” and “survminer” packages. The criterion of P < 0.05 was selected as the filtering threshold. The Least Absolute Shrinkage and Selection Operator (LASSO) was conducted to avoid over-fitting and to delete those tightly correlated genes. Tenfold cross-validation was employed to select the minimal penalty term (λ). The formula of the risk score was constructed as follows:

where Coef _i represents the coefficients and X _i represents the normalized count of each hub genes. The median risk assessment is achieved using “ggrisk” in R. Risk models were used to divide GC patients into high- and low-risk groups based on median risk scores. The “ggplot2” and “survival” packages were used to plot Kaplan-Meier (K-M) survival curves for the overall survival (OS) of different GC subtypes. The time-dependent subject work characteristic (ROC) curve and the area under the ROC curve (AUC) were used to assess the sensitivity and specificity of the model. The AUC value was calculated by the “survival ROC” packages. The “pRRophetic” package was used to predict clinical drug therapy for the different GC subtypes.

To investigate the global protein expression profiles in response to overexpression METTL3, mass spectrometry-based proteomic experiments with TMT labeling were performed by Novogene Biotech Co., Ltd (Beijing, China) to explore the comprehensive changes in proteomic profiles in BGC-823 cells (The efficiency of METTL3 overexpression was shown in Supplementary Figure 1 ). The flowchart of the study is as illutrated in Figure 1 . Principal component analysis (PCA) determined that the METTL3 overexpression and control samples were independently clustered and exhibited significant differences between the two groups ( Figure 2A ). Based on this, a volcano plot was constructed for visualization of DEPs between the two groups, and 173 differerntially expressed proteins were identified (fold-change ≥1.2, P<0.05, 79 up-regulated and 95 down-regulated proteins) following METTL3 overexpression ( Figure 2B ). Then, the top 20 with up- and down-regulated DEPs were selected to construct a heat map for visual analysis ( Figure 2C ). GO is one of the most extensively used tools for classifying biological entities into functionally related groups, and the GO results demonstrated that the downregulated proteins regulated by METTL3 overexpression were mainly enriched in cellular energy metabolism processes, such as ATP metabolism process, mitochondrial dehydrogenase complex, mitochondrial respiratory chain, and NADH dehydrogenase activity especially in intracellular OXPHOS ( Figure 2D , Supplementary Table 1 ). Correspondingly, upregulated proteins were mainly enriched in the pyrimidine nucleoside biosynthetic process, cellular metabolic compound salvage, cadherin binding and other biological processes ( Figure 2E , Supplementary Table 2 ). To explore the role of METTL3-regulated DEPs in OXPHOS, the 173 DEPs were intersected with the gene sets of the OXPHOS pathway obtained from GSEA, and a total of 13 intersecting molecules were screened to generated a heatmap to visualize their expression in GC cells ( Figure 2F ). Interestingly, these proteins were all significantly downregulated in METTL3-overexpressed GC cells compared to the control group. Finally, a PPI network diagram was constructed by using Cytoscape (version 3.8.2), and a strong correlation was identified between METTL3 and low-expressed DEPs ( Figure 2G ). Taken together, these results indicated that the down-regulated proteins regulated by METTL3 are involved in cellular energy metabolism, especially in mediating mitochondria-related oxidative phosphorylation.

To evaluate the impact of METTL3-regulated DEPs on the prognosis of GC patients, firstly, LASSO regression and multivariate Cox analysis were performed on 173 DEPs to determine the best candidate DEPs, and 7 (DCK, TMEM164, ANXA5, OPLAH, NOS3, DNPEP, EPB41L3) independent prognostic DEPs were identified for building risk models, and the coefficients for each DEP in the risk model were -0.03, -0.02, 0.00, -0.01, 0.09, -0.01 and 0.04 respectively ( Supplementary Table 3 ). Risk model was constructed based on the TCGA dataset), and the GSE84437 dataset (as the validation set) to verify the model values. The K-M survival curve showed that the overall survival (OS) of high-risk patients in the train group was significantly lower than that in the low-risk group ( Figure 3A , Supplementary Table 3 ). The ROC curve showed AUC values of 0.691, 0.624 and 0.590 for the probability of survival at 1, 3 and 5 years for GC patients ( Figure 3B ). In addition, as the risk score increased, the number of patients who died increases significantly ( Figures 3C, D ), and we also generated a heatmap to visualize the differential expression of 7 risk factors in high- and low-risk patients ( Figure 3E ). The risk score also significantly distinguished the prognosis of patients in the external independent GEO dataset ( Figures 3F–J ). Finally, gene set enrichment analysis (GSEA) of the 7 DEPs-related to the high- and low-risk groups was conducted to explore the potential mechanisms affecting the prognosis of GC patients. Interestingly, the high-risk group with increased expression levels of ANXA5 and TMEM164 was significantly enriched in energy metabolism-related signaling pathways such as mTORC1, glycolysis, and arachidonic acid metabolism ( Figures 3K-M ). Collectively, these results indicated that METTL3-regulated DEPs might influence the prognosis of GC patients through energy metabolism-related signaling pathways.

Previous studies have demonstrated that METTL3-mediated m6A methylation could regulate gene expression via post-transcriptional alterations (10). To explore the possible impact of m6A modification on the expression of METTL3-induced DEPs, m6A-mRNA & lncRNA epitranscriptomic microarrays were carried out in METTL3 overexpressing BGC-823 cells. The results showed that METTL3 overexpression resulted in a total of 1837 DMGs, with 430 up-regulated and 1407 down-regulated methylated genes ( Figure 4A ). Next, the top 20 methylated genes that were significantly up- and down-regulated were chosen to generate a heatmap ( Figure 4B ). The GO enrichment analysis showed that DMGs regulated by METTL3 were significantly enriched in pathways that responded to oxygen level, carbohydrates and the cellular response to glucose stimulation and so on ( Figure 4C ). KEGG analysis found that the DEGs were predominantly enriched in energy metabolism-related signaling pathways such as MAPK, FOXO, mTOR, and TGF-β ( Figure 4D , Supplementary Table 4 ).

To further validated whether METTL3-mediated m6A methylation is involved in regulating the expression of genes related to energy metabolism, we took the intersection of DEPs and DMGs obtained from proteomic and m6A methylation analyses,and a total of 24 core DEPs were obtained for further analysis ( Figure 4E ). The PPI network constructed with core DEPs displayed the interactions between these core DEPs (containing 24 nodes and 33 edges, P<0.05, Figure 4F ). Then, the relationship between core DEPs and oxidatively phosphorylated molecules was explored, Sankey diagram delineated the correlation between METTL3-regulated core DEPs and oxidatively phosphorylated molecules with correlation coefficients > 0.3 ( Figure 4G ). In short, these results suggest that METTL3-mediated m6A modification is involved in the oxidative phosphorylation pathway.

To further explore the effect of METTL3-regulated m6A modification and protein expression on cellular metabolism and its relationship with the prognosis of GC patients, the mRNA expression and profiling data were downloaded from the TCGA and GEO databases. A consensus clustering algorithm was used to divide GC patients into two subtypes according to the 24 core DEPs. To determine the optimal number of clusters, differences in 24 core DEPs expression clustering stability were assessed using the “Consensus Cluster Plus” package, and when the consensus matrix k=2, there was no crossover between the GC samples of cluster A and B ( Figure 5A , Supplementary Table 5 , Supplementary Figure 1 ). The histogram illustrated that the expression levels of 16 core DEPs (16 of 24 core DEPs above) were significantly different between the two subtypes ( Figure 5B ). Meanwhile, the PCA plot showed significant differences between the two clusters which reflected not only the specificity of the samples in clusters A and B, but also the classification accuracy of the samples ( Figure 5C ). Next, the correlation between the two subtypes and prognostic characteristics in GC patients was evaluated, and the result of the Kaplan-Meier survival curve showed that the OS of cluster B was significantly superior to that of cluster A in GC patients ( Figure 5D ). Finally, a Venn diagram was employed to visualize the relationship between DEGs in the two subtypes and the three major energy metabolism pathways including carbohydrate metabolism, amino acid metabolism, and lipid metabolism. The results showed that the DEGs overlapped with the three major energy metabolism pathways, and a total of 8 lipid metabolism-related, 8 carbon oxide metabolism-related and 6 amino acid metabolism-associated molecules appearing in the intersection ( Figure 5E ). GSEA enrichment analysis found that cluster A samples were significantly enriched in signaling pathways such as TGF-β and MAPK, whereas cluster B samples were significantly enriched in signaling pathways such as pentose phosphate pathway, pyrimidine metabolism, and drug metabolism ( Figure 5F ). Taken together, these results indicated that METTL3 might participate in metabolic processes and influence the prognosis of GC patients by mediating the level of m6A modification and protein expression, suggesting that these differentially expressed genes may serve as novel prognostic markers for GC patients.

Previous studies have reported that METTL3-dependent m6A methylation leads to RNA hypermethylation and participates in the regulation of downstream biological processes (18). We combined the up-regulated DEPs obtained by proteomic analysis and the hypermethylated genes obtained by m6A epitranscriptomic microarray analysis in GC cells, and a total of 5 core highly expressed and hypermethylated molecules were identified ( Figure 6A , Supplementary Table 6 ). Further comprehensive analysis based on TCGA data indicated that the expression level of 3 genes (AVEN, DAZAP2, DNAJB1) were significantly higher in GC tumor tissues than in normal tissues ( Figure 6B ). Next, to investigate whether AVEN, DAZAP2, and DNAJB1 could serve as potential prognostic biomarkers in GC patients by receiver operating characteristic (ROC) curves, the results showed that the AUC values of the three molecules were 0.559, 0.799, and 0.669, respectively, showing a predictive value for clinical prognosis of GC patients ( Figure 6C ). Notably, the boxplot results showed that the expression levels of AVEN and DAZAP2 were significantly correlated with clinicopathological grade and stage of GC patients ( Figure 6D-F ). The correlation between the expression of AVEN, DAZAP2, DNAJB1 and anticancer drugs in GC patients was analyzed, and the results showed that AVEN, DAZAP2, and DNAJB1 were significantly correlated with antitumor drugs (tanspiromycin, procarbazine, amonafide, megestrol acetate, etc.) ( Figure 6G ). Finally, AVEN and DNAJB1 (19, 20), which have DNA damage and tumor formation roles in tumors were selected to verify their protein expression levels in BGC-823 control and METTL3 overexpression cell lines, and the results showed that the AVEN and DNAJB1 level were significantly upregulated by METTL3 expression ( Figure 6H ). Collectively, the above results show that the high degree of m6A modification and increased protein expression mediated by METTL3 are significantly correlated with the prognosis and antitumor drugs of GC patients, signifying that they are promising novel biomarkers and potential clinical therapeutic targets for GC.