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
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).

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
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
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
Figure 3—source data 2

Metabolomics quantification and differential analysis.

https://cdn.elifesciences.org/articles/70283/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).

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
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).

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 5Figure 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
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
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.

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
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.

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.

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 resourceDesignationSource or referenceIdentifiersAdditional information
Cell Line (H. sapiens)WI-38 fibroblastsCoriellAG06814-N
Cell Line (H. sapiens)WI-38 hTERTThis paperimmortalized WI-38 human fibroblasts, methods
Transfected construct (H. sapiens)pCDH-CMV-hTERT-EF1a-puroThis paperLentiviral plasmid, methods
Commercial assay, kitSenescence-β-GalactosidaseCell Signaling Tech.9,860
Commercial assay, kitMycoAlert Mycoplasma DetectionLonzaLT07-218
Commercial assay, kitDirect-zol RNA Miniprep PlusZymo ResearchR2072
Commercial assay, kitChromium Single Cell 3’ v210 x Genomics120,237
Commercial assay, kitChromium Single Cell A Chip10 x Genomics1000009
Commercial assay, kitTagment DNA Enzyme and Bufferillumina20034197
Commercial assay, kitClean-and-Concentrator-5Zymo ResearchD4014
Commercial assay, kitNEBNext High-Fidelity 2 X PCRNEBM0541L
Commercial assay, kitTruSeq Stranded mRNA Libraryillumina20020595
Commercial assayBioanalyzer High Sensitivity DNAAgilent5067–4626
Commercial assay, kitPierce BCA Protein AssayThermo Fisher23,227
antibodyanti-human p16 antibody (Mouse monoclonal)BD BiosciencesRRID:AB_395229(1:250)
antibodyanti-human p21 antibody (Mouse monoclonal)BD BiosciencesRRID:AB_396414(1:250)
Commercial assay, kitPippin Prep 2% 100–600 bpSage ScienceCDF2010
Chemical compound, drugVerteporfinRD Systems1243926
OtherZorbax Extend C18 columnAglient759700–902
OtherSeQuant ZIC-pHILIC columnEMD Millipore150,460
Software, algorithmR (v4.0.3 and 3.6.2)r-project.org/RRID:SCR_001905
Software, algorithmSalmon (v 0.8.2)combine-lab.github.io/salmon/RRID:SCR_017036
Software, algorithmDESeq2 (v1.30.1)bioconductorRRID:SCR_015687
Software, algorithmsva package(v3.38.0)bioconductorRRID:SCR_012836
Software, algorithmfgsea 1.16.0bioconductorRRID:SCR_020938
Software, algorithmCellRanger 3.010 x GenomicsRRID:SCR_017344
Software, algorithmSCTransform (v 0.3.2)Satijalab.org/seurat
Software, algorithmSeurat (v4.0.1.9005)satijalab.org/seuratRRID:SCR_007322
Software, algorithmmonocle3 (v1.0.0)cole-trapnell-lab.github.io/monocle3RRID:SCR_018685
Software, algorithmbowtie2 (v2.3.4.1)bowtie-bio.sourceforge.netbowtie-bio.sourceforge.netRRID:SCR_016368
Software, algorithmsamtools (v1.2)http://ww.htslib.org/RRID:SCR_002105
Software, algorithmPicard (v2.6.4)broadinstitute.github.io/picardRRID:SCR_006525
Software, algorithmmacs2 (v2 2.1.2)hbctraining.github.io/main/RRID:SCR_013291
Software, algorithmGenomicRanges (v1.42.0)bioconductorRRID:SCR_000025
Software, algorithmcutadapt (v2.4)https://github.com/marcelm/cutadaptRRID:SCR_011841
Software, algorithmbcl2fastq (v2.20)IlluminaRRID:SCR_015058
Software, algorithmregioneR v1.22.0bioconductor
Software, algorithmATACseqQC v1.14.4bioconductor
Software, algorithmLIMMA v3.46.0bioconductorRRID:SCR_010943
Software, algorithmQvalue v2.26.0combine-lab.github.io/salmonQvalueRRID:SCR_001073
Software, algorithmLISA v1lisa.cistrome.org

Cell culture

WI-38 maintenance and subculturing for replicative senescence time course-WT and hTERT

Request a detailed protocol

WI-38 cells were obtained from the Coriell Institute (AG06814-N) at PDL 15. The cells were grown in Dulbecco’s modified Eagle’s medium (DMEM, Gibco, 11885084) supplemented with 10% dialyzed fetal bovine serum (FBS, Sigma, F0392) and maintained in an incubator set to 37° C, 5% CO2, and 20% O2. 0.3–0.5 E6 WI-38 cells were seeded and maintained on 10 cm collagen coated plates (Corning, 354450) and split when cells reached about 70% confluence. For WI-38 cells PDL 15 tp PDL 40, it took about 4 days to reach 70% confluence. WI-38s > PDL 40 were slower growing, and therefore split every 5–7 days instead; media was replenished every 3–4 days (supplementary methods). Cells were passaged by washing the cell monolayer with PBS, followed by incubating cells with TrypLE Express (Gibco, 12604013) for 5 min. Media was added to neutralize the TrypLE and the cell suspension was collected into conical tubes. The cell suspension was spun at 200 x g for 5 min and the supernatant was aspirated. The resulting cell pellet was resuspended in media and cells were counted with viability measurements on a ViCell XR Cell Analyzer (Beckman Coulter). Cell numbers taken from the ViCell were used to calculate population doublings (PDLs). The following formula was used to calculate PDLs:

PDL=log2(numberofcellsharvested)log2(numberofcellsseeded)

All cells collected for assays were sampled in triplicate, at 2.5–3 days after seeding, targeting 60–70% confluence to avoid the confounding effects of confluence. Confluence levels were determined with a phase contrast EVOS microscope (ThermoFisher). Slower growing cells (PDL > 40) were sampled at 3.5–4 days after seeding. Cells were tested every month for mycoplasma contamination using the MycoAlert Mycoplasma Detection Kit (Lonza). A smaller subset of three PDL timepoints (PDL 45, PDL 55, PDL 56) were generated in a secondary time course to ensure ‘deep’ senescence (2–4 months after cell cycle cessation) was maintained as previously reported (De Cecco et al., 2019).

WI-38 hTERT cells were maintained using the same conditions as WT WI-38s, targeting 50–70% confluence for cell splitting and sampling alongside WT cells.

Sampling proliferating WI-38 cells at different cell densities confounds assessment of cell cycle state and senescence

Request a detailed protocol

A limited pilot (performed prior to RS time course) experiment was conducted to determine optimal sampling conditions. Using gene expression as a readout for cellular state, we performed RNA-seq on PDL 20, PDL 29, and PDL 39 WI-38 cells. We plated an equal number of cells in triplicate for the three PDLs and cells were collected for RNA extraction after 4 days of growth. We looked for any significant changes in transcript levels between the three samples.

We observed two notable trends. First, we found stronger expression of genes associated with cell cycle and mitosis in the older PDL 39 cells compared to PDL 20 and PDL 29 cells. This was unexpected based on previous observations of slower growth and less active mitosis in cells closer to senescence (Hayflick, 1965). Second, the cells at PDL 20 and PDL 29 expressed genes associated with hypoxic stress (Figure 1—figure supplement 5A).

These data, in tandem with the observation that older cells grow at a slower rate, suggest that the highly proliferative younger cells (PDL 20 and 29) grow to higher cell density faster and therefore encounter earlier in real time, the various normal mitosis-inhibiting stresses, for example hypoxia, contact inhibition, and starvation that trigger cell cycle exit (Dayan et al., 2009).

To understand the effects of growth rates and cell density on gene expression, we performed a cell density study (CD) wherein young PDL cells were sampled at increasing levels of cell density to measure the effects of cell density and proliferation on gene expression (Figure 1A). We observed that based on the seeding density we employed (to minimize the stress of sparse plating) the sampling window we had originally targeted (4 days after seeding) was located in a volatile region with respect to proliferation rate and hypoxic stress using Cyclin B1 (CCNB1) and VegfA (VEGFA) expression as proxies for each, respectively (Figure 1—figure supplement 5B).

Accordingly, we modified the sampling window by choosing a time interval (2.5–3 days) prior to the onset of VEGF-A induction and CCNB1 repression (Figure 1—figure supplement 5B). For older, slower growing cells in the plateauing region of the replicative senescence process (PDL >40) we shifted the sampling window to 3.5–4 days based on visual inspection and cell counts. Ultimately in the actual extended time course, we tracked cell density with a phase contrast EVOS microscope to ensure cellular confluence was similar across all time points during sampling (Figure 2—figure supplement 2).

WI-38 irradiation time course

Request a detailed protocol

WI-38 cells at PDL 20 were plated in DMEM supplemented with 10% FBS with 50,000 cells per well in a 6-well collagen coated plate. Cells were allowed to settle and adhere to plates. All cells were adhered to the TC plate by 2 hr and were subsequently treated with 10 Gy of X-rays (Faxitron CellRad). Cells were sampled between 1 day and 9 days for transcriptome profiling. Each time point was sampled in triplicate. RNA was extracted according to methods below.

WI-38 cell density time course

Request a detailed protocol

0.3 E6 WI-38 cells at PDL 23 were seeded and grown on 10 cm collagen-coated plates in DMEM supplemented with 10% FBS. Cells were grown and sampled intermittently between 1 day and 10 days for transcriptome profiling. Samples at each time point were taken in triplicate. RNA was extracted according to methods below.

SA-β GAL staining

Request a detailed protocol

Cells were stained for senescence associated beta-gal using the Senescence β-Galactosidase Staining Kit (Cell Signaling, 9860) by following the manufacturer’s published protocols exactly.

WI-38 hTERT cells and lentiviral transduction

Request a detailed protocol

293T cells were transfected with a lentiviral target plasmid expressing hTERT (pCDH-CMV-hTERT-EF1a-puro) and lentiviral packaging constructs overnight; 48 hr later, viral supernatant was collected. WI-38 cells were transduced with viral supernatant in the presence of 5 µg/mL polybrene and selected for 7 days with 1 µg/ml of puromycin.

WI-38 + verteporfin treatment

Request a detailed protocol

A total of 200,000 WI-38 cells at PDL 38 were seeded onto each well of a six-well collagen-coated plate and grown with 2 ml of DMEM supplemented with 10% FBS. After 4 days of growth, cells were treated with 5 or 10 µM of verteporfin (RD Systems 5305/10) for 2 hr or DMSO. RNA was extracted after 2 hr of treatment according to methods described below.

Bulk RNA-Seq methods

RNA collection and library preparation

Request a detailed protocol

Total RNA was extracted from cells using the Direct-zol RNA Miniprep Plus kit (Zymo Research R2072) for all bulk RNA-seq time course experiments and verteporfin treatment experiment. Manufacturer’s protocol was followed exactly and in-column DNAase digestion was performed. RNA quality score and concentration was measured using the Fragment Analyzer (Agilent 5200) with the Fragment Analyzer Standard Sense RNA kit (Agilent Technologies DNF-471–0500). All samples required to have a RIN score of >7 for processing. RNA sequencing libraries were prepared as directed using TruSeq Stranded mRNA Library Prep Kit (llumina 20020595), with 1000 ng of input material. Samples were amplified for 12 cycles of PCR with TruSeq RNA CD Index Plate (Illumina) and pooled. 3 nM libraries were loaded across four lanes on the HiSeq 4000 (Illumina).

Read processing and quantification

Request a detailed protocol

Reads generated from the Illumina HiSeq 4000 were demultiplexed with bcl2fastq (version=2.20) based on the barcode sequence of each sample. Average read depth across samples was 50 milllion paired-end reads. Reads were pseudo-aligned and quantified using Salmon (version=0.8.2) by deploying the mapping based mode using a Salmon command ‘index’ with default parameters based on 10 X genomics hg38 transcriptome annotations optimized for single cell RNA-seq (‘refdata-cellranger-GRCh38-3.0.0’, cellRanger version 3) to ensure accurate comparison between bulk and single-cell RNA-seq (Patro et al., 2017). Annotations can be obtained running ‘wget https://cf.10xgenomics.com/supp/cell-exp/refdata-cellranger-GRCh38-3.0.0.tar.gz’.

Differential expression analysis

Request a detailed protocol

DESeq2 (version=1.30.1) was used for differential analysis of the RNA-seq data Love et al., 2014. Wald test was used to estimate fold change and significance using the model:~ time + batch, where time is a numeric variable representing the fraction of time course complete and batch is a categorical variable used only with replicative senescence with the addition of the ‘deep’ senescence time points. DEGs were defined as having FDR adjusted p-values < 0.01.

Batch correction and hierarchical clustering

Request a detailed protocol

For clustering and visualization, we corrected the raw replicative senescence count table by batch before converting to transcripts per million (TPM) with a + 1 pseudocount and combining with radiation induced senescence and cell density samples using Combat-seq (sva R package version=3.38.0) (Zhang et al., 2020). Significant genes from each condition were concatenated to generate a universe of significant change genes used for Figure 1. Each sample was then converted to log2 fold change vs. the mean of the initial time point.

Gene set enrichment analysis

Request a detailed protocol

For gene set enrichment analysis, we downloaded the MSigDB Hallmark gene sets from misgDB website (version 7, gene symbols) Liberzon et al., 2015. For each time course (replicative senescence, radiation-induced senescence, and cell density) we we ranked all genes by the log2 fold change across time generated by DESeq2. Then GSEA was performed on the ranked gene set using the R package ‘fgsea’ version 1.16.0 which is an implementation of GSEA in R (Subramanian et al., 2005; Korotkevich et al., 2019; Mootha et al., 2003). By default, GSEA tests for enrichment of each gene set in each condition in both directions. We report the -log10 Benjamini-Hochberg corrected p-value in Figure 1 and use the normalized enrichment score (NES) to assign direction of the change.

Single-cell 3’ RNASeq methods

Cell collection

Request a detailed protocol

At each time point, singlet hTERT controls and experimental samples were processed with Chromium Single Cell 3’ RNAseq kit V2 (10 x Genomics 120237) through cDNA amplification cleanup, where they were frozen at –20° C. Once all time points were converted to cDNA, the frozen cDNAs were thawed and batched for library construction. The following modifications were made to the process: Reverse transcription reactions were brought up to volume with DMEM +10% FBS instead of water, each emulsion targeted 3000 cells (5200 cells loaded), and cDNA was amplified 12 cycles, with 12 cycles of index PCR. There was no hTERT control for time point 3, and one of the replicates for PDL 25 dropped out during library construction. Remaining samples were pooled at equimolar concentration and sequenced on a HiSeq4000 with the standard 26,8,0,98 run base pairs per read configuration.

Single-cell data processing, normalization, scoring, clustering, and DEG analysis

Request a detailed protocol

The raw single-cell reads were demultiplexed by sample using bcl2fastq. Alignment, cell barcode demultiplexing, transcript quantification and sample merging were carried out with CellRanger 3.0 using the hg38 CellRanger 3.0 gene annotation (‘refdata-cellranger-GRCh38-3.0.0’, cellRanger) with default options and parameters. Filtered cell by gene matrices were normalized using Seurat (version = 4.0.1.9005) and SCTransform (version = 0.3.2) (Hafemeister and Satija, 2019; Stuart et al., 2019). Dimension reductionality was carried out with PCA (n = 50). Clusters were defined using the louvain algorithm resolution = and cells were visualized with UMAP projection McInnes et al., 2018; Becht et al., 2018; Blondel et al., 2008. For cell cycle scoring and phase determination as well as senescence scoring, we applied the Seurat implementation of the scoring function as previously described (Nestorowa et al., 2016).

DESeq2 was used for differential gene expression analysis to identify significantly changing genes within individual clusters as a function of increasing PDL (PDL by mitotic phase) (Figure 3D) using a an input matrix of gene counts by cells per PDL by mitotic phase. Gene counts were summarized across cells for each mitotic phase (PDL 50 removed). To expedite computation, we restricted analysis to highest expressed 8000 genes. For each PDL by cell cycle phase grouping, we required that >15 representative cells must exist to be considered in analysis. To visualize differentially expressed genes, we converted single cell counts to counts per million ( + 1 pseudocount) and averaged across PDL by cell cycle phase and calculated the log2 fold change at each PDL by cell cycle phase against the earliest PDL for that cell cycle phase using a concatenated list of all significantly changing genes generated by DESeq2 from each cell cycle phase.

Pseudotime analysis

Request a detailed protocol

For pseudotime analysis of the single-cell data, we used the R package monocle3 (version=1.0.0) which implements PCA, Leiden clustering, and UMAP prior to partitioning and trajectory analysis. We first isolated WT cells by removing hTERT cells and those not belonging to the primary grouping. To focus trajectory analysis on the replicative senescence progression, we used only the top genes found to significantly change with replicative senescence in bulk RNA-seq for PCA (padj <0.01) and log2 fold change >3 or log2 fold change < –3). In applying the Monocle trajectory analysis: (learn_graph(cds,close_loop = FALSE,learn_graph_control = list(minimal_branch_len = 35, geodesic_distance_ratio = 0.5, euclidean_distance_ratio=1)), we designated cycling cells at the opposite end of UMAP project from the senescent cells as the start point manually (Trapnell et al., 2014; Qiu et al., 2017; Cao et al., 2019; Traag et al., 2019; Levine et al., 2015; McInnes et al., 2018). After establishing a trajectory, we then employed the Monocle3 ‘graph test’ function to isolate genes that significantly change as a function of pseudotime. Pseudotime estimation was output from monocle3 using the learn graph function for building a trajectory. Smoothed pseudotime trajectories use for Figure 7C were calculated for significantly changing genes (with pseudotime) by binning cells across pseudotime into 60 bins and using a cubic spline to estimate expression at each bin. For each gene, the smoothed trajectory was set from 0 (minimal) to (1) maximal expression. Genes were organized with K-median clustering (k=25) using cosine similarity. For visualization the median expression value for each pseudotime bin and cluster was calculated. Genes from each cluster were fed into LISA TF analysis (below).

ATAC-SEQ Methods

ATAC-Seq library preparation and sequencing

Request a detailed protocol

Freshly harvested cells were used for all reactions. Briefly, the cell monolayer was washed with PBS, trypsinized with TrypLE Express (Gibco 12604013), resuspended in media, and cells were pelleted. Cells were counted and 100,000 cells were used in each reaction. Cell lysis, DNA transposition, and library construction was followed from the Omni-ATAC protocol (Corces et al., 2017). Libraries were amplified for 13 total cycles. Sample purification and size selection were performed on the Pippin high throughput size selection platform (Sage Science) using 2% agarose cassettes to isolate nucleosome free fragments that were <300 base pairs in library fragments and thus corresponding to <150 base pairs fragments after accounting for sequencing adaptors. Quality of ATAC-seq libraries were assessed with the Agilent Bioanalyzer 2,100 with DNA high sensitivity chips (Agilent Technologies 5067–4626). Libraries were sequenced on the HiSeq4000 with paired end sequencing using 2 × 150 bp reads (Illumina).

ATAC-Seq data processing, peak calling, and differential accessibility

Request a detailed protocol

We first trimmed the raw fastq files with cutadapt (version=2.4) to remove standard Nextera PE adapters:

  • adapt -a file:$ADAPTER -A file:$ADAPTER -o $SAMPLE.R1.fq -p $SAMPLE.R2.fq –pair-filter=any –minimum-length=30 $R1 $R2

  • Then we aligned with bowtie2 (version version 2.3.4.1) to align the trimmed reads to hg38: bowtie2 -x $INDEX –1 $SAMPLE.R1.fq –2 $SAMPLE.R2.fq –no-mixed –no-discordant -X 1000

After alignment, we used samtools (version=1.2) flags (-f 0 x02 and -q2 0) to filter for only properly paired and high quality reads. PCR duplicates are removed using picard (version=2.6.4) MarkDuplicates. Finally for each bam file, we adjusted the reads ends by Tn5 offset ( + 4 on+ strand, –5 on -strand).

ATAC-seq QC (mitochondrial percent and transcription start site score) and alignment metrics were generated with R package ATACseqQC (version=1.14.4) and multiQC respectively (version=1.11) (Ewels et al., 2016) using TxDb.Hsapiens.UCSC.hg38.knownGene (version = 3.10.0) annotations for transcription start site score calculation.

For peak calling, we created a condition specific peak atlas by pooling all replicates in a specific condition and applied macs2 (version = 2 2.1.2) for peak calling on the pooled bam file with options (-g hs -p 1e-1 –nomodel –shift –37 –extsize 73) (Zhang et al., 2008). In addition, we performed peak calling on each individual replicates as well. Then we performed irreproducible discovery rate analysis on each condition specific peak atlas for each pair of replicates and filter for peaks that are reproducible in at least two replicates (IDR threshold of 0.05). A single accessibility atlas is created by merging condition-specific peak atlas across all conditions (Figure 5—source data 3). Peaks were assigned to nearest gene if it is within 50 kb, otherwise it is annotated as intergenic.

Read count was performed using countOverlaps function from R package GenomicRanges (version = 1.42.0) (Lawrence et al., 2013). We performed quantile normalization of the count matrix using normalize.quantiles function of R package preprocessCore (version=1.52.1).

Limma-voom (Limma version = 3.46.0) was used for differential accessibility analysis (Ritchie et al., 2015). Fold change and fdr adjusted p-values were estimated using moderated t-test statistic based on the model: time+ condition + time:condition. We performed separate tests for time point one versus each of the other time points.

ATAC-Seq chromatin state and NAD/LAD analysis

Request a detailed protocol

To quantify ATAC-seq reads in chromatin states, we retrieved hg38 ENCODE IMR-90 chromatin state labels from (Ernst and Kellis, 2015; Figure 5—source data 1). Next, we quantified coverage for each instance of all 25 chromatin states genome wide using countOverlaps function from R package GenomicRanges (Lawrence et al., 2013). We performed quantile normalization of the chromatin state count matrix using normalize.quantiles function of R package preprocessCore (Figure 5A).

ATAC-seq peaks were assigned to chromatin states using the findOverlapsOfPeaks function from the R package GenomicRanges (Lawrence et al., 2013). We used only significantly (P < 0.001) changing peaks(with senescence) from LIMMA-voom analysis. To simplify the overlap of these two sets of genomic intervals, we took only peaks that fell within (‘inside’) or encompassed a chromatin state annotation (‘inside feature’). Peaks encompassing more than one chromatin state interval were discarded (Figure 5A and B).

For analysis of NAD and LAD domain overlap we collected IMR-90 NAD labels from Dillinger et al., 2017 and IMR-90 LAD labels from Sadaie et al., 2013 (Figure 5—source data 5). For calculating overlap Z-score between genomic intervals sets (e.g. NADs/LADs vs. gene annotation) (TxDb.Hsapiens.UCSC.hg38.knownGene, version = 3.10.0) for Figure S18C we used the overlapPermTest function from the R package regioneR (version = 1.22.0) Gel et al., 2016.

For testing overlaps between NAD/LAD domains and significantly changing genes, we used the overlapPermTest function from the regioneR package using the top induced genes from each bulk RNA-seq condition (RS,RIS, and CD) against a universe of all genes sampled from the same expression distribution to control NAD and LAD bias across different levels of expression. For gene annotations, we used TxDb.Hsapiens.UCSC.hg38.knownGene (version = 3.10.0).

Metabolomics

Extraction of water-soluble metabolites from lung fibroblast cell culture

Request a detailed protocol

Twenty-four hours before metabolite extraction, the medium was aspirated, cells were washed with unconditioned medium and then the medium was replaced. For metabolite extraction, cells were washed once in 37°C warm PBS-buffer immediately followed by the addition of 3.5 mL of freezer-cooled (-20) LC-MS grade 80:20 MeOH/H2O (Sigma Aldrich). The plates were then held at -20 for 2 hr, then harvested with a sterile cell scraper while at –20°C and transferred to -20 cold 5 mL centrifuge tubes (Eppendorf Lo-Bind). After centrifuging the cell-extracts in a 4 centrifuge for 5 min at 2000 x g, the supernatants were transferred into new cold centrifuge tubes and dried under nitrogen at 4. Dried extracts were stored at -20.

LC-MS/MS analysis of cell culture extracts

Request a detailed protocol

Dried supernatants were resuspended in 200 µL of water containing 1 µg/mL of deuterated lysine and deuterated phenylalanine and 250 ng/mL of deuterated succinate (Sigma Aldrich) as internal standards. For negative ion mode, the resuspended samples were diluted 1:4 in water, for positive ion mode, they were diluted 1:4 in acetonitrile. Samples were then centrifuged at 18,000 x g for 5 min, the supernatant was moved to HPLC vials and 5 µL was injected for analysis by LC-MS on Vanquish UPLCs (Thermo Scientific) coupled to Q Exactive Plus mass spectrometers (Thermo Scientific).

For analysis in negative ion mode, separation of compounds was achieved by reverse-phase liquid chromatography on a Zorbax Extend C18 column (150 × 2.1 mm, 1.8 µm particle size [Agilent 759700–902]). Mobile phase A was 10 mM Tributylamine and 15 mM Acetic Acid in 97:3 water:methanol at pH 4.95 and mobile phase B was methanol. Prior to injection, the column was equilibrated in 0% B for 4 minutes. The gradient eluted isocratically in 0% B for 2.5 min, increased to 20% B over 2.5 min, held at 20% B for 2.5 min, increased to 55% B over 5.5 min, increased to 95% B over 2.5 min, maintained at 95% B for 3 min, then decreased to 0% B over 0.5 min, where it was held for 3 min, resulting in a total run time of 26 min. The LC separation was carried out at 30 column temperature using a 200 µL/min flow rate. For MS-analysis, parameters on MS1 were set to 70,000 resolution with an AGC target of 1e6 at a maximum IT of 100ms. The scan range was 70–1050 m/z. MS2 parameters were set to 17,500 resolution at loop count 6, an AGC target of 1e5 at a maximum IT of 50ms, an isolation window of 1 m/z and an underfill ratio of 1%. Dynamic exclusion was set at 20 s with an apex trigger from 3 to 12 s. Stepped collision energies were set to 20, 50, and 100% NCE.

For analysis in positive ion mode, compounds were separated via hydrophilic liquid interaction chromatography (HILIC), using a SeQuant ZIC-pHILIC column (150 × 2.1 mm, 5 µm particle size) Millipore (150460). Mobile phase A consisted of 20 mM ammonium bicarbonate at pH 9.2 in H2O, and mobile phase B was acetonitrile. Prior to injection, the column was equilibrated for 6 minutes in 80% B. The gradient then decreased to 20% B over 20 min, then to 15% B over 2 min, returned to 80% B over 0.5 min and held there for 1.5 min for a total run time of 30 min. The column temperature was 35 with a flow rate of 150 µL/min. For MS-analysis, the MS1 parameters were as described for negative ion mode except the AGC target was 3e6. MS2 parameters were the same with following exceptions: dynamic exclusion was set to 25 s with an apex trigger from 3 to 10 s. Stepped collision energies were set to 20, 40% and 80% NCE.

Metabolomics batch correction and differential analysis

Request a detailed protocol

Protein lysates derived from same plates used for metabolomics collection were assayed for total protein concentration with Pierce BCA Protein Assay Kit (ThermoFisher 23225) for normalizing raw metabolite values. The median protein concentration for each PDL and PDL.ctrl was calculated and divided by the median protein concentration across all samples for either cell line (WT and hTERT) to derive a protein concentration normalization factor. The raw metabolite values for each PDL and hTERT PDL.ctrl were multiplied by the protein concentration normalization factor for that strain and sample. Protein normalized metabolite values were then converted to log2 fold change values using the first time point as a reference (Figure 3—source data 2). To remove sample day batch effects, we took the protein normalized metabolite values and divide each WT metabolite value by the corresponding metabolite value in the temporally paired hTERT sample (e.g. PDL 25 values were divided by PDL 25.ctrl values) and converted to log2 fold change with first sample as reference.

Finally, the normalized log2 ratios were fit to a linear model to test for a linear trend across the time course. For each metabolite, we fit the model yit=β0+β1Tt+ϵit, where i = (1,...,3) indexes the replicates for each time point and t=(1,..,4) indexes the cell passage vector т = (33, 37, 46, 50).

The model was fit with the lm() function in R Version 4.1.2. p-values from the two-sided hypothesis test H0:β1 = 0 were FDR adjusted (across all proteins) using the qvalue function from the qvalue R package (version =2.26.0). If present, asterisks in figures indicate that the FDR adjusted p-value from this test was <.05.

Proteomics methods

Materials and sample preparation, extraction and digestion

Request a detailed protocol

LC-MS grade organic solvents, water, and tandem mass tag (TMT) isobaric reagents were purchased from Thermo Fisher Scientific (Waltham, MA). Trypsin was ordered from Promega Corporation (Madison, WI) and Lys-C from Wako Chemicals USA (Richmond, VA). Sep-Pak C18 cartridges were from Waters Corporation (Milford, MA). Unless otherwise stated, all other chemicals were purchased from Sigma-Aldrich (St. Louis, MO).

At time of sample collection, cells were trypsinized off the monolayer with TrypLE Express and media was used to neutralize the reaction. Cells were pelleted at 200 x g for 5 min and the supernatant was aspirated. Cell pellets were flash frozen in liquid nitrogen and stored at –80°C until all samples were collected and ready for proteomics analysis. Cell pellets were resuspended in 450 µL of lysis buffer (75 mM NaCl, 3% SDS, 50 mM HEPES, pH 8.5) and lysed by passage through a BD PrecisionGlide 21-gauge needle (20 X). The lysate was sonicated for 5 min and then centrifuged (5 min, 16,000 x g) to remove cellular debris and the supernatant was collected.

Proteins were reduced with 5 mM dithiothreitol (DTT) for 30 min at 56 with shaking. Proteins were then alkylated with 15 mM iodoacetamide (IAM) for 30 min at room temperature (RT) in the dark, and excess IAM was quenched with 5 mM DTT for 20 min at RT in the dark. Protein purification was accomplished using a methanol-chloroform precipitation. Briefly, 800 µL methanol, 200 µL chloroform and 600 µL water were sequentially added to 200 µL of cell lysate, with 5-s intervals of vortexing between each addition. The samples were centrifuged for 30 min (16,000 x g at 4°C) to induce phase separation and both the top and bottom layers were removed. The precipitated protein pellet was washed with 600 µL methanol, vortexed briefly, then centrifuged for 10 min (16,000 x g at 4 °C). The methanol layer was removed and protein pellets were dried at RT for 10 minutes. Protein pellets were resuspended in digestion buffer (8 M urea, 50 mM HEPES, pH 8.5). The urea concentration was diluted to 4 M, then proteins were digested with Lys-C overnight (10 ng/µL, 25°C, 16 hr). The urea concentration was further diluted to 1 M and samples were digested with trypsin (5 ng/µL) for 6 hr at 37 ◦C.

Following digestion, peptides were acidified with trifluoroacetic acid (TFA) to a final concentration of 0.5% TFA. Peptides were desalted using Sep-Pak C18 solid-phase extraction (SPE) columns and samples were eluted sequentially, first with 40% acetonitrile (ACN)/0.5% acetic acid and then 80% ACN/0.5% acetic acid. Eluted peptides were dried in a CentriVap Benchtop Vacuum Concentrator (Labconco, Kansas City, MO) running at 30 °C. Peptide concentrations were measured using the Pierce BCA Protein Assay Kit, then 50 µg aliquots of each samples were dried in the CentriVap for further processing.

Tandem mass tag (TMT) labeling

Request a detailed protocol

Dried peptides were resuspended in 50 µL 200 mM HEPES/30% anhydrous ACN, then 200 µg of each TMT tag was added to 50 µg peptides. TMT 131 c was used as the ‘bridge sample’ while the other tags (126, 127 n, 127 c, 128 n, 128 c, 129 n, 129 c, 130 n, 130 c, 131 n) were used to label the individual samples. The TMT reaction was incubated for 1 hour at room temperature with gentle shaking, then quenched with 11 µL 5% hydroxylamine/200 mM HEPES. All samples were acidified to a final concentration of 0.5% TFA. A small amount (4 µL) of each labeled sample was combined and desalted using StageTips to check TMT ratios and labeling efficiency. The TMT-labeled samples were then combined at a 1:1:1:1:1:1:1:1:1:1:1 peptide ratio into 11-plex samples (Rappsilber et al., 2007). The combined samples were desalted using Sep-Pak C18 cartridges and dried under vacuum.

High pH reversed-phase (HPRP) fractionation

Request a detailed protocol

The pooled TMT-labeled peptides were fractionated using high pH reversed-phase liquid chromatography on an Agilent 1260 Infinity HPLC equipped with a diode array detector set at 215, 220, and 254 nm wavelengths (Agilent Technologies, Santa Clara, CA). Peptides were separated on an Agilent ZORBAX Extend-C18 column (4.6 mm x 250 mm, 5 µm particle size) running at 500 µl/min at 25°C. Peptides were eluted with a gradient with initial starting condition of 100% buffer A (5% ACN, 10 mM ammonium bicarbonate) and 0% buffer B (95% ACN, 10 mM ammonium bicarbonate). Buffer B was increased to 35% over 60 min, then ramped up to 100% B in 6 s where it was held for 5 min. Buffer B was then decreased to 0% over 6 s and held for 10 min to re-equilibrate the column to original conditions. The samples were fractionated into 96 fractions, then pooled into 12 fractions as previously described (Huttlin et al., 2010). The fractions were dried under vacuum and resuspended in 5% ACN/5% formic acid (FA) for LC-MS/MS analysis.

Proteomics data acquisition and analysis

Request a detailed protocol

LC-MS/MS Data Acquisition All samples were analyzed by an Orbitrap Fusion Lumos Tribrid mass spectrometer coupled to an EASY-nLC 1200 (Thermo Fisher Scientific). Peptides were separated using a microcapillary column (100 µm x 250 mm long, filled in-house with Maccel C18 AQ resin, 1.8 µm, 120 Ã; Sepax Technologies, Newark, DE) operating at 60 °C with a flow rate of 300 nL/min. Peptides were eluted into the mass spectrometer using a 180 min method, with acetonitrile increasing from 6% to 30% over a 165 min linear gradient in 0.125% formic acid. Mass spectrometry data was collected in data-dependent acquisition (DDA) mode. A high-resolution MS1 scan (500–1200 m/z range, 60,000 resolution, AGC 5 × 105, 100ms max. injection time, RF for S-lens 30) was collected in the Orbitrap, and the top 10 precursors were selected for MS2 and MS3 analysis. Ions were isolated using a 0.5 m/z window for MS2 spectra. The MS2 scan was performed in the quadrupole ion trap (CID, AGC 1 × 104, 30% normalized collision energy, 35ms max. injection time) and the MS3 scan was analyzed in the Orbitrap (HCD, 60,000 resolution, max. AGC 5 × 104, 250ms max. injection time, 50% normalized collision energy). The max. cycle time was set at 5 s. For TMT reporter ion quantification, up to six fragment ions from each MS2 spectra were selected for MS3 analysis using synchronous precursor selection (SPS).

Proteomics data analysis

Request a detailed protocol

The Ting et al., 2011 software pipeline and methods developed in the the Haas and Gygi labs was used to process all proteomics data (Ting et al., 2011). Raw files were converted to mzXML files and searched against a composite human UniProt database containing forward and reverse sequences using the Sequest algorithm. MS/MS spectra were matched with fully tryptic peptides from this composite dataset using a precursor ion tolerance of 20 ppm and a product ion tolerance of 0.6 Da. TMT modification of peptide N-termini and lysine residues ( + 229.162932 Da) and carbamidomethylation of cysteine residues ( + 57.02146 Da) were set as static modifications. Oxidation of methionine residues ( + 15.99492 Da) was set as a variable modification. Peptide spectral matches were filtered to a 1% false discovery rate (FDR) using linear discriminant analysis (LDA) as previously described (Ting et al., 2011). Non-unique peptides that matched to multiple proteins were assigned to proteins that contained the largest number of matched redundant peptides sequences using the principle of Occam’s razor (Ting et al., 2011).

Quantification of TMT reporter ion intensities was performed by extracting the most intense ion within a 0.003 m/z window at the predicted m/z value for each reporter ion. TMT spectra were used for quantification when the sum of the signal-to-noise for all the reporter ions was greater than 200 and the isolation specificity was greater than 0.75.

Peptide level intensities were converted to log2 ratios by dividing each scan by the intensity in the bridge channel. Relative protein abundance, for each sample, was estimated with the posterior mean from a previously described Bayesian model using methods as described by O’Brien et al., 2018 with code located at https://github.com/ColtoCaro/compositionalMS (O’Brien, 2020).

Finally, the relative abundance estimates were fit to a linear model to test for a linear trend across the passages as described by Gaun et al., 2021. For each protein, we fit the model yit=β0+β1Tt+ϵit, where i=(1,...,3)indexes the three replicates for each time point and t=(1,...,7)indexes the cell passage vector T=(18,25,32,33,37,46,50).

The model was fit with the lm() function in R Version 3.6.0. p-Values from the two-sided hypothesis test H0:β1=0 were FDR adjusted (across all proteins) using the p.adjust() function. If present, asterisks in all figures indicate that the FDR adjusted p-value from this test was <0.01.

Regulatory analysis

LISA TF analysis

Request a detailed protocol

Identification of putative regulators of the gene expression changes observed in the bulk RNA-seq experiments and the pseudotime analysis was carried out using (Qin et al., 2020); a recent algorithm built to leverage the vast amount of protein-DNA interactions catalogued via ENCODE. We used the GUI hosted at http://lisa.cistrome.org/ to enter lists of genes derived from K-medians pseudotime clustering clustering. The online GUI has a max input of 500 genes. Accordingly, for clusters containing >500 genes, each gene was ranked by correlation with the cluster median and the top 500 genes were used. LISA output consists of a ranked file of transcription factors and chromatin modifiers with enrichment p-values associated with specific ENCODE experiments. Given that the calculated p-values were derived from gene sets of different sizes and from different number of supporting experiments, we used the ranks as the input parameter for the LISA analysis shown in Figures 6D and 7D. A universe of top factors across clusters was compiled by concatenating the top 5 TFs from each individual LISA output. The rank matrix was then centered across clusters to identify TFs with the highest rank in specific clusters.

ATAC-Seq peak motif enrichment

Request a detailed protocol

We built a binary peak-by-motif matrix where each row is a binary vector representing the presence of 405 motifs we used from CisBP (Catalog of Inferred Sequence Binding Preferences) which is a freely available online database of transcription factor (TF) binding specificities: http://cisbp.ccbr.utoronto.ca (Weirauch et al., 2014). We only used motifs for TFs with measurable RNA-seq in WI-38 cells. In order to characterize transcription factor activity changes during senescence, we implemented both a gene-centric and a gene-indepedent approach.

In the gene-centric approach, we tested enrichment of each of the 405 motifs in peaks that are (1) associated with a gene set of interest (within 50 kb), and (2) significantly associated with senescence (adjusted p-value < 0.001 in PDL 50 vs. PDL 20) as compared to all the other peaks using binomial test.

Alternatively, we trained a ridge logistic regression model using the binary motif matrix as features to distinguish peaks of significantly increased accessibility during senescence from peaks of significantly reduced accessibility in PDL 50 vs. PDL 20. A model trained on two thirds of the data distinguishes the two sets of peaks with AUC=0.67 on held-out peaks. We then trained 10 independent models using all the data to evaluate the coefficients to identify features (motifs) that are most predictive of senescence.

Data availability

Sequencing data have been deposited in GEO under accession code GSE175533. Proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository under accession code EBI PRIDE Code for processing and analyzing data modalities have been deposited at https://github.com/dghendrickson/hayflick, (copy archived at swh:1:rev:e32e36ceb9aed82f15b794b28bbbcee51bf5e9e8). Source data for figures and analysis have been deposited at https://github.com/dghendrickson/hayflick and/or uploaded as source data files.

The following data sets were generated
    1. Chan M
    2. Hendrickson DG
    (2021) NCBI Gene Expression Omnibus
    ID GSE175533. Revisiting the Hayflick Limit: Insights from an Integrated Analysis of Changing Transcripts, Proteins, Metabolites, and Chromatin.
    1. Chan M
    2. Hendrickson DG
    (2021) PRIDE
    ID EBI PRIDE. Revisiting the Hayflick Limit: Insights from an Integrated Analysis of Changing Transcripts, Proteins, Metabolites, and Chromatin.
The following previously published data sets were used
    1. Gregory BD
    2. Gosai S
    3. Crowe EP
    4. Torres CA
    (2014) NCBI Gene Expression Omnibus
    ID GSE58910. Transcriptome analysis of oxdative-stress induced senescence in human astrocytes.
    1. Németh A
    2. Dillinger S
    (2017) NCBI Gene Expression Omnibus
    ID GSE78043. Mapping of nucleolus-associated chromosomal domains during replicative senescence.
    1. Salama RA
    2. Narita M
    (2013) NCBI Gene Expression Omnibus
    ID GSE49341. Redistribution of the Lamin B1 genomic binding profile affects spatial rearrangement of heterochromatic domains and gene expression during senescence.

References

    1. Lodish H
    2. Zipursky SL
    (2001)
    Molecular cell biology
    Biochemistry and Molecular Biology Education: A Bimonthly Publication of the International Union of Biochemistry and Molecular Biology 29:126–133.
    1. Wang C
    2. Zhu X
    3. Feng W
    4. Yu Y
    5. Jeong K
    6. Guo W
    7. Lu Y
    8. Mills GB
    (2016)
    Verteporfin inhibits YAP function through up-regulating 14-3-3σ sequestering YAP in the cytoplasm
    American Journal of Cancer Research 6:27–37.
    1. Zhang K
    2. Rekhter MD
    3. Gordon D
    4. Phan SH
    (1994)
    Myofibroblasts and their role in lung collagen gene expression during pulmonary fibrosis. A combined immunohistochemical and in situ hybridization study
    The American Journal of Pathology 145:114–125.

Decision letter

  1. Weiwei Dang
    Reviewing Editor; Baylor College of Medicine, United States
  2. Kevin Struhl
    Senior Editor; Harvard Medical School, United States
  3. Lynne Cox
    Reviewer
  4. Payel Sen
    Reviewer; National Institute on Aging, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Revisiting the Hayflick Limit: Insights from an Integrated Analysis of Changing Transcripts, Proteins, Metabolites and Chromatin" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Matt Kaeberlein as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Lynne Cox (Reviewer #2); Payel Sen (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) All three reviewers have raised various concerns about robustness of analyses and statistics for the multi-omics data. They recognize that significant amounts of these analyses will need to be redone and some additional analyses were recommended. However, they all agree that these revisions are essential. Please address these concerns as much as possible.

2) Reviewers #2 and #3 pointed out the necessity to provide some biological validations for novel findings from these multi-omics analyses. This will require additional experiments but will greatly improve this manuscript.

3) All three reviewers have identified the lack of description for certain methods and techniques, mislabeling, typographical errors, etc. Please address all these issues.

4) Reviewer #2 pointed out that the proteomics data should be made publicly available at one of the major proteomics data repositories. Reviewer #1 also pointed out that in-house software codes and scripts for data analyses should also be made publicly available.

Reviewer #1 (Recommendations for the authors):

1. There are major concerns about the statistics in the analysis of the dataset.

a. For the differential expression analysis (4.2.3), the authors state: "DEGs were defined as having p-values < 0.001." This is an incorrect method of filtering DEGs as it does not correct for multiple hypothesis testing. This is all the more surprising because DESeq2, used by the authors, natively provides FDR corrected p-values. The data should be re-analyzed after proper multiple hypothesis correction control (native DEseq2 FDR, Bonferroni, or other).

b. Potentially related to the absence of FDR thresholding, the number of DEGs identified-"8969 genes change with increasing PDL … 11652 genes vary with increasing cell density"-seemed, to this reviewer, really high, and likely the result of improper FDR correction. It will be interesting to see how those numbers change with the multiple testing correction.

c. Based on the methods, the other analyses should also always incorporate a multiple testing correction, i.e. differential ATACseq peaks (4.4.3 p < 0.001, no FDR correction), metabolomics (4.5 no information on differential analysis) and proteomics (4.6.5 no information on differential analysis).

d. When carrying out differential analyses, the authors should always avoid filtering on fold change. They do this for instance in section 2.4 "we focused on compounds exhibiting PDL-dependent changes in abundance > log2 0.5.". Fold-change filtering leads to extremely poor FDR control and should not be used post hoc – the only acceptable method to integrate fold change filtering includes a new statistical framework which takes into account fold change for FDR calculations (see PMC2654802). All results should be presented with this criterion removed or re-calculated with a method that takes fold-change into account.

e. Based on figure legends (Figure 1F, 3A, etc) and methods (4.2.5), FDR correction for multiple testing was also not performed for GSEA results. Please recompute everything including the correction, and only report corrected p-values.

f. There are several key issues with the ATAC-seq analysis (4.4.2), which will require a ground-up reanalysis of the data. Quantile normalization of count data and subsequent use of limma for differential accessibility analysis of ATAC-seq (4.4.2, end) is potential problematic. Indeed, limma is not appropriate for count data (it should only be used for continuous data such as microarray) unless the authors mean that they used limma-voom, which is a different algorithm. Please correct this in the revised analyses. In cases where differential depth cannot be corrected for appropriately, the authors should instead use downsampling of mapped reads to match the lowest depth sample. This will remove sequencing depth biases and allow proper count based methods, such as DEseq2, to perform properly.

2. Additional methodological information is needed for long-term reproducibility of analyses.

a. For reproducibility of code and analyses, all analytical code for this study should be deposited in a repository such as github or made available as a Supplementary file.

b. Please include all version numbers for all used software (e.g. R, Trimmomatic, Salmon, LISA, etc.) packages and R packages (e.g. DESeq2, fgsea, etc.) command parameters where relevant (such as for Trimmomatic, if not available elsewhere). The same should be applied to annotation databases, and if a version number doesn't exist, date of access should be provided.

c. All reagents (e.g. DMEM) should include suppliers and catalogue numbers in order to remove any ambiguity as to which product was used.

d. Salmon is a pseudoaligning method that uses a transcriptome reference, not a genome. The methods say that the index was derived from hg38, but how and using which information and software is not clear. In general, the authors should make sure that all key analytical steps are explicitly documented in the methods for long term reproducibility of the study.

3. It is completely unclear which data/plots were generated using the Illumina mRNA kit and which were generated using the Total RNA Kit. Please specify in text and legend.

Reviewer #2 (Recommendations for the authors):

The -omics aspects of the paper are generally well conducted and provide either confirmatory data on pre-existing studies of senescence or add new data (especially on the key FMT pathway and transcription factors involved). Log2FC of +/- 0.5 is somewhat lower than usual cut-offs, so justification is needed where that is selected. No mention is made of quality control steps in RNAseq analysis, but I assume these were conducted? Use of in-house programs for proteomics analysis does mean that others will not be able to conduct the same analyses, though availability of datasets should allow for validation of findings on other data analysis platforms – have these been publicly deposited as links provided are only to RNAseq datasets? From the protocol provided, it is also not clear that the methodology for protein extraction would release chromatin-bound factors – inclusion of a potent nuclease e.g. Benzonase, would ensure that tightly bound proteins are not lost in the spin steps to remove 'cellular debris'.

Figure 2E: there appears to be an aberrant pattern for cells at PDL33 that is present in all parts – is 33 an outlier? What is this cluster?

The cell biology aspects are much less well described making it impossible for others to reproduce the work – in particular, cell seeding density is critical but not even mentioned for longitudinal culture, and Figure 1D images suggest significant cell crowding. Along the senescence trajectory, primary fibroblasts grown in monolayers increase in surface area (and to some extent in volume) and populations at later PDLs undergo contact inhibition at lower cell densities than those at earlier PDLs so it may be the case that some of the 'gradual' changes (especially those associated with cell cycle changes) reflect simply quiescence form contact inhibition rather than steps to senescence. Methodology to count cells and calculate PDL is needed as different labs may conduct this differently. Daily microscopy inspection of cultures to assess appropriate time points for subculturing is essential – simply choosing to split at 4 or 7 days according to PDL does not account for the morphological changes that occur across the replicative lifecourse. It would therefore be helpful to determine to what degree these findings are cell type specific and which changes represent general cell senescence, especially as multiple RNAseq datasets exist on lung and skin fibroblast senescence with various mode of senescence induction. The Methods mentions 'deep senescence' but it is not clear in the main text what datasets are derived from this. Images of cells at different points along the pathway to senescence would be reassuring, and orthogonal validation of senescence markers e.g. p16, p21 by qRT-PCR is needed. It would be extremely helpful if the RNA seq and protein changes reported were cross validated by conventional techniques.

Although studies are conducted on the influence of cell density on -omics changes over a short time course of 10 days, the effect of cell cycle exit in senescence would be better accounted for by comparing with genuinely quiescent cell populations e.g. serum starvation of the hTERT cells eg Figure 3A

Time points used for comparisons are poorly labelled – are hTERT samples taken at daily intervals or at the same time points in long term culture as the RS samples? The text states that the aim is to control for the effects of long-term cell culture, but this is not apparent from later figures. Figure 1D implies that cell population at the same PDL are used for RS and hTERT groups (were cells transduced at PDL20 to allow analysis at PDL25 allowing for drug selection steps?) yet later on time points are simply given numbers e.g. 1,7 in Figure 2B and 2,4,5,6,7 in Figure 3A. It is not even stated if these are days, weeks or sample numbers.

A comparison between early PDL WT WI38 with hTERT WI38 should be conducted (the authors have the relevant data sets) to ensure that the hTERT expression simply leads to cell immortality without impacting on gene expression patterns – this is critical since telomerase has roles other than at the telomere, including in gene expression control, RNA splicing and mitochondrial metabolism. hTERT transduction is very poorly described in the methods section so that others will not be able to reproduce the experiments: "appropriate target plasmid and packaging constructs" are completely uninformative, as is 'selected. with a selection drug". Specific details are needed.

Experiments on radiation damage to induce senescence assay relatively early time points – cell cycle exit is an early event in the DNA damage response but not indicative of senescence per se. True damage-induced senescence onset should be monitored at later time points post-damage (eg from 14 days) to avoid conflating the acute DDR with senescence-specific gene expression patterns. This could account for the large discrepancy between RS and RIS datasets here.

The oxidative phosphorylation data are interesting and consistent with known changes on senescence (though the effect on ROS generation should be considered more). However, the citation in this context of papers that refer to lipid peroxidation – a state of lipid damage – rather than physiological energy generation through β oxidation of fats. Please take care not to conflate energy generating β oxidation and damaging oxidation in the citations.

The NNMT data are fascinating and really provocative – did the authors consider including bisulfite sequencing to assess DNA methylation state in conjunction with the ATAC-seq data? The finding of CTCF upregulation is also interesting in the context of the LAD/NAD data – would 3C/Hi-C give clearer data on overall chromatin configuration? (It would be intersting to see if TADs change significantly, which would be predicted from the data here).

The manuscript has not been fully proof-read which is somewhat frustrating for the reader e.g. the first few references don't even have volume or page numbers, there is incorrect use of upper/lower case in some places (e.g. amp instead of AMP) and there are various typos throughout. Methods section has gaps/omissions e.g. 4.4.3 refers to Figure ??B, 4.5.1 refers to "(vendor?)" and 4.7.2 refers to "(Cite CIS-BP)". Tenses should be consistent.

Most worryingly, very little effort has been made to ensure the Supplementary figures are informative – figure legends are lacking, labels of figure parts are missing e.g. time points/PDLs on Figure S2, Figure S3; the interpretation that EMT (or rather FMT) is a feature of sesncence depends to some extent on the data in Figure S2 but these appear to show that there is no change of the progression from early proliferation to senescence (though hard to tell as PDLs not given)? statements such as "S phase and G2M cells were isolated" – this is wholly misleading as the cells were not isolated physically according to cell cycle stage (eg FACS), but the data were processed post-hoc to determine which cell cycle stage they most closely fit with according to gene expression patterns. Figure S8 legend states cells at PDL20, figure labels cells at PDL25; time points/cell ages are missing for Figure S9- and cannot state mitochondrial functions are already simply by looking at protein levels – safer to say that changes are consistent with altered mitochondrial metabolism. Again, PDL20 is referred to yet throughout the paper, the earliest PDL analysis is at 25. Figure S10 is uninterpretable as no X axis is given, and time points/PDLs are also missing in Figure S11B. Figure S12D – what do the boxes represent? The lack of labels and suitable legend make it impossible to interpret. Why analyse PDL 46 here, as cells are not yet senescent? Is the whole of chromosome 22 shown in part E? Please provide Mb scale – and state what platform was used to view (e.g. IGV?).Figure S14 – GO terms would be a better way of showing the processes than the word clouds used. Figure S15 lacks a legend – myofibroblast markers (from lit) – I presume from literature? Which are the fibrillar collagens (the collagen number matters).

When describing what samples the data shown are compared against, the terminology "the first sample" needs more information – is this PDL25 for RS? But this 'first sample' appears later in metabolomic analysis eg Figure 3E (hence fold changes will not be comparable with proteomics).

Overall this could be a really great paper, bringing together a range of potent -omics techniques to study replicative senescence at high resolution cross time. Most of the paper is well written and provides a clear and sensible discussion of the work carried out, the data and implications of the results. It is somewhat dense (which may be inevitable with multi-omics studies) and could benefit from more accessible explanation of some of the analyses performed (eg readers will be more familiar with tSNE than UMAP – a brief explanation would help here). However, the paper is let down in parts by sloppy presentation, lack of critical experimental detail (meaning the work cannot be reproduced by others), poor design (or one hopes simply poor description) of comparators, an apparent lack of understanding of the nuances of senescent cell behaviour in culture, and by lack of statical validation of the simpler experiments (no n values noted in some cases – e.g. Figure 1D, Figure 5A, C, Figure S3 etc etc). Confirmatory experiments to cross validate top hits from the -omics studies are also required. Mining the existing proteomics data for PTMs should also be possible and may add some value, though not essential here.

Reviewer #3 (Recommendations for the authors):

1. There are some technical points about ATAC-seq that need to be clarified. (a) global chromatin accessibility may be higher in senescent cells. Could the authors share the bioanalyzer profiles of the ATAC-seq libraries from the different timepoints? Are they similar? (b) senescent cells typically produce far more mitochondrial reads. Could the authors confirm that sufficient depth of sequencing was achieved in their ATAC-seq runs? These data including sequencing statistics, peak coordinates, fragment distribution graphs etc. should be included in the Supplementary Section.

2. The ATAC-seq data might benefit from a RepEnrich analysis to identify peaks in repeat elements that are known to be heterochromatinized in proliferating cells and derepressed in senescence. This, in my opinion, is a better way to look at constitutive heterochromatin desilencing compared to NAD/LAD overlap and is also independent of NAD/LAD calls in a different cell line (in this case IMR-90).

3. What are the authors thoughts about the decrease in chromatin accessibility at promoters, enhancers, and gene bodies in their ATAC-seq data (Figure 5B)? Does this imply a global shut down of transcription? Is it reflected in the transcriptomic analyses?

4. While the authors perform an in-depth analyses of transcription factors most likely to drive RS, they haven't really shown the direct binding of these factors to chromatin in RS (and not hTERT, RIS or CD cells) and/or the change in their expression in senescence. Additionally, no functional experiments are performed by knocking down or overexpressing TEAD1 and investigating its effect on senescence. If possible, it would be nice to include this data (i.e., a western blot, genomic binding of TEAD1 and some functional experiments).

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

Author response

Essential revisions:

1) All three reviewers have raised various concerns about robustness of analyses and statistics for the multi-omics data. They recognize that significant amounts of these analyses will need to be redone and some additional analyses were recommended. However, they all agree that these revisions are essential. Please address these concerns as much as possible.

We have preformed a statistical reanalysis of the metabolic and proteomic data sets as per the reviewers requests. In addition we have revamped the ATAC-seq analysis using limma-voom, which is compatible with quantile normalized count data to make significance calls on our ATAC-seq peak data. Lastly, we have employed higher stringency thresholds for many of our other analyses including the differential calls on our bulk RNA-seq and GSEA analysis.

2) Reviewers #2 and #3 pointed out the necessity to provide some biological validations for novel findings from these multi-omics analyses. This will require additional experiments but will greatly improve this manuscript.

We embarked on this project to create a high resolution definition of replicative senescence across multiple modalities wherein the integration of data types would foment novel hypothesis that would not have likely been made looking at any one or two data types (e.g. NNMT→metabolism →chromatin). The effort in collecting, analyzing and integrating all the modalities presented was immense, and we argue that the manuscript stands as a complete piece of work in its own right.

That being said, we wholeheartedly agree that the functional verification of the central regulatory hypothesis would clearly illustrate the utility of this data set.

To that end we have tested the central regulatory hypothesis, that the YAP1/TEAD1 transcription factor activity regulates the EMT/FMT signature in old and senescent WI-38 cells that were isolated in our experiments. Using verteporfin, a molecule that blocks YAP1/TEAD1 activity, we found that inhibition of YAP1/TEAD1 activity robustly down-regulates the published YAP targets in addition to the EMT/FMT signature in older WI-38 cells (Figure 8 and Figure 8—figure supplement 1).

In addition, we have two new analyses based on insightful reviewer questions and observations. First, reviewer number 2 pointed out that the gradual effect of increasing PDL on the transcriptome could arise from earlier onset (prior to sampling) contact inhibition due to the larger size of older cells. To test this hypothesis we scored our single cell data with significant genes that change only during replicative senescence, or only during contact inhibition. We found that only the senescence score increases with PDL suggesting that cell density is not a major confounding factor in our analysis (Figure 2—figure supplement 4).

Second, Reviewer 3 asked if the decrease in accessibility within euchromatin was indicative of a global shut down of transcription. Given that we did observe gene induction and increased accessibility in NADs, we asked if we subdivide the genome into NADs or NAD excluded regions if the chromatin state changes suggest a global shutdown. We found an increase in accessibility with replicative senescence in enhancer and alternate promoter states that reside with NADs (Figure 5—figure supplement 5) This result suggests transcription is not globally inhibited.

In all, these experiments, data, and analyses make the manuscript much stronger and much more compelling in the context of current work at the interface of wound healing, myofibroblast biology, and replicative senescence. We greatly appreciate the reviewers comments and suggestions.

3) All three reviewers have identified the lack of description for certain methods and techniques, mislabeling, typographical errors, etc. Please address all these issues.

We apologize for typos and oversights. We have fleshed out methods and fixed typos.

4) Reviewer #2 pointed out that the proteomics data should be made publicaly available at one of the major proteomics data repositories. Reviewer #1 also pointed out that in-house software codes and scripts for data analyses should also be made publicly available.

The Protemics data has been uploaded to ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD02879.

Additionally all code for figures and data processing can be found here at github repository: https://github.com/dghendrickson/hayflick

Again, we thank the reviewers for their time and thoughts and appreciate their contributions to the revised manuscript.

Reviewer #1 (Recommendations for the authors):

1. There are major concerns about the statistics in the analysis of the dataset.

a. For the differential expression analysis (4.2.3), the authors state: "DEGs were defined as having p-values < 0.001." This is an incorrect method of filtering DEGs as it does not correct for multiple hypothesis testing. This is all the more surprising because DESeq2, used by the authors, natively provides FDR corrected p-values. The data should be re-analyzed after proper multiple hypothesis correction control (native DEseq2 FDR, Bonferroni, or other).

b. Potentially related to the absence of FDR thresholding, the number of DEGs identified-"8969 genes change with increasing PDL … 11652 genes vary with increasing cell density"-seemed, to this reviewer, really high, and likely the result of improper FDR correction. It will be interesting to see how those numbers change with the multiple testing correction.

We have switched to using adjusted p-values generated by DESeq2. We have fixed this oversight in the text as well as methods. The number of significant genes is as follows: 6995 replicative senescence DEGs, 7065 radiation induced senescence DEGs, and 9958 cell density DEGs.

c. Based on the methods, the other analyses should also always incorporate a multiple testing correction, i.e. differential ATACseq peaks (4.4.3 p < 0.001, no FDR correction), metabolomics (4.5 no information on differential analysis) and proteomics (4.6.5 no information on differential analysis).

As per Reviewer 1’s request (see below), we re-analyzed the data with limma-voom and used corrected p-values for differential analysis.

For the metabolites, we have applied a linear model across time for each metabolite in both the hTERT cells and WT WI-38 cells. We now have this analysis described in the methods and include a supplementary data file with the resulting p-values and the FDR adjusted p-values [Storey and Tibshirani, 2003]. We now call out significantly changing metabolites (q-value < 0.05) with asterisks in relevant figures and supplementary figures.

Similarly to the metabolomics, we applied a linear model across time for each protein in both the hTERT cells and WT WI-38 cells. We have now described this analysis in the methods and include a supplementary data file with the resulting p-values and q-values. We also now mark significantly changing metabolites and proteins (adj. p value < 0.01) with asterisks in relevant figures and supplementary figures.

d. When carrying out differential analyses, the authors should always avoid filtering on fold change. They do this for instance in section 2.4 "we focused on compounds exhibiting PDL-dependent changes in abundance > log2 0.5.". Fold-change filtering leads to extremely poor FDR control and should not be used post hoc – the only acceptable method to integrate fold change filtering includes a new statistical framework which takes into account fold change for FDR calculations (see PMC2654802). All results should be presented with this criterion removed or re-calculated with a method that takes fold-change into account.

e. Based on figure legends (Figure 1F, 3A, etc) and methods (4.2.5), FDR correction for multiple testing was also not performed for GSEA results. Please recompute everything including the correction, and only report corrected p-values.

As suggested by the reviewer, the described analysis and significance thresholds have been applied to the metabolomics and proteomics data. The p-values used in the GSEA analysis are Benjamini-Hochberg FDR adjusted pvalues. We have made this clear in text, figure legends and methods. We apologize for the confusion. We thank the reviewer for their attention to detail, and we believe the analyses are now more robust and the conclusions stronger.

f. There are several key issues with the ATAC-seq analysis (4.4.2), which will require a ground-up reanalysis of the data. Quantile normalization of count data and subsequent use of limma for differential accessibility analysis of ATAC-seq (4.4.2, end) is potential problematic. Indeed, limma is not appropriate for count data (it should only be used for continuous data such as microarray) unless the authors mean that they used limma-voom, which is a different algorithm. Please correct this in the revised analyses. In cases where differential depth cannot be corrected for appropriately, the authors should instead use downsampling of mapped reads to match the lowest depth sample. This will remove sequencing depth biases and allow proper count based methods, such as DEseq2, to perform properly.

We have reanalyzed the ATAC-seq data using limma-voom as opposed to limma. In addition, limma-voom, is compatible with quantile normalization as previously described by the voom devloper in Law et al., and others [Shi et al., 2015, Law et al., 2014]. The current analysis using the original limma significance calls on quantile normalized ATAC data have been replaced with the limma-voom significance calls on quantile normalized ATAC-seq data (manuscript Figure 6A and 6B).

2. Additional methodological information is needed for long-term reproducibility of analyses.

a. For reproducibility of code and analyses, all analytical code for this study should be deposited in a repository such as github or made available as a Supplementary file.

b. Please include all version numbers for all used software (e.g. R, Trimmomatic, Salmon, LISA, etc.) packages and R packages (e.g. DESeq2, fgsea, etc.) command parameters where relevant (such as for Trimmomatic, if not available elsewhere). The same should be applied to annotation databases, and if a version number doesn't exist, date of access should be provided.

All code for data processing, figure analysis and plotting, package version history and annotation database version number is now hosted in a github repository:

https://github.com/dghendrickson/hayflick with the appropriate input data files.

Command parameters and annotation database versions and software versions have also been included in methods.

c. All reagents (e.g. DMEM) should include suppliers and catalogue numbers in order to remove any ambiguity as to which product was used.

We have updated the methods section as specified and included a ”Key Reagents” table.

d. Salmon is a pseudoaligning method that uses a transcriptome reference, not a genome. The methods say that the index was derived from hg38, but how and using which information and software is not clear. In general, the authors should make sure that all key analytical steps are explicitly documented in the methods for long term reproducibility of the study.

As mentioned above this information can be located at the Git repository and methods. Specifically we have fleshed out methods sections for our pseudoalignment in manuscript.

3. It is completely unclear which data/plots were generated using the Illumina mRNA kit and which were generated using the Total RNA Kit. Please specify in text and legend.

Reviewer #2 (Recommendations for the authors):

The -omics aspects of the paper are generally well conducted and provide either confirmatory data on pre-existing studies of senescence or add new data (especially on the key FMT pathway and transcription factors involved). Log2FC of +/- 0.5 is somewhat lower than usual cut-offs, so justification is needed where that is selected.

As mentioned in response to reviewer 1, we have revamped the differential testing for metabolomics and proteomics data with higher stringency and FDR adjusted p-values.

No mention is made of quality control steps in RNAseq analysis, but I assume these were conducted?

We have now included a supplementary data file with QC metrics of percent reads aligning to the transcriptome and total number of reads mapped. We included a note in methods that all RNA samples had RIN scores of higher than 7 going into library preparation.

Use of in-house programs for proteomics analysis does mean that others will not be able to conduct the same analyses, though availability of datasets should allow for validation of findings on other data analysis platforms – have these been publicly deposited as links provided are only to RNAseq datasets?

Proteomics methods have been updated to cite the methods used and provide links to code for analysis. Furthermore, data has been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD028799.

From the protocol provided, it is also not clear that the methodology for protein extraction would release chromatin-bound factors – inclusion of a potent nuclease e.g. Benzonase, would ensure that tightly bound proteins are not lost in the spin steps to remove 'cellular debris'.

We understand and agree with the reviewers point. We will take into consideration for future work.

Figure 2E: there appears to be an aberrant pattern for cells at PDL33 that is present in all parts – is 33 an outlier? What is this cluster?

We compiled the genes and ran in GoTerm finder: https://david.ncifcrf.gov/home.jsp and found these enriched terms and keywords (top 10 see below). Having not found a clear cohesive set of enrichments we think this cluster may represent noise.

Phosphoprotein, Acetylation, GO:0005515 protein binding, Nucleus, GO:0005737 cytoplasm, GO:0005829 cytosol, Isopeptide bond, GO:0005634 nucleus, hsa04071:Sphingolipid signaling pathway,Cytoplasm

The cell biology aspects are much less well described making it impossible for others to reproduce the work – in particular, cell seeding density is critical but not even mentioned for longitudinal culture, and Figure 1D images suggest significant cell crowding. Along the senescence trajectory, primary fibroblasts grown in monolayers increase in surface area (and to some extent in volume) and populations at later PDLs undergo contact inhibition at lower cell densities than those at earlier PDLs so it may be the case that some of the 'gradual' changes (especially those associated with cell cycle changes) reflect simply quiescence form contact inhibition rather than steps to senescence. Methodology to count cells and calculate PDL is needed as different labs may conduct this differently. Daily microscopy inspection of cultures to assess appropriate time points for subculturing is essential – simply choosing to split at 4 or 7 days according to PDL does not account for the morphological changes that occur across the replicative lifecourse.

We appreciate the reviewer’s attention to detail. We have added more detailed descriptions to our methods section for cell culture including the equation we used to calculate PDL. We have also added a supplementary figure with representative cell density images at time of harvest. We have also added a supplementary methods document outlining the extensive testing we carried out prior to initiation of the time course and the reasoning behind the sampling windows chosen. In our preliminary experiments testing the effect of cell density on gene expression we found that cells reach contact inhibition quicker at younger PDLs as evidenced by inhibition of cell cycle genes and up-regulation of the hypoxia response (Figure 1—figure supplement 5).

The reviewer raises an insightful point with respect to increased cell size that we had not ruled out with respect to the gradual increase in the senescent gene expression program. Although great care was taken to control cell density this is entirely plausible. Accordingly, we have added the following text, analysis and updated supplementary figure. Below is new in manuscript:

“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 providing additional evidence for a gradual increase in the senescent gene expression program that is independent from changes in cell density.”

The following text is the additional supplementary methods.

“0.1 Sampling proliferating WI-38 cells at different cell densities confounds assessment of cell cycle state and senescence

A limited pilot (performed prior to RS time course) experiment was conducted to determine optimal sampling conditions. Using gene expression as a readout for cellular state, we performed RNA-seq on PDL 20, PDL 29, and PDL 39 WI-38 cells. We plated an equal number of cells in triplicate for the three PDLs and cells were collected for RNA extraction after 4 days of growth. We looked for any significant changes in transcript levels between the three samples.”

We observed two notable trends. First, we found stronger expression of genes associated with cell cycle and mitosis in the older PDL 39 cells compared to PDL 20 and PDL 29 cells. This was unexpected based on previous observations of slower growth and less active mitosis in cells closer to senescence.(Hayflick 1965). Second, the cells at PDL 20 and PDL 29 expressed genes associated with hypoxic stress (Figure 1—figure supplement 5A).

These data, in tandem with the observation that older cells grow at a slower rate, suggest that the highly proliferative younger cells (PDL 20 and 29) grow to higher cell density faster and therefore encounter earlier in real time, the various normal mitosis-inhibiting stresses, e.g. hypoxia, contact inhibition, and starvation that trigger cell cycle exit (Dayan et al., 2009).

To understand the effects of growth rates and cell density on gene expression, we performed a cell density study (CD) wherein young PDL cells were sampled at increasing levels of cell density to measure the effects of cell density and proliferation on gene expression (Figure 1A). We observed that based on the seeding density we employed (to minimize the stress of sparse plating) the sampling window we had originally targeted (4 days after seeding) was located in a volatile region with respect to proliferation rate and hypoxic stress using Cyclin B1 (CCNB1) and VegfA (VEGFA) expression as proxies for each respectively (Figure 1—figure supplement 5B).

Accordingly, we modified the sampling window by choosing a time interval (2.5-3 days) prior to the onset of VEGF-A induction and CCNB1 repression (Figure 1—figure supplement 5B). For older, slower growing cells in the plateauing region of the replicative senescence process (PDL > 40) we shifted the sampling window to 3.5-4 days based on visual inspection and cell counts. Ultimately in the actual extended time course, we tracked cell density with a phase contrast EVOS microscope to ensure cellular confluence was similar across all time points during sampling (Figure 2—figure supplement 5).

It would therefore be helpful to determine to what degree these findings are cell type specific and which changes represent general cell senescence, especially as multiple RNAseq datasets exist on lung and skin fibroblast senescence with various mode of senescence induction.

As mentioned, we do compare our data set to a model compiled by Hernandez-Segura 2017 over multiple skin and lung fibroblast cell lines approaching replicative senescence and find high concordance (manuscript Figure 1D). Although we agree with the reviewer that it would be informative to expand this to a comprehensive analysis of all types of senescence induction, it is outside the scope of the manuscript.

The Methods mentions 'deep senescence' but it is not clear in the main text what datasets are derived from this. Images of cells at different points along the pathway to senescence would be reassuring, and orthogonal validation of senescence markers e.g. p16, p21 by qRT-PCR is needed. It would be extremely helpful if the RNA seq and protein changes reported were cross validated by conventional techniques.

We thank the reviewer for these suggestions. The ’deep senescence’ time course consists of a second time course and three samples at 45, 52, 53 PDLs that were used for bulk RNA-seq. We have added a sample manifest that clearly shows what assays were performed on which samples and time courses (Figure 1—figure supplement 1). We added the second deep senescence time course after De Cecco et al., reported transcriptomic differences between early and deep senescence to ensure we were not missing a crucial aspect of replicative senescence [De Cecco et al., 2019]. The PDL 45 cells served as a bridge sample for batch correction so that we could accurately compare the deep senescence time points to first time course. We have added clarifying language and updated manuscript Figure 1B.

Additionally, as mentioned we now have images of cells at harvest for original time course (Figure 2—figure supplement 5).

For orthogonal validation, we would point out that each different modality of ”omics” experiments validates the other orthogonally. For example, many of the genes that we highlight as induced at the transcriptomic level (e.g. NNMT or TGFB2) are also induced in our proteomics data, our single cell RNA-seq data, and have significant change in their promoter regions in our ATAC-seq data. The strength of this work in part rests with this intrinsic cross-validation across modalities. In addition, all of our experiments were performed with replicates and at adequate depth for NGS experiments.

In addition, a plethora of previous work has demonstrated that qPCR and RNA-seq data are highly correlated [Wu et al., 2014, Shi and He, 2014, Asmann et al., 2009, Griffith et al., 2010]. We argue that the orthogonal validation across multiple modalities is substantial validation and sufficient for the conclusions drawn. In addition we argue superior to qPCR in terms of satisfying orthogonal validation given that RNA-seq and qPCR are both based on reverse transcription and cDNA amplification.

We have collated the RNA-seq and proteomics data for p21, p16, and β-galactodsidase along side NNMT and TGFB2 to demonstrate the overlap between technologies (Figure 1—figure supplement 2). For p16 and p21 we have also included Western blots to confirm as the induction in the proteomics data exhibits high variance.

Although studies are conducted on the influence of cell density on -omics changes over a short time course of 10 days, the effect of cell cycle exit in senescence would be better accounted for by comparing with genuinely quiescent cell populations e.g. serum starvation of the hTERT cells eg Figure 3A

Similar to comparing replicative senescence to other senescence inductions, we agree with the reviewer that a full workup of comparing senescence to many types of cell cycle inhibition would further our understanding about the regulatory logic of cell cycle exit and what may or may not be specific to replicative senescence. However, our aim in using the confluence condition was to account for cell cycle inhibition driven by cell density as we reasoned it was the most likely source of variability in our experimental set up. A comprehensive comparison to multiple types of cell cycle exit, although an important data set, is outside the scope of this work.

Time points used for comparisons are poorly labelled – are hTERT samples taken at daily intervals or at the same time points in long term culture as the RS samples? The text states that the aim is to control for the effects of long-term cell culture, but this is not apparent from later figures. Figure 1D implies that cell population at the same PDL are used for RS and hTERT groups (were cells transduced at PDL20 to allow analysis at PDL25 allowing for drug selection steps?) yet later on time points are simply given numbers e.g. 1,7 in Figure 2B and 2,4,5,6,7 in Figure 3A. It is not even stated if these are days, weeks or sample numbers.

We apologize for the confusing description. We have added clarifying text to emphasize that the hTERT timepoints are matched to WT PDL samples. We have also added a visual dotted line connection in manuscript Figure 1B to highlight that the samples are paired. The transduction to generate the WI-38 hTERT cells happened far in advance of the time course. We have updated the methods to reflect this.

A comparison between early PDL WT WI38 with hTERT WI38 should be conducted (the authors have the relevant data sets) to ensure that the hTERT expression simply leads to cell immortality without impacting on gene expression patterns – this is critical since telomerase has roles other than at the telomere, including in gene expression control, RNA splicing and mitochondrial metabolism. hTERT transduction is very poorly described in the methods section so that others will not be able to reproduce the experiments: "appropriate target plasmid and packaging constructs" are completely uninformative, as is 'selected. with a selection drug". Specific details are needed.

The selection process in generating the hTERT line, has made meaningful comparison to address the specific effect of hTERT over-expression impossible. To test that, we would need an inducible hTERT system to follow cells after hTERT induction. In this study the hTERT cells were only a control for the long term culture effects and not to test the effects of hTERT expression.

We agree with the reviewer that this would be a worthwhile comparison if we had the appropriate experimental set up. We have clarified our methods concerning the generation of the WI-38 hTERT cells to make this clear. The transduction to make the hTERTs occurred far in advance of the time course.

Experiments on radiation damage to induce senescence assay relatively early time points – cell cycle exit is an early event in the DNA damage response but not indicative of senescence per se. True damage-induced senescence onset should be monitored at later time points post-damage (eg from 14 days) to avoid conflating the acute DDR with senescence-specific gene expression patterns. This could account for the large discrepancy between RS and RIS datasets here.

We agree with the reviewer. We have changed the language surrounding comparisons to the radiation induced senescence condition and have added clarifying context where relevant and in the discussion.

The oxidative phosphorylation data are interesting and consistent with known changes on senescence (though the effect on ROS generation should be considered more). However, the citation in this context of papers that refer to lipid peroxidation – a state of lipid damage – rather than physiological energy generation through β oxidation of fats. Please take care not to conflate energy generating β oxidation and damaging oxidation in the citations.

We apologize we are not sure which reference the reviewer is referring to. The only paper with ”lipid peroxidation” in the title we cited was Flor et al., 2017, A signature of enhanced lipid metabolism, lipid peroxidation and aldehyde stress in therapy-induced senescence [Flor et al., 2017]. We cited this paper in the following section summing up our metabolic observations.

Indeed, metabolomic data collected from various types of senescence models has shown increased lipid oxidation, lipid accumulation, TCA up-regulation, and glycolytic alterations [Flor et al., 2017].

We cited Flor et al., in reference to increased lipid oxidation as pertaining to fatty acid β-oxidation as fatty acids are common components of complex lipids. Flor et al., observed increased β-oxidation in Figures 1A,B, C and Figure 4 with therapy induced senescence. To remove any confusion we have changed the text to read ”shown increased fatty acid oxidation, lipid accumulation, TCA up-regulation, and glycolytic alterations”

The NNMT data are fascinating and really provocative – did the authors consider including bisulfite sequencing to assess DNA methylation state in conjunction with the ATAC-seq data? The finding of CTCF upregulation is also interesting in the context of the LAD/NAD data – would 3C/Hi-C give clearer data on overall chromatin configuration? (It would be intersting to see if TADs change significantly, which would be predicted from the data here).

We agree with the reviewer and find this result very exciting. We are in fact working on directly measuring DNA methylation chromatin conformation during replicative senescence but as part of a future publication related to the findings here.

The manuscript has not been fully proof-read which is somewhat frustrating for the reader e.g. the first few references don't even have volume or page numbers, there is incorrect use of upper/lower case in some places (e.g. amp instead of AMP) and there are various typos throughout. Methods section has gaps/omissions e.g. 4.4.3 refers to Figure ??B, 4.5.1 refers to "(vendor?)" and 4.7.2 refers to "(Cite CIS-BP)". Tenses should be consistent.

We apologize and have gone through the manuscript to fix typos.

Most worryingly, very little effort has been made to ensure the Supplementary figures are informative – figure legends are lacking, labels of figure parts are missing e.g. time points/PDLs on Figure S2, Figure S3; the interpretation that EMT (or rather FMT) is a feature of sesncence depends to some extent on the data in Figure S2 but these appear to show that there is no change of the progression from early proliferation to senescence (though hard to tell as PDLs not given)?

We have remade many of the supplementary figures and main figures with more labels. The GSEA enrichment can be hard to interpret with respect to magnitude of change, i.e. even in the face of monotonic increase, the p-value of GSEA enrichment would not necessarily go up as the analysis is based on rank rather than magnitude. We thank the reviewer for calling this out as it is central to the EMT/FMT conclusions drawn. We have updated the figure to include the actual log2 fold changes of the genes driving the enrichment at each PDL to emphasize the progression from early proliferation to senescence (Figure 1—figure supplement 4).

statements such as "S phase and G2M cells were isolated" – this is wholly misleading as the cells were not isolated physically according to cell cycle stage (eg FACS), but the data were processed post-hoc to determine which cell cycle stage they most closely fit with according to gene expression patterns.

We agree with the reviewer and have fixed the language.

Figure S8 legend states cells at PDL20, figure labels cells at PDL25; time points/cell ages are missing for Figure S9- and cannot state mitochondrial functions are already simply by looking at protein levels – safer to say that changes are consistent with altered mitochondrial metabolism. Again, PDL20 is referred to yet throughout the paper, the earliest PDL analysis is at 25.

We have added a label key for the supplementary figure. Reference time points for log2 fold change heatmaps are not shown as values are at or near (when mean subtracted) zero. We have emphasized what is being compared (e.g. vs. PDL20 rather than first timepoint) and what is not being shown in figure legends.

Figure S10 is uninterpretable as no X axis is given, and time points/PDLs are also missing in Figure S11B. Figure S12D – what do the boxes represent? The lack of labels and suitable legend make it impossible to interpret. Why analyse PDL 46 here, as cells are not yet senescent? Is the whole of chromosome 22 shown in part E? Please provide Mb scale – and state what platform was used to view (e.g. IGV?).Figure S14 – GO terms would be a better way of showing the processes than the word clouds used. Figure S15 lacks a legend – myofibroblast markers (from lit) – I presume from literature? Which are the fibrillar collagens (the collagen number matters).

We have added PDL and timepoint labels to all figures. The long arm of Chr22 is shown and that is reflected in figure now.

PDL 46 was used as the changes in epigenome we report begin prior to PDL 46. As the observed heterochromatic opening begins prior to senescence, there was not a requirement to use senescent cells. If these changes were driven by dead cells, then treatment of cells with PMA prior to ATAC-seq would have altered the results in PDL 46.

The ATAC-seq supplementary figure was rendered in R. We have the relevant packages listed and code is now accessible. The heatmap from the pseudotime clustering was made using 800 GO terms that cannot fit into figure. We have now included the data underlying the heatmap for readers who would like to see the GO terms that compose the word clouds.

We have removed Figure S15 in the revised manuscript- However Col5A2 was labeled in the figure.

When describing what samples the data shown are compared against, the terminology "the first sample" needs more information – is this PDL25 for RS? But this 'first sample' appears later in metabolomic analysis eg Figure 3E (hence fold changes will not be comparable with proteomics).

We have added more explicit labels and included a sample manifest (Figure 1—figure supplement 1) to make it easier for readers to see what the time points are available for each modality. It is true that we do not have the exact same time course samples for all modalities e.g. metabolomics and proteomics.

However, we did not, and do not want readers to directly compare fold changes. Even if both metabolomics and proteomics had the exact same reference point, there is not an expectation that any of the modalities should scale similarly in magnitude or time. The aim of the figures is to show change over time from proliferating state to replicative senescence. It was extraordinarily hard to seed and sample the cells required for these experiments. When we lost the PDL20 sample for metabolomics and single cell RNA-seq, we opted not to throw out PDL20 for all other modalities as the change over time with respect to proliferation state was the key parameter and PDL25 and PDL20 are both very much in the proliferative state (manuscript Figure 1B).

Overall this could be a really great paper, bringing together a range of potent -omics techniques to study replicative senescence at high resolution cross time. Most of the paper is well written and provides a clear and sensible discussion of the work carried out, the data and implications of the results. It is somewhat dense (which may be inevitable with multi-omics studies) and could benefit from more accessible explanation of some of the analyses performed (eg readers will be more familiar with tSNE than UMAP – a brief explanation would help here).

We thank the reviewer for praise and appreciate the criticisms which helped improve the manuscript overall. We have added a quick UMAP explanation in the manuscript.

However, the paper is let down in parts by sloppy presentation, lack of critical experimental detail (meaning the work cannot be reproduced by others), poor design (or one hopes simply poor description) of comparators, an apparent lack of understanding of the nuances of senescent cell behaviour in culture, and by lack of statical validation of the simpler experiments (no n values noted in some cases – e.g. Figure 1D, Figure 5A, C, Figure S3 etc etc).

Hopefully, the changes we have made will strengthen the paper. Specifically, we hope that the additional methods and supplementary text concerning cell culture and hTERT cell line generation convince the reviewer that we took great care in our design and cell culture. As mentioned, we specifically chose cell density as a control as cell density is most likely to be altered between sampling periods especially in cells changing in size and growth rate. For example, the cell density experiment combined with the single cell data bolsters our conclusion of ’gradual’ onset senescence and provides evidence beyond that our careful work with respect to sampling density was successful (Figure 2—figure supplement 4).

We have added statistical validation of simpler experiments in manuscript Figure 1. manuscript Figure 5A was an exploratory analysis that was simply a starting point for manuscript Figure 5B, i.e. we would have proceeded with higher resolution analysis even if global percentage changes were not significant. The significance testing for the analysis in manuscript Figure 5C, originally just in text, has been added to the figure.

Confirmatory experiments to cross validate top hits from the -omics studies are also required. Mining the existing proteomics data for PTMs should also be possible and may add some value, though not essential here.

As mentioned, we argue that the different modalities serve as sufficient cross validation and that experiments to confirm changes that are robust across multiple disparate data types would not add to conclusions or invalidate them. For example, an inconsistency from a western blot of NNMT could signify that the purchased antibody has poor binding or specificity since we have available data to validate from the gene level (RNA-seq), protein level (proteomics), and readout of cellular metabolic state (metabolomics). A simple western blot or qPCR would not strengthen the conclusions further. However as mentioned above, we have included p21 and p16 due to variance in proteomics data.

Reviewer #3 (Recommendations for the authors):

1. There are some technical points about ATAC-seq that need to be clarified. (a) global chromatin accessibility may be higher in senescent cells. Could the authors share the bioanalyzer profiles of the ATAC-seq libraries from the different timepoints? Are they similar? (b) senescent cells typically produce far more mitochondrial reads. Could the authors confirm that sufficient depth of sequencing was achieved in their ATAC-seq runs? These data including sequencing statistics, peak coordinates, fragment distribution graphs etc. should be included in the Supplementary Section.

A) We have added a supplementary figure with representative bionalyzer traces (Figure 5—figure supplement 1 and 2). ATACseq libraries from all PDLs sampled share the same fragment length distribution before and after sequencing. We employed a size selection cut-off to enrich for nucleosome free regions.

B) We have added a supplementary figure with mitochondrial percentages and transcription start site enrichment (Figure 9). We do not observe an increase in mitochondrial percentage (less than 2% across all samples) with increasing PDL. We attribute the low mitochondrial percentages to the use of the Corces et al., ”Omni-ATAC” protocol, which claims a 13-fold decrease in reads that map to mtDNA compared to the standard ATAC-seq protocol [Corces et al., 2017].

In addition, we have included a supplementary file with read depth, duplicate read percentages, mitochondrial read percentage, alignment statistics, and TSS score for each sample–all of which pass the ENCODE pipeline standards: https://www.encodeproject.org/atac-seq/.

Significant peak coordinates were included in original manuscript and in revised manuscript. They can also be found in GEO repository along with bedgraphs of alignment and coverage. We included representative fragment length distribution plots which mirror the size selected ATAC-libraries sequenced (Figure 5—figure supplement 1, final pool bioA trace).

2. The ATAC-seq data might benefit from a RepEnrich analysis to identify peaks in repeat elements that are known to be heterochromatinized in proliferating cells and derepressed in senescence. This, in my opinion, is a better way to look at constitutive heterochromatin desilencing compared to NAD/LAD overlap and is also independent of NAD/LAD calls in a different cell line (in this case IMR-90).

We appreciate the reviewer’s suggestion. However, for our analysis we specifically were interested in experimentally validated NAD and LAD domains which contain both gene bodies and heterochromatin as opposed to solely constitutive heterochromatin (see point 3 below). We argue that use of the IMR-90 labels is justified based on the observed decrease in transcription start site coverage in Figure 5—figure supplement 3. These analyses was carried out solely based on genomic gene annotation independent of cell type. That we observe the similar results when using chromatin state labels derived from IMR-90 cells (also fetal lung fibroblasts) independent of gene annotation (manuscript Figure 5A and 5B) suggests that the IMR90 derived labels are accurately reflecting WI-38 chromatin state.

We think that digging deeper into constitutive heterochromatin especially in light of our NNMT findings is an exciting avenue for further research and is part of ongoing follow up work, and we will look into RepEnrich analysis.

3. What are the authors thoughts about the decrease in chromatin accessibility at promoters, enhancers, and gene bodies in their ATAC-seq data (Figure 5B)? Does this imply a global shut down of transcription? Is it reflected in the transcriptomic analyses?

This is an important question. We interpret the results as a winnowing down or specialization of function in senescent cells (e.g. permanent cell cycle exit, etc. ). Although the aggregate pattern is down, we do observe promoters that exhibit increased accessibility resulting in higher transcription of neighboring genes and increased protein levels e.g. NNMT and TGFB2. In fact, in response to this question, we carried out the same analysis used in the manuscript Figure 5B (changing peak accessibility by chromatin state label) and instead divided into peaks overlapping or excluded from NADs (Figure 5—figure supplement 5). We found that ATAC-seq peaks within NADs exhibit an altered chromatin state profile. For NAD exluded peaks, the results largely mirror the results for all peaks (increased accessibility in closed regions and vice versa). Conversely, multiple chromatin states in the enhancer and promoter categories, which decrease in accessibility with PDL in NADs excluded peaks, display increased accessibility with increasing PDL only in NAD overlapping peaks.

Given that we also observe overlap with replicative senescent gene expression from bulk RNA-seq data and NADs (manuscript Figure 5C), we interpret this result as consistent with the hypothesis that transcription is not entirely shutting down rather so much as changing where it is happening. Basically, although widespread, the reduced accessibility in euchromatin observed in aggregate, is not reflective of a pervasive shutdown of transcription across the entire genome.

We appreciate the reviewer’s question. We have added this result as a supplement to the main text.

4. While the authors perform an in-depth analyses of transcription factors most likely to drive RS, they haven't really shown the direct binding of these factors to chromatin in RS (and not hTERT, RIS or CD cells) and/or the change in their expression in senescence. Additionally, no functional experiments are performed by knocking down or overexpressing TEAD1 and investigating its effect on senescence. If possible, it would be nice to include this data (i.e., a western blot, genomic binding of TEAD1 and some functional experiments).

We embarked on this project to create a high resolution definition of replicative senescence across multiple modalities wherein the integration of data types would incite novel hypothesis that would not have likely been made looking at any one or two data types (e.g. NNMT →metabolism →chromatin). The effort in collecting, analyzing and integrating all the modalities presented was extensive, and we argue that the manuscript stands as a complete piece of work in its own right. However, we do agree that the functional verification of the central regulatory hypothesis would clearly illustrate the utility of this data set.

To that end, we have taken late passage PDL 40 WI-38 cells wherein increase in expression of the EMT signature was observed, (Figure 5B) and treated with verteporfin, a molecule that blocks YAP1/TEAD1 interaction as demonstrated by [Mascharak et al., 2021]. We preformed RNA-seq on these treated cells and found robust inhibition of known YAP targets, TGFB2, and the EMT signature that started this line of inquiry (Figure 8—figure supplement 1). These data are consistent with a regulatory role for YAP1/TEAD1 in establishing the EMT/FMT gene expression program in WI-38 cells approaching replicative senescence.

These data make the paper stronger and highlight the utility and relevance of the generated data sets. We thank the reviewer for this suggestion.

References

Asmann et al., 2009 Asmann, Y. W., Klee, E. W., Thompson, E. A., Perez, E. A., Middha, S., Oberg, A. L., Therneau, T. M., Smith, D. I., Poland, G. A., Wieben, E. D., et al. (2009). 3’tag digital gene expression profiling of human brain and universal reference rna using illumina genome analyzer. BMC genomics, 10(1):1–11.

Corces et al., 2017 Corces, M. R., Trevino, A. E., Hamilton, E. G., Greenside, P. G., Sinnott-Armstrong, N. A., Vesuna, S., Satpathy, A. T., Rubin, A. J., Montine, K. S., Wu, B., et al. (2017). An improved atac-seq protocol reduces background and enables interrogation of frozen tissues. Nature methods, 14(10):959–962.

De Cecco et al., 2019 De Cecco, M., Ito, T., Petrashen, A. P., Elias, A. E., Skvir, N. J., Criscione, S. W., Caligiana, A., Brocculi, G., Adney, E. M., Boeke, J. D., et al. (2019). L1 drives ifn in senescent cells and promotes age-associated inflammation. Nature, 566(7742):73–78.

Flor et al., 2017 Flor, A. C., Wolfgeher, D., Wu, D., and Kron, S. J. (2017). A signature of enhanced lipid metabolism, lipid peroxidation and aldehyde stress in therapy-induced senescence. Cell death discovery, 3(1):1–12.

Griffith et al., 2010 Griffith, M., Griffith, O. L., Mwenifumbo, J., Goya, R., Morrissy, A. S., Morin, R. D., Corbett, R., Tang, M. J., Hou, Y.-C., Pugh, T. J., et al. (2010). Alternative expression analysis by rna sequencing. Nature methods, 7(10):843–847.

Hernandez-Segura et al., 2017 Hernandez-Segura, A., de Jong, T. V., Melov, S., Guryev, V., Campisi, J., and Demaria, M. (2017). Unmasking transcriptional heterogeneity in senescent cells. Curr. Biol., 27(17):2652–2660.e4.

Law et al., 2014 Law, C. W., Chen, Y., Shi, W., and Smyth, G. K. (2014). voom: Precision weights unlock linear model analysis tools for rna-seq read counts. Genome biology, 15(2):1–17.

Mascharak et al., 2021 Mascharak, S., Davitt, M. F., Griffin, M., Borrelli, M. R., Moore, A. L., Chen, K., Duoto, B., Chinta, M., Foster, D. S., Shen, A. H., et al. (2021). Preventing engrailed-1 activation in fibroblasts yields wound regeneration without scarring. Science, 372(6540).

Ogrodnik et al., 2019 Ogrodnik, M., Salmonowicz, H., Jurk, D., and Passos, J. F. (2019). Expansion and cell-cycle arrest: common denominators of cellular senescence. Trends in biochemical sciences, 44(12):996–1008.

Shi et al., 2015 Shi, W., Liao, Y., Willis, S. N., Taubenheim, N., Inouye, M., Tarlinton, D. M., Smyth, G. K., Hodgkin, P. D., Nutt, S. L., and Corcoran, L. M. (2015). Transcriptional profiling of mouse b cell terminal differentiation defines a signature for antibody-secreting plasma cells. Nature immunology, 16(6):663–673.

Shi and He, 2014 Shi, Y. and He, M. (2014). Differential gene expression identified by rna-seq and qpcr in two sizes of pearl oyster (pinctada fucata). Gene, 538(2):313–322.

Storey and Tibshirani, 2003 Storey, J. D. and Tibshirani, R. (2003). Statistical significance for genomewide studies. Proceedings of the National Academy of Sciences, 100(16):9440–9445.

Wu et al., 2014 Wu, A. R., Neff, N. F., Kalisky, T., Dalerba, P., Treutlein, B., Rothenberg, M. E., Mburu, F. M., Mantalas, G. L., Sim, S., Clarke, M. F., et al. (2014). Quantitative assessment of single-cell rna-sequencing methods. Nature methods, 11(1):41–46.

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

Article and author information

Author details

  1. Michelle Chan

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Visualization, Writing - original draft, Writing - review and editing
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  2. Han Yuan

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Data curation, Formal analysis, Methodology, Visualization
    Contributed equally with
    Ilya Soifer
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  3. Ilya Soifer

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Software, Visualization
    Contributed equally with
    Han Yuan
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  4. Tobias M Maile

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Formal analysis, Investigation, Methodology
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  5. Rebecca Y Wang

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Investigation, Methodology
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  6. Andrea Ireland

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Investigation, Methodology
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  7. Jonathon J O'Brien

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Formal analysis, Investigation, Methodology
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9660-4797
  8. Jérôme Goudeau

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Investigation, Methodology
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2483-1955
  9. Leanne JG Chan

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Investigation, Methodology
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  10. Twaritha Vijay

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Investigation, Methodology
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  11. Adam Freund

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Project administration, Supervision
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7956-5332
  12. Cynthia Kenyon

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Project administration, Supervision
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3446-2636
  13. Bryson D Bennett

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Project administration, Supervision
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  14. Fiona E McAllister

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Project administration, Supervision
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  15. David R Kelley

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Project administration, Supervision
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  16. Margaret Roy

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Conceptualization, Investigation, Methodology, Project administration, Supervision
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  17. Robert L Cohen

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Conceptualization, Project administration, Supervision
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  18. Arthur D Levinson

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Conceptualization, Project administration, Resources, Software
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
  19. David Botstein

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Conceptualization, Project administration, Resources, Supervision, Writing - review and editing
    For correspondence
    botstein@calicolabs.com
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9499-4883
  20. David G Hendrickson

    Calico Life Sciences, LLC, South San Francisco, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Supervision, Visualization, Writing - original draft, Writing - review and editing
    For correspondence
    dgh@calicolabs.com
    Competing interests
    is affiliated with Calico Life Sciences, LLC. The author has no other competing interests to declare
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1884-5234

Funding

No external funding was received for this work.

Acknowledgements

The authors thank Ari Firestone, Jason Rogers, Nick Bernstein, Jeff Settleman, Rochelle Buffenstein, Emily Stoops, Magdalena Preciado Lopez, and Antoine Roux for their feedback on the experiments, analysis and the manuscript. In addition the authors also thank Leonard Hayflick for guidance and discussion.

Senior Editor

  1. Kevin Struhl, Harvard Medical School, United States

Reviewing Editor

  1. Weiwei Dang, Baylor College of Medicine, United States

Reviewers

  1. Lynne Cox
  2. Payel Sen, National Institute on Aging, United States

Publication history

  1. Preprint posted: May 4, 2021 (view preprint)
  2. Received: May 12, 2021
  3. Accepted: January 31, 2022
  4. Accepted Manuscript published: February 4, 2022 (version 1)
  5. Accepted Manuscript updated: February 8, 2022 (version 2)
  6. Version of Record published: March 18, 2022 (version 3)

Copyright

© 2022, Chan et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 3,189
    Page views
  • 638
    Downloads
  • 2
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Michelle Chan
  2. Han Yuan
  3. Ilya Soifer
  4. Tobias M Maile
  5. Rebecca Y Wang
  6. Andrea Ireland
  7. Jonathon J O'Brien
  8. Jérôme Goudeau
  9. Leanne JG Chan
  10. Twaritha Vijay
  11. Adam Freund
  12. Cynthia Kenyon
  13. Bryson D Bennett
  14. Fiona E McAllister
  15. David R Kelley
  16. Margaret Roy
  17. Robert L Cohen
  18. Arthur D Levinson
  19. David Botstein
  20. David G Hendrickson
(2022)
Novel insights from a multiomics dissection of the Hayflick limit
eLife 11:e70283.
https://doi.org/10.7554/eLife.70283
  1. Further reading

Further reading

    1. Cancer Biology
    2. Chromosomes and Gene Expression
    Arnaud Carrier, Cécile Desjobert ... Paola B Arimondo
    Research Article Updated

    Aberrant DNA methylation is a well-known feature of tumours and has been associated with metastatic melanoma. However, since melanoma cells are highly heterogeneous, it has been challenging to use affected genes to predict tumour aggressiveness, metastatic evolution, and patients’ outcomes. We hypothesized that common aggressive hypermethylation signatures should emerge early in tumorigenesis and should be shared in aggressive cells, independent of the physiological context under which this trait arises. We compared paired melanoma cell lines with the following properties: (i) each pair comprises one aggressive counterpart and its parental cell line and (ii) the aggressive cell lines were each obtained from different host and their environment (human, rat, and mouse), though starting from the same parent cell line. Next, we developed a multi-step genomic pipeline that combines the DNA methylome profile with a chromosome cluster-oriented analysis. A total of 229 differentially hypermethylated genes was commonly found in the aggressive cell lines. Genome localization analysis revealed hypermethylation peaks and clusters, identifying eight hypermethylated gene promoters for validation in tissues from melanoma patients. Five Cytosine-phosphate-Guanine (CpGs) identified in primary melanoma tissues were transformed into a DNA methylation score that can predict survival (log-rank test, p=0.0008). This strategy is potentially universally applicable to other diseases involving DNA methylation alterations.

    1. Chromosomes and Gene Expression
    Desislava P Staneva, Stefan Bresson ... Robin C Allshire
    Research Article

    Kinetoplastids are a highly divergent lineage of eukaryotes with unusual mechanisms for regulating gene expression. We previously surveyed 65 putative chromatin factors in the kinetoplastid Trypanosoma brucei. Our analyses revealed that the predicted histone methyltransferase SET27 and the Chromodomain protein CRD1 are tightly concentrated at RNAPII transcription start regions (TSRs). Here we report that SET27 and CRD1, together with four previously uncharacterized constituents, form the SET27 promoter-associated regulatory complex (SPARC), which is specifically enriched at TSRs. SET27 loss leads to aberrant RNAPII recruitment to promoter sites, accumulation of polyadenylated transcripts upstream of normal transcription start sites, and conversion of some normally unidirectional promoters to bidirectional promoters. Transcriptome analysis in the absence of SET27 revealed upregulated mRNA expression in the vicinity of SPARC peaks within the main body of chromosomes in addition to derepression of genes encoding variant surface glycoproteins (VSGs) located in subtelomeric regions. These analyses uncover a novel chromatin-associated complex required to establish accurate promoter position and directionality.