key: cord-0315629-5f171dyg
authors: Li, Kun; Wu, Yang; Li, Young; Yu, Qiaoni; Tian, Zhigang; Wei, Haiming; Qu, Kun
title: Landscape and dynamics of the transcriptional regulatory network during natural killer cell differentiation
date: 2019-03-10
journal: bioRxiv
DOI: 10.1101/572768
sha: 829ba7eef269c1e5ff6dbf72bfe5a218fc30644e
doc_id: 315629
cord_uid: 5f171dyg

Natural killer (NK) cells are essential in controlling cancer and infection. However, little is known about the dynamics of the transcriptional regulatory machinery during NK cell differentiation. In this study, we applied assay of transposase accessible chromatin with sequencing (ATAC-seq) technique in a self-developed in vitro NK cell differentiation system. Analysis of ATAC-seq data illustrated two distinct transcription factor (TF) clusters that dynamically regulate NK cell differentiation. Moreover, two TFs from the second cluster, FOSL2 and EGR2, were identified as novel essential TFs that control NK cell maturation and function. Knocking down either of these two TFs significantly impacted NK cell transformation. Finally, we constructed a genome-wide transcriptional regulatory network that provides an understanding of the regulatory dynamics during NK cell differentiation.

Introduction 47 significance ( Figure 2B ,10 -8 <P<10 -4 ), keeping cells alive, proliferation and prepare to differentiate. 153

Cluster II comprises 4817 elements, which were highly accessible at the late stage (days 24~35) 154 of NK cell differentiation. GREAT analysis revealed that peaks in cluster II were strongly enriched 155 (P<10 -50 ) for immune-relevant GO terms, such as immune system process, immune response and 156

immune system development among others (Figure 2B) , indicating that the epigenetic states of 157 functional immune cell-specific genes were activated throughout the process and then drove 158 progenitor cells differentiation toward NK cells. Cluster III consisted of only 356 peaks, which 159 showed no enrichment for GO terms and were therefore neglected in downstream analysis. 160

We next examined whether the DNA accessibility chromatin signatures at different stages 161 correlated with those of the corresponding gene expressions. By comparing the ATAC-seq profiles 162 with the genome-wide microarray data during NK cell differentiation, we found that gene loci that 163 gained chromatin accessibility were significantly increased in the gene expression level (P<10 -6 ), 164 while gene loci that lost chromatin accessibility had decreased expression (P <10 -9 ), indicating a 165 high correlation between epigenetic and RNA profiles ( Figure 2C) . 166

Chromatin structure and epigenetic modifications regulate gene expression [34] , however, 167 the chronological order of the dynamics of the chromatin accessibility, gene expression, and cell 168 phenotype have not yet been well studied. Here, we integrated the ATAC-seq signals, microarray 169 profile and percentage of NK cell counts to examine the temporal dynamics of these three features. 170

Interestingly, we noticed that the accessibility of mature NK cell-specific peaks (cluster II in 171 Motif searches of differential ATAC-seq peaks could potentially provide information about 198 the transcriptional regulatory network, but one caveat of this method is unable to distinguish 199 between similar motifs of binding TF family members. Thus, we sought to apply the Genomica's 200 module map algorithm and "TF footprint" analysis to better predict TF occupancy on accessible 201 sites by integrating the ATAC-seq and gene expression microarray profiles (Figure 3) In addition, several novel TFs were also enriched, namely, FOSL2 and EGR2, and are 217 potential regulators of NK cell differentiation. These TFs were assessed along with the other 218 enriched TF families from the motif analysis to identify which gene in the family was expressed 219 or differentially expressed during NK differentiation to further filter out true regulators. At each 220 time point, we plotted both the expression value (colored from red to green to represent high to 221 low expression in the gene microarray profile) and the motif enrichment score (shown by the circle 222 size representing the -log(P-value) of the enrichment) in the same figure (Figure 3C) , and we 223 observed that the known regulators were highly expressed and enriched at different stages, same 224 as the genes FOSL2 and EGR2. By integrating results from the above three analyze ( Figure 3A -225 C), we predict FOSL2 and EGR2 as potential regulators. 226 DNA sequences that are directly occupied by DNA-binding proteins are protected from 227 transposition, and therefore the resulting sequence "footprint" reveals the presence of a DNA-228 binding protein at each site, analogous to DNase digestion footprints. TF "footprint" analysis of 229 our ATAC-seq profile provided further evidence of direct occupancy of a TF candidate on genomic 230 DNA and thus refined the prediction of potential regulators. We illustrated the "footprints" of 2 231 known regulators, STAT5 and T-bet, and observed higher DNA accessibility and deeper 232 "footprints" flanking their motifs in the late compared with the early stage during NK cell 233 differentiation ( Figure 3D ). Similarly, "footprints" of the TFs FOSL2 and EGR2 were also deeper 234 and more accessible at the late stage, suggesting not only that the motifs of these 2 TFs were 235 enriched at stage-specific peaks but that they were most likely physically bound to those accessible 236 chromatin sites, indicating that they are functional regulators of NK differentiation ( Figure 3D) . 237

Overall, the results from footprint analysis were also consistent with that from the HOMER and 238

Genomica's motif enrichment analysis. 239

By combining the TF motif and "footprint" analysis, we can theoretically predict genes 240 that are regulated by any TF of interest. We performed this analysis on the TFs RUNX, T-bet, 241 FOSL2 and EGR2 and integrated the gene expression profiles at each time point (Figure S5B , 242 Table S2 ). We found that genes regulated by each of these TFs were also highly expressed at late 243 stages and were significantly enriched in GO terms of immune system construction and other 244 related functions. Transcription factors ETS1 was reported to drive early stages of NK cell differentiation [53], and we found that there was a FOSL2 binding site in a dynamical accessible 246 sites on the gene body of ETS1, suggesting FOSL2 might regulate NK cell differentiation through 247 ETS1 ( Figure 3E ). Transcription factors GATA3 was found regulating liver-specific NK cells, 248 IFN-γ production, and T-bet expression in mice [54] . Similarly, we found that both FOSL2 and 249 EGR2 bind to the gene body of GATA3 (Figure 3E) , suggesting that the TFs FOSL2 and EGR2 250 may also regulate NK cell differentiation through regulation of GATA3. 

was not affected at the same time ( Figure 4E-F) . These results indicate that knockdown of FOSL2 278 and EGR2 expression, but not viral infection, inhibit NK cell differentiation, suggesting that 279 FOSL2 and EGR2 are necessary for NK cell differentiation. We then tested another important 280 marker of NK cell maturation CD11b [1], and found that it was significantly reduced in GFP-281 positive NK cells, suggesting that FOSL2 and EGR2 might affect NK cell maturation ( Figure 4G) . 282

Overall, we predicted these two TFs FOSL2 and EGR2 as key regulators based on ATAC-seq and 283 gene microarray profile analysis, and then experimentally verified that they were indeed necessary 284 for normal NK cell differentiation. 285

In addition, we attempted to identify the signaling pathways through which FOSL2 and 286 EGR2 were involved in driving NK cell differentiation by performing module map analysis of all 287 the differential peaks. 18 modules were identified according to the patterns of their accessibility 288 (Table S3) We then examined the gene expression profiles of these TFs and found that 120 TFs were 304 expressed in at least one stage during the differentiated process. By applying differential analysis, 305

we obtained 14-26 TFs that were distinctly expressed at each stage (fold change>1.5) (Figure S7) , 306 and defined them as the nodes of the regulatory network. The connections (edges) between any 2 TFs were defined as follows: If the TF A motif is on the promoter of TF B, then we say TF A 308 regulates TF B and draw an arrow from TF A to TF B. Here, only TFs that were expressed at the 309 specific time point were consideration [58] . Using this method, we constructed the transcriptional 310 regulatory network at each time points with both the enrichment (P-value) and expression 311 information for all the relevant TFs (Figure 5A-E) . Interestingly, day 7 specific TFs were densely 312 interconnected at the beginning, and quickly vanished after two weeks ( Figure 5A ). In contrast, 313 the day 35 specific network gradually grew out through the induction of relevant TFs. Many known 314 regulators, such as EOMES, TBX21, ETS1, PRDM1, and GATA3 and also FOLS2 and EGR2 315 were increasingly enriched on the network (Figure 5E ). Similar phenomena were also observed 316 on other networks (Figure 5B-D) . We believe the dynamics of the transcriptional regulatory 317 network explain the increase in the proportion and the differentiation of NK cells. In short, this study provides an epigenomic landscape and dynamics of NK cell differentiation and 328 presents foundational profiles for studying the relationship between chromatin accessibility, gene 329 expression and cell growth during this process (Figure 2) . 330

TFs bind to their motifs and are often obligate nucleosome evictors and the creators of 331 accessible DNA sites, and therefore we can use ATAC-seq to predict critical regulators in NK cell 332 differentiation [26] . By motif analysis in HOMER, we found that several known TFs were enriched 333 at different stages during NK cell differentiation. Similar results were also obtained using 334 Genomica ( Figure 3B) . The discovery of known regulators strongly suggested the reliability of 335 our analysis Furthermore, by integrating results from HOMER, Genomica and motif footprint 336 analysis, we identified two novel TFs, FOSL2 and EGR2 that were essential for NK cell maturation. 337

Knockdown of either of these two TFs significantly inhibited NK cell transformation in the in vitro 338 NK cell differentiation system. Module map analysis suggested these two TFs may regulate NK 339 cell through the JAK-STAT pathway, and therefore further studies of this pathway may facilitate 340 the generation of NK cells and thus promote the NK cell-based immunotherapy. Overall, this study 341 also provided a framework to identify new regulators from chromatin accessible data for NK cell 342

TFs do not usually function alone, they always interact with other molecules to fulfill their 344 unique roles. Hence, we depicted the transcriptional regulatory networks at different stages during 345 NK cell differentiation. In order to construct a stable transcriptional regulatory network, we 346 performed a rigorous screening of TFs to avoid stochastic fluctuations, and integrated both the 347 enrichment (P-value) and expression information for all the relevant TFs. Therefore, even a change 348 of either the enrichment or gene expression cutoff may result in different networks, the most 349 critical TFs to the regulatory process still remain. However, since the screening of TFs mainly rely 350 on the gene microarray, which is not as accurate as RNA-seq, the structure of these predicted 351 network may not as robust as well. 352

From our previous study, we noticed that with a minimal cytokine cocktail, we can generate 353 sufficient number of functional NK cells that express the cytokines necessary for NK cells and and then, the lentivirus particles were concentrated by ultracentrifugation at 50000 g for 2 h at 4°C. 396

Finally, the virus particles were gently resuspended in HBSS and stored at -80°C. After UCB 397 CD34 + cells were cultured with multiple cytokines for 14-18 d, we incubated the lentivirus and the 398 cultured cells with polybrene (5 µg/ml) and centrifuged them at 1000 rpm for 70 min at 10°C.

values less than 0.05 were considered statistically significant. 401

Gene Expression Array. The signal values of the samples were normalized using RMA. 403

Differential analysis was performed using Student's t-test as previously described [19] . The Stage-specific peaks: Each stage (e.g., day 7) consists of genes that were induced (>1.5-fold 417 change) at that stage compared with all the other stages. We defined peaks that were accessible 418 only in one stage as stage-specific peaks, and those that were accessible at all stages as 419 conservative peaks. TFs that regulates NK cell development, the pink box indicates the new TF whose function will 672 be experimentally tested later in Figure 4 . TCF3   TTK   TGIF1   THRA   TGIF2   ZNF711   MYCN   IRF4   MAZ   NFATC1   MYC   MEF2C   TAL1   SOX4  PBX1  RBFOX2   NFE2   SPI1   TGIF2   ERG   TGIF1   GFI1B   MEF2C   THRA   MYC   TTK   MAZ   IRF4  ZNF711   TCF3   TAL1   SOX4   MYCN   PBX1  SPI1   NFATC1   NFE2   TCF4   RBFOX2   MAZ   MEF2C   TGIF1   THRA   IRF4   MYC   ERG   TTK  ZNF711  GFI1B   TGIF2   NFE2   SOX4   TCF4   TAL1   PBX1   NFATC1   SPI1   TCF3   MYCN   RBFOX2   MAZ   TGIF1   ZNF711   MYC   THRA   TTK  GFI1B   IRF4   MEF2C   ERG   TGIF2   PBX1   NFATC1   TCF4   SOX4   TAL1   SPI1  NFE2   TCF3   MYCN   RBFOX2  PBX1  SOX4   TCF4   TAL1   NFATC1   RBFOX2  SPI1   MYCN   TCF3   NFE2   ERG   TGIF2   MAZ   ZNF711  GFI1B   THRA   TGIF1   TTK   MEF2C   MYC   IRF4 

NK 505 cells differentiated from bone marrow, cord blood and peripheral blood stem cells exhibit similar 506 phenotype and functions

Enhancement of 508 human cord blood CD34+ cell-derived NK cell cytotoxicity by dendritic cells

Combined IL-15 and IL-12 drives the generation of CD34(+)-derived natural killer cells with 512 superior maturation and alloreactivity potential following adoptive transfer

Programmed differentiated natural killer 515 cells kill leukemia cells by engaging SLAM family receptors

Differential requirement 517 for the transcription factor PU.1 in the generation of natural killer cells versus B and T cells

The transcription factors T-bet 520 and Eomes control key checkpoints of natural killer cell maturation

Transposition of native 522 chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins 523 and nucleosome position

The epigenetic 525 landscape of T cell exhaustion

Epigenetic 527 stability of exhausted T cells limits durability of reinvigoration by PD-1 blockade

Nfib Promotes 530 Metastasis through a Widespread Increase in Chromatin Accessibility

Individuality and variation of 532 personal regulomes in primary human T cells

Chromatin Accessibility 534 Landscape of Cutaneous T Cell Lymphoma and Dynamic Response to HDAC Inhibitors

Epigenetic control 537 of innate and adaptive immune memory

Progress report on Safe 539 VISITOR: approaching a practical instrument for terahertz security screening. Passive Millimeter-540 Wave Imaging Technology Xiii

Location and cellular stages of natural killer cell development

Human 544 NK cell receptors/markers: a tool to analyze NK cell development, subsets and function

Molecular portraits 547 of human breast tumours

GREAT improves 549 functional interpretation of cis-regulatory regions

Heterochromatin and epigenetic control of gene expression

Runx proteins are involved 553 in regulation of CD122, Ly49 family and IFN-gamma expression during NK cell differentiation

CCAAT/enhancer binding proteins are critical components of the transcriptional regulation of 557 hematopoiesis (Review)

1 expression 559 delineates heterogeneity in primary Th2 cells

PU.1 regulates the commitment 561 of adult hematopoietic progenitors and restricts granulopoiesis

Dynamic regulation of PU.1 expression in 563 multipotent hematopoietic progenitors

The transcription factor 565 PU.1 controls dendritic cell development and Flt3 cytokine receptor expression in a dose-566 dependent manner

Distinctive 568 and indispensable roles of PU.1 in maintenance of hematopoietic stem cells and their 569 differentiation

Mature natural killer cell and lymphoid tissue-571 inducing cell development requires Id2-mediated suppression of E protein activity

A novel Ncr1-Cre 574 mouse reveals the essential role of STAT5 for NK-cell survival and development

RUNX proteins in transcription factor networks that 577 regulate T-cell lineage choice

JASPAR 2016: a major 579 expansion and update of the open-access database of transcription factor binding profiles

The JAK-STAT 582 pathway: impact on human disease and therapeutic intervention

The IRF family transcription factors in 584 immunity and oncogenesis

IRF family of transcription factors as 586 regulators of host defense

Deficiency in 588 the transcription factor interferon regulatory factor (IRF)-2 leads to severely compromised 589 development of natural killer and T helper type 1 cells

IFN Regulatory Factor-2 Deficiency 591 Revealed a Novel Checkpoint Critical for the Generation of Peripheral NK Cells

A novel transcription 594 factor, T-bet, directs Th1 lineage commitment

T-bet 596 regulates the terminal maturation and homeostasis of NK and V alpha 14i NKT cells

Gene 599 deregulation and chronic activation in natural killer cells deficient in the transcription factor ETS1

GATA-3 602 promotes maturation, IFN-gamma production, and liver-specific homing of NK cells

SOX4 is a direct target gene of FRA-605 2 and induces expression of HDAC8 in adult T-cell leukemia/lymphoma

The transcription factors Egr2 608 and Egr3 are essential for the control of inflammation and antigen-induced proliferation of B and 609 T cells

Roles of LAG3 and EGR2 in regulatory T 611 cells

Densely interconnected transcriptional circuits control cell states in human hematopoiesis

Emerging natural killer cell immunotherapies: large-scale ex vivo 616 production of highly potent anticancer effectors

New aspects of natural-killer-cell surveillance 618 and therapy of cancer

Expansion of highly 620 cytotoxic human natural killer cells for cancer cell therapy

ATAC-pipe: general analysis of genome-wide 622 chromatin accessibility

Differential expression analysis for sequence count data

WGCNA: an R package for weighted correlation network analysis. 626

STRING: a web-server to retrieve and display the 628 repeatedly occurring neighbourhood of a gene

Database Resources of the BIG Data Center

The y axis represents the number of TFs that pass the fold change cutoff. 776 C: The ratio of known TFs versus all known TFs that qualify the different fold change cutoffs. 777Known TFs: TFs regulating NK cell development reported in the literature.