Novel insights from a multiomics dissection of the Hayflick limit

Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Data availability
References
Article and author information
Metrics

Abstract

The process wherein dividing cells exhaust proliferative capacity and enter into replicative senescence has become a prominent model for cellular aging in vitro. Despite decades of study, this cellular state is not fully understood in culture and even much less so during aging. Here, we revisit Leonard Hayflick’s original observation of replicative senescence in WI-38 human lung fibroblasts equipped with a battery of modern techniques including RNA-seq, single-cell RNA-seq, proteomics, metabolomics, and ATAC-seq. We find evidence that the transition to a senescent state manifests early, increases gradually, and corresponds to a concomitant global increase in DNA accessibility in nucleolar and lamin associated domains. Furthermore, we demonstrate that senescent WI-38 cells acquire a striking resemblance to myofibroblasts in a process similar to the epithelial to mesenchymal transition (EMT) that is regulated by t YAP1/TEAD1 and TGF-β2. Lastly, we show that verteporfin inhibition of YAP1/TEAD1 activity in aged WI-38 cells robustly attenuates this gene expression program.

Editor's evaluation

This manuscript reports a unique, comprehensive, multi-omic resource for the study of replicative senescence. This resource encompasses temporal metabolomic, proteomic, bulk transcriptomic, single cell transcriptomic, and chromatin accessibility states of fibroblasts as they transition from proliferative to replicatively senescent. Hence, it will be considered a valuable resource by aging researchers.

https://doi.org/10.7554/eLife.70283.sa0

Introduction

Replicative senescence in animal cells growing in vitro was first discovered by Leonard Hayflick. He found that primary human diploid fibroblast cell lines ceased to proliferate after an extended number of serial passages (Hayflick, 1965). Since then, considerable work has been done to describe this phenomenon. A major causal feature of replicative senescence is telomere erosion, a process in which the telomeres gradually shorten with increasing cellular divisions. Eventually, the telomeres become uncapped which triggers a DNA damage response that results in cell cycle exit (Harley et al., 1990). It is understood that this is due to the absence of the telomerase reverse transcriptase (hTERT), the catalytic component of human telomerase which adds telomeric sequences to the ends of chromosomes to maintain telomere length in germ cells and stem cells. As hTERT activity is undetectable in normal human somatic cells, telomere attrition is a common aging phenotype hypothesized to underlie cellular senescence at the organismal level (López-Otín et al., 2013; Meyerson et al., 1997). Replicative senescence of human somatic cell lines in vitro can be avoided by overexpression of hTERT which prevents telomere shortening and confers apparently unlimited replicative capacity (Bodnar et al., 1998).

Beyond growth arrest and telomere shortening, phenotypic changes exhibited in replicatively senescent cells include the senescence associated secretory phenotype (SASP). Proteins of the SASP include proinflammatory cytokines, growth factors, angiogenic factors, and proteases. The SASP has been shown to play a role in paracrine signaling whereby senescent cells can promote local wound healing and/or drive healthy neighboring cells into senescence (Acosta et al., 2013; Coppé et al., 2006; Demaria et al., 2014; Coppé et al., 2008). Replicatively senescent cells also accumulate DNA and protein damage, accumulate lipids, and lose regulatory control of mitochondria and lysosomes (Gorgoulis et al., 2019).

The phenotypic similarity (telomere attrition, epigenetic alterations, mitochondrial dysfunction, and loss of proteostasis) between the cell autonomous aging hallmarks and in vitro senescence has led to the hypothesis that senescent cells in vivo play a causal role in organismal aging and aging-related diseases (Hernandez-Segura et al., 2018; Campisi and d’Adda di Fagagna, 2007; Sedelnikova et al., 2004; Jeyapalan and Sedivy, 2008; Herbig et al., 2006; Childs et al., 2015; López-Otín et al., 2013).

Consistent with this model, several age-related disease states can be directly linked to telomere length or telomerase activity. For example, up to 15% of familial idiopathic pulmonary fibrosis (IPF) cases arise from mutations in telomerase and up to 25% of sporadic cases occur in people with telomere lengths less than the 10th percentile. (Tsakiri et al., 2007; Armanios et al., 2007; Cronkhite et al., 2008; Alder et al., 2008; Stuart et al., 2015; Raghu et al., 2006; Duckworth et al., 2021; Stuart et al., 2014; Dai et al., 2015). Furthermore, the elimination of senescent cells in a number of age-related diseases, such as cardiac fibrosis, pulmonary fibrosis, neurodegenerative diseases, osteoporosis, and metabolic disorders have been argued to alleviate the disease state (Pignolo et al., 2020). Clinical trials for senolytics targeting fibrotic diseases, osteoporosis, frailty, and metabolic syndromes are currently underway (Borghesan et al., 2020).

The mechanism by which senescent cells might contribute to aging phenotypes is currently still unclear (Hernandez-Segura et al., 2017; Hernandez-Segura et al., 2018). In an effort to bring clarity to the replicative senescence process in vitro that in turn could elucidate in vivo function, we revisited and redesigned the original Hayflick experiment. Making use of recent advances in technology we tracked changes in cell state throughout the rplicative lifespan of the original Hayflick WI-38 cell line using bulk RNA-seq, single-cell RNA-seq (scRNA-seq), ATAC-seq, metabolomics and proteomics.

Overall, our data recapitulate many known features of the in vitro senescence process while simultaneously providing novel insight. First, the time resolution of our experiment coupled with single-cell trajectory analysis reveals that senescence is a gradual process that shares transcriptional, proteomic, and metabolomic features with epithelial-mesenchymal transition (EMT). Second, our metabolomic data identifies Nicotinamide N-methyltransferase (NNMT) activity as a potential initiating event in replicative senescence dependent loss of epigenetic silencing. Third, we show that these genomic regions that exhibit increased accessibility with increasing cellular age are concentrated in nucleolar/lamin associated domains and correspond with observed changes in the replicative senescence transcriptome. Lastly, integration across data modalities reveals that senescent WI-38 cells bear a strong resemblance to myofibroblasts. We provide bioinformatic and experimental evidence that the YAP1/TEAD1 transcription factor complex and TGF-β2 signaling are putative regulators of the transition to this state. Together our data suggests that a process similar to fibroblast to myofibroblast transition (FMT; analogous to EMT) is an intrinsic aspect of the replicative senescence phenotype in WI-38 fibroblasts.

Results

Transcriptomic profiles of replicative senescence, radiation-induced senescence, and increasing cellular density

To capture the replicative senescence process with high resolution we designed an experiment to continuously grow and intermittently sample cells from a starter batch of WI-38 cells that had undergone only 20 population doublings (PDL 20) (methods). To distinguish between replicative senescence-dependent changes and those arising from altered cell density and growth rate, we performed a cell density control study (methods). Briefly, early PDL cells were sampled at increasing levels of cell density to measure gene expression changes associated with cellular density and decreasing cell proliferation independent of replicative senescence (Figure 1A). In addition to the cell density control, we also included TERT immortalized WI-38 cells (hTERT) grown in parallel and sampled alongside wild type (WT) cells as a control for long term culturing and day-to-day sampling batch effects (Bodnar et al., 1998). The hTERT cell line was generated in advance of the experiment. Frozen aliquots of the hTERT line were thawed alongside and sampled in parallel with WT WI-38 cells. As expected, hTERT immortalized cells grew at a constant rate and did not slow or cease growth (Figure 1B).

Figure 1 with 5 supplements see all

Download asset Open asset

Expression dynamics of replicative senescence (RS), radiation induced senescence (RIS), increasing cell density (CD) and hTERT WI-38 cells.

(A) Experimental design for the RS, RIS, CD and hTERT experiments. (B) Days in culture vs. population doublings (PDL) for WT WI-38 cells (red) and hTERT immortalized cells (green). Labeled points denote sample collection and are expressed in PDLs for RS time course or PDL controls (PDL.ctrl) for the hTERT time course. Temporally paired samples indicated by vertical dotted lines. PDLs 52 and 53 were collected in a separate ‘deep’ senescence time course. (C) Percent cells staining positive for SA-β gal staining (y-axis) at increasing PDLs (x-axis) for WI-38 cells (red) and hTERT (green). Points represent replicate values. (D) Scatterplot of log2 fold changes in gene expression for PDL 50 vs. PDL 20 (x-axis) vs. senescence log2 fold change derived from a generalized linear model compiling gene expression changes across multiple fibroblast cell lines during replicative senescence (y-axis) (Hernandez-Segura et al., 2017). (E) Hierarchical clustering of significant gene expression changes (FDR adjusted p-value < 0.01) across all conditions (n=3 replicates). Values are log2 fold change vs. the average of the first time point of each condition. Reference time point not shown. (F) Significant (FDR adjusted p-value < 0.01) Gene Set Enrichment Analysis (GSEA) results for RS, RIS, and CD using the MSigDB Hallmarks annotation set. The -log10 p-value is colored by direction of enrichment (red=up, blue=down).

Figure 1—source data 1 RNA-seq quality control metrics.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig1-data1-v3.xlsx
Download elife-70283-fig1-data1-v3.xlsx
Figure 1—source data 2 RNA-seq Gene expression data and differential analysis.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig1-data2-v3.xlsx
Download elife-70283-fig1-data2-v3.xlsx
Figure 1—source data 3 RNA-seq GSEA.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig1-data3-v3.xlsx
Download elife-70283-fig1-data3-v3.xlsx
Figure 1—source data 4 RNA-seq GSEA for indivdual timepoints.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig1-data4-v3.xlsx
Download elife-70283-fig1-data4-v3.xlsx

We also included a radiation induced senescence condition to test for differences between replication and radiation induced senescence and to isolate changes arising from acute DNA damage response. Finally, we sampled proliferating WT and hTERT WI-38 cells at multiple PDLs for RNA-seq, scRNA-seq, proteomics, metabolomics and ATAC-seq until the WT WI-38 cells had reached senescence as measured by the cessation of growth (Figure 1B, Figure 1—figure supplement 1B).

We first examined the bulk RNA-seq data to compare and contrast replicative senescence with radiation induced senescence, cell density, and hTERT cells (quality control metrics located in Figure 1—source data 1). Differential gene expression analysis using RNA-seq revealed 6955 genes change with increasing PDL, 7065 genes change in response to ionizing radiation, and 9958 genes vary with increasing cell density (FDR adjusted p-value < 0.01, Figure 1—source data 2). Notably, the transcriptional changes we observed were consistent with previous studies. Figure 1D shows the strong correlation of the pattern of induced transcription established in our PDL 50 cultures and the fibroblast derived senescence-associated signature compiled across multiple fibroblast senescence gene expression experiments from cell lines derived from different tissue depots (Hernandez-Segura et al., 2017). We also observed induction of senescence-associated-β galactosidase activity (Figure 1C) as well as p16 and p21 in both our RNA-seq, proteomics data, and via western blot (Figure 1—figure supplement 2, Figure 1—figure supplement 3). These data reveal that the cells in our time course display classic features of senescent cells and exhibit transcriptomic changes that are highly correlated with previous studies of replicative senescence. Together, these data strongly suggest that our WI-38 cells successfully reached replicative senescence.

Hierarchical clustering of all significantly changing genes across the four conditions highlighted several important features of our experiment. First, as expected, gene expression in the hTERT-immortalized WI-38 cells remained largely stable. Second, replicative senescence, radiation induced senescence, and cell density exhibited many shared, but also unique, gene expression changes with respect to both the identity of differentially expressed genes and the magnitude of their changes (Figure 1E).

To facilitate a biological interpretation of these shifting transcriptomic landscapes, we applied gene set enrichment analysis (GSEA) using the MSigDB Hallmark annotation sets to learn which general processes are shared and distinct between replicative senescence, radiation induced senescence, and cell density (Figure 1F and Figure 1—source data 3; Liberzon et al., 2015; Subramanian et al., 2005; Mootha et al., 2003). All three perturbations (and not the immortalized cells), exhibited dramatic reductions in expression of genes belonging to S, G2M, and M cell cycle phases, consistent with cessation of cell division (Figure 1F - top cluster). Perhaps driven by this shift, we also found significant overlap in additional enriched annotations among genes induced in replicative senescence, radiation induced senescence and cell density, albeit with some variation in significance. Many of the shared enriched annotations can be categorized as stress responses, for example apoptosis, p53 pathway, inflammatory and interferon responses, and STAT3/5 signaling. In addition, we also observed enrichment of several developmental gene sets including myogenesis, angiogenesis, and adipogenesis.

We observed relatively few gene sets with discordant patterns across the replicative senescence, radiation induced senescence, and cell density experiments (Figure 1E). However, the gene set for epithelial to mesenchymal transition (EMT) was significantly enriched in genes that increased with replicative senescence as opposed to radiation induced senescence (Benjamini-Hochberg corrected p=5.5e-6 vs. Benjamini-Hochberg corrected p=0.12) or cell density wherein the term was actually underrepresented, but not significantly so (Benjamini-Hochberg corrected p=0.4). Although replicative senescence, radiation induced senescence, and cell density are highly similar at the abstracted level of enriched gene sets, it is clear in Figure 1E that there are many induced genes specific to replicative senescence. In addition, we applied GSEA to each individual time point for each condition (Figure 1—figure supplement 4, Figure 1—source data 4) and found the EMT gene set is enriched early and robustly during replicative senescence. Thus, the EMT gene set appears to represent a particularly important aspect of replicative senescence biology.

Single-cell RNA-Seq reveals that replicative senescence is a gradual process

Interestingly, the vast majority of gene expression changes evident by PDL 50 begin to manifest at much earlier PDLs. This observation is consistent with two distinct possibilities: (A) The senescence expression changes accrue early and gradually in the majority of cells without respect to proliferation status, or (B) the bulk RNA-seq profiles are a changing admixture of transcriptionally distinct mitotic, G1, and senescent cells. Both models have substantial support in the literature (Tang et al., 2019; Wiley et al., 2017; Nassrally et al., 2019; Smith and Whitney, 1980; Passos et al., 2007; Smith and Whitney, 1980). To discriminate between these possibilities, we employed single-cell RNA-seq (scRNA-seq) to directly measure percentages of senescent and cycling cells.

Briefly, we collected between 1000 and 2000 cells at increasing PDLs of WI-38s and matched hTERT time points at beginning and end of time course. We applied a modified (see methods) 10 X Genomics protocol to capture and preserve the more fragile senescent cells (Figure 2—figure supplement 1). We profiled a total of 11,000 cells (9000 WT and 2000 hTERT). When aggregated into pseudo-bulk expression profiles (i.e. the summation of all cells), the single cell results are highly correlated with expression changes observed in bulk RNA-seq (r=0.81), (Figure 2—figure supplement 2).

We first classified cells as either S, G2M, or G1 phase using canonical markers for the S and G2M phases as previously described Nestorowa et al., 2016. As expected, based on bulk transcriptomics, the number of cells in either S or G2M decreases with increasing PDL, whereas the cycling cell proportions are stable in the hTERT timepoints (Figure 2A and B).

Figure 2 with 5 supplements see all

Download asset Open asset

Cell cycle exit and distribution on approach to replicative senescence (RS) does not explain gradual increase in the replicative senescence transcriptome.

(A) Individual UMAP projections of WT WI-38 cells; PDL colored by phases of the cell cycle (G1 = green, G2/M = orange, S-phase = purple). (B) Bar graph of cell cycle state percentages defined by transcriptomic score (y-axis) by PDL (x-axis) for WT WI-38 cells (left) and hTERT WI-38 cells (right). (C) UMAP projections of WT and hTERT WI-38 cells by PDL and PDL.ctrl colored by different phases of the cell cycle defined by transcriptomic score. (D) UMAP projection of all WI-38 cells. Increasing PDLs (RS) or PDL control (hTERT) are colored from early (red) to late (blue). (E) The replicative senescence transcriptome manifests early in all phases of cell cycle. Heatmap of hierarchical clustering of gene expression values of differentially expressed genes as aggregated transcriptomic profiles for each cell cycle phase and PDL (left) versus all cells (right). Values are log2 fold change of each PDL against the reference PDL 25 (not shown).

Next, we projected all WT cells and two hTERT time points together to identify broad patterns in replicative aging PDL cells versus cycling hTERT cells. To do this, we employed UMAP projection as opposed to t-SNE as the former is faster and preserves local structure more so than the latter (McInnes et al., 2018). Cells from the two cell lines organize separately, however, share a similar geometry composed of a S/G2M phase roundabout and a large G1 lobe (Figure 2C). Although the overall pattern of WI-38 cells for both cultures was highly similar, we observed progression of cell grouping with PDL that had no concordance with increasing hTERT time points (Figure 2D - red to blue). Specifically, PDL 50 cells (but not the temporally paired hTERT PDL.ctrl 50 time point) organize apart from all other cells within the UMAP projection. To test if these cells were senescent, we scored all cells with a senescent gene expression signature derived from the bulk RNA-seq and found that the highest scoring cells reside within this grouping composed of PDL 50 cells (methods) (Figure 2—figure supplement 4).

To test model A vs. B, we compared the pseudo-bulk gene expression profiles across the cell cycle at each PDL. We reasoned that if replicative senescence were gradual, we would observe the replicative senescence gene expression pattern even in young cycling cells without respect to G1/S/G2M. Indeed, the replicative senescence signature is present in all cell cycle phases even in early PDLs (Figure 2E- left panels vs. far right panel). In addition, bioinformatic segregation and UMAP projection of cells scored as S/G2M phase cells revealed that the dominant source of variance in cycling cells is the PDL (Figure 2—figure supplement 3). Lastly, we did not observe any cells < PDL 50 in close association with the senescent cells in the UMAP projection PDLs (Figure 2C and D).

It has previously been reported and observed here (Figure 2—figure supplement 5) that cell size increases with PDL in primary cell lines (Ogrodnik et al., 2019). We hypothesized that perhaps the higher PDL cells by nature of their larger size, might reach contact inhibition faster than early PDL cells. In this scenario, given the high correlation between the gene expression signatures of senescence and cell density, the increase in senescence signature could be an artifact of cell culture arising from changing morphology. To test this possibility, we scored our single cells again with modified signatures for senescence and cell density composed only of genes uniquely induced in either perturbation (Figure 2—figure supplement 4). We found that the specific senescent signature increases with PDL in WT WI-38 cells in contrast with the specific cell density signature which stays relatively flat. These data provide more evidence supporting a gradual increase in the senescent gene expression program that is independent from changes in cell density.

Together, these results argue in favor of a gradual model (A) of replicative senescence wherein cells ramp up expression of the replicative senescence program with increasing PDL even when still proliferative. We cannot rule out there are a small minority of cells that enter senescence early in the experiment as previously reported (Smith and Whitney, 1980; Passos et al., 2007). However, for the cells that made it to PDL 50, we observed that the senescent gene expression program manifests in their primogenitors. These data are also consistent with the previously reported changes in cellular phenotype (larger cell size and increased cycling time) with increasing PDL (Macieira-Coelho and Azzarone, 1982; Ogrodnik et al., 2019; Neurohr et al., 2019; Absher et al., 1974). Importantly, this observation suggests that aspects of cellular senescence are present in non-senescent cells. This result raises the intriguing possibility that the reported pathological features of senescent or senescent-like cells in vivo might also manifest in cells that are not classically senescent.

Proteomic landscape of replicative senescence WI-38 cells suggests large metabolic alterations

We next turned to our proteomics data generated from the same cultures and time points as above. We obtained high confidence measurements for ~8500 proteins (Figure 3—figure supplement 1, Figure 3—source data 1). Similar to our transcriptional results, the hTERT samples did not exhibit large changes in the proteome with experimental progression (Figure 3—figure supplement 1). Overall, we observed high concordance between transcript and protein levels with very few outliers (r=0.71 for PDL 50 vs. PDL 20, Figure 3—figure supplement 2). We employed GSEA using hallmark annotations and observed a sharp depletion of proliferation and mitosis associated gene sets supporting our previous findings (Figure 3A and Figure 3—source data 1). Likewise, the proteomics data recapitulated enrichment of the EMT annotation set along with the adipogenesis and myogenesis.

Figure 3 with 5 supplements see all

Download asset Open asset

Proteomic and metabolomic changes during replicative senescence.

(A) Significant (Benjamini-Hochberg adjusted pvalue < 0.01) GSEA results for protein changes (n=3) using the MSigDB Hallmarks. Reference time point not shown. Values are -log10 p-value and are colored by the direction (orange=up, green=down). (B) KEGG and GO term enrichment of genes from A for selected annotations. (C) Gene expression changes (n=3) for insulin-independent glucose transporter *SLC2A1* (left) and fatty acid transporter *CD36* (right) in the hTERT, replicative senescence, radiation induced senescence and cell density time courses. (D) Diagram of glycolysis, glycolytic shunts, fatty acid import and oxidation, TCA cycle/oxidative phosphorylation. Metabolites are blue, proteins are green. (E) Metabolite (n=4) and protein changes (n=3) with senescence from D. Average corrected metabolite log2 fold change vs. PDL 25 as reference (not shown) are plotted. Proteomics values are average log2 fold change vs. PDL 20 as reference (not shown). *:Significant changes (FDR adjusted p-value < 0.05 metabolomics, FDR adjusted p < 0.01 proteomics).

Figure 3—source data 1 Proteomics quantification, differential analysis and GSEA.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig3-data1-v3.xlsx
Download elife-70283-fig3-data1-v3.xlsx
Figure 3—source data 2 Metabolomics quantification and differential analysis.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig3-data2-v3.xlsx
Download elife-70283-fig3-data2-v3.xlsx

Multiple enriched sets pointed to replicative senescence-dependent shifts in cellular energy utilization (Figure 3A). Of note, we found that although under-enriched in our replicative senescence transcriptome (Figure 1F), The oxidative phosphorylation hallmark was highly enriched in the proteomics data. We plotted the oxidative phosphorylation genes driving the set enrichment in Figure 3—figure supplement 2. This projection revealed that many genes in this set fall into a quadrant (positive proteomic, negative transcriptomic) suggesting discordant regulation of these genes.

The oxidative phosphorylation hallmark set is composed of multiple closely linked mitochondrial complexes and functions that regulate cellular energy flux, for example TCA cycle, fatty acid oxidation, pyruvate metabolism, ATP synthase etc. We divided the oxidative phosphorylation gene set into into these constituent parts and visualized the changes in proteins with increasing PDL. We observed PDL-dependent increased in proteins from all categories except for mitochondrial assembly and structures. These results are consistent with altered mitochondrial function (Figure 3—figure supplement 3).

Next, we examined the annotations for the highest enriched hallmark sets from Figure 3A and found strong enrichment of multiple KEGG metabolic pathways. (Figure 3B). In addition, both the glucose transporter SLC2A1 and the fatty acid scavenger CD36 transcripts exhibited strong, early, and replicative senescence specific up-regulation suggesting altered cellular energy requirements (Figure 3C). Together these results point to a drastically altered replicative senescence metabolic landscape.

Senescent WI-38 cells exhibit increased utilization of fatty acid metabolism and glycolytic shunts

To generate a metabolic profile for replicative senescence in WI-38 cells, we harvested samples at increasing PDLs alongside input material for all other data types (methods). From our metabolomics data, we identified 285 compounds (Figure 3—figure supplement 4, Figure 3—source data 1). To remove batch effects that were observed in several PDL and the matched hTERT samples, we corrected each PDL with its paired hTERT sample. We then calculated significant changes over time for each metabolite (FDR adjusted p < 0.05, Figure 3—source data 2). To further guide our metabolic analysis, we focused our queries based on the enriched metabolic pathways found in the proteomic analysis. Specifically, we studied the changes to glycolysis, oxidative phosphorylation and fatty acid metabolism (Figure 3D).

Following the glucose and the increased expression of the glucose importer SLC2A1, we looked at the central stem of glycolysis as a potential source of energy fueling the altered senescent metabolism suggested by our proteomics data. (Figure 3D and E). We found that metabolites and enzymes dedicated to pyruvate generation and/or lactate production did not change in a concordant manner. We instead found that most of the glycolytic shunts exhibited upregulation at both the metabolite and protein level (Figure 3E). These results suggest an increase in allocation of glucose for manufacture of various biomolecules and their precursors (glycogen, hexosamines, phospholipids) alongside or in lieu of pyruvate or lactate production (Figure 3E).

Metabolites and enzymes involved in fatty acid import (e.g. CD36) and oxidation for the purpose of energy generation appear strongly up-regulated (Figure 3E). We also saw that metabolites involved in the de novo production of phospholipids via the Kennedy pathway are highly upregulated, specifically phosphotidylethanolamines precursors (Figure 3—figure supplement 5).

From these data, it is clear that WI-38 cells approaching replicative senescence undergo drastic shifts in metabolism. Specifically, we observed increased glucose utilization in glycolytic shunts coupled with an increase in fatty acid import and oxidation. It is possible that the senescent cells are switching to fatty acid oxidation to fuel increased TCA cycling and oxidative phosphorylation as glucose is diverted to macromolecule production. Indeed, metabolomic data collected from various types of senescence models has shown increased fatty acid oxidation, lipid accumulation, TCA up-regulation, and glycolytic alterations (Zwerschke et al., 2003; Ogrodnik et al., 2017; Flor et al., 2017; Unterluggauer et al., 2008; Johmura et al., 2021). Likewise, consensus metabolomic findings in EMT models report increased TCA cycle products, altered lipid metabolism, and activated hexosamine pathway (Hua et al., 2020; Corbet et al., 2020; Lucena et al., 2016). These data provide functional metabolic evidence supporting a connection between replicative senescence, metabolic rewiring, and the increased expression of EMT-related genes.

Nicotinamide N-methyltransferase (NNMT) links nicotinamide adenine dinucleotide (NAD) and methionine metabolism as a putative heterochromatin regulator

The interplay between metabolism and the epigenome is an increasingly recognized mechanism through which nutrient availability shapes chromatin and regulates cell state (Lu and Thompson, 2012). The metabolites nicotinamide adenine dinucleotide (NAD+) and methionine metabolism power the deacetylation and methylation of histone tails and DNA required for maintaining repressive DNA conformations (Lu and Thompson, 2012). In a compelling intersection, perturbation of NAD, methionine, and heterochromatin levels have all previously been reported in multiple aging contexts including replicative senescence (James et al., 2016; López-Otín et al., 2016; Kozieł et al., 2014; Yousefzadeh et al., 2020; Benayoun et al., 2015). In review of our replicative senescence data for the metabolic and proteomic components of the NAD and methionine pathways, we found that one of the largest and earliest changes for any protein (or transcript) was the increased expression of Nicotinamide N-methyltransferase (NNMT) (Figure 4A and B).

Figure 4

Download asset Open asset

Nicotinamide n-methyltransferase (NNMT) links nicotinamide adenine dinucleotide (NAD) and methionine metabolism during replicative senescence.

(A) Gene expression fold changes for the NNMT (top) in the hTERT, replicative senescence (RS), radiation induced senescence (RIS), and cell density (CD) time courses (three replicate average, each point a replicate). Heatmaps of metabolite and protein changes during replicative senescence from B (bottom-replicate average). (B) Metabolic diagram of the methionine and NAD salvage pathways. Metabolites are blue, proteins are green. Shading indicates inferred pathway direction during replicative senescence based on metabolite/protein changes in A. (C) Heatmaps of metabolite (four replicate average) and protein (three replicate average) changes with replicative senescence from B. Metabolomics data was batch corrected against hTERT samples and expressed as log2 fold change vs. PDL 25 as reference (not shown). Proteomics data expressed as log2 fold change vs. PDL 20 as reference (not shown). Significant changes (FDR adjusted p < 0.05 metabolomics, FDR adjusted p-value < 0.01 proteomics) during replicative senescence are denoted with asterisks. (D) Potential model for NNMT and metabolic regulation of chromatin state during replicative senescence.

DNA and histone methylation require abundant levels of the universal methyl donor S-adenosyl methionine (SAM). NNMT not only depletes SAM by catalyzing the removal of its methyl group, it does so by fusing the methyl group to the NAD+ precursor nicotinamide (NAM) resulting in the generation of methyl nicotinamide (MNA) (Ulanovskaya et al., 2013). Thus, NNMT effectively acts as a sink for two of the primary metabolites a cell requires for silencing chromatin and regulating gene expression (Komatsu et al., 2018). We found that MNA levels mirror those of NNMT; increasing early and robustly during replicative senescence progression consistent with high NNMT activity. Both SAM and is methyl-depleted counterpart, S-adenosyl homocysteine (SAH), were depleted with replicative senescence albeit to a lesser extent than the observed increase in MNA. NAMPT and NMNAT1, enzymes in the NAD salvage pathway, were weakly down-regulated at the protein level (Figure 4A and B).

High NNMT expression has been implicated in SAM depletion in embryonic stem cells and cancer-associated fibroblasts (CAFs) (Eckert et al., 2019; Sperber et al., 2015). In both cases, the functional consequences were similar; DNA hypomethylation and decreased capacity to form or maintain silenced loci and heterochromatin. The important point is that the loss of silencing is actively promoted through NNMT activity. One intriguing hypothesis is that increased NNMT activity with replicative senescence promotes loss of silencing and heterochromatin via SAM depletion (Figure 4C).

Increased DNA accessibility and transcription from nucleolar and lamin-associated domains is a dominant feature of the replicative senescence epigenome

To study genome-wide changes in the epigenetic landscape of senescent cells, we collected ATAC-seq data during the replicative senescence time course. ATAC-seq libraries from all PDLs and hTERT PDL controls exhibited similar ATAC-seq fragment size distribution emblematic of nucleosomal patterning, high alignment rates, very low mitochondrial read percentages, and strong enrichment around transcriptional start sites. To enrich for nucleosome free regions of DNA, we applied size selection to our ATAC-seq libraries prior to sequencing (methods, Figure 5—figure supplement 1, Figure 5—figure supplement 2, Figure 5—figure supplement 3 and Figure 1—source data 1). These quality control metrics are indicative of high-quality ATAC-seq libraries at all sampling time points and PDLs.

We next quantified the read distribution across chromatin states previously annotated for another human fetal lung fibroblast cell line (IMR-90, Figure 5—source data 1; Ernst and Kellis, 2015). We categorized 25 distinct states into four broad categories: promoters, enhancers, transcription, and miscellaneous, the last of which is composed largely of heterochromatic and undefined states, for example H3K9me3 and H3K27me3-rich regions, zinc-finger repeats, and quiescent (defined by the absence of histone marks, accessibility, or gene features) which we refer to as ‘undefined’ (Ernst and Kellis, 2015).

In comparing the read distribution across chromatin states from hTERT immortalized cells and WT WI-38 cells, we observed an increase in ATAC-seq reads falling into the undefined chromatin regions at the relative cost to all other states (Figure 5A-left panel). We next plotted the read distribution across chromatin states as a function of increasing time for both cell lines and found these shifts to be associated with increasing PDL in the WT WI-38 cells; hTERT ATAC-seq read proportions across chromatin states remained stable over time (Figure 5A-right panel).

Figure 5 with 5 supplements see all

Download asset Open asset

Increased accessibility within heterochromatin and nucleolar associated domains (NADs) is a dominant feature of the replicative senescence (RS) epigenome.

(A) Percent of ATAC-seq reads falling into four broad chromatin states compiled from the ENCODE IMR-90 25 chromatin state prediction for all samples in WT or hTERT WI-38 cells. Percent of all ATAC-seq reads within two specific states (y-axis) vs PDL (x-axis) for WT (red) and hTERT (green) WI-38 cells (bar height denotes average of 3 replicates). (B) Median log2 fold change in ATAC-seq counts for significantly changing ATAC-seq peaks binned into ENCODE IMR-90 25 chromatin states (row) for each sample (column, n=2 or three replicates). Reference samples are PDL 20 or PDL 20.ctrl for WT and hTERT WI-38 cells respectively (not shown). (C) Observed and expected distribution of overlaps of significantly induced genes (from Figure 1, RNA-seq) with NADs and LADs. Replicative senescence (red), radiation-induced senescence (RIS) (orange), and cell density (CD) (blue). Median expected number of overlaps and observed number of overlaps shown by black and green lines respectively. (D) Average (n=2 or three replicates) log2 fold change distribution (PDL 50 vs. PDL 20) for all ATAC-seq peaks that overlap or are excluded from NADs (left) or LADs (right).

Figure 5—source data 1 ENCODE IMR90 chromatin state labels.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig5-data1-v3.txt
Download elife-70283-fig5-data1-v3.txt
Figure 5—source data 2 ATAC-seq peak coordinates and metadata.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig5-data2-v3.xlsx
Download elife-70283-fig5-data2-v3.xlsx
Figure 5—source data 3 Normalized ATAC-seq coverage in peaks.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig5-data3-v3.zip
Download elife-70283-fig5-data3-v3.zip
Figure 5—source data 4 ATAC-seq fraction of reads in peaks.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig5-data4-v3.txt
Download elife-70283-fig5-data4-v3.txt
Figure 5—source data 5 IMR90 LADs from Salama 2013.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig5-data5-v3.txt
Download elife-70283-fig5-data5-v3.txt
Figure 5—source data 6 IMR90 NADs from Nemeth 2017.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig5-data6-v3.txt
Download elife-70283-fig5-data6-v3.txt

To determine if this shift in accessibility to undefined and heterochromatic states was indicative of increasing noise versus coherent changes in accessibility, we used the ATAC-seq reads to call peaks of accessibility. From the WT and hTERT ATAC-seq data from all time points, we identified 363,470 ATAC-seq peaks (Figure 5—source data 2 and Figure 5—source data 3). We found a slight but significant (~5%, p=0.014,linear regression) linear PDL-dependent decrease in the fraction of reads in peaks during replicative senescence (Figure 5—figure supplement 4A, Figure 5—source data 4). We tested experimentally whether the increased heterochromatic accessibility was an artifact of increasing numbers of dead cells in later PDL samples by repeating the experiment using late (PDL 45) cells with a cross-linking agent (propidium monoazide-PMA) that renders DNA from dead cells inert as previously described (Hendrickson et al., 2018). No significant change in read distribution across chromatin states was observed with addition of PMA (Figure 5—figure supplement 4B).

We then divided the ATAC-seq peaks into each of the 25 discrete chromatin states and calculated the fold change in accessibility for each state and PDL compared to the first time point. We observed that changes in peak accessibility mirrors that for all reads. Again we observed a clear increase in accessibility with PDL across the miscellaneous category with a concomitant decrease in all other states in WT WI-38 cells and not in immortalized hTERTs (Figure 5B).

The undefined state from the miscellaneous category piqued our interest for two reasons, one being that it alone accounted for 20–40% of all ATAC-seq reads and a large 10% increase with senescence. The second point of interest was the observed early onset of change (Figure 5A and B). In looking to define the undefined state, we found that Dillinger et al., 2017 previously reported that the undefined state is largely overlapping with both lamin-associated domains (LADs) and nucleolar-associated domains (NADs). In support of this observation, we too found that the undefined chromatin state is markedly gene poor and overwhelmingly overlapping with both empirically defined LADs and NADs (Z-score 90 and 150) that had previously reported in IMR-90 fibroblasts (Figure 5—figure supplement 5 Figure 5—figure supplement 4C, Figure 5—source data 5; Dillinger et al., 2017; Sadaie et al., 2013). For reference, Figure 5—figure supplement 4E provides a chromosome level view of the significant overlap between replicative senescence accessible ATAC-seq peaks in the undefined state with NADs and LADs compared to all peaks annotated in WI-38 cells.

To control for the fact that the chromatin state annotations we used were derived from a different, albeit closely related cell line, we performed the same analysis using gene annotations which are cell line independent. We found that peaks in intergenic regions increase in accessibility with PDL on average which is consistent with the undefined (and gene poor) chromatin state annotation from IMR-90 cells. These results suggest that the chromatin states defined in IMR-90 cells are capturing the chromatin states in WI-38 cells (Figure 5—figure supplement 4D).

We tested and found that ATAC-seq peaks falling within NADs and LADs exhibit increased accessibility with replicative senescence (p < 10 e-16 for both: Wilcoxon rank-sum test). Notably, this trend was much greater in magnitude for peaks in NADs versus peaks in LADs, with a median log2 fold change of 0.98 for NADs vs. 0.24 for LADs (Figure 5D).

Digging deeper into NADs, we next asked the chromatin residing within NADs might respond differently than the rest of the genome to increasing PDL. We repeated the analysis shown in Figure 5B with two subsets; peaks excluded from NADs or peaks overlapping NADS (Figure 5—figure supplement 5). Dividing these two sets of peaks into the 25 chromatin states revealed that as expected, NAD excluded peaks mirror the pattern observed for significantly changing peaks across the whole genome in the original analysis. Conversely peaks overlapping NADs were found to increase in accessibility in a variety of chromatin states (e.g. poised promoters and many classes of enhancers) that exhibit reduced accessibility across the rest of the genome with increasing PDL.

One interpretation is that although the aggregate peak accessibility data appear to suggest that transcription is globally shutting down, the NAD specific view reveals rather that transcription is induced in a subset of loci. Although widespread, the reduced accessibility in euchromatin we observed in aggregate is not reflective of a pervasive shutdown of transcription across the entire genome. To formally test this hypothesis, we next asked if the shifts in NAD/LAD accessibility correlate with productive transcription and found that indeed, senescence induced genes significantly overlap NAD domains. Radiation and cell density induced genes exhibited weaker (cell density) or nonexistent (radiation) overlap with NADs and LADs. These data suggest that increasing accessibility in LADs and specifically NADs is a prominent feature of the replicative senescence epigenome that is directly connected to the senescent gene expression program (Figure 5C). Together these results highlight a striking increase in accessibility within nucleolar associated DNA that connects changes in the transcription with a global shift in the epigenome.

Transcriptional regulators of the replicative senescence transcriptome and epigenome

To parse out the regulatory logic of the senescence gene expression program, we leveraged our ATAC-seq data to gain insight into which transcription factors regulate replicative senescence accessibility via transcription factor motif analysis. Having cataloged a universe of ATAC-seq peaks with significant changes in accessibility, we next assigned peaks to neighboring genes. Taking the top replicative senescence differentially expressed genes, we searched the proximal ATAC-seq peaks for enriched transcription factor motifs (Figure 6A). We found enriched motifs for TEAD1, CEBP family, SMAD family, AP1 complex, and FOX family transcription factors.

Figure 6 with 1 supplement see all

Download asset Open asset

Master transcriptional regulators of replicative senescence.

(A) Scatter plot of transcription factor motif enrichment in ATAC-seq peaks surrounding significantly induced genes (FDR adjusted p-value < 0.01, log2 FC >0.5) during replicative senescence. The Y-axis is the -log10 p-value; x-axis is the percent of input genes with the transcription factor motif. All points are transcription factors. Specific transcription factors of interest are called out in red. (B) Bar graph of ridge regression coefficient of motif predictive power in model of increasing peak accessibility with replicative senescence. Transcription factors of interest are highlighted in red. (C) Scatterplot of transcription factors enriched for binding in regulatory regions around replicative senescence depleted genes (y-axis, FDR adjusted p-value < 0.01, log2 FC <–0.5) vs. replicative senescence induced genes (x-axis, FDR adjusted p-value < 0.01, log2 FC <–0.5). Curated cell cycle transcription factors are colored in green; transcription factors of interest e.g. EMT/AP1/YAP1/TEAD1 etc. are colored in red. (D) Scatterplot of enriched transcription factors rank for binding enrichment in regulatory regions around senescence induced genes in astrocytes (y-axis) vs. replicative senescence induced genes in WI-38 cells (x-axis). transcription factors with discordant ranks/enrichment are colored in turquoise,transcription factors with concordant ranks/enrichment are colored in red; transcription factors of interest are labeled in red.

We also applied an orthogonal gene-independent methodology for determining which motifs are predictive of replicative senescence induced increases in ATAC-seq peak accessibility. Consistently, we again found TEAD1 motifs to be the most predictive feature (mean coefficient = 0.25 across 10 models-methods). In addition, we also found evidence for FOXE1 and SMAD1 regulation as well as other senescence related transcription factors, for example CEBPB and TP53 (Figure 6B).

Out of the FOX family transcription factors, FOXE1 is unique in that its transcript exhibits one of the most specific (no change in hTERT, cell density, or radiation) and largest increases in replicative senescence gene expression (40–50 x) (Figure 6—figure supplement 1). FOXE1 is annotated as a thyroid-specific transcription factor with putative roles in thyroid development and cancer (Zannini et al., 1997; Kallel et al., 2010; Ding et al., 2019). It has also been reported to regulate two signaling molecules upregulated in the replicative senescence time course; TGFB3 and WNT5A (Figure 1—source data 2; Venza et al., 2011; Pereira et al., 2015). Furthermore, it has been reported that many of the FOX family transcription factors are relatively promiscuous binders of each others canonical motifs. Thus, it is possible that the increased accessibility in peaks with FOX motifs during replicative senescence may be driven by FOXE1 activity Rogers et al., 2019. Despite the specificity and magnitude of replicative senescence induction, FOXE1’s function in replicative senescence is unclear and warrants further investigation.

We returned to our bulk RNA-seq to test for enrichment of protein-DNA binding events (mapped by ENCODE) in regulatory elements proximal to our replicative senescence differentially expressed genes by using Landscape in silico deletion Analysis (LISA) (Qin et al., 2020). Plotting transcription factor enrichment for genes depleted with replicative senescence against genes induced during replicative senescence revealed three broad cohorts of transcription factors; proliferating cell transcription factors, replicative senescence transcription factors, and transcription factors whose binding was enriched around both sets of genes (Figure 6C). As expected, proliferation specific transcription factors are replete with cell-cycle-specific transcription factors for example E2F and RBP family transcription factors. A large portion of the transcription factors exhibiting replicative senescence specificity belong to four categories: inflammation transcription factors (NFKB, CEBPB), AP1 subunits (JUN,JUND,FOSL2), YAP1-TEAD1 components (TEAD1/4, YAP1, WWTR1), and EMT transcription factors (SNAI2, TCF21).

TEAD1 is a member of the TEA domain transcription factors whose functions range across a wide swath of biology depending on context and binding partner. TEAD transcription factors cannot induce gene expression without a cofactor, which is most often YAP1 (yes-associated protein 1) a key downstream effector of Hippo signaling (Azakie et al., 2005; Chen et al., 1994; Benhaddou et al., 2012; Yu et al., 2015; Landin Malt et al., 2012; Vassilev et al., 2001; Zhao et al., 2008; Ma et al., 2019; Piccolo et al., 2014). Consistent with our identification of TEAD1 and YAP1 as replicative senescence regulators, YAP1 activation has been tied to EMT, anti-apoptosis, telomere dysfunction, inflammation, and positive regulation of fatty acid oxidation (Kurppa et al., 2020; Lee et al., 2019; Peng et al., 2017; Zhang et al., 2017; Chakravarti et al., 2020).

Lastly, given that YAP1/TEAD1 activity appears to increase during replicative senescence against the backdrop of an altered epigenomic context, we tested for an interaction between TEAD1 motifs and the increasingly accessible NADs and LADs. We found that not only are TEAD1 sites significantly enriched within NADs and LADs (p < 7.17e-07, and p < 6.12e-03, respectively hypergeometric), but we also discovered a greater-than-additive increase in accessibility with replicative senescence for TEAD1 motifs that occur within a NAD (p < 3.12 e-07). These results suggest that there is a functional connection between TEAD1 activity and NAD domains during replicative senescence. Thus, The changing epigenetic context of replicative senescence may alter TEAD1 binding and target expression induction with increasing PDL.

Collectively, these analyses uncover a common theme amongst putative regulatory transcription factors; Hippo signaling (YAP1/TEAD1), EMT transcription factors, and TGF-β signaling (SNAI2, SMAD activity). SNAI2 has been shown to work in tandem with the YAP1/TEAD1 complex and these pathways often work towards similar biological ends (Tang et al., 2016; Kurppa et al., 2020). Together, these transcription factors are reported as highly involved with proliferation, EMT, ECM production, fibrosis, and apoptosis avoidance (Kim et al., 2019a).

Lastly, given that the LISA results are based on binding events collected from a vast multitude of cell lines, we wondered if the same transcription factors might be found to regulate senescence in a completely different cellular context. To test this, we took significantly induced genes from a senescence model using astrocytes and oxidative stress as the senescence trigger (Crowe et al., 2016). Plotting transcription factor enrichment for genes induced in astrocyte senescence against transcription factor enrichment for WI-38 replicative senescence genes revealed that although there were a substantial number of discordant transcription factors, there was a clear population of transcription factors highly enriched in both senescence models (Figure 6D). Notably, the top concordant transcription factors ranked at the top in both contexts and recapitulate all the previous results for example YAP1/TEAD1, SNAI2, CEBP family transcription factors and AP1 subunits.

scRNA-Seq trajectory analysis resolves WI-38 cells’ approach to replicative senescence

A wealth of recent work in single-cell transcriptomics has demonstrated that ordering single cells in a process-specific trajectory often reveals nuanced timing and dynamics of gene expression that bulk assays cannot capture (Qiu et al., 2017; Trapnell et al., 2014). Mapping this trajectory is frequently referred to as ’pseudotime analysis’. We employed pseudotime analysis to arrange single WI-38 cells along a pathway to senescence setting proliferating cells as the trajectory ‘root‘. As expected, early and late PDL cells concentrated at the beginning and end of the pseudotime trajectory respectively (Figure 7A). We next performed differential expression analysis to identify genes that change significantly over pseudotime. We plotted examples of genes changing early (CENPK–an S-phase cell-cycle-regulated gene), midway (SNAI2–a master regulator of EMT), or late (PAPPA–a prominent SASP factor) in pseudotime in Figure 7B.

Figure 7 with 2 supplements see all

Download asset Open asset

Pseudotime (PS) analysis of WI-38 approach to replicative senescence using single-cell RNA-seq.

(A) UMAP projection of single WI-38 cells collected at increasing PDLs (PDL25-red) to (PDL50-blue) and colored by pseudotime (top). (B) Scatterplot of single-cell gene expression across pseudotime with three genes representative of changes that occur early in pseudotime (top), during the transition phase (middle) and in late pseudotime (bottom). Each point is a cell colored by PDL–PDL25 (red) to PDL50 (blue). x-axis is pseudotime, and y-axis is counts for the gene of interest. Black line is a cubic spline. (C) Hierarchical clustering and heatmap of smoothened gene expression trajectories over 60 pseudotime bins (x-axis) of 25 K-median gene expression clusters (y-axis). Value plotted is the scaled (min expression to max expression; 0–1) median expression of all genes in the cluster. Clusters are divided into three (early, transition, late) pseudotime categories. (D) LISA transcription factor enrichment analysis using the 25 clusters from A divided into the same three (early, transition, late) pseudotime categories. Transcription factors and clusters are further divided on vertical axis into putative groupings based on transcription factor functions and GO term enrichment.

Figure 7—source data 1 Gene expression changes with pseudotime.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig7-data1-v3.txt
Download elife-70283-fig7-data1-v3.txt
Figure 7—source data 2 GO term enrichment in pseudotime clusters.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig7-data2-v3.xlsx
Download elife-70283-fig7-data2-v3.xlsx

We next generated gene expression trajectories for the top (by significance p-value < 0.001) 5000 differentially expressed genes across pseudotime (Figure 7—figure supplement 1, Figure 7—source data 1). To identify the temporal relationship between the biological processes and transcription factors that compose replicative senescence in WI-38 cells, we clustered and ordered pseudotime trajectories (Figure 7—figure supplement 2). Broadly, we classified the pseudotime expression pattern as early, transition, or late based on the maximum median value of all constituent genes for each cluster (Figure 7C). We performed LISA transcription factor and GO enrichment analysis on each cluster and used our findings to assign putative functional labels across pseudotime (Figure 7D, Figure 7—figure supplement 2, Figure 7—source data 1 and Figure 7—source data 2).

Early pseudotime is dominated by the transcription factors (E2F) and GO terms associated with cell cycle progression through the G2M and S phases. Moving down the y-axis deeper into pseudotime, we observe that the next primary cluster of transcription factors and functional annotations exhibited widespread enrichment across all of pseudotime. Furthermore, transcription factors enriched in this cluster are involved with basic cellular functions (e.g. euchromatin maintenance, transcription, growth) and are likely representative of normal WI-38 function in G1 phase.

Moving forward into the transition phase of pseudotime, we observed an enrichment of transcription factors regulating higher order chromatin structure (CTCF) and epigenetic silencing (polycomb group complex). Also in the transition region of pseudotime, the earliest enrichments appear for transcription factors related to EMT (SNAI2), TGF-β signaling (SMAD3), YAP1/TEAD1 activity, and the AP1 complex (JUN, FOS, FOSL2). We observed continued enrichment for these transcription factors and processes throughout the rest of the transition phase and into late pseudotime.

Lastly, late in pseudotime, we observed enrichment of transcription factors and functional annotations related to regulation of inflammatory processes (NFKB1, RELA, ERG1, CEBPB) and changes in cellular morphology (Figure 7D and Figure 7—figure supplement 2). These observations are consistent with observations made with bulk RNA-seq that replicative senescence WI-38 cells exhibit transcriptomic features similar to that observed in TGF-β signaling, EMT, and with YAP1/TEAD1 activity. Collectively, these results present a possible order of operations for replicative senescence progression that highlights an initial cessation of active mitotic cycling, followed by an epigenetic shift that precedes a strong EMT/TGF-β signal before segueing into a pro-inflammatory secretory state.

Senescent WI-38 fibroblasts express canonical myofibroblastmarkers and metabolic features

Given the repeated observations linking replicative senescence with EMT, (Figure 1, Figure 3), TGF-β and YAP1/TEAD1 activity (Figure 6, Figure 7), we considered the possibility that these processes are connected through the fibroblast to myofibroblast transition (FMT), a subtype of EMT in which stressed or injured fibroblasts differentiate into myofibroblasts (Phan, 2008; Piersma et al., 2015; Hinz, 2007). Upon receiving cues mediated by injury or stress (e.g. activated TGF-β), fibroblasts can trans-differentiate into myofibroblasts, whose functions as ‘professional repair cells’ include increased proliferation, migration, apoptosis avoidance, cell and tissue contraction, and ECM/collagen deposition to promote tissue repair and wound closure (Hinz and Lagares, 2020; Gibb et al., 2020).

Previous work has demonstrated that there exists mechanistic and functional association between telomerase inhibition, senescence, and myofibroblasts. Senescence is an integral part of the wound healing processes; upon injury resolution, activation of a senescence-like phenotype prevents unchecked collagen secretion and fibrosis by preventing myofibroblast proliferation and earmarking them for subsequent immune clearance (Demaria et al., 2014; Krizhanovsky et al., 2008; Jun and Lau, 2010; Mellone et al., 2016; Razdan et al., 2018; Liu et al., 2002; Liu et al., 2006).

To further explore this proposition, we retrieved canonical myofibroblast marker genes and direct transcriptional targets of the YAP1/TEAD1 complex to determine to what extent these genes are expressed in WI-38 replicative senescence cells at both the RNA and protein level. First, we examined canonical myofibroblast markers (Hinz and Lagares, 2020), collagen produced by myofibroblasts (Zhang et al., 1994), genes upregulated in myofibroblasts derived from idiopathic pulmonary fibrosis (IPF) patients (Strunz et al., 2020), and effectors and targets of TGF-β signaling in two data modalities: bulk RNA-seq and bulk proteomics (Figure 8A).

Figure 8 with 1 supplement see all

Download asset Open asset

YAP regulation of myofibroblast markers and YAP1/TEAD1 targets during replicative senescence.

(A) RNA-seq and proteomics heatmaps of selected genes based on myofibroblast markers and biology. Values plotted are log2 fold change of each PDL vs. PDL 20 for replicative senescence timecourse and vs. no treatment for late PDL (PDL 40) WI-38 cells treated with the YAP inhibitor verteporfin. Proteomics values are median values for n = 3 replicates. Values for individual RNA-seq replicates (n = 3) are shown. (B) Diagram of collagen processing (left) with metabolites in blue and proteins in green. Collagen specific amino acid derivatives are outlined in red. Heatmap of log2 fold changes for metabolites and proteins involved in collagen processing (right). Proteomics values are median values for n = 3 replicates. Metabolomics values are hTERT batch corrected median values for n = 4 replicates. (C) Heatmaps of YAP1 targets collected by Kurppa et al., 2020. for proteomics and RNA-seq during replicative senescence and in late PDL (PDL 40) WI-38 cells treated with the YAP inhibitor verteporfin. Values plotted are same as in A. (D) Heatmap of log2 fold changes relative to PDL 20 or untreated control of leading edge genes from gene set enrichment analysis driving the hallmark EMT signature in verteporfin-treated WI-38 cells and during replicative senescence. Values for individual RNA-seq replicates (n=3) are shown.

Figure 8—source data 1 RNA-seq data for verteporfin treatment.: https://cdn.elifesciences.org/articles/70283/elife-70283-fig8-data1-v3.xlsx
Download elife-70283-fig8-data1-v3.xlsx

For the majority of genes in the curated myofibroblast panel, expression of transcripts and protein increased with PDL. Importantly, expression of smooth-muscle actin (ACTA2), a classic myofibroblast marker, increases strongly in replicative senescence. The expression of follistatin-like protein (FSTL1), also known to be strongly expressed in smooth muscle, shows a similar pattern, as does fibrillin (FBN1). All three of these are associated with smooth muscle and TGF-β family regulation. In addition we observed an increase in both fibrillar and basal lamina collagens in our data at both the RNA and protein levels.

Collagen processing is a multi-step process requiring the coordination of multiple enzymes and metabolites (Lodish and Zipursky, 2001). Review of the collagen synthesis pathway alongside our metabolomic and proteomic data provide further confirmation that WI-38 replicative senescence cells exhibit altered collagen metabolism (Figure 8B). We observed upregulation of multiple pathway enzymes as well as increased abundance of hydroxyproline, a primary constituent amino acid of collagen protein. In addition, we observe a striking depletion of ascorbate (vitamin C) which is required for proline hydroxylation and is an essential vitamin. It is possible that our observations underestimate the collagen production potential of replicative senescent WI-38 cells as they appear limited in terms of collagen production by the amount of supplemented vitamin C (Boyera et al., 1998).

Mellone et al., 2016 recently reported that although senescent fibroblasts share features with myofibroblasts, this resemblance does not extend to fibrogenic ECM components, e.g. collagen (Mellone et al., 2016). However, it is important to note that Mellone et al., 2016 focused on radiation induced senescence rather than replicative senescence, and we similarly observed less or no induction of many of these same genes in our radiation induced senescence experiment. In addition, it is possible that with more time ( > 10 days), fibrogenic ECM expression could initiate in the irradiated WI-38 cells.

Moving forward in our myofibroblast panel, we observed that expression of TGF-β cytokine, TGF-β1, decreased significantly with replicative senescence in both RNA and protein. However, we see robust induction of the TGF-β isotype 2 (TGF-β2) cytokine with replicative senescence (Figure 8A) as TGF-β1 abundance drops, which indicates a switch in TGF-β isotypes with replicative senescence. It has been shown TGF-β2 is a more potent inducer of the endothelial to mesenchymal transition (EndMT) in vitro compared to TGF-β1 and TGF-β3 in human microvascular endothelial cells, and TGF-β2 may be playing a similar role here in inducing FMT and replicative senescence in WI-38 cells (Sabbineni et al., 2018). The distinct functional roles of different TGF-β isotypes are largely unknown, although both are known to activate the SMAD transcription factors. On the basis of our data, it seems likely that the TGF-β paralog relevant here might not be TGF-β1, but TGF-β2.

Expression of YAP1/TEAD1 targets during replicative senescence

Next, we retrieved a gene set of YAP1/TEAD1 targets assembled by Kurrpa et al. from five separate studies (Kurppa et al., 2020). Taking the intersection of the five YAP1/TEAD1 gene target lists, we only kept genes present in at least two of the studies and plotted the remainder in Figure 8C across the two data types. We arranged the YAP1/TEAD1 targets by the two predominant expression patterns–decrease with replicative senescence (top) and increase with replicative senescence (bottom). As with myofibroblast markers, the data from both modalities is largely concordant.

Interestingly, we noted a striking bifurcation whereby YAP1 targets tend toward either strong down-regulation or strong up-regulation. The down regulated partition is heavily enriched for classic cell cycle regulators such as TOP2A, CDC20, BIRC5, and CDK9 suggesting that this dichotomy in YAP1 activity is heavily influenced by the cell cycle and consistent with recent work (Kim et al., 2019b; de Sousa et al., 2018).

Excitingly, TGF-β2 was one of the two YAP1 targets found in four out of five collected YAP1/TEAD1 data sets, supporting the inference that the TGF-β2 and not TGF-β1 is the relevant paralog and potentially regulated by YAP1/TEAD1 in the replicative senescence context. In addition, another YAP1 target, thrombospondin-1 (THBS1), is also upregulated at both the protein and transcript level with RS; THBS1 is known to act as an activator of TGF-β signaling, and specifically TGF-β2, through proteolytic cleavage of latent TGF-β (Ribeiro et al., 1999).

Inhibiting the YAP1/TEAD1 interaction in WI-38 cells suppresses expression of YAP1/TEAD1 targets and the replicative senescence gene signature

In all, the data reveal that replicative senescence WI-38 cells share multiple defining transcriptomic and proteomic features with myofibroblasts. Furthermore, we found that a subset of YAP1/TEAD1 targets are both induced with replicative senescence, and are principal components of TGF-β signaling (TGF-β2, THBS1). Thus, the YAP1/TEAD1 complex may be acting in convergence with TGF-β signaling with increasing PDL to enact a myofibroblast-like state that we recognize as replicative senescence.

Recently, Mascharak et al., 2021 showed verteporfin treatment of wounds in mice alleviated the fibrotic state during wound healing, reverted the profibrotic transcriptional program, and reduced the myofibroblast population. Verteporfin is also known to be an inhibitor to the formation of the YAP1/TEAD1 complex (Wang et al., 2016). To test for a direct role for YAP1/TEAD1 in transcriptional regulation of YAP target genes, myofibroblast marker expression and the EMT signature during replicative senescence transition, we treated late passage WI-38 PDL 40 cells with verteporfin or DMSO as an untreated control. We utilized RNA-seq to capture the gene expression changes on cells treated with verteporfin. First, to confirm the verteporfin treatment inhibited expression of canonical YAP target genes in WI-38 cells, we checked the expression of CTGF, CYR61, and TGFB2 transcripts (Figure 8—figure supplement 1; Zhao et al., 2008; zhang2009tead; chen2001angiogenic). We found that treatment with 10 µM verteporfin for 2 hours reduced expression of these three YAP target genes by 70% relative to the no treatment control (Figure 8—figure supplement 1, Figure 8—source data 1). We expanded the analysis to the myofibroblast markers and YAP targets used in Figure 8 and observed that verteporfin treatment reduces the gene expression of most myofibroblast markers and collagen and many YAP1 target genes that were shown to increase with replicative senescence (Figure 8A and C).

We applied GSEA analysis using the MSigDB Hallmark annotation to understand what pathways are affected by verteporfin treatment of cells using the RNA-seq data from the verteporfin treated samples relative to non-treated control. Consistent with our hypothesis that YAP1/TEAD1 plays a role in regulating the FMT-like transition of cells undergoing replicative senescence, treatment of late passage cells with verteporfin resulted in negative enrichment of pathways relevant to FMT and TGF-β signaling (Figure 8—figure supplement 1B). We also took genes driving the EMT annotation enrichment during replicative senescence to see if inhibition of YAP1/TEAD1 complex formation aids in reducing EMT pathway gene expression. Compared to the no treatment control, the majority of genes induced during replicative senescence that drive the EMT enrichment exhibit a striking decrease in expression after two hours of treatment. These data suggest that the use of verteporfin to inhibit the YAP1/TEAD1 interaction can impede the FMT-like transition observed in replicative senescent cells (Figure 8—figure supplement 1D).

To gain a global understanding of how inhibiting the YAP1/TEAD1 interaction affects the transcriptomic landscape of late passage WI-38 cells, we compared the gene expression signature of replicative senescence cells to that of verteporfin-treated cells and found that verteporfin treatment globally reverses many of the gene expression changes occurring with replicative senescence (Figure 8—figure supplement 1). This results suggest that YAP1/TEAD1 transcriptional activity is a crucial regulator of a large swath of the replicative senescence transcriptome.

Discussion

The study of replicative senescence in human tissue culture has proven to be an informative model for learning how genetic and environmental factors impact cellular senescence. However, the field has not fully taken advantage of the ‘omics’ revolution. Rekindled by the advent of senolytics, interest in the molecular underpinnings of senescence has burgeoned in recent years as researchers seek to design therapeutic strategies for ablating senescent cells (Amor et al., 2020; Elmore et al., 2018; Wagner and Gil, 2020; Aghajanian et al., 2019). However, many such studies span only one or two systematic data modalities. Here, to fully leverage the power of recent advances in high-dimensional profiling, we revisit the original Hayflick limit in WI-38 lung fibroblasts cells with a battery of assays including RNA-seq, ATAC-seq, scRNA-seq, proteomics, and metabolomics in an effort to capture the defining features of replicative senescence at every step of the central dogma and beyond. Our results are summarized graphically in Figure 9.

Figure 9

Download asset Open asset

Replicative senescence fibroblast to myofibroblast transition (FMT) model.

We divided Replicative senescence progression into three major categories and summarized our results across all data modalities focusing on features in common with myofibroblasts.

One important feature of our study is the number and type of control conditions. From our pilot studies, it was clear that studying replicative senescence in isolation would preclude our ability to know what features of change were specific to replicative senescence versus high cell density and/or DNA-damage-induced growth arrest. Thus, we included cellular density as an alternate method for arresting growth and exposure to ionizing radiation as an alternate methodology for senescence induction and DNA damage response. Lastly, all of our modalities were paired with samples collected from an hTERT immortalized cell line grown in parallel. Importantly, very little changed across all modalities in the hTERT cell line consistent with the changes we highlight in nine as specific to replicative senescence.

The kinetics and precise timing of senescence onset have been obscured by low temporal resolution and ensemble measurements that cannot differentiate between global shifts in gene expression versus changing proportions of senescent cells (Passos et al., 2007; Xu et al., 2013; Victorelli and Passos, 2017; Smith and Whitney, 1980; Tang et al., 2019; Wiley et al., 2017; Nassrally et al., 2019; Smith and Whitney, 1980). Here, with a combination of high time resolution and single-cell RNA-seq, we provide evidence that the early manifestation of the senescent gene expression reflects gradual changes on a per cell basis rather than changing cell proportions. In effect, individual cells ‘show their age’” with increasing PDL long before permanently exiting the cell cycle and transiting fully into the senescent state. The implications of this conclusion extend to organismal aging. For example, the percentage of senescent cells calculated from aging organisms varies greatly depending on the marker/phenotype used (Ogrodnik, 2021). The reported disparities could be explained in part by the use of early versus end stage markers. Likewise, it is possible the reported increase in fibroblast heterogeneity and altered functionality with age is a direct result of cells slowly moving along a spectrum towards replicative senescence (Shin et al., 2020; Mahmoudi et al., 2019). Importantly, the gradual progression suggests that cells need not reach the endpoint to elicit a phenotype. For instance, proliferative fibroblasts isolated from IPF patients exhibited multiple senescent features and phenotypes in addition to accelerated senescence progression (Yanai et al., 2015). Lastly, this phenomenon is not constrained to fibroblasts as we observe that the salient regulatory features of replicative senescence extend to cell types as distant as astrocytes (Figure 6D).

In our data, the pattern of gene expression annotated to EMT as a unique feature of replicative senescence that consistently presents early and robustly at both the RNA, protein, and single cell level. Given the fundamental nature of the EMT transition with respect to cellular function (development, fibrosis, and wound healing), it is not surprising that this hallmark tracks with multiple proteins and gene sets suggesting drastic metabolic rewiring. In our metabolic data, we highlight shifts in carbon and fatty acid utilization that have been reported previously as hallmark metabolic features of EMT. These metabolic changes demonstrate that our observations represent an authentic change in cellular state as opposed to a superficial uptick in a few EMT-related genes.

The data presented above across multiple data modalities to provide a clear connection between senescent cells and myofibroblasts supported by independent observations at the level of DNA, RNA, protein, transcription factor activity and metabolism (Figure 9). In light of these, replicative senescence resembles a specialized subtype of EMT specific to the trans-differentiation of fibroblasts into myofibroblasts (FMT) in response to wound healing (Lombardi et al., 2019; Gibb et al., 2020; Hinz and Lagares, 2020). We hypothesize that during fibrotic disease states and/or age, fibroblasts migrate to sites of micro-injuries. As these proliferating fibroblasts become replicatively aged, they are triggered (by DNA damage or other insults) to rewire their metabolism to induce FMT via active epigenetic reorganization (NNMT-SAM/NAD sink). It is important to note here that the in vitro DNA damage here arises primarily from telomere erosion documented by the observation that our hTERT control cultures do not exhibit the same changes. However, genotoxic stress in vivo may originate from a variety of endogenous and environmental sources e.g. reactive oxygen species, replication stress, chemical exposure etc.

Following increases in DNA accessiblity, expression of newly opened TEAD1/YAP1/SMAD target genes cement FMT transition by promoting fibrosis and a myofibroblast-like, ECM-secreting state. This model is supported by the synergy between TEAD1 motifs and NAD domains we report and reconciles conflicting reports that YAP1/TEAD1 inhibition can both prevent and promote senescence (Fu et al., 2019; Kurppa et al., 2020; Xie et al., 2013; Jia et al., 2018). Basically, the functional consequences of YAP1/TEAD1 inhibition will depend on the epigenetic organization of the cells used.

Further metabolic changes (hexoasmine/collagen synthesis and fatty acid oxidation) then support the new pro-fibrotic state. Finally, end point senescent cells reinforce the senescent state and contribute to neighboring cell progression towards senescence via secretion of inflammatory factors and SASP.

We think that this transition is distinct from classic FMT as the endpoint cells are not proliferative, but instead bear striking resemblance to lingering senescent myofibroblasts that can persist long after wound repair is complete (Hinz and Lagares, 2020). Rather than a privileged or unique state, perhaps replicative senescence is better categorized as a DNA damage mediated path to a potentially common stress-induced endpoint.

The complex interplay between metabolism and epigenetic regulation preclude easy determination of a causal factor in translating DNA damage into the replicative senescence/EMT program in WI-38 cells. Does EMT regulate metabolism or vice versa? Here we present compelling evidence on the side of metabolism. We observe both an early and sharp rise in NNMT expression and activity in addition to a global increase in heterochromatin accessibility. These results are consistent with NNMT’s reported role as global epigenetic regulator through its methylation sink activity (Komatsu et al., 2018; Eckert et al., 2019; Pissios, 2017; Ulanovskaya et al., 2013).

Importantly, the observed shifts in repressed chromatin induced by NNMT are functional and may play a central role in fibroblast biology and stress response in multiple contexts. First, Eckert et al., 2019 recently demonstrated that NNMT activity and the resulting heterochromatin reorganization initiate the expression program of cancer-associated fibroblasts (CAFs) associated with oncogenic stroma in vivo (Eckert et al., 2019). Similar to senescent fibroblasts and myofibroblasts, the defining features of CAFs are increased cytokine production, metabolic rewiring, and ECM alteration and production (Sahai et al., 2020). Second, NNMT is one the most upregulated genes in a TGF-β-mediated in vitro FMT conversion in WI-38 cells (Walker et al., 2019). Finally, the changes in silenced chromatin we observed during replicative senescence overlap with induction of gene expression driving the senescent phenotype (Figure 5D). In fact, pseudotime analysis argues that chromatin reorganization may precedes the FMT induction as evidenced by the loss of polycomb activity that appears prior to enrichment of the EMT and YAP1/TEAD1 transcription factors enrichments (Figure 7D).

The path to replicative senescence process has many in vivo parallels with implications for aging and pathogenesis. After observing induction in WI-38 cells of IPF myofibroblast markers (Figure 8A), we expanded our literature search and found a striking overlap (TGF-β signaling, YAP1 activity, and EMT) between our data and scRNA-seq expression profiles from alveolar epithelial cells collected from IPF patients (Xu et al., 2016). These studies also report a large induction of TGF-β2 relative to TGF-β1 which is consistent with our findings and highlights a clear connection between in vitro replicative senescence and an in vivo disease state.

Another example of such a connection arises from our LISA analysis: one of the top enriched transcription factors, TCF21, has been implicated in atherosclerotic disease progression. Wirka et al., 2019 found that TCF21 promotes the transition of vascular smooth muscle cells into a novel fibroblast-like cell type they dub ‘fibromyocytes’ owing to their possession of both fibroblast and myocyte phenotypes in atherosclerotic lesions in both mice and humans.

In general, the FMT hypothesis provides a conceptual framework and integrative model for linking the multi-modal senescent phenotypes we observed to multiple human age-related diseases. Given the observation that fibrosis and senescence markers correlate with increasing age in multiple tissues, it is possible that FMT might be a widespread phenomenon underlying many age-related pathologies (Yousefzadeh et al., 2020; Idda et al., 2020). Future work harnessing multi-modal single cell technology coupled with relevant in vivo models will aid greatly in determining the exact order of events and physiological import.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers	Additional information
Cell Line (H. sapiens)	WI-38 fibroblasts	Coriell	AG06814-N
Cell Line (H. sapiens)	WI-38 hTERT	This paper	immortalized WI-38 human fibroblasts, methods
Transfected construct (H. sapiens)	pCDH-CMV-hTERT-EF1a-puro	This paper		Lentiviral plasmid, methods
Commercial assay, kit	Senescence-β-Galactosidase	Cell Signaling Tech.	9,860
Commercial assay, kit	MycoAlert Mycoplasma Detection	Lonza	LT07-218
Commercial assay, kit	Direct-zol RNA Miniprep Plus	Zymo Research	R2072
Commercial assay, kit	Chromium Single Cell 3’ v2	10 x Genomics	120,237
Commercial assay, kit	Chromium Single Cell A Chip	10 x Genomics	1000009
Commercial assay, kit	Tagment DNA Enzyme and Buffer	illumina	20034197
Commercial assay, kit	Clean-and-Concentrator-5	Zymo Research	D4014
Commercial assay, kit	NEBNext High-Fidelity 2 X PCR	NEB	M0541L
Commercial assay, kit	TruSeq Stranded mRNA Library	illumina	20020595
Commercial assay	Bioanalyzer High Sensitivity DNA	Agilent	5067–4626
Commercial assay, kit	Pierce BCA Protein Assay	Thermo Fisher	23,227
antibody	anti-human p16 antibody (Mouse monoclonal)	BD Biosciences	RRID:AB_395229	(1:250)
antibody	anti-human p21 antibody (Mouse monoclonal)	BD Biosciences	RRID:AB_396414	(1:250)
Commercial assay, kit	Pippin Prep 2% 100–600 bp	Sage Science	CDF2010
Chemical compound, drug	Verteporfin	RD Systems	1243926
Other	Zorbax Extend C18 column	Aglient	759700–902
Other	SeQuant ZIC-pHILIC column	EMD Millipore	150,460
Software, algorithm	R (v4.0.3 and 3.6.2)	r-project.org/	RRID:SCR_001905
Software, algorithm	Salmon (v 0.8.2)	combine-lab.github.io/salmon/	RRID:SCR_017036
Software, algorithm	DESeq2 (v1.30.1)	bioconductor	RRID:SCR_015687
Software, algorithm	sva package(v3.38.0)	bioconductor	RRID:SCR_012836
Software, algorithm	fgsea 1.16.0	bioconductor	RRID:SCR_020938
Software, algorithm	CellRanger 3.0	10 x Genomics	RRID:SCR_017344
Software, algorithm	SCTransform (v 0.3.2)	Satijalab.org/seurat
Software, algorithm	Seurat (v4.0.1.9005)	satijalab.org/seurat	RRID:SCR_007322
Software, algorithm	monocle3 (v1.0.0)	cole-trapnell-lab.github.io/monocle3	RRID:SCR_018685
Software, algorithm	bowtie2 (v2.3.4.1)	bowtie-bio.sourceforge.net bowtie-bio.sourceforge.net	RRID:SCR_016368
Software, algorithm	samtools (v1.2)	http://ww.htslib.org/	RRID:SCR_002105
Software, algorithm	Picard (v2.6.4)	broadinstitute.github.io/picard	RRID:SCR_006525
Software, algorithm	macs2 (v2 2.1.2)	hbctraining.github.io/main/	RRID:SCR_013291
Software, algorithm	GenomicRanges (v1.42.0)	bioconductor	RRID:SCR_000025
Software, algorithm	cutadapt (v2.4)	https://github.com/marcelm/cutadapt	RRID:SCR_011841
Software, algorithm	bcl2fastq (v2.20)	Illumina	RRID:SCR_015058
Software, algorithm	regioneR v1.22.0	bioconductor
Software, algorithm	ATACseqQC v1.14.4	bioconductor
Software, algorithm	LIMMA v3.46.0	bioconductor	RRID:SCR_010943
Software, algorithm	Qvalue v2.26.0	combine-lab.github.io/salmonQvalue	RRID:SCR_001073
Software, algorithm	LISA v1	lisa.cistrome.org

Share this article

Cite this article

Expression dynamics of replicative senescence (RS), radiation induced senescence (RIS), increasing cell density (CD) and hTERT WI-38 cells.

Figure 1—source data 1

Figure 1—source data 2

Figure 1—source data 3

Figure 1—source data 4

Cell cycle exit and distribution on approach to replicative senescence (RS) does not explain gradual increase in the replicative senescence transcriptome.

Proteomic and metabolomic changes during replicative senescence.

Figure 3—source data 1

Figure 3—source data 2

Nicotinamide n-methyltransferase (NNMT) links nicotinamide adenine dinucleotide (NAD) and methionine metabolism during replicative senescence.

Increased accessibility within heterochromatin and nucleolar associated domains (NADs) is a dominant feature of the replicative senescence (RS) epigenome.

Figure 5—source data 1

Figure 5—source data 2

Figure 5—source data 3

Figure 5—source data 4

Figure 5—source data 5

Figure 5—source data 6

Master transcriptional regulators of replicative senescence.

Pseudotime (PS) analysis of WI-38 approach to replicative senescence using single-cell RNA-seq.

Figure 7—source data 1

Figure 7—source data 2

YAP regulation of myofibroblast markers and YAP1/TEAD1 targets during replicative senescence.

Figure 8—source data 1

Replicative senescence fibroblast to myofibroblast transition (FMT) model.

Author details

Michelle Chan

Contribution

Competing interests

Han Yuan

Contribution

Contributed equally with

Competing interests

Ilya Soifer

Contribution

Contributed equally with

Competing interests

Tobias M Maile

Contribution

Competing interests

Rebecca Y Wang

Contribution

Competing interests

Andrea Ireland

Contribution

Competing interests

Jonathon J O'Brien

Contribution

Competing interests

Jérôme Goudeau

Contribution

Competing interests

Leanne JG Chan

Contribution

Competing interests

Twaritha Vijay

Contribution

Competing interests

Adam Freund

Contribution

Competing interests

Cynthia Kenyon

Contribution

Competing interests

Bryson D Bennett

Contribution

Competing interests

Fiona E McAllister

Contribution

Competing interests

David R Kelley

Contribution

Competing interests

Margaret Roy

Contribution

Competing interests

Robert L Cohen

Contribution

Competing interests