Runx factors launch T-cell and innate lymphoid programs via direct and gene network-based mechanisms

Runx factors launch T-cell and innate

lymphoid programs via direct and

gene network-based mechanisms

Boyoung Shin

, Wen Zhou

1,2,3

, Jue Wang

1,2

, Fan Gao

1,4

, Ellen V. Rothenberg

Division of Biology and Biological Engineering, California Instit

ute of Technology.

Pasadena, CA 91125, USA

Program in Biochemistry and Molecular Biophysics,

California Institut

e of Technology.

Current address: BillionToOne, Menlo Park, CA

Bioinformatics Resource Center, Beckman Instit

ute of California Institute of Technology.

*Corresponding author.

Email: evroth@its.caltech.edu

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

Abstract

Runx factors are essential for lineage spec

ification of various

hematopoietic cells,

including T lymphocytes. However, they

regulate context-specific genes and occupy

distinct genomic regions in different cell types

. Here, we show that dynamic Runx binding

shifts in early T-cell developm

ent are mostly not restricted

by local chromatin state but

regulated by Runx dosage and functi

onal partners. Runx co-fac

tors compete to recruit a

limited pool of Runx factors

in early T-progenitors, and a modest increase in Runx protein

availability at pre-commitm

ent stages causes prematur

e Runx occupancy at post-

commitment binding sites. This results in

striking T-lineage devel

opmental acceleration

by selectively activating T-identity and

innate lymphoid cell programs. These are

collectively regulated by Runx

together with other, Runx-induced

transcription factors that

co-occupy Runx target genes and propagate gene network changes.

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

Introduction

Runx family transcription fa

ctors (Runx1, Runx2, Runx

3, and their cofactor CBF

)

are essential for T cell development from the

earliest steps in the lineage, playing partially

redundant roles

1-6

. However, they are also vital fo

r the establishment of hematopoietic

stem cells in early embryos

and for generation of B cells

and megakaryocytes throughout

life

8-11

, quite different programs.

Early intrathymic T cell development itself spans distinct

regulatory contexts where

Runx factors can implement

stage-unique or continuing

functional inputs at different stag

es. Runx target motifs are consistently highly enriched

around open chromatin sites and lin

eage-specific transcription fa

ctor (TF) binding sites in

multiple hematopoietic lineages

12-19

, suggesting a common contribution to active

enhancers generally. However, we have recent

ly shown that Runx factors occupy

different genomic sites at different stages of ear

ly T cell development, and that this results

in their regulation of

different target genes

. Thus, key questions

are how Runx factors

can accurately guide their contributions to di

stinct cell programs,

whether by intrinsic

DNA-binding sequence specificit

y, epigenetic constraints, or

interactions with other

partner factors. Here, we show how levels

of expression of Runx itself control the

qualitative choices of site occupancy to c

ontrol T cell developmental speed and pathway

choice.

Stages in T cell development are distinguished by changes in chromatin states

and changes in expression of a di

screte set of regulatory factors

20-26

, even though Runx

factors themselves are collectively

active at similar levels throughout

. Driven by Notch

signaling and thymic microenvir

onmental cues, multipotent progenitor cells are converted

to T-lineage committed pro-T cells withi

n the thymus. They pass through CD4

CD8

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

double negative (DN) subs

tages (DN1-4) to CD4

CD8

double-positive (DP) stage before

becoming mature CD4 or CD8 si

ngle-positive (SP) T cells (Fig. 1a). Pro-T cells in DN1

and DN2a stages (“Phase 1”) still resemble

hematopoietic stem and progenitor cells

(HSPC) in regulatory gene expression and chromatin state and can still produce non-T

lineage cells. Definitive T-

lineage commitment normally occurs

in transition from DN2a

to DN2b stages. Up-regulation of T cell identi

ty genes and changes in TF expression and

chromatin states

20, 21, 23

during commitment (DN2b) de

fine “Phase 2”, extending till

successful assembly of T cell receptor

(TCR

) in DN3 stage (Fig. 1a). Runx1 and

Runx3 are crucial for progression through bot

h Phases. Note that although Runx1 and

Runx3 act differently in other contexts

, they appear functionally redundant in pro-T cells

However, they both bind to different genomic

sites and regulate different target genes

from Phase 1 to Phase 2

Profound changes occur in chromatin

looping, accessibility, and histone

modification profiles during T-lineage commitment

22-24

, associated with repression of the

genes associated with progenitor and alternative lineage states and activation of the T

cell identity programs

20, 21, 28, 29

. Multiple TFs also change activity

21, 30, 31

. Runx TFs

themselves can interact physically with multiple

TFs at distinct binding sites, suggesting

that TF cooperativity could be a

major influence on Runx activity

13, 32, 33

. Yet how

important are the Runx factors

themselves in dictating which factor complexes and which

target sites will be active?

Here, we evaluated the chromatin constr

aints on Runx action across the Phase 1

and Phase 2 stages of T cell development and

tested the hypothesis that Runx binding

site shifts depend on competition between Phas

e 1 and Phase 2 partners for a limiting

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

amount of Runx protein. We

found that at modest excess,

when no longer titrated by

Phase 1 partners, Runx factors directed a di

stinctive accelerated form of early T and

innate-like cell development, relieving the need to repress most Phase 1 regulators before

T-lineage regulatory genes could

be upregulated. Both direct (binding site-mediated) and

indirect mechanisms propagated through a

Runx-dependent gene r

egulatory network

drove this acceleration. Thus, Runx factor levels are major timing controllers of the

deployment of the T-cell specif

ication gene regulatory network.

Results

Runx TFs substantially shifted their binding sites during all

stages of T-cell

development

Runx 1 and Runx 3 are functionally redundant

and bind to similar sites in early pro-

T cells, but their site c

hoices are highly stage dependent

. Fig. 1b shows that this not

only distinguishes pro-T cell stages but also

extends to Runx deploy

ment in different

hematopoietic lineage contexts from HSPCs to mature T cells, B lineage and

megakaryocyte-precursor cells

15, 34-38

. Distinct genomic regions showed cell-type specific

Runx occupancies, defining separate regions

preferentially occupied only in HSPC, in B

cell progenitors, in DN1, in DN3,

in DP, or in mature T cells (Fig. 1b, A-F), and regions

occupied in different developmen

tal combinations, with ~12% of

sites shared in all (Fig.

1b, G). This site infidelity of Runx factors contrasted with binding pa

tterns in pro-T cells

of PU.1, a critical Phase 1-specific TF

inherited from bone-marrow progenitor cells, which

showed more similar binding site choices fr

om HSPCs to DN2b pro-T cells (Fig. S1a).

Importantly, each cluster of Runx binding regi

ons from HSPC to mature T cells harbored

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

a distinct set of motifs in addi

tion to the Runx motif, in which motifs for EBF, PU.1, E2A,

TCF1, GATA, JunB, and ETS factors were diffe

rentially enriched in each cluster (Fig.

100

S1b). In our previous r

eport, Runx binding was monitored by ChIP-seq after

101

disuccinimidyl glutarate-assi

sted stabilization crosslinking

, which might have biased our

102

previous results to overrepresent Runx comp

lexes with other protei

ns rather than the

103

distribution of Runx binding preferences themselves. Here, we independently analyzed

104

Runx DNA binding profiles in pro T cells

(DN1 and DN2b/DN3) using cross-linking-

105

independent CUT

&RUN instead

39, 40

(C&R, Cleavage Under Targets & Release Using

106

Nuclease)(see Methods;

Fig. S1c-f shows detailed com

parison of C&R against previous

107

ChIP-seq identified sites). These results suggested that even ex

cluding crosslinking

108

artifacts, Runx factors r

eadily changed their binding sites across all stages of T cell

109

development to interact with distinct cell-type

specific regions which may be occupied by

110

different TF ensembles.

111

112

Runx factors prefer “active” chromati

n compartments but do not follow local

113

chromatin state changes

114

We previously showed that some dynamic

Runx binding shifts during T-lineage

115

commitment occurred in a highly coor

dinated manner, with group appearance or

116

disappearance of multiple Runx occupanc

ies across large genomic domains of 10

-10

117

kb (ref.

). To test whether Runx

TFs were constrained or redirected by large-scale

118

chromatin remodeling during

commitment, the non-promoter

Runx binding sites were

119

categorized into three groups: Phase 1-pr

eferential binding sites (Group 1), Phase 2-

120

preferential binding sites (Group 2), and stably occupied sites

(Group 3) (Fig. 1c). We

121

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

analyzed how Runx binding was correlated with “act

ive” (A) or “inactive” (B) large-scale

122

nucleome compartments

by comparing the principal co

mponent 1 (PC1) values of the

123

previously reported Hi-C correlation matr

ices from ETP (DN1), DN2, and DN3 cells

. All

124

Runx binding sites were preferentially enr

iched within the A comp

artment (84-92%) and

125

were almost depleted from t

he B compartment (3-7%), regardl

ess of Group (Fig. S2a, b;

126

note

Ets1

flanking regions). As

pro-T cells developed from

ETP (DN1) to DN3 stages,

127

most genomic regions remained in the same

compartment (“A-to-

A” 41.9%, “B-to-B”

128

48.1%). Among the minority of genomic regions changing co

mpartment, Runx occupancy

129

tended to follow the active states (Fig. S2a)

. Compartments switch

ing from active to

130

inactive (decreasing PC1 values, A to B tr

end) included more Group 1 (4.75%) than Group

131

2 sites (1.32%), whereas com

partments becoming active (inc

reasing PC1 values, B to A

132

trend) included more Group 2 sites (4.56%) than

Group 1 (2.35%) (Fig. S2a). One striking

133

example of concerted Group 2 site appearance

with a B to A compartm

ent flip was seen

134

in the extended flanking region of

Bcl11b

(Fig. S2b). However, most Runx site shifts

135

occurred within A compartments.

136

Local chromatin states at numerous sites change substantially during the transition

137

from Phase 1 to Phase 2, based on nucleos

ome density [assay of transposase accessible

138

chromatin (ATAC), or DNase acce

ssibility] and histone modifications

22-24

. To test how

139

local chromatin state changes associate with

Runx factor redist

ribution, we coded

140

individual chromatin states across the genom

e from pre-commitment to post-commitment

141

stages using ChromHMM

42, 43

with published datasets for chro

matin state marks in pro-T

142

cells

23, 24, 44

(see Methods). Globally, Runx binding sites were not enriched in genomic

143

regions with repression-associ

ated chromatin states, whether

defined by high levels of

144

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

H3K27me3 or without any active chromatin ma

rks (state 15, 16). Runx binding sites

145

overall were preferentially enriched in open/active chromatin regions (state 1-3, 5-10) or

146

weakly accessible regions harboring H3K4m

e2 marks representing

likely poised regions

147

(state 13)(Fig. S2c). This global bias was

notable because we had verified that C&R could

148

detect Runx binding in closed chromatin at

least as well as ChIP-seq (Fig. S1f), and

149

because Runx factors can work both as repressors and as activators

45-47

150

However, the developmental changes of R

unx binding patterns did not strictly

151

follow developmental changes in local chromati

n states (Fig. S2c). At genomic sites

152

occupied stage-specifically by Runx fact

ors in Phase 1 or Phase 2 (Fig. 1d),

153

developmental shifts in Runx occupancies

could occur without corresponding changes in

154

accessibility of those sites. Of Group 1 (Phas

e 1-specific) sites, only 43.8% were open

155

in a Phase 1-restricted way, and only 21.4%

of the Group 2 (Phase 2-

specific) sites were

156

open selectively in Phase 2. T

herefore, over 50% of Grou

p 1 and about 80% of Group 2

157

sites failed to change local chromatin accessibi

lity as Runx binding changed in the Phase

158

1-Phase 2 transition.

For instance, near

Meis1

multiple Runx occupancies disappeared

159

from DN1 to DN3 (Group 1 peaks), but t

hese sites remained open by ATAC-seq.

160

Conversely, at

Ets1

multiple sites gained Runx occupa

ncies from DN1 to DN3, but these

161

sites had been accessible from DN1 (Fig. 1e, Fig. S2b). Over 1/3 of Group 2 sites

162

remained closed in both Phases (Fig. 1d). T

hus, local chromatin state failed to explain

163

why Runx occupancy was delayed

at Group 2 bindi

ng sites.

164

165

Sequence specific features and partner fact

or interactions distinguish chromatin

166

sites with developmentally

changing Runx occupancies

167

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

Two other possible explanations for site choi

ce shifts could be differences in Runx

168

binding avidity (site affinity times site dens

ity) which could make Group 2 sites highly

169

sensitive to small changes in Runx concentra

tion, or collaborations

with different stage-

170

specific partners

. We evaluated these options by quantitative motif analysis, focusing

171

exclusively on Runx sites that were consistent

ly “open” to minimize

chromatin effects (Fig.

172

1f, g). Runx binding sites mapping to open pr

omoter regions had negl

igible Runx motif

173

frequencies and poor Runx motif quality scores (Fig. 1f). Stably open non-promoter Runx

174

binding regions in Groups 1, 2, and 3

had much higher Runx motif frequencies and motif

175

quality than promoter si

tes, but to different degrees. Consistently occupied Group 3 sites

176

had the highest scores. Both Group 1 and Gr

oup 2 sites showed lower Runx motif

177

frequencies and motif scores than the Group 3 sites, but were

similar to each other. Thus,

178

at non-promoter sites

without chromatin barrier

s, stage-specific redistribution of Runx

179

factors occurred most read

ily between “modest” Runx motif sites without strong

180

advantages for recruiting Runx

factors via DNA recognition

per se

(Fig. 1f).

181

We previously identified dist

inct partner factors for Ru

nx binding in Phase 1 and

182

Phase 2

13, 32

, and found distinct partner

motifs enriched at Runx

sites in different stages

183

(Fig. 1g, Fig. S1b).

De novo

motif enrichment analysis of the open sites confirmed that

184

the Group 1 sites were highly enriched for PU

.1 (ETS subfamily)

motifs whereas the

185

Group 2 sites were highly enriched for E2A (bas

ic helix-loop-helix, bHLH) motifs. Although

186

C&R Runx peaks did not show the extreme

enrichment of ETS motifs seen with ChIP-

187

seq (Fig. S1e), canonical (non-PU.1) ETS fact

or motifs were still enriched, and were

188

found at similar frequencies in all classes of non-

promoter sites (Fig. 1g). Thus, at sites

189

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

that were stably accessible throughout Phases

1 and 2, different en

sembles of TFs might

190

recruit Runx TFs stage-specifically.

191

A central question is whether the T cell

commitment process is switchlike, e.g.

192

whether a single mechanism causes Runx factor

s to shift from Group

1 vs. Group 2 sites.

193

The majority of precommitment-specific bind

ing sites for Runx fa

ctors (Group 1 sites)

194

have been shown to be actual co-binding sites with PU.1

1, 13

. Besides PU.1, in later pro-

195

T cells other TFs such as GATA3 and Bcl11b c

an also collaborate with Runx factors at

196

different sites

32, 33

. Notably, the presence of PU.1 c

an redirect Runx1 occupancy to the

197

PU.1 sites, while depleting Runx1 binding

(“theft”) from alte

rnative Runx sites

. If Runx

198

factor levels are truly limit

ed such that partners have to

compete to recruit Runx to

199

different sites, the tipping of a balance between partners might cause concerted

200

occupancy switches from Gr

oup 1 to Group 2 sites.

201

We hypothesized that if such competition

occurs, it could be overridden if Runx

202

availability were increased. We first tested this hypothesis in a PU.1 “theft” model (Fig.

203

S3). The DN3-like Scid.adh.2C

2 pro-T cell line, represen

ting a Phase 2 state, was

204

retrovirally transduced with exogenous PU.1,

with or without additional exogenous Runx1

205

(Fig. S3a-b). PU.1 activated myeloid mark

ers in the cells with

or without exogenous

206

Runx1 (Fig. S3c) and recruited endogenous Ru

nx1 to a set of new co-occupancy sites

207

with PU.1, most of which had been closed before

(Fig. S3d, PU.1-induced). As previously

208

reported

, without extra Runx1, PU.1 also caused

a loss of Runx1 occupancy from nearly

209

55 % of the normal endogenous Ru

nx binding sites (Fig. S3d,

PU.1-depleted). However,

210

when extra Runx1 was added (OE), although PU.1

was still able to recruit Runx binding

211

to the PU.1-induced sites, occupancy of the

PU.1-depleted sites was fully rescued (Fig.

212

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

S3d). The extra Runx1 also occupied a se

t of novel sites (OE new). These had high

213

quality Runx motifs at high frequency (Fig. S3e)

, but were mostly sequestered in closed

214

chromatin in the normal Scid.adh.2C2. These

results suggest that

either Runx1-PU.1

215

complexes or high-level Runx1 alone could gain

access to normally closed chromatin, but

216

that the ability of PU.1 to remove Runx1

from its default binding sites was based on

217

competitive titration when

Runx1 was limiting.

218

219

Modestly increased Runx1 in Phase 1 pr

ematurely induced developmentally

220

important TFs

221

If titration of potentially competing partner factors affect

s Runx site choice, Runx

222

concentration might affect the T-lineage specific

ation program in early progenitor cells.

223

To test this, we exploited t

he OP9-Delta-like ligand 1 (Dll1)

in vitro

differentiation system

224

as in our previous studies

. Briefly, bone-marrow derived

progenitor cells expressing a

225

Bcl2

transgene and

Bcl11b

-mCitrine reporter were co-cultured with OP9-Dll1 cells and

226

exogenous Runx1 was retrovira

lly delivered to pro-T cells when the progenitor cells were

227

still at DN1 stage (Fig. 2a). Then, we me

asured T-development markers (cKit, CD44,

228

and CD25) and

Bcl11b

-mCitrine expression, normally a

marker for T-lineage commitment

229

(see Fig. 1a, 2a)

. At day 2 after exogenous Runx1 in

troduction (overexpression, OE),

230

total Runx1 protein levels were increased by 2-

3 fold relative to the control conditions,

231

and this increase was stable at day 4 post-infe

ction (Fig. 2b, S3f). Notably, the modest

232

degree of increase was important

for the health of the cells

233

Runx1 OE caused a striking acceleration of

Bcl11b

induction as early as day 2

234

after introducing extra Runx1,

increasing at day 3: ~20% of

control cells but ~50% of

235

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

Runx1 OE DN2 cells showed

Bcl11b

-mCitrine expression (Fig. 2c, d). Abnormally,

236

Bcl11b

-mCitrine was also activated in

some Runx1 OE DN1 cells (CD25

). Furthermore,

237

increased Runx1 levels caus

ed premature appearan

ce of cells resembling DN3 cells

238

(CD44

low

CD25

)(Fig. 2c, d).

239

We examined expression profiles of

developmentally important TFs, TCF1,

240

GATA3, and PU.1. Runx1 OE upregulated TC

F1 and GATA3 protein expression within

241

the cKit

CD25

DN1 (ETP) population at both days 2 and 4 post-infection (Fig. 2e). TCF1

242

and GATA3 levels in control cells only reached

those of the Runx1 OE populations at day

243

4, in the cells that had turned on

CD25 (DN2a cells)(Fig. 2e). PU.1 (

Spi1

) is normally

244

repressed by Runx1 in Phase 2 only

1, 33

. Its expression was not affected by Runx1

245

overexpression at day 2 post-in

fection, but was significantly

downregulated even in the

246

DN1 population at day 4 (Fig. 2e), to

levels lower than in normal CD25

(DN2a) cells.

247

Hence, a mild increase in Runx factor avail

ability in Phase 1 pro-T cells could accelerate

248

aspects of early T-cell development, especially within cKit

CD25

DN1 cells.

249

250

Runx level perturbations in Phase 1 resu

lted in profound changes in single-cell

251

transcriptomes

252

We tested critically the ability of Runx fact

or levels to affect the T-cell specification

253

program as a whole, using single-cell RNA-seq

(scRNA-seq). To identify targets of Runx1

254

OE that were also dependent on normal Runx

levels in controls, we also compared

255

Runx1/Runx3

double knockout (dKO) cells

, measuring the single-

cell transcriptomes of

256

Runx1 OE, control, and

Runx1/Runx3

dKO cells together using the 10X Chromium

257

system. We delivered Runx1-OE vector

or empty-vector

control into

Bcl2

-transgenic

258

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

progenitor cells to test OE, or guide-RNAs (gRNAs) against

Runx1

and

Runx3

or control

259

irrelevant gRNAs into

Cas9;Bcl2

transgenic prethymic progenitor

cells to test dKO. Then

260

these progenitor cells were each co-cultu

red with OP9-Dll1 cells

to two different

261

timepoints (day 2 and day 4 post-infection for

OE; day 3 and day 6 post-infection for dKO)

262

and marked by hashtag oligos before they were pooled and subjected to single cell RNA-

263

seq together (scRNA-seq) (Fig. 3a).

264

From two independent 10X r

uns, we recovered 15,310 cells that successfully

265

passed the quality control criter

ia (see Methods). In a lo

w dimensional transcriptome

266

representation after batch correction, the fi

rst parameters in t-distributed stochastic

267

neighbor embedding (tSNE) and Uniform Manifold

Approximation and Projection (UMAP)

268

reflected cell-cycle phases (Fig

. S4a, b, top panels). Afte

r cell-cycle regression, UMAP2

269

(x axes, Figs. 3b-d) approxim

ately represented real time dev

elopmental progression for

270

normal pro-T cells, while UMAP3 (y axes) reflected perturbation; not

e that cells in each

271

population progressed asynchronously. In c

ontrols, cells with low-UMAP2 values

272

expressed high levels of

DN1 signature genes (

Lmo2, Spi1, Bcl11a, Cd34, Mef2c, Hhex

273

Genes transiently expressed dur

ing DN1 to DN2a transition (

Mycn, Fgf3

) were maximally

274

expressed in control cells at UMAP2-interm

ediate values (Fig. 3b), while DN2 marker

275

genes (

Il2ra, Tcrg-C1, Gata3, Tcf7, Thy1, Cd3g

) were first expressed in control cells at

276

UMAP2-intermediate values and were ma

intained throughout the UMAP2-high cells.

277

Then, mostly at later timepoints, genes

associated with T-lineage commitment and the

278

DN2b stage (

Bcl11b, Ly6d, Lef1, Ets1

) initiated expression in the UMAP2-high control

279

cells (Fig. 3b). Thus, for controls, UMAP2 posit

ions could relate cell states to the normal

280

developmental progression.

281

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

282

Developmental pathway kinetics were sen

sitive to the Runx dosage changes

283

Both Runx1 OE and

Runx1/Runx3

dKO (“Runx dKO”) caused contrasting

284

deviations from the control cell

clusters in the UMAP3 dimensio

n (y axes in Figs. 3b-d).

285

Control cells from all timepoints were conc

entrated at the center of UMAP3, whereas

286

Runx1 OE cells were shifted to lower UMAP

3 values while Runx dKO cells veered to

287

UMAP3-higher values (Fig. 3c). Consistently, Runx perturbation caused cells to form

288

unique clusters (Louvain clustering, Fig. 3d,

S4c), suggesting that Ru

nx factors regulated

289

pathways followed by individual cells rat

her than changing subpopul

ation distributions

290

along the control pathway.

291

Runx1 OE and Runx dKO cells showed eviden

ce for reciprocally shifted landmark

292

gene expression patterns along the UMAP2 axis

as compared to the control group (Fig.

293

3b). Consistent with

faster induction of

Bcl11b

-mCitrine reporter expr

ession in absolute

294

time, Runx1 OE cells upregulated

Bcl11b

at a lower UMAP2 value than control cells. In

295

addition, Runx1 OE cells upregulated various

later-stage genes “prematurely” at lower

296

UMAP2 values and to higher

levels than controls (

Gata3, Tcf7, Cd3g, Tcf12, Ly6d, Lef1

297

Ets1

). However, not all DN2-associated genes

were concurrently induced (e.g.,

Il2ra,

298

Thy1, Tcrg-C1

), nor were all critical DN1 land

mark genes downregulated (e.g., not

Spi1

)

299

in Runx1 OE cells. On the other hand, Runx dKO caused lingering expression of DN1-

300

associated genes (

Lmo2, Spi1, Bcl11a, Cd34, Mef2c

) with markedly impaired

301

upregulation of later stage genes (

Mycn, Fgf3, Il2ra, Tcrg-C1,

Gata3, Tcf7, Thy1, Cd3g,

302

Bcl11b, Ly6d

). Instead, Runx dKO

cells expressed genes associated with non-T cells,

303

such as

Id2,

Cd81

Csf2rb,

and

Ifngr2

, and prolonged expressi

on of HSPC gene

Meis1

304

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

(Fig. 3b, Fig. S4d). Thus

, many developmentally regul

ated genes sensitively responded

305

to perturbations of Runx levels, although

alternative programs were not coherently

306

activated or inhibited together. These results

also indicated that Runx dosage responses

307

occurred within pro-T cells t

hemselves, not only reflecting

enrichment of minority

308

contaminants.

309

310

The common set of genes sensitive to bot

h gain and loss of Runx functions

311

overlapped with essential gen

es for early T cell development

312

We reasoned that the core set of devel

opment-regulating genes that were most

313

likely to be direct Runx targets should be

reciprocally affected by Runx1 OE and

314

Runx1/Runx3

dKO. Differentially expressed genes

(DEGs) affected by Runx1 OE and

315

Runx1/ Runx3

dKO were each similar at both perturb

ation timepoints (Fig. S4e). The

316

global gene expression changes

mediated by Runx dKO vs.

Runx1 OE were inversely

317

correlated at both timepoints, with a core set

making significant, reciprocal responses to

318

both gain- and loss-of Runx functi

ons (Fig.3e, f and S4f, Table S1, S2). These core Runx-

319

activated genes (100 genes) were generally upr

egulated as normal pro-T cells advance

320

to the DN2 and DN3 stages (e.g.,

Cd24a, Hes1, Patz1, Ahr, Myb, Lck, Tcf7, Ly6d, Bcl11b,

321

Lat, Cd3d, Cd3g,

and

Gzma

). In contrast, core genes inhi

bited by Runx factors (46 genes)

322

included ETP signature and non-T genes (e.g.,

Mef2c, Lmo2, Cd34, Pou2f2, Csf1r,

323

Csf2rb, Bcl11a, Id1, Ly6a,

and

Cd81

)(Fig. 3e, f, and S4f).

324

These impressions from landmark genes

were supported globally by Single-

325

sample Gene Set Enrichment Analysis (ssG

SEA)(Fig. 3g). GSEA used curated stage-

326

specific thymocyte gene sets to track devel

opmental progression at different absolute

327

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

times of differentiation, show

n in time-resolved population hist

ograms. Controls started

328

with high ETP (DN1) and low DN2 or DN3 enr

ichment scores and shifted to low ETP

329

(DN1) and high DN2 and DN3 values from day 2

to day 6 post perturbations. Runx dKO

330

cells, at both day 3 and day 6, were more signi

ficantly correlated with ETP signatures and

331

slightly increased myeloid si

gnatures, failing to activate DN2 or DN3 signatures

332

comparably to controls. In contrast, Runx

1 OE at both timepoints showed accelerated

333

loss of associations with ETP and increases

in DN2, DN3 profiles accelerated by about a

334

day relative to controls (Fig. 3g).

335

Thus, single-cell transcriptome analysis suggested that Runx activity preferentially

336

activated genes critical for T-developmental

progression, while inhibiting the genes

337

associated with progenitor and myeloid programs.

338

339

Runx levels controlled T-developmental sp

eed via selective activities on discrete

340

T-identity and lymphoid program modules

341

The fine structure of Runx effects on developmental progression could be

342

measured quantitatively using ps

eudotime. We calculated pseudotime trajectories with

343

Monocle3, defining t

he root cells by high expression of

Flt3

and

Kit

and absence of

Tcf7

344

and

Il2ra

transcripts, based on the phenotype of

the earliest thymus-seeding ETPs

25, 51,

345

. As expected, from day 2 pi to day 6 pi

, cells in control clusters showed gradual

346

pseudotime progression wit

h UMAP2 value (Fig. 4a, b).

Runx dKO displayed slower

347

progression in pseudotime as

compared with control cells, while Runx1 OE markedly

348

accelerated progression (Fig. 4a. b).

349

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint

For control cells, UMAP2 and pseudotime

parameters were correlated, with a

350

highly linear relationship (Pearson correlation score

= 0.92) (Fig. 4c), but Runx dKO and

351

Runx1 OE samples showed weaker linearity t

han the controls (Fig. 4c middle and the

352

right panels). Runx dKO cells exhibited co

nsistently slowed developmental (pseudotime)

353

progression across most UMAP2 values; in

contrast, Runx1 OE cells accelerated

354

development (pseudotime), especially dur

ing a specific low-UMAP2 window (“OE-

355

accelerated-UMAP2-window”: UMAP2: -30 to 5).

This implied a specific early-Phase 1

356

window of opportunity when elevated Ru

nx1 dosage was most

effective.

357

To define the target genes involved, we

compared differentially

expressed genes

358

(DEGs) between control vs. Runx1 OE groups

specifically within t

he stages when they

359

presumably diverge, focusing on cells within

the same “OE-accele

rated-UMAP2-window”

360

(-30 < UMAP2 < 5) (Fig. 4d), approximately corre

sponding to clusters 4, 2, and 0 (control)

361

vs. clusters 3 and 8 (Runx1 OE). Amon

g 411 DEGs, 234 genes were more highly

362

expressed in the OE group than

the control, while 177 were

higher in the control group

363

(Table S3). The genes upregulated by Runx1 OE

in this focused comparison included

364

multiple key T-identit

y genes and driver TFs (

Cd3g, Cd3d, Bcl11b, Tcf7, Lck, Gata3, Myb,

365

Patz1, Hes1

), which were induced not only much earli

er but also at higher levels than in

366

controls. However, Runx activation targets

were not entirely T-lineage specific, as Runx1

367

OE also caused increased expression of genes (

Zbtb16, Nfil3, Clnk, Cd160

) associated

368

with innate lymphoid cell (ILC)

programs, in which Runx3

is more highly expressed

369

These genes are normally transiently activa

ted in DN1-DN2a cells

but repressed during

370

later T cell development by Bcl11b and E-proteins (Fig. 4e, Table S5)

32, 54, 55

371

CC-BY-NC-ND 4.0 International license

available under a

(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint

this version posted November 19, 2022.

;

https://doi.org/10.1101/2022.11.18.517146

doi:

bioRxiv preprint