To characterize gene regulation in early ILC and NK cell development, we performed RNA sequencing (RNA-seq) on FACS sorted LY6D^neg CLPs (35, 36), preNKPs and rNKPs (23–26) ( Figure 1A and Table S1 ). In agreement with preNKPs representing a developmental stage between CLPs and rNKPs (23), principal component analysis (PCA) revealed three distinct groups with the first component (PC1 60%) positioning the related progenitor subsets in the expected hierarchical order ( Figure 1B ). Further in line with this, preNKPs expressed genes otherwise only expressed (≥0.3 TPM in all replicas) in CLPs or rNKPs ( Figure 1C ).

Next, we investigated the expression pattern of genes known to be crucial for NK and ILC development (37, 38) ( Figure 1D ). Many of these genes, including Bcl11b, Tox2, Zbtb16, Tcf7, Rorc, Id2, and IL2rb, were found to be upregulated at the preNKP stage ( Figure 1D ). Looking directly at genes with significant expression changes in the CLP to preNKP ( Figure 1E left) and preNKP to rNKP ( Figure 1E right) transitions, revealed an overall pattern where preNKPs displayed down-regulation of genes expressed at the CLP stage ( Figure 1F , cluster II, and in part I) and up-regulate genes expressed at the rNKP stage ( Figure 1F , cluster III, and in part IV). When annotated (using Metascape), we as expected found that cluster I and II were enriched for B cell lineage associated genes (35, 36). In contrast, genes in cluster III and IV were enriched for genes associated with the NK cell lineage ( Figure S1A ). Hence, suggesting that the B-lineage associated gene program observed in CLPs (35, 36) is shut down in preNKPs. The preNKPs instead adopt a general ILC gene program before a more refined NK gene program is established in the NK-lineage committed rNKPs (23).

To characterize changes in cell-signaling pathways occurring at the developmental transitions, we performed KEGG pathway analysis ( Figure 1G ) on the differentially expressed genes ( Figures 1E, F ). This revealed a significant enrichment of genes involved in cytokine-cytokine receptor interactions as well as in the PI3K-Akt, MAPK, and Rap1 signaling pathways ( Figure 1G ). This is in line with prior observations of the critical involvement of cytokines and downstream signaling for early NK cell development and maturation (22, 39–41). Interestingly, the cytokine, PI3K-Akt, MAPK, and Rap1 pathways all coalesced on the FOXO family. This by either modulating FOXO localization and activity or by altering expression of genes that are direct transcriptional targets of the FOXO family (42–45). We found that FOXO1 and FOXO3 were co-expressed in LY6D^neg CLP, preNKP, and rNKP on the mRNA level ( Figure 1H ). To further validate this observation, we confirmed that FOXO1 and FOXO3 were expressed at the protein level at all three progenitor stags ( Figure 1I ) as well as in committed NK cells from spleen and BM ( Figures S1B, C ). This prompted us to further explore the role of the FOXO family in NK cell development (32, 34).

To understand the role of FOXO1 and FOXO3 in NK cell development, we utilized Vav-iCre (46) to conditionally ablate FOXO1 (FOXO1^ΔVav) (47) and FOXO3 (FOXO3^ΔVav) (48) individually or in combination (FOXO1,3^ΔVav) throughout the hematopoietic system. Littermates lacking Vav-iCre (mainly FOXO1^flox/floxFOXO3^flox/flox animals) were utilized as controls (WT). While residual FOXO1 and FOXO3 proteins were detected in the conditional mice ( Figure S2A ), the floxed DNA binding domains of both genes were found to be very efficiently deleted by Vav-iCre when investigated at the mRNA level ( Figure S2B ).

Neither the loss of FOXO1 or FOXO3 alone resulted in a significant alteration to the number of NK cells in spleen and BM ( Figures 2A, B ). In sharp contrast, NK cell numbers were severely reduced in the spleen (5-10-fold reduced) and clearly decreased in the BM of FOXO1,3^ΔVav animals (2-5-fold reduced) ( Figures 2A, B ). Altogether, this implies that FOXO1 and FOXO3 display functional redundancy in NK cell development and together are critical for the generation or maintenance of normal NK cell numbers.

We next assessed NK cell maturation in the three FOXO-deficient mouse strains based on CD11b and CD27 expression (49, 50). In this scheme, the first stage of NK cell maturation is characterized by low expression of both CD27 and CD11b (CD27^lowCD11b^low). CD27 is then increased (CD27^highCD11b^low) before the subsequent upregulation of CD11b (CD27^highCD11b^high) and finally CD27 being down-regulated again (CD27^lowCD11b^high) at the last step of maturation. The expression of CD27 on CD27^high NK cells was significantly reduced in FOXO1,3^ΔVav as compared to WT ( Figures 2C, D ) but regardless, all four maturation stages could be distinguished ( Figure 2C ). Analysis of the maturation subsets in spleen revealed that FOXO1,3^ΔVav NK cells were less mature compared to NK cells from WT and single knockouts, with an accumulation of the CD27^lowCD11b^low population and a significant increase of the CD27^highCD11b^low population ( Figure 2E ). Interestingly, there was an accumulation of the terminally differentiated mature CD27^lowCD11b^high NK cells in FOXO1^ΔVav mice, that was not observed in FOXO3^ΔVav and FOXO1,3^ΔVav mice ( Figure 2E ). This suggests that in the periphery FOXO1 might act as a brake for FOXO3-driven maturation as its absence only has an effect when FOXO3 is present. In contrast to what was found in the spleen, we observed no accumulation of the immature CD27^lowCD11b^low subset and only a trend towards mature CD27^highCD11b^high NK cells being reduced in the BM of FOXO1,3^ΔVav mice ( Figure 2F ).

To investigate whether the observed phenotype had any relation to FOXO expression, we quantified FOXO protein expression in early (CD11b^low) and late (CD11b^high) NK maturation ( Figure 2G ) in WT mice. In line with the reduced generation of later CD11b^high NK cells in FOXO1,3^ΔVav mice, we found that both FOXO1 and FOXO3 generally displayed higher expression in the more immature CD11b^low NK fraction ( Figure 2G ). Together this suggests that the FOXO proteins to a higher extent influence developmental progression of CD11b^low NK cells and that peripheral NK cell development is more dependent on the FOXO proteins than BM NK cell development. The latter in line with the higher FOXO protein expression observed in splenic NK cells ( Figures S1B, C ).

We next investigated the impact of the FOXO genes on the expression of activating and inhibitory receptors that control signaling and functionality in NK cells (21, 51). Loss of FOXO3 caused no significant changes in the NK receptor repertoire except for a reduction in NK1.1 expression ( Figures 2H, I and Figure S3A, B ). In contrast, the loss of FOXO1 alone was enough to cause significant changes in DNAM-1 and KLRG1 expression ( Figure 2H ). As FOXO1,3^ΔVav mice displayed a similar phenotype as FOXO1^ΔVav mice ( Figure 2H ), this suggests that the perturbation of DNAM-1 and KLRG1 is caused solely by the loss of FOXO1. In addition, FOXO1,3^ΔVav mice uniquely displayed significant changes in the expression of the activating receptors Ly49D, NKp46 and 2B4 ( Figure 2H and Figure S5A ) as well as the inhibitory receptors Ly49A and Ly49G2 ( Figure 2I and Figure S5B ). Hence, the FOXO transcription factors individually or cooperatively influence the expression of specific activating and inhibitory receptors.

As IL-15 is an important cytokine for NK cell development and survival (52), we investigated if the reduced expression of CD122 - a central component of the IL-15 receptor on NK cells (38) - could be a contributing factor to the observed reduction in NK cell numbers in FOXO1,3^ΔVav mice. We found a significant reduction in CD122 expression on both BM and splenic NK cells from FOXO1,3^ΔVav mice ( Figure 3A ) and that there was a clear correlation between NK cell numbers and CD122 expression in peripheral NK cells ( Figure 3B ). Interestingly, dependence on the FOXO proteins differed between BM and splenic NK cells with the CD122 expression on BM NK cells seemingly only relying on FOXO3 ( Figure 3A ). In contrast, the combined deletion was required to affect CD122 expression in splenic NK cells ( Figure 3A ). This, potentially due to the higher expression of FOXO proteins in splenic NK cells as compared to BM NK cells ( Figures S1B, C ).

It has previously been shown that the expression of the activating receptor NKG2D on NK cells is dependent on IL-15 signaling (53–55) and hence it can be utilized as a surrogate marker for IL-15 signaling. Much resembling the expression pattern of CD122 ( Figure 3A ), we found that NKG2D expression was significantly reduced on both BM and splenic NK ells from FOXO1,3^ΔVav mice but that the dependence on the FOXO proteins varied between BM and spleen ( Figure 3C ). Hence, we concluded that the loss of FOXO results in reduced surface expression of CD122 and decreased IL-15 signaling. This strongly suggests that the reduced NK cell numbers in part can be directly contributed to a diminished IL-15 response in FOXO1,3^ΔVav NK cells.

The significant decrease in NK cell numbers ( Figure 2B ) coupled with the mild phenotype in BM NK cell maturation ( Figure 2F ) in FOXO1,3^ΔVav mice, hinted at a defect in NK cell progenitors. To investigate this, we performed phenotypic analysis of the CLP, preNKP, and rNKP compartments ( Figure 4A ). This showed a significant decrease in the numbers of CLPs in FOXO1,3^ΔVav mice but no overt changes to the number of preNKPs ( Figures 4A, B ). In contrast, the downstream CD122⁺ (IL15Rβ) rNKP population was significantly decreased in the FOXO1,3^ΔVav mice ( Figures 4A, B ). No significant changes were observed in preNKP and rNKP numbers in mice lacking only FOXO1 or FOXO3 alone ( Figure S4 ). Hence, the combined loss of FOXO1 and FOXO3 results in seemingly increased generation of preNKPs from CLPs and developmental block at the preNKP to rNKP transition. The latter argues for a significant part of the observed reduction in NK cell numbers – particularly in BM – being due to decreased generation of early NK cell progenitors.

To verify that the observed phenotype was intrinsic to hematopoiesis, we utilized the CD45.1/2 system and adoptively transfered FOXO1,3^ΔVav BM (CD45.2) to irradiated WT (CD45.1) hosts. The number of NK cells generated 12 weeks post transplantation from adoptively transferred FOXO1,3^ΔVav BM cells were markedly reduced in spleen ( Figures S5A, B left), BM ( Figure S5B middle), and blood ( Figure S5B right). The perturbation in NK maturation was also recapitulated with reconstituted FOXO1,3^ΔVav NK cells displaying accumulation of immature (CD27^lowCD11b^low) NK cells and reduced number of more mature (CD27^highCD11b^high) NK cells ( Figure S5C ). Further, we observed a significant reduction in CLP numbers and block at the preNKP to rNKP transition in the progeny of FOXO1,3^ΔVav BM cells ( Figures S5D, E ). Hence, the presence of normal cells in the transplantation setting does not rescue NK cell development from FOXO1,3^ΔVav donor BM cells, supporting the notion of a cell autonomous FOXO requirement in the regulation of early NK progenitors and NK cell maturation.

With the loss of FOXO perturbing the NK developmental pathway already at the level of the CLP ( Figures 5A, B ), we next sought to investigate if we could identify NK lineage related changes at this step of development. To this end, we performed RNA-seq on LY6D^neg CLP (GR1/MAC1/NK1.1/TER119/CD3ε^-CD11C^-LY6C^-IL-7Rα⁺FLT3⁺KIT^lowLY6D^-) from WT, FOXO1^ΔVav mice, FOXO3^ΔVav mice, and FOXO1,3^ΔVav mice. We found 469 genes that displayed significant differential expression between LY6D^neg CLPs from FOXO1,3^ΔVav and WT mice ( Figures 5A, B and Table S2 ). The observed perturbation in gene expression became progressively more distinct with the loss of both FOXO1 and FOXO3 function ( Figure 5A ). This suggests that the FOXO proteins have mainly synergistic functions at this step of development.

In agreement with earlier reports, we found that the Il7ra gene was significantly down-regulated ( Figures 5A–C ) confirming the known positive regulatory role of FOXO in controlling IL-7Ra expression (9, 16). Further, looking specifically at genes previously described to be NK cell signature genes (28), we found that a significant number of these genes displayed altered expression at the LY6D^neg CLP stage ( Figures 5B, C ). Of note, we found that Tcf7, Id2, Il18r1, Il12rβ1, Cxcl10, and Cxcl9 – all genes encoding proteins important for NK cell development, migration, and functions (31, 56–59) - were up-regulated in the FOXO1,3^ΔVav cells. Conversely, looking at the down-regulated genes, we interestingly found that Ets1 – a gene known to be important for the development of NK progenitors and NK maturation (30, 60) - was significantly down-regulated in CLPs ( Figures 5A–C ). The down-regulation of ETS1 was also confirmed at the protein level in CLPs ( Figure S6A ).

We next investigated whether the changes in the expression of known NK-cell-development genes could be related also to changes in the gene regulatory landscape. Using the assay for transposase-accessible chromatin followed by sequencing (ATAC-seq), we identified close to 46,000 open chromatin regions across analyzed LY6D^neg CLP with the vast majority existing both in WT and FOXO1,3^ΔVav cells ( Figure 5D ). Overall, most of the identified open chromatin regions were localized in intergenic and intronic region ( Figure 5E , left). Next, we identified regions with significant differences (adjusted p-value <0.01 and ≥2-fold change in signal amongst peaks identified in ≥2 samples and having >30 reads) in chromatin accessibility between WT and FOXO1,3^ΔVav. This revealed 297 regions that were mainly localized in intergenic and intronic regions ( Figure 5E , right). Hence, this suggests that loss of FOXO in LY6D^neg CLPs results in relatively few but distinct changes to overall chromatin accessibility and that the major effect on gene regulation is via distal elements while promoters remain largely unaffected.

To identify transcription factors where altered binding could cause the alterations in chromatin accessibility, we performed de novo motif enrichment analysis on the peaks with altered chromatin accessibility. This revealed that peaks with decreased accessibility in the FOXO1,3^ΔVav LY6D^neg CLP were enriched for transcription factor binding (TFBS) sites related to the bHLH-, ETS- and, as expected, the FOXO family ( Figure 5F top). Correspondingly, RUNT, IRF, and ETS motifs were found in regions that gained chromatin accessibility ( Figure 5F bottom). To further corroborate that TF binding was altered, we analyzed the Tn5 integration sites (cut-profiles) around the putative TFBS for each of these TF-families. Out of the motifs identified ( Figure 5F ), ETS and IRF produced clear cut-profiles ( Figures 5G and S6B ). This supports that altered ETS and IRF binding directly contribute to the changes in chromatin accessibility while suggesting that the other TFBS are either not used or that TF binding does not produce distinct cut-profiles on these analyzed regions.

Altered binding of a single transcription factor might not cause significant changes to the overall chromatin accessibility at the level of a whole chromatin region as defined by the ATAC-seq peaks. A complementary approach is footprinting analysis, which instead attempts to localize sudden decreases in the number of reads within an open chromatin region to identify individual TF bound regions (61). By means of footprint analysis, we assessed the changes in genome-wide binding of ETS and IRF. We found no major changes in the overall number of IRF associated footprints ( Figure 5H ). This suggests that increased IRF binding is specifically associated with peaks displaying increased chromatin accessibility in FOXO1,3^ΔVav CLPs ( Figures 5F, G ), while the observed decrease Irf4 expression ( Figures 5C, H ) has no major impact on overall IRF-binding.

In contrast, we found a significant decrease in footprints associated with ETS motifs also at the genome-wide level ( Figure 5H ). The fact that Ets1 is the only identified ETS-family member displaying significant changes in expression ( Figures 5C, H ), suggests that the decreased number of ETS-bound regions reflects reduced binding of ETS1. This conclusion is further supported by putative ETS1 target genes (62, 63) including Stat4, Pdk2, Adck3 showing significantly lower expression in FOXO1,3^ΔVav CLPs ( Figure 5I ). Taken together with the reduced expression of ETS1 ( Figures 5B, C and S6A ) and clear ETS cut-profile ( Figure 5G ), this suggests that loss of ETS1 binding in FOXO1,3^ΔVav CLPs contribute to the altered gene regulatory landscape and potentially the reduced capacity to generate NK cells downstream of the CLP (30).

We further looked specifically at TCF7 (31) and NFIL3 (29), as both are known to be critical for BM NK progenitor development. TCF7 expression was increased in CLPs from FOXO1,3^ΔVav mice ( Figures 5B, C ) but no significant change in binding as assessed by footprinting analysis (46 less footprints in the FOXO1,3^ΔVav CLPs) was observed. Hence, TCF7 is potentially controlled by FOXO at the transcriptional level but only has a minor impact in terms of chromatin accessibility at the CLP stage. NFIL3 showed no significant change in expression or overall binding (seven less footprints in the FOXO1,3^ΔVav CLPs). Hence, we find no indication that perturbed TCF7 and NFIL3 activity contribute to the changes observed in FOXO1,3^ΔVav CLPs.

With FOXO1,3^ΔVav mice displaying a block at the preNKP to rNKP transition, we next sought to characterize the transcriptional changes caused by the loss of FOXO in preNKPs. Based on PCA, the FOXO1,3^ΔVav preNKPs overall maintained a preNKP transcriptional profile as compared to its WT counterparts (PC2) with the top loadings of PC2 being CLP and NK-lineage related ( Figures 6A and S7A ). In line with this, the expression of Id2 – which marks the formation of preNKP from CLP (24) – was not altered and generally the expression of ILC related transcription factors (37) was also found to be similar ( Figure 6B ). However, the combined loss of FOXO1 and FOXO3 did cause distinct gene expression changes as observed both by PCA (PC1) ( Figure 6A ) and direct comparison of expression profiles ( Figure 6C and Table S1 ). In relation to normal development, the FOXO1,3^ΔVav preNKPs displayed lower expression of genes commonly expressed throughout the early NK progenitor hierarchy ( Figure 6D , cluster I) and, to a lesser extent, maintained expression of CLP associated genes ( Figure 6D , cluster IV). We additionally found that the FOXO1,3^ΔVav preNKPs failed to properly express genes normally upregulated in the CLP to preNKP transition and further increased in expression in the subsequent preNKP to rNKP transition ( Figure 6D , cluster II). The later gene cluster was as expected enriched for NK lineage associated genes ( Figure S7B ) including amongst others Klrd1 (encoding for CD94) and Tigit shown to function as NK immune checkpoint inhibitors (65). The reduced generation of rNKPs from FOXO1,3^ΔVav preNKPs is hence associated with a failure to properly express an early NK gene program.

Looking specifically at individual genes impacting NK development, we found that TCF7 - in contrast to what was observed in CLPs ( Figures 5B, C ) - was downregulated in FOXO1,3^ΔVav preNKPs while NFIL3 expression remained unaffected ( Figures 6B, C ). Hence, the lower TCF7 expression in preNKPs could potentially contribute to the impaired generation of rNKPs (31). Further, in line with data from the CLP stage, we found that ETS1 displayed a clear trend toward being down-regulated in FOXO1,3^ΔVav preNKPs ( Figures 6B and S7C ). In concordance with the CLP expression and epigenetic data, this argues for the loss of FOXO causing ETS1 downregulation which in turn contributes to the observed impairment generation of NK cells. Further in line with perturbed ETS1 expression contributing to the observed NK phenotype, we found that Ets1 expression was significantly down-regulated also in splenic NK cells ( Figure S7D ).

Interestingly, we also found a significant decrease in expression of zinc finger protein 105 (Zfp105) ( Figure 5B ) – a transcription factor previously shown to be regulated by FOXO1 (15, 66), to augment differentiation toward NK cell lineage (67) and to be a direct target of FOXO1 in CD8 T cells – which might suggest that FOXO-mediated regulation of Zfp105 expression plays a role in NK cell development.

The preNKP represents a heterogenous population which gives rise not only to NK cells but also to the other non-cytotoxic ILC subsets (24, 25). In agreement with this, we observed expression of genes encoding transcription factors linked to development of the non-cytotoxic ILCs, including Tcf7, Tox, Bcl11b, Zbtb16, Rorc, Ets1, Nfil3 (37) in both WT and FOXO1,3^ΔVav preNKP ( Figure 6B ).

We next sought to investigate if we could observe transcriptional changes in the FOXO1,3^ΔVav preNKP indicative of disruptions potentially influencing also the non-cytotoxic ILC subsets. To this end, we utilized the ILC1-3 gene signature published by Robinette et al. (64). Interestingly, we found down regulation of expression across all the three gene signatures ( Figures 6E–G ) but most notably within the ILC1 ( Figure 6E ) and ILC3 ( Figure 6G ) gene signatures.

These transcriptional changes made us speculate that the development of the non-cytotoxic ILC subsets in addition to NK cells could be perturbed in FOXO1,3^ΔVav mice. To investigate this, we performed phenotypic analysis of ILC subsets from BM, spleen, and thymus from WT and FOXO1,3^ΔVav mice ( Figure 6H ). Indeed, we found that the ILC1 population was reduced in BM and spleen ( Figures 6H, I ). ILC3 numbers were similarly reduced in BM and thymus while, in contrast, NKp46^-ILC3s were increased in spleen ( Figures 6H, I ). In contrast, ILC2 subsets were generally increased in the analyzed organs ( Figures 6H, I ). Hence, we concluded that the loss of FOXO causes broad perturbation of the non-cytotoxic ILC subsets. Potentially these changes are due to gene regulatory changes already at the level of the preNKP, meaning that FOXO determine lineage specification of ILC.

To generate animals conditionally lacking FOXO1 and/or FOXO3 throughout the hematopoietic system, we crossed FOXO1^flox/flox (47) and/or FOXO3^flox/flox (48) with Vav^-iCre (46) mice. All alleles were maintained on a C57BL/6 background and mice were predominantly analyzed at 8 to 14 weeks of age. Congenic CD45.1 WT C57BL/6 mice were used as recipients in transplantation experiments. All animal experiments were approved by the local animal ethics committee.

Single-cell suspension of bone marrow, spleen, blood, or thymus were incubated with Fc block (anti-FcγRIII, clone 2.4G2) and subsequently stained with fluorescent antibodies ( Table S3 ) and viability markers (LIVE/DEAD^® Fixable Aqua Dead Cell Stain Kit or propidium iodide, both from Invitrogen). Staining was done at 4°C in PBS with 2% FBS for 20 min. Results were acquired mainly using the BD LSRFortessa™ or BD FACSymphony™ Flow Cytometers (BD Biosciences).

For FACS sorting of progenitor cells, mature cells were depleted using antibodies against TER119, CD3, CD19, GR1, and MAC1 together with sheep anti-rat IgG Dynabeads (Invitrogen) prior to staining with fluorescent antibodies. Cell sorting was performed mainly on a FACSAriaIIu or FACSAria Fusion (BD Biosciences).

Further analysis of FACS data was performed using Flowjo v9.9.6 (TreeStar, Ashland, OR). Corrected mean fluorescence intensities (MFI) were calculated by subtracting the MFI of control sample (stained only with secondary antibody) from the MFI of the FOXO stained sample. Normalized MFIs were calculated by dividing MFIs with the average MFI observed for WT samples in each independent experiment.

To determine in vivo NK lineage output from FOXO1,3^ΔVav hematopoietic stem- and progenitor cells, unfractionated bone marrow (0.2 x 10⁶ WT or 3 x 10⁶ FOXO1,3^ΔVav CD45.2 donor cells), was injected intravenously into irradiated (950cGy) CD45.1 recipients. The number of unfractionated BM cells transplanted were proportional to the frequency of phenotypic hematopoietic stem cells in WT and FOXO1,3^ΔVav mice (data not shown). Part of the FOXO1,3^ΔVav transplanted animals were in addition given 0.2 x 10⁶ unfractionated WT BM cells as support. Reconstitution was analyzed at 12 weeks post-transplantation using flow cytometry.

RNA-seq was done using TotalScript (Epicenter) on RNA prepared from approximately 5000 FACS sorted cells using RNeasy micro (Qiagen) as previously described (78). Libraries were sequenced pair end (2x50 cycles) on the Illumina platform. Reads were mapped to the mouse reference genome (mm10) using STAR v2.3.2b (https://github.com/alexdobin/STAR). Strand-specific reads in exons were quantified using HOMER and assessment of differential gene expression analysis was done using EdgeR (79) on raw read count. Data visualizations were mainly done using ggplot2, pheatmaps, and R base graphics.

ATAC sequencing was performed (using 3000-5000 FACS sorted cells) as previously described (80). Libraries were sequenced pair-end (2x50 cycles) on the Illumina platform (Illumina). Reads were trimmed (using Trim Galore v0.4.1), mapped to the mouse reference genome (mm10) (using Bowtie2 v2.3.3.1) and PCR duplicates removed when making HOMER tag directories (using makeTagDirectory with -tbp 1). Peaks were subsequently identified in sub-nucleosomal reads (read-pairs within 100bp) using HOMER’s findPeaks.pl. Peaks with differential chromatin accessibility were identified using EdgeR on raw read counts in identified peaks. Peaks displaying an adjusted p-value ≤0.01 and ≥2-fold change in read count were considered to have differential chromatin accessibility. Only peaks identified in ≥2 replicas each with >30 reads were considered in the analysis. Annotation and motif enrichment analysis of differential peaks were done using the HOMER’s annotatePeaks.pl and findMotifsGenome.pl with -size given respectively.

To make cut-profiles, the localization of known HOMER transcription factor binding sites (TFBS) belonging to members of the enriched TF family (identified by the motif enrichment in differential ATAC-seq peaks) were localized in the genome using HOMER’s findMotifsGenome.pl. To take into account the position of the Tn5 integration into the genome, custom HOMER tag directories were made off-setting reads on the plus and minus strand with +4 and -5 bases respectively. Read depth centered around TFBS from a specific family were subsequently plotted using HOMER’s annotatePeaks.pl with a -fragLength of 9 (corresponding to the bp covered by Tn5) and -hist 1 (1 bp bins).

Genome-wide footprinting to identify TF binding sites was performed using DNase2TF. In brief, Bowtie2 mapped reads were deduplicated (using Picacard tools’ MarkDuplicates), data from the same population/genotype merged (using SamTools’ merge) and down-sampled to 39 million read-pairs per sample (using Picard Tool’s DownsampleSam.jar). Localization of reads were off set in the.bam file to take into account the Tn5 integration (as described above) using custom scripts. HOMER tag directories and peak finding were done as described above. Peak files and downsampled.bam files were subsequently used as input for DNase2TF (81). Identified footprints with a p-value ≤0.05 were overlapped with the TFBS identified in the mouse reference genome (mm10) using the transfac catalogue as previously described (82). Footprints were associated with TFBS when the center of the TFBS fell within the footprint.

Data visualizations were mainly done using ggplot2, pheatmaps and R base graphics.

NK cells were FACS sorted using BD AriaIII (BD Biosciences). CD117+ BM cells were enriched using CD117 MicroBeads, mouse (Miltenyi). RNA was purified using RNeasy micro kit (Qiagen) and cDNA prepared using MultiScribe Reverse Transcriptase (Life Technologies) or SuperScript II (Life Technologies) in combination with random hexamer priming. qPCR was performed using TaqMan™ Universal PCR Master Mix (Life Technologies) and TaqMan probes against: FOXO1 (Mm00490672_m1), FOXO3 (Mm00490673_m1), hprt (Mm01545399_m1 or Mm00446968_m1) and Ets1 (Mm01175819_m1).

Two hundred to 500 progenitor cells were FACS sorted into lysis solution with DNase I from the Single Cell Lysis Kit (Invitrogen) and samples prepared according to manufacturer’s instructions. RNA-seq libraries were subsequently prepared using the SMART-Seq Stranded Kit (Takara) according to the manufacturer’s instruction. Quality of cDNA library was determined using an Agilent Bioanalyzer according to the manufacturer’s protocol. Libraries were quantified by using the KAPA-SYBR FAST qPCR kit (Roche) and sequenced pair-end (2x75 cycles) on the Illumina NextSeq 500.

Adaptor sequences were trimmed, and low-quality reads removed using Trimmomatic (v.0.36). All sequencing reads aligning (HiSAT2, v.2.1.0) to annotated mouse ribosomal RNA genes were discarded. High-quality and ribosomal RNA depleted sequencing reads were aligned to the genome GRCm38.p6/mm10 genome using HiSAT2. Using sorted bam files (Samtools v.1.10), the number of aligned reads were counted (featurecount in subread package v. 2.0.0). After normalization (TMM: trimmed mean of M-values), a differential gene expression analysis (edgeR v. 3.28.1) was performed. Significant differentially expressed genes were distinguished by a false discovery rate (FDR) <0.05. Data were plotted using ggplot2 (v.2.3.3) in R (v. 3.6.1). Gene ontology analysis was conducted using clusterprofiler (v.3.14.3) with the database org.Mm.eg.db (v.3.10.0) in R. All scripts used for processing of SMART-seq data are deposited on Github: https://github.com/jonasns/NK_FOXO.

Statistical analysis of FACS data was performed using Graphpad Prism version 6 for Mac OSX (Graphpad 83) or R version 3.3.3 (R Development Core 84). Statistics pertaining to RNA-seq and ATAC-seq data were performed as described above.