Prognosis and tumor immune microenvironment of patients with gastric cancer by a novel senescence-related signature.

BACKGROUND: Cellular senescence is a stable state of cell cycle arrest that plays a crucial role in the tumor microenvironment (TME) and cancer progression. Nevertheless, the accurate prognosis of gastric cancer (GC) is complicated to predict due to tumor heterogeneity. The work aimed to build a novel prognostic model in GC.
METHODS: LASSO and Cox regression analysis were constructed to develop a prognostic senescence-related signature. The Gene Expression Omnibus dataset was used for external validation of signature. Afterward, we performed correlation analysis for the risk score and the infiltrating abundance of immune cells, TME scores, drug response, tumor mutational burden (TMB), and immunotherapy efficacy.
RESULTS: Five senescence-related genes (AKR1B1, CTNNAL1, DUSP16, PLA2R1, and ZFP36) were screened to build a signature. The high-risk group had a shorter overall survival, cancer-specific survival, and progression-free survival when compared to the low-risk group. We further constructed a nomogram based on risk score and clinical traits, which can predict the prognosis of GC patients more accurately. Moreover, the risk score was evidently correlated with infiltration of immune cells, TME score, TMB, TIDE score, and chemotherapy sensitivity. Meanwhile, the Kyoto Encyclopedia of Genes and Genomes pathway showed that the PI3K-Akt and Wnt signaling pathway were differentially enriched in the high-risk group.
CONCLUSIONS: The senescence-related signature was an accurate tool to guide the prognosis and might promote the progress of personalized treatment.

1. Introduction

Gastric carcinoma (GC) is the fourth most common cancer, and is the third frequent cancer related mortality globally.[ At present, the main treatment methods for GC include a combination of surgical resection, chemotherapy and targeted therapy, and immunotherapy.[ Despite considerable progress in treatment strategies, the mortality rate of GC is still high.[ Currently, a significant clinical challenge for GC is that the disease is often asymptomatic at an early stage and diagnosed at a later stage, characterized by a poor prognosis and limited treatment options.[ Treatment strategies varied widely, due to the lack of a specific biomarker to assess GC progression. Thus, a reliable biomarker for better understanding of GC and can be utilized for personalized treatment is pressingly needed.

Cellular senescence is a stable state of cell cycle arrest that is penalized in response to various insults, such as deoxyribonucleic acid damage, telomere dysfunction, telomere shortening, and oncogenic stress, resulting in inhibition of potentially dysfunctional, transformed, or aged cells.[ It is characterized by distinct morphological hallmarks, alterations in gene expression and chromatin structure, expression of senescence-associated β-galactosidase, and acquisition of a senescence-associated secretory phenotype.[ Cellular senescence is a double-edged sword that plays an important role in different stages of tumorigenesis, such as tumor initiation, establishment, and escape.[ Oncogene activation can promote cell proliferation and induce senescence.[ Therefore, cellular senescence is the mechanism by which cells are permanently withdrawn from the cell cycle to prevent the progression of benign neoplastic lesions to malignant tumors. However, work in recent decades has demonstrated that senescent cells often also have oncogenic properties. It is mainly because senescent cells remain viable and biologically active over a long period and eventually resume proliferation while signaling their microenvironment with heterotypic.[ Taken together, the senescence status of tumors is heterogeneous. Hence, individualized treatment of GC is an urgent task at present.

In this study, we constructed a senescence-associated signature as a prognostic biomarker. To develop new therapeutic methods for the treatment of GC, we evaluated the impact on immune cell infiltration, clinical outcomes, TMB, immunotherapy, and chemotherapy sensitivity of the signature in GC.

2. Materials and methods

2.1. Data preparation

In this study, a total of 414 GC patients, including complete ribonucleic acid-seq Fragments Per Kilobase Million data and clinical characteristics, were obtained from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). The microarray data set (GSE84437) was obtained from the Gene Expression Omnibus database. TCGA-stomach adenocarcinomas was used for building the prognosis model, and 431 samples from GSE84437 were used to validate the risk model. The samples with missing clinical information and prognosis < 30 days were excluded. This study did not require ethics approval as all data were obtained public databases. Moreover, 279 senescence-related genes were obtained from the CellAge database (https://genomics.senescence.info/cells/signatures.php?).

2.2. Generation and validation of senescence signature

Univariate Cox regression analysis helped screen the potential senescence-related prognostic genes in the TCGA. Following that, LASSO Cox regression analysis was employed to select the most optimal candidates. Subsequently, the multivariable Cox regression analysis was employed to establish prognostic signature. The risk score was estimated using:

In this formula, coef and exp represent correlation coefficient and gene expression levels, respectively. The respective median value was set as the cutoff value and GC patients were classified into low- and high-risk groups. Kaplan–Meier survival curves and time-dependent receiver operating characteristic (ROC) were used to evaluate the prognostic predictive performance of the signature. Principal component analyses (PCA) were performed with “scatterplot3d” R package. The nomogram is drawn based on clinical data and risk scores, and its accuracy is demonstrated by a calibration curve. In addition, in the GSE84437 set, the same methods were employed to examine the accuracy of the signature.

2.3. Assessment of immune landscape in tumor microenvironment (TME)

To assess the link between the signature and tumor-infiltrating immune cells (TIICs) in GC, the CIBERSORT algorithm was employed to quantify 22 TIICs. Based on the ESTIMATE algorithm, we calculated the stromal score, immune score, ESTIMATE score for each sample to quantify TME.

2.4. Predicting chemotherapy drug and immunotherapy response

To access the significance of the signature in predicting the sensitivity to chemotherapy in GC, the “pRRophetic” package was conducted to calculate the half-maximal inhibitory concentration (IC50) of the main chemotherapeutic medications used in the treatment of GC patients. In addition, Tumour Immune Dysfunction and Exclusion (TIDE) algorithm was used to assess the response of immunotherapy.[

2.5. Analysis of somatic mutations and TMB among different risk groups

The number of gene mutations and the type of gene mutation were derived from the somatic mutation data of the TCGA database. We also visualized the top 15 genes with the highest mutation frequency in the low- and high-risk subgroups using the “Maftool” R package.

2.6. Functional analysis

Use the “limma” R package to determine the differentially expressed genes (DEGs) in different risk subgroups. The filter standard is: |LogFC| > 1, and the corrected P-value is < .05. The “clusterProfiler” package was adopted to make gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses on the DEGs.

3. Results

3.1. Establishment of senescence-related signature

In the TCGA set, univariate Cox proportional hazards regression analysis showed that 16 genes were associated with prognosis of GC samples (Fig. 1A). Then, we performed LASSO Cox analysis based on the 16 genes, following 5 senescence-related genes were screened and a risk formula was built using the coefficient from multivariable Cox regression analysis (Fig. 1B–D). Risk score = 0.010 × (AKR1B1 expression value) + 0.018 × (CTNNAL1 expression value) − 0.024 × (DUSP16 expression value) + 0.100 × (PLA2R1 expression value) + 0.001 × (ZFP36 expression value). The risk score of all samples could be calculated according to the risk formula and samples could be classified into low- and high-risk subgroups based upon whether their median risk scores were above or below the median value. In the TCGA set, the survival status and risk sore distribution showed that death risk increased gradually with the increment of risk scores (Fig. 1E). PCA was performed to compare the difference between two risk subgroups based on expression profiles of whole-genome, 279 senescence-related genes, and 5 senescence-related genes of the signature (Fig. 1F–H). The result showed good discriminative performance of the signature (Fig. 1H). Kaplan–Meier analysis demonstrated that the higher risk score was often associated with worse overall survival (OS), cancer-specific survival, and progression-free survival (Fig. 1I–K). The result of ROC curves demonstrated that this signature was accurate and stable to predict the prognosis in GC (Fig. 1L). Meanwhile, after adjusting other clinicopathological variables, the signature remained a powerful and independent factor (Fig. 2A and B).

Figure 1.

Generation of the prognostic signature in the TCGA set. (A) Univariate Cox regression analysis identified 16 senescence-related prognostic genes. (B and C) Lasso regression analysis. (D) The presentation of 5 independent prognosis genes in multivariate Cox regression analysis. (E) The distribution of risk scores and the coherence of survival time and survival status among the two risk subgroups. (F–H) PCA is based on the expression profiles of whole-genome (F), 279 senescence-related genes (G), and 5 senescence-related genes of the signature (H). (I–K) Kaplan–Meier curve of OS (I), CSS (J), and PFS (K) in different risk patients. (L) ROC curve to assess the prognostic value of risk score. CSS = cancer-specific survival, OS = overall survival, PCA = principal component analyses, PFS = progression-free survival, ROC = receiver operating characteristic, TCGA = the Cancer Genome Atlas.

Figure 2.

Establishment and validation of the nomogram. (A and B) Univariate (A) and multivariate analyses (B) indicated the prognostic value of the risk score. (C) Nomogram for predicting 3- and 5-year OS. (D) ROC curve analysis of the nomogram. (E) The calibrate curves are used to verify the nomogram the at 3 and 5 years. ROC = receiver operating characteristic, OS = overall survival.

To better represent the clinical utility of the signature, we built a nomogram combining risk scores with clinicopathological traits, such as age and TNM stage (Fig. 2C). In this nomogram, different clinical information and risk groups corresponded to specific scores, and we could predict 3 and 5-year OS for a patient according to the sum of the scores. The area under the curve values of the ROC curve at 3 and 5 year-OS for the nomogram were 0.700 and 0.729 (Fig. 2D). The calibration curve of the nomogram showed that the constructed nomogram was more consistent with the actual situation (Fig. 2E).

3.2. External validation of the signature

The risk score of samples in GSE84437 validation set was calculated using the formula from the training cohort to verify the robustness of the signature. The distribution of the risk score of each individual revealed that the high-risk had shorter OS, the similar results were received from Kaplan–Meier curves analysis (Fig. 3A and B). PCA was employed to verify the accuracy of categorization. As shown in Figure 3C, there seems significant distinction between the two subgroups. The area under the curve values for the 3-year (0.680) and 5-year (0.702) OS in the validation set showed good sensitivity and specificity (Fig. 3D).

Figure 3.

Verification of the prognostic model in the GSE84437 set. (A) The distribution of risk score and the coherence of survival time and survival status among two risk subgroups. (B) Kaplan–Meier curve of OS in low-risk and high-risk patients. (C) PCA demonstrated overt separation of both subgroups. (D) ROC curve for assessing the predictive accuracy of the signature. OS = overall survival, PCA = principal component analyses, ROC = receiver operating characteristic.

3.3. Correlation between risk score and TME

To better evaluate the correlation between risk scores and TIICs, we employed CIBERSORT algorithm to evaluate the infiltrating level of immune cells and made comprehensive comparisons of two risk groups. According to the Figure 4A and B cells, the high-risk score was negatively linked to infiltration of CD8 + T cells, plasma cells, activated memory CD4 + T cells, and natural killer cells resting, and positively linked to infiltration of B cells naïve, monocytes, macrophages M2, and T cells regulatory (Tregs). Furthermore, the stromal score, immune score, and ESTIMATE score of the low-risk group were overtly higher than those of the high-risk group (Fig. 4B).

Figure 4.

Immune landscape and chemotherapy drug sensitivity of different risk subgroups. (A) Differences in immune cells infiltration between high-risk and low-risk groups. (B) Correlation between the risk score and the TME score. (C) Comparison of TIDE score between two risk groups. (D–I) Comparison of sensitivities to commonly used chemotherapeutics for GC in two risk groups. GC = gastric cancer, TIDE = tumor immune dysfunction and exclusion, TME = tumor microenvironment.

3.4. Immunotherapy response and chemotherapy drug selection

We used the TIDE to assess the likelihood of GC benefiting from immune checkpoint inhibitors therapy. The result indicated that patients in low-risk group led to a statistically significant decrease in TIDE score (P < .05), indicating a better response to immunotherapy (Fig. 4C). In addition, we compared responses to chemotherapeutic drug in two risk groups. We observed that the IC50 of dasatinib, foretinib, lapatinib, and pazopanib was relatively low in the high-risk group, while the IC50 of gefitinib and veliparib were relatively high in the high-risk group (Fig. 4D–I).

3.5. Correlation of risk score with somatic mutation landscapes

We predicted the gene mutation frequency of the high- and low-risk groups. Compared with the low-risk group (92.44%), the mutation rate of the high-risk group (86.47%) was lower (Fig. 5A and B). The top five mutated genes in the two risk groups were TTN, TP53, MUC16, LRP1B, and ARID1A. In addition, we found that the TMB of the high-risk group were evidently lower than those of the low-risk group (Fig. 5C). The Kaplan–Meier survival curves revealed that the OS of the low-TMB group was significantly lower than that in high-TMB group (Fig. 5D). The GC patients with high TMB in low-risk group had significantly longer OS relative to high-risk patients. And the patients with low TMB showed consistent results (Fig. 5E), indicating the prognostic value of the risk score was not affected by TMB.

Figure 5.

Correlation of risk score with somatic mutation landscapes. (A, B) Waterfall chart showing higher mutation frequency genes in high- (A) and low-risk groups (B). (C) Comparison of TMB in the high- and low-risk groups. (D) Survival analysis of high-and-low-TMB group. (E) Survival analysis of GC stratified by TMB and risk score. TMB = tumor mutational burden, GC = gastric cancer.

3.6. Functional analysis

To explore potential biological functions of the signature, we performed GO and KEGG functional enrichment analysis on 651 DEGs and presented significantly enriched biological processes. The GO enrichment results showed that DEGs were significantly enriched in muscle contraction and muscle system process (biological process), collagen-containing extracellular matrix and contractile fiber (cellular component), as well as actin binding and glycosaminoglycan binding (molecular function) (Fig. 6A). The KEGG enrichment results revealed that DEGs prominently enriched in PI3K-Akt signaling pathway, focal adhesion, and Wnt signaling pathway (Fig. 6B). This suggests that the potential biological functions of DEGs are related to the tumor.

Figure 6.

GO and KEGG pathway analysis of the senescence-related signature. (A) The results of GO enrichment analysis. (B) The results of KEGG enrichment analysis. GO = gene ontology, KEGG = Kyoto Encyclopedia of Genes and Genomes.

4. Discussion

To date, Japan has approved two ant-PD-1 inhibitors for GC: nivolumab as a third- or later-line treatment for GC, and pembrolizumab for previously treated patients with microsatellite instability-high tumors.[ Although immunotherapy has achieved unprecedented success, the OS time of individuals varies greatly due to the existence of heterogeneity.[ Senescence is a complex biological process with cell-autonomous and paracrine effects that has a major impact on the TME.[ Senescent tumor cells can culture and shape the tumor microenvironment through senescence-associated secretory phenotype. In the TME, senescent tumor cells are surrounded by tumor cells, stromal cells, and infiltrating immune cells. These immune cells include cells that promote tumor growth and inhibit tumor growth.[ Senescence-associated secretory phenotype exerts anti-tumor and pro-tumor effects by recruiting and activating immune cells and neighboring cells.[ Hence, modeling GC has important implications for deciphering whether molecular determinants of senescence remodel the TME and whether this model has an impact on clinical outcomes and immunotherapy response in GC patients.

In this work, univariate Cox analysis was performed to screen senescence-related genes that were associated with OS. Subsequently, we performed LASSO and multivariate Cox regression analysis in the training set to construct a prognostic model. Based on the result of multivariable Cox regression analysis, a 5-gene (AKR1B1, CTNNAL1, DUSP16, PLA2R1, and ZFP36) signature was built. The signature consists of four risk prognostic genes (AKR1B1, CTNNAL1, PLA2R1, and ZFP36) and one protective prognostic gene (DUSP16). AKR1B1, a member of AKR1 subfamily B, plays a vital role in glucose metabolism and osmoregulation.[ The role of AKR1B1 in cancer is not fully understood, but growing evidence suggests that it has a large impact on cancer progression.[ AKR1B1 can be involved in complex signaling pathways, proteins, and miRNA networks, such as mir-21-mediated mechanisms such as cell cycle, apoptosis, epithelial-mesenchymal transition.[ Wu et al[ revealed that AKR1B1 promoted basal-like breast cancer progression by activating the epithelial-mesenchymal transition program. PLA2R1, a member of the C-type lectin superfamily, has been reported to be a prognostic or diagnostic biomarker for breast cancer.[

The GC patients were then divided into low- and high-risk subgroups. K-M curves revealed that the OS, cancer-specific survival, and progression-free survival rate was markedly higher in the low-risk subgroup. Internal and external validation confirmed the value of the predictive signature. Furthermore, the risk score can be adopted as a prognostic predictor independent of clinical characteristics. To better represent the clinical utility of the signature, we constructed a prognostic nomogram combining risk scores with clinicopathological traits. The KEGG enrichment results revealed that DEGs are prominently enriched in tumor-related signaling pathways. These results suggest that this signature may be effective in assessing patient outcomes, thereby facilitating the implementation of this signature in future clinical practice.

TME is a key mediator of tumor progression and treatment outcome. Its status can reflect the response of tumor patients to immunotherapy.[ Thus, exploring the diversity and complexity of TME can improve the predictive power and clinical guidance of immunotherapy. We assessed the correlation between risk score and TME and showed that the high-risk group presented a typical tumor immunosuppressive microenvironment. This was reflected in that compared with the low-risk group, the infiltration proportion of activated memory CD4 + T cells, CD8 + T cells, plasma cells, and resting natural killer cells in the high-risk group was lower, while the infiltration abundance of macrophages M2 and Tregs was increased. T cells are among the most common immune cells, and T cells can help prevent metastatic progression and are known to be related to patient prognosis.[ This is consistent with our results that the low-risk group had higher levels of CD8 + and CD4 + T cells. M1 and M2 macrophages are one of the major tumor-infiltrating immune cell types, the former typically exerting antitumor effects, including direct-mediated cytotoxicity and antibody-dependent cell-mediated cytotoxicity killing tumor cells; while M2-polarized macrophages, generally considered tumor-associated macrophages, promote tumor cell initiation and metastasis by promoting angiogenesis and lymphangiogenesis and suppressing T cell-mediated antitumor immune responses.[ Tregs suppress the immune activity mainly by reducing the activation of CD8 + T cells, moreover, Tregs inhibit the immune response through its membrane-bound cytotoxic T-lymphocyte-associated protein 4 immune checkpoint molecule.[ In addition, TME scores in the low-risk group were higher than than those in high-risk group.

Research has shown that immunotherapy became one of the most promising methods for GC treatment in recent years.[ In this work, we revealed that low-risk patients presented with lower TIDE scores, indicating better response to immunotherapy. TMB has been proven to be a promising indicator to forecast the response to immune checkpoint inhibitors, which is significantly correlated with immunity.[32] We found that the TMB related to immunotherapy sensitivity was obviously higher in the low-risk subgroup than in the high-risk subgroup. Studies have demonstrated that patients with higher TMB more likely to benefit from immunotherapy.[ These results suggest that patients in the low-expression group may benefit more from immunotherapy. Stratified analysis showed that the signature was not associated with TMB in GC. This means that the signature and TMB represent different aspects of tumor immunobiology and can predict the patient response to immunotherapy independent of TMB. Lastly, GC patient sensitivity to chemotherapy and targeted drugs was next assessed. High-risk patients were more sensitivity to dasatinib, foretinib, lapatinib, and pazopanib, while low-risk patients were more sensitivity to gefitinib and veliparib. Therefore, suitable chemotherapy and targeted drugs can be selected according to different risk groups.

5. Conclusions

In this research, we successfully established a 5-gene senescence-related signature that could be used to classify GC patients. This signature can improve the guidance of the patient personalized treatment.

Authors’ contributions

GLZ and WZ conceptualized and designed this study. GLZ, KCD, and JPL collected the data. GLZ, KCD, and WZ analyzed the data, drew the pictures. GLZ, KCD, and JPL drafted the manuscript, and WZ revised the final manuscript.

Conceptualization: Guanglin Zhang, Wei Zhou.

Data curation: Guanglin Zhang, Kechen Dong, Jianping Liu.

Formal analysis: Guanglin Zhang, Kechen Dong, Wei Zhou.

Investigation: Guanglin Zhang, Kechen Dong.

Methodology: Guanglin Zhang, Kechen Dong.

Project administration: Guanglin Zhang, Kechen Dong.

Resources: Guanglin Zhang, Kechen Dong, Jianping Liu.

Software: Guanglin Zhang, Kechen Dong, Jianping Liu.

Supervision: Wei Zhou.

Validation: Wei Zhou.

Visualization: Wei Zhou.

Writing – original draft: Guanglin Zhang, Kechen Dong, Jianping Liu.

Writing – review & editing: Wei Zhou.