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

Abstract

Natural killer (NK) cells represent the cytotoxic member within the innate lymphoid cell (ILC) family that are important against viral infections and cancer. Although the NK cell emergence from hematopoietic stem and progenitor cells through multiple intermediate stages and the underlying regulatory gene network has been extensively studied in mice, this process is not well characterized in humans. Here, using a temporal in vitro model to reconstruct the developmental trajectory of NK lineage, we identified an ILC-restricted oligopotent stage 3a CD34-CD117+CD161+CD45RA+CD56- progenitor population, that exclusively gave rise to CD56-expressing ILCs in vitro. We also further investigated a previously nonappreciated heterogeneity within the CD56+CD94-NKp44+ subset, phenotypically equivalent to stage 3b population containing both group-1 ILC and RORγt+ ILC3 cells, that could be further separated based on their differential expression of DNAM-1 and CD161 receptors. We confirmed that DNAM-1hi S3b and CD161hiCD117hi ILC3 populations distinctively differed in their expression of effector molecules, cytokine secretion, and cytotoxic activity. Furthermore, analysis of lineage output using DNA-barcode tracing across these stages supported a close developmental relationship between S3b-NK and S4-NK (CD56+CD94+) cells, whereas distant to the ILC3 subset. Cross-referencing gene signatures of culture-derived NK cells and other noncytotoxic ILCs with publicly available data sets validated that these in vitro stages highly resemble transcriptional profiles of respective in vivo ILC counterparts. Finally, by integrating RNA velocity and gene network analysis through single-cell regulatory network inference and clustering we unravel a network of coordinated and highly dynamic regulons driving the cytotoxic NK cell program, as a guide map for future studies on NK cell regulation.

Graphical abstract

Figure 1.

Dynamic changes across Lin⁻CD34⁺ progenitor and CD56⁺ cell compartments over time during in vitro differentiation. (A) Schematic diagram summarizes in vitro generation of NK cells from Lin⁻CD34⁺ progenitors. (B) UMAP created from multiparameter FACS analysis of cells collected from each time point during 4-week culture, with annotated CD34⁺ and CD56⁺ clusters. Bar graphs show proportions of CD34⁺ and CD56⁺ cells across the different time points from 3 individual donors. (C) NK cell in vitro trajectory projected on UMAP embedding separated into 4 clusters based on sequential upregulation of key developmental markers: CD34 at stage1/2, CD161 and CD56 at stage 3a and 3b, and CD94 at stage 4. (D) Changes in proportions of stage1/2, stage 3a, stage 3b, and stage 4 clusters across different culture time points. (E) UMAP shows expression levels of key NK cell receptors upregulated within CD56⁺ cluster during 4-week culture.

Figure 2.

Dissecting the overlapping and distinct phenotypic features between immature stage 3 NK cells and ILC3. (A-C) Multiparameter FACS analysis of cells collected at each week during 4-week in vitro culture. (A) UMAP from multiparameter FACS of combined cell surface markers and intracellular TFs expression showing 4 sequential stages in NK cell trajectory as in Figure 1, with stage 3 phenotype further divided based on bifurcated NKp44 expression. (B) UMAP generated as in panel A showing expression levels of ILC-associated TFs: TBET, EOMES, and RORγt with annotated CD34⁺, CD161⁺, and CD56⁺ clusters in gated regions. (C) Comparisons of levels of NK receptors (DNAM-1 and NKG2D), and ILC3-associated markers (CD117 and CD161) between RORγt⁻ (as S3b NK) and RORγt⁺ (as ILC3) within the CD56⁺CD94⁻NKp44⁺ population by intracellular and cell surface FACS staining after 4 weeks of in vitro culture. Numbers indicate percentages of total CD56⁺CD94⁻NKp44⁺ cells in each quadrant, data obtained from 3 donors.

Figure 3.

TFs and effector molecules expression coupled with the progression of sequential NK cell stages along the in vitro trajectory. (A) Representative FACS gating scheme to define S1/2 (Lin⁻CD34⁺CD161⁻CD56⁻CD94⁻), S3a (Lin⁻CD34⁻CD161⁺CD45RA⁺CD56⁻CD94⁻), S3b-NKp44⁻ (Lin⁻CD34⁻CD56⁺CD94⁻NKp44⁻), S3b-NKp44⁺ (Lin⁻CD34⁻CD56⁺CD94⁻ NKp44⁺DNAM-1^hiCD161^lo), ILC3 (Lin⁻CD34⁻CD56⁺CD94⁻NKp44⁺DNAM-1^loCD161^hi), and S4 (Lin⁻CD34⁻CD56⁺CD94⁺) populations generated in vitro from Lin⁻CD34⁺ progenitors collected from cultures at 2 weeks (S1/2 and S3a) and 4 weeks (S3b-NKp44⁺, S3bNKp44⁻, ILC3, S4). (B) Violin plots display intracellular expression of GZMB and PRF1 within the subsets S3a, S3b-NKp44⁻, S3b-NKp44⁺, S4, and ILC3, from a combined pool of 3 donors. Statistical significance was performed with Kruskal-Wallis test between pair-wise comparisons S4 NK vs all other subsets. (C) Representative FACS histograms of TBET, EOMES, and RORγt expression within the subsets S3a, S3b-NKp44⁻, S3b-NKp44⁺, S4, and ILC3. Values described mean fluorescence intensity (MFI) of each indicated population. (D) Secretion of IFN-γ and IL-22 by in vitro–generated NK-S3b, NK-S4, and ILC3 cells after 4-week culture in response to stimulations with either PMA/ionomycin, IL-12/15/18, or IL-2/1β/23 compared with untreated cells (UT). Statistical significance was performed with 1-way analysis of variance (ANOVA) between all subsets (n = 9). (E) Upregulation of CD107a by in vitro generated NK-S3b (left panels), NK-S4 (middle panels) and ILC3 (right panels) cells after 4-week culture in response to stimulations with either PMA/ionomycin or K562 target cell (with effector-to-target [E:T] ratio of 1:5), compared with UT cells as controls (n = 3 individual donors, with 3 replicates for each donors). (F) Upregulation of CD107a by in vitro–generated NK cells after 2 weeks of further expansion on mbIL21-41BBL-K562 feeder cells, stimulated with K562 target cells. Resting peripheral blood (PB) NK cells were used as reference control and CB-derived in vitro expanded NK cells were separated into NK-S3b and NK-S4 stages. FACS plots show representative CD107a expression and bar plots represent the proportion of CD107a positive NK cells from PB NK cells and in vitro generated expanded NK cells (n = 3 donors). (G) IFN-γ secretion in response to K562 target cell stimulation from the same experiment settings as described in panel F. (F-G) Statistical significance was performed with t test;∗P < .05; ∗∗P < .01; ∗∗∗P < .001; ∗∗∗∗P < .0001; and n.s., not significant.

Figure 3.

TFs and effector molecules expression coupled with the progression of sequential NK cell stages along the in vitro trajectory. (A) Representative FACS gating scheme to define S1/2 (Lin⁻CD34⁺CD161⁻CD56⁻CD94⁻), S3a (Lin⁻CD34⁻CD161⁺CD45RA⁺CD56⁻CD94⁻), S3b-NKp44⁻ (Lin⁻CD34⁻CD56⁺CD94⁻NKp44⁻), S3b-NKp44⁺ (Lin⁻CD34⁻CD56⁺CD94⁻ NKp44⁺DNAM-1^hiCD161^lo), ILC3 (Lin⁻CD34⁻CD56⁺CD94⁻NKp44⁺DNAM-1^loCD161^hi), and S4 (Lin⁻CD34⁻CD56⁺CD94⁺) populations generated in vitro from Lin⁻CD34⁺ progenitors collected from cultures at 2 weeks (S1/2 and S3a) and 4 weeks (S3b-NKp44⁺, S3bNKp44⁻, ILC3, S4). (B) Violin plots display intracellular expression of GZMB and PRF1 within the subsets S3a, S3b-NKp44⁻, S3b-NKp44⁺, S4, and ILC3, from a combined pool of 3 donors. Statistical significance was performed with Kruskal-Wallis test between pair-wise comparisons S4 NK vs all other subsets. (C) Representative FACS histograms of TBET, EOMES, and RORγt expression within the subsets S3a, S3b-NKp44⁻, S3b-NKp44⁺, S4, and ILC3. Values described mean fluorescence intensity (MFI) of each indicated population. (D) Secretion of IFN-γ and IL-22 by in vitro–generated NK-S3b, NK-S4, and ILC3 cells after 4-week culture in response to stimulations with either PMA/ionomycin, IL-12/15/18, or IL-2/1β/23 compared with untreated cells (UT). Statistical significance was performed with 1-way analysis of variance (ANOVA) between all subsets (n = 9). (E) Upregulation of CD107a by in vitro generated NK-S3b (left panels), NK-S4 (middle panels) and ILC3 (right panels) cells after 4-week culture in response to stimulations with either PMA/ionomycin or K562 target cell (with effector-to-target [E:T] ratio of 1:5), compared with UT cells as controls (n = 3 individual donors, with 3 replicates for each donors). (F) Upregulation of CD107a by in vitro–generated NK cells after 2 weeks of further expansion on mbIL21-41BBL-K562 feeder cells, stimulated with K562 target cells. Resting peripheral blood (PB) NK cells were used as reference control and CB-derived in vitro expanded NK cells were separated into NK-S3b and NK-S4 stages. FACS plots show representative CD107a expression and bar plots represent the proportion of CD107a positive NK cells from PB NK cells and in vitro generated expanded NK cells (n = 3 donors). (G) IFN-γ secretion in response to K562 target cell stimulation from the same experiment settings as described in panel F. (F-G) Statistical significance was performed with t test;∗P < .05; ∗∗P < .01; ∗∗∗P < .001; ∗∗∗∗P < .0001; and n.s., not significant.

Figure 4.

Increased NK lineage restriction during progression through sequential NK cell stages along the in vitro developmental trajectory. (A) Schematic diagram summarizes experimental outline: progeny from in vitro cultured Lin⁻CD34⁺ progenitors after 2 weeks were sorted into S1/2 (Lin⁻CD34⁺CD161⁻CD56⁻CD94⁻), S3a (Lin⁻CD34⁻CD161⁺CD45RA⁺CD56⁻CD94⁻), S3b (Lin⁻CD34⁻CD161⁺CD56⁺NKp44⁺⁻CD94⁻), and S4 (Lin⁻CD34⁻CD161⁺CD56⁺CD94⁺) populations, replated, and lineage output analyzed after 2 and 3 weeks. (B) Representative FACS profiles showing gating scheme of different lineage output after 3-week subculture: Monocyte/DC (CD19⁻CD14⁺CD11c⁺), B (CD14⁻CD11c⁻CD19⁺), CD56⁺ ILCs (CD14⁻CD11c⁻CD19⁻CD56⁺), which were further separated into RORγt⁻ NK and RORγt⁺ ILC3 cells. (C) Lineage output generated from CD34⁺, S1/2, S3a, S3b, and S4 populations after 2 and 3 weeks in subculture (3 individual donors, 3 replicates for each donor). (D) Summary of lineage-output scores aggregated across all the donors and time points, compared between sorted populations as in panel C. (E) Lineage output within CD56⁺ cells generated after 2 and 3 weeks from CD34⁺, S1/2, S3a, S3b, and S4 populations categorized as either wells contained both RORγt⁻ and RORγt⁺ (NK_ILC3), or only RORγt⁻ cells (NK_only; 3 individual donors, 3 replicates for each donor). (F) Summary of CD56⁺ ILC lineage-output scores aggregated across all the donors and time points, compared between sorted populations as in panel E. The lineage scoring criteria are described in detail in supplemental Methods. For panels C-F, statistical significance was performed with Kruskal-Wallis test between pairwise comparisons among S3a subset vs the other 4 subsets; ∗P < .05; ∗∗P < .01; ∗∗∗P < .001; ∗∗∗∗P < .0001; and n.s., not significant.

Figure 4.

Increased NK lineage restriction during progression through sequential NK cell stages along the in vitro developmental trajectory. (A) Schematic diagram summarizes experimental outline: progeny from in vitro cultured Lin⁻CD34⁺ progenitors after 2 weeks were sorted into S1/2 (Lin⁻CD34⁺CD161⁻CD56⁻CD94⁻), S3a (Lin⁻CD34⁻CD161⁺CD45RA⁺CD56⁻CD94⁻), S3b (Lin⁻CD34⁻CD161⁺CD56⁺NKp44⁺⁻CD94⁻), and S4 (Lin⁻CD34⁻CD161⁺CD56⁺CD94⁺) populations, replated, and lineage output analyzed after 2 and 3 weeks. (B) Representative FACS profiles showing gating scheme of different lineage output after 3-week subculture: Monocyte/DC (CD19⁻CD14⁺CD11c⁺), B (CD14⁻CD11c⁻CD19⁺), CD56⁺ ILCs (CD14⁻CD11c⁻CD19⁻CD56⁺), which were further separated into RORγt⁻ NK and RORγt⁺ ILC3 cells. (C) Lineage output generated from CD34⁺, S1/2, S3a, S3b, and S4 populations after 2 and 3 weeks in subculture (3 individual donors, 3 replicates for each donor). (D) Summary of lineage-output scores aggregated across all the donors and time points, compared between sorted populations as in panel C. (E) Lineage output within CD56⁺ cells generated after 2 and 3 weeks from CD34⁺, S1/2, S3a, S3b, and S4 populations categorized as either wells contained both RORγt⁻ and RORγt⁺ (NK_ILC3), or only RORγt⁻ cells (NK_only; 3 individual donors, 3 replicates for each donor). (F) Summary of CD56⁺ ILC lineage-output scores aggregated across all the donors and time points, compared between sorted populations as in panel E. The lineage scoring criteria are described in detail in supplemental Methods. For panels C-F, statistical significance was performed with Kruskal-Wallis test between pairwise comparisons among S3a subset vs the other 4 subsets; ∗P < .05; ∗∗P < .01; ∗∗∗P < .001; ∗∗∗∗P < .0001; and n.s., not significant.

Figure 5.

Establishing clonal relationship between NK cell developmental stages by DNA barcode tracing. (A) Schematic summary of experimental outline: Lin⁻CD34⁺ cells transduced with barcode libraries were allowed for in vitro differentiation for 4 weeks followed by barcode recovery from generated progeny through sequencing. (B) Bar graph represents total number of unique clones recovered for each sorted populations after culture. (C) Matrix of Pearson correlations of clonal composition between samples. (D) Correlation plot with dendrogram by Manhattan distance of clonal composition between samples. (E) Principal coordinate analysis displayed variability for overall clonal content in each sample. (F) Heat map illustrates the distribution of the proportions of the top 10 clones recovered from each sample. (G) Triangular plot of clonal bias distribution between CD56⁺ subsets. Each clone is placed according to its generated proportions between 3 populations (NK3, NK4, and ILC3), and clone size is determined by its contribution in the total clones obtained within 1 experiment. Results from 1 representative experiment.

Figure 6.

Single-cell transcriptomic analysis of in vitro–generated CD56 subsets uncovered diverging NK cell and ILC3 signatures highly resembling primary ILC profiles. (A) Schematic summary of experimental outline: cells generated in vitro from Lin⁻CD34⁺ progenitors after 4-week culture were column-purified for CD56⁺ cells and subjected to scRNA-seq. (B) UMAP obtained from scRNA-seq data displayed CD56⁺ cells (cluster 0-5 and cluster 8) and other non-ILC clusters. (C) Violin plots of expression levels of NK cell receptor transcripts (top) and ILC3 associated transcripts (bottom) among CD56⁺ clusters. Statistical significance was performed with Kruskal-Wallis test between pairwise comparisons cluster 4 against the remaining clusters (top panels), or cluster 5 against the remaining clusters (bottom panels); ∗P < .05; ∗∗P < .01; ∗∗∗P < .001; ∗∗∗∗P < .0001; and n.s. not significant. (D) Dot plots showing percentage (as size) and expression level (as color scale) of cytotoxic genes (top) and ILC3-specific genes (bottom) among cells within CD56⁺ clusters. (E) Enrichment score for “human primary bone marrow” gene set from MSigDB database (left) and “NK cell–mediated cytotoxicity” from KEGG pathway (right) across single cells within CD56⁺ clusters. (F) Enrichment score for custom gene modules derived from bulk RNA-seq data from Collins et al, of ex vivo CD56^dim and CD56^br circulating NK cells, and tonsil-derived ILC3 across single cells within CD56⁺ clusters. (G) UMAP obtained from scRNA-seq data displayed cell clusters generated in 2-week in vitro cultures. (H-I) Circle plots depicted cell–cell interactions from Mono-DC (H) and pDC (I) cluster toward NK-ILC cluster, in which projected chord diagrams shown top 5 pairs of ligand-receptor for these cell interactions inferred from CellChat analysis. pDC, plasmacytoid dendritic cell.

Figure 6.

Single-cell transcriptomic analysis of in vitro–generated CD56 subsets uncovered diverging NK cell and ILC3 signatures highly resembling primary ILC profiles. (A) Schematic summary of experimental outline: cells generated in vitro from Lin⁻CD34⁺ progenitors after 4-week culture were column-purified for CD56⁺ cells and subjected to scRNA-seq. (B) UMAP obtained from scRNA-seq data displayed CD56⁺ cells (cluster 0-5 and cluster 8) and other non-ILC clusters. (C) Violin plots of expression levels of NK cell receptor transcripts (top) and ILC3 associated transcripts (bottom) among CD56⁺ clusters. Statistical significance was performed with Kruskal-Wallis test between pairwise comparisons cluster 4 against the remaining clusters (top panels), or cluster 5 against the remaining clusters (bottom panels); ∗P < .05; ∗∗P < .01; ∗∗∗P < .001; ∗∗∗∗P < .0001; and n.s. not significant. (D) Dot plots showing percentage (as size) and expression level (as color scale) of cytotoxic genes (top) and ILC3-specific genes (bottom) among cells within CD56⁺ clusters. (E) Enrichment score for “human primary bone marrow” gene set from MSigDB database (left) and “NK cell–mediated cytotoxicity” from KEGG pathway (right) across single cells within CD56⁺ clusters. (F) Enrichment score for custom gene modules derived from bulk RNA-seq data from Collins et al, of ex vivo CD56^dim and CD56^br circulating NK cells, and tonsil-derived ILC3 across single cells within CD56⁺ clusters. (G) UMAP obtained from scRNA-seq data displayed cell clusters generated in 2-week in vitro cultures. (H-I) Circle plots depicted cell–cell interactions from Mono-DC (H) and pDC (I) cluster toward NK-ILC cluster, in which projected chord diagrams shown top 5 pairs of ligand-receptor for these cell interactions inferred from CellChat analysis. pDC, plasmacytoid dendritic cell.

Figure 7.

A coordinated transcriptional network accompanied NK cell cytotoxic differentiation trajectory diverged from non-cytotoxic ILC identities. (A) Projections of RNA velocity vectors on the CD56⁺ cell UMAP described the 3 major cellular trajectories governed by EOMES, HOBIT, and RORC. (B) RNA velocity latent-time displayed as color scale, indicates suggested regions of root cells (top-right) and end points (for NK trajectory; bottom-right). (C) Expression levels of progenitor-associated genes: MYC, TCF7, BACH2, and α4β7 integrin. (D) Heat map shows dynamic expression of early and late-acting TFs with cells aligned according to NK cell cytotoxic trajectory. (E) Heat map showing selected top 20 master TFs associated with highly dynamic regulons identified via SCENIC, aligned following the RNA velocity latent-time of NK cell trajectory as in panel D. Custom annotation columns indicated: (i) number of genes associated with each regulon (bar plot), (ii) average regulon activity score (area under the curve [AUC] score) respective for each regulon detected (violin plot). (F) Regulon activity (AUC score) among the top 3 regulons: EOMES, IKZF3, and CEBPD as key drivers along the cytotoxic NK cell trajectory (red) compared with alternative trajectory (tissue resident–like NK cell; blue). (G) Network of coordinated genes governed by EOMES identified via SCENIC and the corresponding biological pathways highly enriched for EOMES-driven network.

Figure 7.

A coordinated transcriptional network accompanied NK cell cytotoxic differentiation trajectory diverged from non-cytotoxic ILC identities. (A) Projections of RNA velocity vectors on the CD56⁺ cell UMAP described the 3 major cellular trajectories governed by EOMES, HOBIT, and RORC. (B) RNA velocity latent-time displayed as color scale, indicates suggested regions of root cells (top-right) and end points (for NK trajectory; bottom-right). (C) Expression levels of progenitor-associated genes: MYC, TCF7, BACH2, and α4β7 integrin. (D) Heat map shows dynamic expression of early and late-acting TFs with cells aligned according to NK cell cytotoxic trajectory. (E) Heat map showing selected top 20 master TFs associated with highly dynamic regulons identified via SCENIC, aligned following the RNA velocity latent-time of NK cell trajectory as in panel D. Custom annotation columns indicated: (i) number of genes associated with each regulon (bar plot), (ii) average regulon activity score (area under the curve [AUC] score) respective for each regulon detected (violin plot). (F) Regulon activity (AUC score) among the top 3 regulons: EOMES, IKZF3, and CEBPD as key drivers along the cytotoxic NK cell trajectory (red) compared with alternative trajectory (tissue resident–like NK cell; blue). (G) Network of coordinated genes governed by EOMES identified via SCENIC and the corresponding biological pathways highly enriched for EOMES-driven network.