Senescence-associated morphological profiles (SAMPs): an image-based phenotypic profiling method for evaluating the inter and intra model heterogeneity of senescence

Establishing senescence models

A long-standing limitation within the senescence field is the lack of a universally accepted “marker” of senescence, with the collective assessment of a series of so called “hallmarks” widely accepted as best practice for a reliable classification [9, 19]. Our primary aim within this work, was to explore the intricacies of senescent cell morphology through phenotypic profiling, to understand whether such profiles could constitute a senescence hallmark themselves. An important first step in this process, was to characterise a range of established senescence models via conventional markers. We focused on exploring existing high-content images, as opposed to creating bespoke experiments for our phenotypic profiling, as we believe our approach could be readily implemented by other research groups wishing to extract profiling data from established datasets. This had the advantage of providing an inherent level of variability within the data, as it was generated by independent researchers for unrelated primary research aims. This approach also supported another of our aims: identifying a means of senescent cell classification that does not require resource intense optimisations and which may be applied in most research settings with microscopy capacity. The models employed here comprise four for which we have previously published characterisation data; oncogene-induced senescence (OIS) in IMR90s [20], paracrine senescence via conditioned media from OIS IMR90s [20], replicative senescence in adult human mammary fibroblasts (HMFs) [20, 21] and OIS in human mammary epithelial cells (HMECs) [22]. These image stacks were re-mined for this work, with previous senescence characterisation summarised in Supplementary Figure 1. Furthermore, we have also explored a model of replicative senescence in adult human dermal fibroblasts (HDFs), characterising canonical senescence markers including; loss of cellular proliferation by population doublings and Ki67 staining (Figure 1A and 1B) expression of the tumour suppressor p21 (Figure 1C), appearance of γH2AX DNA damage foci (Figure 1D) and increased SA-β-Gal activity (Figure 1E). As a supporting model, we also characterised UVB-induced senescence in HDFs via assessment of proliferation (Supplementary Figure 2A and 2B) and detection of DNA damage foci (Supplementary Figure 2C). Therefore, in each of the six models, cells were classified as either proliferating or senescent according to conventional senescence hallmarks.

Figure 1. Characterisation of senescence in human dermal fibroblasts (HDFs). (A) Hayflick proliferation curve for human dermal fibroblasts (HDFs) from early proliferation (EP) to deep senescence (DS) through serial cell culture. N = 2–6. (B) Immunofluorescence staining of DAPI (blue) and Ki67 (green) in EP and DS HDFs. N = 2. Scale bar = 50 μm; (C) Immunofluorescence staining of DAPI (blue) and p21 (green) in EP and DS HDFs. N = 2. Scale bar = 50 μm (D) Immunofluorescence staining of DAPI (blue) and γ-H2AX foci (green) in EP and DS HDFs. N = 3. Scale bar = 50 μm. (E) Brightfield assessment of SA-β-Gal (blue) in EP and DS HDFs. N = 2. Scale bar = 100 μm.

Phenotypic profiling of senescent cells – Z-Score profiling

The initial stages of phenotypic profiling comprise the measurement of a large number of discrete image-based features, representing size, texture, spatial and intensity measurements from both nuclei and cells [10]. These are then combined to create a total profile for a particular treatment or cellular state. Here, we utilised two simple and widely used stains (DAPI and Cell Mask) to visualise the nuclear and cellular morphology of proliferating and senescent cells in each of our models. High-content image analysis was then performed, to generate measurements of 62 image-based features (Table 2) for every cell in each model. Means were then generated from both technical and experimental replicates (Figure 2A). Representative images from each model, along with nuclear and cellular masks upon which analysis was based, are summarised in Supplementary Figure 3.

Figure 2. Phenotypic profiling of senescence via Z-score profile generation. (A) Schematic overview of data processing proceeding Z-score generation for each senescence model (P: proliferating condition, S: Senescence condition). (B–E) Z-score profile heatmaps for oncogene-induced senescence (OIS), paracrine senescence (Paracrine), HMF replicative senescence (HMF RS) and HDF replicative senescence (HDF RS) models. Y-axis comprises 62 morphological features (Red = positive modulation, Blue = negative modulation), White = no change). Proliferating conditions: vector induction (Vec_Ind), vector conditioned media (Vec_CM), HMF early proliferating (HMF_EP) and HDF early proliferating (HDF_EP). Senescence conditions: OIS induction (OIS_Ind), OIS conditioned media (OIS_CM), HMF deep senescence (HMF_DS), HDF deep senescence (HDF_DS). (F) Summary Z-score profile heat map of senescence conditions from each senescence model. Y-axis comprises 62 morphological features.

Next, Z-score normalisation was performed in order to compare features that would typically be on vastly different scales and, importantly, between senescence models [23]. This method centres the data to the mean value of the proliferating control in each model and scales based upon the standard deviation of the control. Thus, a Z-score of 1 represents a single standard deviation from the mean of the control condition. This was used within each model to produce a profile of Z-scores for each senescent condition. These were displayed as heatmaps, with positive and negative modulation from the control indicated as red or blue, respectively (Figure 2B–2E). In each of the senescence models, the senescent condition appeared distinct from the control in a number of features. However, the quantity and magnitude of these changes were highly variable between models. The OIS model produced the most striking profile, with potent change from the control condition in 49 out of 62 features. The paracrine senescence profile appeared to have less potent changes, but featured significant alterations in 50 out of 62 features. Both replicative senescence models produced more subtle profiles, with significant alterations in 38 out of 54 features for the senescent HMFs (8 features could not be assessed due to use of an older generation of microscope) and 37 out of 62 features in the senescent HDFs, albeit where changes did occur, they did so potently. This is visualised in the summary heat map (Figure 2F). Of note, the direction of modulation of many features appears consistent between models, opening the possibility that the profiles reflect common senescence-associated changes. However, the Z-score method of normalisation has the disadvantage that the profiles produced are relative to the control condition, thereby masking the “profile” of the proliferating cells. To overcome this limitation, we revised our data processing steps to generate meaningful profiles for both proliferating and senescent conditions.

Phenotypic profiling of senescent cells – standard score profiling

In order to provide a means of data scaling to facilitate comparisons between features and senescence models, the single target data was used to produce “standard scores” [10]. These utilise the mean and standard deviation values from the combined proliferating and senescent single target data within each model, to produce a data set comprising standard scores for each of the 62 features, for every cell in each model (Figure 3A). Median scores for each proliferating and senescent condition were then used to produce heatmap profiles. These supported the observations made using Z-scores, in that the senescent condition in each model was distinct from its proliferating counterpart (Figure 3B–3E). However, whereas some features were indistinguishable from the proliferating control by Z-scores, the standard scores demonstrated modulation in the vast majority of features in each model between the proliferating and senescent states. Similar observations to those of Z-score normalisation were also made in terms of the magnitude of change between models, with the OIS model producing the most potent phenotype and the replicatively senescent HDFs a more subtle one. When all senescence models were compared, hierarchical clustering demonstrated a clear distinction between the senescent conditions and those of the proliferating controls (Figure 3F).

Figure 3. Phenotypic profiling of senescence via standard score profile generation. (A) Schematic overview of data processing proceeding standard-score generation for each senescence model (P: proliferating condition, S: Senescence condition). (B–E) Standard-score profile heatmaps for oncogene-induced senescence (OIS), paracrine senescence (Paracrine), HMF replicative senescence (HMF RS) and HDF replicative senescence (HDF RS) models. Y-axis comprises 62 morphological features (Red = positive modulation, Blue = negative modulation), White = no change). Proliferating conditions: vector induction (Vec_Ind), vector conditioned media (Vec_CM), HMF early proliferating (HMF_EP) and HDF early proliferating (HDF_EP). Senescence conditions: OIS induction (OIS_Ind), OIS conditioned media (OIS_CM), HMF deep senescence (HMF_DS), HDF deep senescence (HDF_DS). (F) Summary standard-score profile heat map and hierarchical clustering of all proliferating and senescence conditions. Y-axis comprises 62 morphological features.

Therefore, it appears that phenotypic profiling may be used to clearly distinguish between senescent and proliferating cells based upon the development of a “senescent-associated morphological profile”. This profile is present within each senescence model but varies in terms of magnitude. Furthermore, these complex profiles seem to consist of both a “core signature” of features that change commonly within senescence, as well a number of idiosyncratic features that reflect the established heterogeneity between inducers of senescence. Alternatively, this could be the result of the different tissues from which the cells were derived or, more likely, a combination of the two. In order to explore this further, we utilised a model of DNA damage-induced senescence in the HDFs via UVB treatment (Supplementary Figure 2A–2C). This led to more potent profiles (both Z-score and standard score) than observed in RS HDFs, supporting the observation from the IMR90s that different inducers of senescence in the same cells produce distinct phenotypic profiles (Supplementary Figure 2D and 2E). In order to emphasise this further, we then applied our analysis pipeline to OIS HMECs, as to this point we had only evaluated fibroblasts (Supplementary Figure 4). As with the IMR90 model, the HMECs were associated with strikingly potent profiles in the OIS condition, both demonstrating the utility of our analysis pipeline in other cell types and further supporting the utility of phenotypic profiling as a powerful tool for senescence characterisation. Importantly, the utilisation of standard scores provided significant benefit over Z-scores, in understanding the nuance between each phenotypic profile. This included the ability to interpret the profile of features in the control condition, allowing comparison of these between models.

Phenotypic profile generation by exploratory factor analysis

EFA model development

Whilst the standard score profiles facilitated the distinction between proliferating and senescent cells in a number of models, the high dimensionality of the data made biological interpretability difficult, as well as any comparison to established senescence-associated features. To overcome this, we adopted a strategy similar to that of Young et al., utilising exploratory factor analysis (EFA) as a method of dimensionality reduction [11]. EFA seeks to determine whether the changes observed in measured features occur through alterations to so called “latent factors” [24, 25]. In the context of HCA, these are assumed to be cellular properties which, whilst not directly assessed, are of potentially more biological relevance than the measured features themselves [11]. Therefore, EFA has the potential to complement HCA by identifying the underlying relationships between features, allowing them to be considered collectively. Here, we utilised EFA to understand whether a select group of biologically interpretable latent factors could provide insight into the morphological phenotypes of senescent cells.

First, the number of latent factors underpinning the full list of 62 features needed to be established. Initially, data from each individual senescence model was used in isolation, to construct four unique EFA models. Scree plots were produced in each instance to determine the number of latent factors. These represent the total variance accounted for by each potential latent factor as an eigenvalue, with the total number of possible latent factors equal to the total number of features (i.e., 62). By employing the widely accepted Kaiser criterion, only factors that had an eigenvalue greater than 1 were retained for factor analysis within each model [26, 27]. This resulted in the following totals of latent factors: OIS induction (8), paracrine senescence (7), RS HMF (7) and RS HDF (9) (Figure 4A). Next, factor loadings were determined for each latent factor. This refers to how closely individual features are associated with each latent factor, and can be considered similar to standard regression coefficients. A stringent loading threshold of 0.5 was selected, so that only features that loaded with a value greater than this were included in each factor. For illustrative purposes, factor loadings for a single factor from the OIS induction EFA model, has been visualised as a polar plot (Figure 4B). Next, the feature composition of individual latent factors in each EFA model was assessed in order to assign potential factor classifications. Interestingly, in each EFA model, the latent factors could be broadly classified according to the following designations 1) Nuclear Size, 2) Cell Size, 3) Cell Intensity, 4) Nuclear Intensity, 5) Nuclear Spacing, 6) Nuclear Shape, 7) Cell Shape 8) Cell Intensity (variance) and 9) Nuclear Intensity (variance), based on the composite features. As EFA represents an iterative tool for investigating relationships within data, as opposed to one applying inflexible statistical cut offs, the total number of latent factors was adjusted to 8 in each EFA model at this point, as initial factor loadings indicated “under-“ and “over-factoring” in the models utilising 7 and 9 latent factors, respectively (see Table 1). By employing this iterative approach, EFA models were constructed for each senescence model, with 7 common factor designations, along with one of either Cell Intensity (variance) or Nuclear Intensity (variance). In order to visualise the composition of factors within each EFA model, as well as the factor classifications themselves, schematic representations were produced (Figure 4C and Supplementary Figure 5).

Figure 4. Exploratory factor analysis (EFA). (A) Scree plots for oncogene-induced senescence (OIS), paracrine senescence (Paracrine), HMF replicative senescence (HMF RS) and HDF replicative senescence (HDF RS) models. Red line indicates eigenvalue = 1. (B) Polar plot of factor loading values for factor 1 from OIS EFA model (designated Nuclear Size). 1–62 refer to features (Table 2). Blue shaded area indicates factor loading threshold of 0.5. (C) Factor loading diagram for OIS EFA model with factor designations. Equivalent diagrams for other senescence models are found in Supplementary Figure 3).

EFA model analysis

Following EFA model development, factor scores were extracted for each cell within each senescence model and median score profiles generated (Figure 5A–5D). Once again, these clearly demonstrated a distinction between the proliferating and senescent condition within each senescence model. However, they also provided insight into the relative importance of each latent factor in determining this distinction. For example, Cell Shape appears to contribute less to the overall phenotype of the paracrine senescent cells, as opposed to the other three models. Conversely, Cell Intensity changes dramatically in the paracrine senescent cells but much less so in the RS HMFs. Therefore, EFA aids the assessment of phenotypic profiles by reducing the overall complexity and directing attention to specific interpretable factors within each model. It is important to note that the factors can only cautiously be compared between models in this way, as whilst similar, the underlying features which contribute to each latent factor (as well as the relative factor loadings) are variable. In order to overcome this, EFA was performed using the combined single target data from all senescence models. Scree plot and eigenvalue assessment suggested that a model comprising 8 latent factors should be constructed, supporting the earlier iterative adjustment to this latent factor threshold (Supplementary Figure 6A). Furthermore, each factor could be identified according to the designations previously established (Supplementary Figure 6B and 6C). Factor scores were again used to construct profiles for each condition, which could now be directly compared, having been based on the same EFA model (Figure 5E). As with the standard score profiles, hierarchical clustering demonstrated a clear distinction between the senescent and proliferating conditions, albeit with the development of highly heterogeneous senescence-associated factor profiles (Figure 5E). This demonstrated that EFA does not preclude phenotypic profiling-based senescence identification despite the large reduction in profile components.

Figure 5. Factor score profiles following EFA. (A–D) Factor score profile heatmaps for oncogene-induced senescence (OIS), paracrine senescence (Paracrine), HMF replicative senescence (HMF RS) and HDF replicative senescence (HDF RS) models. Y-axis comprises scores for 8 latent factors determined by EFA (Red = positive modulation, Blue = negative modulation), White = no change). Proliferating conditions: vector induction (Vec_Ind), vector conditioned media (Vec_CM), HMF early proliferating (HMF_EP) and HDF early proliferating (HDF_EP). Senescence conditions: OIS induction (OIS_Ind), OIS conditioned media (OIS_CM), HMF deep senescence (HMF_DS), HDF deep senescence (HDF_DS). (E) Summary factor score profile heat map and hierarchical clustering of all proliferating and senescence conditions. Y-axis comprises scores for 8 latent factors determined by EFA.

The primary aim of EFA was to identify relationships between features and reduce the complexity of the phenotypic profiles in order to aid interpretation. Whilst we have taken an unbiased and systematic approach, the consistency of latent factor designations between senescence models suggests that extraction of a set of select features could be used as a more straightforward overview of the broader profiles (Supplementary Figure 7A). These representative profiles are less complex to produce and maintain distinction between senescent and proliferating cells both within individual models (Supplementary Figure 7B–7E) and when comparing all models (Supplementary Figure 7F). Such profiles could serve as more accessible (though still valid) representation of senescence-associated morphological profiles, particularly as the majority of the composite features have established literature precedent. Therefore, EFA complements the phenotypic profiling of senescent cells by providing several options for reduced profile complexity and interpretability. However, the reduced factor profiles also limit the potential for more nuanced comparison between models, which we consider a key advantage of the more complex phenotypic profiles. Therefore, recognising a major strength of the initial approach, the highly dimensional nature of the data produced by HCA, we sought to explore the heterogeneity of profiles generated on a single target level.

Exploring the heterogeneity of senescence-associated morphological phenotypes

Both the Z-score and standard score profiles demonstrated that in each model of senescence, a distinct senescence-associated morphological phenotype could be observed. However, each model was characterised by a degree of heterogeneity, both in terms of the magnitude and direction of change in the features observed. We classify this as inter-model heterogeneity, likely representing the distinct pathways engaged in the development of senescence via different inducing stimuli. Another key consideration developing within the field is that of intra-model heterogeneity, the development of distinct populations of senescent cells within a single model. The most well described example of this phenomenon is that of NOTCH1 mediated juxtacrine secondary senescence during OIS induction [28, 29]. Therefore, one would anticipate a degree of heterogeneity within each model, as different populations of cells are induced by distinct stimuli. To explore this, the single target standard scores were explored in a manner similar to that more commonly seen during single-cell RNA sequencing experiments. T-distributed stochastic neighbour embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP) algorithms were employed in order to represent the single target data in an interpretable form. Within these plots, each point represents the high-content, 62-feature profile of a single cell (Figure 6A–6D). Strikingly, in each model, regardless of clustering method, the population of senescent cells was distinct from that of the proliferating control. In the OIS model, two well-defined populations of senescent cells could be observed by both t-SNE and UMAP, aligning with the concept that OIS has both a primary, oncogene-driven component and a second, juxtacrine component, although comprehensive characterisation of NOTCH1 expression would be required to confirm this. The paracrine senescence model had a far more homogenous population of senescent cells, possibly due to the nature of the inducing stimuli being treatment with a SASP comprising a consistent composition. Interestingly, whilst primary senescent clusters were reasonably distinct in both replicative senescence models, this was less apparent than seen in the IMR90 models. This could be attributed to previous observations regarding the complex set of stimuli underpinning replicative senescence, including a degree of paracrine senescence, mitochondrial dysfunction, oxidative damage and telomere shortening within the total population [30]. Furthermore, as evident from the HDF Hayflick curve (Figure 1A), replicative senescence is a far more gradual process than seen in the OIS and paracrine forms of premature senescence. Thus, a greater degree of heterogeneity might be expected, given some cells would likely have reached senescence potentially weeks before others. Overall, this data aims to serve as an alternative visualisation of the morphology profiles, demonstrating variability in the phenotypes at the single cell level, which is not achievable when represented as a heatmap. We suggest that HCA shows good potential in facilitating the opportunity to explore potential heterogeneity present both between and within senescent models, as well as reinforcing our earlier observation that senescent cells are associated with distinct and complex morphological phenotypes.

Figure 6. t-SNE and UMAP profiles for single target morphological profiles. (A–D) t-SNE and UMAP plots for oncogene-induced senescence (OIS), paracrine senescence (Paracrine), HMF replicative senescence (HMF RS) and HDF replicative senescence (HDF RS) models. Proliferating conditions (Blue): vector induction (Vec_Ind), vector conditioned media (Vec_CM), HMF early proliferating (HMF_EP) and HDF early proliferating (HDF_EP). Senescence conditions (Orange): OIS induction (OIS_Ind), OIS conditioned media (OIS_CM), HMF deep senescence (HMF_DS), HDF deep senescence (HDF_DS).

Applying phenotypic profiling to senescence in vivo

Despite the emergence of a variety of hallmarks in vitro, the availability of reliable in vivo markers of senescence remains a challenge for the field. In order to understand whether our observation of in vitro senescence-associated morphological profiles were also detectable in vivo, we assessed the morphology of individual cells within human SK-MEL-103 xenograft tumours following treatment with palbociclib. We utilised p21 staining to define senescent cells within each treatment and confirmed that the percentage of p21 positive (i.e. senescent) cells increased following palbociclib treatment (Figure 7A–7C). When all cells in each treatment condition were assessed collectively, a subtle yet distinct phenotype comprising 27 morphological features (Table 3) distinguished the palbociclib condition from the control, suggesting such profiles could have utility as a prognostic marker of senescence in vivo, or form the basis for a classification tool (Figure 7D). More strikingly, when p21 positive and negative cells were considered separately, within each treatment condition a potent difference in morphological profile was apparent between the senescent and non-senescent cells (via both Z-score and standard score; Figure 7E and 7F, respectively). Hierarchical clustering of standard scores confirmed this distinction, suggesting that the alterations in morphological profiles we observe between senescent and non-senescent cells is not limited to in vitro settings. Importantly we utilised open source software to perform this analysis (QuPath), demonstrating that this methodology could be widely applied within the senescence research community [31].

Figure 7. Assessment of senescence morphology in palbociclib treated human tumour xenografts. (A and B) Analysis mask generation to define p21 positive (red) and negative (blue) nuclei in control and palbociclib treated xenografts. Blue stain = Counterstain, Brown stain = p21. Scale bar = 200 μm. (C) Percentage p21 positivity of nuclei in control and palbociclib treated xenografts. (D) Mean standard score profiles of control and palbociclib samples. Y-axis comprises 27 morphological features (Red = positive modulation, Blue = negative modulation), White = no change). (E and F) Z-score and standard score profiles of p21 positive and negative cells in control and palbociclib treated xenografts. Y-axis comprises 27 morphological features (Red = positive modulation, Blue = negative modulation), White = no change). N = 3.

Cell culture and reagents

Unless indicated, all reagents were sourced from Sigma, UK. IMR90 ER: STOP (vector) and ER:RAS (OIS) fibroblasts were produced as described previously [20, 37] and were kindly gifted by Juan Carlos Acosta (MRC Institute of Genetics & Molecular Medicine, Edinburgh). All cells were maintained in Dulbecco’s Modified Eagles Medium (DMEM; Life Technologies, UK) with additional 10% foetal bovine serum (FBS, Labtech.com, UK) and 2mM L-glutamine (DMEM; Life Technologies, UK). Human mammary fibroblasts (HMFs) were further supplemented with 10 μg/ml bovine pancreas insulin and were kindly provided by Martha Stampfer (Lawrence Berkeley National Laboratory, Berkeley). Tamoxifen inducible ER: RAS (OIS) telomerase-immortalised human mammary epithelial cells (HMEC) were cultured as previously described [22]. Human dermal fibroblasts (HDFs) were anonymously donated by healthy patients (LREC No. 09/HO704/69). Cells were cultured at 37°C/5% CO₂ without the presence of antibiotics. All experimental end points utilised cells at sub-confluent densities.

Antibodies

The following antibodies were used in for immunofluorescence staining: p21 (2947, Cell Signalling, UK; 1:500), Ki67 (NCL-Ki67p, Novocastra, Cell Signalling, UK; 1:1,000), γH2AX (05-636, Upstate Cell Signalling Technologies; 1:2,000), 53BP1 (A300-273A, Bethyl laboratories, 1:200), goat anti-mouse-Alexa Fluor 488 (A32723, Thermo Fisher, UK; 1:500) and donkey anti-rabbit-Alexa Fluor 546 (A10040, Thermo Fisher, UK; 1:500). Immunohistochemistry for p21 was performed using mouse anti-human p21 (M7202, Agilent Dako, US, 1:100), recombinant IgG1+IgG2+IgG3 bridging antibody (ab1334569, Abcam, UK, 1:500) and OmniMap anti-Rb HRP secondary antibody (5269679001, Roche, according to the manufacturer’s instructions).”

Senescence induction

Senescence induction was performed as described previously for both the OIS, paracrine senescence and HMF RS models [20–22]. HDFs were serially cultured to replicative senescence (>200 days) and designated as having reached deep senescence (DS) at passage 39 plus an additional 3 weeks (~280 days). For phenotypic profiling by HCA, early proliferating (EP) HDFs were seeded at 4,000 cells/cm², whilst DS HDFs were seeded at 10,000 cells/cm². For UVB-induced senescence, HDFs were seeded at 5,000 cells/cm² in 10 cm dishes (Starlabs, UK), with the control seeded at 4,000 cells/cm². Cells were incubated for 2 days in standard culture medium, and UV-treated for 5 consecutive days at a dose of 0 mJ/cm² (Control) or 6 mJ/cm² (UVB) (UV-B lamp G8T5E, Sankyo-Denki). During UVB treatment, normal culture medium was replaced with D-PBS (0.9 mM CaCl₂, 0.5 mM MgCl₂-6H₂O; Life Technologies, UK). The cells were the incubated in standard growth medium for a further 2 days before fixation and characterisation.

Immunofluorescence staining and high-content analysis (HCA) microscopy

Immunofluorescence staining and HCA microscopy was performed as described previously [20, 21]. Briefly, following fixation with 3.7% paraformaldehyde, cells were permeabilised with 0.1% Triton X-100 and stained with diamidino-2 phenylindole (DAPI) (Sigma UK, D8417, 1:1,000) and HCS Cell Mask Deep Red (Thermo-Fisher UK, C10046, 1:50,000 – 1:100,000) for 2 hours at room temperature. Cells were then imaged using an IN Cell 2200 automatic microscope (IMR90s and HDFs) or an IN Cell 1000 automatic microscope (HMFs). For senescence marker assessment, cells were blocked for 30 minutes with PBS/BSA following permeabilsation and stained with primary antibody diluted in PBS/BSA for 2 h at room temperature or overnight at 4°C. A PBS/BSA wash step was then used, followed by secondary antibody staining for 2 h at room temperature.

High-content analysis (HCA) and data processing

Images from all senescence models were analysed via InCarta image analysis software (Cytiva). DAPI staining was used to visualise nuclei, whilst Cell Mask staining allowed cellular visualisation. InCarta detection protocols were developed in order to generate nuclear and cellular image masks, upon which data analysis was performed (Supplementary Figure 3). These were consistent within individual experiments, with a level of fine tuning between replicates to ensure optimal mask generation. Image masks were then analysed according to pre-built features (termed measures in the InCarta software), with all features selected apart from those that assess “background” intensities, centre of gravity co-ordinates and major axis angles. These were excluded as they are principally tools for assessing imaging consistency as opposed to biologically profiling cells as opposed to specific cellular properties. This leads to a final list of 62 features, which formed the basis of the phenotypic profiles (Table 2). It is important to note that as we used bespoke HCA software, this is a lower number of total features than is sometimes used for phenotypic profiling, due to a lack of redundancy amongst features [12]. However, phenotypic profiling has been successfully employed with a feature count as low as 35, leading us to believe our profiles to be sufficiently comprehensive [11].

Following HCA, data normalisation was required in order to facilitate comparison between features on different scales. This was initially performed via Z-score normalisation compared to a proliferating control condition using the following equation, as used elsewhere [20, 21, 23]. Z-Score = mean value of three independent experiments for senescent experimental condition – mean value of three independent experiments for control condition / Standard Deviation (SD) of control condition. Z-score profiles were displayed as heat maps, where significant positive or negative change was defined as a Z-score greater than 1.96 and visualised as either red or blue, respectively. Subsequently, single target data was used to generate standard score feature profiles for individual cells from each condition in all models according to the following equation. Standard Score = value of an individual cell – mean value of all targets (proliferating and senescent) / Standard Deviation (SD) of all targets (proliferating and senescent). Standard scores were generated within experimental replicates. Standard scores from all three replicates per condition were then combined and median profiles calculated (to reduce the influence of outliers). These were presented as heat maps, with positive modulation of features represented as red and negative modulation of features represented as blue.

It is important to note that HMF and HMEC image acquisition was performed via an older generation of microscope (INCA 1000), which produces rectangular images as opposed to the square images produced by the INCA2200 that all other models were assessed using. This means that 8 features (Features 33,34,35,36,66,67,68 and 69 – Table 2) could not be assessed in this model. They have been given a consistent value of 0 to allow comparison of the whole profile to other models.

Exploratory factor analysis (EFA)

An iterative approach to EFA was taken as described in detail in the primary results section. EFA was performed in R, using the “psych” and “stats” packages, with EFA models constructed using combined proliferating and senescent single cell data from each senescence model, as well as with all single target data from all models. Scree plots were generated using factor eigenvalues, and the number of latent factors selected based upon the Kaiser criterion (eigenvalues >1), although this was iteratively assessed. Other parameters included extraction via maximum likelihood and an Oblimin oblique rotation method. Factor loadings were determined through use of a 0.5 loading threshold. Factor scores were calculated using Thompson’s regression method (scores parameter set to regression) and median profiles for each condition displayed as heatmaps.

Palbociclib treatment of human tumour xenografts

The experimental procedure was approved by the Ethical Committee for Animal Experimentation at the Parc Cientific de Barcelona and by the Government of Catalunya. Nude mice (Hsd:athymic nude-Foxn1^nu) were purchased from Envigo. The animals were kept under a 12–12 h light-dark cycle and allowed unrestricted access to food and water. 10⁶ SK-MEL-103 cells were injected subcutaneously in the dorsolateral flank of 8-week old male athymic nude mice. When tumours became visible at day 8–10 after injection, the mice were randomly assigned to the control or treated experimental groups. Treated mice received 100 mg/kg palbociclib in 50 mM sodium lactate by oral gavage for seven days. Control mice received vehicle only. Once the treatment was complete, the mice were euthanized. The tumours were extracted, fixed with 10% neutral buffered formalin for 24 hours at 4°C and embedded in paraffin for further processing.

Analysis of palbociclib treatment of human tumour xenografts

Control and palbociclib treated samples were analysed via QuPath open source software [31]. Image type was set as Brightfield (other), with counterstain detected using values 0.651, 0.701 and 0.29 and p21 detected using values 0.269, 0.568 and 0.778. All nuclei were detected via counterstain OD intensity threshold of 0.05 and p21 intensity threshold parameters set at 0.2, 0.4 and 0.6 (1+, 2+ and 3+, respectively) to determine p21 positivity. All features were then extracted apart from those measuring p21 intensity, giving a final list of 27 features (Table 3). Mean standard scores for each feature were then generated using all cells in each treatment condition and presented as heatmaps as described above. Cells identified as p21 positive and negative were then separated and Z-score and standard score profiles generated as described above.

R analysis

All heatmaps were constructed using the “heatmap.2” package, with distance matrices constructed using Euclidean distances. Hierarchical clustering was then performed via the ward.D2 method. T-distributed stochastic neighbour embedding (t-SNE) dimensionality reduction was performed via the “Rtsne” package with a consistent perplexity parameter of 100 used in each case. T-SNE plots were constructed using the “ggplot2” package. Uniform manifold approximation and projection (UMAP) dimensionality reduction was performed using the step_umap function of the “embed” package. UMAP plots were then constructed using the “ggplot2” package.

Senescence-associated beta galactosidase assay

Human dermal fibroblasts were fixed with 0.2% glutaraldehyde for 5 min and incubated at pH 6.0 with 1 mg/mL X-gal (10234923, Fisher Scientific, UK), 150 mM NaCl, 2 mM MgCl₂, 5 mM K₃Fe(CN)₆, 5 mM K₄Fe(CN)₆ and 40 mM NaPi. This was performed without CO2 for 24 hours at 37°C. Cells were then imaged using a Nikon brightfield tissue culture microscope.

Statistical analysis

Statistical assessments were performed using Graphpad Prism 7. Experiments using two groups were assessed via unpaired Students t-tests and multiple groups via ordinary one-way ANOVA with multiple comparisons. P values were represented as: ^*P < 0.05; ^**P < 0.01; ^***P < 0.001; ^****P < 0.0001. Error bars represent SD of 2 independent experiments unless otherwise stated.