The Hypoxic Landscape Stratifies Gastric Cancer Into 3 Subtypes With Distinct M6a Methylation and Tumor Microenvironment Infiltration Characteristics

The interaction between hypoxia and RNA N6-methyladenosine (m6A) is an emerging focus of investigation. However, alterations in m6A modifications at distinct hypoxia levels remain uncharacterized in gastric cancer (GC). Unsupervised hierarchical clustering was performed to stratify samples into different clusters. Differentially expressed gene analysis, univariate Cox proportional hazards regression analysis, and hazard ratio calculations were used to establish an m6A score to quantify m6A regulator modification patterns. After using an algorithm integrating Least absolute shrinkage and selection operator (LASSO) and bootstrapping, we identified the best candidate predictive genes. Thence, we established an m6A-related hypoxia pathway gene prognostic signature and built a nomogram to evaluate its predictive ability. The area under the curve (AUC) value of the nomogram was 0.811, which was higher than that of the risk score (AUC=0.695) and stage (AUC=0.779), suggesting a high credibility of the nomogram. Furthermore, the clinical response of anti-PD-1/CTLA-4 immunotherapy between high- and low-risk patients showed a significant difference. Our study successfully explored a brand-new GC pathological classification based on hypoxia pathway genes and the quantification of m6A modification patterns. Comprehensive immune analysis and validation demonstrated that hypoxia clusters were reliable, and our signature could provide a new approach for clinical decision-making and immunotherapeutic strategies for GC patients.


INTRODUCTION
Gastric cancer (GC) is the fifth most malignant tumor worldwide (1). Greater than 1 million new cases have been identified, and most cases are already advanced at diagnosis, explaining why GC has the third highest number of cancer-related deaths (2). Based on the Lauren/WHO classification and the lymph node metastasis [tumor node metastasis (TNM)] staging of tumors, current treatments exhibit a poor correlation with the molecular pathology of cancer. Despite the development of new pathological classifications, such as The Cancer Genome Atlas (TCGA) subtypes and Asian Cancer Research Group (ACRG) subtypes, the clinical predictive value of these classification systems remains insufficient (3,4). To identify more molecular markers that are closely related to GC progression, accurately predicting developmental trends and providing individualized treatments for patients has become a troublesome point in current relevant research fields.
With greater insight into tumor research, changes in the tumor microenvironment (TME) have drawn more attention, and hypoxia plays an important role in tumorigenesis (5). Hypoxia is one of the characteristics of the microenvironment of solid tumors and one of the greatest obstacles to cancer treatment (6)(7)(8). As a master regulator of cellular adaptation to hypoxia, hypoxia-inducible factor 1 (HIF1) has been proven to extensively regulate the expression of hypoxia genes and hypoxia adaptation-related signal transduction pathways, including EPO, VEGF, iNOS, and other genes to increase oxygen transmission and PDK-1, ALDOA, bcl-2, and other genes to reduce oxygen consumption (5,(9)(10)(11)(12). Another feature of the TME is the change in immune cells, which contributes to maintaining a complex dynamic interaction with tumor cells (13,14). Immunotherapy, especially programmed cell death-1 (PD-1)/PD-1 ligand 1 (PD-L1), and immune checkpoint blockade (ICB), has made remarkable achievements in recent years (15,16). However, sustained clinical responses are only induced in a minority of cancer patients, indicating that more studies on this topic should be performed (17,18).
As the most common RNA modification in eukaryotic cells, N6-methyladenosine (m6A) not only plays a related role in immune regulation but also plays a vital role in the occurrence and development of cancer through various processes, such as proliferation, migration, and invasion (19,20). m6A regulators consist of three types of proteins: "writers" with methyltransferase activity, "erasers" with demethylase activity, and "readers" with m6A binding sites (20)(21)(22). Recent studies have demonstrated that the abnormal m6A modification patterns change the TME and lead to tumor progression, and hypoxia plays a potential role (23,24). Recently, several posttranscriptional modification databases have been established such as the m6AVar and RMBase databases (25,26), which provided important information about m6A-related variants to explore the molecular mechanisms of m6A modification for experimental biologists. Moreover, 2 powerful m6A functional analysis tools ConsRM and m6A2Target (27,28) were also developed. However, the specific mechanisms in GC remain elusive, so a comprehensive analysis of hypoxia and m6A is urgently needed and indispensable.
In this study, we identified three hypoxia pathway subtypes in GC. By correlating hypoxia with m6A modification patterns and defining the m6A score to quantify m6A modification patterns, we ultimately established a robust signature and prognostic nomogram. This study provides information on clinicopathological characteristics and a classification system that are more in line with reality and can be used to guide clinical decision-making. In addition, this study aims to improve GC patient survival.

Data Collection and Preprocessing
GC patients with survival information were retrospectively collected from the Gene Expression Omnibus (GEO, http:// www.ncbi.nlm.nih.gov/geo/) and The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/), and GC samples without clinical data were excluded. In total, 1,673 patients from ten cohorts were enrolled, including The Cancer Genome Atlas-Stomach Adenocarcinoma (TCGA-STAD), GSE13861, GSE26899, GSE26901, GSE57303, ACRG Cohort (GSE62254), Singapore Patient Cohort (GSE15459 and GSE34942), and GSE84437 (GSE84426 and GSE84433). The TCGA-STAD cohort (FPKM normalized) was transformed into the transcripts per kilobase million (TPM) format. For microarray cohorts, the normalized matrix files with expression data and clinical information were directly downloaded and log2 transformed. The remaining cohorts except TCGA-STAD were merged into one cohort, and the "sva" R package was employed to remove batch effects (29). The predictive value of the nomogram was tested using an additional cohort GSE28541. In addition, two immune checkpoint blockade treatment cohorts (IMvigor210 for PD-1 treatment and Nathanson2017 for CTLA-4 treatment) were obtained, and the corresponding normalized data were utilized to determine whether the m6A-related hypoxia signature could be used to screen immunotherapy-sensitive patients. Details are provided in Supplementary Table S1.

Pathway Enrichment Analysis and Single-Sample Gene Set Enrichment Analysis
Using the HALLMARK gene set (downloaded from the MSigDB database v7.1) as the background pathway, gene set variation analysis (GSVA) was performed using the "GSVA" R package to show pathway differences in 3 hypoxiaClusters (32). Moreover, based on immunity-related gene sets reported in a previously published article (33), we employed GSVA enrichment analysis to investigate the distinct response patterns of hypoxiaClusters in innate and adaptive immunity.
Single-sample gene set enrichment analysis (ssGSEA) was performed for other gene sets obtained from previously published studies as follows: the biomarkers of biological processes according to Mariathasan et al. (34), hypoxia biomarkers, T-cell dysfunction, and immunotherapy resistance biomarkers and immunosuppressive cell signatures (Supplementary Table 2).

Immune Cell Infiltration Estimation
For each sample, the ESTIMATE algorithm was adopted to assess the tumor purity and population estimation of stromal and immune cells based on gene expression (35). Twenty-eight different immune cell infiltration patterns, including cells executing antitumor reactivity and cells delivering protumor suppression, were calculated from the gene sets reported in a previous study via ssGSEA (Supplementary Table 2). Furthermore, an additional 22 immune cells calculated by the CIBERSORT deconvolution algorithm, including neutrophils, eosinophils, mast cells, dendritic cells, macrophages, natural killer cells, regulatory T cells (Tregs), B cells, CD4+ T cells, CD8+ T cells, and plasma cells, were assessed to quantify infiltrating pattern heterogeneity (36).

Quantization of the Modification Pattern of m6A Regulators
Based on the m6A score construction method of Shen et al. (22), differentially expressed genes (DEGs) between 3 m6A clusters were extracted using the "limma" R package (37). Univariate Cox regression analysis was performed on DEGs with a p-value <0.05; then, screened genes were employed for the construction of the m6A score and normalized from -1 to 1 to reduce the effect of the gene expression value. Afterwards, we calculated the hazard ratio (HR) of all screened genes and divided them into two groups based on a cut-off score of HR=1, and the m6A score was defined as the difference value of the sums in each group. Tumor mutation burden (TMB) was calculated using "maftools" according to the somatic mutation data acquired from the TCGA database (38). Pearson correlation analysis was employed to reveal the correlation between m6A and TMB.
Subsequently, the distribution differences in somatic mutations between the low and high m6A score groups were analyzed and visualized using a waterfall diagram.

Establishment of an m6A-Related Hypoxia Signature by Machine Learning
Differential expression analysis was performed to screen differentially expressed hypoxiaCluster genes between different hypoxiaClusters. The m6A-related hypoxia pathway genes (MRHPPGs) were defined based on the following criteria: correlation r >0.5 and p-value <0.001 between DEHCGs and the m6A core. Data partitioning and standardization were processed using the classification and regression training (caret) package in the R(caret) package (39). First, the entire dataset was divided into a training and testing cohort at a ratio of 6:4 using a stratified sampling method for each cohort. Then, the function "preProcess" was used to standardize the training cohort and other cohorts based on the parameters calculated in the training cohort. For feature engineering, we used 80% of the samples randomly chosen from the original sample each time. After 1,000 bootstrapping replications, genes with a p-value <0.01 that appeared greater than 900 times in univariate Cox proportional hazards regression analysis were included in further analyses. Furthermore, an algorithm integrating LASSO and bootstrapping was used to identify the best candidate predictive genes (40,41). The optimal candidates that were repeated more than 600 times in 1,000 iterations were determined through 5 cross-validations. Then, stepwise multivariate Cox regression analysis was used to build a prognostic signature, which was determined as follows: risk score = S expression level of gene Xi × Cox coefficient of gene Xi. A survival decision tree was used to show the process of clinical decision-making. The predictive ability of the nomogram for 3-, 5-and 7-year OS was assessed, and external validation was performed using cohort GSE28541. A calibration curve was generated to compare the predicted survival rates with the observed survival rates (42).

Western Blotting
The tumor cells were lysed with RIPA buffer containing a protease inhibitor cocktail. These cells were kept on ice for approximately 1 h and vortexed every 15 min at 12,000 rpm and centrifuged for 15 min. The protein concentration in the lysate supernatant was measured by Bicinchoninic acid (BCA). The whole lysates were diluted to the same concentration, 80 ml of lysates were taken and 20 ml of 5×SDS-PAGE loading buffer were added. The samples were boiled for 15 min. Approximately 10 ml for each sample were loaded when running Sodium dodecyl sulfate-Polyacrylamide gel electrophoresis (SDS-PAGE). The protein was fractionated by 12.5% SDS-PAGE. The protein was transferred to the Polyvinylidene difluoride (PVDF) membrane at 300 V for 1 h in an ice bath. The membrane was blocked with 5% Milk-TBST for 2 h at room temperature. Then, the membrane was probed with primary Abs for glyceraldehyde-3phosphate dehydrogenase (GAPDH) (1:10,000), and HIF-a (1:1,000) overnight at 4°C. HRP-conjugated anti-rabbit IgG (1: 20,000) was used as a secondary Ab. After secondary Ab incubation for 2 h at room temperature, the membrane was washed for 5 times with PBST and then blotted with an Enhanced chemiluminescence (ECL) solution. The blots were imaged in the dark room with an imaging machine.

Statistical Analysis
All data processing was performed using the R 4.0.3 software. For two groups, statistical significance was estimated via unpaired Student's t-tests for normally distributed variables and Wilcoxon rank-sum tests for nonnormally distributed variables. For more than two groups, one-way ANOVA tests and Kruskal-Wallis tests were used (43). The cut-off values of continuous variables, such as OS, were determined using the "survminer" R package. The area under the curve (AUC) of time-dependent receiver operating characteristic (ROC) curves was visualized by the "timeROC" R package (44), and the ROC curve of the immune checkpoint blockade therapy response was assessed using the "pROC" R package (45). Differences with p < 0.05 were considered statistically significant (*p < 0.05; **p < 0.01: ***p < 0.001: ****p < 0.0001).

The Hypoxia Status in Gastric Cancer
The flow chart of this study is shown in (Figure 1). After unsupervised hierarchical clustering, we classified 3 clusters with distinct hypoxia statuses ( Figure 2A). Next, we evaluated how the hypoxia status affected patient prognosis. Both overall survival (OS) and recurrence-free survival (RFS) prognostic analysis for the three major hypoxia statuses demonstrated a particularly prominent survival disadvantage in hypoxiaClusterhigh patients (Figures 2B, C). The expression levels of target genes involved in increased oxygen delivery and reduced oxygen consumption varied among the clusters ( Figures 2D-F), confirming that the different hypoxia clusters exhibited distinct hypoxia statuses. Patients with invasive and Epithelial-Mesenchymal Transition (EMT) subtypes were classified as hypoxiaCluster-high, whereas proliferative and TP53-negative subtypes were classified as hypoxiaCluster-low. We still observed that cancers classified as hypoxiaCluster-high exhibited poorer differentiation and were enriched in the diffuse subtype ( Supplementary Figures 1A-D). In GC, the EMT molecular subtype and diffuse histological type were closely related to a shorter OS. Our hypoxiaCluster classification was consistent with other hypoxia characteristics (Supplementary Figure 1E). These results suggested that there were different hypoxia statuses with a significant prognostic value and that the GC characterized by a high hypoxia state was closely correlated with high malignancy and rapid tumor progression.

TME Landscape in GC Tumors With Distinct Hypoxia Statuses
A significant difference was found: hypoxiaCluster-high, which had the worst outcome, had the highest stromal score, immune score, and ESTIMATE score but had the lowest tumor purity ( Figure 3A). Moreover, the log-rank test revealed that patients with high stromal scores, low immune scores, high ESTIMATE scores, or low tumor purity had a poor prognosis (Supplementary Figures 1I, J). These results suggested that hypoxiaCluster-high might be in a stroma activation state, which is associated with a worse outcome (46). Moreover, hypoxiaCluster-high was prominently associated with high T-cell suppression and exhaustion (Supplementary Figure  1F). We conducted GSVA enrichment analysis to investigate the biological behaviors among these distinct hypoxia clusters. As shown in Supplementary Figure 1G, hypoxiaCluster-high was dramatically enriched in stromal and metastatic activation pathways, such as EMT, angiogenesis, myogenesis, hedgehog signaling, and TNFa signaling via NFkB; hypoxiaCluster-low exhibited enrichment signaling pathways associated with MYC targets V2, MYC targets V1, E2F targets, and the G2 M checkpoint. To our surprise, the subsequent analysis of infiltrating immunocyte populations suggested that hypoxiaCluster-high was significantly enriched in innate immunocytes, including natural killer cells, macrophages, mast cells, MDSCs, and plasmacytoid dendritic cells ( Figure 3B; Supplementary Figure 1H). A previous research reported that the immune-excluded phenotype also exhibited the presence of a great number of immunocytes, but the immunocytes remained in the matrix around the nest of tumor cells rather than penetrating the parenchyma. Stromal activation in the TME is considered to promote T-cell inhibition (15). Moreover, hypoxiaCluster-high was prominently associated with high T-cell suppression and exhaustion (Supplementary Figure 3B). Thus, we conjectured that stromal activation in hypoxiaCluster-high inhibited the antitumor effect of immune cells. In addition, we found that hypoxiaClusterhigh had an increased abundance of immune cell infiltration, including cells performing antitumor functions (e.g., effector memory CD4 T cells, effector memory CD8 T cells, natural killer cells, natural killer T cells, and type 1 T helper cells) and cells executing protumor suppression (e.g., immature dendritic cells, macrophages, MDSCs, neutrophils, plasmacytoid dendritic cells, regulatory T cells, and type 2 T helper cells) ( Figure 3B). Pearson's correlation analysis suggested that the abundances of these two categories of immunocytes were significantly positively associated in the TME ( Figure 3C). This finding indicated the existence of a feedback mechanism in which the antitumor immune response could promote the recruitment or differentiation of cells specialized for immunosuppression. Based on the above inference, we were surprised to confirm that the three hypoxia clusters had dramatically distinct TME cell infiltration features.

The m6A Methylation Modification Patterns Are Distinct Between Hypoxic Conditions
It is generally accepted that m6A methylation modification is involved in diverse biological processes, including dysregulated cell death and proliferation, the degree of tumor malignancy, and immune modifications. Therefore, we similarly classified three m6A methylation modification patterns using the same analysis of hierarchical clustering mentioned above (Supplementary Figure  2A) based on the mRNA expression levels of 21 regulators that presented high heterogeneity (Supplementary Figures 2D, E). We defined these patterns as m6Acluster high, medium, and low, respectively. The Kaplan-Meier survival analysis for the three m6Aclusters demonstrated that m6Acluster low presented a remarkable survival advantage (Supplementary Figures 2B, C). Patients with invasive subtypes, EMT subtypes, or hypoxiaClusterhigh subtypes were also mainly enriched in m6Acluster high (Supplementary Figures 2F, H). To further explore the biological functions affected by m6A modification phenotypes in distinct hypoxia statuses, we performed an unsupervised clustering algorithm based on hypoxia-related genes in the three m6A methylation modification patterns. Analysis indicated that patients with m6Acluster high were mainly concentrated in the hypoxiaCluster-high group ( Figure 4A), which confirmed again that hypoxiaCluster-high was significantly relevant to stromal activation. To further illustrate the potential biological process associated with m6A regulator modification subtypes, we established the m6A score and further tested the relation between the known signatures and the m6A score ( Figure 4B). We observed the distribution differences of somatic mutations between patients with high and low m6A score in the TCGA-STAD cohort. Patients with a low m6A score had more extensive TMB, and the Pearson correlation analysis confirmed that low-m6A-score tumors were significantly negatively related to tumor mutation burden (Supplementary Figure 3A, B). Moreover, the m6Acluster high group exhibited a significantly increased m6A score compared to the other clusters, while the m6Acluster low group showed the lowest m6A score ( Figure 4F). In addition, patients with invasive subtypes, EMT subtypes, or IV stage had the lowest m6A score compared to other corresponding molecular/ histological subtypes (Figures 4C, D, G), which was consistent with previous studies (46). More importantly, this is the first report that a high hypoxia status was associated with a significantly increased m6A score ( Figure 4E). These results showed that the m6A score could also be used to evaluate certain clinical features and was closely linked to hypoxia status.

Construction and Validation of the Prognostic Signature of m6A-Related Hypoxia Pathway Genes
In total, 8 candidate predictive genes were identified ( Figures 5A, B; Table 1). The formula for the risk score is as follows: risk score = 0.091592596 × APOD expression + 0.111440156 × CCN3 expression + 0.211698352 × DACT1 expression + (-0.418107247) × EML1 expression + 0.098618533 × MMP23B expression + 0.230005594 × RBPMS2 expression + 0.06395013 × TAC1 expression + 0.224582981 × TUBB6 expression. The analyses for the biological processes indicated that high-risk scores were significantly associated with increased activation of stromal pathways but presented an immunosuppressive state with decreased immune checkpoints ( Figure 5C). KM curve analysis showed that patients with low risk scores had a better OS in the training cohort, which was consistent with the testing cohort, and the The oxygen transport status between different hypoxiaClusters was analyzed through increasing oxygen delivery and reducing oxygen consumption target genes and correlation. The lines in the boxes represented the median value. (F) Between the two groups, it was revealed that cluster high had the highest increase in oxygen delivery but the lowest reduction in oxygen consumption, while the cluster low had the opposite. PCA, principal component analysis; t-SNE, t-distributed stochastic neighbor embedding; OS, overall survival; RFS, recurrence-free survival. The asterisks represented the statistical p-value (*P < 0.05; **P < 0.01; ***P < 0.001; ****p < 0.0001). The Kruskal-Wallis test was used to compare the statistical difference between three gene clusters.  Figure 5D). The MRHPPGs and HIF-1a expression was significantly elevated in gastric cancer. (Figures 5E, F) For external validation, the prognostic signature also showed a robust predictive ability (Supplementary Figures 3C-G). Moreover, we used the clinicopathological variables and risk score to establish a nomogram quantifying the risk assessment ( Figure 5G). The predicted AUC values were 0.811 and 0.727 in the entire cohort and GSE28541, respectively ( Figure 5H and Supplementary Figure  3H). The calibration curves presented a high credibility of the nomogram ( Figure 5I). To better illustrate the prospect of the clinical application of the MRHPPG signature, a decision tree was used to visualize the stratification level, which displayed significant differences in survival (Supplementary Figure 4). Sankey diagrams clearly depicted that a high risk score was robustly related to other stratification classes with poor prognosis (Supplementary Figure 4). Next, we assessed the predictive value of the MRHPPG signature in the immunotherapeutic cohort. The AUC values of the IMvigor210 cohort and ICB.Nathanson2017 cohort in response to treatment were 0.648 and 0.781, respectively ( Figures 6A, B, D, E, Supplementary Figure 3I). Moreover, we particularly investigated the ability of the risk score to predict the efficacy of anti-PD-L1 and anti-CTLA-4 immunotherapy, suggesting that low-risk patients showed a higher response rate to immunotherapy compared with high-risk patients ( Figures 6C, F and Supplementary Figure 3J).

DISCUSSION
GC is a common malignant tumor with a high recurrence rate (1). Despite great advances in surgery, radiation, and chemotherapy over the past decade, the outcome for advanced GC remains poor (2). The TNM cancer staging system is currently the gold standard for the assessment of the prognosis of cancer, but this system does not consider gene heterogeneity (47). Moreover, since the prognosis of patients with GC varies greatly (48), the establishment of a robust classifier to stratify patients with precise prognosis prediction and risk stratification is urgently needed at present and is essential to

FIGURE 3 | Landscape of the TME between distinct hypoxia status in GC. (A)
The boxplots of the ESTIMATE method were used to explore the TME characteristics among these distinct hypoxia clusters, suggesting that hypoxia-cluster high had the highest stromal score, immune score, and ESTIMATE score but the lowest tumor purity (all P-values < 0.0001). The lines in the boxes represented median value. (B) The heatmap depicted the infiltrating difference of 28 immune cell types in 3 hypoxiaClusters. HypoxiaCluster-high had a higher abundance of immune cell infiltration, including cells performing an anti-tumor function (e.g., effector memory CD4 T cells, effector memory CD8 T cells, natural killer cells, natural killer T cells, and type 1 T helper cells) and cells executing pro-tumor suppression (e.g., immature dendritic cells, macrophages, MDSCs, neutrophils, plasmacytoid dendritic cells, regulatory T cells, and type 2 T helper cells). Moreover, hypoxia-cluster high was significantly rich in innate immunocyte infiltration including natural killer cells, macrophages, mast cells, MDSCs, and plasmacytoid dendritic cells. (C) The correlation between pro-tumor suppression and anti-tumor immunity was analyzed according to 3 hypoxiaClusters, respectively. Pearson's correlation analysis suggested that the abundances of these two categories of immunocytes have a significant positive association in TME (all P-values < 0.001). The asterisks represented the statistical p-value (****p < 0.0001). The Kruskal-Wallis test was used to compare the statistical difference between three gene clusters. maximize the benefits produced by individual therapy and a timely follow-up. A large amount of convincing evidence has suggested that various malignancies are the result of complex interactions between tumor cells and nonmalignant cells in the TME, including fibroblasts, myofibroblasts, endothelial cells, and immune cells (49), which collectively contribute to the formation of a particular niche that promotes tumor growth and metastasis (3,50,51). In this research, using the bioinformatics analysis technology, we identified three distinct hypoxia clusters, which exhibited significantly different m6A methylation modifications and immune cell infiltration by comprehensively mining the public transcriptional data of GC. Moreover, a risk score model was constructed on the basis of 8 MRHPPGs to predict the outcome of GC patients. These findings provide a new perspective for treatment strategies to improve the prognosis and risk stratification of patients by considering the TME characteristics and transcriptomics. As a tumor hallmark, hypoxia (reduced oxygen availability) is caused by an imbalance between increased oxygen consumption and insufficient oxygen supply, and the clinical significance of hypoxia has been widely reported in cancer therapy (52). Although the vigorous metabolism and rapid proliferation of cancer cells can stimulate the formation of a novel vasculature system that is disorderly, only a vascular system with accurate distribution in normal tissue can facilitate the delivery of oxygenated blood (53). Here, we identified three distinct hypoxia clusters with significantly different prognoses, including cluster low with the worst prognosis and cluster high with the best prognosis, using hierarchical clustering analysis based on principal components (removing principal components less than 1%). The immune-inflamed phenotype, which is also referred to as a hot tumor, exhibited abundant immunocyte infiltration in the TME (13,15). Although the immune-excluded phenotype also exhibited the presence of a great number of immunocytes, the immunocytes remained in the matrix around the nest of tumor cells rather than penetrating the parenchyma. The matrix was restricted to the capsule of the tumor, or it penetrated the tumor itself, illustrating that the immune cells were indeed inside the tumor (54). Consistent with the above statement, we revealed that hypoxiaCluster-high presented a significant stroma activation status and was substantially enriched in carcinogenic pathways, including EMT, angiogenesis, the Kirsten rat sarcoma viral oncogene homolog (KRAS) signaling pathway, myogenesis, and TGF beta signaling pathways, which are considered T-cell suppressive. Moreover, this finding further confirmed that the hypoxia cluster high was in an obvious T-cell exhaustion state. Hence, through adequately exploring the characteristics of TME cell infiltration induced by distinct hypoxia states, it was not surprising that hypoxiaClusterhigh had activated innate immunity but the poorest prognosis. m6A methylation is the most common intracellular modification and is ubiquitously present in eukaryotic mRNA (19). Accumulating evidence supports a close link between m6A regulators and hypoxic states. A recent study reported that tumor hypoxia leads to the epigenetic remodeling of m6A (55). Qing et al. (56) reported that HIF-1a-induced YTHDF1 expression was closely related to hypoxia-induced autophagy-related HCC progression. However, the biological function of m6A methylation modification in distinct hypoxia-induced immune states remains unknown. Herein, we defined three m6A subtypes with different clinical outcomes via the same analysis of hierarchical clustering, which further confirmed that m6A methylation dysregulation plays a critical role in the tumorigenesis and progression of various neoplasms. Specifically, m6Acluster high comprised the worst prognosis and was related to the highest hypoxic state; m6Acluster low was associated with the best prognosis and correlated with the lowest hypoxic state. By clarifying m6A gene signatures and establishing the scoring system, we could further precisely assess the effect of m6A modification patterns on GC. Patients with invasive, EMT, m6Acluster high, and IV stage subtypes were significantly associated with a higher m6A score, which demonstrated that the m6A score was a reliable and robust tool for comprehensively evaluating m6A modification patterns  and an independent prognostic biomarker for predicting patient survival in GC. Detailed associations between the m6A score and clinicopathological characteristics were found in our research. Our data also suggested a substantially negative relationship between the m6A score and tumor mutation burden (TMB). Moreover, we found that the high hypoxia cluster had a higher m6A score and that the low hypoxia cluster had a lower m6A score. Based on these results, we hypothesize that the immuneexcluded phenotype of GC patients was accompanied by the activation of the m6A-related hypoxia pathway and the acquisition of other biological abilities, such as EMT and angiogenesis. Previous studies reported that EMT-and TGFbrelated signaling pathway activation led to a weakened transport of T cells into tumors as well as decreased tumor cytotoxicity (34,57). Finally, our study focused on the MRHPPG signature that demonstrated a prognostic value. In the training group, we initially recognized 25 MRHPPGs correlated with a prognosis and established a prognostic signature comprising 8 MRHPPGs via multivariate Cox regression and bootstrap-based univariate analysis with LASSO. Kaplan-Meier analysis suggested that the overall survival of patients with low risk scores was better than that of patients with high risk scores. A dramatically distinct risk score existed between nonresponders and responders, suggesting that we could more accurately predict the GC patients' clinical response to anti-PD-1/CTLA-4 immunotherapy through the MRHPPG risk score. In addition, the analyses of the biological activity of the gene signature indicated that high risk scores were significantly related to lower CD8-positive effector T-cell activity, lower immune checkpoint responses, and higher EMT, further demonstrating that the activation of the m6A-related hypoxia pathway played an important role in immune states, especially in the immuneexcluded phenotype. Next, we built a nomogram to calculate a score representing the OS of GC patients. The calibration plot suggested that the model has a satisfactory fitting curve and better clinical application than the traditional staging system.
Several limitations in this research should be noted. First, several independent external validations were conducted in our research, but it was still difficult to include all of the diverse features of patients from different geographic regions when cases and materials were gathered retrospectively from public databases. Second, the microenvironment features of distinct tumor spatial regions might be different; however, the samples used for analysis were all from the tumor core. Additionally, our study was not completed enough to cover related bioinformatics analysis focusing on m6A RNA modification (e.g., databases like m6AVar and RMBase and functional tools like ConsRM and m6A2Target) (25,28,58). Therefore, further investigations based on well-designed, prospective, multicenter studies are required.

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.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Medical Ethics Committee of the Second Affiliated Hospital of Nanchang University. The patients/ participants provided their written informed consent to participate in this study.