We downloaded various data of STAD patients from the TCGA data portal (https://tcgadata.nci.nih.gov/tcga/) (Hutter and Zenklusen, 2018), including mutation data, gene expression profiles, and clinical data. The STAD probe matrix file (GSE84433 series matrix) and platform file (GPL6947-13512) were also obtained from the GEO database (Barrett et al., 2007). The patient selection criteria for this study was samples of initial gastric adenocarcinoma after primary surgical resection (patients with missing clinical information were not included). To summarize the mutation data, The mutation data were summarized using the Maftools software package (Mayakonda et al., 2018). The relevant parameters of tumor-specific mutant genes were computed to obtain TMB. STAD samples were separated into high and low TMB groups using the median of TMB as the critical value. Then we looked at the relationship between TMB and survival and other clinical factors. Gene expression data from high and low TMB groups were analyzed using limma package (Ritchie et al., 2015). The fold change >2 and the false discovery rate (FDR) < 0.05 were our screening criterion for DEGs.

We used GSEA 4.0.3 software to perform Gene Set Enrichment Analysis (GSEA) and reported the first five Gene Ontology (GO) keywords and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways which enrich most significantly in the high TMB group (compared to the low TMB group) (Subramanian et al., 2005).

The fraction of invading immune cells was computed by applying the CIBERSORT algorithm, and the corresponding violin plot was created. Using the CIBERSORT algorithm with 1,000 permutations and the LM22 signature, an accurate analysis of the immune cells in STAD samples was carried out (Chen et al., 2018).

We performed ssGSEA analysis in TCGA-STAD samples using the GSVA package to determine the immune activity of 23 immune-related genesets, and the correction results ranged between 0 and 1. The heat map (pheatmap R package) and violin plots (ggpubr R package) of the tumor microenvironment were created based on the correlation analysis using the ESTIMATE algorithm between high and low risk group. To clarify the relationship between the risk grouping and immune infiltration, we also used the ssGSEA method to calculate the tumor immune cell infiltration score for each STAD sample in the TCGA database, then used the “Bioconductor Limma” R package to run the differential analysis of immune scoring and immune typing, and finally drew the box plot.

We obtained differentially expressed immune genes by intersecting the list of immune-related genes from the Immunology Database and Analysis Portal (ImmPort Database) (Bhattacharya et al., 2014) with the preceding DEGs. We randomly separated TCGA samples into two groups: a training group and an internal validation group after correlating the expression of differentially immune genes with survival time. Univariate COX analysis (p < 0.01) identified survival-related immune genes in the training group, and using the R package glmnet, the least absolute shrinkage and selection operator (LASSO) (Tibshirani, 1997) excluded the high correlation genes to avoid the model from overfitting. Finally, the best prognostic model was created using stepwise multivariate Cox regression analysis.

The developed prognostic model formula was used to determine the patient riskScores for the training, internal validation, total, and external validation (GSE84433) groups. We divided the patients in the training, internal validation, total, and external validation groups into a high-risk group and a low-risk group using the median riskScore of the training group as the threshold, plotting the survival curve and Receiver Operating Characteristic (ROC) curve of the training, internal validation, total, and external validation groups with R x64 3.6.3 software. Ultimately, univariate and multivariate prognosis analyses were done on the entire group (p < 0.05) to see if the model’s riskScore may be an independent prognostic factor.

The differences in model gene expression were evaluated between the high and low risk groups, with the ggpubr package drawing model gene boxplots to compute the expression differences. To conduct multiple gene survival analysis for all model genes, we used the PROGgeneV2 online program (Goswami and Nakshatri, 2014) and picked the GSE62254 dataset (Cristescu et al., 2015). We created a dynamic nomograph App online based on the aforesaid prognostic model to aid in the rapid calculation of patient survival rates in clinical practice. In order to demonstrate the application method, we randomly selected a low-risk group sample and a high-risk group sample from the training group, and calculated survival probability by inputting gene expression in the dynamic nomograph App.

We entered the cBioportal website (Cerami et al., 2012) and selected the study (Stomach Adenocarcinoma TCGA PanCancer data) to download the mutation status of model genes. Finally, we investigated the association between CNV of the model genes and immune cells infiltration level, as well as the correlation between immune cells infiltration level and survival of STAD patients by using the Tumor Immune Estimation Resource (TIMER) database (Li et al., 2017).

To explore whether the above prognostic model can evaluate the efficacy of immunotherapy for patients, we analyzed a series of immunotherapy biomarkers. We examined the differences between the high and low risk groups by uploading the TCGA sample expression profile into the Tumor immune dysfunction and exclusion (TIDE) database (Jiang et al., 2018) and obtaining TIDE, Microsatellite instability (MSI), Dysfunction, and Exclusion scores for each sample.

Furthermore, 5-years ROC curves were used to evaluate the riskScore’s prognostic prediction efficiency with TIDE and tumor inflammation signature (TIS) scores (Danaher et al., 2018), demonstrating the correctness of the prognostic model we developed.

We downloaded transcriptome data from the CellMiner database (https://discover.nci.nih.gov/cellminer/) and FDA-certified drug sensitivity–related data to clarify the influence of model genes in the prognostic model on drug sensitivity and tolerance. To investigate the relationship between gene expression and drug sensitivity, the Pearson correlation test was used. Then, we used the “pRRophetic” R package to draw a box plot and analyze the differences in drug sensitivity between high and low risk groups.

Human gastric cancer cell lines (BGC823,MKN45) were cultured in RPMI 1640 medium (HyClone, United States)supplemented with 10% fetal bovine serum (FBS, Gibco, United States). These cells were maintained at 37°C under an atmosphere of 5% CO2. GLP2R siRNA were synthesized by RiboBio (Guangzhou, China) and transfected into cells using Lipofectamine 3,000 (Invitrogen, California, United States) according to the manufacturer’s instructions. The cells were cultured in basic RPMI 1640 medium for 12–16 h before the media was for RPMI 1640 medium containing FBS.

RNA was extracted from gastric cancer cell lines (BGC823,MKN45) interfered with si-GLP2R and NC as control. cDNA was synthesized for real-time PCR adopting SYBR Green qPCR mix (Vazyme, China). The primers are listed as following: GLP2R -Forward: TCC​TGG​AAA​TGT​CTC​TGT​ACC​C; GLP2R—Reverse: GGC​GTT​CTC​TAT​CGT​CTG​CC; GAPDH-Forward:GACCACAGTCCATGCCATCA; GAPDH-reverse:GTCAAAGGTGGAGGAGTGGG.

RIPA lysis buffer (Servicebio, China) containing PMSF (Servicebio, China) was used to collect proteins from BGC823 and MKN45 cells. 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) was applied in separating protein samples and a polyvinylidene fluoride membrane (PVDF) membrane (Invitrogen, Carlsbad, United States) was used to transfer the separated protein. The membrane was blocked in 5% skim milk at room temperature for 1 h in a shaker, and then incubated with the primary antibody:GLP2R (Abconlal, CAT# A6602, 1:1,000), and GAPDH (Proteintech, CAT# 60004-1-Ig, 1:20,000) at 4°C overnight and subsequent incubation with the secondary antibody for 1 h. Finally, the indicated proteins were visualized by West Pico plus Chemiluminescent Substrate (Thermo Fisher Scientific, United States).

After transfection with GLP2R siRNA for 48 h, BGC823 and MKN45 cells were cultured in 96-well plates (3,000 cells/well in 100 µl RPMI 1640 medium). The proliferative capacity of treated cells was detected at 24, 48, 72, 96, and 120 h. 10% Cell Counting Kit-8 (CCK8) reagent (Bio-sharp, Hefei, China) was added to each plate according to the kit instructions, and the OD450 value was analyzed by a microplate reader (BioTek, United States). Regarding colony formation experiment, 1,000 cells were seeded in cell culture plates and allowed to grow until visible colonies formed. Then we used methanol to fix clones 15 min, 1% crystal violet to stain clones 20 min, and counted the number of clones (>50 cells).

BGC823 and MKN45 cells were transfected with GLP2R siRNA for 48 h and cultured in 24-well culture plates with 8 mm pore-containing membrane inserts to measure cell migration capacity. 4 × 10⁴ cells were seeded on the upper transwell chambers in 200 µl serum-free culture medium, and 500 µl medium containing 20% FBS was added to the lower chambers. After 24 h incubation, the cells that migrated through membranes were fixed with methanol, stained with 1% crystal violet and counted under light microscope (200×). Additionally, BGC823 and MKN45 cells were cultured in 24-well-plates and scraped with a 10-μl pipette tip. The cells were cultured in RPMI 1640 medium without FBS. Images of wounds were captured at 0, 24, and 48 h, the area of wounds was quantified by ImageJ software (40×).

All data of this study were statistically analyzed by R software 3.6.1 and Prism 9.0. The data were analyzed by two-tailed Student’s t-test and one-way analysis of variance (ANOVA). The difference was considered statistically significant when the p-value was less than 0.05.

We obtained the clinical information and whole-exome sequencing results of STAD patients from the TCGA database (Table 1). Maftools were applied to summarize the mutation data. We further categorized mutations based on the variant effect predictor, among which the frequency of missense mutations is highest (Figure 1A). Among all mutation types, the occurrence frequency of SNP is the highest (Figure 1B). Similarly, the most frequent type of SNV in STAD is C > T transversion (Figure 1C). Figure D shows the mutation of each sample, which is most relevant to TMB (Figure 1D). Figure E also depicts the mutations of the samples (Figure 1E). The first 10 mutated genes are TTN, MUC16, TP53, LRP1B, ARID1A, SYNE1, FAT4, CSMD3, FLG, and PCLO (Figure 1F). The waterfall diagram graphically shows the gene mutations of the samples (Figure 2). The correlation graph visualizes the correlation between two gene mutations; for example, mutations of PIK3CA and ARID1A co-occur, while mutations of PIK3CA and TP53 are mutually exclusive (Figure 3).

TMB was calculated by dividing the non-synonymous protein-coding variance by the genome’s overall sequence length. Based on the median TMB, we classified STAD patients into two groups: high TMB and low TMB (Supplementary File S1). The survival rate of patients in the high TMB group is higher than that of patients in the low TMB group, according to Kaplan-Meier survival analysis (Supplementary Figure S1A). The results show that in STAD, TMB is favorably connected with age and negatively correlated with stage, T, and N stages; female patients have higher TMB than male patients. TMB in female patients is higher than in male individuals (Supplementary Figures S1B–H). Between the high and low TMB groups, 816 DEGs were found (Supplementary File S2).

The top 5 GO keywords (Figures 4A–E) and KEGG pathways (Figures 4F–J) with the highest significant enrichment in the high TMB group (p < 0.05) were obtained using GSEA enrichment analysis comparing the high and low TMB groups. These KEGG pathways and GO keywords are mostly related to genetic material metabolism and DNA repair.

According to studies, the higher the TMB in tumors, the more neoantigens are produced, making tumors more immunogenic (Rizvi et al., 2015). As a result, we investigated the link between TMB and immune markers in STAD. Using the CIBERSORT method, we determined the fraction of invading immune cells. The inference of immune cell types generated by CIBERSORT was relatively reliable at the p < 0.05 threshold. These findings demonstrated that high TMB tumors included large proportions of follicular helper T cells, activated memory CD4+ T-cells, M1 macrophages, M0 macrophages, and Neutrophils. Resting memory CD4+ T cells, regulatory T cells, monocytes, resting dendritic cells, and resting mast cells were all higher in the low TMB group (Figure 5).

From the overlap of immune-related genes and prior DEGs, differentially expressed immune genes were derived (Supplementary File S3). One-time random grouping was used to create the training (Supplementary File S4) and internal validation groups (Supplementary File S5). Based on univariate COX analysis (p < 0.01) for the training group (Supplementary File S6), we applied the LASSO algorithm to deal with the 21 survival-associated immune genes. Following the stepwise multivariate Cox regression analysis, four immune genes were able to develop a predictive model. The riskScore for each patient is calculated as follows: 0.001*APOD +0.005*APOH +0.039*INHA +0.499*GLP2R (Table 2). Each model gene is considered to with high risk attribute.

We split the patients in the training, internal validation, total, and external validation groups into high and low risk groups (Supplementary Files S7, S8, S9, S10) for later survival analysis after computing their riskScores. The survival curves between the high and low risk groups in the training, internal validation, total, and external validation groups are notably different, and the low risk group’s survival rate is significantly greater than the high risk group’s. The accuracy of the prognostic model we developed is dependable, according to the area under the curve (AUC) values of ROC curves (Figure 6). Ultimately, univariate and multivariate prognostic studies (p < 0.05) show that the riskScore derived from the model is a prognostic factor that is independent of other factors (Supplementary Figure S2).

Boxplots of model genes were created to compute the expression level differences between the high and low risk groups (Figures 7A–D). The high-risk group has significantly higher expression of all model genes, which validates the earlier inference that all model genes belong to high risk factors. Furthermore, incorporating all model genes in a multiple gene survival study successfully verifies the efficiency of the prognostic model (Figure 7E). A dynamic nomograph App online (https://u20131050.shinyapps.io/STAD-TMB-Dynamic_nomogram/) was successfully designed to quickly compute patient survival and support clinical decision making in order to improve the clinical translational implications of our study. We take low-risk group sample TCGA-VQ-AA6J and high-risk group sample TCGA-D7-A6F0 as examples to show the calculation results of dynamic nomograph App, which shows that the prediction results are relatively accurate (Supplementary Figure S3).

The cBioportal website provided an overall view of genetic change (Figure 8A) as well as domain mutation plots (Figures 8B–E) for the model genes. These model genes have extremely low mutation frequencies.

According to the TIMER database, CNV of these model genes has little effect on immune cell infiltration level differences (Figures 9A–D), and only Macrophage cell infiltration has a substantial impact on STAD patients’ survival (Figure 9E).

We created a heatmap (Figure 10A) and violin plots (Figures 10B–E) of the tumor microenvironment by analyzing the correlation between risk categorization and tumor microenvironment. In conclusion, high-risk group had higher StromalScore, ImmuneScore, and ESTIMATEScore than low-risk group, while the low-risk group had a higher score of TumorPurity. To investigate the relationship between the risk score and immune status, we used the ssGSEA method to quantify 23 immune cell subsets, discovered that the infiltration of 13 immune cell subsets was significantly higher in the high-risk group than in the low-risk group (p < 0.05). (Figure 10F).

The high risk group had a higher TIDE score, lower MSI score, higher Dysfunction score, and higher Exclusion score, indicating that the high risk group of STAD has increased immune escape potential, resulting in poor immunotherapy efficacy (Figure 11A).

In the end, it can be judged from the 5-years ROC curves that the prognostic model constructed by us has the largest AUC value, so its prognostic prediction efficiency is higher than TIDE score and TIS score (Figure 11B). The above results show that the group with high expression of GLP2R has the highest risk score and the worst prognosis.

By performing a separate drug sensitivity analysis on model genes in the prognostic model, we were able to identify the top 16 drugs with the most statistically significant differences. APOD expression was found to be positively related to the sensitivity of vemurafenib, pd-98059, dabrafenib, hypothemycin, selumetinib, bafetinib, denileukin diftitox ontak and cobimetinib (isomer 1), it indicates that the higher level of APOD expression, the greater sensitivity to the aforementioned drugs, but APOD has a negative correlation with pyrazoloacridine, batracylin, docetaxel and pralatrexate. APOH expression was revealed to be highly linked to the sensitivity of elesclomol, INHA expression was discovered to be strongly tied to the sensitivity of fulvestrant, but inversely correlated with the sensitivity of amonafide. Furthermore, the higher the expression of GLP2R in STAD patients, the greater the patient’s sensitivity to decitabine (Figure 12A). To further improve the clinical value of tumor mutation burden-related prognostic model for the treatment of stomach cancer. We analyzed the commonly used drugs in the clinical treatment of gastric cancer, which include cisplatin, doxorubicin, gemcitabine, lapatinib, and gemcitabine was noticed to be more sensitive in the high-risk group than in the low-risk group (p < 0.001), whereas lapatinib is the opposite (Figures 12B–E).

To evaluate the specific role of GLP2R in gastric cancer, the relative mRNA expression level of GLP2R in the eight gastric cancer cell lines (GES1, BGC-823, MKN45, SNU-216, SGC-7901, MGC-803, AGS and N87) was analyzed. The expression level of GLP2R in the BGC-823, MKN45 cell lines was higher than that in the other cell lines (Figure 13A). We first evaluated the transfection efficiency of the cells by qRT-PCR and Western blot, found that the relative expression level of GLP2R was significantly lower after siRNA 3 transfection (Figures 13B,C). To further confirm the role of GLP2R in proliferation, We performed CCK-8 assays to detect the effect of GLP2R knockdown. After GLP2R silencing, BGC-823, MKN45 cell proliferation significantly decreased compared to control cells (Figure 13D). Colony formation assay also indicated that GLP2R silencing significantly suppressed the growth of BGC-823 cell (Supplementary Figures S4A–B). Transwell and wound healing assays were performed to detect migration, Our results showed that the migration rates of BGC-823, MKN45 cells transfected with siRNA were significantly lower than that of the control-transfected cells (Figures 13E–H). These data suggest that GLP2R knockdown repressed the proliferative and migratory abilities of BGC-823, MKN45 cells.