Analysis of individual cells identifies cell‐to‐cell variability following induction of cellular senescence

Although cellular senescence is marked by strong phenotypic and gene expression changes at the population level, less is known about how individual cells contribute to the overall gene expression pattern. We therefore induced senescence in IMR‐90 human fibroblasts using the clastogen bleomycin (Fig. S1A) (Orjalo et al., 2009) and analyzed individual senescent cells using a commercially available platform from Fluidigm (C1) or manual isolation (MI). To control for growth status, we cultured nonsenescent cells for 3 days in 0.2% FBS to induce quiescence (Fig. S1A). Senescence was confirmed by senescence‐associated beta‐galactosidase (SA‐Bgal) activity (Dimri et al., 1995), markedly reduced EdU incorporation into nuclear DNA (a measure of nuclear DNA replication) (Fig. S1B), and secretion of a major SASP factor, interleukin‐6 (IL‐6) (Fig. S1C).

C1 isolation captured 57 (59%) single senescent cells and 79 (82%) quiescent cells. The different capture efficiencies were like due to the increased size and irregular shape of senescent cells, resulting in several empty wells. In addition, we manually isolated (MI) 163 single senescent and 156 single quiescent cells.

We analyzed the cells for expression (mRNA levels) of 96 genes, including those associated with cell cycle arrest (12 genes), nuclear lamins (four genes), mitochondria and energy production (12 genes), the SASP (46 genes), signal transduction (four genes), and other senescence regulators (11 genes), in addition to seven control genes, the expression of which were expected to be relatively invariant (Data S1). We chose these genes mainly based on prior knowledge of gene expression associated with cellular senescence. Of the 96 genes analyzed, we excluded 25 as either nonamplifying (no or late Ct) or nonspecific (multiple melting curves), resulting in a dataset comprised of 71 gene expression (mRNA) levels in 235 quiescent and 229 senescent cells. These 71 transcripts were detected in at least 30% of the analyzed cells, with many of them being detected in > 60% of the analyzed cells (Fig. S2A). Although several transcripts (e.g., CXCL2, MMP8, IGFBP6) were detected at a higher percentage in the C1 dataset, most transcripts were detected at higher percentage in the MI dataset (Fig. S2B). By comparison, fewer differences were seen in detection percentages between senescent and quiescent datasets, with a few exceptions (Fig. S2C). For example, for LMNB1 gene expression—which declines in senescent cells (Freund et al., 2012)—there were fewer senescent cells in which the transcript was detected compared to transcript detection in quiescent cells. Conversely, TNFRS10C gene expression, which is induced in senescent cells (Coppe et al., 2008), was detected in a higher percentage of senescent cells relative to quiescent cells (Fig. S2C). A few highly studied SASP factors, such as IL‐1B, IL‐6, IL‐7, CXCL1‐3, MMP1, MMP10, and CCL13, were also detected at low percentages (Fig. S2A), limiting the types of analyses that could be performed on these factors.

Our initial analysis revealed few significant differences between the two cell populations primarily because senescent cells had a lower overall transcript abundance, possibly owing to less efficient lysis/amplification. After normalization (Livak et al., 2013), however, it was possible to identify transcriptional profiles that were consistent with transcriptome analyses of bulk cell populations (Fig. 1A). Normalization also resulted in more coordinated expression profiles between the C1 and MI pools of cells (Fig. 1B; Fig. S3A,B). Following normalization, MI cells showed 61 significantly altered genes during senescence, while C1 showed 31, presumably due to the smaller number of successfully amplified genes. Notably, 28 of the altered genes were shared between C1 and MI cells (Fig. 1C).

Using these datasets, we first asked whether senescent and quiescent cells were distinguishable at the single‐cell level by gene expression signatures. Principal component analysis (PCA) allowed us to identify senescent cells by sorting the first and second principal components (consisting of a combination of expression values for 70 genes) (Fig. 2A, right panel). When plotted in two dimensions (X–Y axes), quiescent and senescent populations clustered independently (Fig. 2A, left panel). Senescent cells displayed a larger spread of values, consistent with more variability.

Using the expression of all 70 genes, linear discriminant analysis (LDA) correctly classified > 98% of cells as senescent or quiescent. Furthermore, expression of several genes individually strongly predicted senescence status. The two most predictive genes were TNFRSF10C and CDNK1A (Fig. 2B). TNFRSF10C encodes a secreted decoy receptor that prevents TRAIL‐induced apoptosis (Sheridan et al., 1997). This gene is induced in response to genotoxic stress in a p53‐dependent manner (Sheikh et al., 1999) and is a component of the SASP (identified by its alias TRAILR3) (Coppe et al., 2008). CDKN1A, which encodes p21^Waf1, is a cyclin‐dependent kinase inhibitor (CDKi) and primary mediator of both the transient and permanent (senescence) p53‐induced G1 arrests (El‐Deiry et al., 1993; Dulic et al., 1994). CDKN2A, which encodes p16^INK4a, is another CDKi, and common marker of senescence and aging (Beausejour et al., 2003; Krishnamurthy et al., 2004), but was not assayed in these analyses owing to primer failures. Reduced LMNB1 or HMGB1 expression also strongly predicted senescence. Both gene products are lost from the nuclei of senescent cells in a p53‐dependent manner (Freund et al., 2012; Davalos et al., 2013). As bleomycin causes genotoxic stress and activates p53, CDKN1A, TNFRSF10C, LMNB1, and HMGB1 are likely markers of p53 activation during senescence. Indeed, the combination of CDKN1A, CDKN1B, LMNB1, TNFRSF10C, and CCL3 was sufficient to predict senescence in 97% of cells (P = 9.54E‐39) by LDA (Fig. 2C). Thus, markers of p53 activity are strong indicators of senescence at the single‐cell level.

A few quiescent cells associated more strongly with senescent cells using both PCA and LDA (Fig 2A,C). As even early passage fibroblast populations often contain a small percentage of senescent cells (Dimri et al., 1995), this association was expected. Together, these analyses indicate that senescent cells can be identified at the single‐cell level based upon their gene expression signatures.

There is increased variability in selected mRNA levels among single cardiomyocytes from aged, compared to young, mice, and in mouse fibroblasts treated with H₂O₂ (Bahar et al., 2006). Both age and H₂O₂ increase the number of senescent cells. In agreement with these studies, many of the genes (51%) assayed here displayed significantly increased variability in mRNA levels among single cells following induction of senescence. This increased variability was not universal, however, with a few genes (7%) (such as CYR61 and SOCS3) showing lower variance in senescent compared to quiescent cells (Fig. 3A). Control genes, known to show little change in expression in senescent cells, generally displayed only minor changes in variance. For example, GAPDH displayed a nonsignificant (P > 0.05, Monte Carlo permutation test) decrease in variability, whereas ACTB variability increased slightly (Fig. 3A–C). Interestingly, CDKN1A also showed no significant increases in variance (P > 0.05), despite being induced in senescent cells (Fig. 3B).

Senescence‐associated variability in gene expression could result from genomic rearrangements, as suggested for aged and H₂O₂‐treated cells (Bahar et al., 2006), but also from changes in transcriptional activation or repression in senescent cells. LMNB1 mRNA levels, which decline in senescent cells (Freund et al., 2012), also showed increased variability (Fig. 3B) accompanied, as expected, by fewer transcripts with detectable amplification. Similarly, many SASP factors, such as IL‐8 and IL‐1A, showed strongly increased variability, with a subset of cells displaying a very high level of gene expression (Fig. 3B). This finding suggests the increase in mRNAs encoding SASP factors (Coppe et al., 2008) could result from contributions by only a subset of senescent cells at any time.

Are cells that strongly express SASP factors such as IL‐8 and IL‐1A a subset of senescent cells, or do individual cells express these and other senescence‐associated transcripts in variable quantities? To address these questions, we calculated correlation coefficients (R²) for all genes, eliminating nonsignificant (P > 0.01) correlations (Fig. S4). Two classes of inversely correlated genes emerged from analyses of quiescent cells (arbitrarily labeled ‘Class 1’ and ‘Class 2’) (Fig. 4A, left). From analyses of senescent cells, these classes were largely retained, but most correlation values declined (Fig. 4A, right), with many no longer meeting the P‐value threshold. Pathway analysis of Class 1 genes showed a strong relationship between extracellular matrix regulation, cellular stress responses (including hypoxic responses), and cytokine signaling, whereas Class 2 genes most strongly associated with inflammatory responses (Fig. S5).

Most of the declines in correlation were likely the result of increased variability. However, a subset of genes showed increased correlations (ΔR² > 0.5; R² (sen) > R² (qui)), either directly or inversely (Fig. 4A). Among these genes, there was no relationship between variability in transcript abundance and correlation among paired transcripts (Fig. 4B), suggesting that variability alone cannot account for these increases.

Of all genes that displayed highly altered correlations with senescence, CDKN1A (which was consistently induced in senescent cells; Fig. 3B) was most notable, displaying increased correlations with 25 gene transcripts (Fig. 4C). Furthermore, CDKN1A showed a substantial shift in its correlation patterns, losing correlation with some genes (Class 1) and gaining correlation with others (Class 2) (Fig. 4A).

As many SASP factors are strongly clustered in the genome (Coppe et al., 2010), we asked whether genomic location is linked to expression coregulation. Supporting this possibility, transcripts from the IL‐1 cluster (IL‐1A and IL‐1B, promoters separated by ~50 kb), as well as the CXCL cluster (in order: IL‐8, CXCL1, CXCL5, CXCL3, and CXCL2, separated in total by ~360 kb), went from nonsignificant correlations to stronger, significant direct correlations, suggesting these genes were induced in a coordinated manner (Fig. 4D). By comparison, little to no correlation of expression was observed when the IL‐1 cluster was analyzed against the CXCL cluster (Fig. 4D), which are on different chromosomes. These data suggest that genomic organization can influence gene expression changes in single cells.

Together, our correlation data indicate that, whereas the expression of many genes is coordinated under quiescent conditions, some senescence‐specific gene expression processes appear to be regulated independently of each other.

IMR‐90 human fetal lung fibroblasts at population doubling level 25–30 were cultured in Dulbecco's modified Eagle's medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and penicillin/streptomycin (growth medium). Cells were seeded at 1 x 10⁴ cells/cm² in 175‐cm² flasks, maintained at 37° C, 10% CO_2, and 3% O₂, and were mycoplasma negative. For senescence induction, cells were given growth medium containing 50 μm bleomycin (senescent) or vehicle (PBS—nonsenescent control) in atmospheric oxygen for 3 hrs. The cells were then washed, given growth medium, and returned to 3% O₂. The growth medium was changed every 3 days. On day 7, it was replaced with low serum (0.2% FBS) medium, which induced quiescence in control cells. Cells were harvested by trypsinization for analysis on Day 10 (see Fig. S1A). A subset of cells was plated for senescence assays (below).

Ten days after bleomycin or control treatment, a 12‐well plate was seeded with three wells from each treatment (30 000 cells/well) and allowed to attach overnight. The assay was performed as described (Dimri et al., 1995) using a commercial kit (BioVision, Milpitas, CA, USA).

Quiescent or senescent cells were seeded as for the SA‐Bgal assay and allowed to attach overnight. The following day, the medium was replaced with 500 μL of serum‐free DMEM. Conditioned media were collected 24 hrs later. An ELISA kit (AlphaLISA, Perkin Elmer Biosciences, San Jose, CA, USA) was used to measure IL‐6 in conditioned media according to the manufacturer's instructions; IL‐6 levels were normalized to cell number.

Quiescent or senescent cells were plated into four‐well chamber slides (Labtech, Tampa, FL, USA) and allowed to attach overnight. The next day, the cells were given low serum medium containing 10 μm EdU (Invitrogen, Carlsbad, CA, USA). After 24 hrs, cells were washed with PBS, fixed in 4% formalin, then permeabilized in 0.5% Triton X‐100 for 30 min, and washed twice in ice‐cold PBS. Permeabilized cells were subjected to the EdU reaction according to the manufacturer's instructions (Click‐IT EdU Imaging Kit, Invitrogen).

Cells cultured as described above were trypsinized, pelleted by centrifugation, and resuspended in Krebs–Henseleit buffer. Approximately 20 000 cells were placed in 10 mL of buffer in a 15‐cm dish. Under phase‐contrast microscopy, a small glass needle was used to mouth pipette individual cells into 200‐μL reaction tubes in less than 2 μL volumes. Single cells were immediately frozen on dry ice and maintained at −80° C until analysis.

Single cells were captured in the C1 platform (Fluidigm Inc., South San Francisco, CA, USA) as per the manufacturer's instructions (Fluidigm) and verified as single‐cell captures by microscopic analysis of each nanofluidic chamber. Chambers on the C1 chip containing more than one cell were excluded from subsequent analysis.

Gene expression in single cells was measured by quantitative multiplex PCR following pre‐amplification using the Biomark platform (Fluidigm) according to the manufacturer's instructions. Data were normalized as described (Livak et al., 2013). Briefly, the data adjustment was as follows: Log2 expression data were extracted through the single‐cell analysis package for all PCR plates, including both C1 and manual isolation data. The data were then categorized as deriving from quiescent or senescent cells prior to subsequent analyses.

PCA was performed on qPCR values of 70 genes from 464 individual cells, 235 quiescent cells, and 229 senescent cells (one senescent cell was removed because none of its transcript was detected by qPCR) using the pcaMethods package in R (probability principle component analysis method with missing values imputed by default).

LDA was performed to maximize the separation between the 235 quiescent and 230 senescent cells. To determine which gene subsets were most important in correctly classifying cells as quiescent or senescent, we used a sequential backward selection process (Cotter et al., 2001). Starting with all the data, we used LDA to obtain the maximum separation between quiescent and senescent cells. Then, the loss of a single transcription vector was tested by removing each gene and repeating the LDA analysis with the data from the remaining 70 genes. The gene loss having the smallest impact on the ability to separate quiescent from senescent cells was deleted. This process was repeated until only the most predictive genes for discrimination remained. Because imputation of missing expression values can influence the results of LDA and the backward selection process, three separate imputation schemes were employed, with the final ranking taking each of these gene rankings into account. Missing values were first set to zero, then to the average of all values recorded for that gene, and finally to the average of all quiescent or senescent values.

Class 1 and Class 2 gene lists were submitted to the Database for Annotation, Visualization and Integrated Discovery (DAVID v.6.7) (Huang et al., 2007) for annotation enrichment analysis. The full human genome was used as a background. Pathways and annotations showing significant enrichment following Benjamini multiple testing correction are summarized in Fig. S4: