University of Southern California Dissertations and Theses

A single cell time course of senescence uncovers discrete cell trajectories and transcriptional heterogeneity

(USC Thesis Other)

A single cell time course of senescence uncovers discrete cell trajectories and transcriptional heterogeneity

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content
A SINGLE CELL TIME COURSE OF SENESCENCE UNCOVERS DISCRETE CELL TRAJECTORIES AND
TRANSCRIPTIONAL HETEROGENEITY

by

Serban Ciotlos

A Dissertation Presented To
THE FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(BIOLOGY OF AGING)

May 2023

Copyright 2023 Serban Ciotlos-Ellis
ii

ACKNOWLEDGEMENTS
I am deeply grateful for the countless individuals who have supported me throughout my
graduate school journey. From mentors and instructors at all levels, their guidance and encouragement
have been instrumental in helping me reach where I am today. In particular, I would like to express my
heartfelt appreciation to my Dissertation Committee members, Drs. Simon Melov, Judy Campisi, Lisa
Ellerby, Bérénice Benayoun, and John Tower. Their invaluable feedback and insights have strengthened
my research and propelled me forward in my academic pursuits. I am truly grateful for their useful
advice and suggestions.
In particular, I am deeply grateful to my PhD advisors, Drs. Simon Melov and Judy Campisi, for
their mentorship, scientific expertise, and unwavering support throughout my doctoral studies. Their
guidance has been invaluable in shaping my research and personal development, and I feel fortunate to
have had the opportunity to work in their labs. Simon, thank you for your continuous guidance and
discussions, openness to novel ideas, and push for integrating new and revolutionary techniques in our
lab and into my studies. Your approach and philosophy have allowed to me gain extremely valuable
skills for the future of biology aging. Judy, our discussions about the complex and tricky nature of
cellular senescence were instrumental in pushing the boundaries of my thinking, and feedback and
support on my experiments and studies were critical to developing my studies.
I would like to express my deepest gratitude to my parents for their unwavering support,
patience, listening, and faith in me throughout my academic journey. Your love and encouragement
have sustained me through both the highs and the lows, and I am forever grateful for their presence in
my life and encouragement of my academic and professional pursuits. Mom, Dad, Grandma, Cassie –I
wouldn’t be where I am today without your words, at times gentle and at times bringing reality back
iii

into focus. It was a challenging journey, and I wouldn’t have started or finished it without you, ube.
Nancy, Steve, and Mark – over the last 5 years you have become part of my family, and likewise you
have welcomed me into your families. You’ve listened to my rants and updates, asked great questions,
and your love and support over the years has been instrumental to my sanity and graduation. I’m
excited to continue making memories with all of you.
To my wife Cassie, you’ve basically taken this journey with me since the beginning. Your
presence and love has been a true blessing, and I wouldn’t have made it here without you. I feel like
you’ve earned a PhD-by-osmosis. Our mutual love, respect, and support has been the backbone of my
life. Thank you for putting up with me and the demands and chaotic nature of my studies. I treasure
every day we’ve spent together, and I’m excited for our future together.
I would like to thank all members of the Melov and Campisi labs, past and present, for their
advice and discussions, which have enriched my research and curiosity about aging biology. Your diverse
areas of expertise and perspectives have broadened my scientific horizons and helped me tackle
complex research questions. To my dear friends and colleagues Andrew Cruz, Angelina Holcom, Asia
Davis, Brendan Hughes, Carlos Galicia, Francesco Neri, Jacob Rose, Lauren Wimer, Nicholas Martin,
Ryosuke Doi, Sofiya Galkina, Taekyu Kang, Tyler Hilsabeck – I’ve learned from each of you, appreciated
your feedback and help, and am grateful for your friendship and time. Chandani Limbad and Okhee
Jeon, thank you for your mentorship and guidance at the beginning of my PhD, when I was quite new to
molecular biology – you helped me plan and execute some of my first biology experiments.
Finally, I would like to thank the funding sources that supported my Doctoral studies: NIH U54
AG075932, UG3 CA268105, AG055822, AG061879 (SM), AG051129 (SM).
In the following work, Chapter 1 includes sections of the journal articles entitled “Senolysis
induced by 25-hydroxycholesterol targets CRYAB in multiple cell types”, published in iScience, February
iv

18, 2022, by Limbad, C, Doi, R, Mcgirr, J, Ciotlos, S, Perez, K, Clayton, ZS, Daya, R, Seals, DR Campisi, J,
Melov, S. Also contained in this chapter are sections from the article “Single nuclei profiling identifies
cell specific markers of skeletal muscle aging, frailty, and senescence”, published in Aging, December 13,
2022, by Perez K, Ciotlos S, McGirr J, Limbad C, Doi R, Nederveen JP, Nilsson MI, Winer DA, Evans W,
Tarnopolsky M, Campisi J, Melov S. Finally, a modified version of the article “A Single Cell Time Course of
Senescence Uncovers Discrete Cell Trajectories and Transcriptional Heterogeneity”, currently available
as a preprint on bioRxiv, submitted February 18
th
, 2023, by Ciotlos S, Wimer L, Campisi J, Melov S.

v

TABLE OF CONTENTS
ACKNOWLEDGEMENTS .................................................................................................................. ii
LIST OF TABLES ............................................................................................................................. viii
LIST OF FIGURES ............................................................................................................................ ix
ABSTRACT ........................................................................................................................................ x

CHAPTER 1: INTRODUCTION .......................................................................................................... 1
Mechanisms of Senescence Induction .................................................................................... 1
Senescence initiation, molecular markers, and associated morphological changes ............. 2
The role of SnCs in healthspan and lifespan .......................................................................... 6
The many faces of senescence in different tissues ................................................................. 7
Targeting SnCs in wild type and transgenic mouse models of senescence .......................... 10
Profiling cardiac aging with strain echocardiography ........................................................... 12
Measuring aortic health and aging ....................................................................................... 13
CHAPTER 2: MATERIALS AND METHODS ..................................................................................... 15
In-vivo studies ....................................................................................................................... 15
In-vivo single-cell RNA seq experiments ............................................................................... 15
In-vivo study of effects of 25HC on Doxo-treated mice........................................................ 15
In-vivo testing of 25HC on aged mice ................................................................................... 16
Isolation of FAPs and SCs for single cell RNA-seq ................................................................. 16
FAPs and satellite cells for primary culture .......................................................................... 17
mDFs ..................................................................................................................................... 17
10X single-cell RNA-seq library prep ..................................................................................... 18
10X single-cell RNA-seq analysis ........................................................................................... 18
Cell culture ............................................................................................................................ 19
FAP cells, SCs, and HSMM ..................................................................................................... 19
Human cells ........................................................................................................................... 19
Quantitative RT-PCR (qRT-PCR) ............................................................................................ 19
shRNA .................................................................................................................................... 20
vi

Cell viability assay ................................................................................................................. 20
SA-β-gal assay ....................................................................................................................... 21
Quantification and statistical analysis .................................................................................. 21
RNA extractions from muscle biopsies for bulk analysis ...................................................... 21
Single nuclei RNA extraction and processing of human muscle biopsies ............................. 21
Bulk RNAseq data analysis of human muscle biopsies ......................................................... 22
Single-nuclei 5’ RNAseq data analysis of human muscle biopsies ........................................ 22
Pathway analysis of human muscle biopsies ........................................................................ 23
Spatial transcriptomics of human muscle biopsies ............................................................... 23
Cell culture of HSMMs .......................................................................................................... 24
qRT-PCR of HSMMs ............................................................................................................... 24
IMR90 culture and senescence induction for time course ................................................... 24
10x Single cell library preparation and sequencing for time course .................................... 24
Single cell data analysis for time course ............................................................................... 25
Strain echocardiography measurements .............................................................................. 25
Pulse wave velocity measurements of in vivo aortic stiffness .............................................. 26

CHAPTER 3: A SINGLE CELL TIME COURSE OF SENESCENCE UNCOVERS DISCRETE CELL
TRAJECTORIES AND TRANSCRIPTIONAL HETEROGENEITY .................................................. 27
Abstract ................................................................................................................................. 27
Background ........................................................................................................................... 28
Results ................................................................................................................................... 30
A population of proliferating cells evades acute senescence induction............................... 33
Differentiating between proliferating, senescent, and quiescent cells ................................ 36
Single cell transcriptomic time-course of senescence .......................................................... 42
Validating markers of proliferation and senescence across transcriptomic time course .... 46
Identification of cell trajectories across early, middle, and late phases of senescence ....... 48
Discussion ............................................................................................................................. 64

vii

CHAPTER 4: TREATMENT WITH THE SENOLYTIC 25HC IN OLD ANIMALS IMPROVES MEASURES OF
CARDIOVASCULAR AGING .................................................................................................... 68
Abstract ................................................................................................................................. 68
Results ................................................................................................................................... 69
Heart rate measurements for treated groups in cohort AI3 ................................................ 70
Ejection fraction and GLS improve with senolytic treatments ............................................. 71
Aged animals show higher PWV compared to young control 3MR females ........................ 72
Improvements in aortic health with senolytic treatment..................................................... 72
AI3 cohorts show improvements in aortic health with senolytic treatment ........................ 73
25HC treatment in old 3MR male animals reduces aortic elastic modulus .......................... 73
Discussion ............................................................................................................................. 74

CHAPTER 5: CONCLUSIONS ........................................................................................................... 75

REFERENCES .................................................................................................................................. 77

APPENDICES ................................................................................................................................... 77
APPENDIX A: Seurat_Workflow.R ......................................................................................... 81
APPENDIX B: Seurat_to_Slingshot.R ..................................................................................... 92
APPENDIX C: Slingshot_Tradeseq_workflow.R ..................................................................... 99
APPENDIX D: Analyze_SeuratDEGs_byLineage.R ................................................................ 104

viii

LIST OF TABLES
Table 1.1. Top differentially expressed genes between senescence resistant cells ......................................35
Table 1.2. Pathway analysis of senescence-resistant cells compared to senescent cells .............................36
Table 1.3. Upregulated (UR) lineage markers of early-stage senescence (Pro-D2) ......................................53
Table 1.4. Downregulated (DR) lineage markers of early-stage senescence (Pro-D2) .................................54
Table 1.5. Upregulated (UR) lineage markers of mid-stage senescence (D2-8) ...........................................57
Table 1.6. Downregulated (DR) lineage markers of mid-stage senescence (D2-D8) ....................................58
Table 1.7. Upregulated (UR) lineage markers of late-stage senescence (D8-D12) .......................................61
Table 1.8. Downregulated (DR) lineage markers of late-stage senescence (D8-D12) ..................................62

ix

LIST OF FIGURES
Figure 3.1. Experimental overview and quality control of replicate libraries................................................32
Figure 3.2. Visualization of proliferating cells ................................................................................................35
Figure 3.3. Proliferating, quiescent, and senescent cells form heterogeneous ............................................38
Figure 3.4. Visualizations of senescence time-course ...................................................................................44
Figure 3.5. Normalized expression of markers of proliferation and senescence ..........................................46
Figure 3.6. Senescing cells form distinct transcriptional lineages .................................................................49
Figure 3.7: Differentially regulated pathways in senescence lineages ..........................................................63
Figure 4.1 – Heart rate measurements for treated groups in cohort AI3 ......................................................70
Figure 4.2 – Ejection fraction and GLS measurements improve with senolytic treatments .........................71
Figure 4.3 – Aged animals show higher PWV compared to young control 3MR females .............................72
Figure 4.4 – Longitudinal comparisons show aortic health trend with senolytic treatment .........................72
Figure 4.5 – Longitudinal comparisons show improvements in aortic health with senolytic treatment .....73
Figure 4.6 – 25HC treatment in old 3MR male animals reduces aortic elastic modulus ...............................73

x

ABSTRACT
As one of the main hallmarks of aging, SnCs (SnCs) play a complex role in the aging process. These
permanently grow-arrested cells prevent can suppress the development of cancer; however they also
contribute to age-related tissue dysfunction and chronic inflammation, which are both associated with
many age-related diseases. SnCs are typically studied as endpoints of a complex transformational
process, owing to their frequent maladaptive effects on surrounding tissue and cells. SnCs accumulate
with age, and while they ultimately comprise a small percentage of cells in tissues, they have important
roles in age associated pathologies. Several obstacles remain in understanding the heterogeneous
nature of senescence and formulating potent intervention strategies. One approach targets SnCs and
kills them (“senolytic” approach) and is often driven by a low resolution understanding of SnC identity,
which risks both incomplete clearance and off-target effects. However, any method aiming to efficiently
manipulate SnCs in vivo will need a more accurate targeting method than what is commonly available or
studied, and in parallel a better understanding of cell history.
Cellular senescence is not a singular binary response, but a suite of response trajectories that vary by
multiple parameters including inducer (genotoxic or otherwise acute vs. replicative or other endogenous
stressor) and initial cell state. The heterogeneous nature of SnCs also varies by cell type and species, as
found in our recent publications on 1) the discovery of CRYAB aggregates as a marker of SnCs and the
development of 25HC as a senolytic from associated mouse single cell RNA and protein data, and 2) the
role of aged-associated SnCs in human sarcopenia, assayed at the single nucleus, bulk, and spatial
transcriptomic level. To expand our temporal understanding of senescence and elucidate the
developmental trajectories of SnCs, we performed single-cell RNA sequencing on IMR90 lung fibroblasts
senescencing across a 12-day period. Our analysis revealed substantial heterogeneity in gene expression
within timepoints and across the full time-course. Supervised trajectory inference of the time-course
data uncovered the root-origin and fates of distinct SnC lineages over 3 stages of senescence
xi

development. End and divergence points (clusters) of lineages were compared in terms of regulatory
activity and pathway activation, revealing that senescing lung fibroblasts specialize into different states
or types of senescence. We found that lineages had distinct regulatory and pathway activation profiles;
SnC lineages expressing wound healing signaling, mitochondrial dysfunction, and differential cytokine
signaling can be induced even in a simple in vitro context. These findings confirm the importance of
deeply profiling the nature of SnCs in the relevant biological context of interest, to avoid the paired
threats of poor efficiency and off target effects. Altogether our data provide a novel approach to study
SnC development, identifying SnC subtypes (“senotypes”), and differentiating between SnCs and
quiescent cells (Chapter 3). Applying this approach to known senescence contexts of interest will aid in
identifying key targets and areas for therapeutic intervention in pathologies of aging caused by SnCs.
The cardiovascular impacts of senolytic treatment in old animals, as measured by strain
echocardiography for cardiac health and pulse wave velocity (PWV) for measures of aortic health are
detailed in Chapter 4.

1

CHAPTER 1: INTRODUCTION
In a 1957 publication, the evolutionary biologist George Williams wrote, “It is indeed remarkable that
after a seemingly miraculous feat of morphogenesis, a complex metazoan should be unable to perform
the much simpler task of merely maintaining what is already formed.” [1] We understand this inability to
maintain as aging, manifesting as a gradual impairment of biological function, integrity, and
reproduction, coupled with increased vulnerability to death. It was perhaps premature or dismissive to
label the task of biological maintenance as simple. In the half century that has passed, only more
sensitive technologies for nucleic, proteomic, and organismal-level interrogations have permitted the
field of aging biology to better understand the multi-faceted complexity of aging. In a 2013 review,
Carlos López-Otín et. al. structured our knowledge of aging biology into 9 hallmarks, which pertain to
either genomic, proteomic, metabolic, or cell state deregulations (all interplaying and overlapping) [2].
One such hallmark is cellular senescence, a complex cell state which arises as a result of exogenous and
endogenous stressors, broadly characterized by cell cycle withdrawal, deregulated metabolism,
macromolecular damage, and an altered secretory phenotype (compared to cycling or quiescent cells)
[3-6]. SnCs (SnCs) are implicated in various physiological processes (both adaptive and maladaptive) and
diseases of aging, including Alzheimer’s and Parkinson’s disease, atherosclerosis and cardiovascular
dysfunction, chemotherapy cardiotoxicity, osteoarthritis and osteoporosis; more benign processes
include wound healing, tissue regeneration, embryogenesis, and tumor suppression [7-20].
Mechanisms of Senescence Induction
Originally described by Hayflick et. al. in 1961, cellular senescence was first observed in normal human
fetal fibroblasts (FBs). Cells were observed to cease proliferating after a limited number of divisions in
culture. [21]. While this initial finding was attributed to telomere shortening, the past decades of
research have expanded both our knowledge of senescence inducers as well as the resulting cellular
2

adaptations in this state. Telomeric degeneration occurs as a natural outcome of DNA replication during
proliferation. As a result of incomplete duplication during cell division, these chromosomal regions
shorten with each division (in non-germline cells). Upon reaching a critically short length, the truncated
telomeres are recognized by the cell as DNA damage. “Replicative senescence” is triggered in this event
[5, 22-24]. In contrast, premature or “stress-induced” senescence is triggered by other inducers
including non-telomeric DNA damage, genotoxic drugs such as the chemotherapy drug Doxorubicin
(DOXO), irradiation, oncogenic stress and tumor suppressor loss, oxidative stress, and mitochondrial
dysfunction. The inducer of choice for this proposal is DOXO, an effective cancer drug that has negative
side-effects including cardiotoxicity and senescence induction. In use since 1974, it damages DNA
through several mechanisms including direct intercalation into DNA (blocking synthesis along the
strands), free radical generation, and inhibition of topoisomerase II leading to double strand breaks
(DSBs) [25-29].
Senescence initiation, molecular markers, and associated morphological changes
Two main genetic pathways govern the entrance of a cell into senescence: the p53/p21 (p21/CDKN1A)
and p16Ink4a/retinoblastoma (p16/CDKN2A) protein pathways. Across species and cell types, there is
variability in the main driving pathways (with great relevance to therapeutic approaches); most
senescence inducers activate either or both pathways [30]. For example, telomeric/replicative
senescence in human cells can be modulated by either p16 or p21 signaling, whereas p21 signaling is
preferred in murine cells. Mitochondrial dysfunction and other types of intrinsic cellular stress generally
activate p16 signaling; and DNA damage (including but not limited to telomere attrition) primarily
triggers p21 [31-34]. In cycling cells, retinoblastoma family proteins including RB1, p107/RBL1, and
p130/RBL2 are phosphorylated by cyclin-dependent kinases (CDK) CDK2, 4, and 6. In their
phosphorylated forms, RB proteins are not capable of repressing the E2F transcription factors, which
regulate mitotic genes involved in cell cycle progression from G1 to S-phase [35, 36]. SnCs upregulate
3

and accumulate p21 and/or p16; these in turn inhibit CDK2 and CDK4/6 respectively. In turn RB is
dephosphorylated, begins to inhibit the activity of E2F transcriptions factors, and leads to cell cycle
arrest that is irreversible – even if p53 and p16 are later downregulated [32]. Alongside undergoing cell
cycle arrest, SnCs increase in size and attain a flattened morphology [37], and express senescence-
associated β-galactosidase (SA-βgal; in part due to an underlying increase in lysosomal mass) [38, 39].
SnCs also acquire a senescence-associated secretory phenotype (SASP), which consists of chemokines,
cytokines, growth factors, and matrix remodeling proteins [40-42]. The SASP is a pro-inflammatory
profile, and leads to alterations in tissue microenvironment and interestingly, paracrine transmission of
cell senescence [6, 43-45]. Major components of the SASP include interleukins IL-1, 6, and 8, vascular
endothelial growth factor A (VEGFA), and matrix metalloproteinases such as MMP3. Proinflammatory
cytokine components of the SASP, such as IL6, IL8, IL1α/β, TNFα, and TNFSF7 are generally conserved
across senescence contexts. However, as with senescence markers in general, the composition of the
SASP is fluid, and varies by induction method, cell type, and organism [40, 42, 46-49]. Inducers of
senescence cause permanent genomic damage (DSBs being particularly noxious), activating DNA
damage response (DDR) signaling and leading to DNA damage foci. These foci are initiated by the first
responder ataxia telangiectasia mutated (ATM) gene, which encodes a phosphoinositide-3 kinase. Its
downstream activity recruits p53 binding protein 1 (TP53BP1; involved in repair and checkpoint
activation) and γH2AX, the phosphorylated form of the H2AX histone protein that forms in the DDR
context, to damaged sites. Other important ATM targets include the p53 tumor suppressor, checkpoint
kinase-2 (CHK2), a DDR kinase that promotes growth arrest, and the MRE11-RAD51-NBS1 complex
whose activity ranges from initial DSB processing to repair by non-homologous end joining [6, 32, 50].
Irreparable DSBs lead to the persistent DDR signaling and DNA damage foci, termed DNA segments with
chromatin alterations reinforcing senescence (DNA-SCARs), observed in senescence. Other genomic
alterations include senescence-associated heterochromatin foci (SAHFs) and distension satellites (SADs).
4

The former is characterized by DAPI-dense foci and higher density of nuclear pores, as well as
enrichment of heterochromatin factors H3K9me3, heterochromatin protein 1 and HMGA proteins, and
histone variant macroH2A. Despite a readily visible phenotype, the functional significance is not fully
understood, with data suggesting a limited impact on gene expression [51-53]. SADs are large (up to 3
megabase) genomic reshaping events, where centromeric constitutive heterochromatin becomes
unfolded and decompressed. This process occurs early in the transition to senescence and is well
conserved across mouse and human cell types, making SADs a more reliable marker than SAHFs.
Retrotransposable elements, generally repressed in constitutive heterochromatin, are opened to
transcription and activate the cGAS-STING (a SASP regulator) mediated type-1 interferon response [54-
56]. The nuclear protein lamin b1 (LMNB1) is downregulated in mouse and human SnCs, both in vitro
and in vivo. Loss of expression is independent of ATM and ROS signaling but does occur under both
p53/p21 or p16 signaling. Its conservation across experimental contexts and inducers makes it a robust
marker of senescence [57]. Another important marker of senescence is the loss of nuclear HMGB1, a
protein that facilitates transcription factor (including p53) access to promoters through direct
conformational alterations to DNA [58-60]. It has been validated in human and mouse cells in vitro, and
increased circulation has also been observed in mice. As both an alarmin and a component of damage
associated molecular patterns (DAMPs), HMGB1 normally resides intracellularly and is secreted out of
the cell in stressed states. Upon secretion, it binds to cell surface toll-like receptors (TLRs) and receptors
for advanced glycation end-products (RAGEs), provoking immune and inflammatory responses [61-63]. It
is an early responder to senescence induction, being secreted as quickly as 1 day after x-ray irradiation.
Finally, the secretion of HMGB1 is required for robust secretion SASP factors such as IL-6 and MMP-3
[59]. Figure 1 provides a summary of these hallmarks of senescence [64].

5

To visualize and eliminate SnCs in living mice, the P16-3MR transgenic mouse strain (3MR) was created
[15]. While no universal or specific markers of senescence exist, most cells do initiate their transition to
senescence through the p21 or p16 pathways, as previously mentioned. The 3MR model was created by
insertion of a BAC containing a trimodal reporter fusion protein driven by a copy of the native p16
promoter. Sequencing and analysis performed by myself and our collaborators indicates that this
transgene was inserted into chromosome 13, in an intronic region of a poorly characterized gene. Data
also indicates that two copies of the reporter segment are present on each homolog of chromosome 13.
The approach taken allowed the animal’s genome to main a normal diploid locus for the CDKN2A locus
and its transcript, instead only generating a novel multifunctional protein. The reporter protein is
composed of the functional domains of a synthetic Renilla luciferase (LUC), monomeric red fluorescent
protein, and truncated herpes simplex virus 1 thymidine kinase. Expression of this transcript indicates
activation of the p16 locus and a cellular transition to senescence. The elements of the reporter function

Figure 1. Summary of senescence hallmarks. Cellular senescence is caused by both endogenous (telomere degeneration,
oncogenic stress) and exogenous (genotoxic drugs, irradiation) stressors. Normally cycling or differentiating cells, as well as
transiently quiescent cells, undergo permanent growth arrest upon activation of either or both of the p53/p21 and p16/Rb
pathways. SnCs exhibit a flattened and enlarged morphology, loss of nuclear Lamin B1, an altered secretory profile (SASP),
increased SA-βGal expression, chromatin remodeling events including SAHFs and SADs, altered mitochondrial function,
resistance to apoptosis, and secretion of nuclear HMGB1.
Adapted from Herranz and Gil, 2018

Cycling cell Transition to growth arrest Senescent cell
Senescence induction:
Telomere attrition/replicative stress
Genotoxic drugs
Irradiation
Oncogenic stress
Secretion of HMGB1
6

as follows: LUC allows detection of cells expressing the 3MR transcript in vivo, mRFP permits FACS
sorting of cells from isolate tissues, and HSV-TK allows targeted killing of SnCs by treatment with
Ganciclovir (GCV; HSV-TK converts GCV into a DNA chain terminator, fragmenting mitochondrial DNA
and leading to death by apoptosis) [65, 66]. One drawback to this approach is the potential to miss cells
that senescence through the p21 pathway; work is underway to design a new generation of transgenic
mice which can capture this modality as well. In addition, treatment with DOXO in the current model
leads to senescence induction throughout the entire organism; future mouse models will include tissue-
specific capabilities, which will greatly increase the accuracy of studies. Limitations aside, this powerful
tool has pushed the boundaries of senescence biology and led to the first translation of the field’s
discoveries [13] into human trials for senolytics at Unity Biotechnology.
The role of SnCs in healthspan and lifespan
SnCs play a significant role in the aging process, undergoing alterations that can contribute to various
age-related diseases and organismal health decline. While SnCs can have some beneficial roles, such as
wound healing and tissue repair, their accumulation in tissues over time can negatively impact the aging
process in several ways. SnCs can release pro-inflammatory molecules, known as the senescence-
associated secretory phenotype (SASP). This promotes a low-level chronic inflammation, which
contributes to age-related diseases, such as atherosclerosis, osteoarthritis, and neurodegenerative
disorders. As SnCs accumulate in tissues, they can disrupt the normal functioning of healthy neighboring
cells. This can lead to a decline in tissue repair and regeneration, contributing to frailty and organ
dysfunction. SnCs can interfere with cell signaling, which is essential for maintaining proper tissue
function. This disruption can impact processes like immune system function, cell growth, and tissue
repair. SnCs play an important role in both healthspan and lifespan. Healthspan refers to the period of
life during which an individual is generally healthy, free from serious chronic diseases, and able to
maintain functional independence. Lifespan, on the other hand, refers to the total length of time an
7

individual lives. SnCs can have both beneficial and detrimental effects on healthspan and lifespan,
depending on the context and balance between their accumulation and clearance. In recent years,
research has focused on developing strategies to target SnCs to improve both healthspan and lifespan.
These strategies involve selectively eliminating SnCs (senolysis) or modulating their secretory phenotype
(senomorphics). Preclinical studies in animal models have shown that targeting SnCs can improve
healthspan by delaying the onset of age-related diseases, preserving organ function, and extending
lifespan. Human clinical trials are ongoing to evaluate the potential of senotherapies in improving
healthspan and treating age-related diseases.
The many faces of senescence in different tissues
Although neurons are considered post-mitotic, meaning they do not divide, they can still exhibit some
features of senescence. Neurons may undergo functional decline and show altered gene expression
patterns, reduced synaptic plasticity, and a SASP-like secretory profile. Glial cells, which include
astrocytes, microglia, and oligodendrocytes, play crucial roles in supporting neurons and maintaining
brain homeostasis. Senescence in these cells can have significant consequences for brain health. For
example, senescent astrocytes may show reduced ability to support neurons and maintain the blood-
brain barrier, while senescent microglia may exhibit impaired ability to clear debris and produce an
exaggerated inflammatory response. In the brain, this can contribute to chronic low-grade
inflammation, also known as "inflammaging," which has been implicated in cognitive decline and
neurodegenerative diseases such as Alzheimer's and Parkinson's disease. Senescence in endothelial cells
that form the blood-brain barrier can lead to its dysfunction, allowing potentially harmful substances to
enter the brain and contribute to neuroinflammation and neurodegeneration. These mechanisms link
brain health and function to vascular and cardiac cell senescence.
8

Proper vascularization is also critical to muscle function and health. Sarcopenia is the age-related loss of
skeletal muscle mass, strength, and function, which can lead to frailty, increased risk of falls, and
reduced quality of life. Multiple factors contribute to sarcopenia, and cellular senescence is considered
one of the mechanisms involved in this process. Satellite cells are muscle stem cells that play a vital role
in muscle repair and regeneration. As we age, these cells can become senescent, resulting in a reduced
capacity for muscle repair and regeneration. Senescent satellite cells show decreased proliferative
potential and impaired ability to differentiate into functional muscle cells, leading to a decline in muscle
mass and function. In the context of muscle aging, the SASP can contribute to chronic inflammation and
promote a catabolic environment that favors muscle breakdown over synthesis. This can impair muscle
function and accelerate the progression of sarcopenia. Cellular senescence can disrupt the signaling
pathways that regulate muscle growth and maintenance. For example, senescence can negatively affect
the insulin-like growth factor 1 (IGF-1) and mammalian target of rapamycin (mTOR) signaling pathways,
both of which are essential for maintaining muscle mass and function. As discussed previously, vascular
senescence can impair blood flow and delivery of oxygen and nutrients to muscle tissue. This can
negatively impact muscle function and contribute to sarcopenia.
Vascular senescence refers to the age-related changes and functional decline in the cells that form the
vascular system, which includes blood vessels such as arteries, veins, and capillaries. Vascular
senescence can contribute to cardiovascular diseases, impaired tissue perfusion, and other age-related
disorders. Endothelial cells line the inner surface of blood vessels and play a critical role in maintaining
vascular homeostasis, regulating blood flow, and preventing clot formation. As endothelial cells age and
are exposed to cell stressors in the blood, they can become senescent, exhibiting altered gene
expression and impaired function. Senescent endothelial cells can contribute to atherosclerosis,
endothelial dysfunction, and increased inflammation. Vascular smooth muscle cells (VSMCs) are
responsible for the contractile function of blood vessels, which helps regulate blood pressure and blood
9

flow. Senescent VSMCs can exhibit reduced contractile function, and increased production of pro-
inflammatory molecules. This can lead to the development of vascular diseases, such as atherosclerosis
and arterial stiffness.
According to the American Heart Association, cardiovascular disease (CVD) is the leading cause of death
in the United States. By 2030, CVD (conditions including hypertension, coronary heart disease, heart
failure, stroke, etc.) is projected to reach 40% prevalence, with an economic toll of about $1 trillion
dollars. CVD encompasses a range of cardiac conditions (defective output, vessel occlusions and
ruptures, arrhythmia, inflammation) and underlying mechanisms [67]. Age is the main risk factor for
developing CVD [68-71]. While progress has been made in treating some aspects of CVD in recent
decades, the cellular and molecular mechanisms by which aging compromises CV function remain
underexplored. CV medications are a broadly acting group of compounds, including anticoagulants (ex.
Heparin – decreases the clotting ability of blood, used to prevent clotting during surgery or to treat
heart attacks and anginas), antiplatelet therapies (ex. Aspirin, Plavix – decrease platelet adhesion,
prescribed preventatively against plaque buildup), calcium channel blockers (ex. Amlodipine – blocks
calcium uptake by vascular smooth muscle cells/cardiomyocytes, relaxing vessels and inhibiting
contraction of heart muscle; prescribed for hypertension), cholesterol drugs (ex. Lovastatin – a statin
used to treat high LDL cholesterol, used to prevent atherosclerotic CVD), and vasodilators (ex.
Nitroglycerin – a nitrate that acts to relax blood vessels, rescuing impaired flow of blood and oxygen to
the heart). Much of the CV drug arsenal thus acts on vascular rigidity and blood flow, and indeed a key
process that links advancing age to CVD is the development of arterial dysfunction. With age, oxidative
stress and chronic inflammation cause a gradual impairment of endothelial cell function (reduced
dilation and arterial stiffening) [72-74]. Arterial dysfunction leads to ventricular hypertrophy, systolic
and diastolic dysfunction, atherosclerosis of the coronary artery, and increased heart attack risk [75-78].
10

It has been well established that with age, an increase in SnCs occurs in many human and rodent tissues,
including both vasculature and the heart. Furthermore, the ablation of these SnCs in both progeroid-
modeling and aged mice can significantly ameliorate pathologies of aging including cardiac dysfunction,
sarcopenia, cataracts, and bone degeneration. Apart from natural aging, short-term or acute stressors
can also cause DNA damage that results in senescence and SASP production. As previously mentioned,
one such stressor is DOXO, which is known to induce cardiac and vascular senescence and SASPs. As this
leads to an impaired endothelial functional state, reminiscent to changes seen with aging, it is
commonly used as a model of accelerated CV aging. The persistence of SnCs and the SASP can not only
impair tissue repair, but also lead to further damage and hyperplasia, through the chronic sterile low-
level inflammation which follows.
Targeting SnCs in wild type and transgenic mouse models of senescence
The p16-3MR transgenic mouse model is used to study the selective clearance of SnCs with the drug
ganciclovir, using the expression of p16 as a target. In this model, the p16INK4a promoter is used to
drive the expression of a fusion protein consisting of the herpes simplex virus thymidine kinase (HSV-
TK). This allows for the selective killing of cells that express p16, which are often SnCs that accumulate
with age and contribute to age-related diseases. Treatment with ganciclovir activates the HSV-TK
component of the fusion protein, leading to the death of SnCs that express the p16INK4a gene. This
allows researchers to study the effects of selectively removing SnCs on aging and age-related diseases,
as well as to develop potential therapeutic strategies for targeting SnCs. The HSV-TK component in the
p16-3MR mouse model is an enzyme that can phosphorylate (add a phosphate group) to a variety of
nucleoside analogues, including ganciclovir. Once phosphorylated, ganciclovir is converted into a toxic
compound that interferes with DNA replication and induces apoptosis (programmed cell death) in the
targeted cells. In the p16-3MR mouse model, the HSV-TK component is specifically expressed in cells
that have high levels of p16INK4a, which are often SnCs. Upon treatment with ganciclovir, only these
11

p16INK4a-positive SnCs are targeted for cell death, while normal, healthy cells are not affected. This
selective clearance of SnCs with ganciclovir has been shown to improve health and extend lifespan in
mice, suggesting that targeting SnCs may be a promising strategy for treating age-related diseases in
humans. However, there are also some limitations and potential drawbacks to their use. In some cases,
the expression of the transgene may not accurately reflect the normal pattern of gene expression seen
in humans or other animals. For example, the p16-3MR mouse model uses a modified form of the
p16INK4a gene promoter to drive the expression of the suicide gene. While this may allow for precise
targeting of SnCs, it may not fully reflect the natural regulation of the p16INK4a gene in vivo. The
expression of the transgene may have unintended effects on other cellular processes or functions,
leading to potential off-target effects. For example, the induction of apoptosis in SnCs may have effects
on neighboring cells or tissues that are not senescent, potentially leading to unintended consequences.
To kill SnCs in non-transgenic mice in my studies, the drug ABT-263 (Navitoclax) was used. It is a
senolytic drug that has been shown to selectively eliminate senescent cells. ABT-263 works by inhibiting
the activity of two proteins, Bcl-2 and Bcl-xL, which are members of the B-cell lymphoma 2 (Bcl-2) family
of proteins. These proteins play a key role in regulating programmed cell death, or apoptosis. In
senescent cells, Bcl-2 and Bcl-xL are overexpressed, which helps to protect these cells from apoptosis
and allows them to persist in tissues. By inhibiting the activity of Bcl-2 and Bcl-xL, ABT-263 triggers the
apoptotic pathway in senescent cells, leading to their selective elimination. This process is thought to be
mediated by the activation of pro-apoptotic proteins, such as Bim, which bind to and inhibit the activity
of Bcl-2 and Bcl-xL, thereby promoting cell death. It is important to note that ABT-263 has also been
shown to have some off-target effects and can cause toxicity in certain tissues, such as the platelets.
Therefore, further research is needed to optimize the use of senolytic drugs like ABT-263 and develop
more targeted approaches for eliminating senescent cells.

12

Profiling cardiac aging with strain echocardiography
Strain echocardiography is a technique used in cardiac imaging that allows for the assessment of
myocardial deformation and contractile function of the heart. This technique uses advanced software
algorithms to track the movement of myocardial tissue throughout the cardiac cycle and measure the
amount of deformation, or strain, that occurs in the tissue. Strain echocardiography can be performed
using either 2D or 3D echocardiography. In 2D echocardiography, a series of images are obtained from
different angles, and the software is used to track the movement of specific points on the myocardial
tissue between the images. In 3D echocardiography, a complete volumetric image of the heart is
obtained, and the software is used to track the movement of the entire myocardial tissue volume. Strain
echocardiography allows for the assessment of several different parameters of cardiac function. Global
longitudinal strain (GLS) is one such metric that measures the extent of longitudinal myocardial
deformation during systole, providing information on overall ventricular function. Strain rate measures
the velocity of myocardial deformation, providing information on the rate of contraction and relaxation
of the heart. Twist and torsion measure the rotation and twisting of the heart during the cardiac cycle,
providing information on ventricular function (such as ejection fraction and fractional shortening) and
the interaction between the left and right ventricles. Strain echocardiography can be used to assess
cardiac aging by measuring changes in myocardial deformation and contractile function over time. With
age, the heart undergoes structural and functional changes, including changes in myocardial stiffness,
alterations in the extracellular matrix, and decreased contractile function. These changes can be
detected and quantified using strain echocardiography. Several studies have demonstrated that strain
echocardiography can detect age-related changes in myocardial deformation and contractile function,
even in the absence of clinically apparent cardiac disease. For example, global longitudinal strain (GLS)
has been shown to decrease with age, indicating a reduction in overall ventricular function. In addition,
regional strain analysis has been used to identify specific regions of the heart that are most affected by
13

age-related changes, such as the basal segments of the left ventricle. Strain echocardiography can also
be used to assess the impact of age-related comorbidities on cardiac function, such as hypertension,
diabetes, and atherosclerosis. For example, hypertension can lead to myocardial hypertrophy and
fibrosis, which can be detected using strain echocardiography. Diabetes has been associated with
impaired GLS and abnormalities in twist and torsion, which can also be detected using strain
echocardiography. Overall, strain echocardiography provides a sensitive and non-invasive tool for the
assessment of cardiac aging and the impact of age-related comorbidities on cardiac function. This
information can be used to identify individuals at risk for cardiovascular disease and to guide the
development of targeted interventions to prevent or treat age-related cardiac dysfunction.
Measuring aortic health and aging
The aorta is the largest artery in the body and is responsible for carrying oxygen-rich blood from the
heart to the rest of the body. Aortic stiffness refers to the loss of elasticity in the aortic wall, which
results in reduced compliance and increased resistance to blood flow. This can lead to increased
workload on the heart, increased blood pressure, and an increased risk of cardiovascular disease. Aortic
health changes with age, as the aorta undergoes structural and functional changes over time. The aorta
becomes stiffer with age due to changes in the composition of the extracellular matrix, increased
deposition of collagen and elastin, and increased cross-linking of these proteins. These changes can
result in increased pulse wave velocity, decreased compliance, and increased systolic blood pressure.
Cellular senescence may also play a role in aortic health. Senescent cells accumulate in the vasculature
with age and secrete pro-inflammatory cytokines, growth factors, and extracellular matrix proteins,
which can promote aortic stiffness and contribute to the development of cardiovascular disease.
Senescent cells have been found to be present in aortic tissue from aged mice and humans, and their
elimination has been shown to improve aortic compliance and reduce inflammation in mouse models of
aging. In addition to age and cellular senescence, other factors can contribute to aortic stiffness,
14

including hypertension, diabetes, and atherosclerosis. These conditions can accelerate the development
of aortic stiffness and increase the risk of cardiovascular disease. Overall, aortic stiffness is an important
marker of cardiovascular health and is associated with an increased risk of morbidity and mortality.
Strategies to prevent or treat aortic stiffness, such as lifestyle interventions, pharmacotherapy, or
senolytic therapies, may help to reduce the burden of cardiovascular disease in aging populations. Aortic
pulse wave velocity (PWV) is a measure of aortic stiffness that reflects the speed at which pressure
waves propagate through the arterial tree. It is defined as the distance between two points along the
aorta divided by the time it takes for the pressure wave to travel between those points. Aortic PWV is
typically measured by recording pressure waveforms at two locations along the aorta and determining
the time delay between the waveforms. Aortic PWV increases with age, reflecting the progressive
stiffening of the arterial wall that occurs with aging. In healthy individuals, aortic PWV increases by
about 0.5-1 m/s per decade after age 40, although this rate of increase may be accelerated in the
presence of cardiovascular risk factors or disease. Doppler ultrasound is a non-invasive technique used
to measure aortic PWV. This technique involves placing two ultrasound probes at different locations
along the aorta, typically at the carotid and femoral arteries, and measuring the time delay between the
arrival of the pressure wave at each location. The distance between the two probes can be measured
using a tape measure or a body surface area calculation. The Doppler signal is used to measure the
velocity of blood flow at each location and to calculate the transit time of the pressure wave between
the two sites. The aortic PWV is then calculated by dividing the distance between the two probes by the
transit time. This technique is safe, non-invasive, and can be performed quickly and easily in a clinical
setting. Measurement of aortic PWV provides valuable information about cardiovascular health and is a
predictor of future cardiovascular events. It is often used in research and clinical settings to evaluate the
effectiveness of interventions designed to reduce aortic stiffness, such as lifestyle modifications,
pharmacotherapy, or exercise training.
15

CHAPTER 2: MATERIALS AND METHODS
In-vivo studies
All mice were maintained according to National Institutes of Health guidelines for use of live animals. All
experimental procedures were approved by the Institutional Animal Care and Use Committee at the
Buck Institute (A10203). Gender, genetic background and age of mice are described for specific
experiments below.
In-vivo single-cell RNA seq experiments
Four-month-old male C57BL/6 wildtype (WT) mice were purchased from the Jackson Laboratory, and
injected once intraperitoneally (i.p.) with 10 mg/kg doxorubicin hydrochloride (Doxo) (Tocris Bioscience,
Product# 2252) in endotoxin-free Dulbecco′s PBS (w/o Ca++ and Mg++) (ETF-PBS) (EMD Millipore,
Product TMS-012-A) or only ETF-PBS. After 6 days, the mice were treated with ABT263 in 10% ethanol,
30% PEG400, 60% Phosal 50 PG (vehicle) or vehicle only (Chang et al., 2016). ABT263 or vehicle was
administered by gavage at 50 mg/kg/day for 5 consecutive days per cycle for 2 cycles with 9 days
between cycles. At 14 days after the end of the 2nd ABT263 or vehicle treatment cycle, the mice were
euthanized and skeletal muscles were harvested.
In-vivo study of effects of 25HC on Doxo-treated mice
Male C57BL/6J mice (3–4 months old) were purchased from the Jackson Laboratory. We performed five
i.p. injections with 4 mg/kg of Doxo (Tocris Bioscience, Product# 2252) for each injection over 10 days
(Figure 6A). Therefore, we injected a total of 20 mg/kg of Doxo/mouse over 10 days. We used EFT-PBS
(EMD Millipore, Product# TMS-012-A) as a control. To test effects of 25HC (Tocris Bioscience, Product#
5741) on doxo-treated mice, we prepared three experimental groups (n=15 for each group): 1)
PBS+HβCD vehicle (22.5%), 2) Doxo (20 mg/kg)+vehicle (22.5%), 3) Doxo (20 mg/kg)+25HC (50 mg/kg)
16

(Figure 6A); 25HC was dissolved in 22.5% 2-hydroxypropyl-beta-cyclodextrin (HβCD) (Tocris Bioscience,
Product 0708). On the Day 5 after the last Doxo/PBS injections, the mice were treated with either 25HC
(50 mg/kg) or vehicle by i.p. injections for 7 consecutive days (Figure 6A). 12 hr after the first 25HC or
vehicle treatments, urine was collected (Figure 6A) from each group to test for the senolysis marker 15-
d-PGJ2. 1 wk after the last day of 25HC or vehicle treatments, the mice were euthanized over the period
of 3 days (Figure 6A), and tissues were harvested for gene expression analysis. Senescence markers
p16INK4a, p21Cip1 and CRYAB expression was analyzed on gastrocnemius (Gastroc), soleus (So), tibialis
anterior (TA), visceral fat (VF), kidney (Ki), liver (Li), lung (Lu), heart (He), and skin (Sk) samples by qRT-
PCR.
In-vivo testing of 25HC on aged mice
Two groups of male and female C57BL/6J mice were used. The young group was 5 months old, and the
aged group was 24–26 months old. Six groups were prepared, 3 for male (n=9–15 mice) and 3 for female
(n=10-16 mice/group) mice: 1) Young (5 months) + vehicle, 2) Old (24–26 months) + vehicle, 3) Old (24–
26 months) + 25HC (Figure 7A). Each group was treated with either 25HC (50 mg/kg) or vehicle via i.p.
injections for 5 consecutive days. After 12 h of the first treatments, urine samples were collected in low
protein binding tubes (Figure 7A) to test for the senolysis marker 15-d-PGJ2. Mice were euthanized 2
days after the last injections over a period of 3 days (Figure 7A), and tissues were harvested for gene
expression analysis. Senescence markers and CRYAB were analyzed using mRNA from gastrocnemius
(Gastroc), soleus (So), tibialis anterior (TA), visceral fat (VF), kidney (Ki), liver (Li), lung (Lu), heart (He),
and skin (Sk) samples using qRT-PCR.
Isolation of FAPs and SCs for single cell RNA-seq
Skeletal muscles were harvested from hind limbs of mice from PBS, Doxo or Doxo+ABT groups. Skeletal
muscles from 5 mice were combined for each group. For skeletal muscle dissociation, a kit from Miltenyi
17

Biotech was used (cat#130-098-305). Muscles were chopped into 2-4 mm pieces and enzymatic
digestion was performed as described in the product datasheet using recommended kit enzymes, C
Tubes, (Company: Miltenyi, Cat# 130-096-334), and gentle MACS Octo Dissociator with Heaters (cat#
130-096-427). After digestion, Debris Removal Solution (Company: Miltenyi, cat#130-109-398) was used
to achieve cleaner single-cell suspensions. The suspension was stained with APC-CD31 (Biolegend, cat#
102514), APC-CD45 (Biolegend, cat# 103112), PB-Sca1(Invitrogen, cat# 12-5981-82), FITC-Integrin α7
(MBL life sciences, cat#K0046-4) and Biotin-PDGFRa(Thermo fisher, cat# 13-1401-82), PE/Cy7-
Streptavidin (Thermo fisher, cat# SA1012) for FACS (FACS Aria II, BD Biosciences). FAPs were identified
as Sca1+, PDGFra+, CD31-, CD45- and a7-integrin- and SCs were identified as a7-integrin+, Sca1-,
PDGFra-, CD31-, and CD45- (Joe et al., 2010; Liu et al., 2015).
FAPs and satellite cells for primary culture
Hind limb muscles were digested in DMEM containing 0.5% (w/v) collagenase type II (Worthington) for
90 min at 37 °C with trituration and passed through a 40-μm nylon mesh. Erythrocytes were eliminated
by treating with BD Pharm Lyse (BD Biosciences). Cells were resuspended in a washing buffer of PBS
with 2% FBS and stained with antibodies for 30 min at 4 °C as described (Kamizaki et al., 2017). FAPs and
SCs were isolated from the suspension by FACS (Aria II, BD Biosciences) using anti-CD31, CD45, Sca-1,
PDGFR and VCAM1 antibodies. Antibodies are listed in Table S6. FAPs were identified as Sca1/PDGFR-
positive and CD31-, CD45- and VCAM1-negative cells (Uezumi et al., 2010). SCs were identified as
VCAM1-positive and CD31-, CD45-, Sca-1 and PDGFR-negative cells (Liu et al., 2015).
mDFs
For skin sample collection, mice were euthanized, hair removed from dorsal skin and the area cleaned
with alcohol. A small 1 cm x 1 cm region of skin was removed, and subcutaneous fat removed by scalpel.
Skin was digested in 2–3 ml of 0.25% trypsin at 37º C for 1 h. Hair and epidermis were further removed,
18

and the dermis was minced into 2–3 mm pieces. Digestion of the sample in 5 ml of 1 mg/ml collagenase
(Sigma, Cat#C5138) in HBSS (with Ca++ and Mg++) (Gibco, cat#24020117) was carried out at 370 C on a
shaker for 1 h. Collagenase digestion was stopped by addition of 10 ml 10% FBS DMEM. Filtration of the
homogenate through a 70-micron mesh was performed to remove undigested tissue. Differential
centrifugation and washing were performed at 1200 rpm for 5 min, and the pellet was resuspended in
10% FBS DMEM then seeded in a six-well plate. Media were changed at 5 days after seeding and then
every 3 days.
10X single-cell RNA-seq library prep
FAPs and SCs were isolated by FACS. Samples were prepared for 10X single-cell RNAseq to target 10 K
cells. The 3′ single-cell library preparation was performed by following the 10X kit protocol. QC for cDNA
and libraries were performed using Tapestation, according to the 10X protocol. Sequencing of the
resultant libraries was then performed on Illumina sequencers as per the manufacturer’s protocols.
10X single-cell RNA-seq analysis
Libraries were processed through CellRanger v2.1.0 [1], aligned to mouse genome version mm10. Data
analysis was performed with Seurat version 2.3(Stuart et al., 2019) by integrating libraries from control,
DOXO treated, and DOXO+ABT263 treated animals. Differential gene expression analysis across clusters
was performed between integrated library pairs (Doxo vs. PBS, Doxo+ABT vs Doxo). Clusters with
significant (p < 0.05) expression changes in senescence marker genes between treatment groups were
further analyzed.

19

Cell culture – 25HC
All experiments were performed in incubators at 37 °C, 3% O2 and 5% CO2. Senescence was induced by
treatment with either Doxo (250 Nm for 24 h) (Demaria et al., 2017) or X-irradiation (10 Gy) (Limbad et
al., 2020). For cell culture studies using 25HC, it was used at the indicated concentrations in the media.
25HC was dissolved in 100% EtOH, such that the total concentration of EtOH in the media did not
exceed 0.5%, and vehicle controls used were equivalent volumes of EtOH without 25HC.
FAP cells, SCs, and HSMM
FAPs and SCs were cultured on Matrigel-coated dishes (BD Biosciences) in growth medium (GM)
consisting of DMEM plus 10% FBS, 1% penicillin-streptomycin and 10% MyoCult (STEMCELL
Technologies) for 7 days. Human myogenic cell HSMMs (Lonza) were maintained according to the
manufacturer’s instructions with SkGM-2 medium (Lonza).
Human cells - 25HC
IMR-90 cells were cultured in DMEM with 10% FBS. Cardiac endothelial cells (Sciencell, Catalog 6000,
Lot 15367), liver stellate cells (Sciencell, Catalog 5300, Lot 25907), renal proximal tubule epithelial cells
(Sciencell, Catalog 4100, Lot 19870) and articular chondrocytes (Sciencell, Catalog 4650, Lot 69440) were
cultured in media from Sciencell as per the datasheets.
Quantitative RT-PCR (qRT-PCR)
Total RNA was extracted using either Isolate II RNA mini kit (Bio-52073) or RNeasy Mini Kit (Qiagen),
then reverse transcribed using PrimeScript RT Reagent Kit (Perfect Real Time) (Takara Bio USA, Prod
RR037B). qRT-PCR reactions were performed using TB Green Premix Ex Taq (Tli RNase H Plus) (Takara
Bio USA, Prod RR420L). The amounts of mRNA were normalized relative to those of Actin mRNA. Specific
primer sequences used for qRT-PCR are listed in Table S8.
20

shRNA
Experiments were performed in triplicate. Cells were seeded (2000 cells/well) in 96-well plates. 24 h
later, cells were infected with either a scrambled sequence (control), shCRYAB or shHMOX1 lentiviruses
at 5 MOI in culture medium with 0.1% polybrene. Control and shRNAs were from Mission SHRNA
Custom Lentiviral Particles. Each vial contained 106 TU/mL. For HMOX1, TRCN0000218100 TRC2 and
TRCN0000234075 TRC2 were used. For CRYAB, TRCN0000097213 TRC1 and TRCN0000097210 TRC1
were used. SHC202V TRC2 and SHC002V TRC1 non-targeting controls were used. 18 h after infection,
virus was removed and fresh media added for 48–72 h. Puromycin (2 μg/mL) selection was performed
for 3 days, after which knockdown (KD) of the targeted gene was confirmed by qRT-PCR and Westerns,
and cells were treated with Doxo (250 nM) for 24 h to induce senescence. Cell viability was assessed
every 2–3 days to determine senolysis.
Cell viability assay
Cell viability was assessed using a cell counting kit-8 from Dojindo Molecular Technologies, Inc (Cat
CK04-20). Cells were incubated with culture media + 10% kit cell viability reagent, 150 μL/well in 96-well
plates for 2 h. The supernatant was transferred to a new 96-well plate and absorbance read at 450 nm
as recommended in the product datasheet.
Cell death assay
Cell death was analyzed by the Cytotoxicity LDH Assay Kit (Dojindo), according to manufacturer’s
protocol. Briefly, the supernatant was collected 7 d after removing medium containing Doxo, mixed with
the same volume of working solution, incubated for 30 min at 37 °C in a humidified 5% CO2/3% O2
atmosphere, and then absorbance at 490 nm was measured with a microplate reader.

21

SA- β-gal assay
SA-β-gal staining was performed using the Biovision kit (Prod# K320-250). Cells were plated at 5K/cm2.
SA-β-gal was assessed 24–48 h later, when cells were 60%–70% confluent. Cells were washed with PBS
and fixed for 5 min. After another PBS wash, cells were incubated overnight at 37º C in staining solution.
After washing with PBS, a minimum of 300 cells were counted. Positive (blue) cells were scored as
percentage of total cell number.
Quantification and statistical analysis
Statistical analyses were conducted using unpaired two-tailed Student’s t-test and Dunnett's multiple
comparison test using GraphPad Prism 9 software (GraphPad Software, Inc.). Statistical significance was
assigned for p-values <0.05.
RNA extractions from muscle biopsies for bulk analysis
Muscle biopsies were homogenized using a mortar and pestle with liquid nitrogen. RNA was extracted
from the powdered samples using the RNeasy Fibrous Tissue Mini Kit (Qiagen) and the QIAcube
automatic processor (Qiagen). Integrity and concentration of the RNA was assessed using the
Tapestation 4200 (Agilent Technologies), with the cut-off for acceptable integrity being an RNA integrity
number (RIN) >7. Batch-tag-seq libraries were then produced from the RNA and run on the HiSeq 4000
by the DNA Technologies Core at U.C. Davis.
Single nuclei RNA extraction and processing – human muscle biopsies
Nuclei were isolated from the biopsies using the Singulator instrument (S2 Genomics). The instrument
was primed with cold nuclei isolation and storage buffers (S2 Genomics) and the biopsy was loaded into
the cartridge and covered with the 19.7 mm grinding cap. The “Nuclei_All_Tissues” protocol was used to
isolate nuclei, after which the nuclei were centrifuged for 5 minutes and resuspended in buffer. The
22

nuclei were counted using a Countess II Automated Cell Counter (Invitrogen), centrifuged again for 5
minutes, and resuspended at the proper concentration for use in the Chromium Single Cell 5’ Library
and Gel Bead Kit v1 (10× Genomics). Samples were processed using the Chromium Single Cell A Chip Kit
and Chromium Controller (10× Genomics). Quality control was performed on the Tapestation 4200
(Agilent Technologies). Libraries were sequenced in one lane of a NovaSeqS4 by the U.C. Davis Genomics
Core.
Bulk RNAseq data analysis – human muscle biopsies
Reads were aligned to the human genome using the STAR aligner and GRCh38 as the reference genome.
Counts were computed using the featureCounts function in the subread software. Genes with a total
count of less than 10 were removed from the analysis. PCA of differentially expressed genes were
derived by the DESeq2 library in the R software, with absolute log Fold Change (logFC) >1.5 and False-
Discovery Rate (FDR) <5%. Supplementary Table 3 shows QC details for libraries prepared for each
sample. As a convention, anytime we refer to “A vs. B”, genes with a positive fold-change are
upregulated in “B”.
Single-nuclei 5’ RNAseq data analysis – human muscle biopsies
Reads were mapped to the human genome using Cell Ranger (10× Genomics), and the GRCh38 genome
reference. Cells were removed if they expressed <200 unique genes. Genes not detected in any cell
were removed from subsequent analysis. One sample was discarded due to low quality. Samples in this
study averaged less than 1% mitochondrial counts per sample, but counts were higher in the old
samples (p < 10^–16). Read count normalization, variable feature detection (nfeatures = 2000), scaling,
UMAP (ndim = 10), and differential expression were computed as described in the Seurat package [94].
Clustering was performed by the Louvain algorithm (resolution = 0.2). No batch effect was observed.
Cell types were characterized by a combination of known markers and de novo cluster markers.
23

Pathway analysis – human muscle biopsies
To derive the pathways containing differentially regulated genes, we performed a hyper-geometric test
to assess over-representation. We used the clusterProfiler R package on the database Gene-Ontology
(GO), Kyoto Encyclopedia of Genes and Genomes (KEGG). Gene set enrichment analysis was performed
[95], using log Fold-Change as ranking metric.
Spatial transcriptomics – human muscle biopsies
Paraffin embedded human skeletal muscle biopsies (4 young, 3 old) were profiled using the GeoMx
Digital Spatial Profiler (nanoString), as previously described [96]. The standard GeoMx workflow was
employed [96] 5 μm tissue sections were cut from the fresh-frozen tissues, and mounted on superfrost +
glass slides and used within a week of cutting. Sections were then fixed in 10% NBF for overnight, then
antigen retrieval followed by a proteinase K digestion to make available RNA within the tissue (1 mg/ml).
GeoMx Hu WT probes (>18,000 protein coding genes) were then incubated overnight in a humidified
chamber at 37°C to allow probes to anneal to expressed genes in the tissue. Morphology markers were
then applied to the tissue for 1 hour in a humidity chamber; we used a directly conjugated antibody
(FITC-525 nM) against desmin (for fiber detection), and for nuclei we used Syto83 (Cy3/568 nMM).
Slides containing probes and morphology markers were then loaded into the DSP, and polygonal regions
of interest within individual fibers were then drawn on each sample, with up to 4 tissue sections/slide,
for a total of eighty ROI’s across the 7 samples. ROI counts were normalized by area. Then, PCA,
differential gene expression testing was performed using a linear model on each gene, comparing ROIs
from young and old.

24

Cell culture - HSMMs
Human skeletal muscle myoblasts (HSMMs; Lonza) were maintained at 37°C, 3% O2 and 5% CO2 in
skeletal muscle growth media (SkGM2; Lonza). Differentiation was induced by replacing SkGM2 media
with DMEM supplemented with 2% horse serum. For senescence induction, cells were treated with 50
nM (differentiated cells) or 250 nM (un-differentiated cells) of doxorubicin (Doxo) for 24 hours and then
cultured for 7 days.
qRT-PCR - HSMMs
Total RNA was isolated from HSMMs 7 days after Doxorubicin treatment using the RNeasy mini kit
(Qiagen), and reverse-transcribed using the PrimeScript RT reagent kit (TAKARA Bio Inc.). Expression
levels of the genes of interest were measured by real-time quantitative PCR using a CFX-384 instrument
(BioRad). The sequences of the primer pairs are indicated as shown below. The amount of mRNA was
normalized to that of β-ACTIN mRNA.
IMR90 culture and senescence induction - timecourse
All experiments were performed in incubators at 37 °C, 3% O2 and 5% CO2. Senescence was induced by
treatment with doxo (250 nM for 24 h) (Demaria et al., 2017). Proliferating cells were harvested from
non-synchronized normally cycling cells in 10% FBS. Quiescence was induced by serum starvation (0.2%
FBS) for 3 days, followed by cell collection. Treated cells were harvested at post-induction timepoints
including 12 hours post-doxo, 36 hours, and days 2, 4, 6, 8, 10, and 12.
10x Single cell library preparation and sequencing - timecourse
Harvested cells were prepared as input for 10x Chromium Next GEM Single Cell 3’ Reagent v3.1 kits.
Sample processing and library preparation was performed according to the 10X supplied protocol. To
reduce batch effects, only cell encapsulation and cDNA conversion for each sample was performed on
25

the day of cell harvest. Further processing and library construction were then performed simultaneously
for each sample. QC for cDNA and libraries was performed using the Qubit fluorometer and Agilent
Tapestation, according to the 10X protocol. Sequencing was performed on an S4 flowcell type Illumina
NovaSeq 6000 system.
Single cell data analysis - timecourse
Raw sequencing data was processed with Cellranger count (version 3.0.2). Cellranger output was then
processed in Seurat (version 4.2). Doublets were identified with the DoubletFinder R package (version
2.0). Cells were excluded from the analysis unless they met the following criteria: percent of reads
mapping to MALAT1 < 20, percent reads NEAT1 < 20, nFeature_RNA > 500, nCount_RNA > 2000, percent
mitochondrial reads < 25, percent ribosomal reads < 25. SCTransform was used to perform
normalization, variance stabilization, and feature selection. Percent reads from mitochondrial,
ribosomal, MALAT1, and NEAT1 expression were regressed in a second linear regression. Proliferating
cells were identified by expression of MKI67, CCNA2, TOP2A, MCM2, and PCNA. Further details and R
scripts used for sample integration, cluster identification, DEG analyses, and trajectory fitting and
visualization using Slingshot and tradeSeq, are available in appendices. IPA core analysis was used to
identify differentially regulated pathways and regulatory element activity.
Strain echocardiography measurements
We placed each animal under 1.5% isoflurane anesthesia, in a supine position, and removed hair from
the thorax using depilatory cream. Mice were imaged using a Vevo 3100 high-resolution small animal
ultrasound system (FUJIFILM VisualSonics Inc.). For mice, we used a 40 MHz ultrasound transducer
(MX550D). The probes for mice were positioned to obtain a parasternal long axis and short axis views of
the left ventricle (LV), and to limit rib and lung shadows. The stage was titled to optimize the view of the
LV. B-mode data was acquired across 1000 frames for speckle-tracking analysis. We performed speckle
26

tracking analysis using the Vevo Strain software (Vevo LAB, FUJIFILM VisualSonics Inc.). Three
consecutive cardiac cycles from the long and short axis videos were chosen for analysis using the
anatomical M-mode setting to limit respiratory artifacts in LV wall motion. 20-30 points are placed along
the LV endocardial wall border in one frame within the three isolated cardiac cycles. Points were placed
to be relatively equidistant, avoiding regions where rib or lung shadowing artifacts were observed.
Speckle tracking analysis resulted in values for velocity, displacement, strain, and strain rate at each
point. Wall motion tracking was manually adjusted frame by frame as necessary. P values on multiple
comparisons of cohorts were calculated using ANOVA with Bonferroni post-test. P values of .05 or less
were considered statistically significant.
Pulse wave velocity measurements of in vivo aortic stiffness
Aortic stiffness was measured in vivo using aortic PWV. Mice were placed under light isoflurane
anesthesia (1.5%), positioned supine on a heating pad. Paws were then secured to ECG electrodes. The
two Doppler probes were positioned at the transverse aortic arch and the abdominal aorta. To calculate
aortic PWV, the distance between probes was divided by the time difference between the thoracic and
abdominal aorta impulses.

27

CHAPTER 3: A SINGLE CELL TIME COURSE OF SENESCENCE UNCOVERS DISCRETE CELL TRAJECTORIES
AND TRANSCRIPTIONAL HETEROGENEITY (Ciotlos et al, pre-print 2023, references to my contributions
to Limbad et al 2022, and Perez et al 2022)
Abstract
Senescent cells (SnCs) are typically studied as endpoints of a complex cell state transformation process,
and frequently have maladaptive effects on surrounding tissue and cells. SnCs accumulate with age, and
while they typically comprise a small percentage of cells in tissues, they have important roles in age
associated pathologies. Several obstacles remain in understanding the heterogeneous nature of
senescence, and formulating potent beneficial intervention strategies. One approach targets senescent
cells and kills them (“senolytic” approach), and is often driven by a low resolution understanding of SnC
identity, which risks both incomplete clearance and off-target effects. Cellular senescence is not a
singular binary response, but a suite of response trajectories that vary by multiple parameters, including
the inducer and initial cell state. In order to elucidate the developmental trajectories of SnCs, we
performed single cell RNA sequencing on IMR90 human lung fibroblasts senescencing across a 12 day
time period. Our analysis reveals substantial heterogeneity in gene expression within time points and
across the full time-course. We uncovered unique markers and differentially regulated pathways in cell
populations within each time point. Supervised trajectory inference of the time-course data uncovered
the root-origin and fates of distinct SnC lineages over 3 stages of senescence induction. Altogether our
data provide a novel approach to study SnC development, identifying cell states of interest, and
differentiating between SnCs and quiescent cells. This analysis will aid in identifying key targets for
therapeutic intervention in senescence.

28

Background
Cellular senescence is a complex cell state typically characterized by a permanent cell cycle arrest,
deregulated metabolism, macromolecular damage, and an altered secretory profile relative to cycling or
quiescent cells. Multiple types of stimuli can provoke cellular senescence and a cell type-specific
senescence-associated secretory phenotype (SASP) [1-4]. SnCs are implicated in many physiological
processes (both adaptive and maladaptive) and diseases of aging, including Alzheimer’s and Parkinson’s
disease, atherosclerosis and cardiovascular dysfunction, chemotherapy cardiotoxicity, and
osteoarthritis. Benign processes associated with senescence include wound healing, tissue regeneration,
embryogenesis, and tumor suppression [5-18]. Senescence is triggered by multiple distinct inducers,
including telomeric and non-telomeric DNA damage (from genotoxic drugs such as the chemotherapy
drug doxorubicin), irradiation, oncogene expression, oxidative stress, and mitochondrial dysfunction
[19]. The inducer chosen for specific cell types is a source of potential heterogeneity in the senescence
response. In this study, we induced human fetal lung fibroblasts (IMR90) to senesce using doxorubicin
(doxo), a chemotherapy drug that has multiple deleterious side effects, including cardiotoxicity and
induction of senescence. Two major pathways have been heavily investigated in the context of
senescence: the p53/p21 (p21/CDKN1A) and p16Ink4a/retinoblastoma (p16/CDKN2A) protein pathways.
Across species and cell types, there is variability in the senescence response [20]. Many inducers
activate either or both pathways [21]. SnCs also acquire a SASP, which often consists of chemokines,
cytokines, growth factors, matrix remodeling proteins and bioactive lipids [22-24].
The SASP is often a pro-inflammatory profile that leads to alterations in the tissue microenvironment
and paracrine transmission of cell senescence [2, 25-27] . Major components of several SASPs include
interleukins IL-1, 6, and 8, vascular endothelial growth factor A (VEGFA), and matrix metalloproteinases
such as MMP3. Proinflammatory cytokine components of the SASP, such as IL6, IL8, IL1α/β, TNFα, and
TNFSF7, are generally conserved across senescence contexts. However, as with senescence markers in
29

general, the composition of the SASP is fluid, and varies by induction method, cell type, and organism
[22, 24, 28-31]. Heterogeneity in cellular senescence refers to the variation in the molecular and
functional responses of individual senescent cells, despite their shared growth arrest [20]. This
variability contributes to the diverse roles that SnCs can play in physiological processes and age-related
disorders. The killing or modification of SnCs has been shown to ameliorate age related pathologies and
extend median lifespan in mice [6]. Several approaches are used to accomplish this, including senolytic
molecules, immune modulation, and “senomorphic” compounds which suppress certain aspects of the
SASP without killing the cell [32].
Current therapeutic strategies are based on a limited understanding of senescence, often targeting non-
specific molecules or pathways that are shared broadly across SnCs of different cell types. For example,
resistance to apoptosis is a shared feature of many SnCs. The drug Navitoclax (ABT-263) inhibits the
anti-apoptotic protein BCL2 and induces SnCs to undergo apoptosis [33]. However, this method also
risks eliminating beneficial or neutral cells. Transgenic mouse models of SnC clearance, such as the p16-
3MR line[12], also depend on a broad targeting strategy. In p16-3MR mice, ganciclovir kills any cell
expressing a p16-driven reporter construct. This approach also risks eliminating non-senescent cells
[12]. For example, certain macrophages express p16 as a normal response to certain immune stimuli
[34]. In addition, this model will not target SnCs lacking p16 expression, as the method relies on defining
SCs by p16 expression alone. Previous work on SnC heterogeneity (at the bulk RNA level) showed that
the transcriptional response to senescence induction changes across time, while noting that the ability
to discriminate between senescence subtypes is critical to developing specific therapy targets [20, 35].
To identify specific and time-sensitive markers of SnCs, we used scRNA-seq to survey and compare
proliferating cells, low serum-induced quiescent cells, and senescing cells across 10 time points derived
from a single culture stock. A population of senescence-evading cells was identified across the time
course, which resisted both the initial senescence induction, and remained proliferative despite being
30

cultured in a SASP-enriched environment. These cells distinguished themselves transcriptionally from
both normal proliferating cells and stable SnCs. We observed that p16 and p21 expression was not
limited to SnCs, and identified specific markers for each of the 3 main cell states that were studied.
Finally, we performed pseudotime and trajectory analyses, uncovering transcriptionally distinct lineages
of SnCs and populations engaged in different senescence-associated damage/stress responses, cytokine
signaling patterns, and metabolic activities.
Results
To explore the heterogeneity of cellular senescence both as an end state and a dynamic process, we
cultured human embryonic lung fibroblasts (IMR-90) in three cell states: proliferation, quiescence, and
senescence. To induce quiescence, we cultured cells in low serum media; for senescence induction we
treated cells with doxo for 24h and harvested them 12 hours, 1 day, 2, 4, 6, 8, 10, or 12 days later
(Figure 1A). These cell states and time points were chosen to capture transcriptional changes as cells
transform from a proliferative state to a senescent state. To ensure robustness in our experimental
design, all experiments were run with two biological replicates. Single cells were prepared by
microfluidic encapsulation and libraries were prepared according to the 10x protocol and submitted for
high throughput sequencing at the UC Davis Core Sequencing facility. Quality control included the
removal of cells with low gene and read counts, those that were highly enriched for mitochondrial and
ribosomal gene expression, as well as those that exhibited long non-coding RNAs MALAT1 and NEAT1.
QC exclusion resulted in 53,435 cells from Replicate 1 and 79,074 cells from Replicate 2 (Figure 1B), with
an average sequencing depth of 13,251 and 11,726 reads per cell, respectively. Average gene counts
were 2,972 and 2,966 for Replicate 1 and Replicate 2, respectively. Replicates were visualized together
to assess whether significant differences existed in cell embedding and clustering (Figure1C), which
could indicate either biological or technical variation. Although the replicates do not overlap perfectly,
we did not observe strong batch effects necessitating correction, a process that can reduce biological
31

variation while attempting to correct for technical variation. Cells from each time point do not strongly
segregate by replicate, indicating batch effects did not significantly affect library quality. We
determined that the percent of reads assigned to cells declined at D6, D8, and D10 post-doxo before
slightly recovering at D12 (Figure 1C). This decrease is due to an increase in ambient (non-cell-assigned)
RNA, suggesting that a fraction of senescing cells become fragile 6-10 days following doxo treatment,
and may lyse due to the shear stress that cells experience during microfluidic processing. These cells
were excluded from our analysis. Due to the high degree of homogeneity between replicates, we
combined them for increased power to characterize heterogeneity in a sizable cell population.

32

Figure 1. Experimental overview and quality control of replicate libraries. (A) Proliferating IMR90 cells
treated with 250 nM doxo for 24h are shown at either 12 hours (referred to as “H12”) ,1 (“D1”), 2
(“D2”), 4 (“D4”), 6 (“D6”), 8 (“D8”), 10 (“D10”), or 12 (“D12”) days after doxo treatment. (B) Time point
replicates visualized by UMAP, cells colored by time point. (C) Overlaid replicates visualized by UMAP
and percent of reads assigned to validated cells from each time point replicate library.

33

A population of proliferating cells evades acute senescence induction
Despite treatment with doxo, we observed that some cells were able to “evade” senescence conversion.
We identified 8,449 cells across the 12-day time course which remained proliferative (expressing
MKI67+/CCNA2+/TOP2A+/MCM2+/PCNA+); we visualized these cells at each time point using UMAP
(Figure 2). The presence of proliferating cells suggests that not all cells within a culture are induced to
senescence following 24h of doxo treatment. To further characterize these cells, we identified
differentially expressed genes and upstream regulators, between sequential time points (Table 1). 12
hours after doxo treatment of proliferating cells, the growth factor binding protein IGFBP5 and matrix
metallopeptidase MMP1 were the top up- and down-regulated genes. We identified ILF3, a
transcription factor known to also bind mRNAs of SASP proteins, and TP53, a central activator of the
senescence response, as activated upstream regulators. Conversely, the activity of the transcription
factor E2F3 and methyltransferase DNMT3B was inhibited.
From 12 hours (H12) to 1 day (D1) (H12-D1), we observed increases in growth factor and the SASP
protein GDF15, and decreases in another SASP factor, SFRP1. TP53 gene expression activity continued
to increase from H12 to D1, and histone demethylase KDM5B gene expression changed, consistent with
its epigenetic function. Transcription factors FOXM1, which plays important roles in aging and
senescence, and HSF1 (heat shock regulator whose suppression activates p38–NF-κB-SASP) were
inhibited. Comparing those cells that escape senescent transformation from D2 to D1, we determined
that IGFBP5 is upregulated. In contrast to the previous (D1 vs. H12) comparison, GDF15 was strongly
downregulated by D2. The reduction in GDF15 may indicate an immediate response that is attenuated
with time, or that such cells are able to suppress the need to release GDF15. mRNA levels of JARID2 and
BACH1 also increased. The latter is a known inhibitor of p53 and oxidative stress-induced senescence.
Indeed, transcriptomic regulation by TP53 was inhibited at D2, reversing the trend of the previous time
points. Histone deacetylase HDAC3 function, recently shown to repress p65 NF-κB-mediated induction
34

of the SASP, was also inhibited. Between D2 and D4, expression of the known SASP factor CCL2
increased, while expression of the secreted factor and chaperone CLU declined. These findings
illustrates the unique transcriptional activity that separates the senescence resistant cells from more
robust SnCs at this time point.
Changes in gene expression between the next pair of time points (from D4 to D6) describe a
transcriptomic profile more akin to senescent cells: actin-binding protein FLNA is upregulated and serine
protease inhibitor TFPI2 is downregulated. TFPI2 is an inhibitor that is highly expressed in non-invasive
cells; TFPI2 increases at apoptosis and senescence. FLNA was shown to induce senescence in the
context of hypoxia, by complexing with DRP1 and leading to mitochondrial fragmentation. At this point,
gene expression of the senescence-induced GATA6 transcription factor increased, while androgen
receptor (AR) and HOXD10 activity declined. COL6A1 is the top upregulated gene in senescence
“evaders” (SnEs) from D6 to D8, confirming previous data indicating that its knockdown suppresses cell
proliferation and migration. TXN gene expression declined; its suppression has previously been shown
to induce senescence in normal human fibroblasts. Expression of the transcription factor GLIS1
increases, suggesting a link between senescence evasion and its roles as a facilitator of pluripotency
induction in somatic reprogramming and fibrosis of senescent fibroblasts. Between D8 and D10,
LRRC75A and CENPF declined. CENPF is required for chromosomal segregation and kinetochore
function. CDKN1A/p21 activity increased at this time point, potentially indicating a different role for this
gene in senescence-evading cells and SnCs. FOSB is the top upregulated gene from D10 to D12,
consistent with previous findings linking its expression to increased proliferation. IGFBP7 is secreted by
SnCs and its upregulation inhibits proliferation; in our study this gene is downregulated. Transcription
factors SATB1 and ETV5 are both increased in gene expression at this time point, which supports
previous findings indicating their loss suppresses proliferation.

35

Figure 2: Visualization of proliferating cells. UMAP of 8,449 proliferating cells, drawn from both the
untreated proliferating (“Pro”) time point and senescence resistant (H12-D12) time points, identified
across the full time course.
Table 1. Top differentially expressed genes and transcriptional regulators between senescence
resistant cells from sequential time points. Differentially expressed genes (DEGs) were identified with
Seurat’s FindMarkers function, and regulators were identified with Ingenuity IPA. Empty boxes indicate a
result flagged for bias by Ingenuity’s IPA analysis and are not reported.
Time point
Comparison
Top
Upregulated
Gene
Log
FC
Top
Downregulated
Gene
Log
FC
Top
Activated
Regulators
Top Inhibited
Regulators

H12 vs.
Proliferating
IGFBP5 2.17 MMP1 -1.71 ILF3, TP53
E2F3,
DNMT3B
D1 vs. H12 GDF15 0.99 SFRP1 -0.8 TP53, KDM5B
FOXM1,
HSF1
D2 vs. D1 IGFBP5 1.08 GDF15 -2.91
JARID2,
BACH1
TP53, HDAC3
D4 vs D2 CCL2 1.54 CLU -1.51 E2F3, SF1 TP53, NR3C2
D6 vs D4 FLNA 3.12 TFPI2 -2.94 HLX, GATA6 AR, HOXD10
D8 vs D6 COL6A1 3.94 TXN -3.03 GLIS1
D10 vs D8 LRRC75A 1.73 CENPF -1.16 CDKN1A E2F3
D12 vs D10 FOSB 1.01 IGFBP7 -1.2 SATB1, ETV5

36

Table 2. Pathway analysis of activated and inhibited pathways in senescence-resistant cells compared
to senescent cells. Consistently modulated pathways of both proliferating and senescing cells at all time
points. “Activated” denotes higher pathway activity in proliferating cells while “Inhibited” denotes
lower pathway activity in proliferating cells. Analysis was performed by Ingenuity’s IPA analysis.
Pathway Category Pathway Status
Cell Cycle and
DNA Replication
Cell Cycle Control of Chromosomal Replication Activated
Cyclins and Cell Cycle Regulation Activated
ID1 Signaling Pathway Activated
RAN Signaling Activated
Cell Cycle: G2/M DNA Damage Checkpoint Regulation Inhibited
DNA Repair
Nucleotide Excision Repair Pathway Activated
Base Excision Repair Pathway Activated
Immune Function
Immunogenic Cell Death Signaling Pathway Activated
Fcγ Receptor-mediated Phagocytosis in Macrophages Activated
Granzyme A Signaling Inhibited
Senescence
HMGB1 Signaling Activated
VEGF Signaling Activated

Differentiating between proliferating, senescent, and quiescent cells
We compared the 3 main outcomes from doxo treatment (proliferating, quiescent and senescent cells;
31,444 cells total) to define the transcriptional differences between cell states. Here we report
transcriptional markers for each cell state, while also highlighting heterogeneity within the respective
cell states that may otherwise confound identification. Some cells from the D12 time point retain
37

proliferative capacity (Figure 3A, top), as previously discussed. These cells fall into cluster 12 (Figure 3A,
bottom), and are not included in the results here. A small number of quiescent cells clustered (cluster
14) in close spatial proximity to the main body of untreated proliferating cells (clusters 10 and 11). In the
(untreated) proliferating cell time point, MKI67 was not ubiquitously expressed (found in 45% of cells),
nor were other proliferation markers such as TOP2A (43%), CCNA2 (36%), and MCM2 (25%) (Figure 3B).
Proliferating cell nuclear antigen (PCNA) was expressed in 62% of these cells, as opposed to 55% of QnCs
and 16% of SnCs. In the (untreated) proliferating cell time point, MKI67 was not ubiquitously expressed
(found in 45% of cells), nor were other proliferation markers such as TOP2A (43%), CCNA2 (36%), and
MCM2 (25%) (Figure 3B). Proliferating cell nuclear antigen (PCNA) was expressed in 62% of these cells,
as opposed to 55% of QnCs and 16% of SnCs.
38

Figure 3A)

39

Figure 3B)

40

Figure 3C)

Figure 3D)

41

Figure 3: Cyclin-dependent kinase inhibitors are variably expressed in proliferating, quiescent, and
senescent cells, which form heterogeneous populations identified by non-traditional markers. (A)
Visualization of populations identified in comparisons of proliferating, quiescent, and SnCs; by time
point (top) and cluster (bottom). (B) Scaled expression values of MKI67, a number of CDKN family
genes, and three genes (SH3PXD2A, DDX39A, TMEM243) determined to be specific markers of cell state.
(C) Expression heatmaps of genes identified as being specifically expressed in respective cell states. (D)
Density plot of differentially expressed genes upregulated in senescent (top row), proliferating (middle),
and QnCs (bottom).

The upregulation of cyclin-dependent kinase inhibitors CDKN1A (p21) and CDKN2A (p16) transcripts are
frequently reported to be markers of cell senescence. We compared the expression of these two genes,
and their inhibitors CDKN1B (p27/Kip1), CDKN2B (p15/INK4b), and CDKN1C (P57/Kip2), in our cell states
of interest. In SnCs, the most highly upregulated (by log fold change) of these two genes was CDKN1A,
and was expressed in 92% of these cells. CDKN1A was also expressed in 71% of QnCs and 64% of
proliferating cells, albeit to lower degrees. To our surprise, CDKN2A was most highly upregulated in
QnCs (81%). Interestingly, 43% of proliferating cells also expressed CDKN2A, while only 28% of SnCs had
significant expression. CDKN1B was downregulated in proliferating cells compared to quiescent and
senescent cells, with the latter two cell states showing similar levels of expression. In small numbers of
quiescent and proliferating cells, CDKN2B was slightly up- and downregulated respectively. CDKN1C was
upregulated in both QnCs (67% of cells expressing) and SnCs (39%) compared to proliferating cells
(13%). However, we did not measure protein levels for these genes at the level of the single cell.
The most specific gene expression marker in SnCs was the SH3 and PX domain-containing protein
SH3PXD2A (expressed in 66% of SnCs and only 33% of the other cells. The RNA helicase DDX39A was
both upregulated and a specific marker for proliferating cells (81%, compared to 20% of proliferating
and QnCs). The transmembrane protein TMEM243 was highly specific to QnCs, with 80% of them
expressing the gene compared to only 17% of the remaining cells; other markers of considerable
specificity for each cell state are plotted in Figure 3C. In this 3 way comparison, SnCs upregulated type
42

VI collagen alpha chains COL6A1 and COL6A2, the ECM-associated peroxidase PXDN, secreted ECM
component CCDC80, and myosin light chain kinase MYLK (Figure 3D). Upregulated genes in proliferating
cells included matrix metallopeptidase MMP1, ubiquitin ligase UBE2S, centromere and microtubule
protein CENPF, microtubule component TUBB4B, and pituitary tumor transforming PTTG1. Finally, QnCs
upregulated thioredoxin TXN, acetyltransferase ACAT2, the cell cycle regulator RGCC, the peroxisomal
enzyme IDI1, and monooxygenase MSMO1.
Single cell transcriptomic time-course of senescence
To assess the full extent of heterogeneity within time points and across the time course of
transformation to SnCs, we performed a merged analysis of proliferating (“Pro”) and post-doxo
treatment (H12, D1, D2,...D12) samples. We adjusted clustering parameters in this analysis of 126,163
cells such that the population of senescence evaders ( ~6.5% of all cells) was restricted to one cluster.
Figure 4A groups cells by time point, and the UMAP projection infers an unsynchronized progression
through senescence. Apart from co-clustering of D1 cells, we identified a complex pattern of co-
localization and overlap in cells from different time points. We found 46 (Figure 4B) distinct groups of
cells, including the presence of proliferating cells from the post-doxo treatment (Cluster 33).
Clusters were composed of non-homogenous cell populations from multiple time points (Figure 4C).
Clusters 0, 17, 26, 33, 37, and 44 were dominated by cells from D12, suggesting these are clusters of
stable SnCs. However, many cells from earlier time points also co-cluster with these D12 cells. This
suggests that with senescence induction, some cells quickly (as soon as 12 hours post doxo treatment)
reorient their transcriptional profiles from proliferative potential to senescence. Conversely, some cells
appear to lag in their transformation to SnCs. This lag is illustrated by the presence of D12 cells in
clusters dominated by cells from time points D4/D6/D8 (Clusters 3, 4, 12, 13, 15, 19, 25, 32, 34).
43

Our study shows that the development of senescence is a dynamic process, through which cells do not
progress in lock-step. This phenomenon emphasizes the value of single cell profiling, while raising
concerns about bulk preparation workflows of SnCs at either protein or gene expression levels.
Evaluating the progression of senescence across the experimental timeline provided insight into the
heterogeneity among cell populations at each time point. However, this approach did not facilitate
identifying and tracking cell trajectories across sequential time points. To capture greater
transcriptional specificity in smaller timescales, we carried out integrated analyses of subsets of the time
course.

44

Figure 4A) top, 4B) bottom

45

Figure 4C)
Figure 4: Visualization of senescence time-course. A) Clustering of 126,163 cells across 9 time points,
plotted with the UMAP dimensional reduction technique. Pro = proliferating cells, H12 = 12 hours after
doxo treatment, D1 = 1 day after doxo treatment, D2 = 2 days after doxo treatment, and so on up to
D12. B) Louvain clustering of merged data identifies 46 distinct clusters across the full time-course. C)
Cluster composition plotted by time point.

46

Validating markers of proliferation and senescence across transcriptomic time course
We examined the expression patterns of markers of proliferating cells and commonly reported markers
of senescence. MKI67, TOP2A, and CCNA2 were used as markers of cycling cells. These three genes
were expressed in untreated proliferating cells and were downregulated as quickly as 12 hours after
senescence induction. Furthermore, they were sporadically expressed in subpopulations of cells from
post-doxo treatment time points (Figure 5A).

47

Figure 5: Normalized expression of markers of proliferation and senescence. Darker purple shading
indicates increased expression. A) Genes MKI67 (left), TOP2A (middle), and CCNA2 (right) were
predominantly expressed in clusters of proliferating cells. B) Genes HMGB2 (left), LMNB1 (middle),
EZH2 (right), all commonly negative senescence markers, were also expressed in clusters of proliferating
cells. C) CDKN1A (left), CDKN2A (middle), and GDF15 (right), all commonly positive markers of
senescence, were upregulated to different degrees following doxo treatment.

We also examined expression patterns of negative markers for proliferation HMGB2, LMNB1, and EZH2
(Figure 5B), genes whose proteins are reported downregulated as cells enter senescence. These
markers are not often detected via common bulk RNA analyses, and changes in their expression at the
RNA level have been deliberated in the senescence field. Here, we demonstrate that transcriptional
changes of these markers can be captured using a time course-based single cell approach. We
investigated how the expression of the commonly used positive senescence markers CDKN1A, CDKN2A,
and GDF15 change across the time course employed in this study. CDKN1A/p21 expression was widely
distributed - its low, sporadic presence in proliferating cells may be due to its role in transitory cell cycle
arrest at the G1 checkpoint. As early as 12 hours after senescence induction we observed a transient
increase in its expression, which then increased from D1 to D4. By D12, its expression was lower, which
may be attributed to cells reaching a stable arrest, no longer necessitating p21 overproduction. In
contrast to CDKN1A, CDKN2A (p16 and p14arf) expression was much weaker throughout the time
course, despite an increase at D2 and D4 when compared to proliferating cells. The cytokine GDF15 is an
important senescence-associated secretory phenotype (SASP) molecule, which plays a role in apoptosis
and inflammation signaling. We observed an early upregulation of GDF15 from 12-24 hours (H12 and
D1), before subsiding starting at D2. This finding contrasts with the unchanged CDKN2A levels from H12
to D1. Altogether, these findings serve as confirmation of senescence induction using several commonly
used senescence markers.

48

Identification of cell trajectories across early, middle, and late phases of senescence
We next investigated whether distinct lineages of senescent cells could be identified. We separated the
total time course into three distinct phases: early (proliferating to D2), middle (D2 to D8), and late-stage
senescence (D8 to D12), allowing us to explore each with greater specificity. These temporal stages
were chosen in conjunction with normalization of expression levels of CDKN1A (Figure S3), as well as
changes in ambient RNA fraction (Figure 1). Following processing by the Seurat pipeline , the Slingshot
and tradeSeq R packages were utilized to perform pseudotime analysis and visualize the transcriptional
shifts along identified trajectories. These analysis packages were chosen due to their ability to identify
multiple trajectories in heterogeneous single cell data. Of particular interest were clusters in which a
given lineage diverged from other lineages, and clusters in which the given lineage ended (putatively the
most temporally advanced cells of a lineage). The top up- and down-regulated genes (by log2 fold
change) were identified for the divergence and end clusters, as well as the most specific up and down-
regulated genes (by percent cells expressing a given gene). Early stage senescence is composed of H12,
D1, and D2 cells. All early senescence trajectories commenced in cluster 18, largely composed of
proliferating cells, before splitting off into one of six lineages (Figure 6a).
49

Figure 6a
50

Figure 6b
51

Figure 6c
52

Figure 6d
Figure 6: Senescing cells form distinct transcriptional lineages as they progress from proliferation
through later-stage senescence. Lineage trajectories are identified by number and end cluster. (A) From
53

untreated proliferating cells to the D2 time point, early stage senescence trajectories begin at cluster 18
and mature down to one of six lineages, ending in clusters 19, 16, 21, 15, 13, and 9. (B) Mid-stage
senescence (time points D2 to D8) diverge into seven lineages which start at cluster 17 and end in
clusters 18, 6, 9, 10, 3, 16, and 19. (C) Late-stage senescence (D8 - D12) lineages originated at clusters 18
and 10. Four lineages start from both cluster 18 (lineages ending in clusters 8, 15, 22, 5) and (D) cluster
10 (lineages ending in clusters 14, 21, 19, and 20).
For each lineage, we identified genes with positive expression and specificity to their cluster relative to
others (Table 3), as well as negative changes (Table 4). Positive changes within early-stage senescence
were genes associated with responses to oxidative stress (TP53I3, lineage 1) and inflammatory
responses (CXCL8, lineage 3).
Table 3. Upregulated (UR) lineage markers of early-stage senescence (Pro-D2). Top genes were
determined by the highest Log2 fold change. Specificity was determined by change in the percent of
cells expressing a gene of interest. For both divergence and end clusters, comparisons were performed
by comparing a given cluster against all other clusters of that type.
Lineage Lineage
Divergence
Cluster
Top UR
log2FC gene
Top UR gene
by specificity
Lineage
End
Cluster
Top UR
log2FC
gene
Top UR gene
by specificity
Lineage
1
0 TP53I3 ARHGAP22 13 ACAT2 TUBB2A
Lineage
2
3 C12orf75 IGFBP3 15 GDF15 PSAT1
Lineage
3
12 CXCL8 AREG 21 CXCL8 FOS
Lineage
4
1 IGFBP3 IGFBP3 16 BEX3 TP53TG1
Lineage
5
13 CTGF CHIC2 19 CD63 EMILIN1
Lineage
6
17 MYL9 ACTG2 9 NABP1 CPNE7
54

Table 4. Downregulated (DR) lineage markers of early-stage senescence (Pro-D2). Top Log2FC genes
were determined by highest Log2 fold change, while specificity was determined by change in the
percent of cells expressing a given gene.
Lineage Lineage
Divergence
Cluster
Top DR
log2FC gene
Top DR gene
by specificity
Lineage
End
Cluster
Top DR
log2FC
gene
Top DR gene
by specificity
Lineage
1
0 CXCL8 AREG 13 GDF15 POLR2J3
Lineage
2
3 CLU LRPAP1 15 IGFBP5 IGFBP5
Lineage
3
12 TP53I3 ARHGAP22 21 NABP1 CDC42EP3
Lineage
4
1 COL3A1 SOX4 16 GDF15 DDX5
Lineage
5
13 MMP1 MMP1 19 TXN BLOC1S2
Lineage
6
17 SERPINE2 EMC7 9 GDF15 EMC7

Lineage 1 diverged from the rest at cluster 0, which upregulated the p53 inducible oxidoreductase
TP53I3 and the GTPase activator ARHGAP22 (Table 3). This lineage ended at cluster 13, marked by
expression of acetyltransferase ACAT2 and tubulin TUBB2A. Expression of C12orf75 and IGFBP3 in
cluster 3 marked the divergence of lineage 2, which ended in cluster 15, which expresses growth factor
GDF15 and aminotransferase PSAT1. The SASP factor CXCL8 was also upregulated across lineage 3, from
its divergence point (cluster 12) until its end cluster 21, in which the AP-1 transcription factor
component FOS is expressed. The divergence of lineage 4 occurred in cluster 1, where IGFBP3 was both
the most highly upregulated gene and most specific marker for this cluster. This lineage ended in cluster
55

16, which expressed BEX3, a protein involved in p75NTR-mediated apoptosis. The most specific marker
for this cluster was TP53TG1, a lncRNA with both pro- and anti-oncogenic properties. Lineage 5
expressed CTGF and CHIC2 when diverging, before ending in cluster 19, which expressed the cell surface
protein CD63 and extracellular matrix glycoprotein EMILIN1. The final early-stage lineage was number
6, which began its trajectory by expressing myosin chain MYL9 and smooth muscle actin expressing
ACTG2. This lineage ended in cluster 9, concomitant with upregulation. of NABP1 and CPNE7.
Interestingly, SASP genes including AREG, CLU, CXCL8, MMP1, and GDF15 were downregulated in
different early-stage lineages (Table 4). Our data show that early in senescence, cells separate into
lineages with distinct transcriptional profiles and specific expression of both SASP and non-SASP genes.
Differential regulation of pathways provided insight into functional specializations of some SnC lineages.
We determined that in lineage 1, CXCR4, IL-8, actin cytoskeleton (RHOA, Rho GTPase, and RAC), and
ephrin receptor signaling were all activated. The coronavirus replication pathway was also uniquely
activated in lineage 1, while the virus’ pathogenesis pathway was expressed in lineages 5 and 6 only.
Lineage 2 was marked by increased endoplasmic reticulum stress signaling, p38/MAPK activity, and
inflammatory signaling mediated by IL-10. The unfolded protein response (UPR) activity was highest in
this lineage, as was necroptosis signaling. The population of developing SnCs in lineage 3 was
characterized by IL-17 and hypercytokinemia signaling. Lineage 4 modulates a number of pathways
shared with other lineages, such as EIF2 signaling (also activated in lineages 1-3, inhibited in lineages 5
and 6), oxidative phosphorylation (shared with lineages 1 and 2, inhibited in others) and PPARα/RXRα
activation (shared only with lineage 6). However, no strongly specific or uniquely activated pathways
could be determined. Pathways activated in lineage 5 include osteoarthritis, collagen receptor GP6, and
pulmonary fibrosis signaling. The mitochondrial dysfunction pathway was activated in lineages 3, 5, and
6. The latter two lineages also displayed increased lysosomal activity (CLEAR pathway). Wound healing,
56

one of the beneficial roles of SnCs, was activated in lineages 1, 3, and 5. Finally, in lineage 6, GADD45
(Growth Arrest and DNA Damage-inducible 45 family genes) signaling was uniquely upregulated.
Middle stage senescence was composed of D2, D4, D6, and D8 cells, which were clustered and visualized
separately from other time points or stages of senescence. The UMAP projection shows a gradual
progression of cells from Day 2 to Day 8 after senescence induction. Cells started their trajectories in
cluster 17 (Figure 6b), largely made up of D2 cells, before proceeding into one of seven lineages. As with
early senescence lineages, we observed a large degree of heterogeneity in this sequence of time points.
While the end clusters of lineages 1 (cluster 18), 2 (cluster 6), 3 (cluster 9), and 5 (cluster 3) are
composed mainly of cells from D8, non-terminal clusters had a mixed populations of cells from all time
points. For each lineage, we identified genes with positive expression changes and positive specificity to
their cluster relative to the others (Table 5), as well as negative changes (Table 6).

57

Table 5. Upregulated (UR) lineage markers of mid-stage senescence (D2-8). Top Log2FC genes were
determined by highest Log2 fold change, while specificity was determined by change in the percent of
cells expressing a given gene.
Lineage Lineage
Divergence
Cluster
Top UR
log2FC
gene
Top UR gene
by specificity
Lineage
End
Cluster
Top UR
log2FC
gene
Top UR gene
by specificity
Lineage 1 18 COL6A1 ABCA1 18 N4BP2L2 SNRNP70
Lineage 2 15 MT2A PGAM1 6 FTL SELENOW
Lineage 3 9 CALD1 USP53 9 VMP1 LENG8
Lineage 4 13 ACTA2 ACAT2 10 FDPS ACAT2
Lineage 5 1 GDF15 LRPAP1 3 PTMA MDK
Lineage 6 16 GREM1 GREM1 16 IGFBP3 BDNF
Lineage 7 19 AREG BTG1 19 AREG ZFP36L1

58

Table 6. Downregulated (DR) lineage markers of mid-stage senescence (D2-D8). Top Log2FC genes
were determined by highest Log2 fold change, while specificity was determined by change in the
percent of cells expressing a given gene.
Lineage Lineage
Divergence
Cluster
Top DR
log2FC
gene
Top DR gene
by specificity
Lineage
End
Cluster
Top DR
log2FC
gene
Top DR gene
by specificity
Lineage 1 18 CALD1 USP53 18 TXN DBI
Lineage 2 15 DST N4BP2L2 6 COL1A1 CYR61
Lineage 3 9 COL6A1 ABCA1 9 TXN SUMO2
Lineage 4 13 GREM1 GREM1 10 VMP1 SELENOW
Lineage 5 1 C12orf75 CDC42EP3 3 CCDC80 CCDC80
Lineage 6 16 ACTA2 ACAT2 16 VMP1 TKT
Lineage 7 19 IGFBP3 MYL12A 19 PRDX1 MYL12A

Mid stage lineage 1 cells expressed the collagen gene COL6A1, the ABCA1 cholesterol phospholipid
transporter, transcription repressor complex component N4BP2L2, and small nuclear ribonucleoprotein
SNRNP70. Metallothionein MT2A, glycolytic mutase PGAM1, the iron storage ferritin chain FTL, and
glutathione dependent antioxidant SELENOW were upregulated over lineage 2. Actin-binding protein
CALD1, thioredoxin TXN, ubiquitin-like modifier SUMO2, and ubiquitin peptidase USP53 were
upregulated by cells in cluster 18 / lineage 3. Lineage 4 was defined by upregulation of smooth muscle
actin ACTA2, acetyltransferase ACAT2, and farnesyl diphosphate synthase FDPS. GDF15 was again
upregulated in lineage 5, along with LDL receptor chaperone LRPAP1, antiapoptotic protein PTMA, and
secreted growth factor MDK. Bone morphogenic protein antagonist GREM1 and nerve growth factor
59

BDNF were the strongest markers of lineage 6. Finally, AREG, anti-proliferative factor BTG1, and zinc
finger-like protein ZFP36L1 were upregulated in lineage 7.
Lineage 1 and 3 shared many activated pathways, most significantly the coronavirus pathogenesis
pathway and mitochondrial dysfunction, as well as CLEAR (Coordinated Lysosomal Expression and
Regulation). IL-12 signaling was activated in lineages 1 and 6. Oxidative phosphorylation and neutrophil
extracellular trap signaling were upregulated in lineages 2, 5, and 6. Nucleotide excision repair was
uniquely activated in lineage 2. Lineage 4 was unique in its activation of numerous cholesterol synthesis
pathways, as well as geranylgeranyl diphosphate biosynthesis. Lineages 4 and 5 shared many activated
pathways, including estrogen receptor and immunogenic cell death signaling. Despite differences at the
gene expression level, the number of shared pathways between the two lineages indicate that some
cells may take different paths to the same cell fate. Actin cytoskeleton, RHOA, and integrin signaling
were all activated in lineage 6, indicating a cell lineage undergoing cytoskeletal reorganization and cell
surface remodeling. Il-10 signaling was upregulated in lineages 6 and 2, and lastly IL-6 signaling was
uniquely activated in lineage 7. IL-8 was activated in lineages 4-7 but inhibited in lineages 1-3. From D2
to D8, SnCs evolve into lineages characterized by signaling by specific SASP cytokines, metabolic activity,
and cytoskeletal remodeling. Wound healing signaling was activated in lineages 4-7, continuing a trend
identified in the early-stage lineages.
Late-stage senescence was defined by the presence of D8, D10, and D12 cells. In contrast to early and
middle stages, these cells initiated their lineage start in one of two distinct clusters highly enriched for
D8. The first subset of cells (a) began their trajectories in cluster 18 (Figure 6c) before splitting off into
one of four lineages (lineages 1-4a ending in clusters 8, 15, 22, and 5 respectively). The second subset
also evolved into 4 lineages after starting in cluster 10 (Figure 6d); these were named lineages 1-4b, and
ended in clusters 14, 21, 19, and 20. As with the analyses of previous stages, D8-D12 cells were
embedded and clustered separately from other time points. Lineage 1a upregulated neuronal pentraxin
60

NPTX1 (Table 7), retinol binder RBP1, COL3A1, and procollagen endopeptidase PCOLCE. Follistatin (FST),
squalene epoxidase (SQLE), and again ACAT2 were markers for lineage 2a. Transcription factors KLF6
and SOX4 were upregulated in lineage 3a, alongside secreted Wnt signaling modulator SFRP1 and
histone HIST1H4C. Lineage 4 upregulated myofibroblast marker TPM2, antiapoptotic lncRNA
MTRNR2l12, and IGFBP7. In the second subset of late-stage senescence trajectories, lineage 1b
expressed the metalloproteinase and IGFBP-family cleaving PAPPA when diverging from other lineages.
IGFBP3 and carboxypeptidase CPA4 were upregulated at the end of this lineage. Collagen chain COL6A3
and secreted growth factor PTN were upregulated in lineage 2b, across the entire trajectory from
divergence until the end cluster (cluster 21). Lineage 3b upregulated phosphodiesterase PDE5A and the
lncRNA NORAD, which is produced in response to DNA damage, when diverging at cluster 19. This
lineage also expressed YBX1, a secreted cold shock domain protein that functions as an extracellular
mitogen, and SET, a nuclear protein that inhibits apoptosis and histone acetylation. The final lineage,
4b, was characterized by strong and specific expression of GDF15, unique among the late state lineages.
Pathway analysis revealed significant differences between late-stage lineages from groups (a) and (b).
Specifically activated in lineage 1a (group a) were cytokine storm, collagen receptor GP6, and wound
healing signaling. IL-4 signaling was also unique to lineage 1, whereas IL-12 activation was observed in
both lineage 1a and 3a. Cholesterol signaling was strongly activated in lineage 2a, continuing the
expression pattern from middle stage lineage 4. Lineage 2a also upregulated FC gamma-receptor
mediated phagocytosis and the NRF2-mediated oxidative stress response. Ferroptosis signaling and
CXCR4 activity were increased in lineage 3a, while lineage 4a was characterized by mitochondrial
dysfunction, HER-2 signaling, and PAK (p21-activated kinase) signaling. Oxidative phosphorylation
activity was high in lineages 1-3a and inhibited in 4a. PKR induction of interferons was activated in
lineages 1 and 3 specifically.
61

The second group (b) of late-stage lineages were analyzed separately. We found that lineage 1b
specifically activated the unfolded protein response (UPR), the WNT/Ca pathway, mitochondrial
dysfunction, and HIPPO and lysosomal signaling. Lineages 1b and 2b upregulated IL-4 signaling and the
chaperone mediated autophagy pathway. Interestingly, rheumatoid arthritis signaling was activated
specifically in lineage 2b. In lineage 3b, IL-8 and CXCR4 signaling increased, alongside integrin and
PI3K/AKT signaling. Finally, lineage 4b shared upregulated EIF2 and erythropoietin signaling with lineage
3b, but uniquely activated NRF2-mediated oxidative stress.
Table 7. Upregulated (UR) lineage markers of late-stage senescence (D8-D12). Top Log2FC genes were
determined by highest Log2 fold change, while specificity was determined by change in the percent of
cells expressing a given gene.
Lineage Lineage
Divergence
Cluster
Top UR
log2FC
gene
Top UR gene by
specificity
Lineage
End
Cluster
Top UR
log2FC
gene
Top UR gene
by specificity
Lineage 1a 17 NPTX1 RBP1 8 COL3A1 PCOLCE
Lineage 2a 15 FST SQLE 15 FST ACAT2
Lineage 3a 9 SFRP1 KLF6 22 SOX4 HIST1H4C
Lineage 4a 2 TPM2 MTRNR2L12 5 IGFBP7 IGFBP7
Lineage 1b 14 IGFBP3 PAPPA 14 IGFBP3 CPA4
Lineage 2b 21 COL6A3 PTN 21 COL6A3 PTN
Lineage 3b 19 PDE5A NORAD 19 YBX1 SET
Lineage 4b 20 GDF15 GDF15 20 GDF15 GDF15

62

Table 8. Downregulated (DR) lineage markers of late-stage senescence (D8-D12). Top Log2FC genes
were determined by highest Log2 fold change, while specificity was determined by change in the
percent of cells expressing a given gene.
Lineage Lineage
Divergence
Cluster
Top DR
log2FC
gene
Top DR gene
by specificity
Lineage
End
Cluster
Top DR
log2FC
gene
Top DR gene
by specificity
Lineage 1a 17 FST SQLE 8 CCND1 GAS6
Lineage 2a 15 NPTX1 RBP1 15 TIMP1 IGFBP7
Lineage 3a 9 TIMP1 PRDX1 22 PRDX1 ANXA1
Lineage 4a 2 FOSB FOSB 5 FOSB FOSB
Lineage 1b 14 COL6A3 PTN 14 PTMA PTMA
Lineage 2b 21 IGFBP3 PAPPA 21 IGFBP3 PAPPA
Lineage 3b 19 TUBA1B TUBB4B 19 TIMP1 LRPAP1
Lineage 4b 20 CALD1 IGFBP3 20 SPARC CRIM1

63

Figure 7: Differentially regulated pathways in senescence lineages. Pathway activation z-scores were
calculated in IPA. Pathway activation/inhibition are denoted by positive/negative z-scores. (A) Early
stage senescence is characterized by lineages that differentially express stress response, inflammatory,
and DNA damage related pathways. (B) Pathway activity in middle stage senescence lineages. (C) - (D)
Pathway activity in two branches of late stage senescence lineages.

64

Discussion
Our results reinforce previous findings about the heterogeneity of SnCs in culture, while expanding our
understanding at the single cell level, both within single time points and over a time course of
senescence development. Far from being a monolithic cell state, SnCs undergo dynamic transcriptional
changes and diverge into discrete lineages across the studied time course. Each trajectory represents a
separate senescent outcome following the same induction event by doxo treatment. We identified co-
clustering cells from different time points, indicating that induced cells transition to these distinct
senescent profiles at different speeds. Within these trajectories, differentially expressed genes and
regulated pathways illustrate a wide range of variability in SnC transcriptional activity. Altered pathways,
including chemokine signaling, mitochondrial dysfunction, lysosomal activity, and wound healing, have
previously been associated with SnCs [21, 36-40]. We show that in senescent IMR-90 cells, these
responses are not manifested across all induced cells, but characterize heterogeneous populations of
SnCs, possibly derived from differing induction mechanisms.
Our findings at the single cell RNA level expand on previous work describing the cellular outcomes of
doxo induced senescence [41-44]. While ionizing radiation induces senescence predominantly via
double stranded DNA breaks [45-47], doxo acts through several mechanisms: DNA intercalation and
crosslinking, DNA damage by topoisomerase II inhibition, mitochondrial dysfunction, and free radical
production [48-51]. In contrast to bulk RNAseq, single cell profiling makes it possible to track the
development of SnC populations originating from different cellular damages branching out to multiple
cellular response. We describe these lineages in terms of transcriptional state at root, divergence (from
other lineages), and end clusters, retracing the transcriptional steps SnCs take as they enter growth
arrest. These profiles can be used to identify or target sub-populations of interest. For example, cells
involved in wound healing arise soon after induction, in early lineages 1, 3, and 5 only (Figure 7).
However, by late-stage senescence (Days 8, 10, 12), only one lineage (Late lineage 1a) expressed wound
65

healing signaling. CLEAR (coordinated lysosomal expression and regulation [52]) signaling was also
activated in early stage senescence lineages 5 and 6. It remained activated in middle stage lineages 1
and 3, and late stage lineage 1b, showing that lysosomal activity is upregulated quickly upon senescence
induction and remains active throughout the time course in specific lineages of cells. We found a
complex pattern of interleukin expression across the time course lineages. In early-stage senescence,
lineages 1-3 were characterized by varying degrees of activation of IL-1, 2, 3, 6, 7, 8, 10, and 17.
Lineages 4-6 overall showed inhibition of interleukin signaling, with exceptions in lineage 6 (IL-7
pathway). Interleukin signaling in middle stage senescence was activated sporadically in lineages 4-7,
while lineages 1-3 showed little activity. By late stage, less interleukin signaling pathways were
differentially regulated. For example, late lineage 1a showed increased IL-4, 6, and 12 activity. Together
with other lineages, interleukins 8 and 10 were also activated. Mitochondrial dysfunction was strongly
activated in early-stage lineages 5 and 6, which also showed high CLEAR activity. Interestingly, lineages
5 and 6 also showed inhibited EIF2 signaling in contrast to the activation observed in lineages 1-4. We
observed increased mitochondrial dysfunction signaling in middle-stage lineages 1,3, and 7, which also
activated CLEAR and chaperone-mediated autophagy signaling. Middle-stage lineage 4 was uniquely
characterized by activation of cholesterol biosynthesis, which carried forward specifically into late-stage
lineage 2a. This trajectory thus identifies a subpopulation of senescent cells engaged in altered
cholesterol metabolism. Thus, a complex pattern of pathway modulation describes the development of
senescence into separate lineages.
We identified distinct and dynamic transcriptional profiles for proliferating cells that evade doxo-
induced senescence and resist the impact of SnC conditioned media. This resistance may result from
several factors, including differences in individual cell state and DNA conformation during exposure, as
well as the effectiveness of doxo to enter cells given their spatial position and interactions with
neighboring cells. As IGFBP5 expression has broad biological implications, it is difficult to ascertain its
66

role in the early response (12 hours after doxo treatment) of senescence-escaping cells (Table 1). The
downregulation of MMP1 demonstrates an acute response in extracellular matrix (ECM) remodeling.
LRRC75A (upregulated in D10 vs D8) is a long noncoding RNA associated transcripts of which play
multiple roles in regulating proliferation and senescence. Interestingly, its overexpression in
spontaneous pre-term birth is associated with an increase in p53/p21 signaling and secretion of SASP
factors IL-8 and TNF-α [53]. We found that senescence evaders consistently modulate (either by
activation or inhibition) some pathways at every single time point, when compared to their growth
arrested neighbors (Table 2). Most pathways involved in cell cycle regulation and chromosomal
replication are activated, however some exceptions (ex. G2/M DNA damage checkpoint regulation) are
inhibited. Perhaps most strikingly, both nucleotide excision and base excision repair mechanisms are
consistently activated, suggesting that the proliferative capacity of senescence evaders involves a
sustained DNA repair response. Several immune-related pathways are involved as well. Specifically,
phagocytosis mediated by Fcγ and the immunogenic cell death signaling pathways are activated. In
contrast, granzyme A signaling is repressed, suggesting that senescence evaders are not priming
themselves for caspase-independent apoptosis. Finally, HMGB1 and VEGF signaling are upregulated,
both consistent with a proliferative state. These findings suggest that senescence evading cells walk a
fine line between proliferation, apoptosis, and senescence.
The microfluidics-based single cell processing strategy places some limitations on our data, due to the
enlarged and fragile state of SnCs. Some SnCs grow in size, and the largest of these cells may lyse upon
passing through the channels of the fluidic chips – releasing the contents of the cells into the medium.
Cells also experience shear stresses as they travel through chips, and thus may lyse even before the
encapsulation step. As a result, a caveat of our data is that the identified clusters and reconstructed
lineages are not fully representative of the doxo-induced SnC pool. Enlarged and/or fragile cells that are
not profiled form an understudied population that needs to be studied using different strategies.
67

Solutions to this problem include using nuclei instead of whole cells or using a non-microfluidic
approach such as Split-seq. Single nucleus profiling is attractive for several practical reasons such as
ease of processing and storage. However, the concordance between the transcriptional profile of a
whole cell compared to its nucleus varies by cell type and is not suitable in some cases [54-57].
Furthermore, collecting only the nuclear fraction of RNA prevents studying the fully mature population
of mRNAs in the cytosol. The Split-seq approach to single cell profiling eschews microfluidics altogether.
Instead of droplet encapsulation, fixed and permeabilized cells are used as the vessels for the reverse
transcription of mRNA into cDNA. This plate-based system circumvents cell size limitations and allows
for heterogeneous cells to be processed simultaneously. Combined with gentle pipetting, the lysis of
fragile cells can also be reduced this way, which may provide a more complete picture of the
transcriptomics of senescence.
Our findings increase the understanding of the scale and outcomes of heterogeneity in SnCs,
powered by single cell profiling and the inclusion of the temporal dimension. We show that even within
a mono-culture experimental context, senescent cells display remarkable variability in gene expression.
With the inclusion of the temporal dimension, as well as quiescent and proliferating controls, we
identify new markers for each cell state, and describe a dynamic transformation of proliferating cells to
senescence. We provide transcriptomic profiles of SnC lineages arising from doxo treatment, which
damages cells in a number of ways. These profiles can be used to identify populations of interest that
are characterized by different cellular responses. It can reasonably be expected that SnCs are highly
variable in vivo as well, due to the vast cellular insults continuously encountered by organisms, as well
as the cell type-specific nature of the senescence response. Our study provides a novel, time-resolved
single cell data set of senescence, and describes an analytical approach for teasing apart the various
strands of cellular senescence, laying the foundation for more challenging studies of this cellular
phenomenon in vivo.
68

CHAPTER 4: TREATMENT WITH THE SENOLYTIC 25HC IN OLD ANIMALS IMPROVES MEASURES OF
CARDIOVASCULAR AGING
Abstract
With aging, both cardiac and vascular health and function deteriorate owing to cellular mechanisms
including cellular senescence. SnCs in the heart and aortic vessel have been associated with worsened
measures of cardiovascular health including ejection fraction, fractional shortening, GLS strain, and
aortic pulse wave velocity. The cardiovascular SASP can promote the recruitment of immune cells,
activate fibroblasts, contributing to the development of fibrosis and tissue remodeling in the heart. This
can lead to impaired cardiac function and an increased risk of heart failure. In addition, cellular
senescence may also contribute to the development of age-related cardiac diseases, such as
atherosclerosis and cardiac hypertrophy. Senescent cells can accumulate in the arterial wall, leading to
chronic inflammation and the development of atherosclerotic plaques. In the heart, senescence can
contribute to the development of hypertrophy by promoting inflammation and fibrosis. Furthermore,
chemotherapeutic agents like doxo can also induce cellular senescence in the heart, which may
contribute to cardiotoxicity and an increased risk of heart failure. Anthracycline chemotherapy is a
model of accelerated vascular aging, used to study the effects of chemotherapy on blood vessel
function. Vascular dysfunction in this context is caused by the accumulation of reactive oxygen species
(ROS) and the activation of pro-inflammatory pathways in endothelial cells, leading to oxidative stress
and inflammation. Anthracycline chemotherapy can cause damage to the endothelial cells lining the
blood vessels, leading to impaired vasodilation and increased arterial stiffness. These changes are similar
to the changes seen in aging blood vessels and suggest that anthracycline chemotherapy may accelerate
the aging process of the vasculature. Previous findings indicate that senolytic interventions against age-
related and acutely induced SnCs lead to improvements in cardiac and vascular function. Skeletal muscle
mass and function decline with aging, culminating in sarcopenia, which is linked to an increased burden
69

of senescent cells. Muscle regeneration and maintenance are facilitated by resident mesenchymal
progenitors and muscle stem cells, also known as fibro-adipogenic progenitors (FAPs), and satellite cells
(SCs). Using single cell profiling, we identified CRYAB as a senescence-induced gene and potential target
for senolysis. Chemical inhibitor screening for CRYAB disruption identified 25-hydroxycholesterol (25HC),
an endogenous metabolite of cholesterol biosynthesis, as a potent senolytic. We further confirmed its
efficacy in diverse cell types of both human and mouse origin. 25HC is senolytic in mouse and human
fibroblasts, skeletal muscle cells, and cultured primary human cells from multiple organs. I hypothesized
that its senolytic activity would extend to the cardiovascular system and improve cardiac and vascular
health in aged animals. The study confirmed previous findings from multiple studies on GCV and
cardiovascular SnCs. Furthermore, improvements in ventricular function and aortic health, both critical
to healthy aging, were observed with 25HC treatment when compared to untreated and vehicle treated
animals. These results expand the range of efficiency of 25HC into new tissues and organs, and implicate
the senescence-associated expression of CRYAB in cardiovascular dysfunction. Further studies will
elucidate these findings through molecular and omic approaches.
Results
Three cohorts of aged (27-28 months old) p16-3MR mice were studied; cohorts “AI3” (n=7, 8, 9, and 7 of
25HC, GCV, HBCD, and untreated animals) and “AO1” (n=11, 9, and 11, of 25HC, GCV, and untreated
animals) were composed of females and cohort “AI4b” was composed of males (n=9 vehicle treated,
n=9 25HC treated). Treatments were delivered as one time IP injections. Echocardiography and pulse
wave doppler data was collected for two timepoints: before randomized grouping into treatment
cohorts, and 4 weeks after treatment. I found that 25HC treatment improved ejection fraction (EF) and
strain as measured by speckle-tracking echocardiography. EF is a measure of left ventricular (LV)
contraction, which forces oxygenated blood through the aortic valve for distribution throughout the
body. Impaired EF leads to decreased function either from injury, or maladaptive changes with age, can
70

induce symptoms of cardiovascular disease. Global longitudinal strain (GLS) is a well validated and
reproducible technique for LV longitudinal deformation, and an accurate predictor of cardiovascular
outcomes. Like any ejection-phase parameter, GLS is age-, sex-, and load-dependent. In adults (humans),
GLS <16% is abnormal, GLS >18% is normal, and GLS 16% to 18% is borderline.

Figure 4.1 – Heart rate measurements for 25HC, GCV, untreated, and 25HC Vehicle (HBCD) treated
groups in cohort AI3. No significant differences in heart rate were observed across the AI3 animals.

As animals were under isoflurane anesthesia and not otherwise stressed, cardiac load should not be
increased. We confirmed this by finding no significant differences between heart rate in 25HC, GCV,
untreated, and vehicle treated old 3MR females (Figure 4.1). Improvements in ejection fraction were
identified in longitudinal (before and after treatment) comparisons of the 25HC treatment group.
Furthermore, the senolytic treatments groups had statistically significant improvements in ejection
fraction compared to untreated and vehicle treated animals. 25HC and GCV treatment decreased GLS
longitudinally. 25HC treatment also resulted in significant improvements in GLS compared to untreated
and vehicle treated animals (Figure 4.2). We contextualized our understanding of aortic health through
PWV measurements of young untreated 3MR animals. Compared to the untreated baseline readings of
25HC-pre
25HC-post
GCV-pre
GCV-post
Untreated-pre
Untreated-post
Vehicle25HC-pre
Vehicle25HC-post
0
200
400
600
Heart Rate beats/minute
heart rate, beats per minute
71

cohort AO1, young animals had significantly lower values for PWV. With age, PWV increased from 300-
350 cm/s to up to ~450 cm/s. Senolytics improved PWV compared to the no treatment group with a
single dose. Despite the improvement, the metric did not return to youthful levels (Figure 4.3).

Figure 4.2 – Ejection fraction and GLS measurements improve with senolytic treatments in cohort AI3.
(A) Significant changes in ejection fraction were identified in longitudinal comparisons of the 25HC
treatment group. 25HC and GCV treated groups showed significant improvements in ejection
fraction compared to untreated and vehicle treated animals. (B) 25HC and GCV treatment
decreased GLS in longitudinal comparisons. 25HC treatment also resulted in significant
improvements in GLS compared to untreated and vehicle treated animals.

25HC-pre
25HC-post
GCV-pre
GCV-post
Untreated-pre
Untreated-post
Vehicle25HC-pre
Vehicle25HC-post
30
40
50
60
Ejection Fraction (EF)
pre/post treatment
✱ ✱
✱✱✱
✱✱
✱
✱
25HC-pre
25HC-post
GCV-pre
GCV-post
Untreated-pre
Untreated-post
Vehicle25HC-pre
Vehicle25HC-post
-30
-20
-10
0
Global Longitudinal Strain (GLS)
✱ ✱✱
✱
✱✱✱
✱
72

Figure 4.3 – Aged animals in cohort AO1 show higher PWV compared to young control 3MR females.
Compared to young untreated females, animals in cohort AO1 have higher pulse wave velocities,
indicating a decrease in aortic health.

In cohort AO1, GCV and 25HC treated animals showed significant improvements in PWV compared to
untreated controls (Figure 4.3). However in comparison to AI3, significant longitudinal changes were not
observed, indicating the need for investigation into the variability in response to 25HC treatment (Figure
4.4). Potential reasons for this variability include heterogeneity in life history and aging, as well as timing
of repeated measures post-treatment.

Figure 4.4 – Longitudinal comparisons of AO1 cohorts show an improving trend in aortic health with
senolytic treatment. Untreated control animals showed an increasing but not significant trend in PWV,
between the baseline and timepoint 1 measurements. 25HC and GCV treated groups had decreasing
trends between baseline and treatment timepoint 1.
young
old baseline
GCV
25HC
control
250
300
350
400
450
Young controls vs AO1
PWV (cm/s)
✱ ✱ ✱
✱ ✱
ns
✱ ✱ ✱
73

Figure 4.5 – Longitudinal comparisons of AI3 cohort groups show significant improvements in aortic
health with senolytic treatment. Untreated animals saw a significant increase in PWV between
timepoints, in contrast to vehicle treatment. 25HC and GCV treated animals saw significant
improvements in PWV compared to pre-treatment.

Figure 4.6 – 25HC treatment in old 3MR male animals reduces aortic elastic modulus. Assays carried
out by collaborators in the Seals laboratory at UC Boulder showed a significant improvement in aortic
elastic modulus in old 3MR males following 25HC treatment.

74

Discussion
The elimination of senescent cells can alleviate vasomotor dysfunction in naturally aging mice and mice
with existing atherosclerosis. In addition, removing senescent cells reduces the presence of osteogenesis
markers in advanced intimal plaques, ultimately decreasing the calcification of these plaques. In one
study, the senolytic compound ABT-263 reversed PWV to young control levels, supported by improved
elastic modulus and improved endothelial cell function as measured by endothelial-dependent dilation.
Based on these results, clearance of senescent cells could be a beneficial and complementary approach
to traditional risk factor management in reducing the morbidity and mortality rates associated with
cardiovascular diseases related to aging. In cohort AI3, longitudinal comparisons of GCV and 25HC
treated animals showed significant improvements in PWV compared to pre-treatment (Figure 4.5). In
contrast, the PWVs of untreated animals significantly increased while treatment with the 25HC vehicle
resulted in a non-significant increase. 25HC treatment thus results in improvements in aortic health
comparable to those caused by GCV treatment. In studies done by our collaborators in the Seals lab at
UC Boulder, 25HC treatment led to significant improvement in aortic elastic modulus in old 3MR male,
compared to vehicle treated animals (Figure 4.6). This assay measures the intrinsic mechanical stiffness
of aortic tissue by incremental stress-strain testing via wire myography, and serves as validation of PWV
measurements. GCV-mediated killing of SnCs is limited to the p16-3MR transgenic mouse model,
inherently limiting its usefulness in other mouse models as well as human cells and any translational
studies. In contrast, 25HC is a naturally occurring hydroxycholesterol conserved between humans and
mice, making it an attractive target for development as a broad use senolytic agent. The data presented
here expand the role of 25HC as a compound that can ameliorate aging and senescence associated
dysfunction in the cardiovascular system.

75

CHAPTER 5: CONCLUSIONS
Altogether, the findings described here offer new insights into the development of senescence and our
understanding of SnCs as dynamic and heterogenous cell types. Our studies provide new insights into
the senolytic potential of 25HC, as well as a novel data set and methodology for analyzing senescent
single cells over a time course. From the transcriptional diversity of studied SnC endpoints, a picture
emerges of a cell state that is difficult to encapsulate or define with a single word. It may be more
meaningful to think of SnCs as possessing “senotypes”, which diverge over time from an induction
event. By studying this cellular phenomenon at the single cell level in vitro, using one cell type and one
inducer, it is possible to reduce the problem of SnC heterogeneity to a more manageable dimensionality.
Indeed, despite using only doxo as an inducer, it is possible to show that cells embark on journeys with
ultimately differing destinations. The initial cell state is likely a very important determinant to whether it
will transform into a mitochondrially dysfunctional, or DNA damage driven, or wound healing signaling
SnC. Further studies are needed to capture and understand the senescence outcomes of different initial
cell states – these may consist of cell cycle stage, metabolic activity, or the surrounding biological
context. A given cell’s place in the cell cycle can determine both the amount and conformation of its
DNA, and thus how it would respond to a senescence induction event. Likewise, mitochondrial content
and activity may also play a role in the development of senescence as a result of mitochondrial
dysfunction. While this study was conducted in a relatively simple monoculture experiment, the results
nonetheless provide insight into the transcriptional heterogeneity of SnCs, as well as providing new
marker genes for identifying SnCs compared to quiescent or proliferating versions of the same cell type
– a comparative context which is often ignored. Our study shows that the transformation to senescence
is a process that cells undergo at varying speeds, adding another level of complexity to SnCs – apart
from asking which cells are harmful or otherwise interesting, we need to ask when exactly these cells
appear, in order to tune any intervention to the right timing. These findings are important as identifying
76

genes with expression specificity in SnCs provides new targets for investigation, including senolytic and
senomorphic potential. This applies not only to the 3 mentioned cell states, but also to differentiating
between the identified senotypes. SnCs play an important role in wound healing, and in this study I
showed that cells engaged in wound healing signaling can arise from a pool of doxo treated cells,
alongside others that may be engaged in less benign activities. The development of molecular markers
that allow us to separate (by FACS for example) and further study senotypes is an important avenue of
research that serves to increase our understanding of SnCs and their weaknesses. Apart from simply
providing potential identifiers by which to categorize and investigate SnCs, our studies also serve as a
launching point for the development of more specific SnC-targeting compounds (as shown in our
discovery of 25HC). The experimental and computational approaches described herein also serve as a
model for conducting analogous studies in any cell type of interest, and can be modified and expanded
to deal with a more complex population of SnCs such as that derived from primary tissues. By
incorporating the temporal dimension into a single cell study of senescence, our data provide a new lens
by which to understand this cell state, which plays a critical role in the aging process. With a deeper and
more specific arsenal of tools for engaging with SnCs, we can be hopeful that future therapeutic
interventions will be both more effective and carry reduced risks or off-target interactions.

77

References

1. Bodnar, A.G., et al., Extension of life-span by introduction of telomerase into normal human cells.
Science, 1998. 279(5349): p. 349-52.
2. Campisi, J. and F. d'Adda di Fagagna, Cellular senescence: when bad things happen to good cells.
Nat Rev Mol Cell Biol, 2007. 8(9): p. 729-40.
3. Gorgoulis, V., et al., Cellular Senescence: Defining a Path Forward. Cell, 2019. 179(4): p. 813-827.
4. Serrano, M., et al., Oncogenic ras provokes premature cell senescence associated with
accumulation of p53 and p16INK4a. Cell, 1997. 88(5): p. 593-602.
5. Aguayo-Mazzucato, C., et al., Acceleration of beta Cell Aging Determines Diabetes and Senolysis
Improves Disease Outcomes. Cell Metab, 2019. 30(1): p. 129-142 e4.
6. Baker, D.J., et al., Naturally occurring p16(Ink4a)-positive cells shorten healthy lifespan. Nature,
2016. 530(7589): p. 184-9.
7. Bussian, T.J., et al., Clearance of senescent glial cells prevents tau-dependent pathology and
cognitive decline. Nature, 2018. 562(7728): p. 578-582.
8. Campisi, J., Cellular senescence as a tumor-suppressor mechanism. Trends Cell Biol, 2001.
11(11): p. S27-31.
9. Childs, B.G., et al., Senescent intimal foam cells are deleterious at all stages of atherosclerosis.
Science, 2016. 354(6311): p. 472-477.
10. Chinta, S.J., et al., Cellular Senescence Is Induced by the Environmental Neurotoxin Paraquat and
Contributes to Neuropathology Linked to Parkinson's Disease. Cell Rep, 2018. 22(4): p. 930-940.
11. Demaria, M., et al., Cellular Senescence Promotes Adverse Effects of Chemotherapy and Cancer
Relapse. Cancer Discov, 2017. 7(2): p. 165-176.
12. Demaria, M., et al., An essential role for senescent cells in optimal wound healing through
secretion of PDGF-AA. Dev Cell, 2014. 31(6): p. 722-33.
13. Farr, J.N., et al., Targeting cellular senescence prevents age-related bone loss in mice. Nat Med,
2017. 23(9): p. 1072-1079.
14. Jeon, O.H., et al., Local clearance of senescent cells attenuates the development of post-
traumatic osteoarthritis and creates a pro-regenerative environment. Nat Med, 2017. 23(6): p.
775-781.
15. Menon, R., et al., Placental membrane aging and HMGB1 signaling associated with human
parturition. Aging (Albany NY), 2016. 8(2): p. 216-30.
16. Munoz-Espin, D., et al., Programmed cell senescence during mammalian embryonic
development. Cell, 2013. 155(5): p. 1104-18.
17. Ohtani, N., D.J. Mann, and E. Hara, Cellular senescence: its role in tumor suppression and aging.
Cancer Sci, 2009. 100(5): p. 792-7.
18. Storer, M., et al., Senescence is a developmental mechanism that contributes to embryonic
growth and patterning. Cell, 2013. 155(5): p. 1119-30.
78

19. Kumari, R. and P. Jat, Mechanisms of Cellular Senescence: Cell Cycle Arrest and Senescence
Associated Secretory Phenotype. Front Cell Dev Biol, 2021. 9: p. 645593.
20. Hernandez-Segura, A., et al., Unmasking Transcriptional Heterogeneity in Senescent Cells. Curr
Biol, 2017. 27(17): p. 2652-2660 e4.
21. van Deursen, J.M., The role of senescent cells in ageing. Nature, 2014. 509(7501): p. 439-46.
22. Acosta, J.C., et al., Chemokine signaling via the CXCR2 receptor reinforces senescence. Cell, 2008.
133(6): p. 1006-18.
23. Coppe, J.P., et al., Senescence-associated secretory phenotypes reveal cell-nonautonomous
functions of oncogenic RAS and the p53 tumor suppressor. PLoS Biol, 2008. 6(12): p. 2853-68.
24. Kuilman, T., et al., Oncogene-induced senescence relayed by an interleukin-dependent
inflammatory network. Cell, 2008. 133(6): p. 1019-31.
25. Kuilman, T. and D.S. Peeper, Senescence-messaging secretome: SMS-ing cellular stress. Nat Rev
Cancer, 2009. 9(2): p. 81-94.
26. Nelson, G., et al., A senescent cell bystander effect: senescence-induced senescence. Aging Cell,
2012. 11(2): p. 345-9.
27. Acosta, J.C., et al., A complex secretory program orchestrated by the inflammasome controls
paracrine senescence. Nat Cell Biol, 2013. 15(8): p. 978-90.
28. The gene expression program of prostate fibroblast senescence modulates neoplastic epithelial
cell proliferation through paracrine mechanisms.
29. Rodier, F., et al., Persistent DNA damage signalling triggers senescence-associated inflammatory
cytokine secretion. Nat Cell Biol, 2009. 11(8): p. 973-9.
30. Novakova, Z., et al., Cytokine expression and signaling in drug-induced cellular senescence.
Oncogene, 2010. 29(2): p. 273-84.
31. Aging, cellular senescence, and cancer.
32. Niedernhofer, L.J. and P.D. Robbins, Senotherapeutics for healthy ageing. Nat Rev Drug Discov,
2018. 17(5): p. 377.
33. Chang, J., et al., Clearance of senescent cells by ABT263 rejuvenates aged hematopoietic stem
cells in mice. Nat Med, 2016. 22(1): p. 78-83.
34. Hall, B.M., et al., p16(Ink4a) and senescence-associated beta-galactosidase can be induced in
macrophages as part of a reversible response to physiological stimuli. Aging (Albany NY), 2017.
9(8): p. 1867-1884.
35. Wiley, C.D., et al., Analysis of individual cells identifies cell-to-cell variability following induction
of cellular senescence. Aging Cell, 2017. 16(5): p. 1043-1050.
36. Kim, Y.M., et al., Implications of time-series gene expression profiles of replicative senescence.
Aging Cell, 2013. 12(4): p. 622-34.
37. Schauble, S., et al., Quantitative model of cell cycle arrest and cellular senescence in primary
human fibroblasts. PLoS One, 2012. 7(8): p. e42150.
79

38. Velarde, M.C. and M. Demaria, Targeting Senescent Cells: Possible Implications for Delaying Skin
Aging: A Mini-Review. Gerontology, 2016. 62(5): p. 513-8.
39. Wiley, C.D., et al., Mitochondrial Dysfunction Induces Senescence with a Distinct Secretory
Phenotype. Cell Metab, 2016. 23(2): p. 303-14.
40. Mitochondrial dysfunction in cell senescence and aging.
41. Bientinesi, E., et al., Doxorubicin-induced senescence in normal fibroblasts promotes in vitro
tumour cell growth and invasiveness: The role of Quercetin in modulating these processes. Mech
Ageing Dev, 2022. 206: p. 111689.
42. Eom, Y.W., et al., Two distinct modes of cell death induced by doxorubicin: apoptosis and cell
death through mitotic catastrophe accompanied by senescence-like phenotype. Oncogene, 2005.
24(30): p. 4765-77.
43. Sliwinska, M.A., et al., Induction of senescence with doxorubicin leads to increased genomic
instability of HCT116 cells. Mech Ageing Dev, 2009. 130(1-2): p. 24-32.
44. Yang, M.-Y., et al., Induction of Senescence by Doxorubicin Is Associated with Upregulated Mir-
375 and Induction of Autophagy In CML Cell Line K562. Blood, 2010. 116(21): p. 1843-1843.
45. Adjemian, S., et al., Ionizing radiation results in a mixture of cellular outcomes including mitotic
catastrophe, senescence, methuosis, and iron-dependent cell death. Cell Death Dis, 2020. 11(11):
p. 1003.
46. Chen, Z., et al., Cellular senescence in ionizing radiation (Review). Oncol Rep, 2019. 42(3): p. 883-
894.
47. Li, M., et al., Ionizing Radiation-Induced Cellular Senescence in Normal, Non-transformed Cells
and the Involved DNA Damage Response: A Mini Review. Front Pharmacol, 2018. 9: p. 522.
48. Saleh, T., et al., Therapy-Induced Senescence: An "Old" Friend Becomes the Enemy. Cancers
(Basel), 2020. 12(4).
49. Meredith, A.M. and C.R. Dass, Increasing role of the cancer chemotherapeutic doxorubicin in
cellular metabolism. J Pharm Pharmacol, 2016. 68(6): p. 729-41.
50. Bielak-Zmijewska, A., et al., A comparison of replicative senescence and doxorubicin-induced
premature senescence of vascular smooth muscle cells isolated from human aorta.
Biogerontology, 2014. 15(1): p. 47-64.
51. Mitry, M.A., et al., Accelerated cardiomyocyte senescence contributes to late-onset doxorubicin-
induced cardiotoxicity. Am J Physiol Cell Physiol, 2020. 318(2): p. C380-C391.
52. Palmieri, M., et al., Characterization of the CLEAR network reveals an integrated control of
cellular clearance pathways. Human Molecular Genetics, 2011. 20(19): p. 3852-3866.
53. Jiang, J., et al., Long non-coding RNA SNHG29 regulates cell senescence via p53/p21 signaling in
spontaneous preterm birth. Placenta, 2021. 103: p. 64-71.
54. Lake, B.B., et al., A comparative strategy for single-nucleus and single-cell transcriptomes
confirms accuracy in predicted cell-type expression from nuclear RNA. Sci Rep, 2017. 7(1): p.
6031.
80

55. Chamberlin, J.T., et al., Variable RNA sampling biases mediate concordance of single-cell and
nucleus sequencing across cell types. 2022.
56. Wen, F., et al., Comparison of single‑ nucleus and single‑ cell transcriptomes in hepatocellular
carcinoma tissue. Mol Med Rep, 2022. 26(5).
57. Thrupp, N., et al., Single-Nucleus RNA-Seq Is Not Suitable for Detection of Microglial Activation
Genes in Humans. Cell Rep, 2020. 32(13): p. 108189.

81

APPENDICES

APPENDIX A

Seurat_Workflow.R – R code used to pre-process, quality control, normalize and cluster time course
libraries using the Seurat R package.
.libPaths("/opt/R/4.0.4/lib/R/library")
Sys.unsetenv("GITHUB_PAT")
setwd("/data/scRNA/IMR90_timecourse/")

library(dplyr)
library(purrr)
library(monocle3)
library(Seurat)
library(SeuratWrappers)
library(SeuratDisk)
library(patchwork)
library(ggplot2)
library(patchwork)
library(harmony)
require("Seurat.utils")
library(scCustomize)
library(psupertime)
require(devtools)
library(scanalysis)
library(SingleCellExperiment)
library(Tempora)
library(basicPlotteR)
library(igraph)
library(RColorBrewer)
library(ggrepel)
library(SCpubr)
library(paletteer)
library(glmGamPoi)
library(dittoSeq)
library(shadowtext)
library(magick)
library(ggvennEx)
library(ggbio)

## load doublet cleaned 10x libraries
{

#singlets only Seurat objects
Pro <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_Pro.RDS.gz")
H12 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_H12.RDS.gz")
D1_1 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D1_1.RDS.gz")
D1_2 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D1_2.RDS.gz")
D2_1 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D2_1.RDS.gz")
D2_2 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D2_2.RDS.gz")
D4_1 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D4_1.RDS.gz")
D4_2 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D4_2.RDS.gz")
D6_1 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D6_1.RDS.gz")
D6_2 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D6_2.RDS.gz")
D8_1 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D8_1.RDS.gz")
D8_2 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D8_2.RDS.gz")
D10_1 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D10_1.RDS.gz")
D10_2 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D10_2.RDS.gz")
82

D12_1 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D12_1.RDS.gz")
D12_2 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D12_2.RDS.gz")
Qui_1 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_Qui_1.RDS.gz")
Qui_2 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_Qui_2.RDS.gz")

}

seurlist <- list(Pro, H12, D1_1, D1_2, D2_1, D2_2, D4_1, D4_2, D6_1, D6_2, D8_1, D8_2, D10_1, D10_2, D12_1, D12_2)
saveRDS.compress.in.BG(seurlist, fname = "/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/seurlist_preFilter_preSCT_11-02-
22.RDS")
seurlist <- readRDS("/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/seurlist_preFilter_preSCT_11-02-22.RDS.gz")

seuratsubset_list <- list()
sct_list <- list()
libraryName_list <-
list("Pro","H12","D1_1","D1_2","D2_1","D2_2","D4_1","D4_2","D6_1","D6_2","D8_1","D8_2","D10_1","D10_2","D12_1","D12_2")
dayName_list <- list("Pro","H12","D1","D1","D2","D2","D4","D4","D6","D6","D8","D8","D10","D10","D12","D12")
batch_list <- list("Rep3", "Rep3", "Rep1", "Rep2", "Rep1", "Rep2","Rep1", "Rep2", "Rep1", "Rep2","Rep1", "Rep2", "Rep1", "Rep2","Rep1",
"Rep2")

#using libraryName_list with Quiescent sample commented out, so that the sample will be skipped for timecourse
for (f in seq_along(along.with = libraryName_list)){

print(f)
#Aready run if commented:

seurlist[[f]][["library"]] <- libraryName_list[[f]]
seurlist[[f]][["timepoint"]] <- dayName_list[[f]]
seurlist[[f]][["batch"]] <- batch_list[[f]]
seurlist[[f]][["percent.mito"]] <- PercentageFeatureSet(seurlist[[f]], pattern = "^MT-")
seurlist[[f]][["percent.ribo"]] <- PercentageFeatureSet(seurlist[[f]], pattern = "^RP[SL]")
seurlist[[f]][["percent.MALAT1"]] <- PercentageFeatureSet(seurlist[[f]], pattern = "MALAT1")
seurlist[[f]][["percent.NEAT1"]] <- PercentageFeatureSet(seurlist[[f]], pattern = "NEAT1")

#print(unique(seurlist[[f]]$timepoint))
print(seurlist[[f]])
seuratsubset_list[[f]] <- subset(seurlist[[f]], percent.MALAT1 < 20 & percent.NEAT1 < 20 &
nFeature_RNA > 500 & nCount_RNA > 2000 & percent.mito < 25 & percent.ribo < 25)
print(seuratsubset_list[[f]])

#plotlist_postfilter[[f]] <- (VlnPlot(seuratsubset_list[[f]], features = c("nFeature_RNA", "nCount_RNA", "percent.mito", "percent.ribo",
"nFeature_spliced","nCount_spliced","nFeature_unspliced","nCount_unspliced", "percent.MALAT1","percent.NEAT1"), ncol = 4, pt.size=0.05))
#sct_list[[f]] <- SCTransform(seuratsubset_list[[f]], variable.features.n = 3000,
# vars.to.regress = c("percent.mito", "percent.ribo"), return.only.var.genes = F, verbose = TRUE)
}

saveRDS.compress.in.BG(seuratsubset_list, fname =
"/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/seurlist_postFilter_preSCT_11-02-22.RDS")

################################### Merge all objects -> SCT -> clusters

seuratsubset_list <- readRDS("/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/seurlist_postFilter_preSCT_11-02-22.RDS")

TPs_merged <- merge(seuratsubset_list[[1]], y = c(seuratsubset_list[[2]], seuratsubset_list[[3]], seuratsubset_list[[4]], seuratsubset_list[[5]],
seuratsubset_list[[6]], seuratsubset_list[[7]], seuratsubset_list[[8]], seuratsubset_list[[9]], seuratsubset_list[[10]], seuratsubset_list[[11]],
seuratsubset_list[[12]], seuratsubset_list[[13]], seuratsubset_list[[14]], seuratsubset_list[[15]], seuratsubset_list[[16]]), merge.data = T)
TPs_merged_SCT <- SCTransform(TPs_merged, vst.flavor = "v2", ncells = 15000, variable.features.n = 4000, vars.to.regress = c("percent.mito",
"percent.ribo","percent.MALAT1","percent.NEAT1"), return.only.var.genes = T, verbose = T)

83

saveRDS.compress.in.BG(TPs_merged_SCT, fname =
"/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/seurlist_postFilter_postSCT_11-02-22.RDS")

TPs_merged_SCT <- readRDS("/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/seurlist_postFilter_postSCT_11-02-22.RDS.gz")

Idents(TPs_merged_SCT) <- TPs_merged_SCT$timepoint

ordered_Timepoints <- c("Pro","H12","D1","D2","D4","D6","D8","D10","D12")
TPs_merged_SCT$timepoint <- factor(TPs_merged_SCT$timepoint, levels= ordered_Timepoints)

TPs_merged_SCT <- RunPCA(TPs_merged_SCT, npcs = 100, verbose = F, ndims.print = F)
ElbowPlot(TPs_merged_SCT, ndims = 50)
#DefaultAssay(TPs_merged_SCT)
TPs_merged_SCT <- RunUMAP(TPs_merged_SCT, repulsion.strength = 5, assay = 'SCT', dims = 1:50)
TPs_merged_SCT <- FindNeighbors(TPs_merged_SCT, reduction = "umap", dims = 1:2)
TPs_merged_SCT <- FindClusters(TPs_merged_SCT, resolution = 0.35)

DimPlot(TPs_merged_SCT, group.by = "timepoint", pt.size = 0.001, raster = F, shuffle = T)

####
main_title_size = 18
axis_text_size = 9
axis_title_size = 9
panel_tag_size = 16
####

ambientRNA <- image_read("/data/scRNA/IMR90_timecourse/FractionReadsinCells_barGraph.tif") %>% image_ggplot()
ambientRNA <- ggbitmap("/data/scRNA/IMR90_timecourse/FractionReadsinCells_barGraph.tif")
ambientRNA2 <- ggbio::rescale(ambientRNA, sx = 1.5, sy = 1.5)

Idents(TPs_merged_SCT) <- TPs_merged_SCT$batch
replicate1_plot <- DimPlot(TPs_merged_SCT, group.by = "timepoint", cells = WhichCells(TPs_merged_SCT, idents = "Rep1"), pt.size = 0.001,
raster = F, shuffle = T, combine = T) +
NoLegend() + #comment out and edit legend placement manually as needed
#guides(color=guide_legend(title="Timepoint"), fill = guide_legend(override.aes = list(size=20))) +
ggtitle('Replicate 1') +
labs(tag = "A") +
theme(text = element_text(size = main_title_size, face = "plain"), axis.title=element_text(size=axis_title_size), axis.text =
element_text(size = axis_text_size), plot.tag = element_text(size = panel_tag_size, face = "bold"))

replicate2_plot <- DimPlot(TPs_merged_SCT, group.by = "timepoint", cells = WhichCells(TPs_merged_SCT, idents = "Rep2"), pt.size = 0.001,
raster = F, shuffle = T, combine = T) +
NoLegend() + #comment out and edit legend placement manually as needed
#guides(color=guide_legend(title="Timepoint"), fill = guide_legend(override.aes = list(size=20))) +
ggtitle('Replicate 2') +
#labs(tag = "C") +
theme(text = element_text(size = main_title_size, face = "plain"), axis.title=element_text(size=axis_title_size), axis.text =
element_text(size = axis_text_size), plot.tag = element_text(size = panel_tag_size, face = "bold"))

replicateOverlap_plot <- DimPlot(TPs_merged_SCT, group.by = "batch", cells = WhichCells(TPs_merged_SCT, idents = "Rep3", invert = T), pt.size
= 0.001, raster = F, shuffle = T, combine = T) +
NoLegend() + #comment out and edit legend placement manually as needed
#guides(color=guide_legend(title="Timepoint"), fill = guide_legend(override.aes = list(size=20))) +
ggtitle('Replicates Overlaid') +
labs(tag = "C") +
theme(text = element_text(size = main_title_size, face = "plain"), axis.title=element_text(size=axis_title_size), axis.text =
element_text(size = axis_text_size), plot.tag = element_text(size = panel_tag_size, face = "bold"))

replicateOverlap_plot2 <- DimPlot(TPs_merged_SCT, group.by = "batch", cells = WhichCells(TPs_merged_SCT, idents = "Rep3", invert = T),
pt.size = 0.001, raster = F, shuffle = T, combine = T) +
NoLegend() + #comment out and edit legend placement manually as needed
#guides(color=guide_legend(title="Timepoint"), fill = guide_legend(override.aes = list(size=20))) +
ggtitle('Overlaid') +
84

#labs(tag = "C") +
theme(text = element_text(size = main_title_size, face = "plain"), axis.title=element_text(size=axis_title_size), axis.text =
element_text(size = axis_text_size), plot.tag = element_text(size = panel_tag_size, face = "bold"))

(replicate1_plot | replicate2_plot) / (replicateOverlap_plot | replicateOverlap_plot2)

########################################

axis_title_size = 10
axis_text_size = 9
panel_tag_size = 14

AllTPs_p1 <- DimPlot(TPs_merged_SCT, group.by = "timepoint", cols = paletteer_d("ggthemes::Classic_20", n = 10), shuffle = T, raster = F) +
#NoLegend() + #comment out and edit legend placement manually as needed
#guides(color=guide_legend(title="Timepoint"), fill = guide_legend(override.aes = list(size=20))) +
ggtitle('') +
labs(tag = "A") +
theme(legend.direction = "horizontal", legend.position = "bottom", legend.key.size = unit(0.5,"cm"), legend.justification = "center",
axis.title = element_text(size = axis_title_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold")) +
guides(colour = guide_legend(nrow = 1, override.aes = list(size=7.5))) #force legend to one row, change symbol size

AllTPs_p2 <- DimPlot(TPs_merged_SCT, group.by = "SCT_snn_res.0.35", cols = paletteer_d("ggsci::default_igv", n = 46, direction = 1), pt.size =
0.01, label = T, label.box = T, label.size = 4, raster = F) +
NoLegend() +
ggtitle('') +
labs(tag = "B") +
theme(axis.title = element_text(size = axis_title_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold"))

AllTPs_p1 / AllTPs_p2

#pdf("/data/scRNA/IMR90_timecourse/Figures_final/UMAP_byTimepoint_and_byCluster.pdf")
#P_timepoint / P_cluster
#dev.off()

ordered_Timepoints <- c("Pro","H12","D1","D2","D4","D6","D8","D10","D12")
TPs_merged_SCT$timepoint <- factor(TPs_merged_SCT$timepoint, levels= ordered_Timepoints)

dittoBarPlot(TPs_merged_SCT, var = "timepoint", group.by = "SCT_snn_res.0.35", var.labels.reorder = c(9, 8, 1, 4, 5, 6, 7, 2, 3), color.panel =
paletteer_d("ggthemes::Classic_20", n = 10),
x.reorder =
c(23,36,38,43,18,21,41,14,15,22,26,37,30,9,33,13,25,4,31,42,39,44,35,29,8,24,5,6,12,19,27,2,45,46,3,7,11,16,17,34,28,40,32,20,10,1),
xlab = "Cluster number", x.labels.rotate = T, legend.show = F) + #main = "Cluster composition by Timepoint",
#labs(tag = "A") +
theme(legend.direction = "horizontal", legend.position = "bottom", legend.key.size = unit(0.5,"cm"), legend.justification = "center",
axis.title = element_text(size = axis_title_size, face = "bold"), axis.text = element_text(size = axis_text_size, face = "bold", color =
"black"), plot.tag = element_text(size = panel_tag_size, face = "bold")) +
guides(fill=guide_legend(nrow=1, byrow=TRUE, title = "Timepoint"))

#dittoBarPlot(TPs_merged_SCT, var = "SCT_snn_res.0.35", group.by = "timepoint")
#table(Idents(TPs_merged_SCT), TPs_merged_SCT$SCT_snn_res.0.35)

FeaturePlot_sengenes <- FeaturePlot(TPs_merged_SCT, label = F, label.size = 3, label.color = "black", features = c("MKI67", "TOP2A", "CCNA2",
"PCNA", "HMGB1", "HMGB2", "LMNB1", "EZH2", "CDKN1A", "CDKN2A", "TP53", "GDF15"),
85

keep.scale = "feature", slot = "data", pt.size = .5, cols = c("lightgrey","#481b6d"), raster = T) &
NoLegend() &
#ggtitle('') &
#labs(tag = "B") &
theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag =
element_text(size = panel_tag_size, face = "bold"))#,
#plot.title = element_text(size=15), legend.key.size = unit(2.75, 'mm'),
#legend.position = c(0.85, 0.12)) #legend.position="right",

FeaturePlot_sengenes

saveRDS.compress.in.BG(TPs_merged_SCT, fname =
"/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/TPs_merged_SCT_final.RDS.gz")

TPs_merged_SCT <- readRDS("/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/TPs_merged_SCT_final.RDS.gz")

Idents(TPs_merged_SCT) <- TPs_merged_SCT$timepoint

VlnPlot(TPs_merged_SCT, features = c("MKI67", "CDKN1A"), pt.size = 0) +
stat_summary(fun.y = median, geom='point', size = 25, colour = "black", shape = 95)

# Proliferating cells analysis

TPs_merged_SCT <- readRDS("/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/TPs_merged_SCT_final.RDS.gz")

prolif1 <- subset(TPs_merged_SCT, (MKI67 > 0 & CCNA2 > 0) | (MKI67 > 0 & TOP2A > 0) | (TOP2A > 0 & CCNA2 > 0) | (MCM2 > 0 & MKI67 > 0) |
(MCM2 > 0 & TOP2A > 0) | (MCM2 > 0 & CCNA2 > 0) |
(PCNA > 0 & MKI67 > 0) | (PCNA > 0 & TOP2A > 0) | (PCNA > 0 & CCNA2 > 0) | (PCNA > 0 & MCM2 > 0))

#MKI67 CCNA2 TOP2A MCM2 PCNA
#DefaultAssay(TPs_merged_SCT)

#prolifInvert <- subset(TPs_merged_SCT, invert = T, (MKI67 > 0 & CCNA2 > 0) | (MKI67 > 0 & TOP2A > 0) | (TOP2A > 0 & CCNA2 > 0) | (MCM2 >
0 & MKI67 > 0) | (MCM2 > 0 & TOP2A > 0) | (MCM2 > 0 & CCNA2 > 0) )

sct_prolif <- SCTransform(prolif1, vst.flavor = "v2", vars.to.regress = c("percent.mito", "percent.ribo","percent.MALAT1","percent.NEAT1"),
return.only.var.genes = F, verbose = T)

sct_prolif <- RunPCA(sct_prolif, npcs = 100)
ElbowPlot(sct_prolif, ndims = 100)
sct_prolif <- RunUMAP(sct_prolif, assay = 'SCT', dims = 1:30, min.dist = 0.5, spread = 2.5)
sct_prolif <- FindNeighbors(sct_prolif, reduction = "umap", dims = 1:2)
sct_prolif <- FindClusters(sct_prolif, resolution = 0.35)

main_title_size = 18
axis_text_size = 15
axis_title_size = 15
panel_tag_size = 24

prolif_timepoints_p1 <- DimPlot(sct_prolif, pt.size = 0.05, group.by = "timepoint", cols = paletteer_d("ggthemes::Classic_20", n = 10), label = F)
+ #labs(tag = "A") +
NoLegend() + ggtitle('') + #ylim(min_y, max_y) +
theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold"))

prolif_timepoints_p1

#prolif_clusters_p2 <- DimPlot(sct_prolif, pt.size = 0.05, group.by = "SCT_snn_res.0.1", cols = paletteer_d("ggsci::default_igv", n = 46, direction
= 1), label = T, label.box = T, label.size = 4) +
86

# labs(tag = "B") +
# NoLegend() + ggtitle('') + #ylim(min_y, max_y) +
# theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold"))

#compare prolif vs non-prolif in each TP
{

Prolif_cells <- WhichCells(prolif1)
TPs_merged_SCT$prolif <- ifelse(colnames(TPs_merged_SCT) %in% Prolif_cells, "P", "N")

Idents(TPs_merged_SCT) <- TPs_merged_SCT$timepoint
H12_PvsNp <- FindMarkers(TPs_merged_SCT, ident.1 = "P", ident.2 = "N", group.by = "prolif", subset.ident = "H12", assay = "SCT", min.pct =
0.05, only.pos = F)
H12_PvsNp$cluster <- "H12_PvsNp"
H12_PvsNp$gene <- rownames(H12_PvsNp)

D1_PvsNp <- FindMarkers(TPs_merged_SCT, ident.1 = "P", ident.2 = "N", group.by = "prolif", subset.ident = "D1", assay = "SCT", min.pct = 0.05,
only.pos = F)
D1_PvsNp$cluster <- "D1_PvsNp"
D1_PvsNp$gene <- rownames(D1_PvsNp)

D2_PvsNp <- FindMarkers(TPs_merged_SCT, ident.1 = "P", ident.2 = "N", group.by = "prolif", subset.ident = "D2", assay = "SCT", min.pct = 0.05,
only.pos = F)
D2_PvsNp$cluster <- "D2_PvsNp"
D2_PvsNp$gene <- rownames(D2_PvsNp)

D4_PvsNp <- FindMarkers(TPs_merged_SCT, ident.1 = "P", ident.2 = "N", group.by = "prolif", subset.ident = "D4", assay = "SCT", min.pct = 0.05,
only.pos = F)
D4_PvsNp$cluster <- "D4_PvsNp"
D4_PvsNp$gene <- rownames(D4_PvsNp)

D6_PvsNp <- FindMarkers(TPs_merged_SCT, ident.1 = "P", ident.2 = "N", group.by = "prolif", subset.ident = "D6", assay = "SCT", min.pct = 0.05,
only.pos = F)
D6_PvsNp$cluster <- "D6_PvsNp"
D6_PvsNp$gene <- rownames(D6_PvsNp)

D8_PvsNp <- FindMarkers(TPs_merged_SCT, ident.1 = "P", ident.2 = "N", group.by = "prolif", subset.ident = "D8", assay = "SCT", min.pct = 0.05,
only.pos = F)
D8_PvsNp$cluster <- "D8_PvsNp"
D8_PvsNp$gene <- rownames(D8_PvsNp)

D10_PvsNp <- FindMarkers(TPs_merged_SCT, ident.1 = "P", ident.2 = "N", group.by = "prolif", subset.ident = "D10", assay = "SCT", min.pct =
0.05, only.pos = F)
D10_PvsNp$cluster <- "D10_PvsNp"
D10_PvsNp$gene <- rownames(D10_PvsNp)

D12_PvsNp <- FindMarkers(TPs_merged_SCT, ident.1 = "P", ident.2 = "N", group.by = "prolif", subset.ident = "D12", assay = "SCT", min.pct =
0.05, only.pos = F)
D12_PvsNp$cluster <- "D12_PvsNp"
D12_PvsNp$gene <- rownames(D12_PvsNp)

all_DEGs <- rbind(H12_PvsNp, D1_PvsNp, D2_PvsNp, D4_PvsNp, D6_PvsNp, D8_PvsNp, D10_PvsNp, D12_PvsNp)
all_DEGs$p_val_adj = p.adjust(all_DEGs$p_val, method='BH')
all_DEGs$delta <- abs(as.double(all_DEGs$pct.1)-as.double(all_DEGs$pct.2))
all_DEGs$minor <- ifelse(all_DEGs$pct.1 < 0.10 | all_DEGs$pct.2 < 0.10, 1, 0)
all_DEGs <- all_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
all_DEGs <- all_DEGs[all_DEGs$p_val_adj<0.05,]

write.table(all_DEGs, "/data/scRNA/IMR90_timecourse/DEGs_ProVSNonPro_allTPs.tsv", quote = F, row.names = F, sep = "\t")

}

87

Idents(sct_prolif) <- sct_prolif$timepoint

proSubset_H12vsPro <- FindMarkers(sct_prolif, ident.1 = "H12", ident.2 = "Pro", assay = "SCT", min.pct = 0.05, only.pos = F)
proSubset_H12vsPro$cluster <- "H12vsPro"
proSubset_H12vsPro$gene <- rownames(proSubset_H12vsPro)
proSubset_D1vsH12 <- FindMarkers(sct_prolif, ident.1 = "D1", ident.2 = "H12", assay = "SCT", min.pct = 0.05, only.pos = F)
proSubset_D1vsH12$cluster <- "D1vsH12"
proSubset_D1vsH12$gene <- rownames(proSubset_D1vsH12)
proSubset_D2vsD1 <- FindMarkers(sct_prolif, ident.1 = "D2", ident.2 = "D1", assay = "SCT", min.pct = 0.05, only.pos = F)
proSubset_D2vsD1$cluster <- "D2vsD1"
proSubset_D2vsD1$gene <- rownames(proSubset_D2vsD1)
proSubset_D4vsD2 <- FindMarkers(sct_prolif, ident.1 = "D4", ident.2 = "D2", assay = "SCT", min.pct = 0.05, only.pos = F)
proSubset_D4vsD2$cluster <- "D4vsD2"
proSubset_D4vsD2$gene <- rownames(proSubset_D4vsD2)
proSubset_D6vsD4 <- FindMarkers(sct_prolif, ident.1 = "D6", ident.2 = "D4", assay = "SCT", min.pct = 0.05, only.pos = F)
proSubset_D6vsD4$cluster <- "D6vsD4"
proSubset_D6vsD4$gene <- rownames(proSubset_D6vsD4)
proSubset_D8vsD6 <- FindMarkers(sct_prolif, ident.1 = "D8", ident.2 = "D6", assay = "SCT", min.pct = 0.05, only.pos = F)
proSubset_D8vsD6$cluster <- "D8vsD6"
proSubset_D8vsD6$gene <- rownames(proSubset_D8vsD6)
proSubset_D10vsD8 <- FindMarkers(sct_prolif, ident.1 = "D10", ident.2 = "D8", assay = "SCT", min.pct = 0.05, only.pos = F)
proSubset_D10vsD8$cluster <- "D10vsD8"
proSubset_D10vsD8$gene <- rownames(proSubset_D10vsD8)
proSubset_D12vsD10 <- FindMarkers(sct_prolif, ident.1 = "D12", ident.2 = "D10", assay = "SCT", min.pct = 0.05, only.pos = F)
proSubset_D12vsD10$cluster <- "D12vsD10"
proSubset_D12vsD10$gene <- rownames(proSubset_D12vsD10)
proSubset_D12vsPro <- FindMarkers(sct_prolif, ident.1 = "D12", ident.2 = "Pro", assay = "SCT", min.pct = 0.05, only.pos = F)
proSubset_D12vsPro$cluster <- "D12vsPro"
proSubset_D12vsPro$gene <- rownames(proSubset_D12vsPro)

all_DEGs <-
rbind(proSubset_H12vsPro,proSubset_D1vsH12,proSubset_D2vsD1,proSubset_D4vsD2,proSubset_D6vsD4,proSubset_D8vsD6,proSubset_D10v
sD8,proSubset_D12vsD10,proSubset_D12vsPro)
all_DEGs$p_val_adj = p.adjust(all_DEGs$p_val, method='BH')
all_DEGs$delta <- abs(as.double(all_DEGs$pct.1)-as.double(all_DEGs$pct.2))
all_DEGs$minor <- ifelse(all_DEGs$pct.1 < 0.10 | all_DEGs$pct.2 < 0.10, 1, 0)
all_DEGs <- all_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
all_DEGs <- all_DEGs[all_DEGs$p_val_adj<0.05,]

write.table(arrange(all_DEGs,desc(cluster),desc(avg_log2FC)),
"/data/scRNA/IMR90_timecourse/ProliferatingCells_merged_DEGs_sequentialTPs.tsv", quote = F, row.names = F, sep = "\t")

####################### proliferating vs quiescent vs d12 senescent
Pro <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_Pro.RDS.gz")
D12_1 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D12_1.RDS.gz")
D12_2 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_D12_2.RDS.gz")
Qui_1 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_Qui_1.RDS.gz")
Qui_2 <- readRDS("/data/scRNA/IMR90_timecourse/Singlets_initial_seuratObjects/Singlets_seuratObj_Qui_2.RDS.gz")

Pro$timepoint <- "Pro"
D12_1$timepoint <- "D12"
D12_2$timepoint <- "D12"
Qui_1$timepoint <- "Qui"
Qui_2$timepoint <- "Qui"

MainStates_merged <- merge(Pro, y = c(D12_1, D12_2, Qui_1, Qui_2), add.cell.ids = c("Pro", "D12_1", "D12_2", "Qui_1", "Qui_2"), merge.data
= T)

ordered_Timepoints <- c("Pro","D12","Qui")
MainStates_merged$timepoint <- factor(MainStates_merged$timepoint, levels= ordered_Timepoints)

MainStates_merged$library <- ''
88

MainStates_merged$library[grep(pattern="Pro_", colnames(MainStates_merged))] <- "Pro"
MainStates_merged$library[grep(pattern="D12_1", colnames(MainStates_merged))] <- "D12_1"
MainStates_merged$library[grep(pattern="D12_2", colnames(MainStates_merged))] <- "D12_2"
MainStates_merged$library[grep(pattern="Qui_1", colnames(MainStates_merged))] <- "Qui_1"
MainStates_merged$library[grep(pattern="Qui_2", colnames(MainStates_merged))] <- "Qui_2"
MainStates_merged[["percent.mito"]] <- PercentageFeatureSet(MainStates_merged, pattern = "^MT-")
MainStates_merged[["percent.ribo"]] <- PercentageFeatureSet(MainStates_merged, pattern = "^RP[SL]")
MainStates_merged[["percent.MALAT1"]] <- PercentageFeatureSet(MainStates_merged, pattern = "MALAT1")
MainStates_merged[["percent.NEAT1"]] <- PercentageFeatureSet(MainStates_merged, pattern = "NEAT1")

MainStates_subset <- subset(MainStates_merged, percent.MALAT1 <= 20 & percent.NEAT1 <= 20 &
nFeature_RNA > 500 & nCount_RNA > 2000 & percent.mito < 25 & percent.ribo < 25)

D12_proliferating <- WhichCells(subset(MainStates_subset, subset = timepoint == "D12" & ((MKI67 > 0 & CCNA2 > 0) | (MKI67 > 0 & TOP2A > 0)
| (TOP2A > 0 & CCNA2 > 0) | (MCM2 > 0 & MKI67 > 0) | (MCM2 > 0 & TOP2A > 0) | (MCM2 > 0 & CCNA2 > 0) |
(PCNA > 0 & MKI67 > 0) | (PCNA > 0 & TOP2A > 0) | (PCNA > 0 & CCNA2 > 0) | (PCNA > 0 & MCM2 > 0))))

##### rename D12 proliferating and D12 senescent
Idents(MainStates_subset) <- "timepoint"
MainStates_subset <- SetIdent(MainStates_subset, cells = D12_proliferating, value = "D12 proliferating")
MainStates_subset$timepoint <- Idents(MainStates_subset)

D12_senescent <- WhichCells(subset(MainStates_subset, subset = timepoint == "D12"))
MainStates_subset <- SetIdent(MainStates_subset, cells = D12_senescent, value = "D12 senescent")
#make sure the new identity name gets assigned to the timepoint metadata
MainStates_subset$timepoint <- Idents(MainStates_subset)

saveRDS.compress.in.BG(MainStates_merged, fname = "/data/scRNA/IMR90_timecourse/MainStates_merged_12-20-22.RDS.gz")

MainStates_subset <- SCTransform(MainStates_subset, vst.flavor = "v2", ncells = 15000, variable.features.n = 4000,
vars.to.regress = c("percent.mito", "percent.ribo","percent.MALAT1","percent.NEAT1"),
return.only.var.genes = T, verbose = T)

#MainStates_subset$timepoint[grep(pattern="Qui_", colnames(MainStates_subset))] <- "Qui"
MainStates_subset <- RunPCA(MainStates_subset, npcs = 100, verbose = F, ndims.print = F)
ElbowPlot(MainStates_subset, ndims = 50)
MainStates_subset <- RunUMAP(MainStates_subset, min.dist = 0.25, dims = 1:30, verbose = T)
MainStates_subset <- FindNeighbors(MainStates_subset, reduction = "umap", dims = 1:2, verbose = T)
MainStates_subset <- FindClusters(MainStates_subset, resolution = 0.125, verbose = T)

ordered_Timepoints <- c("Pro","D12 senescent", "D12_proliferating","Qui")
MainStates_merged$timepoint <- factor(MainStates_merged$timepoint, levels= ordered_Timepoints)

DimPlot(MainStates_subset, group.by = "timepoint")

main_title_size = 18
axis_text_size = 9
axis_title_size = 9

MainStates_timepoint_plotLegend <- DimPlot(MainStates_subset, group.by = "timepoint") +
ggtitle('') +
theme(text = element_text(size = main_title_size, face = "plain"), axis.title=element_text(size=axis_title_size), axis.text = element_text(size =
axis_text_size)) +
theme(legend.position = "bottom", legend.direction = "horizontal")
#theme(legend.position = "bottom", legend.direction = "horizontal", axis.text.x=element_text(size=14), axis.title.x = element_text(size=14),
axis.title.y = element_text(size=14))

MainStates_cluster_plot <- DimPlot(MainStates_subset, group.by = "SCT_snn_res.0.15", label = T, label.box = T) +
89

ggtitle('') +
NoLegend() +
theme(text = element_text(size = main_title_size, face = "plain"), axis.title=element_text(size=axis_title_size), axis.text = element_text(size =
axis_text_size))
(MainStates_timepoint_plotLegend / MainStates_cluster_plot)

saveRDS.compress.in.BG(MainStates_subset, fname = "/data/scRNA/IMR90_timecourse/MainStates_merged_1-1-22.RDS.gz")

MainStates_subset <- readRDS("/data/scRNA/IMR90_timecourse/MainStates_merged_12-12-22.RDS.gz")

################# Find DEGs
Idents(MainStates_subset)
MainStates_ProMarkers <- FindMarkers(MainStates_subset, ident.1 = "Pro", ident.2 = c("D12 senescent","Qui"), assay = "SCT", min.pct = 0.05,
only.pos = F)
MainStates_ProMarkers$cluster <- "Pro"
MainStates_ProMarkers$gene <- rownames(MainStates_ProMarkers)

#MainStates_ProMarkers2 <- FindMarkers(MainStates_subset, ident.1 = "Pro", ident.2 = c("D12","Qui","D12_proliferating"), assay = "SCT",
min.pct = 0.05, only.pos = F)
#MainStates_ProMarkers2$cluster <- "Pro2"
#MainStates_ProMarkers2$gene <- rownames(MainStates_ProMarkers2)

MainStates_D12Markers <- FindMarkers(MainStates_subset, ident.1 = "D12 senescent", ident.2 = c("Pro","Qui"), assay = "SCT", min.pct = 0.05,
only.pos = F)
MainStates_D12Markers$cluster <- "D12 senescent"
MainStates_D12Markers$gene <- rownames(MainStates_D12Markers)

MainStates_QuiMarkers <- FindMarkers(MainStates_subset, ident.1 = "Qui", ident.2 = c("Pro","D12 senescent"), assay = "SCT", min.pct = 0.05,
only.pos = F)
MainStates_QuiMarkers$cluster <- "Qui"
MainStates_QuiMarkers$gene <- rownames(MainStates_QuiMarkers)

#MainStates_D12proMarkers <- FindMarkers(MainStates_subset, ident.1 = "D12 proliferating", ident.2 = c("D12 senescent","Qui"), assay =
"SCT", min.pct = 0.05, only.pos = F)
#MainStates_D12proMarkers$cluster <- "D12pro"
#MainStates_D12proMarkers$gene <- rownames(MainStates_D12proMarkers)

#MainStates_all_DEGs <- rbind(MainStates_ProMarkers, MainStates_D12Markers, MainStates_QuiMarkers, MainStates_D12proMarkers)
MainStates_all_DEGs <- rbind(MainStates_ProMarkers, MainStates_D12Markers, MainStates_QuiMarkers)
MainStates_all_DEGs$p_val_adj = p.adjust(MainStates_all_DEGs$p_val, method='BH')
MainStates_all_DEGs$delta <- abs(as.double(MainStates_all_DEGs$pct.1)-as.double(MainStates_all_DEGs$pct.2))
MainStates_all_DEGs$minor <- ifelse(MainStates_all_DEGs$pct.1 < 0.10 | MainStates_all_DEGs$pct.2 < 0.10, 1, 0)
MainStates_all_DEGs <- MainStates_all_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
MainStates_all_DEGs <- MainStates_all_DEGs[MainStates_all_DEGs$p_val_adj<0.05,]
MainStates_all_DEGs$avg_log2FC <- round(MainStates_all_DEGs$avg_log2FC, 2)
#MainStates_all_DEGs$p_val_adj <- round(MainStates_all_DEGs$p_val_adj, 2)
MainStates_all_DEGs$pct.1 <- round(MainStates_all_DEGs$pct.1, 2)
MainStates_all_DEGs$pct.2 <- round(MainStates_all_DEGs$pct.2, 2)
MainStates_all_DEGs$delta <- round(MainStates_all_DEGs$delta, 2)

write.table(arrange(MainStates_all_DEGs,desc(cluster),desc(avg_log2FC)), "/data/scRNA/IMR90_timecourse/MainStates_DEGs.tsv", quote = F,
row.names = F, sep = "\t")

#old figures
{
#SCpubr::do_DotPlot(MainStates_subset, features = c("TXN","ACAT2","RGCC","IDI1", "MSMO1",
# "NEAT1","COL6A1","MEG3","PXDN","CCDC80",
# "MMP1","UBE2S", "TUBB4B","CENPF","PTTG1"),
# group.by = "timepoint", dot.scale = 8, flip = T, rotate_x_labels = T, assay = "SCT", ylab = "Timepoint",plot.title = "Top Markers by
Log2 Fold Change")
#
#MainStates_EnrichmentHeatmap_plot <- SCpubr::do_EnrichmentHeatmap(MainStates_subset, list_genes =
c("SH3PXD2A","PIEZO2","CADM1","CCND2","SEMA5A",
90

# "TMEM243","HOXB2","UNC50","SDHD","MFAP4",
# "DDX39A","CDK2AP1","TUBG1","CCND3","GALE"),
# group.by = "timepoint", split.horizontal = F, column_names_rot = 45)
}

inverse_D12proliferating <- WhichCells(subset(MainStates_subset, subset = timepoint == "D12 senescent" | timepoint == "Qui" | timepoint ==
"Pro"))

subset(MainStates_subset, subset = timepoint == "Pro")

#ordered_Timepoints <- c("D12","Qui","Pro")
#MainStates_subset$timepoint <- factor(MainStates_subset$timepoint, levels= ordered_Timepoints)

DoHeatmap(MainStates_subset, cells = inverse_D12proliferating,
features = c("COL6A1","PXDN","CCDC80", "COL6A2", "MYLK", "SPARC", "FOSB","MFGE8","DST","COL5A1", #D12 senescent
"MMP1","UBE2S", "CENPF","TUBB4B","PTTG1", "HIST1H4C","CKS2","TUBA1C","CCNB1", "UBE2C", #Proliferating
"TXN","ACAT2","RGCC","IDI1", "MSMO1","CTHRC1","FDPS","FABP3","GLRX","EEF1A1"), #Quiescent
group.by = "timepoint", assay = "SCT",
angle = 0, group.bar = F, draw.lines = T, lines.width = 100) + #scale_fill_viridis(option = "B")
scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(n = 11, name = "RdBu"))) +
guides(color=FALSE) +
guides(fill=guide_legend(title="Scaled Expression")) +
theme(axis.text = element_text(face = "bold", colour = "black"))

Mainstates_DEGbyTP_densityPlot <- Nebulosa::plot_density(subset(MainStates_subset, subset = timepoint == "D12 senescent" | timepoint ==
"Qui" | timepoint == "Pro"),
features = c("COL6A1","PXDN","CCDC80", "COL6A2", "MYLK",
"MMP1","UBE2S", "CENPF","TUBB4B","PTTG1",
"TXN","ACAT2","RGCC","IDI1", "MSMO1"), size = 0.01, combine = T) &
NoLegend() &
theme(plot.title = element_text(size = main_title_size, face = "plain", hjust = 0.5), axis.title=element_text(size=axis_title_size),
axis.text = element_text(size = axis_text_size))

Mainstates_DEGbyTP_densityPlot + plot_layout(ncol = 5)

DoHeatmap(MainStates_subset, cells = inverse_D12proliferating, group.colors = c("white","white","white","white"),
features = c("SH3PXD2A","CADM1","CCND2", "NFIX", "PCSK7", "COL18A1", "PIEZO2","MRVI1","SEMA5A","AP2B1", #D12 senescent
"DDX39A","CDK2AP1", "TUBG1","FOSL1","CCND3", "GALE","CENPW","GGH","FKBP4", "PSMD9", #Proliferating
"TMEM243","HOXB2","UNC50","SNHG25", "MFAP4","SDHD","PPIL3","FABP3","CDKN2A","SVBP"), #Quiescent
group.by = "timepoint", assay = "SCT",
angle = 0, group.bar = T, draw.lines = T, lines.width = 100, ) + #scale_fill_viridis(option = "B")
scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(n = 11, name = "RdBu"))) +
guides(color=FALSE) +
guides(fill=guide_legend(title="Scaled Expression")) +
theme(axis.text = element_text(face = "plain", colour = "black"))

DefaultAssay(MainStates_subset)

MainStates_featurePlot <- FeaturePlot(MainStates_subset, cols = c("grey","darkblue"), pt.size = 0.01, keep.scale = "feature",
features = c("MKI67", "CDKN1A","CDKN2A","CDKN1B","CDKN2B","CDKN1C","SH3PXD2A","DDX39A","TMEM243"), label =
F, combine = T) &
NoLegend() &
theme(text = element_text(size = main_title_size, face = "plain"), axis.title=element_text(size=axis_title_size), axis.text =
element_text(size = axis_text_size))

MainStates_featurePlot <- FeaturePlot(MainStates_subset, cols = c("grey","darkblue"), pt.size = 0.01, keep.scale = "feature",
features = c("MKI67", "CDKN1A","CDKN2A","CDKN1B","CDKN2B","CDKN1C"), label = F, combine = F) +
NoLegend()
#scale_colour_gradientn(colours = rev(brewer.pal(n = 5, name = "RdBu")))
#scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(n = 11, name = "RdBu")) )

91

(MainStates_featurePlot[[1]] | MainStates_featurePlot[[2]] | MainStates_featurePlot[[3]]) /
(MainStates_featurePlot[[4]] | MainStates_featurePlot[[5]] | MainStates_featurePlot[[6]])

Idents(MainStates_subset) <- MainStates_subset$seurat_clusters
MainStates_allMarkers <- FindAllMarkers(MainStates_subset, assay = "SCT", min.pct = 0.05)
MainStates_allMarkers$p_val_adj = p.adjust(MainStates_allMarkers$p_val, method='BH')
MainStates_allMarkers$delta <- abs(as.double(MainStates_allMarkers$pct.1)-as.double(MainStates_allMarkers$pct.2))
MainStates_allMarkers$minor <- ifelse(MainStates_allMarkers$pct.1 < 0.10 | MainStates_allMarkers$pct.2 < 0.10, 1, 0)
MainStates_allMarkers <- MainStates_allMarkers[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
MainStates_allMarkers <- MainStates_allMarkers[MainStates_allMarkers$p_val_adj<0.05,]
MainStates_allMarkers$avg_log2FC <- round(MainStates_allMarkers$avg_log2FC, 2)
MainStates_allMarkers$p_val_adj <- round(MainStates_allMarkers$p_val_adj, 2)
MainStates_allMarkers$pct.1 <- round(MainStates_allMarkers$pct.1, 2)
MainStates_allMarkers$pct.2 <- round(MainStates_allMarkers$pct.2, 2)
MainStates_allMarkers$delta <- round(MainStates_allMarkers$delta, 2)

write.table(arrange(MainStates_allMarkers,desc(cluster),desc(avg_log2FC)), "/data/scRNA/IMR90_timecourse/MainStates_DEGsClusters.tsv",
quote = F, row.names = F, sep = "\t")

Idents(MainStates_subset) <- MainStates_subset$seurat_clusters
MainStates_pro <- subset(MainStates_subset, timepoint == "Pro")
table(Idents(MainStates_pro))
MainStates_allMarkers_Pro <- FindAllMarkers(MainStates_pro, assay = "SCT", min.pct = 0.05, verbose = T, min.cells.group = 50)
MainStates_allMarkers_Pro$p_val_adj = round(p.adjust(MainStates_allMarkers_Pro$p_val, method='BH'), 2)
MainStates_allMarkers_Pro$delta <- round(abs(as.double(MainStates_allMarkers_Pro$pct.1)-as.double(MainStates_allMarkers_Pro$pct.2)), 2)
MainStates_allMarkers_Pro$minor <- ifelse(MainStates_allMarkers_Pro$pct.1 < 0.10 | MainStates_allMarkers_Pro$pct.2 < 0.10, 1, 0)
MainStates_allMarkers_Pro <- MainStates_allMarkers_Pro[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
MainStates_allMarkers_Pro <- MainStates_allMarkers_Pro[MainStates_allMarkers_Pro$p_val_adj<0.05,]
MainStates_allMarkers_Pro$avg_log2FC <- round(MainStates_allMarkers_Pro$avg_log2FC, 2)
MainStates_allMarkers_Pro$pct.1 <- round(MainStates_allMarkers_Pro$pct.1, 2)
MainStates_allMarkers_Pro$pct.2 <- round(MainStates_allMarkers_Pro$pct.2, 2)

write.table(arrange(MainStates_allMarkers_Pro,desc(cluster),desc(avg_log2FC)),
"/data/scRNA/IMR90_timecourse/MainStates_allClusterMarkers_withinPro.tsv", quote = F, row.names = F, sep = "\t")

Idents(MainStates_subset) <- MainStates_subset$seurat_clusters
MainStates_qui <- subset(MainStates_subset, timepoint == "Qui")
table(Idents(MainStates_qui))
MainStates_allMarkers_Qui <- FindAllMarkers(MainStates_qui, assay = "SCT", min.pct = 0.05, verbose = T, min.cells.group = 50)
MainStates_allMarkers_Qui$p_val_adj = round(p.adjust(MainStates_allMarkers_Qui$p_val, method='BH'), 2)
MainStates_allMarkers_Qui$delta <- round(abs(as.double(MainStates_allMarkers_Qui$pct.1)-as.double(MainStates_allMarkers_Qui$pct.2)), 2)
MainStates_allMarkers_Qui$minor <- ifelse(MainStates_allMarkers_Qui$pct.1 < 0.10 | MainStates_allMarkers_Qui$pct.2 < 0.10, 1, 0)
MainStates_allMarkers_Qui <- MainStates_allMarkers_Qui[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
MainStates_allMarkers_Qui <- MainStates_allMarkers_Qui[MainStates_allMarkers_Qui$p_val_adj<0.05,]
MainStates_allMarkers_Qui$avg_log2FC <- round(MainStates_allMarkers_Qui$avg_log2FC, 2)
MainStates_allMarkers_Qui$pct.1 <- round(MainStates_allMarkers_Qui$pct.1, 2)
MainStates_allMarkers_Qui$pct.2 <- round(MainStates_allMarkers_Qui$pct.2, 2)

write.table(arrange(MainStates_allMarkers_Qui,desc(cluster),desc(avg_log2FC)),
"/data/scRNA/IMR90_timecourse/MainStates_allClusterMarkers_withinQui.tsv", quote = F, row.names = F, sep = "\t")

Idents(MainStates_subset) <- MainStates_subset$seurat_clusters
MainStates_D12sen <- subset(MainStates_subset, timepoint == "D12 senescent")
table(Idents(MainStates_D12sen))
MainStates_allMarkers_D12sen <- FindAllMarkers(MainStates_D12sen, assay = "SCT", min.pct = 0.05, verbose = T, min.cells.group = 50)
MainStates_allMarkers_D12sen$p_val_adj = round(p.adjust(MainStates_allMarkers_D12sen$p_val, method='BH'), 2)
MainStates_allMarkers_D12sen$delta <- round(abs(as.double(MainStates_allMarkers_D12sen$pct.1)-
as.double(MainStates_allMarkers_D12sen$pct.2)), 2)
MainStates_allMarkers_D12sen$minor <- ifelse(MainStates_allMarkers_D12sen$pct.1 < 0.10 | MainStates_allMarkers_D12sen$pct.2 < 0.10, 1,
0)
MainStates_allMarkers_D12sen <- MainStates_allMarkers_D12sen[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta",
"minor")]
92

MainStates_allMarkers_D12sen <- MainStates_allMarkers_D12sen[MainStates_allMarkers_D12sen$p_val_adj<0.05,]
MainStates_allMarkers_D12sen$avg_log2FC <- round(MainStates_allMarkers_D12sen$avg_log2FC, 2)
MainStates_allMarkers_D12sen$pct.1 <- round(MainStates_allMarkers_D12sen$pct.1, 2)
MainStates_allMarkers_D12sen$pct.2 <- round(MainStates_allMarkers_D12sen$pct.2, 2)

write.table(arrange(MainStates_allMarkers_D12sen,desc(cluster),desc(avg_log2FC)),
"/data/scRNA/IMR90_timecourse/MainStates_allClusterMarkers_withinD12sen.tsv", quote = F, row.names = F, sep = "\t")

APPENDIX B
Seurat_to_Slingshot.R – R code used to create slingshot objects, calculate pseudotime, lineages, and
curves.
.libPaths("/opt/R/4.0.4/lib/R/library")
setwd("/data/scRNA/IMR90_timecourse/")
#Sys.unsetenv("GITHUB_PAT")

library(SingleCellExperiment)
library(Tempora)
library(Seurat)
library(Seurat.utils)
library(SeuratWrappers)
library(Matrix)
library(dyndimred)
library(Seurat.utils)
library(diffusionMap)
library(slingshot)
BiocParallel::register(BiocParallel::SerialParam())
library(tradeSeq)
library(paletteer)
library(hsslda)
library(khroma)
library(palettesForR)
library(pals)
library(randomcoloR)
library(RColorBrewer)
library(igraph)
library(psupertime)
library(scanalysis)
library(dplyr)
library(tidyr)
library(tibble)
library(patchwork)
library(grDevices)

# the ':::' operator allows access to unexported functions!!

TPs_merged_SCT <- readRDS("/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/TPs_merged_SCT_final.RDS.gz")

DimPlot(TPs_merged_SCT, group.by = "timepoint", label = T)

TPs_merged_SCT$orig.ident <- NULL
TPs_merged_SCT$percent.MALAT1 <- NULL
TPs_merged_SCT$percent.NEAT1 <- NULL
TPs_merged_SCT$seurat_clusters <- NULL

#DefaultAssay(TPs_merged_SCT) <- 'RNA'
#TPs_merged_SCT_diet <- DietSeurat(TPs_merged_SCT, assays = "RNA")
93

list_TPs_merged_SCT <- SplitObject(TPs_merged_SCT, split.by = "timepoint")

#RPCA

# Pro-H12-D1-D2 list
#list_for_rpca <- list(list_TPs_merged_SCT[[1]], list_TPs_merged_SCT[[2]], list_TPs_merged_SCT[[3]],list_TPs_merged_SCT[[4]])

# D2-D4-D6-D8
#list_for_rpca <- list(list_TPs_merged_SCT[[5]], list_TPs_merged_SCT[[6]], list_TPs_merged_SCT[[7]],list_TPs_merged_SCT[[8]])

# D8-D10-D12
list_for_rpca <- list(list_TPs_merged_SCT[[7]], list_TPs_merged_SCT[[8]], list_TPs_merged_SCT[[9]])

list_for_rpca <- lapply(X = list_for_rpca, FUN = SCTransform, method = "glmGamPoi")
features <- SelectIntegrationFeatures(object.list = list_for_rpca, nfeatures = 5000)
list_for_rpca <- PrepSCTIntegration(object.list = list_for_rpca, anchor.features = features)
list_for_rpca <- lapply(X = list_for_rpca, FUN = RunPCA, features = features)

#k30_anchors <- FindIntegrationAnchors(object.list = list_for_rpca, normalization.method = "SCT", anchor.features = features, dims = 1:50,
reduction = "rpca", k.anchor = 30)
k10_anchors <- FindIntegrationAnchors(object.list = list_for_rpca, normalization.method = "SCT", anchor.features = features, dims = 1:50,
reduction = "rpca", k.anchor = 10)
#k5_anchors <- FindIntegrationAnchors(object.list = list_for_rpca, normalization.method = "SCT", anchor.features = features, dims = 1:50,
reduction = "rpca", k.anchor = 5)

#k30_rpca_sct <- IntegrateData(anchorset = k30_anchors, normalization.method = "SCT", dims = 1:30, verbose = T)
k10_rpca_int <- IntegrateData(anchorset = k10_anchors, normalization.method = "SCT", dims = 1:75, verbose = T)
#k10_rpca_int <- IntegrateData(anchorset = k10_anchors, normalization.method = "SCT", dims = 1:100, verbose = T)

#k5_rpca_sct <- IntegrateData(anchorset = k5_anchors, normalization.method = "SCT", dims = 1:50, verbose = T)

#k30_rpca_sct <- RunPCA(k30_rpca_sct, verbose = T, npcs = 100)
k10_rpca_int <- RunPCA(k10_rpca_int, verbose = T, npcs = 100)
#k5_rpca_sct <- RunPCA(k5_rpca_sct, verbose = T, npcs = 100)

#k30_rpca_sct <- RunUMAP(k30_rpca_sct, reduction = "pca", dims = 1:20, verbose = T)
#k10_rpca_int <- RunUMAP(k10_rpca_int, repulsion.strength = 10, min.dist = 1, reduction = "pca", dims = 1:50, verbose = T)
k10_rpca_int <- RunUMAP(k10_rpca_int, repulsion.strength = 15, min.dist = 0.4, reduction = "pca", dims = 1:100, verbose = T)
k10_rpca_int <- RunUMAP(k10_rpca_int, repulsion.strength = 5, min.dist = 1, reduction = "pca", dims = 1:75, verbose = T)
#k5_rpca_sct <- RunUMAP(k5_rpca_sct, repulsion.strength = 5, min.dist = 1, reduction = "pca", dims = 1:75, verbose = T)
#DefaultAssay(k10_rpca_int) <- "SCT"

DimPlot(k10_rpca_int, reduction = "umap", group.by = "timepoint")
DimPlot(k10_rpca_int, reduction = "umap", group.by = "timepoint") | FeaturePlot(k10_rpca_int, features = "MKI67")
DimPlot(k5_rpca_sct, reduction = "umap", group.by = "timepoint") | FeaturePlot(k5_rpca_sct, features = "MKI67")
#saveRDS.compress.in.BG(k10_rpca_int, fname = "Pro-H12-D1-D2_k10_rpca_int.RDS.gz")
k10_rpca_int2 <- readRDS("/data/scRNA/IMR90_timecourse/Pro-H12-D1-D2_k10_rpca_int.RDS.gz")

#k10_rpca_int <- IntegrateData(anchorset = k10_anchors, normalization.method = "SCT", dims = 1:30, verbose = T)
#k10_rpca_int <- RunPCA(k10_rpca_int, verbose = T, npcs = 100)
k10_rpca_int <- RunUMAP(k10_rpca_int, repulsion.strength = 5, min.dist = 0.5, reduction = "pca", dims = 1:75, verbose = T)
DimPlot(k10_rpca_int, reduction = "umap", group.by = "timepoint")
saveRDS.compress.in.BG(k10_rpca_int, fname = "Pro-H12-D1-D2_k10_rpca_int.RDS.gz")

#k10_rpca_int <- IntegrateData(anchorset = k10_anchors, normalization.method = "SCT", dims = 1:50, verbose = T)
#k10_rpca_int <- RunPCA(k10_rpca_int, verbose = T, npcs = 100)
#k10_rpca_int <- RunUMAP(k10_rpca_int, repulsion.strength = 5, min.dist = 1, reduction = "pca", dims = 1:75, verbose = T)
saveRDS.compress.in.BG(k10_rpca_int, fname = "D2-D4-D6-D8_k10_rpca_int.RDS.gz")

#k10_rpca_int <- IntegrateData(anchorset = k10_anchors, normalization.method = "SCT", dims = 1:75, verbose = T)
#k10_rpca_int <- RunPCA(k10_rpca_int, verbose = T, npcs = 100)
#k10_rpca_int <- RunUMAP(k10_rpca_int, repulsion.strength = 10, min.dist = 0.2, reduction = "pca", dims = 1:50, verbose = T)
seu_D8_to_D12 <- RunUMAP(k10_rpca_int, repulsion.strength = 7.5, min.dist = 0.2, reduction = "pca", dims = 1:50, verbose = T)
saveRDS.compress.in.BG(k10_rpca_int, fname = "/data/scRNA/IMR90_timecourse/D8-D10-D12_k10_rpca_int.RDS.gz")

Idents(seu_D8_to_D12) <- seu_D8_to_D12$timepoint
94

ordered_Timepoints <- c("D8", "D10", "D12")
seu_D8_to_D12$timepoint <- factor(seu_D8_to_D12$timepoint, levels= ordered_Timepoints)
seu_D8_to_D12 <- RunUMAP(seu_D8_to_D12, repulsion.strength = 7.5, min.dist = 0.2, reduction = "pca", dims = 1:50, verbose = T)
seu_D8_to_D12 <- FindNeighbors(seu_D8_to_D12, k.param = 50, reduction = "umap", dims = 1:2, verbose = T)
seu_D8_to_D12 <- FindClusters(seu_D8_to_D12, resolution = 0.15, verbose = T)
DimPlot(seu_D8_to_D12, group.by = "timepoint", reduction = "umap",pt.size = 0.01, label = F)
saveRDS.compress.in.BG(seu_D8_to_D12, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_forLateTrajectoryPlotting.RDS.gz")

##########

seu_Pro_to_D2 <- readRDS("/data/scRNA/IMR90_timecourse/Pro-H12-D1-D2_k10_rpca_int.RDS.gz")
seu_D2_to_D8 <- readRDS("/data/scRNA/IMR90_timecourse/D2-D4-D6-D8_k10_rpca_int.RDS.gz")
seu_D8_to_D12 <- readRDS("/data/scRNA/IMR90_timecourse/D8-D10-D12_k10_rpca_int.RDS.gz")

#used to optimize on 10/13/22
seu_D8_to_D12 <- RunUMAP(seu_D8_to_D12, repulsion.strength = 7.5, min.dist = 0.2, reduction = "pca", dims = 1:50, verbose = T)
DimPlot(seu_D8_to_D12, group.by = "timepoint", reduction = "umap",pt.size = 0.1, label = F)

Idents(this_subset) <- this_subset$timepoint
ordered_Timepoints <- c("Pro","H12","D1","D2","D4","D6","D8","D10","D12")
ordered_Timepoints <- c("Pro","H12","D1","D2")
ordered_Timepoints <- c("D2", "D4", "D6", "D8")
ordered_Timepoints <- c("D8", "D10", "D12")

this_subset$timepoint <- factor(this_subset$timepoint, levels= ordered_Timepoints)

DefaultAssay(this_subset) <- "integrated"
set.seed(1)

this_subset <- FindNeighbors(this_subset, k.param = 25, reduction = "umap", dims = 1:2, verbose = T)
this_subset <- FindNeighbors(this_subset, k.param = 50, reduction = "umap", dims = 1:2, verbose = T)

# Use 0.2 for pro to d2
this_subset <- FindClusters(this_subset, resolution = 0.2, verbose = T)
# Use 0.3 for d2 to d8
this_subset <- FindClusters(this_subset, resolution = 0.25, verbose = T)
# Use 0.15 for d8 to d12
this_subset <- FindClusters(this_subset, resolution = 0.15, verbose = T)

DimPlot(this_subset, group.by = "timepoint", reduction = "umap",pt.size = 0.1, label = F)
#DimPlot(this_subset, group.by = "integrated_snn_res.0.15", reduction = "umap",pt.size = 0.1, label = TRUE, repel = TRUE)

# use for pro to d2
#DimPlot(this_subset, group.by = "timepoint", reduction = "umap",pt.size = 0.1, label = F) |
#DimPlot(this_subset, group.by = "integrated_snn_res.0.2", reduction = "umap",pt.size = 0.1, label = TRUE, repel = TRUE)

# use for d2 to d8
#DimPlot(this_subset, group.by = "timepoint", reduction = "umap",pt.size = 0.1, label = F) |
#DimPlot(this_subset, group.by = "integrated_snn_res.0.25", reduction = "umap",pt.size = 0.1, label = TRUE, repel = TRUE)

#DimPlot(this_subset, split.by = "timepoint", reduction = "umap",pt.size = 0.1, label = F)

# Use 0.2 for pro to d2
#clustering <- this_subset@meta.data$integrated_snn_res.0.2
# Use 0.3 for d2 to d8
#clustering <- this_subset@meta.data$integrated_snn_res.0.25

dimred <- this_subset@reductions$umap@cell.embeddings

set.seed(1)
# use for pro to d2
#lineages <- as.SlingshotDataSet(getLineages(data = dimred, clusterLabels = clustering, start.clus = '18'))
#lineages

# use for d2 to d8
lineages <- as.SlingshotDataSet(getLineages(data = dimred, clusterLabels = clustering, start.clus = '17'))
95

lineages

curves <- as.SlingshotDataSet(getCurves(lineages, approx_points = 300, thresh = 100, stretch = 0.8, allow.breaks = T, shrink = 1))
curves

#saveRDS.compress.in.BG(this_subset, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_this_subset.RDS.gz")
#saveRDS.compress.in.BG(dimred, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_dimred.RDS.gz")
#saveRDS.compress.in.BG(clustering, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_clustering.RDS.gz")
#saveRDS.compress.in.BG(curves, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_curves.RDS.gz")
#saveRDS.compress.in.BG(lineages, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_lineages.RDS.gz")

# OLD this_subset <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_this_subset.RDS.gz")
#this_subset$droplet_type <- this_subset_remerged$droplet_type
#this_subset
#table(this_subset$droplet_type)
#table(this_subset$timepoint)
#this_subset_singlets <- subset(this_subset, subset = droplet_type == "Singlet")
#this_subset_singlets
saveRDS.compress.in.BG(this_subset, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_allCells.RDS.gz")
saveRDS.compress.in.BG(this_subset_singlets, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_singletsOnly.RDS.gz")

#dimred <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_dimred.RDS.gz")
#clustering <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_clustering.RDS.gz")
#curves <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_curves.RDS.gz")
#lineages <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_lineages.RDS.gz")

# OLD this_subset <- readRDS("/data/scRNA/IMR90_timecourse/seu_D2_to_D8_this_subset.RDS.gz")
#this_subset$droplet_type <- this_subset_remerged$droplet_type
#this_subset_singlets <- subset(this_subset, subset = droplet_type == "Singlet")
#table(this_subset$droplet_type)
#table(this_subset_singlets$droplet_type)
saveRDS.compress.in.BG(this_subset, fname = "/data/scRNA/IMR90_timecourse/seu_D2_to_D8_allCells.RDS.gz")
saveRDS.compress.in.BG(this_subset_singlets, fname = "/data/scRNA/IMR90_timecourse/seu_D2_to_D8_singletsOnly.RDS.gz")

#dimred <- readRDS("/data/scRNA/IMR90_timecourse/seu_D2_to_D8_dimred.RDS.gz")
#clustering <- readRDS("/data/scRNA/IMR90_timecourse/seu_D2_to_D8_clustering.RDS.gz")
#curves <- readRDS("/data/scRNA/IMR90_timecourse/seu_D2_to_D8_curves.RDS.gz")
#lineages <- readRDS("/data/scRNA/IMR90_timecourse/seu_D2_to_D8_lineages.RDS.gz")

#saveRDS.compress.in.BG(TPs_merged_SCT, fname =
"/data/scRNA/IMR90_timecourse/intermediate_seuratObjects/TPs_merged_SCT_final.RDS")
#this_subset <- AddMetaData(this_subset, metadata = as.data.frame(TPs_merged_SCT$droplet_type), col.name =
"droplet_type")#TPs_merged_SCT$droplet_type, $droplet_type

#these two are D8-D12-subset1 and D8-D12-subset2:
this_subset <- subset1
this_subset <- subset2

table(this_subset$timepoint)
colnames(this_subset)
this_subset_split <- SplitObject(this_subset, split.by = "library")

length(this_subset_split)

for(x in seq_along(along.with = this_subset_split)){

#x=1
print(x)
print(head(colnames(this_subset_split[[x]])))
table(this_subset_split[[x]]$timepoint)
this_subset_split[[x]]$editCellName <- colnames(this_subset_split[[x]])
96

head(this_subset_split[[x]]$editCellName)
this_subset_split[[x]]$editCellName <- substr(this_subset_split[[x]]$editCellName,start = 1,stop = 18) #18 = length of wrong cell name
(AAACCCAAGTCTTCCC-1_6) without the "_6" ending
head(this_subset_split[[x]]$editCellName)
this_subset_split[[x]] <- RenameCells(this_subset_split[[x]], new.names = this_subset_split[[x]]$editCellName)
this_subset_split[[x]] <- RenameCells(this_subset_split[[x]], add.cell.id = unique(this_subset_split[[x]]$library))
print(head(colnames(this_subset_split[[x]])))
this_subset_split[[x]]

}

table(merge_for_droplettype$timepoint)
#colnames(merge_for_droplettype)

droplet_key <- merge_for_droplettype$droplet_type
this_subset_remerged <- AddMetaData(this_subset_remerged, metadata = droplet_key, col.name = "droplet_type")
this_subset_remerged$droplet_type <- replace_na(this_subset_remerged$droplet_type, "Doublet")
table(this_subset_remerged$droplet_type)

this_subset$droplet_type <- this_subset_remerged$droplet_type
this_subset_singlets <- subset(this_subset, subset = droplet_type == "Singlet")
table(this_subset$droplet_type)
table(this_subset_singlets$droplet_type)

saveRDS.compress.in.BG(this_subset, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s1_allCells.RDS.gz")
saveRDS.compress.in.BG(this_subset_singlets, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s1_singletsOnly.RDS.gz")
saveRDS.compress.in.BG(this_subset, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s2_allCells.RDS.gz")
saveRDS.compress.in.BG(this_subset_singlets, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s2_singletsOnly.RDS.gz")

#Plot the lineages
par(mfrow=c(1,1))
par(bg = "white")

#saveRDS.compress.in.BG(pal1, fname = "/data/scRNA/IMR90_timecourse/pal1.RDS")
#saveRDS.compress.in.BG(pal2, fname = "/data/scRNA/IMR90_timecourse/pal2.RDS")
#saveRDS.compress.in.BG(pal3, fname = "/data/scRNA/IMR90_timecourse/pal3.RDS")
#saveRDS.compress.in.BG(pal4, fname = "/data/scRNA/IMR90_timecourse/pal4.RDS")

#pal1 <- distinctColorPalette(k = 40, altCol = F, runTsne = FALSE)
pal1 <- paletteer_d(`"palettesForR::Bears"`, n = 45, direction = -1)

plot(dimred, col = pal1[clustering], cex = 0.5, pch = 16)
lines(lineages, lwd = 1, col = 'black')
for(i in levels(clustering)){
text( mean(dimred[clustering==i,1],)-0.15,
mean(dimred[clustering==i,2],)+.55, labels = i, cex = 1.75, col = "black", offset = 1, pos = 3) }

#plot(dimred, col = 'white', cex = 0.3, pch = 16)
plot(dimred, col = pal1[clustering], cex = 0, pch = 16)
lines(curves, lwd = 3, col = 'black')
for(i in levels(clustering)){
text( mean(dimred[clustering==i,1],)-0.25,
mean(dimred[clustering==i,2],)-0.9, labels = i, cex = 1.75, col = "black", offset = 1, pos = 3) }

# Split d8 to d12 into two distinct lineages
#

DimPlot(this_subset, group.by = "timepoint", reduction = "umap",pt.size = 0.1, label = F) |
DimPlot(this_subset, group.by = "integrated_snn_res.0.15", reduction = "umap",pt.size = 0.1, label = TRUE, repel = TRUE)
Idents(this_subset) <- this_subset$integrated_snn_res.0.15

97

#subset1 <- subset(this_subset, subset = integrated_snn_res.0.15 %in% c("10", "7", "2", "6", "9", "13", "16",
# "18", "17", "19", "3", "15", "12"))

subset1 <- subset(this_subset, subset = integrated_snn_res.0.15 %in% c("18", "13", "2", "5", "4", "9", "11",
"22", "16", "17", "25", "12", "8", "15"))

table(subset1$integrated_snn_res.0.15)

#subset2 <- subset(this_subset, subset = integrated_snn_res.0.15 %in% c("11", "8", "21", "24", "1", "20", "4",
# "14", "22", "5", "0"))

subset2 <- subset(this_subset, subset = integrated_snn_res.0.15 %in% c("10", "6", "20", "7", "1", "19",
"24", "0", "3", "14", "21"))

table(subset2$integrated_snn_res.0.15)

# Use these for d8 to d12
DimPlot(subset1, group.by = "integrated_snn_res.0.15", reduction = "umap",pt.size = 0.5, label = TRUE, repel = TRUE) |
DimPlot(subset2, group.by = "integrated_snn_res.0.15", reduction = "umap",pt.size = 0.5, label = TRUE, repel = TRUE)

dimred1 <- subset1@reductions$umap@cell.embeddings
dimred2 <- subset2@reductions$umap@cell.embeddings

# Use 0.15 for d8 to d12
clustering1 <- subset1@meta.data$integrated_snn_res.0.15
clustering2 <- subset2@meta.data$integrated_snn_res.0.15

# Use these for d8 to d12
lineages1 <- as.SlingshotDataSet(getLineages(data = dimred1, clusterLabels = clustering1, start.clus = '18')) # 2 is starting cluster for D2_t
lineages2 <- as.SlingshotDataSet(getLineages(data = dimred2, clusterLabels = clustering2, start.clus = '10')) # 2 is starting cluster for D2_t

curves1 <- as.SlingshotDataSet(getCurves(lineages1, approx_points = 500, thresh = 0.01, stretch = 0.8, allow.breaks = T, shrink = 0.99))
curves2 <- as.SlingshotDataSet(getCurves(lineages2, approx_points = 500, thresh = 0.01, stretch = 0.8, allow.breaks = T, shrink = 0.99))

#use these for subset1 and subset2 of d8 to d12

par(mfrow=c(1,2))
#plot(dimred, col = 'white', cex = 0.3, pch = 16)
plot(dimred1, col = pal1[clustering1], cex = 0.3, pch = 16)
lines(curves1, lwd = 3, col = 'black')
for(i in levels(clustering1)){
text( mean(dimred1[clustering1==i,1],)-0.25,
mean(dimred1[clustering1==i,2],)-0.9, labels = i, cex = 1.75, col = "black", offset = 1, pos = 3) }
#plot(dimred, col = 'white', cex = 0.3, pch = 16)

plot(dimred2, col = pal1[clustering2], cex = 0.3, pch = 16)
lines(curves2, lwd = 3, col = 'black')
for(i in levels(clustering2)){
text( mean(dimred2[clustering2==i,1],)-0.25,
mean(dimred2[clustering2==i,2],)-0.9, labels = i, cex = 1.75, col = "black", offset = 1, pos = 3) }

counts1 <- as.matrix(subset1@assays$SCT@counts[VariableFeatures(this_subset), ])
counts2 <- as.matrix(subset2@assays$SCT@counts[VariableFeatures(this_subset), ])
filt1 <- rowSums(counts1 > 5) > ncol(counts1)/100
filt2 <- rowSums(counts2 > 5) > ncol(counts2)/100
counts1 <- counts1[filt1, ]
counts2 <- counts2[filt2, ]

saveRDS.compress.in.BG(subset1, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_subset1.RDS.gz")
saveRDS.compress.in.BG(subset2, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_subset2.RDS.gz")
saveRDS.compress.in.BG(dimred1, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_dimred1.RDS.gz")
98

saveRDS.compress.in.BG(dimred2, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_dimred2.RDS.gz")
saveRDS.compress.in.BG(clustering1, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_clustering1.RDS.gz")
saveRDS.compress.in.BG(clustering2, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_clustering2.RDS.gz")
saveRDS.compress.in.BG(lineages1, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_lineages1.RDS.gz")
saveRDS.compress.in.BG(lineages2, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_lineages2.RDS.gz")

saveRDS.compress.in.BG(counts1, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_counts1.RDS.gz")
saveRDS.compress.in.BG(counts2, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_counts2.RDS.gz")
saveRDS.compress.in.BG(curves1, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_curves1.RDS.gz")
saveRDS.compress.in.BG(curves2, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_curves2.RDS.gz")

subset1 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12_subset1.RDS.gz")
subset2 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12_subset2.RDS.gz")

######## Tempora
########

{

#counts <- readRDS("/data/scRNA/IMR90_timecourse/counts_09-30-22.RDS")
#curves <- readRDS("/data/scRNA/IMR90_timecourse/curves_09-30-22.RDS")
seu_Pro_to_D2 <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2.RDS.gz")

TPs_tempora <- ImportSeuratObject(seu_Pro_to_D2, clusters = "SCT_snn_res.0.175", timepoints = "timepoint", assayType = "SCT",
timepoint_order = c("Pro","H12","D1","D2"))

saveRDS.compress.in.BG(TPs_tempora, fname = "/data/scRNA/IMR90_timecourse/TPs_tempora_09-22-22.RDS")
#use this as the original object to run all Calculate PW Profiles further below
TPs_tempora <- readRDS("/data/scRNA/IMR90_timecourse/TPs_tempora_09-22-22.RDS")

TPs_tempora_PWs <- CalculatePWProfiles(TPs_tempora, gmt_path =
"/data/scRNA/IMR90_timecourse/Human_GO_AllPathways_no_GO_iea_August_03_2022_symbol.gmt",
method="gsva", min.sz = 25, max.sz = 30, parallel.sz = 0)
TPs_tempora_DrugBank <- CalculatePWProfiles(TPs_tempora, gmt_path =
"/data/scRNA/IMR90_timecourse/Human_DrugBank_all_symbol.gmt",
method = "gsva", min.sz = 10, max.sz = 1000, parallel.sz = 0)
TPs_tempora_miRs <- CalculatePWProfiles(TPs_tempora, gmt_path =
"/data/scRNA/IMR90_timecourse/Human_miRs_MSigdb_August_03_2022_symbol.gmt",
method = "gsva", min.sz = 10, max.sz = 1000, parallel.sz = 0)
TPs_tempora_TFs <- CalculatePWProfiles(TPs_tempora, gmt_path =
"/data/scRNA/IMR90_timecourse/Human_TranscriptionFactors_MSigdb_August_03_2022_symbol.gmt",
method = "gsva", min.sz = 10, max.sz = 1000, parallel.sz = 0)

#TPs_tempora_list <- list(TPs_tempora_PWs, TPs_tempora_DrugBank, TPs_tempora_miRs, TPs_tempora_TFs)
TPs_tempora_list <- list(TPs_tempora_PWs)
saveRDS.compress.in.BG(TPs_tempora_list, fname = "/data/scRNA/IMR90_timecourse/initial_TPs_tempora_list_09-22-22.RDS")
for (n in seq_along(along.with = TPs_tempora_list)){

print(n)
TPs_tempora_list[[n]] <- BuildTrajectory(TPs_tempora_list[[n]], n_pcs = 11)
TPs_tempora_list[[n]]@cluster.metadata$label <- TPs_tempora_list[[n]]@cluster.metadata$Id
TPs_tempora_list[[n]] <- invisible(IdentifyVaryingPWs(TPs_tempora_list[[n]], pval_threshold = 0.05))
#pdf("/data/scRNA/IMR90_timecourse/test.PDF", width = 10)
TPs_tempora_list[[n]] <- invisible(PlotTrajectory(TPs_tempora_list[[n]]))
#dev.off()
}

pdf("/data/scRNA/IMR90_timecourse/TPs_tempora_plotVaryingPWs_AllPathways_09-30-22.pdf")
PlotVaryingPWs(TPs_tempora_list[[1]])
dev.off()

}
99

APPENDIX C
Slingshot_Tradeseq_workflow.R – R code used to plot and visualize lineages and trajectories across
stages of senescence.

.libPaths("/opt/R/4.0.4/lib/R/library")
setwd("/data/scRNA/IMR90_timecourse/")
#Sys.unsetenv("GITHUB_PAT")

library(SingleCellExperiment)
library(Tempora)
library(Seurat)
library(Seurat.utils)
library(SeuratWrappers)
library(Matrix)
library(dyndimred)
library(Seurat.utils)
library(diffusionMap)
library(slingshot)
library(tradeSeq)
library(paletteer)
library(hsslda)
library(khroma)
library(palettesForR)
library(pals)
library(randomcoloR)
library(RColorBrewer)
library(igraph)
library(psupertime)
library(scanalysis)
library(dplyr)
library(tidyr)
library(tibble)
library(patchwork)
library(grDevices)
library(hues)

#BiocParallel::register(BiocParallel::SerialParam())

seu_Pro_to_D2 <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_singletsOnly.RDS.gz")
seu_D2_to_D8 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D2_to_D8_singletsOnly.RDS.gz")
seu_D8_to_D12_s1 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s1_singletsOnly.RDS.gz")
seu_D8_to_D12_s2 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s2_singletsOnly.RDS.gz")

seu_D8_to_D12_s1 <- RunUMAP(seu_D8_to_D12_s1, repulsion.strength = 10, min.dist = 0.2, reduction = "pca", dims = 1:50, verbose = T)
seu_D8_to_D12_s1 <- FindNeighbors(seu_D8_to_D12_s1, k.param = 20, reduction = "umap", dims = 1:2, verbose = T)
seu_D8_to_D12_s1 <- FindClusters(seu_D8_to_D12_s1, resolution = 0.15, verbose = T)

#DimPlot(seu_Pro_to_D2, group.by = "timepoint", shuffle = T) | DimPlot(seu_Pro_to_D2, group.by = "batch", split.by = "timepoint", shuffle = T)
#DimPlot(seu_D2_to_D8, group.by = "timepoint", shuffle = T) | DimPlot(seu_D2_to_D8, group.by = "batch", split.by = "timepoint", shuffle = T)
DimPlot(seu_D8_to_D12_s1, group.by = "timepoint", shuffle = T)# | DimPlot(seu_D8_to_D12_s1, group.by = "batch", split.by = "timepoint",
shuffle = T)
#DimPlot(seu_D8_to_D12_s2, group.by = "timepoint", shuffle = T) | DimPlot(seu_D8_to_D12_s2, group.by = "batch", split.by = "timepoint",
shuffle = T)

#Pro_to_D2_clustering <- seu_Pro_to_D2@meta.data$integrated_snn_res.0.2
#D2_to_D8_clustering <- seu_D2_to_D8@meta.data$integrated_snn_res.0.25
#D8_to_D12_s1_clustering <- seu_D8_to_D12_s1@meta.data$integrated_snn_res.0.15
100

#D8_to_D12_s2_clustering <- seu_D8_to_D12_s2@meta.data$integrated_snn_res.0.15
#Pro_to_D2_dimred <- seu_Pro_to_D2@reductions$umap@cell.embeddings
#D2_to_D8_dimred <- seu_D2_to_D8@reductions$umap@cell.embeddings
#D8_to_D12_s1_dimred <- seu_D8_to_D12_s1@reductions$umap@cell.embeddings
#D8_to_D12_s2_dimred <- seu_D8_to_D12_s2@reductions$umap@cell.embeddings

Pro_to_D2_curves <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_curves_singlets.RDS.gz")
D2_to_D8_curves <- readRDS("/data/scRNA/IMR90_timecourse/seu_D2_to_D8_curves_singlets.RDS.gz")
D8_to_D12_s1_curves <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s1_curves_singlets.RDS.gz")
D8_to_D12_s2_curves <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s2_curves_singlets.RDS.gz")

set.seed(1)

##### use for pro to d2
DimPlot(seu_Pro_to_D2, group.by = "timepoint", label = T)

seu_Pro_to_D2_clustering <- seu_Pro_to_D2@meta.data$integrated_snn_res.0.2
seu_Pro_to_D2_dimred <- seu_Pro_to_D2@reductions$umap@cell.embeddings
#saveRDS.compress.in.BG(seu_Pro_to_D2_clustering, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_clustering_singlets.RDS")
#saveRDS.compress.in.BG(seu_Pro_to_D2_dimred, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_dimred_singlets.RDS")

seu_Pro_to_D2_lineages <- as.SlingshotDataSet(getLineages(data = seu_Pro_to_D2_dimred, clusterLabels = seu_Pro_to_D2_clustering,
start.clus = '18'))
seu_Pro_to_D2_lineages
seu_Pro_to_D2_curves <- as.SlingshotDataSet(getCurves(seu_Pro_to_D2_lineages, approx_points = 150, stretch = 0, allow.breaks = T, shrink =
0.99))
seu_Pro_to_D2_curves
saveRDS.compress.in.BG(seu_Pro_to_D2_curves, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_curves_singlets.RDS")

seu_Pro_to_D2_counts <- as.matrix(seu_Pro_to_D2@assays$SCT@counts[VariableFeatures(seu_Pro_to_D2), ])
filt <- rowSums(seu_Pro_to_D2_counts > 5) > ncol(seu_Pro_to_D2_counts)/1000
sum(filt)
seu_Pro_to_D2_counts <- seu_Pro_to_D2_counts[filt, ]
saveRDS.compress.in.BG(seu_Pro_to_D2_counts, fname = "/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_counts_singlets.RDS")

#
# use for d2 to d8

DimPlot(seu_D2_to_D8, group.by = "timepoint", label = T)
DimPlot(seu_D2_to_D8, group.by = "integrated_snn_res.0.25", label = T)

seu_D2_to_D8_clustering <- seu_D2_to_D8@meta.data$integrated_snn_res.0.25
seu_D2_to_D8_dimred <- seu_D2_to_D8@reductions$umap@cell.embeddings
saveRDS.compress.in.BG(seu_D2_to_D8_clustering, fname = "/data/scRNA/IMR90_timecourse/seu_D2_to_D8_clustering_singlets.RDS")
saveRDS.compress.in.BG(seu_D2_to_D8_dimred, fname = "/data/scRNA/IMR90_timecourse/seu_D2_to_D8_dimred_singlets.RDS")

seu_D2_to_D8_lineages <- as.SlingshotDataSet(getLineages(data = seu_D2_to_D8_dimred, clusterLabels = seu_D2_to_D8_clustering, start.clus
= '17'))
seu_D2_to_D8_lineages
seu_D2_to_D8_curves <- as.SlingshotDataSet(getCurves(seu_D2_to_D8_lineages, approx_points = 150, stretch = 0, allow.breaks = T, shrink =
0.99))
seu_D2_to_D8_curves
saveRDS.compress.in.BG(seu_D2_to_D8_curves, fname = "/data/scRNA/IMR90_timecourse/seu_D2_to_D8_curves_singlets.RDS")

seu_D2_to_D8_counts <- as.matrix(seu_D2_to_D8@assays$SCT@counts[VariableFeatures(seu_D2_to_D8), ])
filt <- rowSums(seu_D2_to_D8_counts > 5) > ncol(seu_D2_to_D8_counts)/5000
sum(filt)
seu_D2_to_D8_counts <- seu_D2_to_D8_counts[filt, ]
saveRDS.compress.in.BG(seu_D2_to_D8_counts, fname = "/data/scRNA/IMR90_timecourse/seu_D2_to_D8_counts_singlets.RDS")

101

##################################################################################################################
# use for d8 to d12 s1

DimPlot(seu_D8_to_D12_s1, group.by = "integrated_snn_res.0.15", label = T, label.size = 10)
seu_D8_to_D12_s1_clustering <- seu_D8_to_D12_s1@meta.data$integrated_snn_res.0.15
seu_D8_to_D12_s1_dimred <- seu_D8_to_D12_s1@reductions$umap@cell.embeddings
#saveRDS.compress.in.BG(seu_D8_to_D12_s1_clustering, fname =
"/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s1_clustering_singlets.RDS")
#saveRDS.compress.in.BG(seu_D8_to_D12_s1_dimred, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s1_dimred_singlets.RDS")

seu_D8_to_D12_s1_lineages <- as.SlingshotDataSet(getLineages(data = seu_D8_to_D12_s1_dimred, clusterLabels =
seu_D8_to_D12_s1_clustering, start.clus = '18'))
seu_D8_to_D12_s1_lineages
seu_D8_to_D12_s1_curves <- as.SlingshotDataSet(getCurves(seu_D8_to_D12_s1_lineages, approx_points = 150, stretch = 0, allow.breaks = T,
shrink = 0.99))
seu_D8_to_D12_s1_curves
saveRDS.compress.in.BG(seu_D8_to_D12_s1_curves, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s1_curves_singlets.RDS")

seu_D8_to_D12_s1_counts <- as.matrix(seu_D8_to_D12_s1@assays$SCT@counts[VariableFeatures(seu_D8_to_D12_s1), ])
filt <- rowSums(seu_D8_to_D12_s1_counts > 5) > ncol(seu_D8_to_D12_s1_counts)/50000
sum(filt)
seu_D8_to_D12_s1_counts <- seu_D8_to_D12_s1_counts[filt, ]
saveRDS.compress.in.BG(seu_D8_to_D12_s1_counts, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s1_counts_singlets.RDS")

#
# use for d8 to d12 s2

DimPlot(seu_D8_to_D12_s2, group.by = "integrated_snn_res.0.15", label = T)
seu_D8_to_D12_s2_clustering <- seu_D8_to_D12_s2@meta.data$integrated_snn_res.0.15
seu_D8_to_D12_s2_dimred <- seu_D8_to_D12_s2@reductions$umap@cell.embeddings
#saveRDS.compress.in.BG(seu_D8_to_D12_s2_clustering, fname =
"/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s2_clustering_singlets.RDS")
#saveRDS.compress.in.BG(seu_D8_to_D12_s2_dimred, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s2_dimred_singlets.RDS")

seu_D8_to_D12_s2_lineages <- as.SlingshotDataSet(getLineages(data = seu_D8_to_D12_s2_dimred, clusterLabels =
seu_D8_to_D12_s2_clustering, start.clus = '10'))
seu_D8_to_D12_s2_lineages
seu_D8_to_D12_s2_curves <- as.SlingshotDataSet(getCurves(seu_D8_to_D12_s2_lineages, approx_points = 150, stretch = 0, allow.breaks = T,
shrink = 0.99))
seu_D8_to_D12_s2_curves
saveRDS.compress.in.BG(seu_D8_to_D12_s1_curves, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s1_curves_singlets.RDS")

seu_D8_to_D12_s2_counts <- as.matrix(seu_D8_to_D12_s2@assays$SCT@counts[VariableFeatures(seu_D8_to_D12_s2), ])
filt <- rowSums(seu_D8_to_D12_s2_counts > 5) > ncol(seu_D8_to_D12_s2_counts)/20000
sum(filt)
seu_D8_to_D12_s2_counts <- seu_D8_to_D12_s2_counts[filt, ]
saveRDS.compress.in.BG(seu_D8_to_D12_s2_counts, fname = "/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s2_counts_singlets.RDS")

#
# Plot slingshot/tradeseq workflow lineages

axis_label_size = 10
axis_text_size = 8
panel_tag_size = 16

hex_codes <- iwanthue(n = 22)
hex_codes[21] <- "#989C4C"

#
# Pro to D2

Pro_to_D2_curves_DF <- slingCurves(seu_Pro_to_D2_curves, as.df = TRUE) %>% arrange(Order)

102

pal1 <- paletteer_d(`"palettesForR::Bears"`, n = 45, direction = -1)
plot(seu_Pro_to_D2_dimred, col = pal1[seu_Pro_to_D2_clustering], cex = 1, pch = 16)
lines(seu_Pro_to_D2_lineages, lwd = 3, col = 'black')
#for(i in levels(seu_Pro_to_D2_clustering)){
# text( mean(seu_Pro_to_D2_dimred[seu_Pro_to_D2_clustering==i,1],)-0.25,
# mean(seu_Pro_to_D2_dimred[seu_Pro_to_D2_clustering==i,2],)-0.9, labels = i, cex = 1.75, col = "black", offset = 1, pos = 3) }

seu_Pro_to_D2_p1 <- DimPlot(seu_Pro_to_D2, group.by = "timepoint", label = F, pt.size = 0.01) +
ggtitle('') +
NoLegend() +
theme(legend.position = "bottom") +
theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold")) +
labs(tag = "A")

#seu_Pro_to_D2_dimred
#seu_Pro_to_D2_clustering

seu_Pro_to_D2_p2 <- ggplot(as.data.frame(seu_Pro_to_D2_dimred), aes(x = UMAP_1, y = UMAP_2)) +
geom_point(aes(color = seu_Pro_to_D2_clustering), shape = 20, size = 0.1) +
scale_color_manual(values = hex_codes) +
#lines(seu_Pro_to_D2_lineages, lwd = 3, col = 'black') +
geom_path(data = Pro_to_D2_curves_DF, aes(group = Lineage), size = 1.5, arrow = arrow(length = unit(0.25,"cm"), ends = "last", type =
"closed")) +
theme_classic() +
theme(legend.position="none") +
theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold")) +
labs(tag = "B")

(seu_Pro_to_D2_p1 / seu_Pro_to_D2_p2)

#
# D2 to D8

D2_to_D8_curves_DF <- slingCurves(seu_D2_to_D8_curves, as.df = TRUE) %>% arrange(Lineage)

# Subset to deal with weird glitch - cells from wrong cluster added at the end; these cells actually belong to the origin cluster.
# Maybe during lineage/curve building, it's trying to even out cell numbers and grabs the closest/lowest numbered cell.
D2_to_D8_curves_DF_ss1 <- subset(D2_to_D8_curves_DF, D2_to_D8_curves_DF$Lineage %in% c(2, 4, 5, 6) |
(D2_to_D8_curves_DF$Lineage == 7 & D2_to_D8_curves_DF$Order <= 105) |
(D2_to_D8_curves_DF$Lineage == 1 & D2_to_D8_curves_DF$Order <= 136) |
(D2_to_D8_curves_DF$Lineage == 3 & D2_to_D8_curves_DF$Order <= 137))

seu_D2_to_D8_p1 <- DimPlot(seu_D2_to_D8, group.by = "timepoint", label = F, pt.size = 0.01) +
ggtitle('') +
NoLegend() +
theme(legend.position = "bottom") +
theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold")) +
labs(tag = "A")

#seu_D2_to_D8_clustering
#seu_D2_to_D8_dimred

seu_D2_to_D8_p2 <- ggplot(as.data.frame(seu_D2_to_D8_dimred), aes(x = UMAP_1, y = UMAP_2)) +
geom_point(aes(color = seu_D2_to_D8_clustering), shape = 20, size = 0.1) +
scale_color_manual(values = hex_codes) +
geom_path(data = D2_to_D8_curves_DF_ss1, aes(group = Lineage), size = 1.5, arrow = arrow(length = unit(0.25,"cm"), ends = "last", type
= "closed")) +
theme_classic() +
theme(legend.position="none") +
103

theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold")) +
labs(tag = "B")

(seu_D2_to_D8_p1 / seu_D2_to_D8_p2)

#
# D8 to D12, s1

D8_to_D12_s1_curves_DF <- slingCurves(seu_D8_to_D12_s1_curves, as.df = TRUE) %>% arrange(Lineage)

seu_D8_to_D12_s1_p1 <- DimPlot(seu_D8_to_D12_s1, group.by = "timepoint", label = F, pt.size = 0.01) +
ggtitle('') +
NoLegend() +
theme(legend.position = "bottom") +
theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold")) +
labs(tag = "A")

seu_D8_to_D12_s1_dimred
seu_D8_to_D12_s1_clustering

seu_D8_to_D12_s1_p2 <- ggplot(as.data.frame(seu_D8_to_D12_s1_dimred), aes(x = UMAP_1, y = UMAP_2)) +
geom_point(aes(color = seu_D8_to_D12_s1_clustering), shape = 20, size = 0.1) +
scale_color_manual(values = hex_codes) +
geom_path(data = D8_to_D12_s1_curves_DF, aes(group = Lineage), size = 1.5, arrow = arrow(length = unit(0.25,"cm"), ends = "last",
type = "closed")) +
theme_classic() +
theme(legend.position="none") +
theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold")) +
labs(tag = "B")

(seu_D8_to_D12_s1_p1 | seu_D8_to_D12_s1_p2)

#
# D8 to D12, s2

D8_to_D12_s2_curves_DF <- slingCurves(seu_D8_to_D12_s2_curves, as.df = TRUE) %>% arrange(Lineage)

seu_D8_to_D12_s2_p1 <- DimPlot(seu_D8_to_D12_s2, group.by = "timepoint", label = F, pt.size = 0.01) +
ggtitle('') +
NoLegend() +
theme(legend.position = "bottom") +
theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold")) +
labs(tag = "A")

seu_D8_to_D12_s2_dimred
seu_D8_to_D12_s2_clustering

seu_D8_to_D12_s2_p2 <- ggplot(as.data.frame(seu_D8_to_D12_s2_dimred), aes(x = UMAP_1, y = UMAP_2)) +
geom_point(aes(color = seu_D8_to_D12_s2_clustering), shape = 20, size = 0.1) +
scale_color_manual(values = hex_codes) +
geom_path(data = D8_to_D12_s2_curves_DF, aes(group = Lineage), size = 1.5, arrow = arrow(length = unit(0.25,"cm"), ends = "last",
type = "closed")) +
theme_classic() +
theme(legend.position="none") +
104

theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold")) +
labs(tag = "B")

(seu_D8_to_D12_s2_p1 | seu_D8_to_D12_s2_p2)
(seu_Pro_to_D2_p1 / seu_Pro_to_D2_p2)
(seu_D2_to_D8_p1 / seu_D2_to_D8_p2)
(seu_D8_to_D12_s1_p1 / seu_D8_to_D12_s1_p2)
(seu_D8_to_D12_s2_p1 / seu_D8_to_D12_s2_p2)

seu_D8_to_D12 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12_forLateTrajectoryPlotting.RDS.gz")
seu_D8_to_D12_Full <- DimPlot(seu_D8_to_D12, group.by = "timepoint", reduction = "umap",pt.size = 0.01, label = F) +
#ggtitle('Late Senescence Subpopulations') +
ggtitle('') +
NoLegend() +
theme(text = element_text(size = axis_label_size), axis.text = element_text(size = axis_text_size), plot.tag = element_text(size =
panel_tag_size, face = "bold")) +
labs(tag = "A")

seu_D8_to_D12_Full / (seu_D8_to_D12_s1_p2 | seu_D8_to_D12_s2_p2) / (D8_to_D12_s1_p3 | D8_to_D12_s2_p3)
seu_D8_to_D12_Full / (seu_D8_to_D12_s1_p2 | seu_D8_to_D12_s2_p2)

APPENDIX D
Analyze_SeuratDEGs_byLineage.R – R code for analyzing gene expression changes and cluster
composition after generating trajectories with Slingshot and Tradeseq, based on initial Seurat
clustering.
.libPaths("/opt/R/4.0.4/lib/R/library")
setwd("/data/scRNA/IMR90_timecourse/")
#Sys.unsetenv("GITHUB_PAT")

library(SingleCellExperiment)
library(Tempora)
library(Seurat)
library(Seurat.utils)
library(SeuratWrappers)
library(Matrix)
library(dyndimred)
library(Seurat.utils)
library(diffusionMap)
library(slingshot)
BiocParallel::register(BiocParallel::SerialParam())
library(tradeSeq)
library(paletteer)
library(hsslda)
library(khroma)
library(palettesForR)
library(pals)
library(randomcoloR)
library(RColorBrewer)
library(igraph)
library(psupertime)
library(scanalysis)
library(dplyr)
library(tidyr)
library(tibble)
library(patchwork)

#####################################################################################################################
############################################################
#

105

#get top/bottom row
head_and_tail <- function(x, n = 1) {
rbind(
head(x, n),
tail(x, n)
)
}

#ProD2 <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_this_subset.RDS.gz")
#DefaultAssay(ProD2) <- "SCT"
#ProD2 <- PrepSCTFindMarkers(ProD2)
ProD2 <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_singletsOnly_prepSCTfindmarkers.RDS.gz")
ProD2curves <- readRDS("/data/scRNA/IMR90_timecourse/seu_Pro_to_D2_curves_singlets.RDS.gz")

DimPlot(ProD2, group.by = "seurat_clusters", label = T, label.box = T)
ProD2curves #lineage end clusters: 9, 13, 15, 16, 19, 21

ProD2_clusterlist <- list(9, 13, 15, 16, 19, 21)
Idents(ProD2) <- ProD2$seurat_clusters
DefaultAssay(ProD2) <- "SCT"

ProD2_c13_markers <- FindMarkers(ProD2, ident.1 = "13", ident.2 = c("15","21","16","19", "9"), min.pct = 0.01, only.pos = F, verbose = T)
#exclude cluster 9 (lineage 6) bc shortest pseudotime lineage
ProD2_c13_markers$cluster <- "L1c13"
ProD2_c13_markers$gene <- rownames(ProD2_c13_markers)

ProD2_c15_markers <- FindMarkers(ProD2, ident.1 = "15", ident.2 = c("13","21","16","19", "9"), min.pct = 0.01, only.pos = F, verbose = T)
#exclude cluster 9 (lineage 6) bc shortest pseudotime lineage
ProD2_c15_markers$cluster <- "L2c15"
ProD2_c15_markers$gene <- rownames(ProD2_c15_markers)

ProD2_c21_markers <- FindMarkers(ProD2, ident.1 = "21", ident.2 = c("13","15","16","19", "9"), min.pct = 0.01, only.pos = F, verbose = T)
#exclude cluster 9 (lineage 6) bc shortest pseudotime lineage
ProD2_c21_markers$cluster <- "L3c21"
ProD2_c21_markers$gene <- rownames(ProD2_c21_markers)

ProD2_c16_markers <- FindMarkers(ProD2, ident.1 = "16", ident.2 = c("13","15","21","19", "9"), min.pct = 0.01, only.pos = F, verbose = T)
#exclude cluster 9 (lineage 6) bc shortest pseudotime lineage
ProD2_c16_markers$cluster <- "L4c16"
ProD2_c16_markers$gene <- rownames(ProD2_c16_markers)

ProD2_c19_markers <- FindMarkers(ProD2, ident.1 = "19", ident.2 = c("13","15","21","16", "9"), min.pct = 0.01, only.pos = F, verbose = T)
#exclude cluster 9 (lineage 6) bc shortest pseudotime lineage
ProD2_c19_markers$cluster <- "L5c19"
ProD2_c19_markers$gene <- rownames(ProD2_c19_markers)

ProD2_c9_markers <- FindMarkers(ProD2, ident.1 = "9", ident.2 = c("13","15","21","16","19"), min.pct = 0.01, only.pos = F, verbose = T)
ProD2_c9_markers$cluster <- "L6c9"
ProD2_c9_markers$gene <- rownames(ProD2_c9_markers)

#rbind(head_and_tail(ProD2_c9_markers), head_and_tail(ProD2_c13_markers), head_and_tail(ProD2_c15_markers),
head_and_tail(ProD2_c21_markers), head_and_tail(ProD2_c16_markers), head_and_tail(ProD2_c19_markers))

ProD2_lineageEnd_DEGs <- rbind(ProD2_c9_markers, ProD2_c13_markers, ProD2_c15_markers, ProD2_c21_markers, ProD2_c16_markers,
ProD2_c19_markers)
ProD2_lineageEnd_DEGs$p_val_adj = p.adjust(ProD2_lineageEnd_DEGs$p_val, method='BH')
ProD2_lineageEnd_DEGs$delta <- abs(as.double(ProD2_lineageEnd_DEGs$pct.1)-as.double(ProD2_lineageEnd_DEGs$pct.2))
ProD2_lineageEnd_DEGs$minor <- ifelse(ProD2_lineageEnd_DEGs$pct.1 < 0.10 | ProD2_lineageEnd_DEGs$pct.2 < 0.10, 1, 0)
ProD2_lineageEnd_DEGs <- ProD2_lineageEnd_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
ProD2_lineageEnd_DEGs <- ProD2_lineageEnd_DEGs[ProD2_lineageEnd_DEGs$p_val_adj<0.05,]
106

#sorted by log2FC
#ProD2_lineageEnd_DEGs <- ProD2_lineageEnd_DEGs[order(ProD2_lineageEnd_DEGs[,1], -ProD2_lineageEnd_DEGs[,3],
ProD2_lineageEnd_DEGs[,4] ),]
#write.table(ProD2_lineageEnd_DEGs, "/data/scRNA/IMR90_timecourse/ProD2_lineageEnd_DEGs.tsv", quote = F, sep = "\t", row.names = F)

#Lineage1: 18 20 14 11 2 5 8 10 7 0 19
#Lineage2: 18 20 14 11 2 5 8 10 3 6 16
#Lineage3: 18 20 14 11 2 5 8 10 7 12 21
#Lineage4: 18 20 14 11 2 5 8 1 4 15
#Lineage5: 18 20 14 11 2 13
#Lineage6: 18 20 14 11 17 9

ProD2_L1split_0vs12_markers <- FindMarkers(ProD2, ident.1 = "0", ident.2 = c("12"), min.pct = 0.01, only.pos = F, verbose = T)
ProD2_L1split_0vs12_markers$cluster <- "L1c0"
ProD2_L1split_0vs12_markers$gene <- rownames(ProD2_L1split_0vs12_markers)

ProD2_L2split_3vs7_markers <- FindMarkers(ProD2, ident.1 = "3", ident.2 = c("7"), min.pct = 0.01, only.pos = F, verbose = T)
ProD2_L2split_3vs7_markers$cluster <- "L2c3"
ProD2_L2split_3vs7_markers$gene <- rownames(ProD2_L2split_3vs7_markers)

ProD2_L3split_12vs0_markers <- FindMarkers(ProD2, ident.1 = "12", ident.2 = c("0"), min.pct = 0.01, only.pos = F, verbose = T)
ProD2_L3split_12vs0_markers$cluster <- "L3c12"
ProD2_L3split_12vs0_markers$gene <- rownames(ProD2_L3split_12vs0_markers)

ProD2_L4split_1vs10_markers <- FindMarkers(ProD2, ident.1 = "1", ident.2 = c("10"), min.pct = 0.01, only.pos = F, verbose = T)
ProD2_L4split_1vs10_markers$cluster <- "L4c1"
ProD2_L4split_1vs10_markers$gene <- rownames(ProD2_L4split_1vs10_markers)

ProD2_L5split_13vs5_markers <- FindMarkers(ProD2, ident.1 = "13", ident.2 = c("5"), min.pct = 0.01, only.pos = F, verbose = T)
ProD2_L5split_13vs5_markers$cluster <- "L5c13"
ProD2_L5split_13vs5_markers$gene <- rownames(ProD2_L5split_13vs5_markers)

ProD2_L6split_17vs2_markers <- FindMarkers(ProD2, ident.1 = "17", ident.2 = c("2"), min.pct = 0.01, only.pos = F, verbose = T)
ProD2_L6split_17vs2_markers$cluster <- "L6c17"
ProD2_L6split_17vs2_markers$gene <- rownames(ProD2_L6split_17vs2_markers)

ProD2_lineageSwitch_DEGs <- rbind(ProD2_L6split_17vs2_markers, ProD2_L5split_13vs5_markers, ProD2_L4split_1vs10_markers,
ProD2_L2split_3vs7_markers, ProD2_L3split_12vs0_markers, ProD2_L1split_0vs12_markers)
ProD2_lineageSwitch_DEGs$p_val_adj = p.adjust(ProD2_lineageSwitch_DEGs$p_val, method='BH')
ProD2_lineageSwitch_DEGs$delta <- abs(as.double(ProD2_lineageSwitch_DEGs$pct.1)-as.double(ProD2_lineageSwitch_DEGs$pct.2))
ProD2_lineageSwitch_DEGs$minor <- ifelse(ProD2_lineageSwitch_DEGs$pct.1 < 0.10 | ProD2_lineageSwitch_DEGs$pct.2 < 0.10, 1, 0)
ProD2_lineageSwitch_DEGs <- ProD2_lineageSwitch_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
ProD2_lineageSwitch_DEGs <- ProD2_lineageSwitch_DEGs[ProD2_lineageSwitch_DEGs$p_val_adj<0.05,]
write.table(ProD2_lineageSwitch_DEGs, "/data/scRNA/IMR90_timecourse/ProD2_lineageSwitch_DEGs.tsv", quote = F, sep = "\t", row.names =
F)

#upregulated, sorted by log2FC
t1 <- ProD2_lineageEnd_DEGs[order(ProD2_lineageEnd_DEGs[,1], -ProD2_lineageEnd_DEGs[,3]), ]
write.table(t1, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/ProD2_LinEnd_upregulated_FCsorted.tsv", quote = F, sep = "\t", row.names =
F)
#downregulated, sorted by log2FC
t2 <- ProD2_lineageEnd_DEGs[order(ProD2_lineageEnd_DEGs[,1], ProD2_lineageEnd_DEGs[,3]), ]
write.table(t2, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/ProD2_LinEnd_downregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by delta
t3 <- subset(ProD2_lineageEnd_DEGs[order(ProD2_lineageEnd_DEGs[,1], -ProD2_lineageEnd_DEGs[,7], -ProD2_lineageEnd_DEGs[,4]), ],
avg_log2FC > 0)
write.table(t3, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/ProD2_LinEnd_upregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by delta
t4 <- subset(ProD2_lineageEnd_DEGs[order(ProD2_lineageEnd_DEGs[,1], -ProD2_lineageEnd_DEGs[,7], -ProD2_lineageEnd_DEGs[,4]), ],
avg_log2FC < 0)
107

write.table(t4, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/ProD2_LinEnd_downregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by log2FC
t5 <- ProD2_lineageSwitch_DEGs[order(ProD2_lineageSwitch_DEGs[,1], -ProD2_lineageSwitch_DEGs[,3]), ]
write.table(t5, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/ProD2_LinSwitch_upregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by log2FC
t6 <- ProD2_lineageSwitch_DEGs[order(ProD2_lineageSwitch_DEGs[,1], ProD2_lineageSwitch_DEGs[,3]), ]
write.table(t6, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/ProD2_LinSwitch_downregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by delta
t7 <- subset(ProD2_lineageSwitch_DEGs[order(ProD2_lineageSwitch_DEGs[,1], -ProD2_lineageSwitch_DEGs[,7], -
ProD2_lineageSwitch_DEGs[,4]), ], avg_log2FC > 0)
write.table(t7, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/ProD2_LinSwitch_upregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by delta
t8 <- subset(ProD2_lineageSwitch_DEGs[order(ProD2_lineageSwitch_DEGs[,1], -ProD2_lineageSwitch_DEGs[,7], -
ProD2_lineageSwitch_DEGs[,4]), ], avg_log2FC < 0)
write.table(t8, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/ProD2_LinSwitch_downregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)

#

D2D8 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D2_to_D8_singletsOnly_prepSCTfindmarkers.RDS.gz")
D2D8curves <- readRDS("/data/scRNA/IMR90_timecourse/seu_D2_to_D8_curves_singlets.RDS.gz")

#lineages: 7
#Lineage1: 17 8 4 0 2 11 7 18
#Lineage2: 17 8 4 0 2 15 12 6
#Lineage3: 17 8 4 0 2 11 7 9
#Lineage4: 17 8 4 14 13 10
#Lineage5: 17 8 1 5 3
#Lineage6: 17 8 4 14 16
#Lineage7: 17 8 4 19
#
#

DimPlot(D2D8, group.by = "seurat_clusters", label = T, label.box = T)
D2D8curves
D2D8_clusterlist <- list("18", "6", "9", "10", "3", "16", "19")
Idents(D2D8) <- D2D8$seurat_clusters
DefaultAssay(D2D8) <- "SCT"

D2D8_c18_markers <- FindMarkers(D2D8, ident.1 = "18", ident.2 = c("6","9","10", "3", "16", "19"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_c18_markers$cluster <- "L1c18"
D2D8_c18_markers$gene <- rownames(D2D8_c18_markers)

D2D8_c6_markers <- FindMarkers(D2D8, ident.1 = "6", ident.2 = c("18","9","10", "3", "16", "19"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_c6_markers$cluster <- "L2c6"
D2D8_c6_markers$gene <- rownames(D2D8_c6_markers)

D2D8_c9_markers <- FindMarkers(D2D8, ident.1 = "9", ident.2 = c("6","18","10", "3", "16", "19"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_c9_markers$cluster <- "L3c9"
D2D8_c9_markers$gene <- rownames(D2D8_c9_markers)

D2D8_c10_markers <- FindMarkers(D2D8, ident.1 = "10", ident.2 = c("6","9","18", "3", "16", "19"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_c10_markers$cluster <- "L4c10"
D2D8_c10_markers$gene <- rownames(D2D8_c10_markers)

D2D8_c3_markers <- FindMarkers(D2D8, ident.1 = "3", ident.2 = c("18", "6", "9", "10", "16", "19"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_c3_markers$cluster <- "L5c3"
D2D8_c3_markers$gene <- rownames(D2D8_c3_markers)

D2D8_c16_markers <- FindMarkers(D2D8, ident.1 = "16", ident.2 = c("18", "6", "9", "10", "3", "19"), min.pct = 0.01, only.pos = F, verbose = T)
108

D2D8_c16_markers$cluster <- "L6c16"
D2D8_c16_markers$gene <- rownames(D2D8_c16_markers)

D2D8_c19_markers <- FindMarkers(D2D8, ident.1 = "19", ident.2 = c("18", "6", "9", "10", "3", "16"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_c19_markers$cluster <- "L7c19"
D2D8_c19_markers$gene <- rownames(D2D8_c19_markers)

D2D8_lineageEnd_DEGs <-
rbind(D2D8_c18_markers,D2D8_c6_markers,D2D8_c9_markers,D2D8_c10_markers,D2D8_c3_markers,D2D8_c16_markers,D2D8_c19_markers
)

D2D8_lineageEnd_DEGs$p_val_adj = p.adjust(D2D8_lineageEnd_DEGs$p_val, method='BH')
D2D8_lineageEnd_DEGs$delta <- abs(as.double(D2D8_lineageEnd_DEGs$pct.1)-as.double(D2D8_lineageEnd_DEGs$pct.2))
D2D8_lineageEnd_DEGs$minor <- ifelse(D2D8_lineageEnd_DEGs$pct.1 < 0.10 | D2D8_lineageEnd_DEGs$pct.2 < 0.10, 1, 0)
D2D8_lineageEnd_DEGs <- D2D8_lineageEnd_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
D2D8_lineageEnd_DEGs <- D2D8_lineageEnd_DEGs[D2D8_lineageEnd_DEGs$p_val_adj<0.05,]
#write.table(D2D8_lineageEnd_DEGs, "/data/scRNA/IMR90_timecourse/D2D8_lineageEnd_DEGs.tsv", quote = F, sep = "\t", row.names = F)

#7 lineages
#Lineage1: 17 8 4 0 2 11 7 18
#Lineage2: 17 8 4 0 2 15 12 6
#Lineage3: 17 8 4 0 2 11 7 9
#Lineage4: 17 8 4 14 13 10
#Lineage5: 17 8 1 5 3
#Lineage6: 17 8 4 14 16
#Lineage7: 17 8 4 19

D2D8_L1split_18vs9_markers <- FindMarkers(D2D8, ident.1 = "18", ident.2 = c("9"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_L1split_18vs9_markers$cluster <- "L1c18"
D2D8_L1split_18vs9_markers$gene <- rownames(D2D8_L1split_18vs9_markers)

D2D8_L2split_15vs11_markers <- FindMarkers(D2D8, ident.1 = "15", ident.2 = c("11"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_L2split_15vs11_markers$cluster <- "L2c15"
D2D8_L2split_15vs11_markers$gene <- rownames(D2D8_L2split_15vs11_markers)

D2D8_L3split_9vs18_markers <- FindMarkers(D2D8, ident.1 = "9", ident.2 = c("18"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_L3split_9vs18_markers$cluster <- "L3c9"
D2D8_L3split_9vs18_markers$gene <- rownames(D2D8_L3split_9vs18_markers)

D2D8_L4split_13vs16_markers <- FindMarkers(D2D8, ident.1 = "13", ident.2 = c("16"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_L4split_13vs16_markers$cluster <- "L4c13"
D2D8_L4split_13vs16_markers$gene <- rownames(D2D8_L4split_13vs16_markers)

D2D8_L5split_1vs4_markers <- FindMarkers(D2D8, ident.1 = "1", ident.2 = c("4"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_L5split_1vs4_markers$cluster <- "L5c1"
D2D8_L5split_1vs4_markers$gene <- rownames(D2D8_L5split_1vs4_markers)

D2D8_L6split_16vs13_markers <- FindMarkers(D2D8, ident.1 = "16", ident.2 = c("13"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_L6split_16vs13_markers$cluster <- "L6c16"
D2D8_L6split_16vs13_markers$gene <- rownames(D2D8_L6split_16vs13_markers)

D2D8_L7split_19vsOth_markers <- FindMarkers(D2D8, ident.1 = "19", ident.2 = c("0","14"), min.pct = 0.01, only.pos = F, verbose = T)
D2D8_L7split_19vsOth_markers$cluster <- "L7c19"
D2D8_L7split_19vsOth_markers$gene <- rownames(D2D8_L7split_19vsOth_markers)

D2D8_lineageSwitch_DEGs <-
rbind(D2D8_L1split_18vs9_markers,D2D8_L3split_9vs18_markers,D2D8_L2split_15vs11_markers,D2D8_L7split_19vsOth_markers,D2D8_L4spli
t_13vs16_markers,D2D8_L6split_16vs13_markers,D2D8_L5split_1vs4_markers)
D2D8_lineageSwitch_DEGs$p_val_adj = p.adjust(D2D8_lineageSwitch_DEGs$p_val, method='BH')
D2D8_lineageSwitch_DEGs$delta <- abs(as.double(D2D8_lineageSwitch_DEGs$pct.1)-as.double(D2D8_lineageSwitch_DEGs$pct.2))
D2D8_lineageSwitch_DEGs$minor <- ifelse(D2D8_lineageSwitch_DEGs$pct.1 < 0.10 | D2D8_lineageSwitch_DEGs$pct.2 < 0.10, 1, 0)
D2D8_lineageSwitch_DEGs <- D2D8_lineageSwitch_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
109

D2D8_lineageSwitch_DEGs <- D2D8_lineageSwitch_DEGs[D2D8_lineageSwitch_DEGs$p_val_adj<0.05,]
#write.table(D2D8_lineageSwitch_DEGs, "/data/scRNA/IMR90_timecourse/D2D8_lineageSwitch_DEGs.tsv", quote = F, sep = "\t", row.names =
F)

#upregulated, sorted by log2FC
t1 <- D2D8_lineageEnd_DEGs[order(D2D8_lineageEnd_DEGs[,1], -D2D8_lineageEnd_DEGs[,3]), ]
write.table(t1, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D2D8_LinEnd_upregulated_FCsorted.tsv", quote = F, sep = "\t", row.names =
F)
#downregulated, sorted by log2FC
t2 <- D2D8_lineageEnd_DEGs[order(D2D8_lineageEnd_DEGs[,1], D2D8_lineageEnd_DEGs[,3]), ]
write.table(t2, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D2D8_LinEnd_downregulated_FCsorted.tsv", quote = F, sep = "\t", row.names
= F)

#upregulated, sorted by delta
t3 <- subset(D2D8_lineageEnd_DEGs[order(D2D8_lineageEnd_DEGs[,1], -D2D8_lineageEnd_DEGs[,7], -D2D8_lineageEnd_DEGs[,4]), ],
avg_log2FC > 0)
write.table(t3, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D2D8_LinEnd_upregulated_Deltasorted.tsv", quote = F, sep = "\t", row.names
= F)
#downregulated, sorted by delta
t4 <- subset(D2D8_lineageEnd_DEGs[order(D2D8_lineageEnd_DEGs[,1], -D2D8_lineageEnd_DEGs[,7], -D2D8_lineageEnd_DEGs[,4]), ],
avg_log2FC < 0)
write.table(t4, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D2D8_LinEnd_downregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by log2FC
t5 <- D2D8_lineageSwitch_DEGs[order(D2D8_lineageSwitch_DEGs[,1], -D2D8_lineageSwitch_DEGs[,3]), ]
write.table(t5, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D2D8_LinSwitch_upregulated_FCsorted.tsv", quote = F, sep = "\t", row.names
= F)
#downregulated, sorted by log2FC
t6 <- D2D8_lineageSwitch_DEGs[order(D2D8_lineageSwitch_DEGs[,1], D2D8_lineageSwitch_DEGs[,3]), ]
write.table(t6, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D2D8_LinSwitch_downregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by delta
t7 <- subset(D2D8_lineageSwitch_DEGs[order(D2D8_lineageSwitch_DEGs[,1], -D2D8_lineageSwitch_DEGs[,7], -D2D8_lineageSwitch_DEGs[,4]),
], avg_log2FC > 0)
write.table(t7, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D2D8_LinSwitch_upregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by delta
t8 <- subset(D2D8_lineageSwitch_DEGs[order(D2D8_lineageSwitch_DEGs[,1], -D2D8_lineageSwitch_DEGs[,7], -D2D8_lineageSwitch_DEGs[,4]),
], avg_log2FC < 0)
write.table(t8, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D2D8_LinSwitch_downregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)

#

D8D12s1 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12s1_singletsOnly_prepSCTfindmarkers.RDS.gz")
curves1 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s1_curves_singlets.RDS.gz")

DimPlot(D8D12s1, group.by = "seurat_clusters", label = T, label.box = T)
Idents(D8D12s1) <- D8D12s1$seurat_clusters
DefaultAssay(D8D12s1) <- "SCT"

#Lineage1: 18 13 4 16 25 12 17 8
#Lineage2: 18 13 4 16 25 12 15
#Lineage3: 18 13 4 9 11 22
#Lineage4: 18 13 4 2 5

D8D12s1_c8_markers <- FindMarkers(D8D12s1, ident.1 = "8", ident.2 = c("15","22","5"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s1_c8_markers$cluster <- "L1c8"
D8D12s1_c8_markers$gene <- rownames(D8D12s1_c8_markers)

D8D12s1_c15_markers <- FindMarkers(D8D12s1, ident.1 = "15", ident.2 = c("8","22","5"), min.pct = 0.01, only.pos = F, verbose = T)
110

D8D12s1_c15_markers$cluster <- "L2c15"
D8D12s1_c15_markers$gene <- rownames(D8D12s1_c15_markers)

D8D12s1_c22_markers <- FindMarkers(D8D12s1, ident.1 = "22", ident.2 = c("15","8","5"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s1_c22_markers$cluster <- "L3c22"
D8D12s1_c22_markers$gene <- rownames(D8D12s1_c22_markers)

D8D12s1_c5_markers <- FindMarkers(D8D12s1, ident.1 = "5", ident.2 = c("15","22","8"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s1_c5_markers$cluster <- "L4c5"
D8D12s1_c5_markers$gene <- rownames(D8D12s1_c5_markers)

D8D12s1_lineageEnd_DEGs <- rbind(D8D12s1_c8_markers, D8D12s1_c15_markers, D8D12s1_c22_markers, D8D12s1_c5_markers)

D8D12s1_lineageEnd_DEGs$p_val_adj = p.adjust(D8D12s1_lineageEnd_DEGs$p_val, method='BH')
D8D12s1_lineageEnd_DEGs$delta <- abs(as.double(D8D12s1_lineageEnd_DEGs$pct.1)-as.double(D8D12s1_lineageEnd_DEGs$pct.2))
D8D12s1_lineageEnd_DEGs$minor <- ifelse(D8D12s1_lineageEnd_DEGs$pct.1 < 0.10 | D8D12s1_lineageEnd_DEGs$pct.2 < 0.10, 1, 0)
D8D12s1_lineageEnd_DEGs <- D8D12s1_lineageEnd_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
D8D12s1_lineageEnd_DEGs <- D8D12s1_lineageEnd_DEGs[D8D12s1_lineageEnd_DEGs$p_val_adj<0.05,]
#write.table(D8D12s1_lineageEnd_DEGs, "/data/scRNA/IMR90_timecourse/D8D12s1_lineageEnd_DEGs.tsv", quote = F, sep = "\t", row.names
= F)

#Lineage1: 18 13 4 16 25 12 17 8
#Lineage2: 18 13 4 16 25 12 15
#Lineage3: 18 13 4 9 11 22
#Lineage4: 18 13 4 2 5

D8D12s1_L1split_17vs15_markers <- FindMarkers(D8D12s1, ident.1 = "17", ident.2 = c("15"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s1_L1split_17vs15_markers$cluster <- "L1c17"
D8D12s1_L1split_17vs15_markers$gene <- rownames(D8D12s1_L1split_17vs15_markers)

D8D12s1_L2split_15vs17_markers <- FindMarkers(D8D12s1, ident.1 = "15", ident.2 = c("17"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s1_L2split_15vs17_markers$cluster <- "L2c15"
D8D12s1_L2split_15vs17_markers$gene <- rownames(D8D12s1_L2split_15vs17_markers)

D8D12s1_L3split_9vs16_markers <- FindMarkers(D8D12s1, ident.1 = "9", ident.2 = c("16", "2"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s1_L3split_9vs16_markers$cluster <- "L3c9"
D8D12s1_L3split_9vs16_markers$gene <- rownames(D8D12s1_L3split_9vs16_markers)

D8D12s1_L4split_2vs9_markers <- FindMarkers(D8D12s1, ident.1 = "2", ident.2 = c("9", "16"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s1_L4split_2vs9_markers$cluster <- "L4c2"
D8D12s1_L4split_2vs9_markers$gene <- rownames(D8D12s1_L4split_2vs9_markers)

D8D12s1_lineageSwitch_DEGs <-
rbind(D8D12s1_L1split_17vs15_markers,D8D12s1_L2split_15vs17_markers,D8D12s1_L3split_9vs16_markers,D8D12s1_L4split_2vs9_markers)

D8D12s1_lineageSwitch_DEGs$p_val_adj = p.adjust(D8D12s1_lineageSwitch_DEGs$p_val, method='BH')
D8D12s1_lineageSwitch_DEGs$delta <- abs(as.double(D8D12s1_lineageSwitch_DEGs$pct.1)-as.double(D8D12s1_lineageSwitch_DEGs$pct.2))
D8D12s1_lineageSwitch_DEGs$minor <- ifelse(D8D12s1_lineageSwitch_DEGs$pct.1 < 0.10 | D8D12s1_lineageSwitch_DEGs$pct.2 < 0.10, 1, 0)
D8D12s1_lineageSwitch_DEGs <- D8D12s1_lineageSwitch_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta",
"minor")]
D8D12s1_lineageSwitch_DEGs <- D8D12s1_lineageSwitch_DEGs[D8D12s1_lineageSwitch_DEGs$p_val_adj<0.05,]
#write.table(D8D12s1_lineageSwitch_DEGs, "/data/scRNA/IMR90_timecourse/D8D12s1_lineageSwitch_DEGs.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by log2FC
t1 <- D8D12s1_lineageEnd_DEGs[order(D8D12s1_lineageEnd_DEGs[,1], -D8D12s1_lineageEnd_DEGs[,3]), ]
write.table(t1, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s1_LinEnd_upregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by log2FC
t2 <- D8D12s1_lineageEnd_DEGs[order(D8D12s1_lineageEnd_DEGs[,1], D8D12s1_lineageEnd_DEGs[,3]), ]
111

write.table(t2, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s1_LinEnd_downregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by delta
t3 <- subset(D8D12s1_lineageEnd_DEGs[order(D8D12s1_lineageEnd_DEGs[,1], -D8D12s1_lineageEnd_DEGs[,7], -
D8D12s1_lineageEnd_DEGs[,4]), ], avg_log2FC > 0)
write.table(t3, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s1_LinEnd_upregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by delta
t4 <- subset(D8D12s1_lineageEnd_DEGs[order(D8D12s1_lineageEnd_DEGs[,1], -D8D12s1_lineageEnd_DEGs[,7], -
D8D12s1_lineageEnd_DEGs[,4]), ], avg_log2FC < 0)
write.table(t4, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s1_LinEnd_downregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by log2FC
t5 <- D8D12s1_lineageSwitch_DEGs[order(D8D12s1_lineageSwitch_DEGs[,1], -D8D12s1_lineageSwitch_DEGs[,3]), ]
write.table(t5, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s1_LinSwitch_upregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by log2FC
t6 <- D8D12s1_lineageSwitch_DEGs[order(D8D12s1_lineageSwitch_DEGs[,1], D8D12s1_lineageSwitch_DEGs[,3]), ]
write.table(t6, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s1_LinSwitch_downregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by delta
t7 <- subset(D8D12s1_lineageSwitch_DEGs[order(D8D12s1_lineageSwitch_DEGs[,1], -D8D12s1_lineageSwitch_DEGs[,7], -
D8D12s1_lineageSwitch_DEGs[,4]), ], avg_log2FC > 0)
write.table(t7, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s1_LinSwitch_upregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by delta
t8 <- subset(D8D12s1_lineageSwitch_DEGs[order(D8D12s1_lineageSwitch_DEGs[,1], -D8D12s1_lineageSwitch_DEGs[,7], -
D8D12s1_lineageSwitch_DEGs[,4]), ], avg_log2FC < 0)
write.table(t8, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s1_LinSwitch_downregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)

#

D8D12s2 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12s2_singletsOnly_prepSCTfindmarkers.RDS.gz")
curves2 <- readRDS("/data/scRNA/IMR90_timecourse/seu_D8_to_D12_s2_curves_singlets.RDS.gz")

DimPlot(D8D12s2, group.by = "seurat_clusters", label = T, label.box = T)
Idents(D8D12s2) <- D8D12s2$seurat_clusters
DefaultAssay(D8D12s2) <- "SCT"

#Lineage1: 10 6 7 1 24 0 3 14
#Lineage2: 10 6 7 1 24 0 3 21
#Lineage3: 10 6 7 1 19
#Lineage4: 10 6 20

D8D12s2_c14_markers <- FindMarkers(D8D12s2, ident.1 = "14", ident.2 = c("21","19","20"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s2_c14_markers$cluster <- "L1c14"
D8D12s2_c14_markers$gene <- rownames(D8D12s2_c14_markers)

D8D12s2_c21_markers <- FindMarkers(D8D12s2, ident.1 = "21", ident.2 = c("14","19","20"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s2_c21_markers$cluster <- "L2c21"
D8D12s2_c21_markers$gene <- rownames(D8D12s2_c21_markers)

D8D12s2_c19_markers <- FindMarkers(D8D12s2, ident.1 = "19", ident.2 = c("21","14","20"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s2_c19_markers$cluster <- "L3c19"
D8D12s2_c19_markers$gene <- rownames(D8D12s2_c19_markers)

D8D12s2_c20_markers <- FindMarkers(D8D12s2, ident.1 = "20", ident.2 = c("21","19","14"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s2_c20_markers$cluster <- "L4c20"
D8D12s2_c20_markers$gene <- rownames(D8D12s2_c20_markers)

112

D8D12s2_lineageEnd_DEGs <- rbind(D8D12s2_c14_markers,D8D12s2_c21_markers,D8D12s2_c19_markers,D8D12s2_c20_markers)

D8D12s2_lineageEnd_DEGs$p_val_adj = p.adjust(D8D12s2_lineageEnd_DEGs$p_val, method='BH')
D8D12s2_lineageEnd_DEGs$delta <- abs(as.double(D8D12s2_lineageEnd_DEGs$pct.1)-as.double(D8D12s2_lineageEnd_DEGs$pct.2))
D8D12s2_lineageEnd_DEGs$minor <- ifelse(D8D12s2_lineageEnd_DEGs$pct.1 < 0.10 | D8D12s2_lineageEnd_DEGs$pct.2 < 0.10, 1, 0)
D8D12s2_lineageEnd_DEGs <- D8D12s2_lineageEnd_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta", "minor")]
D8D12s2_lineageEnd_DEGs <- D8D12s2_lineageEnd_DEGs[D8D12s2_lineageEnd_DEGs$p_val_adj<0.05,]
#write.table(D8D12s2_lineageEnd_DEGs, "/data/scRNA/IMR90_timecourse/D8D12s2_lineageEnd_DEGs.tsv", quote = F, sep = "\t", row.names
= F)

#Lineage1: 10 6 7 1 24 0 3 14
#Lineage2: 10 6 7 1 24 0 3 21
#Lineage3: 10 6 7 1 19
#Lineage4: 10 6 20

D8D12s2_L1split_14vs21_markers <- FindMarkers(D8D12s2, ident.1 = "14", ident.2 = c("21"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s2_L1split_14vs21_markers$cluster <- "L1c14"
D8D12s2_L1split_14vs21_markers$gene <- rownames(D8D12s2_L1split_14vs21_markers)

D8D12s2_L2split_21vs14_markers <- FindMarkers(D8D12s2, ident.1 = "21", ident.2 = c("14"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s2_L2split_21vs14_markers$cluster <- "L2c21"
D8D12s2_L2split_21vs14_markers$gene <- rownames(D8D12s2_L2split_21vs14_markers)

D8D12s2_L3split_19vs24_markers <- FindMarkers(D8D12s2, ident.1 = "19", ident.2 = c("24"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s2_L3split_19vs24_markers$cluster <- "L3c19"
D8D12s2_L3split_19vs24_markers$gene <- rownames(D8D12s2_L3split_19vs24_markers)

D8D12s2_L4split_20vs7_markers <- FindMarkers(D8D12s2, ident.1 = "20", ident.2 = c("7"), min.pct = 0.01, only.pos = F, verbose = T)
D8D12s2_L4split_20vs7_markers$cluster <- "L4c20"
D8D12s2_L4split_20vs7_markers$gene <- rownames(D8D12s2_L4split_20vs7_markers)

D8D12s2_lineageSwitch_DEGs <-
rbind(D8D12s2_L1split_14vs21_markers,D8D12s2_L2split_21vs14_markers,D8D12s2_L3split_19vs24_markers,D8D12s2_L4split_20vs7_marker
s)

D8D12s2_lineageSwitch_DEGs$p_val_adj = p.adjust(D8D12s2_lineageSwitch_DEGs$p_val, method='BH')
D8D12s2_lineageSwitch_DEGs$delta <- abs(as.double(D8D12s2_lineageSwitch_DEGs$pct.1)-as.double(D8D12s2_lineageSwitch_DEGs$pct.2))
D8D12s2_lineageSwitch_DEGs$minor <- ifelse(D8D12s2_lineageSwitch_DEGs$pct.1 < 0.10 | D8D12s2_lineageSwitch_DEGs$pct.2 < 0.10, 1, 0)
D8D12s2_lineageSwitch_DEGs <- D8D12s2_lineageSwitch_DEGs[,c("cluster","gene","avg_log2FC", "p_val_adj", "pct.1","pct.2","delta",
"minor")]
D8D12s2_lineageSwitch_DEGs <- D8D12s2_lineageSwitch_DEGs[D8D12s2_lineageSwitch_DEGs$p_val_adj<0.05,]
#write.table(D8D12s2_lineageSwitch_DEGs, "/data/scRNA/IMR90_timecourse/D8D12s2_lineageSwitch_DEGs.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by log2FC
t1 <- D8D12s2_lineageEnd_DEGs[order(D8D12s2_lineageEnd_DEGs[,1], -D8D12s2_lineageEnd_DEGs[,3]), ]
write.table(t1, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s2_LinEnd_upregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by log2FC
t2 <- D8D12s2_lineageEnd_DEGs[order(D8D12s2_lineageEnd_DEGs[,1], D8D12s2_lineageEnd_DEGs[,3]), ]
write.table(t2, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s2_LinEnd_downregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by delta
t3 <- subset(D8D12s2_lineageEnd_DEGs[order(D8D12s2_lineageEnd_DEGs[,1], -D8D12s2_lineageEnd_DEGs[,7], -
D8D12s2_lineageEnd_DEGs[,4]), ], avg_log2FC > 0)
write.table(t3, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s2_LinEnd_upregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by delta
t4 <- subset(D8D12s2_lineageEnd_DEGs[order(D8D12s2_lineageEnd_DEGs[,1], -D8D12s2_lineageEnd_DEGs[,7], -
D8D12s2_lineageEnd_DEGs[,4]), ], avg_log2FC < 0)
113

write.table(t4, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s2_LinEnd_downregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by log2FC
t5 <- D8D12s2_lineageSwitch_DEGs[order(D8D12s2_lineageSwitch_DEGs[,1], -D8D12s2_lineageSwitch_DEGs[,3]), ]
write.table(t5, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s2_LinSwitch_upregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by log2FC
t6 <- D8D12s2_lineageSwitch_DEGs[order(D8D12s2_lineageSwitch_DEGs[,1], D8D12s2_lineageSwitch_DEGs[,3]), ]
write.table(t6, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s2_LinSwitch_downregulated_FCsorted.tsv", quote = F, sep = "\t",
row.names = F)

#upregulated, sorted by delta
t7 <- subset(D8D12s2_lineageSwitch_DEGs[order(D8D12s2_lineageSwitch_DEGs[,1], -D8D12s2_lineageSwitch_DEGs[,7], -
D8D12s2_lineageSwitch_DEGs[,4]), ], avg_log2FC > 0)
write.table(t7, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s2_LinSwitch_upregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)
#downregulated, sorted by delta
t8 <- subset(D8D12s2_lineageSwitch_DEGs[order(D8D12s2_lineageSwitch_DEGs[,1], -D8D12s2_lineageSwitch_DEGs[,7], -
D8D12s2_lineageSwitch_DEGs[,4]), ], avg_log2FC < 0)
write.table(t8, "/data/scRNA/IMR90_timecourse/DEGs_Lineages/D8D12s2_LinSwitch_downregulated_Deltasorted.tsv", quote = F, sep = "\t",
row.names = F)

Abstract (if available)

Abstract As one of the main hallmarks of aging, SnCs (SnCs) play a complex role in the aging process. These permanently grow-arrested cells prevent can suppress the development of cancer; however they also contribute to age-related tissue dysfunction and chronic inflammation, which are both associated with many age-related diseases. SnCs are typically studied as endpoints of a complex transformational process, owing to their frequent maladaptive effects on surrounding tissue and cells. SnCs accumulate with age, and while they ultimately comprise a small percentage of cells in tissues, they have important roles in age associated pathologies.
To expand our temporal understanding of senescence and elucidate the developmental trajectories of SnCs, we performed single-cell RNA sequencing on IMR90 lung fibroblasts senescencing across a 12-day period. Our analysis revealed substantial heterogeneity in gene expression within timepoints and across the full time-course. Supervised trajectory inference of the time-course data uncovered the root-origin and fates of distinct SnC lineages over 3 stages of senescence development. We found that lineages had distinct regulatory and pathway activation profiles. These findings confirm the importance of deeply profiling the nature of SnCs in the relevant biological context of interest, to avoid the paired threats of poor efficiency and off target effects. Altogether our data provide a novel approach to study SnC development, identifying SnC subtypes (“senotypes”), and differentiating between SnCs and quiescent cells. Applying this approach to known senescence contexts of interest will aid in identifying key targets and areas for therapeutic intervention in pathologies of aging caused by SnCs.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Investigating brain aging and neurodegenerative diseases through omics data

PDF

Natural divergence of traits in species of mice reveals novel molecular mechanisms of cellular senescence

PDF

Identifying and measuring cell-intrinsic and cell-extrinsic factors influencing aging

PDF

Longitudinal assessment of neural stem-cell aging

PDF

Characterization of senescent cell heterogeneity using cell culture models

PDF

Metabolomic and proteomic approaches to understanding senescence and aging in mammary epithelial cells

PDF

Peripheral blood mononuclear cell capture and sequencing: optimization of a droplet based capture method and its applications

PDF

The overlap between mTOR signaling, rapamycin and cellular senescence

PDF

Fibroblastic connective tissue cells: the blastema stem cells and source of large-scale chondrogenesis in the regenerating lizard tail

PDF

The impact of genetic pleiotropy on heterogeneity in the developing forebrain and hindbrain

PDF

Single cell gene expression analysis of KSHV infection in three-dimensional oral epithelial tissue cultures

PDF

Modeling neurodegenerative diseases using induced pluripotent stem cells and identifying therapeutic targets

PDF

Transcriptomic maturation of developing human cone precursors in fetal and 3D hESC-derived tissues

PDF

Determining the epigenetic contribution of basal cell identity in cystic fibrosis

PDF

Investigating hematopoietic stem cells via multiscale modeling, single-cell genomics, and gene regulatory network inference

PDF

Data-driven learning for dynamical systems in biology

PDF

Deconvolution of circulating tumor cell heterogeneity and implications for aggressive variant prostate cancer

PDF

Transcription errors induce proteotoxic stress in mammalian cells

PDF

Pten deletion in adult pancreatic beta-cells induces cell proliferation and G1/S cell cycle progression

PDF

A generalization of the accumulation curve in species sampling and its applications to high‐throughput sequencing

Asset Metadata

Creator Ciotlos, Serban (author)

Core Title A single cell time course of senescence uncovers discrete cell trajectories and transcriptional heterogeneity

School Leonard Davis School of Gerontology

Degree Doctor of Philosophy

Degree Program Biology of Aging

Degree Conferral Date 2023-05

Publication Date 05/05/2023

Defense Date 03/29/2023

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag aging,bioinformatics,cellular senescence,Computational Biology,fibroblast,IMR90,longevity,OAI-PMH Harvest,proliferation,pseudotime,quiescence,scRNAseq,senescence,single cell profiling,single cell RNA sequencing,trajectory,trajectory inference

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Melov, Simon (committee chair), Benayoun, Berenice (committee member), Campisi, Judith (committee member), Ellerby, Lisa (committee member), Tower, John (committee member)

Creator Email serbanciotlos@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC113099918

Unique identifier UC113099918

Identifier etd-CiotlosSer-11791.pdf (filename)

Legacy Identifier etd-CiotlosSer-11791

Document Type Dissertation

Format theses (aat)

Rights Ciotlos, Serban

Internet Media Type application/pdf

Type texts

Source 20230505-usctheses-batch-1038 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu