# Single-cell profiling reveals the trajectories of natural killer cell differentiation in bone marrow and a stress signature induced by acute myeloid leukemia

Adeline Crinier, Pierre-Yves Dumas, Bertrand Escaliere, Christelle

Piperoglou, Laurine Gil, Arnaud Villacreces, Frederic Vely, Zoran Ivanovic, Pierre Milpied, Emilie Narni-Mancinelli, et al.

► To cite this version:

Adeline Crinier, Pierre-Yves Dumas, Bertrand Escaliere, Christelle Piperoglou, Laurine Gil, et al. Single-cell profiling reveals the trajectories of natural killer cell differentiation in bone marrow and a stress signature induced by acute myeloid leukemia. Cellular and molecular immunology, 2020, 10.1038/s41423-020-00574-8. hal-03099198

HAL Id: hal-03099198

https://amu.hal.science/hal-03099198

Submitted on 6 Jan 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

**ARTICLE** **OPEN**

# Single-cell profiling reveals the trajectories of natural killer cell differentiation in bone marrow and a stress signature induced by acute myeloid leukemia

Adeline Crinier<fnref n="1" />, Pierre-Yves Dumas²,³,⁴, Bertrand Escalière<fnref n="1" />, Christelle Piperoglou⁵, Laurine Gil<fnref n="1" />, Arnaud Villacreces³,⁴, Frédéric Vély<fnref n="1" />,⁵, Zoran Ivanovic⁴,⁶, Pierre Milpied<fnref n="1" />, Émilie Narni-Mancinelli<fnref n="1" /> and Éric Vivier<fnref n="1" />,⁵,⁷

Natural killer (NK) cells are innate cytotoxic lymphoid cells (ILCs) involved in the killing of infected and tumor cells. Among human and mouse NK cells from the spleen and blood, we previously identified by single-cell RNA sequencing (scRNAseq) two similar major subsets, NK1 and NK2. Using the same technology, we report here the identification, by single-cell RNA sequencing (scRNAseq), of three NK cell subpopulations in human bone marrow. Pseudotime analysis identified a subset of resident CD56 bright NK cells, NK0 cells, as the precursor of both circulating CD56 dim NK1-like NK cells and CD56 bright NK2-like NK cells in human bone marrow and spleen under physiological conditions. Transcriptomic profiles of bone marrow NK cells from patients with acute myeloid leukemia (AML) exhibited stress-induced repression of NK cell effector functions, highlighting the profound impact of this disease on NK cell heterogeneity. Bone marrow NK cells from AML patients exhibited reduced levels of CD160, but the CD160 high group had a significantly higher survival rate.

**Keywords:** NK cells; scRNASeq; AML

Cellular & Molecular Immunology _###########_; https://doi.org/10.1038/s41423-020-00574-8

## INTRODUCTION

Natural killer (NK) cells are large granular lymphocytes in the ILC family. NK cells are endowed with the capacity to kill stressed cells, such as infected cells and tumor cells, without specific prior activation, in humans and mice.<fnref n="1" /> NK cells express an array of inhibitory and activating receptors, the engagement of which regulates NK cell activation. Inhibitory receptors are essential for sensing decreases in or total absence of constitutively expressed self MHC-I molecules on target cells. Decreases in MHC-I expression reduce the strength of the inhibitory signals delivered to NK cells, rendering these cells more prone to activation.²,³ NK cell activation results from the engagement of activating receptors, such as the activating isoforms of Ly49 and KIRs, natural cytotoxicity receptors (NCRs), SLAM (signaling lymphocyte activating molecule)-related receptors, NKG2D and CD16, which can induce NK cell activation by initiating different signaling pathways.⁴ The NCR family is composed of three molecules: NKp30 (NCR3, CD337) and NKp44 (NCR2, CD336), which are expressed in humans, and NKp46 (NCR1, CD335), which is expressed in all mammals and is highly conserved between humans and mice. NKp46 is expressed mostly by NK cells and ILC1 cells but is also present on a small population of

T lymphocytes and a subset of ILC3 cells (NCR⁺ ILC3 cells) in mucosa.⁵⁻⁷

NK cells are present in blood and in primary and secondary lymphoid organs, including the spleen, bone marrow, liver, lymph nodes, lungs, tonsils, skin, uterus, and gut.⁸ The NK cell compartment consists of several subsets differing in maturation status. Human NK cells are defined as CD3⁻ CD14⁻ CD19⁻ CD56⁺ lymphocytes. They can be classified into subsets on the basis of the intensity of cell surface CD56 expression. Two principal NK cell subsets are found in the bloodstream in healthy individuals: CD56 dim CD16⁺ NK cells and CD56 bright CD16⁻ NK cells.⁹⁻¹¹ The phenotypic differences between these cells are associated with functional differences: CD56 bright CD16⁻ NK cells are less cytotoxic than CD56 dim CD16⁺ cells but produce greater amounts of cytokines in response to exposure to environmental cues, such as interleukin (IL)-12 and IL-18.⁹ We recently used high-throughput scRNAseq to investigate the heterogeneity of the NK cell compartment in the blood and spleen of humans and mice.¹² We identified two subsets conserved in human and mouse blood and spleen: NK1 and NK2 cells. In humans, the CD56 dim NK1 subset had high levels of cytotoxic activity, whereas the transcriptome of CD56 bright NK2 cells was enriched in gene

<footnote n="1">¹Aix Marseille University, CNRS, INSERM, Centre d'Immunologie de Marseille-Luminy, Marseille, France; ²CHU Bordeaux, Service d'Hématologie Clinique et de Thérapie Cellulaire, Bordeaux, France; ³Bordeaux University, Bordeaux, France; ⁴Institut National de la Santé et de la Recherche Médicale, U1035 Bordeaux, France; ⁵APHM, Hôpital de la Timone, Marseille-Immunopôle, Marseille, France; ⁶Établissement Français du Sang Nouvelle Aquitaine, Bordeaux, France and ⁷Innate Pharma Research Laboratories, Innate Pharma, Marseille, France</footnote>

<footnote>Correspondence: Émilie Narni-Mancinelli (narni@ciml.univ-mrs.fr) or Éric Vivier (vivier@ciml.univ-mrs.fr)</footnote>

<footnote>Statement of Significance: Single-cell profiling of bone marrow NK cells identified three subpopulations in all healthy donors: a CD56 bright subpopulation identified as the precursor of both of the other subpopulations, circulating CD56 dim NK1-like NK cells and CD56 bright NK2-like cells. Acute myeloid leukemia had a profound impact on NK cell heterogeneity, with stress-induced repression of effector functions.</footnote>
ontology (GO) terms relating to the secretion of cytokines and chemokines, consistent with previous findings. We also found two tissue-specific human NK subsets resident in the spleen, which we named hNK_Sp3 and hNK_Sp4. hNK_Sp3 cells resembled CD56<sup>bright</sup> NK2 cells, whereas hNK_Sp4 cells seemed to be more mature than hNK_Sp3 cells and more closely resembled the NK1 subset.

Several studies have suggested that CD56<sup>bright</sup> CD16<sup>-</sup> cells can develop into CD56<sup>dim</sup> CD16<sup>+</sup> NK cells. Indeed, CD56<sup>bright</sup> NK cells constitute the main NK cell population shortly after hematopoietic stem cell transplantation and subsequently disappear as CD56<sup>dim</sup> cell levels rise, 3 months after transplantation.<sup>13</sup> Purified CD56<sup>bright</sup> NK cells cocultured in vitro with synovial fibroblasts display downregulation of CD56 and a phenotype consistent with that of CD56<sup>dim</sup> NK cells.<sup>14</sup> In vivo, adoptive transfer of CD56<sup>bright</sup> NK cells into immunodeficient NOD-SCID recipient mice led to a CD56<sup>dim</sup>CD16<sup>+</sup> phenotype in the blood, spleen, and lymph nodes 10 days later.<sup>14</sup> Finally, CD56<sup>bright</sup> NK cells have been shown to have longer telomere repeats than CD56<sup>dim</sup> NK cells, suggesting that they are less mature than CD56<sup>dim</sup> NK cells.<sup>14</sup>

In both humans and mice, the expansion of a particular population of NK cells, called “adaptive NK” or “memory NK” cells, has been described during cytomegalovirus (CMV) infection.<sup>15</sup> These adaptive NK cells persist long after the contraction of the anti-CMV immune response and are capable of more powerful effector functions upon reinfection. Adaptive NK cells have also been described following the activation of NK cells by cytokines<sup>16</sup> or exposure to hantavirus infection,<sup>17</sup> herpes simplex virus 2 infection<sup>18</sup> or influenza vaccination in humans.<sup>19</sup> Adaptive NK cells express the cell surface activating receptor NKG2C in humans and also express the maturation marker CD57, and they lack many transcription factors and signaling proteins, including FCRγ, PLZF, Siglec-7, EAT-2, and SYK.<sup>20-22</sup>

Acute myeloid leukemia (AML) is a hematological malignancy characterized by bone marrow invasion of proliferative, clonal immature myeloid cells.<sup>23</sup> Blast cells replace normal hematopoietic cells and exploit the normal bone marrow microenvironment to survive and to support disease expansion, also affecting immune cells. NK cells from AML patients display low levels of cell surface NCR expression, which is associated with impaired cytotoxicity.<sup>24,25</sup> Conversely, higher levels of NKp30 and NKp46 on the cell surface have been linked to better outcomes.<sup>26,27</sup> Downregulation of other activating receptors, including DNAM-1, NKG2C and 2B4<sup>28-30</sup> and upregulation of inhibitory receptors, such as NKG2A,<sup>28,31</sup> on the surface of NK cells from AML patients also contribute to the low cytotoxicity and IFN-γ secretion levels of these cells.<sup>32</sup>

Here, we used unsupervised scRNAseq and unsupervised graph-based clustering to investigate human bone marrow NK cell heterogeneity under physiological conditions and at AML diagnosis. As previously found for the blood and spleen, healthy bone marrow contained NK1 and NK2 cells and an additional CD56<sup>bright</sup>-like resident NK cell subset, NK0 cells, from which NK1 and NK2 cells may differentiate. Comparison of transcriptomic profiles between NK cells from healthy and AML-derived bone marrow further revealed that AML profoundly modified the bone marrow NK cell compartment. Thus, AML affects NK cells at the pathological site.

## RESULTS

High-throughput scRNAseq identifies adaptive NK cells in human bone marrow

We investigated the heterogeneity of human bone marrow NK cells with the 10x Genomics high-throughput droplet-based scRNAseq pipeline, which can be used for unbiased transcriptomic characterization of thousands of cells simultaneously.<sup>33</sup> We analyzed ~24,000 flow cytometry-sorted NK cells, defined by their SSC/FSC profile and as CD3<sup>-</sup>CD14<sup>-</sup>CD19<sup>-</sup>CD45<sup>+</sup>CD56<sup>+</sup> cells,

from the bone marrow of eight healthy donors. We verified that the transcriptomic profile of the sorted cells contained the core NK cell signature previously defined<sup>12</sup> (data not shown). We confirmed that Lin<sup>-</sup>CD56<sup>+</sup> helper-like ILCs were absent from this population via module score analysis comparing single-cell human NK cell and helper-like ILC gene signatures, as previously described.<sup>10,34</sup> This evaluation resulted in the removal of only ten contaminating helper-like ILCs from the ~24,000 cells analyzed (Fig. 1A). We analyzed bone marrow NK cell heterogeneity by projecting cells into two dimensions by uniform manifold approximation and projection (UMAP) analysis. UMAP analysis revealed the clustering of human bone marrow NK cells into four different subsets, hereafter referred to as hNK_Bm1, hNK_Bm2, hNK_Bm3, and hNK_Bm4 (Fig. 1B). The three main NK cell subsets (hNK_Bm1, hNK_Bm2, and hNK_Bm3) were common to all samples, whereas the hNK_Bm4 subset was absent from some donors (Fig. 1C,D). We nevertheless checked the accuracy of hNK_Bm4 cell identification by a machine learning approach (data not shown).

We first focused on hNK_Bm4 cells and identified 139 differentially expressed genes distinguishing this cluster from the other human bone marrow NK cell subsets: 88 of these genes were downregulated and 51 were upregulated in hNK_Bm4 cells (Fig. 1E and Supplementary Table S1A). We found that compared with the other human bone marrow NK cell subsets, hNK_Bm4 cells had lower expression levels of the *KLRB1* and *CD7* genes, encoding activating receptors; *KLRC1* (encoding NKG2A) and *SIGLEC7*, encoding inhibitory receptors; *IL2RB*, encoding the IL-2/IL-15Rβ (CD122) surface receptor essential for the development and differentiation of NK cells;<sup>35</sup> and *FCER1G* (FCγR) (Fig. 1E). By contrast, hNK_Bm4 cells had higher expression levels of genes encoding NK cell-associated effector proteins (CCL5, GZMM, and GZMH)<sup>9,10</sup> and *ZEB2*, encoding a transcriptional regulator involved in terminal NK cell maturation.<sup>36</sup> GO term analysis revealed that the transcriptome of hNK_Bm4 cells was specifically enriched in NK cell effector functions, such as positive regulation of interferon (IFN)-γ production, NK cell-mediated cytotoxicity and granzyme-mediated apoptotic signaling (Fig. 1F). The hNK_Bm4 subset also displayed upregulation of several genes relating to the class II major histocompatibility complex (MHC), T cell costimulation, and the regulation of T cell activation and proliferation (Fig. 1F). The hNK_Bm4 subset was identified in three of the five donors seropositive for human cytomegalovirus (HCMV) (Supplementary Table S2A). *FCER1G* downregulation (Fig. 1E,G) and GO term enrichment for antigen presentation and T cell activation (Fig. 1E) are other properties exhibited by adaptive NK cells.<sup>37-40</sup> We therefore investigated whether the hNK_Bm4 subset could be considered to correspond to adaptive NK cells. We established canonical and adaptive NK cell gene signatures by analyzing the transcriptomic profiles of human blood CD3<sup>-</sup>CD56<sup>dim</sup>CD57<sup>+</sup>NKG2C<sup>+</sup> adaptive NK cells and CD3<sup>-</sup>CD56<sup>dim</sup>CD57<sup>-</sup>NKG2C<sup>-</sup> canonical NK cells from a public genome-wide bulk RNAseq dataset<sup>41</sup> (Supplementary Fig. S1). Unsupervised hierarchical clustering (Supplementary Fig. S1A) and principal component analysis (PCA) (Supplementary Fig. S1B) segregated the samples into two groups corresponding to adaptive and canonical NK cells. A bilateral comparison of the two populations identified 912 genes among the total gene set as differentially expressed (Supplementary Table S1B). These genes included *KLRC2* (encoding NKG2C) (Supplementary Fig. S1C), which was overexpressed in adaptive NK cells. The expression of this gene is associated with the adaptive phenotype of NK cells during cytomegalovirus infection.<sup>42</sup> The top 200 genes in both of these gene signatures intersected broadly with the genes detected in our scRNAseq dataset, which included 173 of the top 200 driver genes in the adaptive gene signature and 174 of the top 200 driver genes in the canonical gene signature (Supplementary Fig. S1D), indicating comparability. Module score analysis comparing the adaptive and
<figure><img src="image_3.png" /><figcaption>**Fig. 1** A Single-cell module score<sup>10,34</sup> for each of the bone marrow NK cell samples sorted from eight different healthy donors for the NK and ILC gene signatures defined by<sup>10</sup> and,<sup>8</sup> respectively. Cells are color-coded according to donor origin. **B** UMAP of healthy human bone marrow NK cells from eight donors. Cells are color-coded according to the defined subset. **C** UMAP of 23,684 healthy human bone marrow NK cells from eight donors, identifying four subsets of NK cells. Cells are color-coded according to donor origin. **D** Percentage of each healthy bone marrow NK cell subset across the eight donors. Subsets are color-coded as above. **E** Scatter plot of the genes differentially expressed between NK cells from the hNK_Bm4 cluster and NK cells from the other three clusters. Genes displaying significant differential expression (p < 0.05) are represented by yellow dots, and selected genes are highlighted. **F** Selected Gene Ontology terms displaying enrichment in hNK_Bm4 cells. The Benjamini–Hochberg-corrected −log₁₀ p value calculated by a hypergeometric test is reported. The black dotted line indicates the threshold of significance set at −log₁₀(0.05). **G** Feature plot showing the relative expression levels of the FCER1G gene in each of the 23,762 healthy human bone marrow NK cells. **H** Module score for the CD56<sup>dim</sup> CD57<sup>+</sup> NKG2C<sup>+</sup> adaptive and CD3<sup>−</sup> CD56<sup>dim</sup> CD57<sup>−</sup> NKG2C<sup>−</sup> canonical NK cell gene expression profiles obtained from reanalysis of the data of,<sup>41</sup> for each of the human healthy bone marrow NK cells, at the single-cell level. **I** Violin plots showing the distribution of module scores for CD56<sup>dim</sup> CD57<sup>+</sup> NKG2C<sup>+</sup> adaptive cells and for each blood NK cell grouped by subset. Statistical significance was assessed by the Kruskal–Wallis test with Dunn's multiple comparison test, and p values were adjusted by the Benjamini-Hochberg method. \*p value < 0.05, \*\*p value < 0.01, \*\*\*p value < 0.001, \*\*\*\*p value < 0.0001</figcaption></figure>
canonical gene signatures of each of the bone marrow NK cell subsets revealed that the gene expression profile of hNK_Bm4 cells most closely resembled that of adaptive NK cells (Fig. 1H). Violin plots of the module scores among the subsets also revealed significant enrichment of the gene signature of adaptive NK cells in hNK_Bm4 cells (Fig. 1I). The hNK_Bm4 subset therefore had a transcriptomic profile different from that of conventional NK cells but closely resembling that of adaptive NK cells, suggesting that adaptive NK cells may populate the bone marrow of some HCMV-seropositive individuals.

Three NK cell subsets are present in human bone marrow under physiological conditions

For further characterization of the heterogeneity of bone marrow NK cells among donors, we removed the hNK_Bm4 adaptive NK cell-like NK cell subset, which was absent from some donors, from the analysis. Unsupervised hierarchical clustering (data not shown) and UMAP analysis of the remaining ~23,000 human bone marrow NK cells again clustered the cells into three different subsets: hNK_Bm1, hNK_Bm2, and hNK_Bm3 (Fig. 2A). These three subsets were present in all samples (Fig. 1C, D). Significant differential expression between the hNK_Bm1, hNK_Bm2, and hNK_Bm3 subsets was observed for 255 of the genes composing the three different transcriptomic signatures of bone marrow NK cells (Fig. 2B, Supplementary Table S1C). No significant individual-specific phenotype or batch effect was observed in the unsupervised hierarchical clustering analysis (data not shown), and the NK subsets from a given individual did not cluster together on the PCA score plots (Fig. 2C). Instead, PC1 and PC2, accounting for 42% and 22% of the variance, respectively, and unsupervised hierarchical clustering (data not shown) segregated the eight donors into three different groups corresponding to the hNK_Bm1, hNK_Bm2, and hNK_Bm3 subsets (Fig. 2C). PC1 separated the hNK_Bm1, hNK_Bm2, and hNK_Bm3 subsets. The hNK_Bm3 subset was also separated from the hNK_Bm1 and hNK_Bm2 subsets by PC2. The GNLY and GZMB genes associated with cell cytotoxicity were identified as drivers for hNK_Bm1 along PC1. CD160, encoding an NK cell activating receptor, and CCL3 and CCL4, encoding chemokines, were identified as drivers of the hNK_Bm2 subset, whereas ZFP36L2 (which is involved in lymphocyte hematopoiesis43,44) and DUSP2 (a regulator of the ERK pathway) were drivers of the hNK_Bm3 subset along PC1 and PC2 (Fig. 2D).

We then analyzed the top ten genes encoding secreted proteins, cell membrane markers and transcription factors expressed in hNK_Bm1, hNK_Bm2, and hNK_Bm3 cells (Fig. 2E). hNK_Bm1 cells displayed upregulation of genes associated with cytotoxic factors (FGFBP2, GZMB, GZMH, and PRF1) and overexpressed FCGR3A (encoding CD16), which is characteristic of the CD56<sup>dim</sup> NK cell subset;38-40 moreover, these cells overexpressed S1PR5, which encodes a sphingosine 1-phosphate receptor promoting NK cell egress from the lymph nodes and bone marrow45 (Fig. 2E). The hNK_Bm2 subset transcriptome displayed enrichment with genes encoding soluble factors associated with NK cell effector functions (CCL3, CCL4, XCL1, XCL2, GZMK, and CCL3L1), amphiregulin (AREG), the activating receptor encoded by CD160, the activation marker encoded by CD69, and the chemokine receptor encoded by CXCR6, which is involved in ILC egress from the bone marrow46 (Fig. 2E). Like hNK_Bm2 cells, the hNK_Bm3 subset displayed higher expression levels of GZMK, XCL1, XCL2, and AREG, but these cells also displayed upregulation of the LTB, SELL (encoding CD62L) and CD44 genes encoding adhesion molecules associated with homing and anchorage in the bone marrow47,48 (Fig. 2E). Biological process enrichment analysis revealed specific enrichment of the hNK_Bm1 subset transcriptome in the FcγR signaling pathway, the granzyme-mediated apoptotic signaling pathway and cytolysis (Fig. 2F). By contrast, the transcriptomes of both the hNK_Bm2 and hNK_Bm3 subsets

displayed enrichment in chemokine and cytokine responses, cytokine production and cytokine-mediated signaling. The hNK_Bm2 subset transcriptome was more strongly enriched in responses, and the hNK_Bm3 subset transcriptome was more specifically enriched in IL-2 production. These results thus reveal the presence of three different subsets of human bone marrow NK cells under physiological conditions.

The three human bone marrow NK cell subsets encompass an NK1-CD56<sup>dim</sup>-like subset, an NK2-CD56<sup>bright</sup>-like subset and an additional CD56<sup>bright</sup>-like subset

We then compared the transcriptomic signatures of the three human bone marrow NK cell subsets with those of the blood and splenic NK1 and NK2 subsets, which resembled CD56<sup>dim</sup> and CD56<sup>bright</sup> cells, respectively, and two additional spleen-specific hNK_Sp3 and hNK_Sp4 human NK cell populations, as previously described.<sup>10,12</sup> Gene signature heatmaps (Fig. 3A–C) and module score analysis (Supplementary Fig. S2A–C) revealed that hNK_Bm1 cells corresponded to NK1, the related splenic hNK_Sp1 and CD56<sup>dim</sup> NK cells, consistent with the pattern of FCGR3A expression. hNK_Bm2 corresponded to NK2, the related hNK_Sp2 and CD56<sup>bright</sup> NK cells (Fig. 3A, B and Supplementary Fig. S2A–C). Further analysis revealed that the hNK_Bm3 subset shared part of the NK2 gene signature but was most similar to the minor hNK_Sp3 subset residing in the spleen (Fig. 3A, C and Supplementary Fig. S2C). No hNK_Sp4 subset signature was found in bone marrow NK cells (Fig. 3C).

Flow cytometric analysis of bone marrow NK cells revealed that the expression of CD160 and CD52 (CAMPATH-1) on the surface of CD56<sup>bright</sup> cells was mutually exclusive, further discriminating between hNK_Bm2 and hNK_Bm3, respectively, as previously reported for hNK_Sp2 and hNK_Sp3 in the spleen<sup>12</sup> (Fig. 3D). These analyses showed that the human bone marrow NK cell compartment consists mostly of CD56<sup>dim</sup> NK cells (~50%), with a few CD56<sup>bright</sup>-like CD160<sup>+</sup> CD52<sup>-</sup> NK cells (~10%) and CD56<sup>bright</sup> CD160<sup>-</sup> CD52<sup>+</sup> NK cells (~5%) (Fig. 3E). The protein levels of perforin, KIR2DL2/2DL3 and KIR2DL1/2DS1 were higher in CD56<sup>dim</sup>-like cells than in CD56<sup>bright</sup>-like CD160<sup>+</sup> CD52<sup>-</sup> NK cells and CD56<sup>bright</sup> CD160<sup>-</sup> CD52<sup>+</sup> NK cells, confirming our transcriptomic data (Fig. 3F, G). The hNK_Bm2 subset was characterized by higher levels of CD122, CD27, and CD69, whereas hNK_Bm3 cells had higher levels of CD127 (Fig. 3F, G).

Together, these results demonstrate that the hNK_Bm1 and hNK_Bm2 subsets resemble the NK1-CD56<sup>dim</sup> and NK2-CD56<sup>bright</sup> subsets, respectively, which also reside in the blood and spleen in humans. The third subset, hNK_Bm3, most closely resembled the hNK_Sp3 population resident in the spleen and absent from the bloodstream.

The NK0 subset has features consistent with a role as the precursor of both the NK2 and NK1 subsets

We then investigated whether the four NK cell subsets identified in the bone marrow reflected differences in maturation status. KLRC1 was expressed at higher levels in the hNK_Bm2 and hNK_Bm3 subsets, whereas expression of the KIR2DL3 and KIR3DL1 killer cell immunoglobulin-like receptor (KIR) genes was restricted to the hNK_Bm1 subset (Fig. 4A–C). These results suggest that the hNK_Bm1 subset is more mature than the hNK_Bm2 and hNK_Bm3 subsets because a decrease in cell surface expression of NKG2A and a concomitant increase in KIR expression have been associated with human NK cell differentiation.<sup>49</sup> We then searched for a developmental path across bone marrow NK subsets by applying the pseudotime algorithm Monocle DDRTree, which computationally orders the transcriptomic profiles of cells along a trajectory without prior information about their clustering.<sup>50</sup> The pseudotime algorithm ordered the cells along a trajectory segregating our ~24,000 bone marrow NK cells into two branches and four clusters. These four clusters
A.

<figure><img src="image_4.png" /></figure>

C.

<figure><img src="image_5.png" /></figure>

E.

<table><thead><tr><th></th><th>Top 10</th><th>Secreted</th><th>Cell membrane</th><th>Transcription factors</th></tr></thead><tbody><tr><td>hNK_Bm1</td><td>FGFBP2<br>GZMB<br>SPON2<br>FCGR3A<br>LGALS1<br>GZMH<br>PRF1<br>LAIR2<br>S100A4<br>IGFBP7</td><td>FGFBP2<br>GZMB<br>SPON2<br>LGALS1<br>GZMH<br>PRF1<br>LAIR2<br>IGFBP7<br>PRSS23<br>GNLY</td><td>FCGR3A<br>S100A6<br>SYNE2<br>EMP3<br>CD247<br>CD99<br>KIR2DL3<br>S1PR5<br>KIR2DL1<br>FCRL6</td><td>HOPX<br>GTF3C1</td></tr><tr><td>hNK_Bm2</td><td>CCL3<br>XCL1<br>GZMK<br>XLC2<br>ZNF331<br>TOX2<br>PIK3R1<br>NR4A2<br>CD160<br>AC245014.3</td><td>CCL3<br>XCL1<br>GZMK<br>XCL2<br>CCL3L1<br>SERPINE2<br>FAM177A1<br>ABHD15<br>AREG<br>LEPROTL1</td><td>CD160<br>IL2RB<br>CD27<br>CXCR6<br>CD69<br>ADGRG3<br>SPRY1<br>ALOX5AP<br>TIGIT<br>IFRD1</td><td>ZNF331<br>TOX2<br>NR4A2<br>FOS<br>IRF8<br>TSC22D3<br>KLF6<br>ELF1<br>JUND<br>STAT4</td></tr><tr><td>hNK_Bm3</td><td>SELL<br>GZMK<br>LTB<br>CD44<br>IL7R<br>COTL1<br>CXCR3<br>CAPG<br>XCL1<br>SPTSSB</td><td>GZMK<br>LTB<br>CAPG<br>XCL1<br>SPTSSB<br>XCL2<br>AREG</td><td>SELL<br>CD44<br>IL7R<br>CXCR3<br>SPTSSB<br>GPR183<br>CD52<br>CD3E<br>FXYD5<br>TMEM123</td><td>TCF7<br>BHLHE40<br>FOS<br>MAFF<br>JUN</td></tr></tbody></table>

B.

<figure><img src="image_6.png" /></figure>

D.

<figure><img src="image_7.png" /></figure>

F.

<figure><img src="image_8.png" /></figure>

**Fig. 2** A UMAP of 23,033 healthy human bone marrow NK cells from eight donors, without hNK_Bm4, identifying three subsets of NK cells. Cells are color-coded according to the defined subset. B Heatmap of the 255 genes tested with the Wilcoxon rank sum test separating the 23,033 healthy human bone marrow NK cells into subsets. Cells are plotted in columns according to subset source. Genes are shown in rows and ranked by adjusted *p* values < 0.05. Gene expression is color-coded on a scale based on the z-score distribution, from −2.5 (purple) to 2.5 (yellow). The squares indicate specific transcriptomic signatures of human bone marrow NK cell subsets. C Principal component analysis on the three healthy human bone marrow NK cell subsets in each sample based on the mean expression levels of genes with differential expression. D Driver genes for each cell subset, accounting for 20% of the total variance for each PC. E Top ten most highly expressed genes among the total gene set and among genes encoding secreted proteins, cell membrane markers and transcription factors with significant differences among the three healthy human bone marrow NK cell subsets. Gene symbols and annotations were retrieved from public databases. Cell membrane protein-encoding genes are color-coded in orange, transcription factor-encoding genes in dark blue, secreted protein-encoding genes in red and genes encoding other proteins in gray. Genes are ranked by *p* value. F Selected Gene Ontology terms displaying enrichment in the three healthy human bone marrow NK cell subsets. The Benjamini–Hochberg-corrected −log₁₀ *p* value calculated by a hypergeometric test is reported. The black dotted line indicates the threshold of significance set at −log₁₀(0.05)
<figure><img src="image_9.png" /><figcaption>**Fig. 3** A Heatmap showing the hNK1 and hNK2 gene expression profiles of each of the three healthy human bone marrow NK cell subsets. B Heatmap showing the hNK_Sp1, hNK_Sp2, hNK_Sp3, and hNK_Sp4 gene expression profiles of each of the three healthy human bone marrow NK cell subsets. C Flow cytometric analysis of healthy bone marrow NK cells, showing CD52, CD160, and CD56 cell surface expression in one representative healthy donor. D Frequencies of the indicated cell populations among total bone marrow NK cells from five healthy donors. E FACS profiles of the indicated cell surface markers in the CD56<sup>dim</sup>, CD56<sup>bright</sup> CD160<sup>+</sup> CD52<sup>+</sup>, and CD56<sup>bright</sup> CD160<sup>-</sup> CD52<sup>+</sup> cell subsets in one representative healthy donor. F Mean fluorescence intensity of the cell surface markers shown in E for each of the five healthy donors. Statistical significance was assessed by Friedmann's analysis with Dunn's post hoc test on paired subset data, and *p* values were adjusted by the Benjamini-Hochberg method. The error bars indicate the mean (±SD) values. \* *p* value < 0.05, \*\* *p* value < 0.01, \*\*\* *p* value < 0.001, \*\*\*\* *p* value < 0.0001</figcaption></figure>
A.

<figure><img src="image_10.png" /></figure>

C.

<figure><img src="image_11.png" /></figure>

D.

<figure><img src="image_12.png" /></figure>

E.

<figure><img src="image_13.png" /></figure>

H.

<figure><img src="image_14.png" /></figure>

I.

<figure><img src="image_15.png" /></figure>

B.

<figure><img src="image_16.png" /></figure>

F.

<figure><img src="image_17.png" /></figure>

G.

<figure><img src="image_18.png" /></figure>

J.

<figure><img src="image_19.png" /></figure>

K.

<figure><img src="image_20.png" /></figure>

**Fig. 4** A Feature plot of the relative expression level of the *KIR2DL1* gene in each of the 23,066 healthy human bone marrow NK cells. **B** Feature plot of the relative expression level of the *KIR2DL3* gene in each of the 23,066 healthy human bone marrow NK cells. **C** Feature plot of the relative expression level of the *KLRC1* gene in each of the 23,066 healthy human bone marrow NK cells. **D** Pseudotime inference of the bone marrow NK cell trajectory in one representative healthy donor from the eight studied donors, with cells color-coded according to the corresponding bone marrow NK cell subset identified by UMAP. **E** Module score for the CD34<sup>+</sup> CD38<sup>+</sup> CD123<sup>-</sup> CD45RA<sup>+</sup> CD7<sup>+</sup> CD10<sup>+</sup> CD127<sup>-</sup> bone marrow NKP gene expression profile obtained by reanalysis of the data of<sup>51</sup> for one representative healthy donor by the pseudotime algorithm. **F** Module score for the CD3<sup>-</sup> CD56<sup>+</sup> NKp46<sup>+</sup> bone marrow mature NK cell gene expression profile obtained by reanalysis of the data of<sup>51</sup> for one representative healthy donor by the pseudotime algorithm. **G** Violin plot showing the distribution of the module scores for CD3<sup>-</sup> CD56<sup>+</sup> NKp46<sup>+</sup> mature NK cells among each human healthy bone marrow NK cell subset, grouped by subset. **H** Pseudotime inference by the pseudotime algorithm from an starting point selected as the most distal cell of the hNK_Bm3 subset in one healthy donor representative of the eight studied donors. The scale indicates the maturation state, from dark blue (least mature) to yellow (most mature). **I** Pseudotime inference of the splenic NK cell trajectory in one healthy donor representative of the three studied donors, with cells color-coded according to the corresponding splenic NK subset identified by UMAP, as previously described by.<sup>12</sup> **J** Violin plot showing the distribution of the module scores for CD3<sup>-</sup> CD56<sup>+</sup> NKp46<sup>+</sup> mature NK cells among each human healthy splenic NK cell subset, grouped by subset. **K** Pseudotime inference by the pseudotime algorithm from an starting point selected as the most distal cell of the hNK_Sp3 subset in one healthy donor representative of the three studied donors. The scale indicates the maturation state, from dark blue (least mature) to yellow (most mature)
overlapped perfectly with the hNK_Bm1, hNK_Bm2, hNK_Bm3, and hNK_Bm4 subsets defined by UMAP, providing additional support for our previous conclusions. The hNK_Bm3 subset was located at the intersection of two branches linking the hNK_Bm1-hNK_Bm4-hNK_Bm3 and hNK_Bm2-hNK_Bm3 subsets (Fig. 4D).

We then tried to define the starting point of the putative developmental trajectory by comparing our dataset with an NK cell precursor (NKP) gene signature. We examined the transcriptomic profiles of bone marrow Lin⁻CD34⁺CD38⁺CD123⁺CD45RA⁺CD7⁺CD10⁺CD127⁻NKP and CD3⁻CD56⁺NKp46⁺mature NK cells from a public genome-wide RNAseq dataset.⁵¹ Unsupervised hierarchical clustering separated the samples into two main branches on the basis of a precursor or a mature NK cell phenotype, suggesting that these two populations can be differentiated on the basis of a specific transcriptomic pattern (Supplementary Fig. S3A). PCA confirmed this finding, as mature NK cells clustered away from progenitors along a combined PC1 and PC2 axis, accounting for 30.7% and 16.6% of the total variance, respectively (Supplementary Fig. S3B). Bilateral comparison of the two populations identified 448 genes from the total gene set as differentially expressed, with 265 upregulated in NK cells and 183 upregulated in progenitors (Supplementary Fig. S3C, Supplementary Table S1D). Mature NK cells had increased levels of CD160, FCGR3A, GZMB, and CCL3 expression, whereas NK cell progenitors displayed upregulation of TCF4 and RAG2, which are expressed during NK cell development⁵² (Supplementary Fig. S3C). These gene signatures broadly intersected with the genes detected by our scRNAseq analysis, with 263 genes from the mature NK cell signature and 160 genes from the progenitor gene signature found in our dataset (Supplementary Fig. S3D). Module score analyses comparing the NKP and mature NK cell gene signatures among bone marrow NK cell subsets were performed on the developmental trajectory defined by the pseudotime algorithm (Fig. 4E, F). These analyses revealed that the hNK_Bm3 subset displayed specific enrichment for the NKP signature, whereas the hNK_Bm1 and hNK_Bm4 subsets were enriched for the mature NK cell profile (Fig. 4E, F). Violin plots of the NK gene signature module scores also showed that the hNK_Bm3 subset was the least mature (Fig. 4G). The hNK_Bm3 subset was, therefore, defined as the starting point of the pseudotime developmental trajectory (Fig. 4H).

Pseudotime analysis with Monocle DDRTree revealed a specific developmental trajectory of human bone marrow NK cells, with the CD56<sup>bright</sup>-like hNK_Bm3 minor subset differentiating into the NK1-CD56<sup>dim</sup>-like hNK_Bm1 subset and a distinct NK2-CD56<sup>bright</sup>-like hNK_Bm2 subset (Fig. 4H). The NK1-CD56<sup>dim</sup>-like hNK_Bm1 and NK2-CD56<sup>bright</sup>-like hNK_Bm2 subsets followed different developmental pathways. hNK_Bm4, the adaptive NK cell subset, appeared to arise at a later time point from the hNK_Bm1 subset (Fig. 4H). Similar results were obtained with another pseudotime algorithm (PAGA Tree from Dynverse), providing additional support for our putative differentiation trajectory (data not shown).

In addition to the bone marrow, we reanalyzed our previous splenic dataset¹² with the recently developed pseudotime algorithm Monocle DDRTree. Pseudotime analysis ordered the splenic NK cells into two branches and four clusters corresponding to the hNK_Sp1, hNK_Sp2, hNK_Sp3, and hNK_Sp4 subsets (Fig. 4I). The hNK_Sp3 subset was identified as the least mature, with the hNK_Sp2, hNK_Sp1, and hNK_Sp4 subsets sequentially more mature (Fig. 4J). Thus, in the spleen, as in the bone marrow, a minor CD56<sup>bright</sup>-like hNK_Sp3 subset differentiates into the NK1-CD56<sup>dim</sup>-like hNK_Sp1 subset and a distinct NK2-CD56<sup>bright</sup>-like hNK_Sp2 subset (Fig. 4K). Thus, hereafter, we refer to the CD56<sup>bright</sup>-like hNK_Bm3 and CD56<sup>bright</sup>-like hNK_Sp3 subsets as NK0 (Supplementary Fig. S4). The hNK_Sp4 subset, an additional minor CD56<sup>dim</sup> splenic NK cell subset for which no corresponding subset was found in the bone marrow, seemed to differentiate from the NK1-CD56<sup>dim</sup>-like hNK_Sp1 subset.

Acute myeloid leukemia affects the transcriptome of human bone marrow NK cells

Acute myeloid leukemia (AML) is a specific bone marrow disease known to be associated with strong impairment of NK cell function in the periphery.²⁷ We therefore investigated the impact of AML on the heterogeneity of bone marrow NK cells. We analyzed ~15,000 flow cytometry-sorted bone marrow NK cells, defined as CD3⁻CD14⁻CD19⁻CD34⁻CD56⁺ cells, obtained from eight AML patients at diagnosis (Supplementary Table S2B). Module score analysis comparing human NK cell and helper-like ILC gene signatures resulted in the removal of 39 ILC-like cells (Supplementary Fig. S5A). UMAP analysis on the remaining cells revealed the existence of several NK cell subsets (Fig. 5A), but their assignment on the basis of donor origin showed that these clusters had a strong donor phenotype (Fig. 5B). There was no conserved NK cell subset in the AML patients (data not shown). Moreover, there was no conserved NK cell subset when AML patients were considered according to the same French-American-British classification of AML (data not shown). We then investigated whether module score analysis could be used to compare the specific signatures of each of the hNK_BM subsets to this AML dataset (Supplementary Fig. S5B-E). Module score analysis was performed on the UMAP representation of total NK cells from the AML patients (Supplementary Fig. S5B-E). A high level of transcriptomic heterogeneity was observed at the patient level for both subset composition and percentage. Based on protein levels, we detected heterogeneity in the percentages of the CD56<sup>bright</sup> and CD56<sup>dim</sup> subsets among the AML samples (Fig. 5C). Thus, AML strongly affected the transcriptomic profile of bone marrow NK cells, making it impossible to analyze their heterogeneity with classical annotation tools.

We therefore further investigated the effects of AML on bone marrow NK cells by analyzing a mixture of bone marrow NK cells from healthy donors and AML patients. UMAP analysis of the ~38,700 bone marrow NK cells revealed that NK cells clustered differentially according to the clinical status of the donor (Fig. 5D). In total, 197 genes were differentially expressed between bone marrow NK cells from healthy donors and those from donors with AML: 90 genes were upregulated in NK cells from healthy donors, and 107 genes were upregulated in NK cells from donors with AML (Figs. 5E, Supplementary Figs. S6 and S7 and Table S1E). Bone marrow NK cells from AML patients had higher expression levels of interferon-induced genes (IFI44L, IFI6, IFIT3, and IFI44), HLA molecule-encoding genes (HLA-DPB1, HLA-DPA1, HLA-DRB5, and HLA-DRB1), ZEB2, and KLF2, which encode transcription factors controlling NK cell maturation³⁶ and survival,⁵³ respectively. By contrast, NK cells from the bone marrow of healthy donors displayed higher expression levels of genes encoding NK cell effector molecules (CCL3, CCL4 and GZMM) and genes encoding surface receptors (CD16 (FCGR3A), CD161 (KLRB1), NKG2A (KLRC1), CD300A, FcγR (FCER1G), KIR2DL3, KIR2DL1, and CD160) (Fig. 5E and Supplementary Fig. S6 and Table S1E). The higher level of CD160 expression in bone marrow NK cells from healthy donors than in bone marrow cells from donors with donors with AML was confirmed by flow cytometric assessment of protein levels (Fig. 5F). We then further investigated the potential role of CD160 in AML pathogenesis and progression by studying the clinical outcomes of cancer patients in the TCGA database. The survival rates of patients with CD160<sup>high</sup> AML were much higher than those of patients with CD160<sup>low</sup> cancer, suggesting that CD160 is an antitumor biomarker in AML (Fig. 5G). In addition, bone marrow NK cells from healthy donors displayed stronger expression of genes involved in regulating NK cell effector functions (DOCK2, EVL, LCK, RAP1B, LAT, JAK1, CORO1A, and FGR<sup>54-61</sup>), PTPN6 (encoding SHIP-1), and ITGB2, which has been reported to promote the NK cell response<sup>62,63</sup> (Fig. 5E and Supplementary Fig. S7 and Table S1E).
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>**B.** Single-cell profiling reveals the trajectories of natural killer cell...
A Crinier et al.
15,005 cells - 1,030 mean genes/cell</figcaption></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

9

**Fig. 5** A UMAP of 15,005 human bone marrow NK cells from eight donors with AML. Cells are color-coded according to the defined subset. B UMAP of 15,005 human bone marrow NK cells from eight donors with AML. Cells are color-coded according to donor origin. C Frequencies of CD56^bright^ and CD56^dim^ cells among CD3^−^CD56^+^ or CD3^−^CD14^−^CD19^−^CD56^+^ cells from seven healthy donors and 22 donors with AML from the pooled data of two independent experiments. D UMAP of 38,691 human bone marrow NK cells from eight healthy donors and eight donors with AML. Cells are color-coded according to the donor's clinical status. E Scatter plot of the genes differentially expressed between healthy and AML NK cells. Genes displaying significant differential expression (*p* value < 0.05) are represented by yellow dots, and selected genes are highlighted. F Frequencies of CD160-expressing cells among the total CD3^−^CD14^−^CD19^−^CD56^+^ NK cell population from five healthy donors and seven donors with AML. Statistical significance was assessed by Friedmann's analysis with Dunn's post hoc tests on paired subset data, and *p* values were adjusted by the Benjamini-Hochberg method. The error bars indicate the mean (±SD) values. \* *p* value < 0.05, \*\* *p* value < 0.01, \*\*\* *p* value < 0.001, \*\*\*\* *p* value < 0.0001. G Kaplan-Meier overall survival curves for AML patients from TCGA stratified by the CD160 expression level. The optimal cutoff for patient stratification was obtained with a Cox proportional hazards model, and the *p* value was calculated by the log-rank test. CD160^high^ group (≥7.427, *n* = 78); CD160^low^ group (<7.427, *n* = 83)
Biological process enrichment analysis revealed enrichment in the GO terms "response to cytokine" and "type I interferon signaling pathways" in bone marrow NK cells from AML patients (Fig. 6A) relative to bone marrow NK cells isolated from healthy individuals. By contrast, NK cells purified from healthy bone

marrow displayed enrichment in the GO terms "NK cell-mediated cytotoxicity", "FcγR signaling pathway", "exocytosis", "response to IL-12" and "IL-15-mediated signaling pathways". Collectively, these results suggest that bone marrow NK cells from healthy donors have a more activated phenotype than those isolated from AML

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

**Fig. 6** A Selected Gene Ontology terms displaying enrichment in healthy and AML samples. The Benjamini–Hochberg-corrected –log₁₀ *p* value calculated by a hypergeometric test is reported. The black dotted line indicates the threshold of significance set at –log₁₀(0.05). B Frequencies of CD107⁺ and IFN-γ⁺ cells among CD3⁻ CD56⁺ NK cells from five healthy donors and five donors with AML from the pooled data of two independent experiments. Statistical significance was assessed by Friedmann's analysis with Dunn's post hoc tests on paired subset data, and *p* values were adjusted by the Benjamini–Hochberg method. The error bars indicate the mean (±SD) values. \* *p* value < 0.05, \*\* *p* value < 0.01, \*\*\* *p* value < 0.001, \*\*\*\* *p* value < 0.0001. C Heatmap of selected Gene Ontology terms related to gene expression, the stress response and metabolism. The Benjamini–Hochberg-corrected –log₁₀ *p* values calculated by a hypergeometric test are reported on a scale from 0 to 6. The red star indicates the threshold of significance set at –log₁₀(0.05)
patients, which displayed a strong type I IFN response signature. Indeed, upon simulation with phorbol myristate acetate and ionomycin, higher frequencies of IFN-γ-producing NK cells were detected in healthy donors than in AML patients, consistent with the above findings (Fig. 6B). We extended our analysis of NK cell transcriptomic profiles further by assessing the enrichment in GO terms related to cell metabolism (Fig. 6C). Relative to bone marrow NK cells from healthy individuals, bone marrow NK cells from AML patients displayed enrichment in GO terms associated with the regulation of gene expression through the processes of transcription, splicing, and epigenetic modifications; the regulation of macromolecules; and metabolic processes and stress-related signaling. These results highlight the impact of AML on bone marrow NK cells.

## DISCUSSION

The main objective of this study was to address, from an unbiased transcriptome-wide perspective, the heterogeneity of NK cells in human bone marrow under physiological conditions and at the time of diagnosis of acute myeloid leukemia, a hematologic malignancy of the bone marrow.

We report here four different clusters of Lin⁻CD56⁺ cells in bone marrow—hNK_Bm1, hNK_Bm2, hNK_Bm3, and hNK_Bm4— the last of which was not present in all donors. The unsupervised hierarchical clustering results were further validated by machine learning algorithms, supporting the robustness of our analysis. By contrast, scRNAseq profiling of Lin⁻CD7⁺ cells from healthy human bone marrow identified six NK cell clusters,⁶⁴ with CD7 expressed by both progenitor and mature ILCs.⁶⁵ Our analysis of Lin⁻CD56⁺ cells was thus expected to yield different results. However, similarities may exist between the two studies, as the ‘mature’ and ‘terminal’ populations characterized by PRF1, GZMB, S100A6, FCGR3A, LGALS1, and SPON2 expression⁶⁴ resembled the hNK_Bm1 subset; the “active” and “transitional” populations expressing CCL3, XCL2, PI3KR1, IRFD1, and FOS⁶⁴ resembled the hNK_Bm2 subset; and the CD56bright subset⁶⁴ resembled the CD56bright hNK_Bm3 subset, with the expression of IL7R, CD44, LTB, XCL1, GZMK, SELL, and XCL2.

We identified three subsets of NK cells, including the conserved NK1-CD56dim-like and NK-CD56bright-like subsets, in the bone marrow of all healthy humans. Both the NK1 and NK2 subsets are, therefore, present in the blood, spleen, and bone marrow. This finding is consistent with the origination of NK cells in the bone marrow and with their release into the bloodstream and continual recirculation throughout the body, leading to the population of several organs, including the spleen.⁹,¹¹ Consistent with our previous findings,¹² we observed that the NK1-CD56dim subset had a transcriptomic profile enriched in markers of cytotoxic activity, whereas the transcriptomic profile of NK2-CD56bright cells was more heavily enriched in genes involved in the response to cytokines and chemokines. Preliminary data indicated that NK1 and NK2 cells were also present in nonlymphoid organs, such as the duodenum and lungs, but these findings require further exploration.

We detected an additional NK0 CD56bright-like subset, originally identified as hNK_Bm3, which resembled the hNK_Sp3 subset residing in the human spleen.¹² NK0 cells were not found in the blood, explaining the previous lack of detection of this subset. The cells of this population expressed high surface levels of CD56 (CD56bright) and had a transcriptomic signature similar to that of the CD56bright NK subset but could be distinguished from NK2 cells on the basis of their higher levels of CD52 and CD127 expression and lower levels of CD160 expression. Additional experiments are necessary to clarify the real relationship between the hNK_Bm3 subset and the conventional CD56bright NK subset. We found that the transcriptome of NK0 cells displayed specific enrichment with genes associated with the NK cell precursor (NKP)

signature. The pseudotime algorithm Monocle DDRTree, which infers developmental trajectories from scRNAseq data, showed that bone marrow NK0 (hNK_Bm3) cells can differentiate into both NK1 and NK2 cells (Supplementary Fig. S4). Our data thus show that a minor subset of the CD56bright NK cell population, the subpopulation of tissue-resident CD56bright CD127⁺ CD160⁻ CD52⁺ cells (referred to hereafter as NK0 cells) can give rise to other NK2/CD56bright CD160⁺ CD52⁻ cells and to NK1/CD56dim- perforin⁺high cells under physiological conditions. These results are consistent with those of previous studies showing that CD56bright NK cells can differentiate into CD56dim NK cells¹³ and provide the first demonstration, at the transcriptomic level, of such differentiation in humans in vivo. NK1 and NK2 cells appeared at the end of the two branches on the pseudotime trajectory, suggesting that these two populations are not developmentally related. This observation may explain why not all CD56bright cells differentiate into CD56dim cells.¹⁴ However, upon cytokine-induced activation,⁶⁶ coculture with synovial fibroblasts or transfer into NOD/SCID mice,¹⁴ CD56bright NK2-like cells can differentiate into CD56dim NK cells. These results thus suggest that under certain circumstances, circulating CD56bright NK2-like NK cells in peripheral blood can also give rise to CD56dim-like cells. However, it remains to be determined whether this phenomenon naturally occurs without external manipulation and to what extent the CD56dim-like subset generated resembles the conventional CD56dim subset.

Similar analyses of the human spleen showed that splenic NK0 (hNK_Sp3) cells can also differentiate into both NK1 and NK2 cells. It was long supposed that NK cells develop exclusively within the bone marrow, but it is now clear that there is enrichment of some NK cell precursors in extramedullary tissues.⁶⁷ The biological relevance of extramedullary NK cell differentiation and the mechanisms governing this process remains to be determined.

In this study, we also used public datasets containing transcriptomic profiles of canonical and adaptive NK cells⁵¹ to define the first transcriptomic signature of human adaptive NK cells. Based on transcriptomic data, we detected adaptive NK cells in the bone marrow of human CMV-seropositive individuals. The presence of adaptive NK cells was correlated with CMV seropositivity, but not all CMV-seropositive individuals have circulating adaptive NK cells.⁴⁰ Thus, adaptive NK cells appear to be present in the bone marrow, in addition to the blood and spleen.

A population of tissue-resident NK cells, ItNK cells, has recently been described.⁶⁸ These cells are absent from the blood but can account for ~29% of NK cells in the bone marrow, ~45% of splenic NK cells and ~56% of NK cells in lymph nodes.⁶⁸ ItNK cells have a CD69⁺ CXCR6⁺ cell surface phenotype but disparate levels of CD56 and CD16,⁶⁸ suggesting that this population, at least the CD56bright and CD56dim subsets, is heterogeneous. The transcriptomes of the hNK_Bm3 and hNK_Sp3 subsets, corresponding to our NK0 tissue-resident progenitors, were not enriched in the transcriptomic signature of ItNK cells defined by bulk RNAseq on flow cytometry-sorted CD69⁺ CXCR6⁺ cells.⁶⁹ Given that in RNAseq, the gene expression profile of the considered population is averaged and that ItNK cells have a heterogeneous surface phenotype, these cells may also have a heterogeneous transcriptomic profile, preventing their identification at the transcriptomic level by scRNAseq.

We tried to explore the heterogeneity of the bone marrow NK cell compartment in AML patients at diagnosis. We found that each sample had a unique specific transcriptomic profile, preventing subset assignment. A similar patient-specific pattern was also recently observed for the expression of NKp30 and CD27 on the surface of NK cells from AML patients.⁷⁰ We were
nevertheless able to extract a transcriptomic signature of NK cells under AML disease conditions relative to physiological conditions. Bone marrow NK cells from AML patients had a transcriptomic profile enriched in genes involved in responses to cytokine and type I interferon signaling pathways. Consistent with the well-known impairment of NK cell function during AML<sup>24,25,27</sup> NK cells isolated from healthy individuals had a transcriptomic profile enriched in genes involved in cell cytotoxicity, the FcγR receptor signaling pathway, exocytosis, and the response to IL-12 and IL-15 signaling pathways, indicating a more activated NK cell phenotype for these cells than for those isolated from AML patients. Indeed, higher percentages of bone marrow NK cells from healthy donors produced IFN-γ upon PMA/ionomycin stimulation than did than bone marrow NK cells from AML patients. In addition, the impairment of NK cells from AML patients at the site of disease onset was consistent with the capacity of AML blasts to affect NK cells in vitro.<sup>31,71,72</sup> The exact mechanism by which tumor cells affect NK cells remains to be determined.

Bone marrow from AML cancer patients also exhibited down-regulation of mRNA expression of CD160, whose protein expression was also found to be downregulated. Moreover, among AML patients, higher CD160 levels correlated with better survival, suggesting that CD160 could be a marker of interest in AML. This glycoprotein activates the effector functions of NK cells.<sup>73</sup> However, the biological relevance of CD160 signaling at the surface of NK cells in the context of AML awaits further investigation.

Overall, our study identified an NK0 CD56<sup>bright</sup>-like subset as the precursor of the NK1-CD56<sup>dim</sup> and NK2-CD56<sup>bright</sup> subsets. In addition, our findings regarding AML revealed that NK cells undergo donor-specific cancer-associated editing and provided further evidence of profound NK cell impairment in the bone marrow at diagnosis.

## EXPERIMENTAL MODEL DETAILS

### Samples

Frozen healthy bone marrow samples were recovered from filters used for allogeneic hematopoietic graft preparation. The use of such cells for scientific research is approved by the French authorities in accordance with the regulations in force (authorization DC 2018-3143). Frozen bone marrow samples obtained from patients with AML at diagnosis were selected from an anonymized database registered with the "Commission Nationale de l'Informatique et des Libertés" (the French data protection agency; authorization no. FbP1089790#) and provided by the "Centre de Ressources Biologiques Cancer, Bordeaux Biothèques Santé" at Bordeaux University Hospital. Written informed consent was obtained from each patient for the use of his/her biological samples for research, in accordance with the Declaration of Helsinki. Detailed information about the donors is provided in Supplementary Table S2. This study was performed in accordance with French law (Art. L.1243-1 and Art. L.1245-2 of the French Public Health Code).

## DETAILED METHODS

### Cell preparation and NK cell enrichment

Samples were maintained in liquid nitrogen, transported on dry ice and stored either in liquid nitrogen or at -80 °C if processing was planned shortly after arrival. On the day of processing, samples were thawed and put in RPMI + 20% FCS medium containing DNase I (Roche). Cells were then washed with DPBS + 5% FCS + 2 mM EDTA. Cell viability was assessed with Trypan Blue dye. Cells were then transferred into a 5 mL FACS tube and washed with FACS buffer. NK cells were then enriched by magnetic labeling of the contaminating cells with a Miltenyi

Biotec NK cell isolation kit, leading to specific enrichment of unlabeled NK cells.

### NK cell sorting by flow cytometry

After enrichment, cells were washed and incubated at a dilution of 1/50 with mouse serum (Sigma-Aldrich) for 10 min at 4 °C. Cells were stained with anti-CD3 (clone UCHT1, FITC, Beckman Coulter), anti-CD14 (clone RMO52, FITC, Beckman Coulter), anti-CD19 (clone J3-119, FITC, Beckman Coulter) anti-CD45 (clone HI30, APC-Cy7, Biolegend) and anti-CD56 (clone N901(NKH-1), PE, Beckman Coulter) antibodies for 30 min at 4 °C. We added an anti-CD34 (clone 581, FITC, BD Biosciences) antibody to the lineage to sort AML samples. Cells were then washed twice with DPBS and incubated in DPBS with a dead cell marker for 10 min at 4 °C. Cells were washed again and sorted in an Influx or Melody cell sorter from BD (Becton Dickinson, San Diego, USA).

### Single-cell RNA sequencing

After sorting, cells were washed in DPBS + 0.04% BSA, as recommended by the 10x Genomics sample preparation protocol, and kept on ice until counting was performed. We used a 10x Genomics Chromium single-cell 3' v2 kit and protocol to prepare the libraries. HalioDx (Marseille, France) performed RNA sequencing on the NextSeq 500 platform with a sequencing depth of at least 50,000 reads per cell.

### Flow cytometric assessment of the cell surface phenotype

Frozen bone marrow from healthy donors and a set of AML patients different from that used for scRNAseq analysis was thawed, treated with DNAse I (Roche), and stained the following day with antibodies against CD122 (clone Mik-β3, BV786, BD Biosciences), CD127 (clone HIL-7R-M21, Pe-Cy7, BD Biosciences), CD14 (clone M5E2, BUV737, BD Biosciences), CD16 (clone 3G8, BUV737, BD Biosciences), CD160 (clone BY55, Alexa Fluor® 488, BD Biosciences), CD19 (clone SJ25C1, BUV737, BD Biosciences), CD27 (clone M-T271, PerCP-CyTM™5.5, BD Biosciences), CD3 (clone UCHT1, BUV496, BD Biosciences), CD45 (clone 2D1, APC-H7, BD Biosciences), CD52 (clone 4C8, Alexa Fluor® 647, BD Biosciences), CD56 (clone B159, Pe-Cy7, BD Biosciences), CD69 (clone FN50, BV605, BD Biosciences), perforin (clone δG9, PE-CF594, BD Biosciences), KIR2DL2/L3/S2 (clone GL183, PE-Cy5, Beckman Coulter), CD45 (clone HI30, APC-Cy7, Biolegend), and KIR2DL1/S1 (clone 11PB6, PE-Vio770, Miltenyi Biotec), and the dead cell marker (BV510, Invitrogen). Samples were analyzed in a Canto II cytometer from BD (Becton Dickinson, San Diego, USA) or on a Fortessa cytometer with UV configuration from BD (Becton Dickinson, San Diego, USA). NK cells were defined as CD3⁻CD56⁺ or CD3⁻CD14⁻CD19⁻CD56⁺ cells. Statistical significance was assessed by Friedmann's analysis with Dunn's post hoc test on paired subset data or with Dunn's test, and *p* values were adjusted by the Benjamini-Hochberg method. \* *p* value < 0.05, \*\* *p* value < 0.01, \*\*\* *p* value < 0.001.

### Flow cytometry to assess NK cell functionality

Frozen bone marrow from healthy donors and a set of AML patients different from that used for scRNAseq analysis was thawed, treated with recombinant DNAse I (150 U, Roche) and then dispensed the following day into U-bottom 96-well plates at a density of 5 × 10⁵ PBMCs/well. Cells were incubated for 4 h at 37 °C in the presence of Golgi Stop (1/1500; Becton Dickinson), anti-CD107 antibodies (anti-CD107a clone H4A3 and anti-CD107b clone H4B4, FITC, Becton Dickinson), PMA (2 ng/mL) and ionomycin (500 ng/mL). Cells were then washed in PBS supplemented with 2% FCS and 1 mM EDTA and stained by incubation for 30 min at 4 °C in buffer containing an anti-CD3 antibody (clone SK7, PerCP-Cy5.5, Becton Dickinson), an anti-CD56 antibody (Clone N901, APC, Beckman Coulter), and 2% normal mouse serum (Sigma Aldrich). Cells were fixed and permeabilized and were then
incubated with 2% paraformaldehyde (PermWash; Becton Dickinson) for 10 min in the dark at room temperature. They were then stained by incubation for 30 min at 4 °C with an anti-IFNγ antibody (clone 4S-B3, PE, Becton Dickinson). Samples were analyzed in a Canto II cytometer from BD (Becton Dickinson, San Diego, USA). NK cells were defined as CD3<sup>-</sup> CD56<sup>+</sup> cells. Statistical significance was assessed by Dunn's test, and *p* values were adjusted by the Benjamini-Hochberg method. \* *p* value < 0.05, \*\* *p* value < 0.01, \*\*\* *p* value < 0.001, \*\*\*\* *p* value < 0.0001.

# QUANTIFICATION AND STATISTICAL ANALYSIS

## Preprocessing of samples

Raw FASTQ files were processed with CellRanger software (v3.0.0 or 3.0.1), which performs alignment, filtering, barcode counting and UMI counting. CellRanger software was used to align reads with the GRch38 genome. Low-quality cells were excluded in an initial quality-control step, which removed genes expressed in fewer than three cells and cells expressing fewer than 200 genes. Cells with less than 10% ribosomal genes in samples from healthy controls and less than 15% ribosomal genes in samples from patients with AML were excluded. We also excluded cells with a UMI count of more than three median absolute deviations from the median value. In total, 23,850 healthy and 15,046 AML cells were retained for further analysis. Library size normalization was performed on the gene expression values for each cell barcode after UMI collapsing, with scaling by the total number of transcripts and multiplication by 10,000. The data were then log-transformed before further downstream analysis with Seurat.<sup>74</sup>

## Sample analysis

We first considered each donor separately. We selected genes with a high variance using the FindVariableGenes function with default parameters. We then reduced the dimensionality of our data by PCA and selected the number of variable PCs retained by random sampling. Cells were clustered with Seurat's FindClusters function with the Louvain algorithm. We tested a range of cluster resolution parameters and selected the parameter validated by the machine learning approach of Valentine Svensson used in one of our previous studies ([http://www.nxn.se/valent/2018/3/5/actionable-scrna-seq-clusters](http://www.nxn.se/valent/2018/3/5/actionable-scrna-seq-clusters)).<sup>12</sup> We removed cells assigned to clusters accounting for <2% of the cells to avoid doublets, and we also removed a small cluster of cells from two donors that were clustered away from the other clusters and probably resulted from ILC contamination (87 cells from one donor and 32 cells from the other donor). These cells had a high ILC score (as shown below). We also removed a donor-specific subset of 37 cells not found in seven other donors. We checked for the absence of any remaining contamination with the SingleR package.<sup>75</sup> For visualization, we applied RunTsne and RunUmap to the cell loadings of the previously selected PCs with default parameters to view the cells in two dimensions. We identified cluster markers with FindAll-Markers, with the parameter only.pos set to TRUE, to obtain only upregulated genes as markers of a cluster relative to all other cells. Marker genes were defined as genes with an adjusted *p* value of < 0.05 by the nonparametric Wilcoxon rank-sum test.

## Pooled sample analysis

We pooled samples with the CellRanger aggregate function and subjected them to normalization to equalize the read depth between libraries. We performed regression on the runs to overcome a batch effect that divided human healthy bone marrow NK cell samples into two categories on the basis of processing. We created three different pools of samples: a pool of healthy donor bone marrow NK cells, a pool of bone marrow NK cells from donors with AML and a pool combining these other pools. We then analyzed the pooled samples as described above.

## Assessment of NK cell purity

We checked for ILC contamination in these pools by scoring our cells with the AddModuleScore function (Seurat) and the gene signatures described by<sup>34</sup> for ILCs and<sup>10</sup> for NK cells. The ILC signature was obtained from<sup>34</sup> by selecting the genes from the NK vs ILC1, 2, and 3 comparison that were upregulated in ILCs and for which the *p* value was <0.05. The NK cell gene signature was obtained from<sup>10</sup> by merging the list of CD56<sup>bright</sup> CD16<sup>-</sup> vs CD56<sup>dim</sup> CD16<sup>+</sup> upregulated gene with the list of CD56<sup>bright</sup> CD16<sup>-</sup> vs CD56<sup>dim</sup> CD16<sup>+</sup> downregulated genes, both of which were sorted with an FC of >2. This step resulted in the additional removal of 10 cells from the bone marrow of healthy donors and 41 from the bone marrow of donors with AML.

## Unsupervised hierarchical clustering

We calculated the mean expression levels of genes with differential expression across clusters and performed hierarchical clustering with these values. For all samples, the Euclidean distance was used for samples, genes and clusters.

## Principal component analysis

PC clustering analysis was performed with the Ade4 package on the mean expression level of all genes among clusters. The gene loadings for PC1 and PC2 corresponded to the genes making the largest contributions, which accounted for 20% of the total information for each PC.

## Heatmap

Heatmaps were generated from the scaled expression (log-normalized UMI counts) values for the discriminating gene sets defining each subset, with an adjusted *p* value of <0.05 by the nonparametric Wilcoxon rank-sum test. The color scale is based on the z-score distribution.

## Gene annotations

Cell membrane protein, secreted protein and transcription factor annotations were retrieved from public databases (UniProt, MGI, and NCBI for mice and UniProt, GeneCards, and The Human Protein Atlas for humans). Genes encoding transcription factors were defined as such only if "transcription factor activity" was found among the GO annotations for the gene to prevent the confusion of cofactors and coregulators with transcription factors.

## GO enrichment analysis

We performed GO enrichment analysis with BiomaRt<sup>76</sup> and the GOStats package.<sup>77</sup> Enrichment scores (*p* values) for selected GO annotations were calculated by a hypergeometric statistical test with a significance threshold of 0.05. The data were plotted as the -log10 *p* values after Benjamini-Hochberg correction. The significance threshold was set at -log10(0.05).

## Comparison with known NK cell subsets

Module scores for CD56<sup>bright</sup>, CD56<sup>dim</sup>, hNK1, hNK2, hNK_Sp1, hNK_Sp2, and hNK_Sp3 gene expression profiles, as defined by<sup>10</sup> and<sup>12</sup> were determined with the AddModuleScore function of Seurat for each of our NK cells at the single-cell level. Briefly, the mean expression level of each gene in the defined expression profiles was calculated for each cell, and the aggregated expression for the control gene sets was then subtracted. All analyzed genes were binned on the basis of the mean expression level, and the control genes were randomly selected from each bin. Violin plots were used to assess the distribution of module scores for each NK cell grouped by subset. Statistical significance was assessed by the Wilcoxon rank sum test with continuity correction or by the Kruskal-Wallis test with Dunn's multiple comparisons test, with *p* value adjustment by the
Benjamini–Hochberg method. *p value < 0.05, **p value < 0.01, ***p value < 0.001, ****p value < 0.0001.

## Pseudotime analysis

For bone marrow cells, we performed pseudotime analysis with Monocle 3 (v0.1.3), and the cells were filtered with Seurat, as described above, for each sample individually. We imported Seurat clustering and NK/NKP scoring metadata into the analysis. The starting point of the trajectory was defined as the endpoint of the branch with the highest NKP signature score.

For splenic cells, we used our previously published dataset¹² and again performed pseudotime analysis with Monocle 3 (v0.1.3). The starting point of the trajectory was defined as the endpoint of the branch with the lowest mature NK cell score, and cluster assignments were retrieved from our previous analysis.

## Definition of the NK cell vs. NKP gene signature

Affymetrix CEL files were analyzed with the limma package (v3.34.9) with RMA normalization. Genes with an adjusted *p* value of <0.01 were considered to be differentially expressed.

Definition of the canonical vs. adaptive NK cell gene signature RNA-seq analysis was performed with the DESeq2 package (v1.26.0). Genes with a log2 fold change in expression of >1 and an adjusted *p* value of <0.05 were considered to be differentially expressed.

## Analysis of TCGA data

Analysis of TCGA data was performed with Xena, developed by the University of California Santa Cruz.⁷⁸ The TCGA dataset for acute myeloid leukemia was selected, and the overall survival of patients stratified by the intensity of CD160 expression was analyzed via the Kaplan–Meier method. The optimal cutoff value for patient stratification is given, and the *p* value indicated in the plot was calculated with the log-rank test.

## DATA AND SOFTWARE AVAILABILITY

All sequencing data will be deposited in the NCBI GEO repository and will be accessible in GEO. The schematic representation in Fig. S4 was designed with BioRender.

## ACKNOWLEDGEMENTS

The E.V. laboratory at CIML and Assistance-Publique des Hôpitaux de Marseille is supported by funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (TILC, grant agreement no. 694502 and Minfla-TILC, grant agreement no. 875102 - Minfla-Tilc); the Agence Nationale de la Recherche, including the PIONEER Project (ANR-17-RHUS-0007); MSDAvenir; Innate Pharma; and institutional grants to the CIML (INSERM, CNRS, and Aix-Marseille University) and to Marseille Immunopole. This work was partly supported by a grant from the Fondation de France (DDS Cancer 2017), and P.-Y. D. received a Fondation de France research fellowship. We thank the Centre de Ressources Biologiques Cancer, Bordeaux Biothèques Santé (BB-0033-00036) at Bordeaux University Hospital for providing biological material. We thank the Bioinformatics and Genomics Core Facilities of Centre d'Immunologie de Marseille-Luminy for helpful discussions and comments.

## AUTHOR CONTRIBUTIONS

E.V. conceived the project. A.C., P.-Y.D., E.N.M. and E.V. designed the experiments. P.-Y.D. and F.V. provided access to human samples. A.C., P.-Y.D., C.P., P.M., L.G., and A.V. performed the experiments. A.C. performed the 10x Genomics protocol to generate the libraries. B.E. designed the bioinformatic algorithms for data analysis. B.E. analyzed the scRNAseq data with the help of A.C. A.C., P.-Y.D., E.N.M., and E.V. analyzed the data. A.C., P.-Y.D., E.N.M., and E.V. wrote the manuscript.

## ADDITIONAL INFORMATION

The online version of this article (https://doi.org/10.1038/s41423-020-00574-8) contains supplementary material.

**Competing interests:** E.V. is a cofounder and employee of Innate Pharma. The other authors declare no competing interests.

## REFERENCES

1. Vivier, E. et al. Innate or adaptive immunity? The example of natural killer cells. *Science* **331**, 44–49 (2011).

2. Moretta, A. et al. Activating receptors and coreceptors involved in human natural killer cell-mediated cytolysis. *Annu Rev. Immunol.* **19**, 197–223 (2001).

3. Narni-Mancinelli, E., Ugolini, S. & Vivier, E. Tuning the threshold of natural killer cell responses. *Curr. Opin. Immunol.* **25**, 53–58 (2013).

4. Guia, S., Fenis, A., Vivier, E. & Narni-Mancinelli, E. Activating and inhibitory receptors expressed on innate lymphoid cells. *Semin Immunopathol.* **40**, 331–341 (2018).

5. Spits, H. et al. Innate lymphoid cells-a proposal for uniform nomenclature. *Nat. Rev. Immunol.* **13**, 145–149 (2013).

6. Walzer, T., Jaeger, S., Chaix, J. & Vivier, E. Natural killer cells: from CD3(-)Nkp46(+) to post-genomics meta-analyses. *Curr. Opin. Immunol.* **19**, 365–372 (2007b).

7. Yu, J. et al. Nkp46 identifies an NKT cell subset susceptible to leukemic transformation in mouse and human. *J. Clin. Invest* **121**, 1456–1470 (2011).

8. Bjorkstrom, N. K., Ljunggren, H. G. & Michaelsson, J. Emerging insights into natural killer cells in human peripheral tissues. *Nat. Rev. Immunol.* **16**, 310–320 (2016).

9. Cooper, M. A., Fehniger, T. A. & Caligiuri, M. A. The biology of human natural killer-cell subsets. *Trends Immunol.* **22**, 633–640 (2001).

10. Hanna, J. et al. Novel insights on human NK cells’ immunological modalities revealed by gene expression profiling. *J. Immunol.* **173**, 6547–6563 (2004).

11. Yu, J., Freud, A. G. & Caligiuri, M. A. Location and cellular stages of natural killer cell development. *Trends Immunol.* **34**, 573–582 (2013).

12. Crinier, A. et al. High-dimensional single-cell analysis identifies organ-specific signatures and conserved NK cell subsets in humans and mice. *Immunity* **49**, 971–986 e975 (2018).

13. Dulphy, N. et al. An unusual CD56(bright) CD16(low) NK cell subset dominates the early posttransplant period following HLA-matched hematopoietic stem cell transplantation. *J. Immunol.* **181**, 2227–2237 (2008).

14. Chan, A. et al. CD56(bright) human NK cells differentiate into CD56dim cells: role of contact with peripheral fibroblasts. *J. Immunol.* **179**, 89–94 (2007).

15. Lopez-Verges, S. et al. Expansion of a unique CD57(+)-NKG2Chi natural killer cell subset during acute human cytomegalovirus infection. *Proc. Natl Acad. Sci. USA* **108**, 14725–14732 (2011).

16. Cooper, M. A. et al. Cytokine-induced memory-like natural killer cells. *Proc. Natl Acad. Sci. USA* **106**, 1915–1919 (2009).

17. Bjorkstrom, N. K. et al. Rapid expansion and long-term persistence of elevated NK cell numbers in humans infected with hantavirus. *J. Exp. Med.* **208**, 13–21 (2011).

18. Abdul-Careem, M. F. et al. Critical role of natural killer cells in lung immunopathology during influenza infection in mice. *J. Infect. Dis.* **206**, 167–177 (2012).

19. Dou, Y. et al. Influenza vaccine induces intracellular immune memory of human NK cells. *PLoS One* **10**, e0121258 (2015).

20. Cichocki, F. et al. CD56dimCD57+NKG2C+ NK cell expansion is associated with reduced leukemia relapse after reduced intensity HCT. *Leukemia* **30**, 456–463 (2016).

21. Lee, J. et al. Epigenetic modification and antibody-dependent expansion of memory-like NK cells in human cytomegalovirus-infected individuals. *Immunity* **42**, 431–442 (2015).

22. Tesi, B., Schlums, H., Cichocki, F. & Bryceson, Y. T. Epigenetic regulation of adaptive NK cell diversification. *Trends Immunol.* **37**, 451–461 (2016).

23. Short, N. J. et al. Advances in the treatment of acute myeloid leukemia: new drugs and new challenges. *Cancer Discov.* **10**, 506–525 (2020).

24. Costello, R. T. et al. Defective expression and function of natural killer cell-triggering receptors in patients with acute myeloid leukemia. *Blood* **99**, 3661–3667 (2002).

25. Fauriat, C. et al. Deficient expression of NCR in NK cells from acute myeloid leukemia: evolution during leukemia treatment and impact of leukemia cells in NCRdull phenotype induction. *Blood* **109**, 323–330 (2007).

26. Chretien, A. S. et al. NKp46 expression on NK cells as a prognostic and predictive biomarker for response to allo-SCT in patients with AML. *Oncoimmunology* **6**, e1307491 (2017a).

27. Chretien, A. S. et al. NKp30 expression is a prognostic immune biomarker for stratification of patients with intermediate-risk acute myeloid leukemia. *Onco-target* **8**, 49548–49563 (2017c).
28. Chretien, A. S. et al. Natural killer defective maturation is associated with adverse clinical outcome in patients with acute myeloid leukemia. *Front Immunol.* **8**, 573 (2017b).

29. Pende, D. et al. Analysis of the receptor-ligand interactions in the natural killer-mediated lysis of freshly isolated myeloid or lymphoblastic leukemias: evidence for the involvement of the Poliovirus receptor (CD155) and Nectin-2 (CD112). *Blood* **105**, 2066–2073 (2005).

30. Sanchez-Correa, B. et al. Decreased expression of DNAM-1 on NK cells from acute myeloid leukemia patients. *Immunol. Cell Biol.* **90**, 109–115 (2012).

31. Stringaris, K. et al. Leukemia-induced phenotypic and functional defects in natural killer cells predict failure to achieve remission in acute myeloid leukemia. *Haematologica* **99**, 836–847 (2014).

32. Chiossone, L., Dumas, P. Y., Vienne, M. & Vivier, E. Natural killer cells and other innate lymphoid cells in cancer. *Nat. Rev. Immunol.* **18**, 671–688 (2018).

33. Zheng, G. X. et al. Massively parallel digital transcriptional profiling of single cells. *Nat. Commun.* **8**, 14049 (2017).

34. Bjorklund, A. K. et al. The heterogeneity of human CD127(+) innate lymphoid cells revealed by single-cell RNA sequencing. *Nat. Immunol.* **17**, 451–460 (2016).

35. Huntington, N. D. et al. IL-15 trans-presentation promotes human NK cell development and differentiation in vivo. *J. Exp. Med* **206**, 25–34 (2009).

36. van Helden, M. J. et al. Terminal NK cell maturation is controlled by concerted actions of T-bet and Zeb2 and is essential for melanoma rejection. *J. Exp. Med* **212**, 2015–2025 (2015).

37. Costa-Garcia, M. et al. Human cytomegalovirus antigen presentation by HLA-DR+ NKG2C+ adaptive NK cells specifically activates polyfunctional effector memory CD4+ T lymphocytes. *Front Immunol.* **10**, 687 (2019).

38. Liu, L. L. et al. Critical role of CD2 co-stimulation in adaptive natural killer cell responses revealed in NKG2C-deficient humans. *Cell Rep.* **15**, 1088–1099 (2016).

39. Liu, L. L. et al. Harnessing adaptive natural killer cells in cancer immunotherapy. *Mol. Oncol.* **9**, 1904–1917 (2015).

40. Schlums, H. et al. Cytomegalovirus infection drives adaptive epigenetic diversification of NK cells with altered signaling and effector function. *Immunity* **42**, 443–456 (2015).

41. Cichocki, F. et al. ARID5B regulates metabolic programming in human adaptive NK cells. *J. Exp. Med.* **215**, 2379–2395 (2018).

42. Beziat, V. et al. NK cell responses to cytomegalovirus infection lead to stable imprints in the human KIR repertoire and involve activating KIRs. *Blood* **121**, 2678–2688 (2013).

43. Stumpo, D. J. et al. Targeted disruption of Zfp36L2, encoding a CCCH tandem zinc finger RNA-binding protein, results in defective hematopoiesis. *Blood* **114**, 2401–2410 (2009).

44. Zhang, L. et al. ZFP36L2 is required for self-renewal of early burst-forming unit erythroid progenitors. *Nature* **499**, 92–96 (2013).

45. Walzer, T. et al. Natural killer cell trafficking in vivo requires a dedicated sphingosine 1-phosphate receptor. *Nat. Immunol.* **8**, 1337–1344 (2007a).

46. Chea, S. et al. CXCR6 expression is important for retention and circulation of ILC precursors. *Mediators Inflamm.* **2015**, 368427 (2015).

47. Burgess, M. et al. CD62L as a therapeutic target in chronic lymphocytic leukemia. *Clin. Cancer Res.* **19**, 5675–5685 (2013).

48. Avigdor, A. et al. CD44 and hyaluronic acid cooperate with SDF-1 in the trafficking of human CD34+ stem/progenitor cells to bone marrow. *Blood* **103**, 2981–2989 (2004).

49. Beziat, V., Descours, B., Parizot, C., Debre, P. & Vieillard, V. NK cell terminal differentiation: correlated stepwise decrease of NKG2A and acquisition of KIRs. *PLoS One* **5**, e11966 (2010).

50. Trapnell, C. et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. *Nat. Biotechnol.* **32**, 381–386 (2014).

51. Renoux, V. M. et al. Identification of a human natural killer cell lineage-restricted progenitor in fetal and adult tissues. *Immunity* **43**, 394–407 (2015).

52. Leavy, O. Natural killer cells: RAG keeps natural killers fit. *Nat. Rev. Immunol.* **14**, 716–717 (2014).

53. Rabacal, W. et al. Transcription factor KLF2 regulates homeostatic NK cell proliferation and survival. *Proc. Natl Acad. Sci. USA* **113**, 5370–5375 (2016).

54. Celis-Gutierrez, J. et al. Dok1 and Dok2 proteins regulate natural killer cell development and function. *EMBO J.* **33**, 1928–1940 (2014).

55. Huang, L., Zhu, P., Xia, P. & Fan, Z. WASH has a critical role in NK cell cytotoxicity through Lck-mediated phosphorylation. *Cell Death Dis.* **7**, e2301 (2016).

56. Jevremovic, D. et al. Cutting edge: a role for the adaptor protein LAT in human NK cell-mediated cytotoxicity. *J. Immunol.* **162**, 2453–2456 (1999).

57. Awasthi, A. et al. Rap1b facilitates NK cell functions via IQGAP1-mediated signalosomes. *J. Exp. Med* **207**, 1923–1938 (2010).

58. Mace, E. M. & Orange, J. S. Lytic immune synapse function requires filamentous actin deconstruction by Coronin 1A. *Proc. Natl Acad. Sci. USA* **111**, 6708–6713 (2014).

59. Sasanuma, H. et al. Dual function for the adaptor MIST in IFN-gamma production by NK and CD4+NKT cells regulated by the Src kinase Fgr. *Blood* **107**, 3647–3655 (2006).

60. Wilton, K. M., Overlee, B. L. & Billadeau, D. D. NKG2D-DAP10 signaling recruits EVL to the cytotoxic synapse to generate F-actin and promote NK cell cytotoxicity. *J. Cell Sci.* **133**, 1–12 (2019).

61. Witalisz-Siepracka, A. et al. Loss of JAK1 drives innate immune deficiency. *Front. Immunol.* **9**, 3108 (2018).

62. Crozat, K. et al. Impact of beta2 integrin deficiency on mouse natural killer cell development and function. *Blood* **117**, 2874–2882 (2011).

63. Viant, C. et al. SHP-1-mediated inhibitory signals promote responsiveness and anti-tumour functions of natural killer cells. *Nat. Commun.* **5**, 5108 (2014).

64. Yang, C. et al. Heterogeneity of human bone marrow and blood natural killer cells defined by single-cell transcriptome. *Nat. Commun.* **10**, 3931 (2019).

65. Lim, A. I. et al. Systemic human ILC precursors provide a substrate for tissue ILC differentiation. *Cell* **168**, 1086–1100 e1010 (2017).

66. Romagnani, C. et al. CD56brightCD16- killer Ig-like receptor- NK cells display longer telomeres and acquire features of CD56dim NK cells upon activation. *J. Immunol.* **178**, 4947–4955 (2007).

67. Freud, A. G., Yu, J. & Caligiuri, M. A. Human natural killer cell development in secondary lymphoid tissues. *Semin Immunol.* **26**, 132–137 (2014).

68. Lugthart, G. et al. Human lymphoid tissues harbor a distinct CD69+CXCR6+ NK cell population. *J. Immunol.* **197**, 78–84 (2016).

69. Melsen, J. E., Lugthart, G., Lankester, A. C. & Schilham, M. W. Human circulating and tissue-resident CD56(bright) natural killer cell populations. *Front Immunol.* **7**, 262 (2016).

70. Bruck, O. et al. Immune profiles in acute myeloid leukemia bone marrow associate with patient age, T-cell receptor clonality, and survival. *Blood Adv.* **4**, 274–286 (2020).

71. Khaznadar, Z. et al. Acute myeloid leukemia impairs natural killer cells through the formation of a deficient cytotoxic immunological synapse. *Eur. J. Immunol.* **44**, 3068–3080 (2014).

72. Scoville, S. D. et al. Human AML activates the aryl hydrocarbon receptor pathway to impair NK cell development and function. *Blood* **132**, 1792–1804 (2018).

73. Le Bouteiller, P. et al. CD160: A unique activating NK cell receptor. *Immunol. Lett.* **138**, 93–96 (2011).

74. Satija, R., Farrell, J. A., Gennert, D., Schier, A. F. & Regev, A. Spatial reconstruction of single-cell gene expression data. *Nat. Biotechnol.* **33**, 495–502 (2015).

75. Aran, D. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. *Nat. Immunol.* **20**, 163–172 (2019).

76. Durinck, S., Spellman, P. T., Birney, E. & Huber, W. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. *Nat. Protoc.* **4**, 1184–1191 (2009).

77. Falcon, S. & Gentleman, R. Using GOstats to test gene lists for GO term association. *Bioinformatics* **23**, 257–258 (2007).

78. Goldman, M. J. et al. A user guide for the online exploration and visualization of PCAWG data. *Nat. Commun.* **11**, 3400–3409 (2020).

<figure><img src="image_31.png" /></figure>

**Open Access** This article is licensed under a Creative Commons

Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020