Different types of aging cells marked by p16 and p21 in human and mouse tissues

Our further analysis revealed the presence of three subpopulations: (1) cells exclusively expressing p21; (2) cells exclusively expressing p16; and (3) cells expressing both p21 and p16 (Fig. 1C,D). The unbiased trajectory analysis using scVelo and PhyloVelo (describing the rate of gene expression change at a given time point) showed that p21+ and p16+ cells do not share a common developmental pathway, indicating that neither cell type arises from the other. This suggests that these populations are distinct and do not transition into one another (Wang et al, 2024b) (Fig. 1E). To gain deeper insights into the composition of the associated secretome within these identified subpopulations, we leveraged our previously established SenMayo gene set (Saul et al, 2022), which has been demonstrated to be commonly regulated across various age-related transcriptome datasets and is primarily composed of SASP-related factors (Brizio et al, 2024; Carapeto et al, 2024; Suryadevara et al, 2024). Interestingly, the variability of expressed SASP factors among p16+ cells was substantially higher, as demonstrated by their standard error, compared to that of p21+ cells. This suggests that the secretory phenotype in p21+ cells may be a more conserved or consistent mechanism—or that p21+ cells are more homogeneous (Fig. 1F). We observed a distinct secretory profile for p21+ and p16+ cells, marked by limited overlap between the two populations. Notably, only two genes, Cxcl16 and Plaur, were expressed in both subpopulations (Fig. 1G). Likewise, an expanded analysis of 1989 secreted protein genes (as listed by the Human Protein Atlas (Uhlén et al, 2015)), which included the secreted proteins plus SenMayo, revealed limited overlap between p21+ and p16+ cells. Only four genes, Col8a2, Plaur, Cxcl16 and Fbln5 were shared between the two subpopulations (Fig. EV1C). To eliminate potential cell type-specific biases between p21+ and p16+ cells, we focused on microglia and oligodendrocytes, as these represent the cell types with the highest proportion of p21+ and p16+ cells (Fig. 1B). We found that p21+ and p16+ microglia or oligodendrocytes had distinct secretory profiles, with only partial overlap in microglia, mirroring the differences seen in the overall population. This suggests that these variations are primarily driven by p16 or p21 expression rather than intrinsic cell type differences (Fig. EV1D,E). To further explore the differences between p21+ and p16+ cells and refine the RNA velocity analysis, we conducted a biased pseudotime analysis using Monocle3 (Cao et al, 2019), following the initial unbiased PhyloVelo approach. This allowed us to trace the dynamic trajectories of these subpopulations within the brain dataset, providing a more detailed understanding of their developmental pathways. Surprisingly, the results showed that the trajectory dynamics were independent of p16 or p21 expression status, with no evidence of direct lineage progression from p16+ to p21+ cells (Fig. EV1F,G).

To validate our initial observations, we analyzed an independent brain dataset from eight young and eight old whole mouse brains (excluding hindbrain (Ximerakis et al, 2019)). As observed in the first dataset, p21+ and p16+ cells showed minimal overlap and displayed distinct secretory profiles (Fig. EV2A–E). We also found that the majority of secretory markers were associated with p21+ cells and closely mirrored the patterns identified in the first dataset (Fig. EV2F). Consistent with our previous results, p16 expression significantly increased with age in the brain, whereas p21 upregulation was confined to specific cell clusters (Fig. EV2G,H). However, it is important to highlight that the first dataset was derived from isolated hippocampal tissue, whereas the second dataset included the entire brain. This distinction may account for some of the observed differences in cell distribution and secretory profiles between the two datasets.

Overall, these data emphasize the diversity as well as variability in p21+ and p16+ associated secretome profiles among cells with core features of senescence in the brain.

In our examination of all three tissues (Fig. 2), we observed a consistent pattern of distinct cell subpopulations that expressed either p21 or p16, with only a rare subset co-expressing both markers (Fig. 2B,F,J). Similar to our findings in the brain, RNA velocity and pseudotime analysis revealed that there was no evidence of a common ancestor between p21+ and p16+ cells, reflecting distinct developmental trajectories for each subpopulation. Furthermore, p16+ cells exhibited greater heterogeneity compared to p21+ cells (Fig. 2C,G,K).

In all three tissues, the p21+ cell population was considerably more abundant than the p16+ population. Using the SenMayo gene set, we found that the pattern observed in the brain was consistent across other tissues: while p16+ cells uniquely expressed a small subset of SenMayo factors, p21+ cells displayed a broader range of SenMayo genes (Fig. 2D,H,L). Although a small number of genes were co-expressed by cells positive for both p16 and p21, these were not consistently shared among the three tissues. Indeed, our investigations suggest that there is a distinctiveness in the composition of secretory profiles for each tissue, irrespective of the expression of p21 and p16. If we extend our analysis to a larger set of genes, including an extensive list from the Human Protein Atlas (Uhlén et al, 2015) plus SenMayo, we observe a consistent pattern across murine tissues: p21+ cells are associated with a larger secretome-associated transcriptional profile, with limited overlap in secreted gene expression between p21⁺ and p16⁺ cells (Fig. EV3A).

We confirmed the generality of our findings by performing the same analysis on separate independent scRNA-seq datasets from muscle (Fig. EV3B–J). Additionally, we assessed whether our analysis of the p16+ vs. p21+ secretome-associated transcriptional profile remained consistent when focusing on specific cell types. We found that a similar distribution of the p16+ vs. p21+ secretome was evident in fibro-adipogenic progenitors (FAPs) within skeletal muscle (Fig. EV3), in Lepr+ MSCs in bone (Fig. EV4A–J), and in endothelial cells in the liver (Fig. EV5A–J).

To determine if the two distinct cell populations are also present at the protein level in bone tissue, we analyzed data from mass cytometry by time-of-flight (CyTOF) on bones from young and aged mice (Doolittle et al, 2023). Similar to the scRNA-seq analysis, we excluded cells positive for Ki67 and CD45. Our analysis revealed distinct subsets of cells expressing either p21 or p16 exclusively, along with a rare subgroup co-expressing both markers at the protein level (Fig. EV6A–C). However, it is important to note that CyTOF relies on targeted panels of antibodies, which limited our analysis of secretory profile components. Most importantly, the lack of available SASP components makes it impossible to directly compare these proteomic data to our transcriptomic analyses. Nevertheless, we identified proteins such as Serpine1 exclusively expressed in p16+ cells and IL-6 exclusively expressed in p21+ cells, suggesting the existence of distinct secretory transcriptional profiles (Fig. EV6A–C). We further analyzed murine spleen (Fig. EV7A–E) and kidney (Fig. EV7F–J) using the Calico dataset (Kimmel et al, 2019). Consistent with our previous observations, p16+ and p21+ cells remained distinct populations in these tissues, each characterized by a unique secretome profile.

To determine whether the differences observed between p21+ and p16+ cell populations during aging are specific to age-related senescence or also occur in disease contexts, we compared our results with a mouse model of metabolic dysfunction-associated steatohepatitis (MASH) (Bendixen et al, 2024), a condition characterized by elevated senescence markers (Yashaswini et al, 2024) (Fig. EV8A–C). Similar to an aging liver, MASH showed minimal overlap between p21+ and p16+ cells (Fig. EV8B). However, the secretome transcriptional profiles associated with p21 and p16 in MASH were distinct from those observed in the aging liver (Fig. EV8D,E), suggesting that the composition of the secretory profile is influenced by the specific stimuli driving senescence in vivo.

Overall, our findings support the hypothesis that p21- and p16-expressing cells represent distinct senotypes with unique secretory profiles that vary depending on the tissue type and the senescence-inducing stimuli. This is in line with previous studies in which overexpression of p21 or p16 in mouse embryonic fibroblasts led to distinct secretome profiles (Sturmlechner et al, 2021), showing only partial overlap (Fig. EV9A,B).

To further substantiate our mRNA-based findings, we performed immunofluorescence imaging of sun-protected aged human skin (Fig. 3E), which confirmed at the protein level, distinct populations of p16⁺ and p21⁺ cells (Fig. 3F).

We independently confirmed these findings using another human skin dataset (Solé-Boldo et al, 2020) (Fig. EV10E–J). Because all our previous observations were derived from scRNA-seq, we sought to validate these findings using an independent transcriptomic approach that also preserves spatial information. This was particularly important, as cell dissociation required for scRNA-seq may result in the loss of fragile senescent cells. Applying MERFISH to human skin confirmed that p16⁺ and p21⁺ cells represent distinct populations across different skin cell types, each with a divergent SASP profile and only minimal overlap (Fig. EV10K–M).

A similar pattern was observed in the human lung. Here we utilized a recently published scRNA-seq dataset from healthy human lungs ranging from 21 to 78 years of age (Jia et al, 2023). Within our analysis, we identified 18 distinct cell populations (Fig. 3F). After excluding Ki67+, CD45^high and S-phase cells (Fig. EV11A–C), we observed that a significant majority of cells displayed p21 expression, while only a small number of cells exhibited p16 or the combination of p16 and p21 (Fig. 3G). There was no common ancestor, as demonstrated by RNA-velocity, and p16+ cells exhibited a heterogeneous secretome expression (Fig. 3H). Notably, cells that were p16+ demonstrated limited expression of SenMayo components, with Serpine1 being the exclusive gene expressed by this subgroup. In contrast, p21+ cells exhibited the expression of multiple SenMayo components (Fig. 3I). In addition, we analyzed alveolar type II (AT2) cells and again found distinct secretome transcriptional profiles associated with p21 and p16 (Fig. EV11D).

In an independent experiment, we analyzed human diseased lung tissue (COPD and pulmonary fibrosis/IPF) by immunofluorescence, which confirmed at the protein level the presence of distinct p16⁺ and p21⁺ cell populations (Fig. 3J).

We independently validated our findings in a murine lung dataset from the Calico dataset (Kimmel et al, 2019) (Fig. EV11E–J). When we broaden our examination to include a wider range of genes which involves a comprehensive list of secreted proteins (Uhlén et al, 2019) alongside SenMayo (Saul et al, 2022), we observed that there are relatively few secreted protein genes that are common between p21 and p16 expressing cells in both human and murine lung (Fig. EV11A,I,J). This further underscores the functional specificity inherent in these subpopulations.

In the case of p21+ cells, we observed that Icam1 was consistently expressed across all six analyzed tissues, while several other factors, including Jun, Edn1, Igfbp6, Tnfrsf1a, Cxcl16, Ets2, Fas, Igfbp4, Plaur, and Ccl2 were expressed in eight or more of the 14 analyzed tissues (Fig. 4D). Combining p16 (red) and p21 (blue) markers that are expressed in more than one (p16) and five (p21) tissues, a "core" secretome profile comprised of Icam1, Jun, Igfbp6, Tnfrsf1a, Edn1, Cxcl16, Igfbp4, Fas, and Plaur is established. While Icam1 may exhibit a stronger association with p21, Igfbp4 and Igfbp6 are prominently expressed by both cell types, suggesting these as "common" p16/p21-associated secretory factors (Fig. 4E).

It should be noted that SenMayo mainly consists of secreted genes (c.a. 73%), while SenSig just has 265 secreted genes (c.a. 15%, based on the Human Protein Atlas), and CellAge has 78 secreted genes (c.a. 15%, based on the Human Protein Atlas). In our search for a common p21+ and p16+ secretome signature, we broadened our analysis to encompass a greater number of secreted factors by using SenSig (Cherry et al, 2023). This approach enabled us to identify 29 genes expressed in p16+ cells that were consistently present in (at least) four of the analyzed tissues (Fig. 4F). In p21+ cells, we identified 24 distinct genes expressed in a minimum of 11 tissues, and a single gene, Cdkn1a, was expressed in all analyzed tissues (Fig. 4G). Upon combining all factors to identify a "common" senescence-associated phenotype, we observed that Cdkn1a, Btg2, Bcl2l11, Klf3, Atf3, and Ier5 remain abundant primarily in p21+ cells. Additionally, Cdkn2b is expressed in four out of five tissues by p16+ cells (Fig. 4H). Using CellAge (Avelar et al, 2020) as reference gene set, the p16-associated genes were (next to Cdkn2a) Cdkn2b and Lum in 6 of the 14 datasets (Fig. 4I). The p21-associated genes were (next to Cdkn1a) Btg2, Icam1, Mdm2, Pnrc1, Tap1 and Tspan13, as expressed in 12 of the 14 datasets (Fig. 4J). The "core" secretory phenotype was (after Cdkn2a and Cdkn1a) mostly consisting of Btg2, Stx12 and Tmbim1 in 13 of the 14 datasets (Fig. 4K).

To further validate our findings, we leveraged the p21-Cre mouse model (Wang et al, 2024a) in which p21-positive cells were GFP-tagged in murine visceral fat tissue (Fig. EV12A,B). Using this model, we confirmed that GFP/p21-positive cells exhibited a significant enrichment of the previously identified p21-associated SenMayo, SenSig, and CellAge scores in GFP/p21-positive cells (Fig. EV12C–G).

This heterogeneity is further reinforced when exploring the transcription factors that control gene expression in p16+ and p21+ cells. To achieve this, we employed SCENIC, a computational method that allows the prediction of interactions between transcription factors and target genes based on single-cell RNA-seq data (Aibar et al, 2017). SCENIC analyses reveal that transcription factors regulating gene expression in p16 (Fig. EV13A,B) and p21 (Fig. EV13C,D) positive cells were mostly tissue-specific, with minimal overlap between tissues. This implies that the secretory phenotype is not exclusively specific to p16+ or p21+ cells; rather, these cell subtypes are further transcriptionally regulated by factors specific to the tissue in which they are situated.

In summary, the diversity in cellular communication and transcriptional regulation between p16+ and p21+ cells within the unique environments of the six analyzed tissues underscores their heterogeneity and suggests potentially distinct functions.

The scRNA-seq data were aligned and quantified using the 10X Genomics Cell Ranger Software Suite (v6.1.1) against the murine reference genome (mm10) and human reference genome (hg19). The Seurat package (v4.3.0.1 and 5.0.0) (Butler et al, 2018; Stuart et al, 2019) was used to perform integrated analyses of single cells. Genes expressed in <3 cells and cells that expressed <200 genes and >20% mitochondria genes were excluded from downstream analysis in each sample. The datasets were SCTransform-normalized, and the top 3000 highly variable genes across cells were selected. The datasets were integrated based on anchors identified between datasets before principal component analysis (PCA) was performed for linear dimensional reduction. After normalization and scaling, a shared nearest neighbor (SNN) graph was constructed to identify clusters in the low-dimensional space (top 30 statistically significant principal components, PCs). An unbiased clustering according to the recommendations of the Seurat package was used, if not provided by the authors. The cell types were assigned according to the authors' recommendations or provided metadata. The alluvial plots were designed with the ggalluvial package (v0.12.5). For the DimPlots, the RNA slot was used, and every value above 0 was counted as "positive". After normalization and scaling, we proceeded with the differentially expressed genes analyses.

The differentially expressed markers were identified by the FindMarkers function (ident.1 was specified, and differences calculated to all other clusters) and the Wilcoxon signed-rank test. We used the SenMayo gene set (n = 125) or secreted proteins, obtained by the human protein atlas, augmented by the SenMayo secreted SASP factors (n = 1989 (Uhlén et al, 2019)) to select the SASP factors or secreted proteins. The cytoscape bubble plots were designed with cytoscape (v3.1.0), and the size of each bubble is proportional to the avg_log2FC. The circular plots were designed with ggplot2 (v3.4.4). The spider plots were designed with the package fmsb (v0.7.5). The circle size is proportional to the log2FC compared to all other clusters, while the color codes the respective tissue. The bars show the median of all tissues in which the respective gene is upregulated.

RNA Velocity analyses were performed using scVelo as an unbiased measurement of velocity. Here, unsliced and spliced mRNAs are used to predict cellular states, given that the unspliced counts precede the spliced counts (La Manno et al, 2018). After the direction of the velocity had been obtained, a depiction of the trajectory has been done via PhyloVelo, more closely translating the transcriptional dynamics of splicing kinetics with a likelihood-based, more dynamical model (Wang et al, 2024b). In some tissues, these analyses are further amended by a graphical representation with monocle3 (Cao et al, 2019).

Regarding the heterogeneity analyses, the standard error was calculated as standard deviation, divided by the square root of the sample size for each of the SenMayo SASP factors, separated by p16- or p21-positivity.

For the intercellular communication heatmap, we first calculated the respective intercellular communication via CellChat (v1.6.1). The centrality values were extracted and used for the central heatmap, for which we used GraphPad Prism (Version 9.0). The probability of the inferred communication was used for the sankey-network. The depiction of the sankey-network was done via the networkD3 package (v0.4) and exported via jsonlite (v1.8.7). Since muscle had too little p16+ cells for the CellChat analyses, these were excluded.

For the SCENIC analyses, we used the standard settings (v.1.3.1). For plotting the most important factors, the relative activity for each transcription factor above 1 was chosen for p21+ and P16+ cells, respectively, per tissue. Since the liver has very few p16+ cells for proper calculation with SCENIC, these were excluded. For the pie charts, the sum of transcription factors sharing a tissue was calculated, and the respective percentage is demonstrated within the plot.

The provided fcs-files were read into R by the flowCore package (v2.8.0) and transformed into a Seurat object. All subsequent analyses followed the standard Seurat procedure as described above.

The fastq files were mapped to the murine genome (mm10), and analysis was performed using the DESeq2 package (v1.38.3) as previously described (Saul et al, 2022). Significantly differentially regulated genes were selected by a Benjamini–Hochberg adjusted p value <0.05 and log2-fold changes above 0.5 or below −0.5.

Formalin-fixed, paraffin-embedded (FFPE) tissue sections (5 µm) were prepared from three sun-protected aged skin biopsies from individuals aged 73, 79, and 83 years, and two lung tissue sections from 75-year-old individuals diagnosed with IPF and with COPD, respectively.

Sections were deparaffinized in Histoclear (twice, 5 min each) and rehydrated through a graded ethanol series: 100% ethanol (twice, 5 min each), 90% ethanol (twice, 5 min each), 70% ethanol (1 × 5 min), followed by two washes in distilled water (5 min each). Antigen retrieval was performed by heating the slides to boiling in citrate buffer (pH 6.0; Agilent-Dako, S236984) for 10 min. Slides were allowed to cool at room temperature for 30 min and then rinsed twice in PBS (5 min each).

Tissue sections were blocked with 1% BSA and 1:60 normal goat serum (NGS; Agilent-Dako, X090710-8) for 30 min at room temperature, followed by overnight incubation with primary antibodies at 4 °C. The next day, sections were washed in PBS and incubated with appropriate secondary antibodies for 1 h at room temperature. After three additional PBS washes, slides were mounted using ProLong™ Gold Antifade Mountant with DAPI (Invitrogen).

The primary antibodies used were anti-p16 (725-4793) mouse monoclonal (dilution as provided by kit, VENTANA) and anti-p21 (2947S) (1:200, Cell Signalling).

A punch biopsy was collected from the upper thigh of a 66-year-old individual. Following collection, the sample was cleaned and the surrounding adipose tissue was carefully removed. The tissue was then centered in a plastic embedding mold, embedded in OCT compound, and rapidly frozen by immersion into super-cooled isopentane, which had been pre-chilled in a metal container placed within a liquid nitrogen bath. Using forceps, the mold was carefully lowered into the isopentane–liquid nitrogen bath without full submersion. Once the OCT was fully solidified, the tissue blocks were stored at –80 °C until shipment to the manufacturer.

MERFISH was performed by the manufacturer (Vizgen) using the pre-designed ImmunOncology gene panel in combination with a custom-designed probe targeting CDKN2A (ENST00000579755). The ImmunOncology panel included 550 genes encompassing markers of cellular senescence, cell-cycle regulation, inflammation, immune cell types and others. All steps of the MERFISH workflow, including tissue preparation, hybridization, imaging, and data acquisition, were conducted by Vizgen according to their standardized protocols. The resulting data were exported using the MERSCOPE Visualizer software, version 2.3 (Vizgen, Inc.).

The spatial coordinates were extracted and transferred into a UMAP using the Seurat package (v 5.0.0).

Statistical analyses were performed using either GraphPad Prism (Version 9.0) or R version 4.2.0. A p value <0.05 (two-tailed) was considered statistically significant.

Human sun-protected skin tissue samples were obtained from participants at Mayo Clinic and from the Baltimore Longitudinal Study on Aging (BLSA). At Mayo Clinic, skin biopsies were collected from consented participants under IRB protocol 21-008529, in accordance with the International Conference on Harmonization, Good Clinical Practice guidelines, and the Declaration of Helsinki. Samples from BLSA participants were collected under a protocol approved by the Institutional Review Board of the Intramural Research Program of the National Institutes of Health (protocol #03AG0325). Human lung tissues were obtained incidentally during clinically indicated thoracic surgeries under IRB-approved protocols (Mayo IRB 08-002518), with informed consent obtained prior to collection. All samples were de-identified before analysis.