Scoring of senescence signalling in multiple human tumour gene expression datasets, identification of a correlation between senescence score and drug toxicity in the NCI60 panel and a pro-inflammatory signature correlating with survival advantage in peritoneal mesothelioma

Development of senescence scoring approach

Bypass of senescence is a requirement for full transformation. However, many cancer cells retain the ability to undergo senescence in response to genotoxic insults [3, 4]. Therefore, tumours must retain extant senescence pathways, though it is unclear how their expression is changed during the process of transformation in different cells and whether these expression levels might relate to the propensity to undergo cell cycle arrest in response to appropriate therapeutic intervention. In this study we have explored the use of a senescence scoring approach based on expression profiling of well-defined senescence markers as a means to evaluate latent senescence pathways in relevant samples. Published literature was mined for gene expression signatures and biomarkers of senescence (Table 1) and two sub-signatures were established. DAS biomarkers represent genes with known roles in senescence from literature involving DNA damage and chromatin responses [29–31]. Additionally, a modified secretory senescence (mSS) signature was created by combination of senescence messaging secretome/senescence associated secretory phenotype (SMS/SASP) signatures with a signature of 4 genes representing DNA damage and telomere dysfunction which increase in expression/activity in senescent cells and are detectable in serum of ageing human patients [11, 12, 14, 15]. At present our measure of the senescence response relies on relative expression of individual biomarkers. We propose that with an ever increasing number of genes associated with the senescence response it should be viewed as a new ontological process consisting of a complex network of gene expression and signalling pathways. To demonstrate this point Figure 1 shows known functional relations between the genes in Table 1, which are present in the MetaCore interactions database from GeneGo Inc [32]. Additionally, the figure shows representative GO processes for each gene, which further illustrates the cooperation of diverse processes of damage and inflammation in senescence.

To assign senescence scores, microarray data was imported into BRB ArrayTools [33] and filtered to show only those genes associated with senescence. Data was normalised to a median value across all arrays to make all data relative to a median value of 1. Therefore genes with expression greater than 1 were up-regulated and less than 1 down-regulated. Each gene was scored as senescent or not depending on whether expression was in a pro-senescence direction (see Table 1) with reference to the median expression level. For example, expression of the proliferation-associated gene KI67 below median would be marked "senescent" but above median expression would not. The number of genes marked "senescent" was then calculated for each sample and expressed as a percentage of the total number of genes expressed in each particular signature, all biomarkers in Table 1, DAS biomarkers only and mSS biomarkers only, to give an overall senescence score for each signature. By this approach, we hypothesised that the level of pro-senescence signalling in each pathway can be measured. Although very simple mathematically, taking the percentage of genes still signalling in a senescent manner provides a single continuous senescence measurement in human tumours incorporating the sum contributions of many pathways, something which has not previously been reported. In addition the scoring workflow could be used as a hypothesis driving tool to investigate the effect of modulating the expression of different genes in the senescence signalling network or indeed applied to other biological scenarios by integrating alternative expression signatures. We have utilised this scoring system to perform investigations of senescence signalling in various settings (Figure 2).

Senescence signalling pathways induced by radiation therapy

Our first main aim was to establish proof of concept as to whether our scoring approach can provide a reliable readout in settings in which senescence would be predicted to occur. We therefore investigated whether senescence signalling was altered in tumours and normal tissue of the colon before and after radiotherapy in the publicly available dataset GSE15781 (Figure 2). Examination of all biomarkers showed that tumours had a higher level of senescence signalling compared to their respective normal tissues (score 69% and 56%, respectively) (Figure 2A). Furthermore, irradiated tumours have higher median levels of pro-senescence signalling than their non-irradiated counterparts (score: 74% in irradiated tumours compared with 69% in non-irradiated). These results suggest that the scoring system is able to capture increased senescence signalling resulting from radiotherapy.

As well as measuring overall senescence signalling, we assigned scores for the sub-signatures DAS and mSS. Irradiation of tumours causes high levels of DNA damage with the aim of causing cell necrosis or apoptosis [34]. Figure 2B shows that irradiated normal tissue has increased median DAS signalling compared to non-irradiated normal, however the overall spread of the data is lower in the irradiated normal tissue (50%-87% in non-irradiated compared to 37%-75% in irradiated normal). Median DAS signalling in adenocarcinoma is unchanged by irradiation, although the minimum scores are at a higher level in irradiated tumours than their non-irradiated counterparts (62% in irradiated compared to 50% in non-irradiated adenocarcinoma), suggesting a trend in favour of higher DAS signalling in irradiated tumours. Interestingly, irradiation did not increase mSS signalling in tumours (score: 73% pre- and post-irradiation). In contrast, normal tissues showed a median increase of 6.7% following irradiation (Figure 2C). Radiotherapy therefore appears to specifically activate the DNA damage and chromatin aspects of senescence whilst not altering levels of the secretory senescence signalling pathway in tumours. This finding may present novel opportunities for improved senescence induction through therapeutic activation of secretory pathways in these tumours.

Senescence scoring in drug treated breast cancer cells

To assess the performance of the scoring approach in analysis of drug-induced accelerated senescence, we examined the public data set GSE1647 in which the authors investigated drug-specific toxicity-related expression profiles in a panel of breast cancer cell lines treated with chemotherapy drugs. Senescence was not directly assessed in this study, though each drug was administered at its IC50 for each cell line (assessed by a mitochondrial dye conversion assay). Because the set contains few replicates for each treatment, for statistical robustness we have here pooled data for the 3 cancer cell lines ZR-75-1, ME16C, and MCF7 treated for 24 hours with either doxorubicin or 5-FU, both of which have previously been observed to induce the accelerated senescence phenotype in tissue culture [35, 36].

Across all cell lines, the median basal senescence score is comparatively low (20%). As observed with radiotherapy of colon adenocarcinoma, treatment with either 5-FU or doxorubicin results in induction of the overall senescence score (Figure 2A, 5-FU score: 35%; doxorubicin score: 33%). DAS markers have basal expression of 14% in untreated cells. As expected, these genes are strongly induced, resulting in median scores of 45.5% and 41% for 5-FU and doxorubicin, respectively (Figure 2B). The mSS component is also induced by 24 hour drug treatments, though more modestly than the DAS component, rising to 31% and 28% from a basal level of 25%. Thus, both radiotherapy and chemotherapy induce multiple senescence markers as expected, though mSS appears to have a greater contribution in the chemotherapeutic setting.

Senescence signalling dynamics in human mesenchymal stem cells

Having shown that the scoring approach reports increased senescence signalling under therapeutically relevant conditions of genotoxic stress, we next investigated replicative senescence of human mesenchymal stem cells (hMSCs). We scored the senescence signalling signatures at passages 2 to 11 in hMSCs from the same patient (GSE9593). Prior to passage 11, at which time hMSCs were completely senescent, the authors of this study observed continuous changes in the cellular phenotype including gradual increases in cell size, autofluoresence, SA-bGal staining and osteogenic differentiation potential, while adipogenic differentiation potential slowly decreased [26]. These results are consistent with the model of partially stochastic dynamics in replicative senescence onset that has previously been proposed [37–39].

Our analysis of all markers (Figure 2A) is consistent with these observations, indicating that senescence gradually increases with passage. A plateau between passages 6 and 8 is evident when all markers are scored (score 79%). Strikingly, however, examination of the DAS and mSS signatures shows that this plateau is in fact a timed change in signalling type. Secretory senescence signalling remains level at 75% until passage 6 at which point a slight reduction to 69% is observed. At this time point there is a corresponding sharp increase in DAS signalling (75% to 92%). A subsequent spike in secretory (mSS) expression (score 81%) occurs at passage 7, followed by a return to base levels (Figure 2C). Examination of DAS markers at this time point (Figure 2B) shows a transient decrease in their expression (score 75% compared with 92% at passage 6) before returning to maximal levels. This time course corresponds with a transient change in the global expression profile of these cells observed between passages 6-8 in the original paper [26].

The dynamics of telomere shortening were not addressed in the original study. However, P7 occurs immediately before the onset of the main plateau phase of growth, suggesting that a significant fraction of cells may have undergone telomere dysfunction. DAS signalling spikes at P6, prior to the peak in mSS signalling at P7, raising the intriguing possibility that a threshold level of damage signalling triggers cells to communicate cellular stress to the surrounding microenvironment through the mSS pathway [40]. Notably, the spike in secretory signalling also correlates with a transient increase in expression of glycoprotein GPNMB/osteoactivin which was validated by QPCR in the original study [26]. GPNMB is an osteoclastic factor also implicated in endothelial cell adhesion and regulation of pro-inflammatory macrophage signalling, supporting the idea that signalling to the microenvironment is indeed increased at this time point [41, 42]. To our knowledge this is the first time such a temporal switch in senescence signalling has been documented, however the differential regulation of secretory and other senescence pathways during replicative senescence is also concurrent with our other recent findings in human T-cells during senescence bypass [43].

Senescence scoring in a Melanoma progression dataset

The above results suggest that the scoring system effectively captures increased activity of senescence pathways during both accelerated and physiological senescence. The expression levels of extant senescence pathways are likely to be altered by the processes of transformation and tumour progression and may reflect the route to immortalisation in individual tumour types. However, these questions have not yet been substantially investigated. We therefore examined the changes in senescence signalling during tumour progression in a melanoma dataset (GSE4587). Figure 2 shows levels of senescence signalling in benign nevi (score 55%, for all markers). Interestingly, an increased senescence score was found in melanoma when compared to benign nevi, regardless of the type of senescence signalling explored (DAS score: 58% versus 50%; mSS score 75% versus 63%). Although transformation has been achieved in these melanomas some aspects of the senescence program still appear to be active. The specific identities of induced senescence genes during tumour progression may point to the route of immortalisation for specific tumour types.

Comparison of primary melanoma and metastasis using all markers shows a slight increase in metastasis (melanoma score: 70%; metastasis score: 71.4%) (Figure 2A). However exploration of DAS and mSS reveals striking differences between the two tumour stages. Primary melanomas have a high contribution from mSS signalling and low levels of expression of the DAS biomarkers (75% and 58% respectively), while metastases have the opposite expression pattern (DAS score 75%; mSS score 69%) (Figure 2B and 2C). These results may suggest that the proinflammatory microenvironment produced by the secretory senescence signalling pathway is advantageous to a primary tumour in processes such as angiogensis. The low DAS senescence scores in primary tumours may reflect low levels of signalling shortly after immortalization and senescence bypass. After which, genomic instability may increase during the process of tumour progression leading to higher levels of DAS senescence signalling pathways in metastasis. Alternatively, DAS related signalling pathways may be induced during bypass of anoikis. In summary, these data suggest that DAS signalling increases during tumour progression. However these increases are clearly insufficient to induce cell cycle arrest.

Correlation between senescence score and drug toxicity in the NCI60 panel

The results above provide proof of concept that senescence scores obtained using these signatures conform with expectations given current understanding of replicative and accelerated senescence. To extend these results we investigated the possibility that latent senescence pathways may affect cellular responses to drugging. We first applied the scoring approach to gene expression data for the NCI60 cell line panel. Since extensive pharmacological characterisation data is available in addition to the expression profile in the dataset GSE5846 it is an ideal test set to examine the relationship between senescence signalling and drug response. Figure 3A-C show scores for all biomarkers (A), DAS (B) and mSS (C) in this panel with cell lines divided according to tissue of origin. For quality assurance purposes we also re-scored our own stock of two of the NCI60 cell lines, DU145 and HT29 and compared these results with those of the public dataset (figure 3D). Although we found somewhat lower scores overall, the trend toward higher expression of DAS than mSS is preserved in our hands in both lines (GSE5846 DAS scores: 75%, both cell lines; GSE5846 mSS scores: 56.3%, both cell lines; our cells DAS scores: 45.5% and 40.9%; our cells mSS scores: 33.3%, both lines).

Each tissue of origin had an individual scoring range profile. Comparison of scoring revealed clear divergence between DAS and mSS in some tissues. However in general DAS scores were higher than mSS. For example, leukaemia and colorectal cell lines generally scored more highly for expression of DAS markers than mSS (median scores: 66.7% (leukaemia), 75% (colorectal), DAS; 53.1% (leukaemia) and 50% (colorectal), mSS). Conversely CNS cell lines showed slightly higher mSS expression (mSS median score: 71.9%; DAS median score: 70.8%). These differences may reflect common routes to tumorigenesis and/or microenvironmental interactions in these tumour types. For example, increased secretion of pro-inflammatory and chemotactic molecules might be expected to select against blood-borne tumours whereas in solid tumours of immune privileged sites such as CNS, these may assist in processes such angiogenesis.

To explore the relationship between senescence signalling and drug induced toxicities in this panel, we analysed GI₅₀ data for ~1500 compounds, tested at least four times against each of the cell lines, as described in a previous study by Scherf et al [44]. One finding presented in the original study, is a 1376 gene signature capable of clustering the cell lines according to their drug sensitivity patterns assessed by GI₅₀. To determine if senescence signalling affects drug sensitivities, we performed regression analyses comparing the GI₅₀ scores of each drug across the cell panel with senescence scores for DAS and mSS. Examples of significant regressions using the DAS and mSS signatures are shown in figures 4A and 4B, respectively, for compounds NSC300288 (p = 0.0015, slope = -3.34) and NSC638279 (p = 0.0002, slope = -4.04). Examination of regression trends across the compound set revealed 78 and 328 compounds with unique significant relationships (regression p < 0.05) between growth inhibition and DAS or mSS expression respectively, with a further 5 compounds displaying a relationship to both signatures (figure 4C). To explore the contribution of each signature to drug resistance or sensitivity we examined direction and slope of each unique significant regression for DAS or mSS (Figure 4D shows the mean value across all significant compounds unique for each signature). This analysis revealed that cell lines with high DAS and high mSS are generally more sensitive to compound induced growth inhibition, as shown by the overall negative regression slopes (average DAS slope = -2.41; average mSS slope = -1.93). These results suggest that latent expression levels of these pathways may potentially have a significant impact on response to some therapies.

Compound activity prediction for significant DAS/mSS related compounds

To extend our analysis of the relation between DAS/mSS and drug sensitivity into prediction of target class sensitivities, we undertook activity modelling for the significant compounds using graph-theory connectivity indices to build a decision tree model based on chemical similarity with compounds of known activities in the classes GPCR agonist inhibition, GPCR antagonist inhibition, kinase inhibition, protease inhibition, PDE inhibition, ligand-gated ion channel inhibition and nuclear hormone receptor inhibition. We compared the fraction of assigned predicted activities for each target class in the DAS/mSS related compounds with those of the entire test compound set to obtain fold enrichment distributions as shown in figure 4E.

Protease inhibitor pharmacophores were enriched in the DAS related compounds compared with their representation in the test set (1.05-fold). Interestingly, since DAS apparently confers sensitivity to these compounds, these results may suggest that pharmacophores corresponding with existing protease targets may be therapeutically favourable in high DAS contexts. In all other drug categories mSS showed higher enrichment than DAS, where PDE inhibitor, Ion channel and GPCR agonist chemotypes were the most enriched in mSS related compounds (1.32-, 1.36- and 1.2-fold respectively). Since mSS also correlates with increased toxicity, intervention strategies based on such compounds may efficiently target mSS expressing tumours.

Differential senescence signalling patterns in mesenchymal tumours

Finally, to better understand the contribution of senescence signalling to the biology of mesenchymal tumours, we applied the approach to our own mesenchymally derived tumour gene expression dataset [28] (previously published under GSE17118 [45] (Figure 5). Considering all markers, the median scores were identical across all tumour types (46% (Figure 5A). Dissection of individual senescence signalling pathways showed that DAS signalling is lowest in MPNST (42%) and mesotheliomas (46%) and higher in liposarcoma (score: 50%) (Figure 5B). In contrast, mSS expression is higher in mesothelioma (score: 50%) than either MPNST or liposarcomas (44% and 47%, respectively) (Figure 5C).

These results suggest that even small changes in senescence signalling pathways can be detected using this scoring method and allow for hypothesis driven investigations. For example differences in extant/latent senescence signalling could affect factors such as patient prognosis and response to therapy.

Supervised Hierarchical clustering using the mSS signature highlights heterogeneous subgroups within mesenchymal tumours and a prognostic signature in peritoneal mesotheliomas

Given the essential nature of the establishment of a secreted inflammatory network of signalling in senescence induction [46] and our observations of differential expression patterns of secretory senescence in mesenchymal tumours and hMSCs, we used the mSS signature to further explore senescence signalling patterns in this pathway in mesenchymal tumours. After normalisation of gene expression array data, we examined the expression patterns of only the mSS signature within each tumour type using supervised hierarchical clustering (Figure 6A, and data not shown).

Analysis of GO processes for these prognostic gene groups revealed that group 1 genes have functions in immune activation (p = 4.1E-10), chemotaxis (p = 1.6E-11), inflammation (p = 1.2E-8) and negative regulation of cell proliferation (p = 3E-16). This is concurrent with our previous findings in other tumour types, that patients with high levels of local inflammation and immune invasion at the tumour site have improved survival [47]. Furthermore, the genes in group 2 are involved with the regulation of cell growth (p = 3.4E-16), proliferation (p = 2E-13) and cell motion (p = 1.3E-14). The increased expression of genes involved in cellular proliferation and cell motion may be indicative of more aggressive tumours leading to decreased survival rate.

We performed the same analysis on the other two mesenchymal tumour types. Similar gene groups also split MPNST into 2 tumour groupings while liposarcomas were divided into three groups with high group 1 expression (A), low group 1 expression (B) and low expression of both groups (C) (data not shown). However, examination of survival revealed no significant difference between the groups in these tumours (Figures 6C and 6D). Nevertheless, these expression groups may reflect biological factors which might warrant further investigation to improve understanding of the underlying biology and role of senescence in these subgroups. The same analysis was also carried out using the DAS signature, however no clear groupings were produced by hierarchical clustering that would facilitate Kaplan-Meier survival analysis.

Senescence scoring of gene expression data

Publicly available or in-house generated gene expression microarray data files were imported into BRB ArrayTools (BRB-ArrayTools developed by Dr. Richard Simon and BRB-ArrayTools Development Team [33]) and filtered to show only values associated with the senescence signature gene list in Table 1. Filtered data was then exported to Microsoft Excel (Microsoft Corp.) for further analysis. Data was normalised to a median value calculated across all arrays in that experiment so that data >1 was up-regulated and <1 down-regulated relative to a median of 1. Genes were scored as "senescent" or "non-senescent" depending on their gene expression levels relative to the median and their known expression pattern during senescence signalling (see Table 1). For example, expression of the proliferation-associated gene KI67 below median would be marked "senescent" but above median expression would not; likewise expression of the cyclin dependent kinase inhibitor p21 above median would be marked "senescent" but below median expression would not. The number of genes signalling in a "senescent" direction was then summed and represented as a percentage of the total number of genes expressed from the senescence signature. For example, when 20 of the 28 genes in Table 1 were expressed and 16 genes were signalling in a pro-senescent manner the score was 80%.

To calculate DAS and mSS scores data was sorted in Microsoft Excel and split into the two signature lists. Scores were calculated as a percentage of only DAS or mSS expressed genes respectively.

Median, maximum and minimum values were plotted using Microsoft Excel or Minitab (v15, Minitab Inc, Coventry, UK).

Correlations of senescence signature rankings

To assess the correlations of senescence scores given to each individual mesenchymal tumour by each signature the scores for each tumour were placed into Minitab (v15, Minitab Inc, Coventry, UK) and rankings for each signature score and Pearson correlations calculated. Regression analysis comparing NCI60 GI₅₀ data with senescence score was performed using the statistical package of Perl and exporting results to Excel.

Chemoinformatic prediction of compound activities

The dataset was pre-processed to remove couterions, assign physiological pH-appropriate atomic charges and add hydrogen atoms. Fifty 2D molecular descriptors were computed for each compound including Petitjean graph theory descriptors, Weiner path and polarity descriptors, GCUT descriptors, Kier and Hall connectivity indices, Zagreb index, Balaban indices, Van der Waals surface area and volume descriptors and topological polar surface area [51]. The descriptors were also computed for a training dataset of 3,537 publicly-available known active compounds collated from the PubChem database and the ChEMBL database (437 kinase inhibitors, 781 protease inhibitors, 806 GPCR agonists, 737 GPCR antagonists, 100 PDE inhibitors, 297 nuclear hormone receptor inhibitors and 379 ligand-gated ion channel inhibitors). This database of known actives was input to a binary decision tree algorithm, to derive a 7-category decision tree using 38/50 descriptors which correctly classified 62% of the dataset. The decision tree was subsequently used to predict activities of the NCI60 test dataset. All modelling conducted within the MOE software suite (Chemical Computing Group Inc, 1255 University Street, Montreal, Quebec, Canada).

RNA extraction and Gene expression microarray analysis of Mesenchymal tumours

A subset of 24 liposarcoma samples and 20 peritoneal mesotheliomas from a larger study population previously described [52, 53] and 16 MPNST were used for gene expression analysis. All samples were obtained under local ethical approval. RNA extraction and microarray hybridisation was performed as described in [28, 52, 53]. Fully MIAME compliant data files have previously been published under accession number GSE17118 [45].

Hierarchical Clustering and Survival Analysis

Hierarchical clustering of each tumour type was performed using Spearman correlation, merging branches with similarity less than 0.01 in Genespring GX (v7.3). Kaplan-Meier plots of survival and Log Rank (Cox Regression) tests adjusting for age and sex were performed on all resulting first level clusters of samples in all 3 tumour types using SPSS (v15, SPSS Inc, Chicago, Illinois, USA).