Research Paper
Single-Cell Transcriptomics Reveals Global Markers of Transcriptional Diversity across Different Forms of Cellular Senescence
Authors
Shane A. Evans,1 Yee Voan Teo,2 Samuel J. Hinthorn,1 Kelly Clark,1 Takahiro Ito,2 John M. Sedivy,2 and Nicola Neretti1,2,*
1Center for Computational Molecular Biology, Brown University, Providence, RI, USA
2Department of Molecular Biology, Cell Biology and Biochemistry, Brown University, Providence, RI, USA
*Corresponding author: nicola_neretti@brown.edu
DOI:https://doi.org/10.59368/agingbio.20230008
Full Article | PDF | Supplementary
Abstract
Cellular senescence (CS) is a state of irreversible cell cycle arrest, and the accumulation of senescent cells contributes to age-associated organismal decline. The detrimental effects of CS are due to the senescence-associated secretory phenotype (SASP), an array of signaling molecules and growth factors secreted by senescent cells that contribute to the sterile inflammation associated with aging tissues. Recent studies, both in vivo and in vitro, have highlighted the heterogeneous nature of the senescence phenotype. Single-cell transcriptomics has revealed that oncogene-induced senescence (OIS) is characterized by the presence of subpopulations of cells expressing different SASP profiles. We have generated a comprehensive dataset via single-cell transcriptional profiling of genetically homogenous clonal cell lines from different forms of senescence, including OIS, replicative senescence, and DNA damage-induced senescence. We identified subpopulations of cells that are common to all three major forms of senescence and show that the expression profiles of these subpopulations are driven by markers formerly identified in individual forms of senescence. These common signatures are characterized by chromatin modifiers, inflammation, extracellular matrix remodeling, and ribosomal protein gene expression (measured at the RNA level). The expression patterns of these subpopulations recapitulate primary and juxtacrine secondary senescence, a phenomenon where a pre-existing (primary) senescent cell induces senescence in a neighboring (secondary) cell through cell-to-cell contact. Hence, our results demonstrate that the formation of juxtacrine secondary populations of cells is common to multiple types of senescence and occurs in competition with primary senescence. Finally, we show that these subpopulations show differential susceptibility to the senolytic agent Navitoclax, suggesting that senolytic agents targeting the apoptotic pathways may be clearing only a subset of senescent cells based on their inflammatory profiles.
Introduction
Cellular senescence (CS) is a programmed stress response that leads to a cell’s permanent exit from the cell cycle and can be induced by a variety of factors including telomere attrition, oncogene activation, oxidative stress, and DNA damaging agents1–3. Although CS comes in different forms, an established senescent pathway involves a persistent DNA-damage response that leads to the activation of the tumor suppressor protein 53 (P53), which, in turn, activates the cyclin-dependent kinase inhibitor gene CDKN1A. The translated protein encoded by the CDKN1A gene, p21, holds the cell in cell cycle arrest until upregulation of the cyclin-dependent kinase inhibitor gene CDKN2A, which encodes the p16 protein that maintains the cell in an irreversible senescent state1,4. During senescence, the cell undergoes global epigenetic changes including dramatic chromatin alterations and increased expression of an array of extracellular remodeling proteins, growth factors, and inflammatory molecules such as interleukins and interferons that compose the senescence-associated secretory phenotype (SASP)4–7. SASP leads to inflammation disrupting the tissue microenvironment, reinforces the senescent phenotype by contributing to the cell cycle arrest, and can induce paracrine senescence in normal cells8,9.
Studies describing the heterogeneity of senescent cells have shown that there are different forms of SASP. For example, one form is dominated by transforming growth factor beta (TGF-beta) signaling10–12. TGF-beta-dominated profiles are characterized by extracellular matrix remodeling, collagen deposition, and by the expression of growth factors such as the connective tissue growth factor (CTGF)10,13,14. Another form of SASP, however, shows a more pro-inflammatory profile with higher levels of interleukins such as IL1A, IL1B, IL6, and other genes regulated by Nuclear factor kappa-light-chain-enhancer of activated B cells (NFKB)10–12,15. In oncogene-induced senescence (OIS), the prevalence of these profile changes as cells persist in the senescent state: TGF-beta signaling is typically higher in the earlier stages of senescence, whereas in the later stages senescent cells transition to a pro-inflammatory phenotype10.
Not only do transcriptional profiles change over the course of time, but they can also characterize different types of senescent states. For instance, TGF-beta-dominated SASP profiles are prevalent in notch-induced senescence (NIS) 10–12, which occurs when a primary senescent cell (or a pre-existing senescent cell) makes direct contact and activates the notch signaling pathway in a neighbor cell, causing the spread of senescence through juxtacrine signaling and thereby creating a secondary senescent cell. OIS cells can act as the primary senescence source, and display more pro-inflammatory SASP profiles with higher levels of interleukins and NFKB regulated genes in contrasts with the anti-inflammatory TGF-beta-dominated profiles of NIS cells11.
CS is an inherently heterogeneous state, as it can be cell type- and insult-dependent1,4,16. However, most of the data describing senescence have been collected using bulk sequencing technologies that measure average gene expression across large heterogeneous pools of cells and is blind to cell-to-cell transcriptional variability. With the advent of single-cell transcriptomics, it has become clear that distinct subpopulations contribute to the population-wide average11,17. Despite the limited single-cell transcriptional data currently available for senescent cells, it is becoming increasingly clear that senescent cells are subject to significant transcriptional diversity10,16–18. Consistent with this transcriptional heterogeneity, subpopulations of senescent cells have been observed in single-cell transcriptomic studies from fibroblasts and endothelial cell lines11,17,18. For example, using single-cell RNA-seq, Teo et al showed that two subpopulations coexist in the OIS cultures11. One population possessed the familiar OIS transcriptional profile (primary senescence), whereas the second population, which was also senescent, possessed SASP profiles that were growth factor rich and dominated by TGF-beta signaling, and was composed of NIS cells (secondary senescence).
The senescent phenotype is characterized by chromatin modifications, DNA-damage response pathways, and SASP profiles1,3. However, an in-depth description of how these different aspects vary across senescent cell populations is still lacking17,19. Here, we used single-cell transcriptomics to identify subpopulation of senescent cells in genetically homogenous clonal cell lines with a inducible HRAS:G12V transgene that was activated only in OIS cultures. This strategy ensures that the subpopulations we observe are not a caused by variability of the number and genomic location of the transgene constructs, and the consequent variability in the expression levels of the transgene across cells. This experimental design allowed us to conduct a novel, unified analysis of senescent cell heterogeneity across the major forms of senescence.
Materials and Methods
Single-cell RNA sequencing
Human diploid fibroblast cells were trypsinized and centrifuged at 500 rcf for 10 minutes. Cells were resuspended in cold phosphate-buffered saline and passed through a 40 μM Flowmi Cell Strainer. Cells were then counted and loaded onto the 10× chromium using V2 chemistry of 10× genomics’ 3-prime single-cell reagents with version 2 chemistry. Libraries were prepared according to the manufacturer protocol and sequenced on a Hi-Seq platform at GeneWiz with manufacture recommend sequencing specifications. Senescent libraries were subjected to two rounds of sequencing. All single-cell RNA-seq datasets generated in this work have been deposited in the Gene Expression Omnibus database with accession number GSE173879.
scRNA-seq data processing and filtering
Cell-specific barcodes were error corrected and identified from fastq files, and data were aligned to the hg19 reference genome using CellRanger V2.1 Command Line tools. For measuring expression of neo selectable marker, a second round of alignment was conducted to the hg19 reference genome along with the neomycin sequence from pLNCX2. Secondary analysis was conducted using Seurat R package version 319,20. Data generated from each cell population (growing, replicative senescence [RS], OIS, and DNA damage-induced senescence [DDIS] cells) were filtered separately. Stringent filtering methods were applied using parameters described in the literature such as the number of genes detected, number of unique transcripts detected, percent mitochondrial genes detected, and percent ribosomal protein RNA detected19,20. To filter cells, we refrained from regressing out the effects of the abovementioned quality control metrics when implementing clustering protocols. Then, we were able to cluster the majority of “low-quality” cells separately from the good-quality cells. Low-quality clusters were then discarded.
Clustering scRNA-seq data
Filtered datasets were merged together, scaled, and re-normalized using Seurat19,20. In this case, the number of genes detected, number of unique transcripts detected, and percent of mitochondrial genes detected were all scaled out using Seurat’s ScaleData method with the var.to.regress parameter, thereby ensuring that any conclusions drawn from downstream analysis were driven by variable gene expression and not technical factors such as cell-to-cell variability in sequencing depth. Next, principal component analysis (PCA) was conducted, and, depending on the dataset, the top 20–30 principal components were fed into Seurat’s shared nearest neighbor-based Louvain clustering algorithm.
Differential expression analysis
Differential expression was performed using Seurat with model-based analysis of single cell transcriptomics (MAST) methodology21,22. Cluster 1 was compared to cluster 2 in each type of senescence.
Ingenuity pathway analysis (IPA)
The computed log fold changes of significantly differentially expressed genes (false discovery rate less than 0.05) were used in IPA by Qiagen. Upstream regulators were plotted according to their z score and q value of their overlap using GGplot2 R package. Q values for upstream regulators were computed using p.adjust function in base R.
Gene set enrichment analysis (GSEA)
Computed log fold changes of genes differentially expressed between clusters 1 and 2 were used in a pre-ranked GSEA analysis. Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were downloaded from MSigDB23. GO terms, included custom SASP lists, were based on the existing literature24–30. KEGG pathways for primary and secondary senescence were generated by Teo et al and were added to the list of KEGG pathways11. Bar charts were generated using ggplot2, and tables were generated in Microsoft Excel.
Navitoclax experiments and analysis
Senescent cells were treated with 1 μM of Navitoclax for three days. For OIS and DDIS samples, 4OHT and Etoposide were removed from the growth media before Navitoclax treatment. Dimethyl sulfoxide (DMSO)-treated cells served as controls. Cells recovered in regular growth media supplemented with DMSO overnight before being harvested for scRNA-seq. Clusters were identified in DMSO-treated controls. p values were calculated using the chisq.test function in R.
Retroviral infections of ERT2-HRAS:G12V
The 293T cells were co-transfected with plasmid DNAs of a retroviral vector and the helper vectors using FuGENE HD (Promega). Medium was collected 24, 36, and 48 hours later for infection of LF1 cells. The vector backbone was clonetech pQCXIN (Addgene Cat Number 631514). Clonal cell lines were generated through serial dilution in 500 μg/mL of G418. Colonies generated from single cells were selected and further propagated.
Cell culture and senescence induction
We sequenced proliferating female human diploid fibroblast cells (LF1 cells) along with populations of LF1 cells that were induced into RS, OIS, and DDIS. LF1 cells were obtained from the Sedivy lab. All cell populations were generated with a clonal cell line that possessed an ERT2-HRAS:G12V transgene. For proliferating, populations cells were passaged in regular growth media and harvested at 60%–80% confluency. Meanwhile, replicative senescent cells were passaged in regular growth media until replicative exhaustion. DDIS was induced with addition of etoposide at 40 μM for three weeks to regular growth media. OIS was induced by adding 4-Hydroxytamoxifen (4-OHT) to regular growth media, at which point the cells underwent a hyperproliferative phase before senescing after six days. Regular growth media consisted of Ham’s F10, 15% Fetal Bovine Serum (FBS), and 1× of penicillin streptomycin and glutamine. Cells were maintained at 37°C at 5% CO2 and 2.5% O2.
Results
Single-cell transcriptional profiling of clonal lines
We performed single-cell RNA sequencing using the 10× chromium microfluidics platform to study cell-to-cell gene expression heterogeneity in RS, OIS, and DDIS. All cell populations were generated using female human diploid fibroblast cells (LF1 cells) and were derived from a clonal cell line that possessed a 4-OHT inducible HRAS:G12V transgene ( Fig. 1A). Although subpopulations of senescent cells had been identified in OIS cultures, little is known about the cell-to-cell diversity in other forms of senescence. We hypothesized that the diversity found in OIS cultures is a phenomenon common across multiple forms of senescence. To test this hypothesis, we sequenced different types of senescence in different cell culture conditions. For full list of culture conditions sequenced for this study, see Table 1. For results in the main figures in this paper, we focus on the three of the datasets comprising RS, OIS, and DDIS cells. Analysis of the remaining datasets, which revealed strikingly similar diversity, is shown in the supplemental figures. Specifically, for OIS cultures we have samples that were in 4-OHT, whereas others were removed from 4-OHT prior to sample collection. This helped us to determine how 4-OHT affected diversity. For RS, we sequenced two clonal cells lines, referred to as RS and RS Rep2. This let us show that the observed heterogeneity was not due to the use of a specific clonal cell line but instead could be recapitulated across multiple clonal lines. Moreover, for OIS samples we also performed 5′ end sequencing. This let us show that the heterogeneity we discovered in our datasets is not due to a 3′ sequencing bias.
(A) Schematic showing generation of clonal cell lines with HRAS:G12V transgene and overall experimental design. (B) Uniform Manifold Approximation and Projection (UMAP) plots showing all datasets generated for this study (left) and their predicted stage in the cell cycle (right). (C) Violin plots showing downregulation of cell cycle genes (used to generate G2 and S scores) and chromatin modifiers HMGB and HMGB2 in senescent cells. Meanwhile, there is an upregulation of CDKN2A and CDKN1A in senescent cultures.
| Dataset Name | Condition | 5′ or 3′ | Clonal Line | Used in Main Text or Supplementary Figures |
|---|---|---|---|---|
| DDIS | Etoposide treated for three weeks. | 3′ | 1 | Main text |
| Etoposide was removed and replaced with DMSO for the last three days | ||||
| OIS | 4-OHT treated for three weeks. | 3′ | 1 | Main text |
| 4-OHT was removed and replaced with DMSO for the last three days. | ||||
| RS | Regular growth media | 3′ | 1 | Main text |
| RS Rep-2 | Regular growth media | 3′ | 2 | Supplementary |
| OIS + Nav | 4-OHT treated for three weeks. | 3′ | 1 | Supplementary |
| 4-OHT was removed and replaced with Navitoclax dissolved in DMSO for the last three days. | ||||
| OIS 5’ | 4-OHT treated for three weeks. | 5′ | 1 | Supplementary |
| 4-OHT was removed and replaced with DMSO for the last three days. | ||||
| OIS + Nav 5’ | 4-OHT treated for three weeks. | 5′ | 1 | Supplementary |
| 4-OHT was removed and replaced with Navitoclax dissolved in DMSO for the last three days. | ||||
| DDIS + Nav | Etoposide treated for three weeks. | 3′ | 1 | Supplementary |
| Etoposide was removed and replaced with Navitoclax dissolved in DMSO for the last three days. | ||||
| OIS in 4-OHT | 4-OHT treated for one month. 4-OHT was in the media up until cells were harvested. | 3′ | 1 | Supplementary |
| Growing | Regular growth media | 3′ | 1 and 2 | Main text |
| 4-OHT, 4-Hydroxytamoxifen; DDIS, DNA damage-induced senescence; OIS, oncogene-induced senescence; RS, replicative senescence. | ||||
Because all of our datasets were generated with cells that were infected with an HRAS:G12V transgene, controlling for variable expression of this transgene is important if we are to draw conclusions regarding the formation of subpopulations in senescence. Without clonal cell lines, the transgene would be expressed at varying levels across the cells, and because HRAS is upstream of many important molecular pathways that are heavily implicated in senescence, then this would have been a confounding factor in our experiments. Our approach contrasts previous studies which conducted scRNA-seq on non-clonal cell lines that variably express the transgene11. Moreover, working with clonal cell lines ensured that any variability in gene expression that we observed was a product of the senescent phenotype and not due to a pre-existing heterogeneity in the proliferating cell populations.
For all types of senescence, we noticed the presence of a large populations of cells that showed signs of low-quality data. These included higher than average expression of ribosomal protein genes in combination with the low number of genes being expressed and low number of unique molecular identifiers (UMIs). In addition, we saw a different and smaller subpopulation of cells with the high number of mitochondrial reads. Our proliferating control cells did not possess these subpopulations in large numbers. Taken together, these observations suggest that senescent cells are more fragile to the microfluidics used to generate these libraries. For all datasets, we filtered out these low-quality droplets and the remainder of our analysis focused on the higher-quality cells which we retained. See Materials and Methods section and Supplemental Figure 7, for a detailed explanation of our filtering strategy. RS, OIS, and DDIS cultures showed clear senescent gene expression patterns including upregulation of CDKN2A and downregulation of cell cycle genes, HMGB1, and HMGB2 (Figs. 1C, 2H,J).
(A) Schematic showing the process for analyzing each senescent dataset separately. Oncogene-induced senescence (OIS) data were chosen for this illustration, but the same process was used for all datasets. (B–D) UMAP plots showing DNA damage-induced senescence (DDIS), OIS, and replicative senescence (RS) clusters. (E–G) The corresponding cluster trees. (H–J) Violin plots showing similar levels of senescence markers for each senescent subpopulation (ns: p > 0.05; *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ****p ≤ 0.0001).
All datasets were merged using the Seurat Bioconductor Package19,20. Our filtered dataset included a total of 6108 cells split across multiple types of senescence and conditions. We projected individual cells’ transcriptional profiles onto two dimensions using the Uniform Manifold Approximation and Projection reduction and observed a clear separation between senescent and growing cells. Computing a cell cycle score with Seurat shows that the proliferating population is composed of cells in the S, G2/Mitosis, and G1/G0 phases of the cell cycle. We noticed that the different types of senescent cells are positioned more closely to control cells in the G1/G0 phase, showing that the senescent state as a distinct and final cell cycle phase (Fig. 1B).
Clustering reveals two subpopulations in clonal cell lines
We combined each senescent dataset along with an equal-sized subpopulation of control cells (Fig. 2A). Clustering was performed on the resulting datasets (see Materials and Methods section), and in all cases, we could identify subpopulations of senescent cells. For each type of senescence, we call the subpopulations “-1” and “-2” (Fig. 2B–G). We analyzed the expression of the transgene between the two clusters by aligning to the amino 3′-glycosyl phosphotransferase (neo) selectable marker gene. We saw that the transgene is much more evenly expressed between the subpopulations we observe in cultures generated from clonal cell lines when compared to previous published experiments which did not use clonal lines (Supplemental Fig. 6). In addition, when combined, the cluster 2 cells tend to be more similar to each other across the three forms of senescence than to cells in cluster 1 (Supplemental Fig. 9).
Subpopulations are distinguished by TGF-beta signaling, DNA-damage response, and inflammatory pathways
We performed a differential expression analysis between clusters 1 and 2 in each type of senescence using the Seurat package with MAST methodology21,22. The computed log fold changes of genes passing a false discovery rate of 0.05 were used to identify potential upstream regulators via the IPA (Fig. 3). We plotted the upstream regulators according to two dimensions. On the x-axis are the computed z scores, which show the direction of regulation. On the y-axis is the negative logarithmic transformation of their false discovery rate, which determines the upstream regulators whose gene sets show a statistically significant overlap with the list of differentially expressed genes in the data.
(A) Violin plots showing expression of marker genes for primary and secondary senescence for clusters 1 and 2 in DDIS, OIS, and RS cells. HMGA1 is known to be higher in primary OIS (Teo et al, 2019), and the remaining genes are known to be higher in secondary notch-induced senescence. (B) Volcano plots showing predicted upstream regulators (Ingenuity Pathway Analysis) for cluster 1 (red) and cluster 2 (blue) for each type of senescence. Regulators are plotted according to their z score on the x-axis, which shows if they regulate cluster 1 or cluster 2. Positive values indicate an upstream regulator for cluster 1 and negative for cluster 2. The “−log q value of overlap” is plotted on the y-axis. This is the negative logarithmic transformation of the false discovery rate, which determines the upstream regulators whose gene sets show a statistically significant overlap with the list of differentially expressed genes in the data. Gray lines indicate significance values. Upstream regulators located above the horizontal line and to the left of the vertical left gray line or right of the vertical right gray line have significant z scores and false discovery rates. (C) Heatmap showing z scores for the upstream regulators for each type of senescence. Positive values indicate an upstream regulator for cluster 1 and negative for cluster 2 (ns: p > 0.05; *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ****p ≤ 0.0001).
Using a z score cutoff 1.64 and −log q-value cutoff of 1.301, we were able to identify statistically significant upstream regulators. We noticed that the predicted upstream regulators were very similar for each type of senescence. More specifically, cluster 1 showed upstream regulators TP53, a marker for DNA damage13,31–33. We also saw Interferon term IFNA2, suggesting a more pro-inflammatory profile30. Meanwhile, cluster 2 showed upstream regulators corresponding to TGF-beta signaling and extracellular matrix remodeling. These regulators included TGF-beta receptors, TGF-beta1–3, and several Suppressor of Mothers against Decapentaplegic (SMAD) proteins that are transducers for TGF-beta signaling. We also performed a comparison analysis for each type of senescence and noticed a strong similarity in the z scores for the upstream regulators between all types of senescence that we sequenced (Fig. 3). We also confirmed the existence of these subpopulations in all other senescent datasets (see Supplemental Fig. 1, for senescent markers, supplemental tables provided in Excel sheet for GSEA and KEGG terms, and Supplemental Fig. 8 for IPA upstream regulators).
Next, we took the log fold changes between clusters 1 and 2 for each type of senescence and conducted a pre-ranked GSEA with GO terms and KEGG pathways24–29. We plotted the normalized enrichment scores for statistically significant (p < 0.05) GO terms and KEGG pathways. We noticed that subpopulations in cluster 1 showed enrichment for terms related to DNA-damage response pathways such as telomere organization, mismatch repair, base excision repair, and DNA-ligation. Moreover, there was also an enrichment for SASP and inflammation pathways including our custom SASP lists (see Materials and Methods section, for details on custom SASP lists), inflammatory response, cytokine interactions, viral life cycle, neutrophil migration, and viral transcription (Tables 2–4).
| GSEA/KEGG | NES | p value |
|---|---|---|
| viral_life_cycle | 2 | 0 |
| Up_primary_secondary_(KEGG) | 1.9 | 0 |
| intrinsic_apoptotic_signaling_by_DNA_damage | 1.8 | 0 |
| GO_SHELLY_NATURE_SASP | 1.7 | 0.006 |
| p53_signaling | 1.7 | 0 |
| GO_CAMPISI_V2 | 1.6 | 0.027 |
| chemokine_signaling | 1.6 | 0.024 |
| positive_inflammatory_response | 1.4 | 0.049 |
| apoptotic_signaling_pathway | 1.3 | 0.018 |
| cytokine_receptor_interaction_(KEGG) | 1.3 | 0.032 |
| regulation_of_actin_cytoskeleton_(KEGG) | −1.6 | 0 |
| cell_adhesion_molecules_(KEGG) | −1.6 | 0.004 |
| actin_filament_based_process | −1.7 | 0 |
| actin_filament_organization | −1.7 | 0 |
| cell_junction_assembly | −1.7 | 0 |
| negative_regulation_of_cytokine | −1.7 | 0.002 |
| positive_wound_healing | −1.7 | 0 |
| regulation_of_wound_healing | −1.7 | 0 |
| response_to_wounding | −1.7 | 0 |
| wound_healing | −1.7 | 0 |
| actomyosin_structure_organization | −1.8 | 0 |
| cell_junction_organization | −1.8 | 0 |
| cell_matrix_adhesion | −1.8 | 0 |
| cell_substrate_adhesion | −1.8 | 0 |
| substrate_adhesion_cell_spreading | −1.8 | 0 |
| collagen_fibril_organization | −1.9 | 0 |
| extracellular_structure_organization | −1.9 | 0 |
| ECM_receptor_interaction_(KEGG) | −1.9 | 0 |
| Up_secondary_primary_(KEGG) | −2.2 | 0 |
| GSEA results for cluster 1 versus cluster 2 for OIS samples are shown. Positive normalized enrichment scores (NESs) are pathways/gene sets enriched in cluster 1, and negative NES are enriched in cluster 2. | ||
| GSEA/KEGG | NES | p value |
|---|---|---|
| viral_life_cycle | 2 | 0 |
| IL1_secretion | 1.7 | 0.012 |
| interferon_response | 1.6 | 0.02 |
| NIK_NF_KappaB_signaling | 1.6 | 0 |
| IL1B_production | 1.5 | 0.021 |
| IL1_production | 1.5 | 0.031 |
| innate_immune_activation | 1.4 | 0.011 |
| DNA_damage_signaling | 1.4 | 0.037 |
| KANNAN_TP53_TARGETS_UP | 1.4 | 0.019 |
| positive_innate_immune_response | 1.3 | 0.026 |
| p53_signaling | 1.3 | 0.021 |
| TNF_signaling | 1.3 | 0.028 |
| Up_primary_secondary_(KEGG) | −1 | 0.389 |
| NK_cell_mediated_cytotoxicity_(KEGG) | −1.3 | 0.039 |
| response_to_TGFbeta | −1.5 | 0 |
| LABBE_WNT3A_TARGETS_UP | −1.5 | 0 |
| TGFbeta_production | −1.6 | 0.003 |
| response_to_wounding | −1.6 | 0 |
| actin_filament_based_movement | −1.7 | 0 |
| actin_filament_organization | −1.7 | 0 |
| negative_regulation_of_cytokine | −1.7 | 0.001 |
| regulation_of_wound_healing | −1.7 | 0 |
| substrate_adhesion_cell_spreading | −1.7 | 0 |
| wound_healing | −1.7 | 0 |
| regulation_of_actin_cytoskeleton_(KEGG) | −1.7 | 0 |
| actin_filament_based_process | −1.8 | 0 |
| cell_junction_assembly | −1.8 | 0 |
| cell_junction_organization | −1.8 | 0 |
| cell_substrate_adhesion | −1.8 | 0 |
| cell_adhesion_molecules_(KEGG) | −1.8 | 0 |
| Up_secondary_primary_(KEGG) | −1.8 | 0 |
| extracellular_structure_organization | −1.9 | 0 |
| ECM_receptor_interaction_(KEGG) | −1.9 | 0 |
| GSEA results for cluster 1 versus cluster 2 for RS samples are shown. Positive normalized enrichment scores (NESs) are pathways/gene sets enriched in cluster 1, and negative NES are enriched in cluster 2. | ||
| GSEA/KEGG | NES | p value |
|---|---|---|
| viral_life_cycle | 2 | 0 |
| intrinsic_apoptotic_pathway | 1.9 | 0 |
| p53_signaling | 1.9 | 0 |
| DNA_damage_signaling | 1.9 | 0 |
| cellular_response_to_UV | 1.8 | 0 |
| intrinsic_apoptotic_by_DNA_damage_and_p53 | 1.8 | 0.006 |
| intrinsic_apoptotic_pathway_by_p53 | 1.7 | 0.003 |
| response_to_UV | 1.6 | 0.009 |
| chemokine_signaling | 1.5 | 0.031 |
| stat_cascade | 1.5 | 0.011 |
| SMAD_signaling | 1.5 | 0.045 |
| apoptotic_signaling_pathway | 1.4 | 0 |
| Up_primary_secondary_(KEGG) | 1.3 | 0 |
| TGF_beta_signaling_pathway_(KEGG) | −1.5 | 0.004 |
| positive_substrate_adhesion_cell_spreading | −1.6 | 0.003 |
| LABBE_WNT3A_TARGETS_UP | −1.6 | 0 |
| cell_matrix_adhesion | −1.7 | 0 |
| negative_regulation_of_cytokine | −1.7 | 0.001 |
| non_canonical_Wnt_signaling | −1.7 | 0 |
| positive_regulation_of_wound_healing | −1.7 | 0 |
| response_to_TGFbeta | −1.7 | 0 |
| cell_adhesion_molecules_(KEGG) | −1.7 | 0.001 |
| actin_filament_based_process | −1.8 | 0 |
| actin_filament_organization | −1.8 | 0 |
| cell_junction_assembly | −1.8 | 0 |
| cell_junction_organization | −1.8 | 0 |
| response_to_wounding | −1.8 | 0 |
| actomyosin_structure_organization | −1.9 | 0 |
| cell_substrate_adhesion | −1.9 | 0 |
| regulation_of_actin_cytoskeleton_(KEGG) | −1.9 | 0 |
| extracellular_structure_organization | −2.1 | 0 |
| Up_secondary_primary_(KEGG) | −2.1 | 0 |
| GSEA results for cluster 1 versus cluster 2 for DDIS samples are shown. Positive normalized enrichment scores (NESs) are pathways/gene sets enriched in cluster 1, and negative NES are enriched in cluster 2. | ||
In contrast, cells in cluster 2 showed higher levels of extra cellular matrix activity including terms related to integrin signaling pathways, extracellular organization, cell adhesion, collagen organization, Major Histocompatibility Complex (MHC) class 1 antigen presentation, and cell junction organization. Cluster 2 was also enriched for terms related to TGF-beta signaling such as regulation of TGF-beta, SMAD protein signaling, Wingless-Type (WNT) signaling, and anti-inflammatory gene activity (Tables 2–4). Collagen deposition, extra cellular matrix remodeling, and WNT signaling are all regulated or co-activated by TGF-beta signaling13,14,31,34. Therefore, this is consistent with the expression profile of cluster 2 cells as being dominated by TGF-beta signaling pathways. These GO terms and KEGG pathways were also verified in all other datasets we collected (Supplemental Tables provided in Excel sheet).
Subpopulations resemble primary and secondary senescent cells
We investigated whether there were any differences in the expression levels of important chromatin modifiers. We found that High Mobility Group AT-Hook 1 (HMGA1) is consistently upregulated in cluster 1 (Fig. 3A). The HMGA1 gene encodes a highly abundant chromatin-associated protein that has been shown to organize chromatin architecture in senescence and is critical for the onset of OIS. Moreover, previous studies have shown that HMGA1 is expressed higher in OIS when compared with NIS, which is characterized by increased levels of TGF-beta signaling. This difference in HMGA1 expression is known to be responsible for many of the differences in chromatin architecture that exist between OIS and NIS cells, and many of these differences are predictive of gene expression10–12. These studies comparing NIS and OIS, along with the fact that we see TGF-beta signaling and elevated HMGA1 expression arising from two distinct subpopulations in our scRNA-seq data, suggests that this HMGA1 expression is elevated only in the senescent cells where TGF-beta expression is low and levels of DNA damage are high, which we would expect from a primary senescent cell. This suggests that cells induced into RS, OIS, and DDIS develop subpopulations that resemble primary and secondary senescence, which before was only demonstrated in OIS.
Because HMGA1 is differentially expressed between OIS and NIS cells, we decided to compare our data with the single-cell transcriptional profiles from Teo et al11, as they showed that OIS cultures generated from non-clonal cell lines are composed of primary and secondary senescent subpopulations. For all four datasets (RS, OIS, DDIS, and Teo et al), we see differential expression of collagen genes, TGF-beta response genes, and HMGA1 which are markers for primary and secondary senescence (see Fig. 3A, for our data, and Supplemental Fig. 4, for Teo et al data). Moreover, the expression of HMGA1 is anti-correlated with TGF-beta signaling. We then ran a KEGG analysis with terms generated by Teo et al relating to primary and secondary senescent phenotypes, we see that our subpopulations are enriched for these terms, referred to as “UP_primary” and “UP_secondary,” respectively (see Table 1, for main text data, and Supplemental Tables provided in Excel sheet for all other datasets). Like with our own data, we conducted a GSEA analysis for the data generated by Teo et al and identified very similar GO terms and KEGG pathways. Cluster 1 is enriched for SASP, inflammation, and DNA-damage response pathways. Meanwhile, cluster 2 is enriched for extracellular matrix remodeling, collagen deposition, and WNT signaling (see Supplemental Tables provided in Excel sheet).
Subpopulations show differential sensitivity to the senolytic agent Navitoclax
We wanted to determine if we could characterize the subpopulation of cells that are more resistant to treatment with BH3 mimetics. We treated DDIS and OIS cells with 1 μM of the BH3 mimetic Navitoclax for three days35. Cells were then harvested and run on the 10× chromium platform to generate single-cell libraries along with DMSO-treated controls; see Materials and Methods section and Figure 4A, for Navitoclax experimental details. Data were aligned and filtered as previously described. We identified the clusters 1 and 2 in Navitoclax-treated samples as well (Fig. 4B–G). We observed that Navitoclax preferentially induced apoptosis in DDIS and OIS cells in cluster 2, the subpopulation enriched for a secondary senescent phenotype, TGF-beta signaling, extracellular matrix remodeling, and lower levels of DNA damage (Fig. 5). Moreover, we also saw this phenomenon in OIS cells that were sequenced from the 5′ end, although this dataset did not reach statistical significance (Supplemental Fig. 3). This suggests that TGF-beta signaling in senescence may sensitize cells to drug-induced apoptosis, emphasizing the translational importance of these subpopulations. In general, these data demonstrate that for the development on senolytic agents, it is important to consider the inflammatory profile of the target cell.
(A) Schematic showing experimental design for the Navitoclax experiments conducted on OIS and DDIS cells. (B–D) UMAP plots, cluster trees, and violin plots of primary and secondary senescence for DDIS cells treated with Navitoclax. (E–G) UMAP plots, cluster trees, and violin plots of primary and secondary senescence for OIS cells treated with Navitoclax. (H) Bar plots showing that Navitoclax preferentially induces apoptosis in cluster 2 cells for OIS and DDIS cells (ns: p > 0.05; *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ****p ≤ 0.0001).
(A) Model explaining the observed expression profiles in our datasets. (B) Violin plots showing the expression of collagen (middle) and ribosomal genes (left) in cluster 1 and cluster 2 for RS, OIS, and DDIS cells. (Right) Scatter plot showing that cells expressing higher levels of ribosomal reads express lower levels of collagen reads. Cells are colored according to their cluster 1 and cluster 2 assignments (ns: p > 0.05; *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ****p ≤ 0.0001).
Subpopulations differentially express collagen and ribosomal protein genes
Another feature we noticed in our analysis is that cells belonging to cluster 1 significantly accumulate ribosomal protein transcripts. Altered ribosome biogenesis is implicated heavily in senescence36, and we extend these findings further by showing that it is a process differentially regulated between subpopulations. We plotted cells according to their expression of ribosomal genes on the x-axis and expression of collagen genes on the y-axis and saw a strong anticorrelation for all datasets, suggesting that cluster 2 and secondary senescent cells are characterized by the high levels of collagen genes, and cluster 1 cells express high levels of ribosomal protein genes (Fig. 5).
Subpopulations expression profiles are identified in senescent endothelial cells
Because we identified these subpopulations in multiple forms on senescence in our clonal cell lines in addition to the IMR90 experiments conducted by Teo et al, we wanted to verify if these diversity profiles existed in another cell type. We downloaded single-cell transcriptomic datasets generated from Human Umbilical Vein Endothelial Cells (HUVEC) cell lines18. These data were collected from cells as they transitioned into RS. This means that there were many more subpopulations in the HUVEC dataset, and it was difficult to distinguish senescent form pre-senescent cells. However, we were still able to show major features of the diversity we saw in our own data.
We wanted to see if HUVEC cells expressing higher level of HMGA1 expressed low levels of collagen. We plotted cells according to their expression of HMGA1 (markers for cluster 1 and primary senescence) on the x-axis and expression of collagen genes (marker of cluster 2 and secondary senescence) on the y-axis. We saw an anticorrelation between these two features. We also plotted cells according to the number of reads mapping to ribosomal protein genes (which we saw higher in cluster 1) on the x-axis and saw a strong anticorrelation with cells that expressed high level of collagen reads. Moreover, reads mapping to ribosomal protein genes also anticorrelates with cells that have the high number of reads mapping to collagen genes. See Supplemental Figure 5, for this analysis.
Discussion
Our analysis compares subpopulations of senescent cells across different forms of senescence. We observe that DDIS, RS, and OIS cells are composed of two subpopulations. For our experiments, we used clonal cell lines that were infected with a 4-OHT HRAS:G12V inducible transgene that was only activated in OIS cultures. The fact that we have discovered subpopulations of senescent cells being generated from clonal cell lines suggests that the formation of these subpopulations is an inherent property of the senescent phenotype, and not due to a pre-existing heterogeneity that was present in proliferating cells. Therefore, these results offer additional insight into the formation of senescent cell subpopulations when compared with pre-existing scRNA-seq studies that have been generated from non-clonal cell lines. For instance, scRNA-seq studies on OIS samples that were infected with an inducible HRAS construct and then selected for with antibiotics yields a population of cells with varying amounts of HRAS expression that influences the diversity of the final senescent population.
It is important to mention that due to the experimental conditions, our cell yield in the various forms of senescence was different, which affected the sensitivity with which we could detect differential expression between the two clusters. Although this did not compromise our ability to detect common signatures, some of the genes that appear to be specific to one form of senescence and not the others might be the result of false negatives due to low cell numbers and statistical power.
For all forms of senescence, we have identified a subpopulation of cells that shows higher levels of extracellular matrix remodeling genes and genes that are regulated by TGF-beta such as collagen genes. Meanwhile, the other subpopulation shows higher levels of DNA damage, pro-inflammatory profiles, and HMGA1 expression, and a strong accumulation of ribosomal protein gene transcripts. These molecular signatures are reminiscent of what is seen in senescent cells inducing notch senescence in neighboring cells10–12. Consistently with this model, we see that the subpopulations we observe show similar gene expression patterns as described in the study conducted by Teo et al11, where a primary population of cells was driven into the senescent state through the activation of the HRAS oncogene, and, in turn, induced neighboring cells into senescence through the notch signaling pathway (juxtacrine-induced senescence), giving rise to the second cluster of cells characterized by the activation of the TGF-beta signaling pathway11. Our study extends this work by showing that these subpopulations appear to be a universal feature of senescence as they consistently appear in clonal cell lines and in all the main forms of senescence. One possible explanation is that that ability for OIS cells to induce secondary senescence though the notch signaling pathway is also shared by RS and DDIS cells (Fig. 5A).
Another possible explanation of the formation of these subpopulations is that senescent cells tend to express TGF-beta early on in their senescent lifespan, while expressing higher levels of inflammatory and interferon response genes later in their lifespan10. Our analysis could be capturing a transition from a senescence phase dominated by TGF-beta signaling to a phase displaying a more pro-inflammatory profile. If this is the case, then HMGA1 offers an interesting candidate whose role in this transition could be further investigated. As mentioned previously, HMGA1 has been shown to be downregulated in secondary senescent cells induced through a notch signaling pathway when compared with their primary OIS counterparts. An important difference between OIS and NIS is the expression of HMGA112,37,38. HMGA1 is a highly abundant chromatin-associated protein that is an essential component of the senescent chromatin architecture and is critical for the onset of the senescence program in OIS. OIS typically upregulates HMGA1, and this creates dramatically different chromatin structures compared to NIS cells12. Even though HMGA1 is a very important protein for the execution of the senescence program and the organization of the senescent chromatin, HMGA1’s role in senescence and its opposition of anti-inflammatory and TGF-beta-dominated SASP profiles has almost entirety been explored in the context of OIS and NIS. NIS is also characterized by high levels of TGF-beta signaling, which is consistent with a model where HMGA1 is anticorrelated with the expression of TGF-beta pathways10–12,37. Another interesting possibility is that the heterogeneity we are seeing is a result of differences in levels of DNA damage. We noticed that cells in cluster 1 are enriched in pathways relating to the DNA-damage response. Therefore, higher levels of DNA damage could predispose a cell to displaying the molecular signatures that we see in cluster 1 such as increased SASP, lower TGF-beta signaling, and higher levels of HMGA1. Our study sets the stage for important questions. For instance, how does HMGA1 interact with DNA-damage response pathways and TGF-beta signaling pathways in CS?
We also see that these subpopulations are translationally important by showing that treatment with the commonly used senolytic agent, Navitoclax, preferentially kills cells in cluster 2. A possible explanation for this is that the higher levels of TGF-beta signaling, which is a known affecter of the apoptotic pathways39, sensitizes cells to the apoptotic effects of Navitoclax35. Another explanation is that cells in cluster 1 are at a later stage of senescence, and therefore are more resistant to apoptosis and senolytic drugs. It is also possible that cells in cluster 1 are indeed primary senescent cells, and Navitoclax resistance is an inherent property of primary senescence. In any of these cases, this study points to the translational importance of senescent cell heterogeneity.
Acknowledgments
The N.N. laboratory was supported by Institutional Development Award (IDeA) grant P20GM109035 (Center for Computational Biology of Human Disease) from National Institutes of Health (NIH) National Institute of General Medical Sciences (NIGMS), grant 1R01AG050582-01A1 from NIH National Institute on Aging (NIA), and grant 1 UG3CA268202 from NIH National Cancer Institute (NCI) (SenNet Consortium). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or any other organization. Figure models were created with BioRender.com.
Declaration of Interests
The authors declare no competing conflict of interests.
Supplemental Information
Supplemental information can be found online at Supplementary.
References