Predictive value of N6-methyladenosine (m6A)-related genes for prognosis and correlation with tumor microenvironment in gastric cancer

To the Editor: Gastric cancer (GC) poses a great threat to human health worldwide. N6-methyladenosine (m6A) RNAmethylation, the most abundant post-transcriptional RNA modification, plays a crucial role in RNA splicing, export, processing, translation, and decay. The tumor microenvironment (TME), an aggregation of tumor cells and surrounding tumor-related non-tumor cells, innately modulates tumor progression. However, the prognostic value ofm6A regulator-related genes and their relationship with the TME have not been largely explored in GC. The aim of our study was to screen prognostic genes related to m6A modification in GC patients in The Cancer Genome Atlas (TCGA) database, establish an m6A regulatorassociated signature for predicting clinical outcome, and validate the performance in TCGA and Gene Expression Omnibus (GEO) cohorts. The relationship between m6A and immune infiltration was also assessed in GC. With the underlyingmechanismuncovered,m6A-relatedgenescould become predictive biomarkers for GC and provide deeper insights into the TME and immunotherapy.

The RNA expression profile and clinical characteristic data of patients with GC were downloaded from the TCGA (http://portal.gdc.cancer.gov) and GEO (https:// www.ncbi.nlm.nih.gov/geo/) databases. After excluding samples with overall survival (OS) <30 days, the data for 334 GC patients from the TCGA database and 433 patients in the GSE84437 dataset were downloaded on March 22, 2021. This study was approved by the Ethics Committee of China-Japan Friendship Hospital (No. 2018-116-K85-1). The procedures in this study were in accordance with the ethical standards of the Responsible Committee on Human Experimentation and with the Helsinki Declaration of 1975, as revised in 2000. An informed consent statement was not necessary because all data were available to the public. R version 4.0.3 and SPSS version 26.0 were utilized for data analysis. P < 0.05 was considered statistically significant.
Clinical characteristics, including age, gender, grade, tumornode-metastasis (TNM) stage, and the expression level of programmed cell death-ligand 1 (PD-L1), were compared between different clusters. The TME scores for each sample, including the immune score and stromal score, were inferred by the "ESTIMATE" R package to quantify the infiltration of immune and stromal cells. Based on RNA profiles, the "CIBERSORT" package was used to quantify the infiltration levels of 22 immune cells. TME scores and immune infiltration patterns were analyzed between the two clusters.
The TCGA cohort was randomly classified into training set and testing set. The prognostic signature was obtained using the least absolute shrinkage and selection operator (LASSO) regression analysis based on the expression levels of the prognostic m6A-associated genes in the TCGA training set. Kaplan-Meier curves were generated for survival analysis. Receiver operating characteristic curves were utilized to assess the predictive efficiency of the prognostic model according to the area under the curve (AUC). Univariate and multivariate Cox regression analyses were performed to identify independent prognostic factors. Furthermore, the GSE84437 dataset was utilized as an independent external dataset to validate the efficiency of the m6A regulator-associated signature.
Based on the 280 genes strongly correlated with the expression levels of the 23 m6A regulators in the TCGA cohort, 19 m6A regulator-related genes were identified to be associated with the prognosis of GC patients using univariate Cox analysis [Supplementary Figure 1A, http:// links.lww.com/CM9/B125]. Compared with normal tissue, only the expression levels of PJA2 and FBXL7 were decreased in GC, while the other 17 genes were upregulated in GC. Overall, m6A regulator-related genes might play crucial roles in GC carcinogenesis.
After performing consensus clustering, all 334 GC patients in the TCGA cohort were assigned to the following two clusters: cluster 1 (n = 204) and cluster 2 (n = 130). Compared to cluster 1, cluster 2 tended to include older patients with earlier TNM stages. In cluster 2, PJA2 and FBXL7 were downregulated, and the other 17 genes were upregulated. Moreover, cluster 2 was characterized by better OS, higher PD-L1 expression, lower immune score, lower stromal score [Supplementary Figure 1B To investigate the prognostic value of m6A-related genes, we randomly assigned the patients in the TCGA cohort to training set (n = 202) and testing set (n = 132) in a 6:4 ratio. The baseline characteristics of these two sets were not significantly different [Supplementary Table 1  The risk score was higher in younger patients ( 65 years), in cluster 1, in the high immune score group, and in the high stromal score group, respectively [Supplementary Figure 1K-1N The GC patients in the GSE84437 cohort were classified into high-risk and low-risk groups according to the median risk score. Patients with high-risk scores had a worse prognosis than those with low-risk scores. The AUC for this signature was 0.665. The high-risk score was associated with a high immune score, high stromal score [Supplementary Figure Table 2, http://links.lww.com/ CM9/B37]. Overall, these findings confirmed that the five-gene m6A-related signature was an applicable prognostic indicator and was associated with the TME score and immune cell infiltration levels in GC patients.
Accumulating evidence has revealed that m6A regulation might be involved in the occurrence and progression of GC. [3,4] In the present study, a five-gene prognostic signature was constructed based on the expression of m6A regulator-related genes for clinical outcome prediction. High-risk patients had a poorer prognosis, higher immune score, higher stromal score, and larger immunosuppressive cell infiltration. The performance of the signature was successfully validated in TCGA and GEO cohorts.
There are still some limitations in our study. First, largescale clinical samples in prospective cohorts should be used to validate the robustness of the five-gene signature. Moreover, the biological functions of the m6A regulator-associated genes should be confirmed by further in vitro and in vivo studies to explain the mechanisms involved in m6A modification and the alteration of the TME in GC as well as to assess their potential as immunotherapy targets in the future.