Identification of senescence as a hallmark of gastric cancer
To elucidate the senescence character of GC, we comprehensively explored the transcriptomic data in multiple datasets of GC. In TCGA-STAD cohort, the results of gene set enrichment analysis (GSEA) showed that those biological pathways indicative of cellular senescence, including cellular senescence (NES=1.63, p<0.01), DNA damage/telomere stress induced senescence (NES=1.98, p<0.01), and oxidative stress induced senescence (NES=1.58, p<0.01), were significantly enriched in the GC over normal gastric tissues (Fig. 2A). And Gene Ontology (GO) about cellular senescence was also significantly enriched in GC compared with normal tissues (Fig. S1A). Similarly, consistent results for senescence were observed in the another independent cohort (Fig. 2B, Fig. S1B). Moreover, the gene set variation analysis between GC and their paired adjacent normal tissues in 3 datasets showed GC qualified with markedly activation in multifarious senescence pathways and GO terms (Fig. 2C). These results suggested senescence as a prominent hallmark of GC.
To investigate the genes impacting senescence in GC, we firstly screened the differently expressed genes (DEGs) between GC and no-tumor adjacent tissues in TCGA-STAD dataset. It was shown that 3225 DEGs were identified, including 1455 up-regulated and 1770 down-regulated genes (Fig. S1C, D), and the up-regulated genes were markedly enriched in the cellular senescence pathway (Fig. S1E). Then, we acquired senescence genes based on all the gene sets associated with senescence in Molecular Signatures Database 19 (Table S1). Combined the DEGs with senescence genes, 103 senescence-associated genes in GC were identified. GO and pathway enrichment analysis of these genes was conducted, showing significantly enriched biological processes comprising aging, cell cycle and cellular senescence (Fig. 2D, E). As shown in Fig. 2F, the interaction network between top five enriched pathways and genes encompassed the specific marker for cellular senescence (CDKN2A) and immune-associated genes (e.g. CXCL8 and IL6).
Single cell profiling of senescence within tumour microenvironment in gastric cancer
To further dissect the senescence characteristics within GC tissues at single cell level, we investigated the single cell transcriptome of 19042 cells from normal and cancerous gastric tissues (Fig. 3A, B, Fig. S2A). In agreement with the results from the bulk analysis above, more cells with significant activation of cellular senescence pathway were identified in GC as compared with those in normal gastric tissues (Fig. 3C, Fig. S2B). Interestingly, in GC, certain cell types consisting of endothelial cells, enteroendocrine cells, proliferative cells and macrophages as well as fibroblasts, but not the cancer cells, are highly enriched in various senescence pathways (Fig. 3D), suggesting these cell types, rather than the cancer cells, may primarily contribute to the senescence features in GC.
To identify the key genes associated with senescence in GC, we integrated significantly different expression genes (DEGs), the survival-associated genes and senescence genes derived from senescence gene sets, and finally established the core senescence genes (CSGs) of GC (Fig. 3E, F). At single cell level, CSGs were expressed with higher proportions in cells of TME such as endothelial cells, fibroblasts and macrophages (Fig. 3G), which is in line with above pathway analysis.
Furthermore, we assessed the incidence of copy number variations (CNV) and somatic mutations of 17 CSGs in GC using genomic data. It was indicated that COL1A2 gene exhibited the highest mutation frequency followed by ABCA8, while the investigation of CNV alteration frequency showed that a prevalent CNV amplification occurred in COL1A2, C7, ABCA8, as well as SERPINE1 gene (Fig. 3H and Fig. S2C, D). Although the gene-expression changes were associated with the genomic transformation, especially the alteration of CNV (Fig. S2C,E), the complexity of genomics and heterogeneity of tumor led to the incomplete consistency in the landscape of genomic and expressional alterations, highlighting that expression characteristics of CSGs played a considerably pivotal role in GC.
The identification of senescence subtypes based on the core senescence genes
In order to decipher the underlying biological characteristics of the senescence phenotype, the unsupervised hierarchical clustering analysis was performed firstly in the discovery cohort of TCGA-STAD based on the expression of CSGs. Two senescence-associated clusters were identified (Fig. 4A), and principal component analysis (PCA) indicated the significant distinction existed on the transcriptional profile between the two subtypes (Fig. 4B, Fig. S3A). Except for EZH2, TNFAIP2 and IL1A, the remaining genes exhibited higher expression in cluster 1 than those in cluster 2 (Fig. S3B). Furthermore, functional enrichment analysis showed that cluster 2 displays significantly higher activation in senescence pathways, suggesting cluster 2 was associated with senescence subtype (Fig. 4C). Survival analysis showed that cluster 2 was associated with longer overall survival (OS) (Fig. 4D), which was independently validated in a pooled cohort of 1484 patients (Fig. S3C,D). Interestingly, the higher mutation density was observed in cluster 2 in comparison with cluster 1 in TCGA data set (Fig. 4E).
Moreover, the senescence subtypes were parsed in the clinical and molecular features in TCGA cohort. The results showed that majority of patients with early clinical stage (stage I) were in cluster 2 (Fig. S4A). And no significant distribution difference was found between groups stratified by the age of patients (Fig. S4B), which was probably resulted from the little age gap among patients (Table S2). We also found that tumors in cluster 1 presented poorer differentiation and were enriched in the diffuse and signet-ring histological classification (Fig. S4C), suggesting cluster 1 was associated with the more aggressive subtype. With regards to the molecular subtypes of TCGA project 3, the genomically stable (GS) subtype was mainly concentrated in the cluster 1, while the subtypes of microsatellite instability (MSI) and Epstein-Barr virus (EBV)-positive substantially appeared in the cluster 2 (Fig. S4D). In consistent with the molecular subtype, cluster 2 was higher in the proportion of MSI tumors by the MSI status analysis and the clonal deletion score for quantificational chromosomal instability (CIN) than those in cluster 1 (Fig. S4E,F).
Additionally, we analyzed the canonical biological processes associated with TME, indicating that the cluster 2 displayed activation of DNA damage response (DDR), antigen processing and presentation (APM), which may result from the higher mutation load (Fig. 4F)20. Furthermore, we then explored the characteristics of the immune cell infiltration in 1484 GC patients in the validation cohort 1. The results showed that cluster 1 showed significant increases in the infiltration of anti-tumor immune cells such as the activated CD8+ T cells and central memory CD4+ T cells, as well as pro-tumor immune cells like the regulatory T cells (Treg), and myeloid-derived suppressor cells (MDSC). By contrast, cluster 2 appeared high infiltration of activated CD4+ T cells, CD56 bright natural killer cells (NK), and memory B cells (Fig. 4G). These results indicated that the senescence subtypes were associated with distinct features of TME.
Construction of a senescence scoring model
The senescence phenotype acted non-negligible roles in shaping the TME and markedly affected the clinical prognosis. However, these analyses were only based on the patient population and could not accurately predict the senescence character in individual patients. Therefore, we sought to construct a scoring system to quantify the senescence level in tumour tissues for individual patient with GC. In the TCGA-STAD cohort, we constructed a senescence scoring model, termed “senescore”, by using 6 senescence signature genes (ADH1B, IL1A, SERPINE1, SPARC, EZH2 and TNFAIP2) identified with the LASSO Cox regression algorithm (Fig. S5A, S5B). The survival analysis demonstrated patients with high senescore were correlated with better outcomes in TCGA-STAD cohort (log-rank test, P<0.0001; median OS 19.6 months versus 56.2 months, Fig. 5A and Table S3). Furthermore, when the senescore was evaluated as a continuous variable with the multivariate Cox regression model, the senescore was determined to be an independent and robust prognostic factor (HR, 0.14; 95%CI, 0.065–0.32; P<0.001; Fig. 5B). More importantly, when patients were stratified with TNM clinical stages, the senescore successfully separated patients with distinct clinical prognosis at the same disease stage (Fig. 5C). These findings suggested that the senescore reserved its prognostic relevance even after classic clinicopathologic prognostic features have been taken into accounts. Therefore, for better clinical utility, we integrated senescore with the age and clinical stages, two independent prognostic factors identified above, to generate a comprehensive nomogram (Fig. 5D). The prognostic capacity of the nomogram was also demonstrated by the area under the curve (AUC) of the time-dependent receiver operating characteristic (ROC) curve with a value of 0.794 at 5 years (Fig. S5C). Furthermore, the calibration curves of the model for the possibility of OS at 3 years and 5 years closely approximated the observed estimates, manifesting the accurate predictive ability (Fig. S5D,S5E).
To better illustrate the feature of senescore, the overview of the connection among the senescence subtypes, molecular subtypes, immune subtypes and the senescore of individual patients in TCGA cohort was depicted by the alluvial diagram (Fig. 5E). And the correlations between the senescore and varying biological processes suggested that senescore was positively associated with the activation of DDR, senescence, and DNA damage repair pathway (Fig. 5F). The senescore was significantly higher in cluster 2 subtype than that in cluster 1 (Fig. 5G). Whereas for the established molecular subtypes, the GS subtype had the lowest median senescore, while the EBV subtype and the MSI subtype showed significantly correlated with higher senescore (Fig. 5H). Similarly, the senescore was also associated with the MSI status, corresponding microsatellite instability-high (MSI-H) to the highest senescore (Fig. S5F). Furthermore, the senescore was positively correlated with the mutation density (Fig. 5I). Additionally, senescore signature was associated with histological classification, showing patients with intestinal GC had highest senescore (Fig. S5G). These findings suggested that the senescore was closely associated with clinical characteristics of GC patients and may serve as a prognostic biomarker.
Validation of senescore to effectively predict patient outcomes
To confirm the robustness of the senescore in different populations, we further explored ACRG cohort that was provided with the comprehensive clinical information. Similar with the TCGA cohort, the survival analysis revealed that the high senescore was associated with the longer OS as well as the RFS (p < 0.001; hazard ratio, 0.12 [95%CI, 0.04-0.33]) (Fig. 6A). And the prognostic nomogram demonstrated an impressive predicting ability for OS (1-year AUC = 0.81, 3-year AUC = 0.77, 5-year AUC = 0.75) as indicated by the time-dependent ROC curve analysis (Fig. 6B). Moreover, the significant differences of the senescore were also observed among molecular subtypes in ACRG cohort. Rather, the MSI subtype showed the highest senescore, while the MSS/EMT subtype had the lowest senscore (Fig. 6C). In addition, the histological classification also displayed the obvious differences of the senescore values (Fig. S6A), which was in consistent with the results in TCGA-STAD cohort. Additional validation of the senescore model was further confirmed in the validation cohort1 (n=1484) for OS and in the validation cohort2 (n=999) for RFS, as well as all GC patients in our study (Fig. 6D,E, Fig. S6B). Next, the prognostic power of senescore was examined in a wide spectrum of gastrointestinal tumours, and significantly validated in colorectal cancer and oesophageal cancer, respectively, suggesting senescore signature may be conserved in other types of cancer (Fig. S6C-G).
Senescore predicts the efficacies of chemotherapy and immunotherapy
Considering that cellular senescence could be induced by DNA damage drug, we explored whether senescore, the senescence scoring model, could predict the treatment efficiency of chemotherapy. In ACRG cohorts with available adjuvant chemotherapy information, patients were divided into high and low senescore group and the difference in OS and RFS was independently assessed. Adjuvant chemotherapy was found to markedly improve the overall survival rate and recurrence-free survival rate in patients with low senescore (P = 0.00037), while patients with the high senescore showed only a moderate benefit from adjuvant chemotherapy (P=0.08, Fig. 6F, S7A). Furthermore, patients with GC in another three cohorts of adjuvant chemotherapy were categorized into high- and low-senescore group. And significant survival difference was observed in the low-senescore group, but not in the high-senescore group (Fig. 6G, S7B), which was consistent with the results in ACRG cohort. These result suggests that the patients with low-senesecore were more inclined to gain benefits from the chemotherapy.
However, in contrast to chemotherapy, patients with high senescore, especially the highest quartile senescore, exhibited significantly prolonged OS in an independent anti-PD-L1 therapy cohort (n=348) (Fig. S7C, S7D) 21. In addition, patients responsive to anti-PD-L1 therapy displayed higher senesore (Fig. S7E) and senesore in this immunotherapy cohort was positively related with tumor mutation load (Fig. S7F), suggesting patients with high senescore was associated with better immunotherapy efficiency. Similarly, patients with high senescore in a small anti-PD1/PD-L1 cohort (n = 27) 22 displayed longer progress free survival (PFS) (Fig. 6H), as well as higher overall response (50% vs 21%) compared with low senesceore group (Fig. 6I). Similar findings were independently confirmed in a cohort of large size (n = 399), in which the patients of urothelial tumors 21 and melanoma 23 receiving anti-PD1/PD-L1 therapy were combined, where the high-senescore group displayed a significantly longer OS (Fig. 6J), and a higher proportion of overall response (31% vs 19%) (Fig. 6K). These results indicated that senescore might serve as a potential biomarker helping the clinicians to select appropriate patients for either chemotherapy or immunotherapy.
Validation of senescore at translational level
To further validate our model at translational level for the purpose of clinical pathology, a sophisticated method of multiplex immunofluorescence histochemistry, which allows simultaneous detection of multiple target proteins, was employed to analyse the protein expression of the 6 senescence signature genes identified above on a GC tissue microarray (TMA, Fig. 7A). Among these senescence signatures, four molecules, each by their own, were associated with clinical prognosis (Fig. S8A-F). Consistent with transcriptional analysis above, higher senescore in TMA cohort was significantly associated with longer OS (HR = 0.28, 95%CI: 0.17-0.45; p < 0.0001) (Fig. 7B). Furthermore, the ROC curve was conducted to evaluate the sensitivity and specificity of the senescore, and it illustrated that the AUC values was 0.882, 0.854 and 0.878 for 1, 3 and 5 years, respectively (Fig. 7C). In addition, ROC curve analysis also demonstrated the forceful predictive ability of the prognostic nomogram integrating senescore with clinical stage and age (Fig. 7D). More importantly, the senescore displayed stronger predictive value for OS in comparison with TNM stages, while the nomogram model combining the clinical stage and age with senescore only slightly elevated the predictive accuracy of the patients’ prognosis as compared with senescore alone (Fig. 7E). Collectively, senescore is promising in predicting clinical prognosis by feat of the immunofluorescence histochemistry strategy.