A temporal developmental map separates human NK cells from noncytotoxic ILCs through clonal and single-cell analysis

Natural killer (NK) cells are innate lymphoid cells (ILCs) that mediate host protection against viral infections and malignancy. Human NK cells are known to develop from multipotent bone marrow–residing progenitors that progressively transition through multiple intermediate stages into the secondary lymphoid tissues, including the tonsils, and finally enter the peripheral blood as CD56^bright and CD56^dim NK cell subsets with fully developed effector functions.¹ This dynamic process is collectively proposed as a linear model of NK cell development, and can be divided into distinct stages based on combined expression of cytokine and chemokine receptors, along with a vast array of activating and inhibitory NK cell receptors.^2,3 Briefly, stage 1 and 2 consist of bone marrow–resident multipotent CD34⁺ progenitors that can give rise to both myeloid and lymphoid compartments.^4,5 Acquisition of the receptor IL1R1 and the integrin complex α4β7, along with C-KIT expression within CD34⁺ cells marks an ILC-restricted progenitor population.^6-8 The Lin⁻CD117⁺CD127⁺CD94⁻ stage 3 cells present both in the peripheral blood as well as in tonsils, have a potential to generate different ILC subsets, including NK cells.^9-11 The stage 4 NK cells upregulate CD94 receptor, and can be further separated based on the surface expression of NKp80 into S4a (NKp80⁻) subset and S4b (NKp80⁺) population that represents the NK lineage commitment step.¹² Final stages of NK cell maturation according to this model involve the acquisition of the Fcγ receptor CD16 (stage 5) and CD57 that marks a terminal differentiation status (stage 6). In addition, the presence of several unique NK cell progenitor subsets across diverse tissues and during ontogeny may contribute to the diversity of the NK cell repertoire that mediates different functions at local tissues, depending on the specific immunity context.^13-15 

The emerging role of the other ILC subsets that not only share partially overlapping developmental pathway with NK cells, but also carry multiple NK lineage-specific markers, poses a challenge in ongoing efforts to dissect their respective distinct differentiation routes. In particular, the canonical human NK cell marker CD56 is known to be expressed on human ILC1 and ILC3 subsets in vitro and in vivo, together with other nonexclusive NK cell markers such as IL7R, KIT, CD161, and NKp44 (NCR2).^11,16-18 Consequently, the extent to which developmental pathway of NK lineage concurs or diverge from that of noncytotoxic ILC subset remains largely unresolved, and the heterogeneity of in vitro generated CD56⁺ ILCs requires more detailed dissection.

To address this, we present here a comprehensive trajectory of these diverse CD56⁺ ILCs in vitro that establishes the critical NK lineage differentiation checkpoints and how they are overlapped with developmental stages of noncytotoxic ILC subsets. Leveraged from a previously well-established in vitro coculture system, we further explored the clonal diversity and relationship between the NK cell and CD56⁺ ILC subsets by applying a clonal DNA-barcoding technique. Finally, through in silico transcriptomic analysis using RNA velocity and single-cell regulatory network inference and clustering (SCENIC) tool, we provide an in-depth layer of the transcription factor (TF) network driving the NK cell cytotoxic program. Ultimately, the established framework of NK cell development could provide a reference guide map to facilitate novel NK cell expansion strategies, as well as to better understand dysregulated NK cell development during disease.

NK cell differentiation on OP9 stroma was performed as previously described.¹³ Briefly, cord blood (CB) Lin⁻CD34⁺ cells were plated on preseeded OP9 cell layer and cultured for 4 weeks with a combination of cytokines (interleukin-3 [IL-3], first week; IL-7, IL-15, stem cell factor, and FLT3L), the medium was half-replaced after every consecutive week. For the analysis of NK cell generation kinetics, cell samples were collected at consecutive weeks and kept frozen until all time points were obtained for simultaneous phenotypic analysis by fluorescence-activated cell sorting (FACS). In replating experiments, cells generated in OP9 cocultures after 2 weeks were FACS sorted for the corresponding S1/2, S3a, S3b, S3b-NKp44⁺, S3b-NKp44⁻, and S4 phenotypes (see details in supplemental Methods) and placed on new OP9 stroma layers for an additional 2 to 3 weeks before progeny output was determined by FACS analysis.

All FACS experiments for phenotypic analysis and cell sorting were performed according to standard protocol with 30 minutes incubation at 4°C with appropriate fluorochrome-conjugated antibody mixes in FACS buffer (phosphate-buffered saline, 1% fetal bovine serum, 1 mM EDTA) in the presence of Fc blocking reagent. For the intracellular staining, an additional fixation/permeabilization step was included using BD CytoFix/CytoPerm kit for intracellular antigens; and BD Transcription Factor Buffer Set for TF detection. All analyses were performed on BD Fortessa-X20 and FACS sorting was done on BD Aria-III. FCS files were analyzed with FlowJo software (v10.6).

In vitro–generated CD56⁺ cells after 4 weeks were purified using NK cell isolation kit (Miltenyi), and further cultured on irradiated K562 feeder cells with membrane-bound IL-21 and 41BB-ligand,¹⁹ in the presence of IL-2 and IL-15 for a further 2-week expansion (see details in supplemental Methods). After this, cells were counted and analyzed using cytometry by time of flight (CyTOF) panel on a Helios platform (Standard BioTools). CyTOF data analysis was performed with cytometry data analysis tools and diffcyt R packages.^20,21 

FACS-sorted S1/2, S3a, S3b, S3b-NKp44⁺, S3b-NKp44⁻, and S4 populations were collected at different time points during in vitro culture, and total RNA was isolated with RNeasy-Plus Micro kit (Qiagen). Synthesis of complementary DNA was performed using SuperScript IV RT kit (Invitrogen) and Taqman quantitative polymerase chain reaction (PCR) gene expression assay was carried out on QuantStudio1 instrument (Applied Biosystems). Gene expression was calculated using ΔΔCt method and normalized against glyceraldehyde-3-phosphate dehydrogenase or hypoxanthine phosphoribosyltransferase.

Cells generated in vitro after 4 weeks of culture on OP9 stroma were functionally evaluated based on cytokine secretion (interferon gamma [IFN-γ] and IL-22) and in degranulation assay (CD107a), as previously described.^13,22 Briefly, CD56⁺ cells were directly sorted into U-bottom 96-well plates containing either medium alone (untreated), or with phorbol 12-myristate 13-acetate (PMA)/ionomycin, IL-12/IL-15/IL-18, or IL-2/IL-1β/IL-23 (all cytokines at 50 ng/mL), or with K562 target cells. At 6 hours after stimulation, effector cytokines (IFN-γ and IL-22) were detected by intracellular staining with brefeldin A, and surface CD107a expression was measured by FACS in the presence of GolgiStop (BD Biosciences).

Newborn and adult NOD/SCIDγc^null (NSG) mice irradiated with 1.0 or 3.5 Gy, respectively, injected with 500 to 10 000 CB Lin⁻CD34⁺ cells, as previously described,^13,23 were analyzed for the human cell engraftment at 16 weeks after transplantation. The NSG mice carrying humanized ossicles²⁴ irradiated with 2.5 Gy were injected with 1 × 10⁶ CB Lin⁻CD34⁺ cells per ossicle and generated human hematopoietic progeny was analyzed at 3 and 8 weeks later (see details in supplemental Methods).

The CellTag barcode library was obtained from Addgene, further amplified, and packaged into lentiviral particles, as previously described.²⁵ For barcode labeling, CB-isolated Lin⁻CD34⁺ cells were transduced with low multiplicity of infection pooled library lentiviral vectors (transduction efficiency of <10%). The transduced GFP^hiLin⁻CD34⁺ cells were sorted, and ∼570 cells (equivalent to 10% of barcode library diversity) were plated onto OP9 layers for in vitro differentiation. After 4 weeks, the generated progeny was purified via FACS into B cells, ILC3, S3b NK, or S4 NK cells, and barcodes were recovered by total RNA isolation, followed by PCR amplification of barcode libraries and sequencing with Illumina MiSeq platform. Sequenced reads from each sample were further analyzed with the R package barcodetrackR²⁶ (see details in the supplemental Methods).

Single-cell RNA (scRNA) sequencing (scRNA-seq) of CD56⁺ cells generated in vitro in 4-week cultures was performed on 10X Genomics platform using the Chromium Single Cell 3′ Reagent Kit version 3.1 protocol. Unique molecular index–barcoded library was amplified and sequenced with Illumina NovaSeq6000 and obtained sequence reads were processed using CellRanger pipeline. Sample preprocessing, quality control, and cell clustering were performed following the Scanpy pipeline.²⁷ Cross-referencing scRNA-seq with publicly available gene sets (MSigDB and Kyoto Encyclopedia of Genes and Genomes [KEGG]) and ILC gene signatures (derived from a bulk RNA-seq data set from Collins et al²⁸) were carried out with the Scanpy function “sc.tl.score_genes,” under default parameters. For RNA velocity and trajectory construction, spliced vs unspliced transcript matrix counts were quantified with Alevin-fry tool²⁹ and used as input for sc-Velo³⁰ with “dynamical” mode. Construction of regulatory gene networks from the preprocessed scRNA-seq data was carried out using pySCENIC³¹ workflow. For preserving the consensus of the obtained regulons, the workflow was performed twice and retained only regulons commonly detected between these runs for downstream analysis. Subsequently, pseudotime-directed cell trajectories were built with Slingshot,³² and testing for highly dynamic regulons across the examined trajectories was carried out with TradeSeq package.³³ CellChat package³⁴ was applied to infer cellular interactions between non-ILC/NK cells and differentiating NK cells from the scRNA-seq data set.

All statistical tests were performed either with GraphPad Prism or natively in the corresponding R packages, as indicated in the figure legends. Statistical significance was denoted as ∗P < .05, ∗∗P < .01, and ∗∗∗P < .001, and n.s., not significant.

Hematopoietic stem and progenitor cells were collected and isolated from human umbilical CB after informed donor consent with the approval from the ethical review board at Lund University.

Although previous in vitro studies mainly focused on NK cell and other ILC subsets generated at the end point of culture, the kinetics of this dynamic process remained poorly defined. To address this, we performed serial sampling of the progeny generated at consecutive weeks during OP9 coculture and applied multiparameter FACS analysis of key markers defining early NK cell differentiation and NK cell–receptor diversity (Figure 1A; supplemental Figure 1A). Temporally combined uniform manifold approximation and projection (UMAP) created from this approach clearly depicted the gradual decrease of a CD34⁺ cell cluster with time, along with an accumulation of a CD56⁺ cluster (Figure 1B). Next, we superimposed key stage-defining markers proposed in the NK cell linear developmental model onto this UMAP. This further analysis revealed that the in vitro NK cell trajectory can be broadly divided into 4 stages starting with CD34-expressing progenitors (stage S1/2), followed by the upregulation of CD161 (stage S3a), CD56 (stage S3b), and subsequently CD94 (stage S4) (Figure 1C). The transition between these stages was highly consistent between donors (Figure 1D). The emergence of the NK-receptor repertoire diversity also complied with the sequential in vitro differentiation timeline, among which were DNAM-1 at a very early stage, followed by NKp46 and NKG2 members (NKG2A, CD94, and NKG2D). The mature NK cell phenotype, equivalent to stage 5 NK cells in vivo with CD16 and KIR expression was also observed at the end of the UMAP trajectory, highlighting the phenotypic compatibility between in vitro stages and the proposed in vivo model (Figure 1E).

Because the conventional NK cell marker, CD56, was also reported to be expressed on other noncytotoxic human ILC subsets, including NCR⁺ ILC3s and intraepithelial ILC1s,^17,35 we sought to unravel the heterogeneity within in vitro–generated CD56⁺ by integrating the expression of lineage-specific TFs. Consistent with previous results (Figure 1C), the major 4 stages could be identified based on the cell surface markers alone, except for the NKp44 that further separated the CD56⁺CD94⁻ fraction into 2 smaller clusters (Figure 2A). Although the CD56⁺CD94⁻NKp44⁺ phenotype has been used to define human ILC3s, overlaying the ILC3-specific TF RORγt expression onto this population clearly demonstrated that only a minor fraction of these cells represents bona fide RORγt^hi ILC3 identity (Figure 2A-B). As expected, the EOMES^hi cells within the main CD56⁺ cluster were mostly restricted to the S4 phenotype (Figure 2B). To search for the putative surface marker that labels true ILC3 identity and provides a better separation from the stage 3b NK subset, we explored the expression of NK cell–specific receptors (DNAM-1 and NKG2D) and ILC3-associated markers (Figure 2C; supplemental Figure 1E). Only CD161 and CD117 expression was highly correlated with the RORγt^hi ILC3s, whereas other reported ILC3-specific markers, including CD200R1, CXCR6, IL23R, and IL1R1, did not effectively separate ILC3s from the S3b NK cells. As expected, DNAM-1 and NKG2D, known to be specific for early NK cell differentiation stages, were strongly expressed among the RORγt⁻ S3b NK cell subset but weakly by the RORγt^hi ILC3s generated in vitro.

To validate our new strategy to separate different in vitro NK cell developmental stages from the ILC3 subset, we assessed the presence of effector cytotoxic molecules granzyme-B (GZMB) and perforin-1 (PRF1), as well as the TFs TBET, EOMES, and RORγt (Figure 3A). The expression of GZMB and PRF1 steadily increased following sequential stages of NK cell maturation yet remained at a very low level in the CD117^hiCD161^hiDNAM-1^lo ILC3 cells (Figure 3A-B). Furthermore, the progression of NK cell maturation from the S3a toward the S4 stage was accompanied by an increased expression of EOMES and TBET, whereas high level of RORγt was largely restricted to the ILC3 subset (Figure 3C; supplemental Figure 2A). Interestingly, a small S3a cell fraction displayed a high RORγt level (Figure 3C), together with the high expression of NFIL3 and ID2 transcripts (supplemental Figure 2A), suggesting the possibility that S3a population represents an ICLP-like progenitor capable of giving rise to multiple ILC subsets.

Next, we assessed the function of the S3b NK, S4 NK, and ILC3 subsets, by comparing effector cytokine secretion. The IFN-γ release was most robust in the S4 NK cells, whereas lower in the S3b NK subset and barely detectable in ILC3s (Figure 3D). In contrast, as expected, the ILC3 population displayed the highest IL-22 secretion (Figure 3D). Similarly, we evaluated cytotoxic activity among these subsets in degranulation assays. In line with the previous results, the S4 NK cells displayed the highest degranulation response, followed by the S3b NK cells, whereas the ILC3s remained hyporesponsive (Figure 3E; supplemental Figure 2D).

Because the S3b NK cells contained both NKp44⁺ and NKp44⁻ subsets, we also compared their transcriptional profiles and function. Consistent with the similar levels of GZMB and PRF1 expression (Figure 3B; supplemental Figure 2B), both S3b-NKp44⁺ and S3b-NKp44⁻ cells upregulated CD107a after stimulations, with slightly higher levels in S3b-NKp44⁺ population than in the S3b-NKp44⁻ fraction (supplemental Figure 2C). Similarly, both S3b-NKp44⁺ and 3b-NKp44⁻ cells were effectively triggered to secret IFN-γ. Interestingly, the S3b-NKp44⁺ subset produced more IFN-γ after activation with PMA/ionomycin than the S3b-NKp44⁻ population but not in response to K562-target cells. Nevertheless, the S4-NK cells were the most robust responders both in degranulation and IFN-γ production compared with either of the S3b NK subsets (supplemental Figure 2C). In line with their NK lineage transcriptional profile (supplemental Figure 2A), both S3b-NKp44⁺ and 3b-NKp44⁻ cells generated NK cell output with similar efficiency (supplemental Figure 4E). Interestingly, the generated progeny mostly maintained their respective NKp44 phenotype during the further subcultures (supplemental Figure 4E).

Together, these results corroborate the expected functional properties of the assigned NK cell stages and ILC3 subset in vitro.

We also tested the feasibility and potency of in vitro–generated NK cells after secondary expansion on mbIL21-41BBL-K562 feeder cells (supplemental Figure 3A). Importantly, after 2 weeks, the in vitro–generated CD56⁺ cells expanded 10-fold to 15-fold (supplemental Figure 3B) and maintained their cytokine production as well as cytotoxic activity at comparable level to primary peripheral NK cells (Figure 3F-G). Notably, even after secondary expansion, stage 4 NK cells still retained stronger effector function than the S3b NK stage (Figure 3F-G; supplemental Figure 2E). In addition, comprehensive immunophenotypic CyTOF analysis suggested that expanded NK cells from either CB or peripheral blood phenotypically more resembled each other than the resting peripheral NK cells (supplemental Figure 3C-D). These findings supported that NK cells derived from CB could serve as an equivalent resource to the peripheral blood for large-scale expansion.

To address the multilineage potency of the in vitro NK developmental stages, we performed subculturing experiments using FACS-purified S1/2, S3a, S3b, and S4 populations and tested the lineage composition of their progeny generated after an additional 2- or 3-week culture (Figure 4A; supplemental Figure 4A-B). As expected, the CD34⁺ and S1/2 subsets produced multilineage outputs. In contrast, the S3a population showed the most robust generation of CD56⁺ progeny but virtually lacked B and myeloid lineage potential, suggesting that this stage indeed reflected an ILC-restricted commitment step (Figure 4B-C). Likewise, the S3b and S4 populations showed no non-CD56⁺ outputs (Figure 4B-D). Next, we investigate when the RORγt^hiILC3-generating potential is lost along the NK cell stages, to further pinpoint the definite NK-lineage restriction step. The CD34⁺, S1/2, and S3a populations were effectively generating both RORγt^hi and RORγt^lo CD56⁺ILCs, whereas the S3b and S4 populations gradually lost the ILC3-generation output, supporting that these populations are progressively restricted to the NK lineage (Figure 4E-F). Overall, these data support a hierarchy of NK cell differentiation starting from the multipotent S1/2 stage, which is followed by the loss of non-ILC potential and CD56⁺ restriction at the S3a stage, and thereafter increased acquisition of NK cell identity from the S3b to S4 NK cell stages.

We also investigated whether generation of NK cells from Lin⁻CD34⁺ progenitors in vivo undergoes the equivalent stages as during in vitro culture. Indeed, at 16 weeks after transplantation, the S3a, S3b, and S4 NK cell developmental stages were present in the bone marrow of NSG mice engrafted with human cells (supplemental Figure 5A-D). Furthermore, in NSG mice carrying humanized ossicles,²⁴ all stages of human NK cell differentiation were clearly detectable at 3 and 8 weeks after transplantation (supplemental Figure 5E-H).

Together these data demonstrate that the proposed in vitro human NK cell developmental model is compatible with their in vivo differentiation in immune-deficient mice.

By introducing a novel clonal tracing tool³⁶ into our in vitro system, we aimed to resolve the developmental divergence between NK and ILC3 lineages. For this, barcode-labeled CD34⁺ progenitors were cultured for 4 weeks, followed by barcode recovery from S3b NK cell, S4 NK cell, and ILC3 subsets (Figure 5A-B; supplemental Figure 6A-B), together with B cells used as independent controls. Pearson correlation between barcode compositions across these different stages indicated an expected high relationship between B-cell stages, as well as between NK cell stages (Figure 5C). In contrast, the ILC3 subset showed poor correlation with the B cells and moderately higher correlation with the NK3 and NK4 samples, suggesting that ILC3s shared closer clonal composition with the NK-lineage. Indeed, hierarchical clustering with Manhattan distance depicted 2 branches: B-cell subsets were grouped together, whereas the other branch grouped NK cell stages together and placed ILC3 in further distance (Figure 5D). This was in line with the overall pattern of clonal diversity in principle coordinate analysis among the samples (Figure 5E). In addition, a heat map of lineage contribution from the top 10 clones recovered in each subset clearly reflected the broad segregation between the lineages and among each NK/ILC subsets (Figure 5F). Remarkably, within ILC-restricted clones, we observed a largely balanced contribution between S3b NK and S4 NK subsets among each clone whereas the ILC3-biased clones were much more diverse (Figure 5G). This indicates that, in general, clonal bias was only observed between ILC3 and NK cell subsets but not within the NK cell stages. Similar results were consistently obtained from an additional 2 independent clonal tracing experiments (supplemental Figure 6C-G).

To provide an in-depth transcriptional landscape of CD56⁺ cellular diversity at the single-cell resolution, we performed scRNA-seq of CD56⁺ cells generated after 4 weeks, enriched using column purification (Figure 6A). The majority of cells were CD56^hi (NCAM-1) and could be further separated into 7 clusters (Figure 6B), with some residual non-CD56 clusters that were subsequently excluded from downstream analysis (supplemental Figure 7A-C). Notably, cluster 4 expressed the highest CD94 (KLRD1) level, together with KLRK1 (NKG2D), indicating a mature NK cell profile corresponding to the stage 4 NK (Figure 6C). In line with this, high expression of cytotoxic granule transcripts including GZMs family and PRF were also observed in this cluster. In contrast, cluster 5 distinctively expressed CD117 (KIT), CD161 (KLRB1), and NKp44 (NCR2), together with other ILC3-associated genes (RORC, AHR, IL23R, IL1R1, and TNFSF13B) and generally was devoid of cytotoxicity-related genes (Figure 6C-D). Altogether, this suggested that cluster 5 represented ILC3 identity.

To determine the degree of conservation of gene signatures between in vitro–generated NK cell stages and in vivo counterparts, we leveraged the publicly available databases including MSigBD and KEGG pathways and RNA-seq profiles from multiple human NK/ILC subsets isolated from primary tissues.²⁸ In line with the high expression of cytotoxic genes, cluster 4 displayed strong enrichment scores for both human bone marrow NK cell gene set and NK cell–mediated cytotoxic pathway, supporting the resemblance of cluster 4 gene profile to that of primary cytotoxic NK cells (Figure 6E). Among the ex vivo NK/ILC custom-defined modules (supplemental Figure 7D), cluster 4 showed the highest enrichment score for the CD56^dim module, whereas cluster 0 displayed a partial enrichment for CD56^br module, and cluster 5 was strongly enriched for the primary ILC3 profile (Figure 6F; supplemental Figure 7E).

Induced pluripotent stem cells (iPSC) are another important source of NK cells for clinical applications; therefore, we also compared the publicly available profiles of iPSC-derived NK cells³⁷ to our scRNA-seq data to determine the compatibility between diverse NK cell differentiation platforms. Bulk RNA-seq profile from the Goldenson data set contained both iSPC-derived and CB CD34-derived NK cells, in which the iPSCs-NK cell signature strongly mapped to cluster 4, whereas the CD34 CB–derived NK cells more resembled cluster 0, suggesting differences in NK cell maturation levels (Figure 6D-F; supplemental Figure 8A-B). Overall, we could identify the NK cell signatures from different sources of origin, which supports the compatibility of our model with NK cell development under different experimental settings.

In summary, scRNA-seq data of in vitro–generated CD56⁺ cells supported the identification of 3 major cell clusters with transcriptional profiles that highly recapitulated primary ILC gene profiles, with cluster 4 most strongly associated with cytotoxic NK cells, with cluster 5 corresponding to primary ILC3.

To explore the multilineage complexity and cellular interplay within the in vitro cultures, we also investigated potential interactions between the non-ILC/NK lineages and differentiating NK cell populations using CellChat analysis tool.³⁴ In addition to known communications between NK cells and different dendritic cell (DC) subsets,^38-41 we further corroborated the tightly linked interactions between different DC subsets (both myeloid dendritic cell and plasmacytoid dendritic cell) and NK cells mediated by CD161, NKG2D, and CD94 receptors with their respective CLEC2B and HLA-E ligands (Figure 6G-I; supplemental Figure 8C-D). These findings highlighted the intercellular connections during in vitro differentiation that may provide critical support to shape proper NK cell maturation, therefore warranting further explorations.

Finally, we explored the transcriptional regulation involved in the developmental changes at different stages during in vitro NK cell differentiation by using computational analysis and in silico gene network interference using RNA velocity and SCENIC.^30,31 Interestingly, the generated RNA velocity vectors revealed 3 diverging trajectories initiated from cluster 3 and progressing toward the respective end points, which were distinctively defined over the gradient of TF EOMES, HOBIT, and RORC (Figure 7A; supplemental Figure 9A). Indeed, a progenitor-like signature with enhanced levels of BACH2, TCF7, MYC, and α4β7-integrin in cluster 3 suggests that these cells might represent a highly proliferative precursor population (Figure 7B-C). In concert with this, we observed a highly dynamic shift in all the ILC-associated TFs starting from early acting ILC TFs (NFIL3, GATA3, TCF7, and TBX21), a delayed yet stable ID2 level, and late-acting mature NK cell–specific TFs (IKZF3, ETS1, TOX, and EOMES) along the NK cell trajectory continuum (Figure 7D).

Lastly, to unravel this global transcriptional change specifying the NK cell trajectory, we undertook a system-biology approach described in the SCENIC framework,⁴² which captures coordinately regulated TFs and associates them with respective downstream target genes based on curated knowledge of consensus motifs and ENCODE databases. Through this analysis, we found that the most dynamically regulated TFs aligned with NK cell cytotoxic trajectory, among which were ETS1, IKZF3, EOMES, CEBPD, and MAF (Figure 7E-F; supplemental Figure 9B-C). Finally, pathway enrichment analysis of the captured EOMES regulons pointed toward a series of consensus signatures for NK cells, including NK cell hallmark from Human Cell Atlas, type-1 immunological responses, and NK cell–mediated cytotoxicity pathway (Figure 7G; supplemental Figure 9D-F). Altogether, these data strengthen the close link between the identified EOMES regulon to the physiological function of primary human NK cells.

In this study, we explored a comprehensive developmental pathway along the NK lineage, with specific focus on dissecting the diverging points, as well as phenotypic and transcriptional heterogeneity within in vitro generated CD56⁺ cells, which remains incompletely resolved. Through multifaceted approaches, we demonstrated a coordinated and sequential acquisition of stage-specific markers, including CD161, CD56, CD94, and ultimately CD16. These in vitro differentiation stages were similarly identified after transplantation of human CD34⁺ progenitors into immune-deficient NSG mice.^13,23,24 Importantly, this pattern closely recapitulates the proposed linear model of human NK cell development in vivo.^2,43 Of particular interest, we identified a transient S3a population, with the phenotype that mirrored the previously reported ILCP precursor found in various human tissues.^9-11,44 This intermediate S3a population indeed demonstrated a robust capacity to generate CD56-expressing NK/ILC1 and ILC3 subsets.

In attempt to faithfully separate bona fide NK cells from CD56⁺ noncytotoxic ILCs, we screened multiple putative primary ILC markers and showed, here, that CD161 and CD117 expression was highly linked to RORγt levels, whereas other reported ILC3-associated markers^45-47 including CD200R1, IL1R1, CXCR6, or IL23R did not improve separation. This discrepancy could be explained by the differences in the levels of cell surface receptor expression between cultured cells and their respective counterparts isolated from primary tissues.

Results from clonal bias analysis among the CD56⁺ ILC subsets suggested that in most cases, a few dominant clones contributed to the majority of NK lineage output, whereas ILC3-biased clones yielded only low cellular output, consistent with previous studies applying similar experimental conditions.⁴⁸ These findings could be explained by the presence in the culture of NK cell–promoting factor IL-15, which may confer a proliferative advantage to NK cell-biased clones over ILC3-biased clones.^48-50 

One major concern when using culture system to model human ILC development is to what extent in vitro generated cells mirror their respective in vivo counterparts. We address this question by both functional assays and cross referencing our RNA-seq analysis to the multiple data sets from primary ILC subsets. Combined, these results support that in vitro generated subsets indeed highly resemble their in vivo counterparts. Interestingly, within the cluster of cytotoxic NK cells, we also observed a high enrichment score for a gene signature ascribed to human ieILC1, in line with the similarities in gene expression between human NK cells and ILC1.^17,48,51 Although our data consolidated gene signatures between human NK cell and other group-1 ILC subsets, in accordance with previous works, further studies are clearly warranted to fully resolve the complexity of group-1 ILCs and the molecular signals driving their respective development.

We also compared our scRNA-seq profiles with the gene expression data sets of iPSC-derived NK cells.³⁷ In line with this study,³⁷ the CD34 CB–derived NK cells displayed the immature phenotype compared to the iPSCs-generated NK cell gene signature that closely mapped to highly mature NK cell profile. Together, this supports the compatibility of our model with NK cell development in different settings.

Finally, with the aim to establish a global transcriptional regulatory scheme governing the NK lineage trajectory, we undertook a combinatorial approach of integrating RNA velocity and gene network inference using SCENIC, which unraveled a dynamically coordinated wave of key TFs driving NK cell cytotoxic programming. This global switch in phenotypes and gene regulatory program was reflected by an initiation of progenitor-like root cell population, displaying high level of MYC, BACH2, and TCF7 followed by a continuum of transitional stages characterized by mature NK cell regulators such as EOMES, IKZF3 (AIOLOS), XBP1, and ETS1.^28,52-56 Together, our data suggest a consistent framework of gene regulation that drives NK cell identity and functions in both primary tissues and in vitro cultures. In addition, our in vitro model allowed us to investigate interactions between DCs and developing ILC/NK cells and identify a previously unappreciated role of CD161, NKG2D, and CD94 receptors in this process.

Collectively, this framework broadens the insights in transcriptional regulation of NK cells and provides a tool for future studies on the development of new strategies for in vitro NK cell generation and expansion, as well as understanding disease-induced NK cell defect.

The authors thank Anna Fossum, Zhi Ma, and Isabel Hidalgo at Lund Stem Cell Center FACS core facility for technical support. The authors thank Phuong Cao Thi Ngoc and Stefan Lang for critical bioinformatics support, and Marie Magnusson for support in mouse experiments.

This work was supported by grants from the Swedish Cancer Society (E.S. and K.-J.M.), the Swedish Research Council (E.S. and K.-J.M.), the Alfred Österlunds Foundation (E.S.), and the Fysiografiska Foundation (D.N.V.), The Swedish Children's Cancer Society, Sweden's Innovation Agency, and the Karolinska Institutet (K.-J.M.). The work was further supported by the Research Council of Norway, Center of Excellence: Precision Immunotherapy Alliance, the Norwegian Cancer Society, The South-Eastern Norway Regional Health Authority, EU H2020-MSCA Research and Innovation program, Knut & Alice Wallenberg Foundation, Swedish Foundation for Strategic Research, the National Cancer Institute (K.-J.M.).

Contribution: D.N.V. designed the study, performed experiments with contributions from O.Y., M.K., O.K., E.C., and G.T.-D.; K.H.d.L. performed single-cell RNA-seq and bioinformatics analysis; S.H.R. designed and provided critical resource for RNA-seq; S. Soneji performed bioinformatics analysis; H.L. and S. Scheding contributed to humanized ossicle experiments; D.B. and K.-J.M. provided critical input to study design and discussion; E.S. conceived, designed, and supervised the work; and D.N.V. and E.S. analyzed and interpreted data, prepared the figures, and wrote the manuscript with input from all co-authors.

Conflict-of-interest disclosure: K.-J.M. is a consultant for and received research grants from Fate Therapeutics Inc; is a member of the advisory board at Vycellix; and receives research support from Oncopeptides. All relations have been approved by the University of Oslo, Norway, and Karolinska Institute in Stockholm, Sweden. The remaining authors declare no competing financial interests.

Correspondence: Ewa Sitnicka, Lund Stem Cell Center, Lund University, BMC B12, Klinikgatan 26, 22184 Lund, Sweden; email: ewa.sitnicka_quinn@med.lu.se.