RNA sequencing data, clinical information, and somatic mutation data of TCGA-STAD (n = 370) were acquired from the GDC (https://portal.gdc.cancer.gov/). Counts data were normalized to CPM values through edgeR package (Robinson et al., 2010). Microarray expression profiling and clinical information were gathered from the GSE84426 (n = 76) and GSE84433 (n = 357) (Yoon et al., 2020) on the Illumina platform, used for external verification.

Totally, 449 DNA damage repair genes were gathered from the Molecular Signatures Database (http://www.broadinstitute.org/msigdb), which were listed in Supplementary Table S1 (Liberzon et al., 2015).

Unsupervised clustering approach based upon Euclidean and Ward’s linkage was adopted for determining molecular subtypes in accordance with the transcriptional levels of DNA damage repair genes. ConsensusClusterPlus package was implemented for identifying the optimal number of clusters according to consensus cumulative distribution function (CDF) across TCGA-STAD patients (Wilkerson and Hayes, 2010). Principal component analysis (PCA) was conducted for proving the distribution difference between subtypes. Kaplan–Meier (K-M) curves were plotted for comparing overall survival (OS) of distinct subtypes, followed by log-rank test. Then, conventional clinicopathological parameters were compared between subtypes. The subtyping classification was externally verified in the GSE84426 dataset.

The up- or downregulated Gene Ontology (GO) or Kyoto Encyclopedia of Genes and Genomes (KEGG) terms in two groups were analyzed through adopting GSEA software (Subramanian et al., 2005). The reference gene sets were acquired from the GO (Ashburner et al., 2000) and KEGG (Kanehisa and Goto, 2000) databases.

Somatic mutation data from TCGA-STAD dataset were acquired utilizing TCGAbiolinks package (Colaprico et al., 2016), and mutation annotation format was analyzed and summarized with maftools package (Mayakonda et al., 2018). TMB was also computed, which was defined as the somatic mutation number per megabase of interrogated genomic sequence (Sha et al., 2020).

The fraction scores of 22 immune cell types were inferred utilizing CIBERSORT computational method (Newman et al., 2015).

Genes with differential expression were explored between the DNA damage repair-based subtypes utilizing edgeR package. The criteria were set as |fold change (FC)|≥2 and adjusted p < 0.05.

For revealing the prognostic implication of DNA damage repair, prognosis-related DNA damage repair-relevant genes with p < 0.01 were selected for least absolute shrinkage and selection operator (LASSO) analysis utilizing glmnet package (Engebretsen and Bohlin, 2019). Genes with non-zero coefficients were chosen with ten-fold cross-validation. TCGA-STAD samples were randomly classified as training and test datasets with a ratio of 1:1. Meanwhile, GSE84433 dataset was adopted as external verification. Thereafter, a DNA damage repair-relevant gene signature was generated, following the formula: R i s k S c o r e = ∑ i C o e f f i c i e n t   o f   g e n e   i * e x p r e s s i o n   o f   g e n e   i . All cases were classified as low- and high-risk groups through the median RiskScore. K-M curves of OS were conducted, and 1-, 3-, and 5-year receiver operating characteristic (ROC) curves were plotted with timeROC package. Uni- and multivariate Cox regression approaches were utilized for examining whether the RiskScore acted as an independent prognostic parameter.

Immunotherapy response was inferred through Tumor Immune Dysfunction and Exclusion (TIDE) algorithm based upon two major tumor immune escape mechanisms: triggering T cell dysfunction in tumor tissue with highly infiltrated cytotoxic T lymphocytes (CTLs) and preventing the infiltration of T cells into tumor tissue with lowly infiltrated CTLs (Jiang et al., 2018). High TIDE score indicates a greater possibility of anti-tumor immune evasion, thus exhibits a low immunotherapy response. Transcriptome data and survival information of patients who received immunotherapy were gathered from the GSE78220 (Hugo et al., 2016), GSE91061 (Riaz et al., 2017) and IMvigor210 (Necchi et al., 2017) cohorts.

Expression matrix and drug processing information were obtained from the Cancer Genome Project 2014, and IC50 values of anti-cancer drugs were estimated with pRRophetic package (Geeleher et al., 2014).

Human gastric mucosa epithelial cell line GES-1 and human gastric cancer cell lines HGC-27, MKN-28 and AGS were provided by Procell (China). All cells were cultivated in RPMI-1640 medium (Gibco, United States) plus 10% fetal bovine serum, 100 IU/mL penicillin and 100 mg/mL streptomycin at 37°C with 5% CO₂.

Total RNA was extracted utilizing Trizol (Solarbio, China), with 5 µg for cDNA synthesis. RT-qPCR was executed through SYBR Premix Ex Taq kit (Takara, China) together with Step-one real-time PCR system (ABI, United States). The mRNA level was standardized to Tubulin. The primers were synthesized by Sangon (China) (Table 1). Relative mRNA level was computed based upon 2^−ΔΔCT approach.

RIPA lysis (Solarbio) was adopted for extracting total protein, and protein samples electrophoresed onto polyacrylamide gels were loaded onto PVDF membranes that were then blocked utilizing 5% BSA for 1 hour. The membranes were probed with primary antibody of PAPPA2 (PA5-21046; ThermoFisher, United States), MPO (1/1000; ab208670; Abcam, United States), MAGEA11 (1/500; ab96236), DEPP1 (1/1000; ab230977), CPZ (15944-1-AP; Proteintech, China), COLEC12 (1/1000; ab278081), CYTL1 (1/500; ab129767) or Tubulin (1/5000; ab7291) at 4°C in a shaker overnight. Incubation with HRP goat anti-rabbit IgG (1/10,000; ab288151) for 1 hour, followed by development. ImageJ software was adopted for grayscale analysis.

Small interfering RNAs (siRNAs) targeting CYTL1 (RiboBio) were transfected into cells based upon Lipofectamine RNAiMax (Life Technologies) in accordance with the manufacturer’s protocols.

In accordance with the manufacturer’s protocols (RiboBio), intracellular ROS was measured utilizing ROS fluorescent probe (Dihydroethidium). Under a fluorescence microscope (Olympus, Japan), images were investigated and photographed.

Cells were planted into 6-well plates. Thin scratches were created utilizing a sterile pipette tip. Under an inverted microscope (Olympus), images were photographed immediately (0 h) and marked the 6-well plate so that the same field could be located again. After incubation for 24 h, the culture medium was removed and the cells were washed for removing surrounding cellular debris.

All statistical analyses were executed through R language (version 4.2.1). Two groups with normally distributed variables were compared with Student’s t-test, with Wilcoxon test for non-normally distributed variables. One-way ANOVA was conducted for comparing ≥3 groups. Pearson or Spearman correlation test was utilized for estimating the association between variables. p < 0.05 was regarded as statistical significance, as denoted: *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.

The current study gathered 449 DNA damage repair genes, and according to their expression values, TCGA-STAD samples were classified as DNA damage repair activating and inhibitory subtypes called cluster1 and cluster2 through adopting unsupervised clustering method (Figures 1A–C). In contrast to cluster2, most DNA damage repair genes displayed high expression in cluster1, demonstrating the DNA damage repair activating and inhibitory status of cluster1 and cluster2 (Figure 1D). PCA proved the prominent difference in two subtypes based upon the transcription levels of DNA damage repair genes (Figure 1E). In Figure 1F, we investigated the difference in OS outcomes between subtypes, with better OS for cluster1. This DNA damage repair-based subtyping was externally verified in the GSE84426 cohort. As expected, cluster1 exhibited the notable advantage in OS compared with cluster2 (Figure 1G). Afterwards, we assessed clinicopathological traits of two subtypes in the TCGA-STAD cohort. TCGA project has classified gastric cancer into four major subtypes. Among them, microsatellite instability (MSI) subtype often comprises genetic and/or epigenetic silence of mismatch repair genes. As shown in Figure 1H, cluster1 had the higher ratios of high MSI (MSI-H), and low MSI (MSI-L) as well as the lower ratios of microsatellite stable (MSS) in comparison to cluster2. Patients in cluster1 had older age in contrast to those in cluster2 (Figure 1I). In addition, it was found that higher ratios of male cases and more advanced histological grade were observed in cluster2 (Figures 1J, K). However, we did not observe the notable differences in pathological stage, and TNM stage between subtypes (Figures 1L–O).

Next, mechanisms underlying the two DNA damage repair-based subtypes were analyzed through GSEA. Consequently, Cluster1 exhibited higher activity of immune-related pathways (allograft rejection, asthma, and intestinal immune network for IgA production), with higher activity of linoleic acid metabolism, steroid biosynthesis and Fanconi anemia pathway in Cluster2 (Figure 2A). Thereafter, relative enrichment levels of DNA damage repair pathways were compared between subtypes. Intriguingly, homologous recombination and Fanconi anemia pathway were mainly enriched in Cluster2 (Figure 2B).

We gather somatic mutation data from the TCGA-STAD dataset, and computed TMB score. Figure 2C depicts the distribution of TMB across TCGA-STAD samples, with the median value of 2.16/MB. In contrast to cluster2, higher TMB score was observed in cluster1 (Figure 2D). SYNE2, RYR1, TTN, NEB, DOCK3, AHNAK2, EPHA5, DNAH8, FAT3, and SSPO occurred more frequent mutations in cluster1 (Figures 2E–G). In addition, we observed the co-occurrence of mutated genes in each subtype. Intriguingly, there was more extensive co-occurrence of mutation in cluster1 in contrast to cluster2 (Figures 2H, I).

Then, we measured the transcriptional levels of immunomodulators (chemokines, receptors, MHC I, MHC II, immunostimulators, and immunoinhibitors) across TCGA-STAD. Notably, chemokines (CCL14, CCL2, CCL11, etc.), and immunostimulators (CD48, CD28, LTA, etc.) had higher levels in cluster2, while most MHC molecules displayed higher activity in cluster1 (Figure 3A). Through adopting CIBERSORT computational method, we computed the relative fractions of 22 immune cell types across TCGA-STAD (Figure 3B). In addition, the heterogeneity in immune cells between subtypes was assessed. In contrast to cluster1, cluster2 exhibited higher fractions of memory and naïve B cells, M2 macrophages, resting mast cells, and monocytes as well as lower fractions of M0 and M1 macrophages, activated and resting NK cells, activated memory CD4⁺ T cells, and follicular helper T cells (Figures 3C, D).

Totally, we determined 470 up- and 1256 downregulated genes in cluster1 versus cluster2 based upon the criteria of |FC|≥2 and adjusted p < 0.05 (Figures 4A, B). Among them, 43 were DNA damage repair genes, with notable upregulation in cluster1 (Figure 4C). Thereafter, the current study evaluated their prognostic implication in gastric cancer. With p ≤ 0.01, 138 DNA damage repair-relevant genes were significantly linked to OS of gastric cancer (Table 2), which were input into LASSO analysis for variable selection. TCGA-STAD samples were randomly and equally separated into training and test datasets. When the lambda value was 0.0693, and the regression coefficient was not equal to 0, seven DNA damage repair-relevant genes COLEC12, CPZ, CYTL1, DEPP1, MAGEA11, MPO, and PAPPA2 were finally selected (Figures 4D, E). Figure 4F depicts the univariate cox regression results of above genes. All of them acted as risky factors of gastric cancer OS. Based upon them, a DNA damage repair-relevant signature was generated, following the formula: RiskScore = 0.0103488452696626 * COLEC12 expression +0.0682893864499146 * CPZ expression + 0.0361252660190149 * CYTL1 expression + 0.04553168636357 * DEPP1 expression +0.0457508199288027 * MAGEA11 expression +0.00340309980118412 * MPO expression + 0.00350839302933862 * PAPPA2 expression. Figure 4G shows the distribution of RiskScore across TCGA-STAD. With the increase of RiskScore, the transcriptional level of COLEC12, CPZ, CYTL1, DEPP1, MAGEA11, MPO, and PAPPA2 gradually increased (Figure 4H). Based upon the median RiskScore, TCGA-STAD samples were classified as low- and high-risk groups. We observed less dead and recurred/progressed cases in low-than high-risk group (Figures 4I, J). In the training dataset, low-risk cases exhibited better OS outcomes (Figure 4K). ROC curves proved the significant superiority of this DNA damage repair-relevant RiskScore in predicting long-term OS outcomes with AUC at 5-year survival >0.8 (Figure 4L).

The predictive performance of the DNA damage repair-relevant signature was further verified. Both in the test and entire datasets, low-risk cases displayed more favorable OS outcomes in contrast to high-risk cases, with excellent performance in prediction of long-term survival (Figures 5A–D). The GSE84433 cohort was adopted for external verification, and the prognostic significance of this RiskScore were proven (Figure 5E). Higher RiskScore was observed in more advanced histological grade (Figure 5F), pathological stage (Figure 5G), and T, N, M stage (Figures 5H–J) across TCGA-STAD. The similar findings were found in the GSE84433 cohort (Figures 5K, L).

To assess the sensitivity of the DNA damage repair-relevant signature in predicting patient survival, we conducted subgroup survival analysis. Gastric cancer patients were stratified into distinct subgroups according to conventional clinical parameters: sex, grade, stage, T, N, and M (Figures 6A–L). As a result, in each subgroup, high-risk patients presented poorer survival outcomes versus low-risk patients, demonstrating the excellent sensitivity of the signature in survival prediction.

From uni- and multivariate cox regression results, the DNA damage repair-relevant signature together with age were independent risky factors of OS outcomes (Figures 7A, B). Molecular mechanisms underlying the signature were assessed through GSEA. For biological process, high RiskScore was positively linked to negative regulation of cartilage development and chondrocyte differentiation, and pulmonary valve morphogenesis, with negative relationships to centriole replication, vesicle cargo loading and NLS-bearing protein import into nucleus (Figure 7C). For cellular component, higher activity of clathrin-coated endocytic vesicle membrane, parallel fiber to purkinje cell synapse, anchored component of external side of plasma membrane was observed in high-risk cases, with lower activity of endoplasmic reticulum exit site, cohesion complex and cornified envelope (Figure 7D). For molecular function, high RiskScore exhibited positive relationships with metalloendopeptidase inhibitor activity, fibronectin binding, extracellular matrix structural constituent conferring compression resistance, with negative correlations to taste receptor activity, bicarbonate transmembrane transporter activity and bitter taste receptor activity (Figure 7E). In addition, we observed higher enrichment levels of immune-related pathways (glycosaminoglycan biosynthesis-chondroitin sulfate dermatan sulfate, allograft rejection and asthma) in high-risk cases, with lower enrichment levels of nitrogen metabolism, homologous recombination and Fanconi anemia pathway (Figure 7F). Also, the RiskScore was negatively correlated to DNA damage repair pathways (nucleotide excision repair, mismatch repair, homologous recombination and Fanconi anemia pathway) (Figure 7G).

Next, relationships of the DNA damage repair-relevant signature with tumor microenvironment were assessed across TCGA-STAD. As illustrated in Figure 8A, High-risk cases exhibited higher transcriptional levels of immunostimulators (CD276, ENTPD1, TNFSF18, TNFSF4, LTA, etc.), and chemokines (CCL14, CCL17, CCL22, etc.) in contrast to low-risk cases. Higher infiltration of naïve B cells, M2 macrophages, resting mast cells, monocytes, as well as lower infiltration of M0 macrophages, resting NK cells, activated memory CD4 T cells were observed in the high-risk samples (Figures 8B, C). Above findings unveiled the heterogeneity in tumor microenvironment between low- and high-risk groups.

TIDE computational approach was adopted for inferring immunotherapeutic response. As a result, lower TIDE score was investigated in high-risk patients (Figure 8D), and this subpopulation had higher proportions of responders to immunotherapy (Figure 8E). Three independent immunotherapy cohorts: GSE91061, GSE78220, and IMvigor210 were collected for further evaluating the efficacy of the DNA damage repair-relevant signature in predicting immunotherapeutic response. It was shown that high-risk patients had better survival outcomes in comparison to low-risk patients after immunotherapy (Figures 8F–H), demonstrating that high-risk patients had the higher possibility to benefit from immunotherapy. In addition, low-risk group displayed lower IC50 scores of AZ628, Bortezomib, CHIR.99021, Cyclopamine, and GW843682X, as well as higher IC50 scores of AZD8055, Bosutinib, RO.3306, Sunitinib, and VX.702 compared with high-risk group (Figures 8I, J).

The DNA damage repair-relevant genes were further experimentally verified. In contrast to gastric mucosa epithelial cell line GES-1, transcriptional levels of PAPPA2, MPO, MAGEA11, DEPP1, CPZ, and COLEC12 were higher in gastric cancer cell lines HGC-27, MKN-28 and AGS, with lower transcriptional level of CYTL1 (Figures 9A–G). The consistent results were observed at the protein level (Figures 9H–O).

Among the DNA damage repair-relevant genes, the biological significance of CYTL1 in gastric cancer remains indistinct. CYTL1 expression was notably silenced by its specific siRNAs both in HGC-27, and MKN-28 gastric cancer cells (Figures 10A, B). According to ROS fluorescence probe detection, we found that CYTL1-silenced HGC-27, and MKN-28 cells presented the higher accumulation of ROS (Figures 10C–F). In addition, migrative capacity of HGC-27, and MKN-28 cells was remarkably impaired by CYTL1 knockdown (Figures 10G–J).