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

MATERIALS AND METHODS

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 .

Unsupervised Hierarchical Clustering Reveals Distinct Characteristics of Different Clusters

We systematically collected a set of hypoxia-related genes (https://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_HYPOXIA) and a total of 23 m6A regulators, including 8 writers (METTL3, METTL14, WTAP, RBM15, RBM15B, ZC3H13, CBLL1, and VIRMA), 2 erasers (ALKBH5 and FTO), and 13 readers (IGF2BP1/2/3, YTHDF1/2/3, YTHDC1/2, FMR1, ELAVL1, HNRNPC, HNRNPA2B1, and LRPPRC) (22, 30). To group hypoxia clusters, we performed principal component analysis (PCA) for data reduction. According to the Kaiser–Harris criterion, principal components <1% were considered noise and removed. After calculating the Euclidean distance, the ten-combined cohort was grouped using unsupervised hierarchical clustering with the “ward. D2” linkage criterion. Target genes for increasing oxygen delivery and reducing oxygen consumption were obtained from the Hypoxia-inducible factor 1 (HIF-1) signaling pathway to explore the hypoxic status between different clusters. Similarly, the hierarchical clustering method using the “ward.D” linkage criterion divides patients into high, medium, and low clusters according to m6A regulators. Furthermore, the results were visualized in clustering dendrograms, PCA, and t-distributed stochastic neighbor embedding (t-SNE) figures, and Kaplan Meier (KM) curves were employed to show the trends of overall survival (OS) and recurrence-free survival (RFS). Additionally, we explored the connections between the subtypes defined above and previous molecular stratifications of GC via a percentage stacking diagram (3, 4, 31).

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

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 ( ). 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 = Σ 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).

Cell Culture and Quantitative Real-Time PCR

The human gastric epithelial cell line GES-1 and GC cell lines HGC-27, NCI-N87 were obtained from Shanghai Anwei Biotechnology Co., LTD, China. GES-1 and NCI-N87 were cultured in an Roswell Park Memorial Institute (RPMI) 1640 (Gibco, US) medium, and HGC-27 was cultured in a Dulbecco's modified eagle medium (DMEM) (Gibco, US) medium with 10% fetal bovine serum (FBS; Gibco) and 1% penicillin/streptomycin in a humidified atmosphere of 5% CO2 and 20% O2 at 37°C. Total RNA was extracted using FastPure Cell/Tissue Total RNA Isolation Kit V2 (Vazyme, RC112-01) and reverse-transcribed into cDNA using HiScript III All-in-one RT SuperMix Perfect for qPCR (Vazyme, R333-01). Real-time PCR was performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Q711-02). The primer pairs used in qRT‐PCR were as follows: APOD 5′‐AATCGAAGGTGAAGCCACCC‐3′ (forward) and 5′-GTGCCGATGGCATAAACCAG‐3′ (reverse); CCN3 5′‐AGGCAGAGTTTCAGTGCTCC‐3′ (forward) and 5′‐TGCAGGTCCCAATGACCATC‐3′ (reverse); DACT1 5′‐TTGAACTGTTTGAGGCGAAGAG‐3′ (forward) and 5′‐ACTGAACACCGAGTTAGAGGAAT‐3′ (reverse); EML1 5′‐CAGTTCTGCAACGATGACAGC‐3′ (forward) and 5′‐GCCGAACCACATCAGCTAGAG‐3′ (reverse); MMP23B 5′‐TAGGCTTCTACCCGATCAACC‐3′ (forward) and 5′‐CGCTGTCGTCGAAGTGGAT‐3′ (reverse); RBPMS2 5′‐AAGACAGCCTGTTGGTTTTGT‐3′ (forward) and 5′‐CGAATACCGTTCAGCGCATT‐3′ (reverse); TUBB6 5′‐TGGTGGACTTAGAGCCAGG‐3′ (forward) and 5′‐CCCTTTCGCCCAGTTGTTC‐3′ (reverse).

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 μl of lysates were taken and 20 μl of 5×SDS-PAGE loading buffer were added. The samples were boiled for 15 min. Approximately 10 μl 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-3-phosphate dehydrogenase (GAPDH) (1:10,000), and HIF-α (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).

RESULTS

The Hypoxia Status in Gastric Cancer

The flow chart of this study is shown in ( ). After unsupervised hierarchical clustering, we classified 3 clusters with distinct hypoxia statuses ( ). 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 hypoxiaCluster-high patients ( ). The expression levels of target genes involved in increased oxygen delivery and reduced oxygen consumption varied among the clusters ( ), 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 ( ). 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 ( ). 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.

Figure 1

Flow chart.

Figure 2

Three hypoxia types with distinct prognosis characteristic and oxygen transport status. (A) Identification of hypoxiaClusters by unsupervised hierarchical clustering analysis. PCA and t-SNE analysis supported to divide patients into 3 hypoxiaClusters. (B, C) Kaplan–Meier curves were plotted to demonstrate the difference of prognosis by OS and RFS. (D, E) 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.

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 ( ). 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 ( ). 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 ( ). We conducted GSVA enrichment analysis to investigate the biological behaviors among these distinct hypoxia clusters. As shown in , hypoxiaCluster-high was dramatically enriched in stromal and metastatic activation pathways, such as EMT, angiogenesis, myogenesis, hedgehog signaling, and TNFα signaling via NFκB; 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 ( ; ). 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 ( ). Thus, we conjectured that stromal activation in hypoxiaCluster-high inhibited the antitumor effect of immune cells. In addition, we found that hypoxiaCluster-high 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) ( ). Pearson’s correlation analysis suggested that the abundances of these two categories of immunocytes were significantly positively associated in the TME ( ). 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.

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.

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 ( ) based on the mRNA expression levels of 21 regulators that presented high heterogeneity ( ). 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 ( ). Patients with invasive subtypes, EMT subtypes, or hypoxiaCluster-high subtypes were also mainly enriched in m6Acluster high ( ). 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 ( ), 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 ( ). 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 ( ). 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 ( ). In addition, patients with invasive subtypes, EMT subtypes, or IV stage had the lowest m6A score compared to other corresponding molecular/histological subtypes ( ), 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 ( ). These results showed that the m6A score could also be used to evaluate certain clinical features and was closely linked to hypoxia status.

Figure 4

Comprehensive analysis of distinct hypoxia status between 3 m6A methylation modification patterns. (A) The heatmap of gene set variation analysis (GSVA) was drawn on the basis of 25 immunity-related gene sets. m6Acluster high, leukocyte migration, myeloid leukocyte cytokine immune response, and mastocyte activation were mainly concentrated in the hypoxia-cluster high, which suggested hypoxia-cluster high might be in a immunosuppressed and stromal activation state. (B) The correlation analysis between m6Ascore and the known signatures. Orange indicates an R value > 0, blue indicates an R value < 0. (C–F) Calculation of the m6Ascore in different cohorts and the correlation in different subtypes, including (C) Singapore patient cohort (GSE15459 and GSE34942), (D) ACRG cohort (GSE62254), (E) hypoxiaCluster, (F) m6Acluster, (G) Stage. The upper and lower ends of the boxes meant the interquartile range of values. The lines in the boxes represented the median value. MSI, microsatellite instability. The Kruskal–Wallis test was used to compare the statistical difference between three gene clusters.

Construction and Validation of the Prognostic Signature of m6A-Related Hypoxia Pathway Genes

In total, 8 candidate predictive genes were identified ( ; ). 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 ( ). 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 entire cohort served as the validation cohort ( ). The prognostic accuracy of the risk score in the entire set was assessed; the areas under the ROC curve (AUCs) were 0.66 (62.94%–69.17%, 95% CI), 0.68 (65.18%–71.42%, 95% CI), and 0.67 (63.74%–70.74%, 95% CI) at 3, 5, and 7 years, respectively ( ). The MRHPPGs and HIF-1α expression was significantly elevated in gastric cancer. ( ) For external validation, the prognostic signature also showed a robust predictive ability ( ). Moreover, we used the clinicopathological variables and risk score to establish a nomogram quantifying the risk assessment ( ). The predicted AUC values were 0.811 and 0.727 in the entire cohort and GSE28541, respectively ( and ). The calibration curves presented a high credibility of the nomogram ( ). 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 ( ). Sankey diagrams clearly depicted that a high risk score was robustly related to other stratification classes with poor prognosis ( ). 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 ( , ). 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 ( and ).

Figure 5

Construction and validation of signature and nomogram. (A) The best candidate predictive genes were selected to construct the signature according to the frequency of genes by Lasso-bootstrapping that was defined more than 600 times in 1,000 times replications. (B) The coefficient value of the 8 selected genes came from the stepwise multivariate Cox regression analysis. (C) The difference of core biological processes was evaluated between high- and low-risk patients that were divided according to the optimal cut-off value. The lines in the boxes represented the median value. (D) The area under the ROC curves (AUC) of the entire data were 0.66 (62.94%–69.17%, 95%CI), 0.68 (65.18%–71.42%, 95%CI), and 0.67 (63.74%–70.74%, 95%CI) at 3, 5, and 7 years, respectively. (E) qRT-PCR results of MRHPPG expression level in GSE-1 and GC cell lines. (Data are presented as mean ± SD., **P < 0.01). (F) Western blot analyses of HIF-α protein levels in total cell lysates from paired clinical specimens of normal (N) and tumor (T) tissues from 12 patients with GC. (G) A survival nomogram of 3-, 5-, and 7-year OS was drawn in view of the risk score and stage. (H) AUC for risk score, stage, and nomogram reached at 0.659, 0.779, and 0.811, respectively, which mean a high accuracy to be a reliable method. (I) The calibration curve of nomogram showed a favorable result in 3, 5, and 7 years. (***P < 0.001,****p < 0.0001).

Table 1

Information of 8 m6A-related hypoxia pathway genes.

id	coef	HR	HR.95L	HR.95H	pvalue
APOD	0.091593	1.095918	0.953976	1.258979	0.195593
CCN3	0.11144	1.117887	0.992896	1.258612	0.065457
DACT1	0.211698	1.235775	1.074642	1.421068	0.002979
EML1	-0.41811	0.658292	0.550826	0.786723	4.27E-06
MMP23B	0.098619	1.103645	0.983001	1.239096	0.094981
RBPMS2	0.230006	1.258607	1.070618	1.479604	0.005325
TAC1	0.06395	1.066039	0.957504	1.186878	0.24309
TUBB6	0.224583	1.251801	1.066052	1.469914	0.006135

Figure 6

Validation of signature through two immune-checkpoint blockade treatment cohorts. (A, D) Kaplan–Meier curves were plotted in (A) IMvigor210 cohort and (D) Nathanson2017 cohort to confirm the credibility of the signature. (B, E) the AUC in IMvigor210 cohort and ICB.Nathanson2017 cohort were 0.648 and 0.781, respectively. (C, F) The response of immune-checkpoint blockade treatment could be connected with high- and low-risk scores, especially in anti-CTLA-4 (χ2 test). SD, stable disease; PD, progressive disease; CR, complete response; PR, partial response.

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 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 hypoxiaCluster-high 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-1α-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 immune-excluded 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 TGFβ-related 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 immune-excluded 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/ .

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.

Author Contributions

ZZ, HL and Z-KN conceived and designed this study. Z-KN and H-KT collected and assembled the data. C-GH and JL drafted the manuscript. Z-LY and H-NZ revised the manuscript. All authors read and approved the final version of the manuscript.

Funding

This study was supported by the National Natural Science Foundation of China (Grant Number: 81860433 and 82103645), Training Plan for Academic and Technical Young Leaders of Major Disciplines in Jiangxi Province (Grant Number: 20204BCJ23021), the Natural Science Youth Foundation of Jiangxi Province (Grant Numbers: 20192BAB215036), the Key Technology Research and Development Program of Jiangxi Province (Grant Number: 20202BBG73024), the Foundation for Fostering Young Scholar of Nanchang University (Grant Number: PY201822) and Science and Technology plan of Jiangxi Provincial Health and Family Planning Commission (Grant Number: 20195072) and Chinese Medicine Foundation of Jiangxi Provincial Health Commission Administration (Grant Number: 2018B038).

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.

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.