Identification of Novel Prognostic Biomarkers That are Associated with Immune Microenvironment Based on GABA-Related Molecular Subtypes in Gastric Cancer

Introduction

Gastric cancer (GC) is a type of malignancy with a worldwide prevalence of not less than one million in recent years, resulting in over 780,000 deaths, and is the third leading cause of cancer-associated death.^1,2 In recent years, the incidence of GC has decreased. This has been made possible by screening for Helicobacter pylori and improving diets. However, the incidence of GC is increasing in younger age groups.³ Despite advances in the current management of GC, including radiation, chemotherapy, and surgery, the 5-year overall survival of patients with advanced GC remains below 30%.^4,5 Poor prognosis and high mortality in GC are closely related to lymph node metastasis, delayed diagnosis, and limited treatment options.⁶ Therefore, the development of new prognostic markers is very helpful for early diagnosis and GC prognosis.

Peritoneal metastasis is common in patients with GC, and approximately 55% of patients with distant metastatic GC experienced it.^7,8 Once peritoneal metastases are diagnosed, patients have a survival time of only 4 months, which is significantly lower compared to GC patients without peritoneal metastases.⁹ Gamma-aminobutyric acid (GABA) is a neurotransmitter that is important for the balance between neuronal excitation and inhibition in the mammalian central nervous system.¹⁰ It has been shown that GABA deficiency is closely associated with several neurological diseases, including mania, Alzheimer’s disease, Parkinson’s disease depression, and Huntington’s chorea.^11–13 Recent studies revealed that GABA is involved in the progression of cancer. GABA is implicated in the migration, differentiation, and proliferation of cancer cells.¹⁴ Tumor-derived GABA is involved in activating β-catenin signaling and contributed to the growth of tumors.¹⁵ GABA is involved in prostate cancer progression via inducing gastrin-releasing peptide generation in neuroendocrine cells.¹⁶ GABA receptor is a potential therapeutic target and prognostic marker for non-small cell lung cancer.¹⁷ Many pathways have been associated with metastasis in cancer, including GABA-related pathways. GABA treatment increased the risk of metastasis in breast cancer.¹⁸ Nembutal is a GABA-receptor agonist, which has been to inhibit colon cancer metastasis.¹⁹ GABA signaling promoted breast cancer metastasis by mediating the ERK1/2 pathway.²⁰ However, there are few studies on GABA-related pathways involved in GC prognosis.

The advancement of gene sequencing technology and the accessibility of public databases have facilitated the gathering of gene expression data, thereby accelerating the discovery of prognostic biomarkers.^21,22 In this study, our team aimed to utilize bioinformatics methods to pinpoint diagnostic biomarkers for gastric cancer. Firstly, GSVA algorithm and Kaplan-Meier were used to assess the prognosis value of GABA-related pathways. Next, GC patients were divided into two subtypes based on the median score of Reactome GABA receptor activation (RGRA). Subsequently, GSEA, functional enrichment analysis, and immune infiltration analysis were performed between the two subgroups. Finally, differentially expressed analysis, and weighted gene co-expression network analysis (WGCNA) were used to identify hub GABA-related genes. The prognosis and expression of hub gene were then analyzed and validated in the TCGA database, GEO database, and clinical samples. Overall, we found that there were two subtypes of GC with different prognoses. To gain more insights into the relationship between subtypes and immunity, the study examined the differences in immune cell infiltration and signaling pathways between the two subtypes. Additionally, ATP1A2 was identified as a potential prognostic biomarker that showed a high correlation with immune cell infiltration in GC. These findings provide a theoretical basis for further studies on the role of ATP1A2 in GC carcinogenesis and progression.

Materials and Methods

Data Collection and Processing

The RNA sequencing expression profile and corresponding clinical characteristics of 353 GC patients were downloaded from The Cancer Genome Atlas (TCGA) database, and this dataset was used as the training cohort. The GSE13911, GSE62254, and GSE84437 datasets and corresponding clinical characteristics were downloaded from the Gene Expression Omnibus (GEO) database, and these datasets were used as the validation cohorts. Probe names were converted to the corresponding gene symbols and data batches were normalized before data analysis. Furthermore, subsequent analyses only included GC patients with complete survival data. A total of five GABA-related pathways and their corresponding genes were downloaded from Molecular Signature Database (https://www.gsea-msigdb.org/gsea/msigdb).

Identification of Prognosis-Related GABA Pathways in GC

Gene Set Variation Analysis (GSVA), is a method that unsupervised assesses the variation in pathway activity across a sample population. By pre-entering a selected set of genes, GSVA is used to assess potential changes in the activity of the pathway in each sample.²³ In the present study, the GSVA algorithm was applied to calculate the score of GABA-related gene sets in each GC sample. To evaluate the association between GABA-related pathways and overall survival in GC patients, a Kaplan-Meier survival analysis was then performed. A P-value less than 0.05 was considered a significant difference.

Differentially Expressed Genes (DEGs) Analysis

Based on the median score of Reactome GABA receptor activation (RGRA), GC patients were divided into low- and high-RGRA subgroups. RGRA-related genes (RRGs) between both subgroups, and DEGs between the normal and tumor groups were identified using the “limma” package. The thresholds were adjusted p < 0.05 and llogFCl≥ 0.5. The heatmap and volcano map was plotted using the ‘heatmap’ and ‘ggplot2’ packages.

Enrichment Analysis

We performed Gene Set Enrichment Analysis (GSEA) by uploading the gene expression profile data to GSEA software (https://www.gsea-msigdb.org/g2sea/index.jsp) to evaluate the primary signaling pathways relevant to these subgroups. In addition, we performed Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the DEGs using the “clusterProfiler” package. A term with p < 0.05 was considered statistically significant. The results of the top 10 terms were visualized using the “ggplot2” package.

Immune Infiltration Analysis Between the Subgroups

In the present study, a comprehensive assessment of immune infiltrate cell levels in each sample was performed using CIBERSORT, MCPCounter, and TIMER algorithms.^24–26 The results of were visualized using the “ggplot2” package. The ESTIMATE software was applied to compare the ESTIMATE score, Immune score, and Stromal score between the two subgroups. The R package “ggplot2” was used to perform Spearman correlation analysis.

Construction of Weighted Gene Co-Expression Network

We established a co-expression network based on the gene expression profiles of the TCGA dataset using the WGCNA tool.²⁷ The co-expression modules were then identified based on the dynamic cutting algorithm by applying “PickSoftThreshold” to test the optimal soft threshold. To test which important module had the highest correlation with the RGRA score, the correlation heatmap of modules and sample features was established. As hub genes, genes were selected with a module membership (MM) value > 0.8 and a gene significance (GS) value > 0.2.²⁸

Construction of Protein-Protein Interaction (PPI) Network

The interrelationship of differentially expressed RRGs (DERGs) was assessed using the Search Tool for the Retrieval of Interacting Genes (STRING) and the PPI network was constructed.²⁹ The visualization of the results of the PPI network was performed with the Cytoscape software (3.8.2). The core genes from the DERGs were identified by the cytoHubba plug-in of Cytoscape software.³⁰

Survival Analysis

Gene Expression Profiling Interaction Analysis (GEPIA) is a database that provides information on the survival results of genes using data from TCGA. Survival analyses were conducted using GEPIA for the association between key gene expression levels and prognostic outcomes.³¹ Kaplan-Meier plotter could assess the prognostic value of biomarkers for multiple tumors. The results of the survival analysis were validated using the Kaplan-Meier plotter and GEO datasets (GSE62254 and GSE84437).

Cell Culture

The human normal gastric cell line (GES-1) and GC cell lines (BGC-823, MGC-803, and BGC-803) were obtained from Shanghai, China. The cells were cultured in RPMI-1640 medium (Life technologies, Shanghai, China) and maintained at 37 °C in a humidified atmosphere containing 5% CO2.

Measurement of the mRNA Expression of ATP1A2 in GC

From 2021 to 2022, eight pairs of cancer and para-cancer samples were collected from GC patients undergoing surgery at the Department of Gastroenterology, Guangdong Provincial People’s Hospital. The informed consent was obtained from patients. Approval for this research was from the Ethics Committee of Guangdong Provincial People’s Hospital. TRIzol reagent (Pufei, Shanghai, China) was used to extract total RNA from the tissues. DNA was prepared using Transcriptor reverse transcriptase (ThermoFisher Scientific, USA). The gene expression level was quantified by SYBR green PCR Master Mix (ThermoFisher Scientific, USA) using a Real-Time PCR System (ThermoFisher Scientific, USA). Table S1 presented the primers used in this study.

Statistical Analysis

All analyses have been carried out using R software (3.4.1). Comparisons between the two groups were made using the Wilcoxon test. Statistical significance was defined as p < 0.05.

Results

Identification of Prognosis-Related GABA Pathways

Figure 1 showed the detailed analysis process of this study. First, the GSVA algorithm was applied to calculate the score of GABA-related gene sets in each GC sample. Then, the results of the Kaplan-Meier survival analysis indicated that a high level of RGRA was associated with a poor prognosis in GC patients (Figure 2B, p = 0.03). However, the other four GABA-related gene sets (Figure 2A and C–E) did not predict prognosis in GC patients.

Figure 1 Flow chart of the present study.

Figure 2 Identification of prognosis-related GABA pathways. Association of BIOCARTA GABA PATHWAT (A), REACTOME GABA RECEPTOR ACTIVATION (B), REACTOME GABA SYNTHESIS RELEASE REUPTAKE AND DEGRADATION (C), WP GABA METABOLISM AKA GHB (D), WP GABA RECEPTOR SIGNALING (E) with overall survival of GC.

GSEA Between the Low- and High-RGRA Subgroups

Based on the median score of RGRA, GC patients were divided into low- and high-RGRA subgroups. A total of 136 RRGs between the subgroups were identified, of which the top 50 genes were shown in a heatmap (Figure 3A). To explore the biological functions that were enriched in the two subgroups, GSEA was carried out. As shown in Figures 3B and C, protein-protein interactions at synapses (NES = 2.94, p.adj < 0.001), GABA B receptor activation (NES = 2.53, p.adj < 0.001), leukocyte transendothelial migration (NES = 2.16, p.adj < 0.001), CD8 TCR pathway (NES = 2.42, p.adj < 0.001), B cell receptor signaling pathway (NES = 2.33, p.adj < 0.001), T cell receptor signaling pathway (NES = 2.35, p.adj < 0.001), chemokine signaling pathway (NES = 2.84, p.adj < 0.001), fceri mediated NF KB activation (NES = 3.01, p.adj < 0.001), signaling by the B cell receptor BCR (NES = 3.06, p.adj < 0.001), and neuronal system (NES = 3.17, p.adj < 0.001) were differentially enriched in the two subgroups.

Figure 3 Identification of DEGs and GSEA. (A) The heatmap presented the top 50 DEGs between low- and high-RGRA subgroups. Bubble diagrams (B) and multipeaked maps (C) represented the GSEA results.

The Landscape of Immune Cell Infiltrations in the Two Subgroups

CIBERSORT, MCPCounter, and TIMER algorithms indicated there was a significant difference in the levels of immune cells between the two subgroups. Compared with the high-RGRA group, the low-RGRA group exhibited a lower proportion of B cell naïve (p < 0.01), T cells regulatory (p < 0.05), monocytes (p < 0.01), dendritic cells resting (p < 0.01), mast cells resting (p < 0.05), and eosinophils (p < 0.05), whereas the level of T cells CD4 memory activated was relatively higher (p < 0.001) (Figure 4A). Compared with the high-RGRA group, the low-RGRA group had a lower proportion of T cells (p < 0.001), CD8 T cells (p < 0.05), cytotoxic lymphocytes (p < 0.01), B lineage (p < 0.001), monocytic lineage (p < 0.001), myeloid dendritic cells (p < 0.001), neutrophils (p < 0.001), endothelial cells (p < 0.001), and fibroblasts (p < 0.001) (Figure 4B). Compared with the high-RGRA group, the low-RGRA group had a lower proportion of B cell (p < 0.001), T cell CD4 (p < 0.001), T cell CD8 (p < 0.001), neutrophil (p < 0.001), macrophage (p < 0.001), and dendritic cells (p < 0.001) (Figure 4C).

Figure 4 Characterization of tumor immune microenvironment between the low- and high-RGRA subgroups. CIBERSORT (A), MCPCounter (B), and TIMER (C) algorithms were used to assess the immune infiltrate cell levels. *p < 0.05, **p < 0.01, ***p < 0.001.

Identification of DERGs and Functional Enrichment Analysis

To identify differential genes associated with RGRA, we performed differential expression analysis between two RGRA-associated subgroups. A total of 136 RRGs were identified between the high- and low-RGRA subgroups (Figure 5A). Meanwhile, a total of 6619 DEGs were screened between the normal and tumor groups (Figure 5B). By using the Venn tool, we found 73 DEGs belonging to RRGs (Figure 5C). In addition, GO and KEGG pathways were carried out to further assess the function of 73 DEGs. Cell-cell signaling, regulation of system processes, and chemical synaptic transmission were significantly enriched in biological processes (BP) (Figure 5D). The neuron part, synapse, and somatodendritic compartment were remarkably enriched in cellular components (CC) (Figure 5E). Potassium ion transmembrane transporter activity, monovalent inorganic cation transmembrane transporter activity, and ion-gated channel activity were significantly enriched in molecular function (Figure 5F). Neuroactive ligand-receptor interaction, cAMP signaling pathway, and cell adhesion molecules were significantly enriched in the KEGG pathway (Figure 5G).

Figure 5 Identification of DERGs and functional enrichment analysis. (A) The volcano plot depicts the RRGs between the low- and high-RGRA subgroups. (B) The volcano plot depicts the DEGs between the normal and tumor groups. (C) Identification of intersection genes between RRGs and DEGs. The GO-BP (D), GO-CC (E), and GO-MF (F) results were presented by bubble diagrams. The KEGG (G) results were presented by loop diagrams.

Identification of Key Modules Related to RGRA

To further identify hub RGRA-associated genes, we performed WGCNA. Figure 6A and B showed that we have chosen β = 5 to ensure a scale-free network and perform WGCNA. Then, 64 modules were identified with different colors (Figure 6C). Among these modules, the bule (r = 0.6, p = 3.3e-38) and paleturquoise (r = 0.6, p = 1.5e-38) modules had the highest correlation with RGRA. Therefore, the genes in the two modules were the subject of a hub gene identification.

Figure 6 Identification of key module related to RGRA. (A and B) To ensure the scale-free distribution, the optimal soft threshold was selected. (C) Heatmap for the correlation between RGRA and modules.

Identification and Validation of Prognostic Genes

Figure 7A illustrated the PPI network diagram for 73 DERGs. To screen genes with a high degree of association within the key modules, the thresholds of MM > 0.8 and GS > 0.2 were used to identify an important gene (ATP1A2) for further analysis (Figure 7B). The survival analysis by GEPIA showed that the ATP1A2 expression was significantly associated with overall survival in GC patients (p = 0.012, Figure 7C). Compared with the normal group, the expression of ATP1A2 was significantly downregulated in the GC group (p < 0.05, Figure 7D). Furthermore, these results were further validated by external datasets and clinical samples. In the GSE13911 dataset (Figure 8A) and clinical samples (Figure 8B), ATP1A2 expression was significantly downregulated in the tumor group. In the Human Protein Atlas (HPA), the protein expression of ATPA2 was also significantly downregulated in the tumor group (Figure S1A and B). In the Kaplan-Meier plotter (p = 1.4e-06, Figure 8C), GSE62254 (p = 3e-04, Figure 8D), and GSE84437 (p = 0.0172, Figure 8E), ATP1A2 expression was associated with overall survival in GC patients.

Figure 7 Identification of prognostic genes. (A) PPI network of DERGs. (B) Overlap of top 10 DERGs and key module genes. Survival analysis (C) and expression analysis (D) of ATP1A2 genes in GC by GEPIA. *p < 0.05.

Figure 8 Validation of ATP1A2 gene by independent datasets. Validation of ATP1A2 expression by GSE13911 dataset (A) and clinical samples (B). Validation of ATP1A2 prognostic value by Kaplan-Meier plotter (C), GSE62254 (D), and GSE84437 (E). **p < 0.01, ***p < 0.001.

The relationship between ATP1A2 and clinical features was also investigated. ATP1A2 expression was lower in the over 65 than in the under 65 (p < 0.01, Figure S2A). ATP1A2 expression was lower in the G1&G2 grade than in the G3 grade (p < 0.001, Figure S2B). ATP1A2 expression was lower in stage I than in stage II and stage III (p < 0.01, Figure S2C). ATP1A2 expression was lower in the T1&T2 stage than in the T3&T4 stage (p < 0.01, Figure S2D). No statistical difference was found between ATP1A2 expression and gender (Figure S2E), pathologic N stage (Figure S2F), and pathologic M stage (Figure S2G). All these results showed that ATP1A2 was associated with clinical features in GC patients.

Relationship Between ATP1A2 and Immune Cell Infiltration

In the present study, the level of infiltrating immune cells in GC was assessed through GSEA and ESTIMATE algorithms. Compared with the high-ATP1A2 group, the low-ATP1A2 group exhibited a lower proportion of B cells, CD8 T cells, cytotoxic cells, DC, eosinophils, iDC, macrophages, mast cells, NK cells, pDC, T cells, Tcm, Tem, TFH, Tgd, and Th1 cells, whereas the level of Th17 cells and Th2 cells was relatively higher (Figure 9A). The expression of ATP1A2 was positively correlated with the level of B cells, CD8 T cells, cytotoxic cells, DC, eosinophils, iDC, macrophages, mast cells, NK cells, pDC, T cells, Tcm, Tem, TFH, Tgd, and Th1 cells; while negatively correlated with the level of Th17 cells and Th2 cells (Figure 9B). As shown in Figure 9C, compared with the high-ATP1A2 group, the low-ATP1A2 group exhibited a lower of the stromal score, immune score, and ESTIMATE score.

Figure 9 Relationship between ATP1A2 and immune cell infiltration. (A) The histogram represents the different levels of immune infiltration between the low- and high-ATP1A2 subgroups. (B) Correlation analysis between ATP1A2 expression and immune cell infiltration. (C) The histogram represents the different levels of Stromal score, Immune score, and ESTIMATE score between the low- and high-ATP1A2 subgroups. *p < 0.05, **p < 0.01, ***p < 0.001. Abbreviation: ns, not significant.

Validation of ATP1A2 Gene Expression by Cell Experiments

In order to verify the results of our bioinformatics analysis, we further conducted cell experiments. The qRT-PCR results indicated that the levels of ATP1A2 expression were down-regulated in AGS, BGC-823, and MGC-803 cells compared to GES-1 cells (p < 0.01) (Figure 10).

Figure 10 Analysis of ATP1A2 gene expression level in GC cell lines. **p < 0.01 and ***p < 0.001.

Discussion

In recent years, the treatment of GC is based on surgical resection with adjuvant treatments such as targeted therapy, gene therapy, chemotherapy, and radiotherapy, which has only led to a significant decrease in the incidence of gastric cancer.^5,32,33 Despite this, GC is still the second most common cause of cancer-related deaths in China.³⁴ Developing new prognostic markers and therapeutic targets will therefore help diagnose and predict GC. The origin of GABA remains unclear, but it has been found in breast, ovarian, stomach, and glioma cancer.³⁵ In addition, GABA also involves in the biological process of tumor cells including migration, differentiation, and proliferation.¹⁴ Moreover, GABA level is a predictor of mortality in patients with colon and lung cancer.¹⁵ Studies have shown that GABA accumulation promotes immune suppression and tumor cell proliferation and that inhibiting the GABA pathway may improve the effectiveness of immunotherapies against tumors.³⁶ Therefore, genes associated with the GABA pathway should receive attention, and we hypothesized that molecular subtypes based on the GABA pathway could guide the diagnosis of GC and improve its survival rate.

In the present study, the GSVA algorithm was used to calculate the score of GABA-related gene sets in each GC sample. The results indicated that a high level of RGRA was associated with a poor prognosis in GC patients. Subsequently, based on the median score of RGRA, GC patients were divided into low- and high-RGRA subgroups. GSEA results showed that GABA B receptor activation, B cell receptor signaling pathway, T cell receptor signaling pathway, and chemokine signaling pathway were differentially enriched in these two subgroups. A previous study has indicated that GABA participated in the proliferation of the gastric cancer cell line.³⁷ In GC, CXC chemokines and their receptors regulated tumor formation by regulating tumor metastasis, invasion, growth, and transformation, and indirectly by regulating tumor-leucocyte interactions and angiogenesis.³⁸ The immune microenvironment also plays an important role in the progression of chronic inflammation to GC.³⁹ In addition, CIBERSORT, MCPCounter, and TIMER algorithms were used to characterize the immune infiltration landscape between the two subtypes. Our results showed that the high-RGRA subtype had a worse prognosis and an activated immune microenvironment characterized by high infiltration levels of immune cells. Substantially reduced DC, macrophages, T cells, and CD8 T cells, especially in the high-RGRA subtype with a poor prognosis. Their results were consistent with previous studies. The prognosis of GC patients was influenced by DC subsets and CD4/CD8+ cells.⁴⁰ The M2 phenotype of tumor-associated macrophages could predict the prognosis in GC.⁴¹ The presence of CD8+ tumor-infiltrating lymphocytes was an independent prognostic marker for GC patients.⁴² Thus, the involvement of the GABA pathway in the pathogenesis of GC may be mediated through the regulation of immune responses between the two molecular subtypes.

ATP1A2 is a gene encoding the catalytic subunits of Na/K-ATPases that play a key role in the basic electrophysiology of neurons.⁴³ Recent literatures have shown that there is abnormal expression of ATP1A2 in some type of cancers. This abnormal expression is thought to play a role in the development and onset of cancer.^44,45 In addition, ATP1A2 suppressed invasion, migration, apoptosis, and proliferation in prostate cancer by inactivating the Smad/TGF-β pathway.⁴⁶ In our study, the results of multiple databases showed that ATP1A2 was down-regulated in GC patients. We also found that a higher level of ATP1A2 was associated with shorter overall survival time in GC patients, which was consistent with the effect of ATP1A2 on prognosis in ovarian cancer and bladder cancer.^45,47 In addition, this result was validated in other independent datasets. These results indicated that ATP1A2 has potential as a tumor promoter and its high expression was significantly associated with poor prognosis in patients with GC.

Tumor microenvironment is composed of immune cells, endothelial cells, fibroblasts, and stromal cells. It served as an indicator of patients’ prognosis and the effectiveness of immunotherapy.^48,49 A recent study has shown that a high immune score is an independent predictor of poor overall survival in patients with GC.⁵⁰ Compared to the high-ATP1A2 group, the low-ATP1A2 group displayed a decrease in stromal score, immune score, and ESTIMATE score. Our study presented a new angle on the relationship between ATP1A2 and the prognosis of GC patients. It proposed that the correlation between ATP1A2 and immune score implied that ATP1A2 expression could impact the infiltration of immune cells. Additionally, another important result of the present study was that ATP1A2 expression was positively correlated with the level of immune cells, including CD8 T cells, cytotoxic cells, plasmacytoid dendritic cells (pDC), eosinophils, macrophages, mast cells, and T cells. A recent study has revealed that IL-10+ tumor-associated macrophages may contribute to an immunoevasive microenvironment, leading to a poor prognosis and reduced response to fluorouracil-based adjuvant chemotherapy.⁵¹ The presence of a higher level of M2 phenotype of tumor-associated macrophages has been linked to poor prognosis in patients with GC.⁴¹ A previous study indicated that the presence of mast cells in the tumor area can result in a higher incidence of tumor vascularization, leading to higher rates of tumor growth and metastasis.⁵² The abundance of intratumoral PD-1+CD8+ T cells is a predictor of an unfavorable prognosis in GC.⁵³ It has been reported that pDC in tumor tissue and peripheral blood could predict poor outcome in patients with GC.⁵⁴ These findings further showed that ATP1A2 may have an impact on the prognosis of GC patients, in part through the influence of immune cell infiltration in the tumor. Collectively, our research indicated that increased expression of ATP1A2 contributed to immune cell infiltration and was associated with a negative prognosis in patients with GC, which provided insight into the correlation between ATP1A2 expression and patient outcomes in GC.

Conclusion

Overall, two subtypes of GC were characterized based on RGRA-related genes and found that different subtypes of GC exhibited different prognoses. Further insights into the relationship between RGRA and immunity were gained by examining differences in immune cell infiltration and signaling pathways between the RGRA-related subtypes of GC. Furthermore, ATP1A2 was identified as a potential prognostic biomarker that was highly correlated with immune cell infiltration in GC. This provided a theoretical basis for further studies on the role of ATP1A2 in GC carcinogenesis and progression.