High-Dimensional Single-Cell Analysis Identifies Organ-Specific Signatures and Conserved NK Cell Subsets in Humans and Mice.

Natural killer (NK) cells are innate lymphoid cells (ILCs) involved in antimicrobial and antitumoral responses. Several NK cell subsets have been reported in humans and mice, but their heterogeneity across organs and species remains poorly characterized. We assessed the diversity of human and mouse NK cells by single-cell RNA sequencing on thousands of individual cells isolated from spleen and blood. Unbiased transcriptional clustering revealed two distinct signatures differentiating between splenic and blood NK cells. This analysis at single-cell resolution identified three subpopulations in mouse spleen and four in human spleen, and two subsets each in mouse and human blood. A comparison of transcriptomic profiles within and between species highlighted the similarity of the two major subsets, NK1 and NK2, across organs and species. This unbiased approach provides insight into the biology of NK cells and establishes a rationale for the translation of mouse studies to human physiology and disease.

Introduction

Natural killer (NK) cells are cytotoxic innate lymphoid cells (ILCs) that produce proinflammatory cytokines, chemokines, and growth factors, such as interferon (IFN)-γ, CCL3, and GM-CSF (Vivier et al., 2011, Vivier et al., 2018). NK cell effector activities are tightly controlled by inhibitory and activating signals originating from cell surface receptors (Bryceson et al., 2006). Several subsets of NK cells with different effector functions or maturation states have been described in both humans and mice, based on the relative abundances of specific cell-surface proteins. In humans, NK cells are subdivided into different populations based on the surface markers CD16 and CD56. The two main blood populations in healthy individuals are the cytotoxic CD56dimCD16+ NK cells and the CD56brightCD16− NK cells, which are less cytotoxic but produce larger amounts of cytokines upon exposure to environmental cues, such as interleukin (IL)-12 and IL-18 (Freud et al., 2017, Spits et al., 2016). CD56brightCD16− NK cells are thought to be at an intermediate stage on the way to becoming CD56dimCD16+ cells (Romagnani et al., 2007). Several NK maturation steps have been described in mice, based on the relative amounts of CD27 and CD11b. CD27+CD11b– cells are immature NK cells. They mature into double-positive CD27+CD11b+ cells and, finally, into CD27−CD11b+ NK cells. This developmental program is associated with the acquisition of NK cell effector functions (Chiossone et al., 2009, Hayakawa and Smyth, 2006, Kim et al., 2002).

Human NK cells have been studied mostly in the blood compartment, at least partly for practical and ethical reasons. By contrast, most studies in mice have been performed on splenic NK cells. Despite these experimental differences, it is becoming clear that NK cells are highly diverse in both species (Shi et al., 2011, Wilk and Blish, 2018). Nevertheless, key unanswered questions remain. Are the NK cells in the bloodstream functionally similar to those in tissues? How relevant are comparisons of NK cells from the two species? Can information gleaned from studies in mice be translated into clinically useful knowledge for treating human diseases? Efforts have been made to profile NK cells. However, such studies were limited by the technical approaches used, such as bulk transcriptomic analysis (Bezman et al., 2012, Robinette et al., 2015), which is dependent on data averaging, or mass cytometry, which is biased by the choice and number of molecules included in the assay (Horowitz et al., 2013, Simoni et al., 2017). Tonsil NK cells have been recently profiled by single-cell RNA-seq analysis (scRNA-seq), but it remains unclear whether these cells are representative of blood NK cells (Björklund et al., 2016). There is, therefore, a need to define NK cell heterogeneity at the single-cell level in an unsupervised pan-genomic analysis across organs and species.

We used high-throughput scRNA-seq to decipher the similarities and differences between human and mouse and between blood and splenic NK cells. scRNA-seq is a state-of-the-art method for characterizing the gene expression profiles of thousands of cells simultaneously at the single-cell level (Zheng et al., 2017). It is particularly useful for deciphering cellular complexity and diversity. We used unsupervised hierarchical clustering to identify organ-specific NK cell gene signatures in both humans and mice. We also distinguished multiple NK cell subsets in the blood and spleen of both species and predicted functional differences between the different subsets of NK cells in humans and mice. Comparisons between organs within or between species revealed similarities in specific NK cell subpopulations. These data will also serve as a valuable resource for translating the results of mouse NK cell studies into information improving our understanding of the role of these cells in human physiology and diseases.

Results

scRNA-Seq Profiling of Splenic and Blood Human and Mouse NK Cells Identifies a Species-Specific NK Cell Gene Signature

We investigated the NK cell profiles in the paired blood and spleen of brain-dead human donors and mice, with a Drop-Seq-based scRNA-seq technology (10× Genomics, Figures S1A and S1B and STAR Methods) (Zheng et al., 2017). We analyzed ∼4,000 NK cells from each of the mouse blood and spleen and from human spleen samples, and ∼3,000 human blood NK cells. scRNA-seq detected a mean of 914 genes per cell in the human dataset and 919 in the mouse dataset, resulting in a total of 16,806 human genes and 15,144 mouse genes detected in all cells.

We analyzed the specific NK cell signature of the cells recovered. NKp46 is a cell surface receptor detected on subsets of ILC1s, ILC3s, and T cells, but overall recognized as a marker of the NK cell lineage in humans and mice (Walzer et al., 2007). We generated with BioGPS a list of genes with an expression pattern strongly correlated (correlation coefficient ≥ 0.9) with that of the gene encoding NKp46 (NCR1 or Ncr1). For each species, we cross-checked this list of genes against the known NK cell-specific genes from available datasets (Bezman et al., 2012, Björklund et al., 2016). The intersect between the resulting list for each species and our own scRNA-seq dataset was then used to define NK cell-specific transcriptomic signatures. This approach identified 13 genes as defining a robust consensus NK cell signature in mice and humans (Figure S1C). In mice, the 13 NK cell-defining genes encoded activating receptors (Klrb1c [NK1.1], Klrk1 [NKG2D], and Ncr1 [NKp46]), inhibitory receptors (Klrg1 [KLRG1] and Klrb1b), the transcription factor eomesodermin (Eomes), the cytolytic proteins perforin (Prf1) and granzyme A (Gzma), the IL-18 coreceptor (Il18rap), and a sphingosine 1-phosphate receptor promoting NK cell egress from lymph nodes and bone marrow (S1pr5 [S1P5]) (Figure S1C). The 13 NK cell-defining genes in humans encoded activating receptors (CD160, KLRC3 [NKG2E], CD244 [2B4], KLRF1 [NKp80]), an inhibitory receptor (KLRC1 [NKG2A]), the antimicrobial protein granulysin (GNLY), the cytolytic protein perforin (PRF1), an inflammatory cytokine (XCL2), and the IL-18 coreceptor (IL18RAP). PRF1 and IL18RAP were the only genes common to both mouse and human NK cell gene signatures (Figure S1C). These results show that even if no gene could be assigned as NK cell specific, the combination of 13 genes in humans and mice defines a robust NK cell transcriptomic signature.

Mouse NK Cells Have an Organ-Specific Transcriptomic Profile, Indicative of a More Active Phenotype in the Spleen than in Blood

Projection of cells onto two dimensions in a t-distributed stochastic neighbor-embedding (t-SNE) analysis provides a visual representation of cell clustering based on their transcriptomic profile. In a t-SNE analysis of 8,118 mouse cells, splenic and blood mouse NK cells segregated into two different organ-specific clusters. Within the spleen and blood, NK cells from the same group of mice did not cluster together, leading to the conclusion that there was no individual-specific phenotype and no overwhelming technical batch effect (Figure 1A). We found that 164 genes (101 spleen-specific and 63 blood-specific) displayed significant differences in expression between splenic and blood NK cells (Table S1). The expression profiles of these genes resulted in a heatmap that clearly separated splenic and blood NK cells (Figure 1B). We identified the ten genes with the lowest p values and those encoding secreted proteins, cell membrane markers, and transcription factors (Figure 1C).

Figure 1

Mouse NK Cells Have an Organ-Specific Transcriptomic Profile

(A) t-SNE plot of 8,118 NK cells from mouse blood (3,936) and spleen (4,182).

(B) Heatmap of the 164 total genes distinguishing between mouse blood (63) from splenic (101) NK cells tested with a Wilcoxon rank-sum test. Cells are plotted in columns, by organ source, and genes are shown in rows, ranked by adjusted p value < 0.05. Gene expression is color coded with a scale based on z-score distribution, from −2 (purple) to 2 (yellow). Squares identify the mouse organ-specific NK cell transcriptomic signature.

(C) Top ten expressed total genes and top ten expressed genes encoding secreted proteins, cell membrane markers, and transcription factors differentiating significantly between spleen and blood NK cells. Genes are ranked by p value.

(D) Selected Gene Ontology terms. Benjamini and Hochberg-corrected −log10 p values from hypergeometric tests. The black dotted line represents the significance threshold set at −log10(0.05).

The top ten genes expressed more strongly in mouse splenic than blood NK cells included several encoding transcription factors (Jun, Fos, Fosb, and Klf2), two chemokine-encoding genes (Ccl3 and Ccl4), and four genes involved in regulatory protein or transcription factor activity (Nfkbia, Ppp1r15a, Pim1, and Ier5) (Figure 1C). The genes expressed more strongly in splenic NK cells than in blood NK cells included genes encoding well-known NK cell effector proteins (Ifng, Prf1, Ccl3, Ccl4, and Xcl1) and the activating receptors (Ncr1, Klrk1, and Klrb1c), the common cytokine receptor polypeptide (Il2rg), and an activation marker (Cd69). NK cells from mouse blood had a very different gene expression profile. They differentially expressed Tgfb1 (encoding cytokine transforming grown factor β1) and Ctla2a (encoding a negative regulator of the inflammatory response in activated T cells), the gene encoding a subunit of the IFN-γ receptor (Ifngr1), and genes encoding proteins involved in cell signaling: Nrarp (encoding a protein involved in Notch signaling), Dusp2 (encoding a regulator of the ERK pathway), and Arhgap45 (encoding a Rho guanosine triphosphatase activating protein). These data suggested that splenic NK cells have a more activated phenotype than blood NK cells.

Gene ontology (GO) enrichment analysis indicated that mouse blood NK cells were specifically enriched in genes associated with the Notch signaling pathway (Figure 1D). By contrast, splenic NK cells displayed an enrichment in many biological process terms, such as response to stress, response to stimulus, defense response, signal transduction, and regulation of gene expression, consistent with the more strongly activated phenotype predicted from analysis of the top-ranking genes in the various categories and their association with NK cell activity (Figures 1C and 1D). Consistently, splenic NK cells reacted more strongly than their paired blood NK cell samples upon stimulation (Figure S2).

High-Throughput scRNA-Seq Identifies Three Subsets of Mouse Splenic NK Cells

To assess mouse NK cell heterogeneity within the spleen, we performed unsupervised hierarchical clustering on the 4,182 mouse splenic NK cells (data not shown). NK cells did not cluster on the basis of sample, but into three different subsets for each sample, which we named mNK_Sp1 to 3. A representative t-SNE plot of 1,439 cells from one sample shows these three subsets (Figure 2A). We identified 159 genes with significantly different expression patterns between mouse splenic NK cell subsets (Table S2), and a heatmap of these genes clearly separated these subsets on the basis of differences in their transcriptomic profiles (Figure 2B).

Figure 2

High-Throughput scRNA-Seq Identifies Three Mouse Splenic NK Subsets

(A) t-SNE plot of 1,439 mouse splenic NK cells representative of one sample.

(B) Heatmap of the 159 genes tested with a Wilcoxon rank sum test separating the 4,182 splenic NK cells into subsets. White lines separate the samples. Squares identify specific transcriptomic signatures of mouse splenic NK cell subsets.

(C) Left: PCA for the three mouse splenic NK cell subsets from each sample based on the mean expression of the genes with variable expression. Right panel: driving genes for each cell subset accounting for 40% of total information in each PC.

(D) Top ten genes significantly distinguishing between the three splenic mouse NK cells in a representative sample.

(E) Selected Gene Ontology terms.

(F) Left: Module scores of CD27−CD11b+ and CD27+CD11b− gene expression programs defined by Chiossone et al. (2009) for each spleen NK cell at the single-cell level. Middle and right: Violin plots represent the distribution of module score for CD27−CD11b+ (middle) and CD27+CD11b− (right) cells for each spleen NK cell, grouped by subsets. Kruskal-Wallis with Dunn’s multiple-comparison tests and Benjamini-Hochberg adjusted p values. Error bars indicated the mean (±SD). ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

We identified the sets of genes driving the segregation of the mouse spleen NK cell subsets by performing PC analysis (PCA) (Figure 2C). PC1 and PC2 encompassed 76% of the information and segregated the samples into three well-defined groups, with each set of samples represented in each of the groups (Figure 2C, left). PC2 separated mNK_Sp1 and mNK_Sp2 whereas both PC1 and PC2 contributed to the separation of mNK_Sp3 from the others. These results suggested that mNK_Sp3 differ from both mNK_Sp1 and mNK_Sp2, which appear to be more closely related to each other. Mapping the genes onto these groups identified four driving genes for mNK_Sp1: Cma1 (encoding a chymotryptic serine proteinase), Ly6c2 (encoding a cell membrane protein), Lgals1 (encoding a galectin), and Gzmb. mNK_Sp2 had three driving genes: Ctla2a, Xcl1, and Cd7 (encoding a cell surface receptor potentially involved in NK cell activation). mNK_Sp3 was defined by five genes: Ccl3, Ccl4, Pim1, Nfkbia, and Nr4a1 (encoding a member of the nuclear receptor family of transcription factors) (Figures 2C, right, and S3). Four of these five driving genes for mNK_Sp3 were among the top ten genes displaying the strongest preferential expression in splenic rather than blood NK cells: Ccl3, Ccl4, Nfkbia, and Pim1 (Figure 1C). As the mNK_Sp3 subset did not appear to be the largest of the spleen NK cell population (Figure 2A), this overlap indicates that mNK_Sp3 drives the splenic transcriptional profile.

We analyzed the top ten genes expressed in mNK_Sp1, mNK_Sp2, and mNK_Sp3, together with the top ten expressed genes encoding secreted proteins, cell membrane markers, and transcription factors (Figure 2D). Prf1 and Gzmb (encoding proteins with cytolytic activity) were differentially expressed in the mNK_Sp1 subset, which was also characterized by the expression of Itgam (Cd11b), Itgb2, Klrg1, and S1pr5, encoding cell membrane proteins. At the protein level, this population can be visualized by the expression of CD11b and KLRG1, a surface marker upregulated with maturation (Huntington et al., 2007) (Figures S4A and S4B). mNK_Sp2 cells were characterized by the expression of Xcl1, Ctla2a, Ltb (encoding effector proteins), and Cd7. mNK_Sp2 cells were also characterized by a higher Cd27 expression. This population was defined by flow cytometry as expressing CD27, CD28, and CD90 (Thy-1) (Figures S4A and S4B). An analysis of biological processes for mNK_Sp1 cells revealed specific enrichment in cytolysis and leukocyte migration, two processes involved in inflammatory responses. mNK_Sp2 cells were enriched in lymphocyte activation, cell adhesion, and the regulation of leukocyte migration (Figure 2E).

Consistent with the PC analysis (Figure 2C), the mNK_Sp3 subset displayed a pattern of gene expression regulation different from those of the other subsets. mNK_Sp3 cells appeared to be engaged in complex transcriptional regulation, as indicated by higher expression of several genes encoding proteins involved in the NF-κB pathway: Nfkbia, Nfkbid, Nfkbiz, and Nr4a1 (Figure 2D). mNK_Sp3 cells also expressed genes involved in cell survival and proliferation (Pim1) and growth and apoptosis (Gadd45b). The mNK_Sp3 cell subset was the subset displaying the highest expression of Ccl3 and Ccl4, encoding NK cell inflammatory chemokines. Unfortunately, no antibodies for cell surface markers were available yet for the discrimination of this population by flow cytometry. Consistent with its gene expression profile, the mNK_Sp3 subset was characterized by numerous biological processes, including the regulation of gene expression, response to stimulus, response to stress, and signal transduction (Figure 2E). As for the genes driving the mNK_Sp3 subset, the biological processes overlapped with those of the total mouse spleen NK cell population. This subset of NK cells thus makes a major contribution to the tissue-related differences in transcriptional profiles between the NK cells of the spleen and those of the blood (Figure 1D).

Using module scores calculated at the single-cell level, we compared the gene expression profiles of the three NK subsets from mouse spleen with those of CD27+CD11b− (immature) and CD27−CD11b+ (mature) NK cells, as previously defined (Chiossone et al., 2009). mNK_Sp1 cells, which had a higher Itgam (Cd11b) expression than the other two subsets, were characterized by high scores for the CD27−CD11b+ NK cell gene signature (Figure 2F, left). The genes strongly expressed in both CD27−CD11b+ cells and mNK_Sp1 cells were Cma1, Gzmb, Klrg1, Irf8, and Prf1. Conversely, mNK_Sp2 cells, which had a higher Cd27 expression than the other two subsets, were characterized by high scores for the CD27+CD11b− NK cell gene signature (Figure 2F, left). The genes in common between CD27+CD11b− cells and mNK_Sp2 were Ctl2a2, Emb, Xcl1, Ltb, and Cxcr3. This analysis also revealed a continuum between the mNK_Sp1 and mNK_Sp2 subsets with the greatest similarity to previously defined populations (Figure 2F, left). Violin plots of each module score between subsets revealed a significant enrichment in the gene signature of CD27−CD11b+ cells for mNK_Sp1 cells (Figure 2F, middle) and in the signature of CD27+CD11b− cells for mNK_Sp2 cells (Figure 2F, right). mNK_Sp3 cells were different from both these known subsets, consistent with the results of the PCA (Figures 2F and 2C), suggesting that this subset represents a previously unknown subset of mouse NK cells.

We also performed a pan-genomic, organ-matched module score analysis on our mouse splenic NK cells, comparing them with splenic ILC1 signatures in mice (Robinette et al., 2015; Figure S4E). Only a few cells might have corresponded to ILC1s and did not interfere with our analysis.

High-Throughput scRNA-Seq Identifies Two Splenic-like NK Cell Subsets in Mouse Blood

Unsupervised hierarchical clustering analysis of the variable genes detected in 3,936 paired mouse blood NK cells identified two populations in each sample, which we named mNK_Bl1 and mNK_Bl2, as shown by a representative t-SNE analysis on 1,223 cells from one sample (Figure 3A and data not shown). Different transcriptomic signatures can be inferred for these subsets from the 87 genes displaying significant differences in expression between the two subsets in all samples (Table S2 and Figure 3B).

Figure 3

High-Throughput scRNA-Seq Identifies Two Splenic-like NK Cell Subsets in Mouse Blood

(A) t-SNE plot of 1,223 mouse blood NK cells representative of one sample.

(B) Heatmap of the 87 genes tested with a Wilcoxon rank sum test separating the 3,936 mouse blood-specific NK cells into subsets. Squares identify the specific transcriptomic signatures of mouse blood NK cell subsets.

(C) Left: PCA for the two mouse blood NK cell subsets. Right: driving genes.

(D) Top ten genes significantly differentiating between the two blood mouse NK cell subsets in a representative sample.

(E) Selected Gene Ontology terms.

(F) Left: Module scores for CD27−CD11b+ and CD27+CD11b− gene expression programs. Middle and right: Violin plot showing the distribution of module scores for CD27−CD11b+ (middle) and CD27+CD11b− (right) cells, for each blood NK cell, grouped by subset. Wilcoxon rank sum test with continuity correction. Error bars indicated the mean (±SD). ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

(G) Left: PCA on the two blood and three splenic mouse NK cell subsets. Right: driving genes.

(H) Heatmap of the 428 genes tested with a Wilcoxon rank sum test separating the 8,118 mouse blood and splenic NK cells into subsets. Squares identify mouse-specific transcriptomic signatures of NK cell subsets.

PCA revealed that 92% of the information was encompassed by PC1 and PC2 (Figure 3C, left). mNK_Bl1 and mNK_Bl2 cells differed from each other on PC1, with Gzmb, Cma1, Irf8, and Ly6c2 identified as the driving genes for mNK_Bl1 cells and Ctla2a and Emb as the driving genes for mNK_Bl2 cells (Figure 3C, right, and S3). Gzmb, Cma1, Ly6c2, and Ctla2a were the genes characterized as driving the differences between splenic subsets (Figure 2C).

An analysis of the top ten differentially expressed genes and the ten top differentially expressed genes encoding secreted proteins, cell membrane proteins, and transcription factors showed that mNK_Bl1 cells expressed genes encoding the NK cell effector proteins Gzmb and Prf1, the cell membrane proteins Klrg1, Ly6c2, Klra9, Emp3, and Itgb2, and the transcriptional regulator Zeb2. All of these genes were also differentially expressed by mNK_Sp1 cells. Similarities between the mNK_Sp2 and mNK_Bl2 populations were evident in the top ten expressed genes, for which six genes were found to be common to both subsets (Ctla2a, Xcl1, Ltb, Cd27, Cd7, and Emb). Biological process analysis showed the mNK_Bl1 subset to be enriched in genes associated with cytolysis, defense response, response to stress, and response to virus, whereas the mNK_Bl2 subset was enriched in genes involved in lymphocyte activation, lymphocyte differentiation, and cell adhesion (Figure 3E). These results are similar to those obtained for the spleen subsets mNK_Sp1 and mNK_Sp2, respectively.

We calculated single-cell module scores for comparison of our two mouse blood NK subsets with the gene signatures of CD27+CD11b− (immature) and CD27−CD11b+ (mature) NK cells defined in a previous study (Chiossone et al., 2009). We found that the gene signature of mNK_Bl1 cells resembled that of CD27−CD11b+ NK cells and that the gene signature of mNK_Bl2 was similar to that of CD27+CD11b− NK cells (Figure 3F, left). As observed for the NK cells from mouse spleen, there was a continuum between the mNK_Bl1 and mNK_Bl2 populations most closely resembling the previously defined cell subsets (Figure 3F, left). CD11b and KLRG1 also identified CD27−CD11b+ cells in the blood, whereas CD27, CD28, and CD90 were associated with the CD27+CD11b− population (Figures S4C and S4D). Violin plots of each module score revealed a significant enrichment in the expression program of CD27−CD11b+ cells for mNK_Bl1 cells (Figure 3F, middle) and in that of CD27+CD11b− cells for the NK_Bl2 cells (Figure 3F, right). Consistent with this observation, we detected a higher expression of Cd27 in the mNK_Bl2 subset than in the mNK_Bl1 subset.

These analyses thus indicated similarities between mNK_Sp1 and mNK_Bl1 cells and between mNK_Sp2 and mNK_Bl2 cells. We tested this hypothesis by performing unsupervised hierarchical clustering analysis on all 8,118 mouse NK cells (data not shown). NK cells formed five subsets that could be grouped by organ source, consistent with an organ-specific phenotype of transcriptomic profiles. A PC analysis of all samples revealed that PC1 and PC2 together encompassed 80% of the information (Figure 3G, left). PC1 distinguished blood from splenic NK cells and mNK_Sp3 cells from the other subsets. The genes driving PC1 (Ccl3, Ccl4, Pim1, and Nfkbia) were identified as drivers of the splenic phenotype, as well as the difference between mNK_Sp3 cells and the other two splenic subsets (Figure 3G, right). PC2 distinguished the CD27−CD11b+ (mNK_Sp1 and mNK_Bl1) population from the CD27+CD11b− (mNK_Sp2 and mNK_Bl2) population in both spleen and blood. Ctla2a and Cd7 were drivers of the CD27+CD11b− phenotype, whereas Gmzb, Lgals1, Cma1, Ly6c2, and Irf8 were drivers of the CD27−CD11b+ phenotype (Figures 2C and 3C, right). A heatmap representation of the 8,118 mouse NK cells and the 428 differentially expressed genes between subsets and tissues clearly showed similarities between the transcriptomic signatures of mNK_Bl1 and mNK_Sp1 cells and between those of mNK_Bl2 and mNK_Sp2 cells (Figure 3H). The mNK_Sp3 subset was not represented in the blood. Overall, these results indicate that, despite the tissue-specific phenotype of mouse NK cells, the clustering of at least two of the splenic subsets was based on a transcriptomic signature similar to that underlying the clustering of blood NK cells.

Human NK Cells Have an Organ-Specific Transcriptomic Profile, Indicative of a More Active Phenotype in the Spleen than in Blood

In sharp contrast with the mouse data, scRNA-seq analysis of 7,403 human blood and spleen NK cells revealed that the human samples had a strong donor phenotype, probably reflecting differences in genetic origin and history of immunological challenges between individuals (Figure 4A). However, there was still a clear separation between splenic and blood NK cells, indicating that NK cells had an organ-specific phenotype in humans as well (Figure 4A).

Figure 4

Human NK Cells Have an Organ-Specific Transcriptomic Profile

(A) t-SNE plot of 7,403 human NK cells from blood (3,200) and spleen (4,203).

(B) Heatmap of the 294 total genes distinguishing human blood (106) from splenic (188) NK cells tested with a Wilcoxon rank-sum test.

(C) Top ten genes significantly differentiating spleen from blood NK cells.

(D) Selected Gene Ontology terms.

We identified 294 genes (188 spleen-specific and 106 blood-specific) displaying significant differential expression between human splenic and blood NK cells (Table S1). A heatmap representation of their expression separated splenic and blood NK cells (Figure 4B). We investigated the organ-specific differences in human NK cell populations, by comparing the top ten expressed genes in organ-specific NK cells and the top ten genes encoding secreted proteins, cell membrane proteins, and transcription factors (Figure 4C). The top ten expressed genes in human spleen NK cells included several genes encoding anti-inflammatory proteins, such as a subunit of the IKK complex (NFKBIA) and the products of the glucocorticoid-responsive genes TSC22D3 and ANXA1. Genes encoding anti-apoptotic proteins (BCL2A1 and MCL1) were also in the top ten genes expressed in human spleen NK cells. hNK_Sp NK cells were also characterized by stronger expression of genes encoding the secreted proteins CCL3, XCL1, IFNG, and GZMK, the proinflammatory cytokine IL-32 (an inducer of TNF-α and IL-6 by myeloid cells), and the cell-membrane proteins CD69, CD94 (KLRD1), ICAM1, and CD161 (KLRB1). Enrichment analysis for biological processes indicated that the gene expression profile of hNK_Sp was associated with a high NK functional activity, as represented by GO terms such as defense response, response to stimulus, response to cytokine, response to stress, and signal transduction (Figure 4D). By contrast, the most strongly expressed genes in human blood NK cells included several genes encoding “housekeeping” proteins, such as two members of the S100 family (S100A4 and S100A6, encoding proteins that increase in abundance with the differentiation of NK cells [Scheiter et al., 2013]), LDHA (encoding lactate dehydrogenase), and SLC7A5 (encoding a transporter for neutral amino acids). The gene expression profile of human blood NK cells was, thus, not enriched in any biological process (Figure 4D). Nevertheless, hNK_Bl cells strongly expressed GNLY, TGFB1, FGFBP2 (encoding a fibroblast growth factor-binding protein known to be secreted by cytotoxic lymphocytes), and HMGB1 (encoding a protein associated with CXCL12 and promoting the recruitment of inflammatory cells).

High-Throughput scRNA-Seq Identifies Four Subsets of Human Splenic NK Cells

Given the heterogeneity of the NK cells isolated from each of the donors, we performed unsupervised hierarchical clustering on individual data for each donor (data not shown). We identified four subsets of human splenic NK cells (hNK_Sp1 to 4) regardless of the donor analyzed. As shown in the representative t-SNE plot, hNK_Sp1 cells constituted the largest group, and hNK_Sp4 cells, the smallest (Figure 5A). As the hNK_Sp4 subset was quite small, we tested its identification by a machine-learning approach (STAR Methods). The 4th subset was reproducibly found to be different from subset 1 in all three donors (Figure S5). We then used the 270 genes identified as differentially expressed within the splenic subsets of one donor and analyzed their expression in the other samples (Table S3). These genes formed different subset-specific NK cells transcriptomic signatures in all human spleen samples (Figure 5B).

Figure 5

High-Throughput scRNA-Seq Identifies Four Human Splenic NK Subsets

(A) t-SNE plot of 1,321 human splenic NK cells representative of one sample.

(B) Heatmap of the 270 genes from one representative individual tested with a Wilcoxon rank sum test discriminating the human splenic NK cells into subsets. Squares identify specific transcriptomic signatures of the human splenic NK cell subsets.

(C) Left: PCA on the four human splenic NK cell subsets. Right: driving genes.

(D) Top ten genes significantly differentiating between the four spleen NK cell subsets in one representative sample.

(E) Selected Gene Ontology terms. hNK_Sp1 and hNK_Sp4 gene signatures are grouped.

(F) Left: Module scores of CD56bright and CD56dim gene expression programs defined by Hanna et al. (2004) for each of the blood NK cells at the single-cell level. Middle and right: Violin plots representing the distribution of module score for CD56dim (middle) CD56bright (right) for each blood NK cell, grouped by subset. Kruskal-Wallis with Dunn’s multiple-comparison test and Benjamini-Hochberg adjusted p values. Error bars indicated the mean (±SD). ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

We then analyzed the genes differentially expressed between the four subsets. PC1 and PC2 encompassed 70% of the information and distinguished between three groups of human spleen NK cells (Figure 5C, left). PC2 separated the hNK_Sp2 and hNK_Sp3 subsets. The hNK_Sp1 and hNK_Sp4 subsets were grouped together and separated from the hNK_Sp2 and hNK_Sp3 subsets by PC1 and PC2. Projecting the genes on the first two axes of the PCA identified LMNA (encoding nuclear lamin proteins displaying a higher expression in CD56bright cells) and CD52 as the genes driving the hNK_Sp3 subset (Figures 5C, right and S6). RGS1 and CD160 were the driving genes for the hNK_Sp2 subset. FCGR3A (encoding CD16, one of the two surface proteins used to identify human NK cell populations) was identified as one of the driving genes of the mixed cluster including hNK_Sp1 and hNK_Sp4 cells, together with other genes related to NK cell function, such as GZMB, GZMA, PRF1, CCL3, CCL4, GNLY, and FGFBP2 (Figures 5D and S6). hNK_Sp1 and hNK_Sp4 cells appeared to be related in the PCA analysis, but the tSNE and heatmap analyses revealed differences in gene expression that were sufficient to separate these cells into different clusters.

We compared the top ten differentially expressed genes in the four human splenic subsets and the top ten genes encoding secreted proteins, cell membrane proteins, and transcription factors (Figure 5D). The hNK_Sp1 and hNK_Sp4 subsets were characterized by a high expression of the effector protein-encoding genes FGFBP2, GZMB, GNLY, CCL4, and PRF1, and these cells expressed FCGR3A. The differential expression of these genes as specific markers was significant for hNK_Sp1 cells but not for the hNK_Sp4 subset. hNK_Sp2 cells expressed genes encoding the secreted molecules GZMK, XCL1, XCL2, LIF, and GM-CSF (CSF2) and the transcription factors IRF8, GATA3, FOSB, and NFKB1 more strongly than the other subsets. The hNK_Sp3 subset was characterized by a high of expression of LTB and GZMK (encoding secreted proteins), GATA3 and TCF7 (encoding transcription factors), and IL7R, CD52, and CD44 (encoding cell-surface proteins). Finally, hNK_Sp4 cells had specifically a higher expression of TXNIP (encoding a protein inhibiting the regulation of cellular redox signaling leading to oxidative stress), ARRDC3 (encoding a regulator of the cell-surface expression of the adrenergic receptor), and GIMAP7 (encoding a guanosine triphosphatase).

Biological process enrichment analysis for grouped hNK_Sp1 and hNK_Sp4 cells versus the hNK_Sp2 and hNK_Sp3 subsets identified multiple processes for each human spleen NK subset (Figure 5E). The hNK_Sp1 and hNK_Sp4 subsets had gene profiles displaying specific biological process enrichment in cytolysis; whereas hNK_Sp2 cells displayed specific enrichment in biological pathways associated with cytokine production, response to cytokine, immune system process, and defense response. Both hNK_Sp2 and hNK_Sp3 cells were enriched in the response to cytokines. Overall, human splenic NK cells, like mouse splenic NK cells, appeared to be highly functional.

We then compared the gene expression profiles of our four human spleen NK subsets with those of CD56bright and CD56dim NK cells (Hanna et al., 2004). An analysis of module scores calculated at the single-cell level showed that hNK_Sp1 and hNK_Sp4 cells were more like CD56dim NK cells (Figure 5F, left), which is consistent with FCGR3A being one of their driving genes (Figure 5C). The genes strongly expressed in both hNK_Sp1 and CD56dim NK cells were GZMB, FGFBP2, GZMH, CCL4, CCL3, and S1PR5. hNK_Sp2 and 3 cells were more like CD56bright NK cells (Figure 5F, left), and the genes strongly expressed in all three of these populations were GZMK, COTL1, and GATA3. As for the other NK cell subsets, we observed a continuum between the cells most like CD56bright NK cells and those most like CD56dim NK cells (Figure 5F, left). Violin plots of each module score revealed a significant enrichment in the gene signature of CD56dim cells for hNK_Sp1 and hNK_Sp4 cells (Figure 5F, middle) and in the signature of CD56bright cells for hNK_Sp2 and hNK_Sp3 cells (Figure 5F, right). These data, therefore, suggest that hNK_Sp1 cells resemble CD56dim cells and that hNK_Sp2 resemble CD56bright cells, whereas hNK_Sp3 and hNK_Sp4 cells may represent minor subsets of CD56bright and CD56dim cells, respectively. Although CD160 has been originally described as being expressed only in CD56dim NK cells in the blood, flow cytometry analysis of splenic CD56bright NK cells distinguished between hNK_Sp2 and 3 on the basis of the mutually exclusive expression of CD160 and CD52 (CAMPATH-1), respectively, as predicted from our transcriptomic data (Figures S7A and S7B).

We also performed a pan-genomic module score analysis on our human splenic NK cells, comparing them with tonsil NK and ILC signatures in humans (Björklund et al., 2016). Only 10 of the 4,203 NK cells analyzed had a signature more closely resembling the ILC signature than the NK cell signature (Figure S7E). Therefore, as with our analysis of mouse splenic NK cells, only a few cells might have corresponded to ILC1s and did not interfere with our analysis.

High-Throughput scRNA-Seq Identifies Two Splenic-like NK Cell Subsets in Human Blood

As for the human spleen samples, we performed an unsupervised hierarchical analysis of the 1,166 blood NK cells from one donor (data not shown). This analysis identified two blood NK cell subsets, which we named hNK_Bl1 and hNK_Bl2. These two subsets were represented on a t-SNE plot (Figure 6A). Fifty-five genes (33 more strongly expressed in hNK_Bl1 and 22 in hNK_Bl2) displayed significant differential expression between the two subsets and were used as markers for all the human blood NK cell samples (Table S3). These genes produced distinctive transcriptomic signatures of NK cell subsets in human blood (Figure 6B).

Figure 6

High-Throughput scRNA-Seq Identifies Two Splenic-like NK Cell Subsets in Human Blood

(A) t-SNE plot of 1,166 human blood NK cells representative of one sample.

(B) Heatmap of the 55 genes from one representative individual tested with a Wilcoxon rank-sum test discriminating the total 3,200 human blood NK cells into subsets. Squares identify specific transcriptomic signatures of human blood NK cell subsets.

(C) Left: PCA of the two human blood NK cell subsets. Right: driving genes.

(D) Top ten genes significantly differentiating between the two blood NK cell subsets in one representative sample.

(E) Selected Gene Ontology terms.

(F) Left: Module scores of CD56bright and CD56dim gene expression programs. Middle and right: Violin plots representing the distribution of module scores for CD56dim (middle) and CD56bright (right) cells for each blood NK cell, grouped by subset. Wilcoxon rank-sum test with continuity correction. Error bars indicated the mean (±SD). ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

(G) Left: PCA on the two blood and three splenic NK cell subsets. Right: driver genes.

(H) Heatmap of the 290 genes from one representative individual tested with a Wilcoxon rank-sum test separating the total 7,403 human blood and splenic NK cells into subsets. Squares identify specific transcriptomic signatures of the human NK cell subsets.

PC analysis did not differentiate between NK cell subsets from human blood (Figure 6C, left). However, PC1 expressed most of this difference (52%), tending to segregate the cell subsets, whereas PC2 (31%), which separated donors, did not. Blood NK cell populations therefore appeared to be more heterogeneous than paired splenic NK cell populations from the same donor. No obvious clustering was observed between subsets from different donors, but CCL4, CCL3, GZMA, GZMB, PRF1, LGALS1, SPON2, CYBA, CFL1, and S100A4 were identified as the driver genes of PC1. CCL4, CCL3, GZMA, GZMB, PRF1, and S100A4 were also driver genes separating hNK_Sp1 and hNK_Sp4 cells from the hNK_Sp2 and hNK_Sp3 cells. These data indicate that these seven genes are the drivers of human NK cell subpopulations (Figures 5C, 6C, and S6). We compared the top ten differentially expressed genes in these two human blood NK cell subsets, and the top ten genes encoding secreted proteins, cell membrane proteins, and transcription factors (Figure 6D). Like hNK_Sp1 cells, hNK_Bl1 cells displayed a higher expression of FGFBP2, GZMB, SPON2, CST7, LAIR2, PRF1, and FCGR3A. Like hNK_Sp2 cells, the hNK_Bl2 subset displayed differential expression of GZMK, XCL1, and COTL1. Biological process analysis showed an enrichment in cytolysis in hNK_Bl1 cells, whereas the hNK_Bl2 subset displayed significant enrichment in the response to cytokine and regulation of cytokine production functions (Figure 6E).

Module score analysis revealed that the gene expression profile of hNK_Bl1 cells was similar to that of CD56dim NK cells (Figure 6F, left; Hanna et al., 2004). The genes identified as strongly expressed in both CD56dim NK cells and hNK_Bl1 cells were GZMB, FGFBP2, GZMH, and FCGR3A. The gene expression profile of hNK_Sp2 cells was similar to that of CD56bright NK cells (Figure 6F, left). GZMK, LTB, COTL1, XCL1, and CD44 were strongly expressed in both subsets. We noted a continuum with intermediate expression profiles between the cells most like CD56bright NK cells and those most like CD56dim NK cells (Figure 6F, left). Violin plots of each module score between subsets revealed a significant enrichment in the gene expression program of CD56dim cells for hNK_Bl1 cells (Figure 6F, middle) and in the gene expression program of CD56bright cells for hNK_Bl2 cells (Figure 6F, right). At the protein level, differential cell surface densities of CD99 on human blood NK cells, with a higher expression on CD56dim NK cells, were consistent with our transcriptomic analysis (Figures S7C and S7D).

For the identification of genes driving the similarities between hNK_Bl1 and hNK_Sp1 subsets and between hNK_Bl2 and hNK_Sp2 subsets, we performed PC analysis on all human blood and spleen samples. PC1 separated CD56bright-like cells (hNK_Bl2 and hNK_Sp2) from CD56dim-like cells (hNK_Bl1 and hNK_Sp1), regardless of organ source, whereas PC2 distinguished between blood and splenic NK cells (Figure 6G, left and right). Heatmap representation of the 7,403 human blood and splenic NK cells across the 249 differentially expressed genes in the NK subsets highlighted the correspondence between hNK_Bl1 and hNK_Sp1 and between hNK_Bl2 and hNK_Sp2 cells (Figure 6H, upper red and blue rectangles). This analysis also revealed the tissue-specific signatures of these subsets (Figure 6H, lower red and blue rectangles). By contrast, hNK_Sp3 and hNK_Sp4 cells did not resemble any of the blood NK subsets (Figure 6H) and thus appear to be unique to the spleen.

High-Throughput scRNA-Seq Reveals Transcriptomic Signatures Common to Organs and Species

We investigated possible functional relationships between the organ-specific populations of NK cells from mice and humans, by comparing organ-specific orthologs for each organ (Figures 1 and 4). We identified, among the mouse and human genes, 37 orthologous genes preferentially expressed in spleen NK cells and 8 orthologous genes preferentially expressed in blood NK cells (Figure 7A, right). The 37 spleen NK cell genes included those encoding secreted effector proteins (PRF1, IFNG, CCL3, and XCL1), proteins involved in transcriptional regulation (e.g., JUN, KLF6, and FOSB), and membrane proteins (e.g., KLRC1, ICAM1, and CD69). One of the eight genes with orthologs in the blood NK cells of the two species encoded a secreted protein (TGFB1), four encoded transcriptional regulators (MXD4, CEBPB, ARID5A, and HHEX), and the others encoded lipid-binding adaptor proteins (PLEKHA5) and proteins involved in posttranslational protein regulation (PTPN7 and RNF125).

Figure 7

High-Throughput scRNA-Seq Reveals Transcriptomic Signatures Common to Organs and Species

(A) Spleen-specific (left) and blood-specific (right) transcriptomic signature common to human and mouse NK cells.

(B) Subset-specific transcriptomic signature common to human and mouse splenic NK cells.

(C) Subset-specific transcriptomic signature common to human and mouse blood NK cells.

(B and C) Underlined genes are common to the blood and spleen subsets.

(D) Schematic representation of mouse and human spleen and blood NK cells subsets based on scRNA-seq analysis.

We compared NK cell heterogeneity between mice and humans by performing pairwise comparisons of the organ-specific orthologs for each NK subset (Figures 7B and 7C). For spleen NK cells, subset 1 from both species (mNK_Sp1 and hNK_Sp1) were related and had a 9-gene signature in common, with a higher expression than the other subsets for genes encoding the NK cell effector proteins GZMB, PRF1, and LGALS1, genes encoding cell membrane proteins such as EMP3, ITGB2, and S1P5, and the gene encoding the transcriptional regulator ZEB2, which is involved in terminal NK cell maturation (Figure 7B). Likewise, human and mouse splenic subset 2 (mNKSp2 and hNK_Sp2) had a 6-gene signature in common, including XCL1, CD7, and CD160. With the exception of JUN between mNK_Sp3 and hNK_Sp4 cells, no orthologous gene with differential expression related the human splenic subsets 3 and 4 with any of the mouse splenic subsets (Figure 7B). For blood NK cells, there were fewer orthologous genes as compared to the spleen, but hNK_Bl1 cells resembled mNK_Bl1 cells, through a common 7-gene signature including GZMB, PRF1, EMP3, ITGB2, and ZEB2 (Figure 7C). hNK_Bl2 and mNK_Bl2 cells had a common signature of only three genes: XCL1, LTB, and RPL36A (encoding a ribosomal protein) (Figure 7C). A 5-gene signature (GZMB, PRF1, EMP3, ITGB2, and ZEB2) for cytolytic subset 1 and XCL1 for subset 2 were found to be shared from both spleen and blood (underlined in Figures 7B and 7C).

Finally, we summarized the correspondence between subset 1 of NK cells in the mouse (CD27−CD11b+) and in human (CD56dim) and between subset 2 of NK cells in the mouse (CD27+CD11b−) and in human (CD56bright) (Figure 7D).

Discussion

The main objective of this work was to assess, from an unbiased transcriptome-wide perspective, the similarities and differences between human and mouse NK cells and between blood and spleen NK cells. We defined the heterogeneity of human and mouse NK cells through scRNA-seq on thousands of individual NK cells. The unsupervised clustering distinguished splenic from blood NK cells in both species. Further unsupervised analysis in each organ revealed three subsets in mouse spleen, four subsets in human spleen, and two subsets each in mouse and human blood. Moreover, these analyses defined two main NK cell subsets across organs and species, which we refer to here as NK1 and NK2.

Both NK1 and NK2 were present in the blood and spleen. This suggests that splenic subsets may originate from blood subsets, because blood NK cells are thought to recirculate continually throughout the body (Lysakova-Devine and O’Farrelly, 2014). We identified several genes and cell surface proteins that defined the NK1 and NK2 subsets. GO term enrichment analysis supported the cytotoxic activity of the NK1 subset consistent with earlier studies showing that blood human CD56dimCD27lo NK cells have strong cytotoxic activity (Vossen et al., 2008). The NK2 subset was characterized by the conserved expression of a single gene across organs and species, XCL1, encoding an chemoattractant produced by NK cells that recruits dendritic cells to promote tumor control (Chiossone et al., 2018, Bottcher et al., 2018).

Human and mouse splenic NK cells shared 37 orthologous genes, including genes involved in tuning NK cell responses, such as those encoding the effector proteins CCL3, IFNG, PRF1, and XCL1, the cell membrane proteins CD160, CD69, and NKG2A (KLRC1), and the transcription factors KLF2 and KLF6. By contrast, blood NK cells from humans and mice had only eight differentially expressed orthologous genes in common, including genes encoding proteins involved in immune regulation, such as the protein tyrosine phosphatase PTPN7, and TGFB1. Of note, tissue-resident human NK cells have been reported to express CD69 (Lugthart et al., 2016), which was part of the conserved splenic NK cell-specific gene signature.

Mouse blood NK cells were enriched in a few biological processes and no enrichment in any specific process was noted for human blood NK cells. By contrast, splenic NK cells from both species were enriched in response to stimulus, response to stress, defense response, signal transduction, and cell activation. Therefore, despite the presence of similar major NK cell subsets in the spleen and blood, splenic NK cells appeared to be functionally more active than blood NK cells in both species. Consistently, mouse splenic NK cells displayed greater reactivity than paired blood NK cells upon stimulation. Thus, our organ-specific gene signature of NK cells raises the possibility that organ-specific factors shape NK cell phenotype and functions.

In addition, we identified species-specific subsets in the spleen—two in humans and one in mice. In humans, each of the additional splenic subsets had some of their gene signatures in common with NK1 and NK2 cells. The hNK_Sp3 subset was more like CD56bright cells (NK2), with significantly a higher expression of the gene encoding the homing receptor CD62L (SELL). By contrast, the hNK_Sp4 subset had many driving genes in common with CD56dim cells (NK1). Given these similarities, we hypothesize that these two human splenic subsets represent minor CD56bright and CD56dim subsets, respectively. In mouse spleen, we identified a third minor population of NK cells different from the CD27+CD11b−-like and CD27−CD11b+-like cells. mNK_Sp3 was less represented in the spleen, but its gene signature contributed greatly to the differences between organs. These three spleen-specific subsets constitute previously undefined NK cell subsets in mice and humans and provide additional evidence of the efficacy of scRNA-seq approaches for studying cell heterogeneity. Further characterization based on surface marker phenotyping, proteomics analysis, and functional assays is required to define the maturation and functional states of these cells.

Of note, we could not identify CD27+CD11b+ double-positive (DP) cells, a transient maturation stage, by scRNA-seq. Each DP cell could potentially cluster with either the CD11b− or CD27− subsets, depending on its position in the trajectory of NK cell maturation. A recent study challenged the typical cell-surface phenotype of adaptive NK cells (Hammer et al., 2018). A separate study would thus be required to analyze the heterogeneity of adaptive NK cells at the single-cell resolution.

ILC1s exhibit common attributes with NK cells (Spits et al., 2016). Our organ-matched module score analysis showed that only a few cells corresponded to ILC1s. They did not interfere with our analysis.

Overall, our study strengthens the dichotomy between the CD27−D11b+ (NK1) and CD27+CD11b− (NK2) subsets in mice and the CD56dim (NK1) and CD56bright (NK2) subsets in humans, showing the correspondence between these subsets across species and documents the presence of previously unknown populations in the mouse and human spleen. This unbiased analysis of NK cell transcriptomic signatures also reveals spleen- and blood-specific NK cell signatures common in both species, highlighting the importance of the organ of origin in the definition of a cell population. Beyond, this scRNA-seq analysis should provide a useful resource for the translation of mouse NK cell studies into information about the role of these cells in human physiology and disease.

STAR★Methods

Key Resources Table

Contact for Reagent and Resource Sharing

Further information and requests should be directed to and will be fulfilled by Lead Contact, Eric Vivier (vivier@ciml.univ-mrs.fr).

Experimental Model Details

Control subjects

Peripheral blood and spleen biopsy specimens were obtained from brain-dead beating-heart donors eligible for liver donation, at the end of the graft harvesting procedure. The hepatic vasculature was flushed, removing most of circulating NK cells from tissues. Organ harvesting was performed by the surgeons of the team of Prof. Jean Hardwigsen from the General Surgery and Hepatic Transplantation Department of Timone Hospital for Adults (Marseille, France). The donors were two women (donors 1 and 3) of 73 and 72 years of age, both of whom died from hemorrhagic stroke, and one 50-year-old man (donor 2) who died from sudden cardiac death.

Donor 1 had antibodies consistent with past EBV, CMV and toxoplasmosis infections, and suffered from type 2 diabetes mellitus and high blood pressure, for which she was treated with metformin and irbestan. Donor 2 had antibodies consistent with past EBV, CMV and toxoplasmosis infections and was on no current treatment. Donor 3 had antibodies consistent with past EBV and, CMV infection, suffered from a stroke in 2015, and had type 2 diabetes mellitus, hypercholesterolemia, and high blood pressure, for which she was treated with metformin, pravastatin and atenolol. She also took aspirin for rheumatism. This study was performed in accordance with French rules (Art. L.1243-1 and Art. L.1245-2 of the French Public Health Code). The study was approved by the French Biomedicine Agency (Agence de la biomedicine) under approval number PFS16-009. All samples were stored in IGL-1 medium, at 4°C, and were processed as soon as they were received.

Mice

All mice were obtained from Janvier and were maintained in pathogen-free facilities at CIML. Seven-week-old male C57BL/6j mice were used. Five mice were pooled per sample. All animal experiments were performed with the approval of the Animal Ethics Committees of the CIML and in accordance with European laws.

Method Details

Cell preparation and flow cytometry for scRNaseq

Human blood samples were centrifuged on a Ficoll gradient to obtain a preparation enriched in peripheral blood mononuclear cells (PBMCs), which was then frozen in 90% decomplemented SVF + 10% DMSO, as previously described (Tomasello et al., 2012). Human spleen biopsy specimens were weighed, washed with PBS and mechanically disrupted. The preparation was enriched in splenocytes on a Ficoll gradient and frozen. Mouse blood and spleen samples were obtained by intra-cardiac puncture and during surgery, respectively. PBMCs were obtained by lympholyte-based centrifugation followed by red blood cell lysis. Mouse spleens were weighed, washed with DPBS and mechanically disrupted. Cells were filtered and centrifuged and subjected to red blood cell lysis. Samples 1 and 3 from human donors were thawed and processed on the same day as samples 1 and 3 were harvested from mice. Samples 2 (mice and human donors) were processed on another day.

NK cell sorting

Briefly, individual human samples were thawed, washed and incubated for 10 min at 4°C in DPBS + 5% FCS + 2 mM EDTA with 2% normal mouse serum. Cells were then stained with anti-CD3, -CD14, - CD19, -CD45 and anti-CD56 antibodies for 30 min at 4°C. Cells were washed twice with DPBS and incubated in DPBS with dead cell marker for 10 min at 4°C. Freshly prepared mouse NK cells were first incubated for 10 min at 4°C in DPBS + 5% FCS + 2 mM EDTA with 1:50 mouse Fc Block and then stained for 45 min at 4°C with anti-CD3, -CD19, -CD45.2, -NK1.1 and anti-NKp46 antibodies. Cells were washed twice with DPBS and incubated in DPBS with dead cell marker for 10 min at 4°C. Mouse and human cells were washed and sorted on an Influx cell sorter from BD (Becton Dickinson, San Diego, USA).

Flow cytometry to assess cell-surface phenotype

Mouse blood was collected by flushing the whole vasculature. Matched spleens were then harvested. PBMCs were obtained by lympholyte-based centrifugation followed by red blood cell lysis. Mouse spleens were weighed, washed with DPBS and mechanically disrupted. Cells were filtered and centrifuged and red blood cells were lysed. Staining was performed on the same day, with acquisition on a LSR Fortessa from BD (Becton Dickinson, San Diego, USA). Frozen PBMCs from healthy human individuals and frozen splenic cells from brain-dead donors were thawed and stained on the same day, with acquisition on a LSR II cytometer with UV-configuration from BD (Becton Dickinson, San Diego, USA). Statistical significances were calculated using Friedmann analysis with Dunn post hoc tests on paired subsets measurements, or a Dunn test and p-values were adjusted with the Benjamini-Hochberg method.

∗ p-value < 0.05, ∗∗ p-value < 0.01, ∗∗∗ p-value < 0.001, ∗∗∗∗ p-value < 0.0001.

Single-cell RNA sequencing

After sorting, cells were washed in DPBS + 0.04% BSA, as recommended by the sample preparation protocol of 10x Genomics, and kept on ice before counting. Each mouse sample was pooled with a human sample from the same organ, in a 1:1 ratio for accurate determination of the doublet rate. We used the 10x Genomics Chromium single-cell 3′ v2 kit and protocol to prepare the libraries. HalioDx (Marseille, France) performed the RNA sequencing on a NextSeq500 machine with a sequencing depth of at least 50,000 reads per cell.

NK cell activation

Paired mouse splenocytes and blood cells were stimulated with 5,000 units of IL-2 for 5 days. Cells were then cocultured for 4 hours with YAC-1 in the presence of monensin and brefeldin and CD107a, and surface-stained and processed for the detection of intracellular IFN-γ. 2HB Immulon plates were coated with antibodies against NKp46 or NK1.1 at a concentration of 10 μg/mL. Cells were incubated for four hours in the presence of monensin, brefeldin and CD107a, surface-stained and processed for intracellular IFN-γ detection.

Quantification and Statistical Analysis

Preprocessing of mouse and human samples

Raw FASTQ files were processed with CellRanger software (v1.3.1), which performs alignment, filtering, barcode counting and UMI counting. Reads were first simultaneously aligned with the mm10 and GRch38 genomes, with CellRanger software (v1.3.1), for their assignment as murine or human sequences. Species assignments were improved by individual alignment with the mm10 and GRch38 sequences, resulting in two clusters, one for each species. Finally, we took the intersection between the corresponding clusters of the two alignments and we removed cells identified as doublets (within a range of 1.5 to 3.8%). Low-quality cells were excluded in an initial quality-control (QC) step, by removing genes expressed in fewer than three cells, cells with less than 200 genes expressed, and cells expressing more than 2500 genes. Cells with more than 0.05% of mitochondrial-associated genes among their expressed genes were also removed. In total, 8,193 mouse cells and 8,471 human cells were retained for further analysis.

Library-size normalization was performed on the UMI-collapsed gene expression values for each cell barcode, by scaling by the total number of transcripts and multiplying by 10,000. The data were then log-transformed before further downstream analysis with Seurat (Satija et al., 2015).

Mouse sample analysis

We first selected genes with a high variance, by using the FindVariableGenes function with log-mean expression values between 0.0125 and 8 and a dispersion (variance/mean) between 0.5 and 8. This selected 1,433, 1,759, and 1,081 genes for blood samples and 1,223, 1,109, and 1,610 genes for spleen samples 1, 2, and 3, respectively. We then reduced the dimensionality of our data by principle component analysis (PCA) and we identified by random sampling 20 significant PCs for each sample, except one, the spleen of individual 3, for which 10 PCs were selected, with the PCElbowPlot function. Cells were clustered with Seurat’s FindClusters function. In total, 8,118 of the total of 8,193 cells were assigned to NK cell subsets. A range of cluster resolution parameters were tested: 0.7 was retained for all but one sample, the spleen of sample 3, for which a value of 0.5 was applied. The resulting clusters were reproducible and there were no subdivisions into clusters of very small numbers of cells. We used the percentage of doublets as a minimum threshold for cluster retention. For visualization, we applied RunTsne to the cell loadings of the previously selected PCs with default perplexity, to view the cells in two dimensions. We identified cluster markers with FindMarkers, with the parameter only.pos set to TRUE, to obtain only upregulated genes as markers of a relative to all other cells. Marker genes were defined as genes with an adjusted p-value < 0.05 tested with a non-parametric Wilcoxon rank-sum test.

Human sample analysis

We selected the sample with the best balance between the number of cells and the number of genes as a gold standard for each tissue. We performed the analysis described above on this sample (233 markers, 20 PCs and resolution of 0.6 for the spleen sample and 55 markers, 20 PCs and a resolution of 0.4 for the blood sample). We used the markers identified from both the blood and spleen to perform clustering on other human samples. Subsets were assigned by comparison with the first sample on the basis of the expressed driving genes from PCA.

We removed CCL3, CCL4, CCL3L3, and CCL4L2 from the list of genes displaying variable expression in blood samples, because they were expressed higher in only one sample, completely driving PC1. We removed 947 of the 1498 cells (from the cells expressing highly these genes) from the same blood sample, because they appeared to be highly patient-specific (sample 3). We removed 87 of the 1607 cells from the second blood sample, because they were found to be substantially different from the other cells in the sample on the basis of t-SNE representation and PCA. The cells removed were probably contaminating T cells. We retained 7403 of the 8471 total cells for the human NK cell subsets.

Pooled sample analysis

For the identification of tissue-specific markers, we pooled samples for each species together with the CellRanger aggregate function and normalization, to equalize read depth between libraries. We then analyzed pooled samples as described above.

Assessment of NK cell purity

We checked the purity of the NK cell preparation analyzed, by checking our transcriptomic results for genes other than those specific to the NK cell lineage, such as those specific for hematopoietic stem cells (human: CD34; mouse: Sca1), myeloid cells (human: CD14, ITGAX; mouse: Itgax, Adgre1), granulocytes (mouse: Ly6G), T cells (human: TRAC, CD8B; mouse: Cd3, Cd4), and B cells (human: CD19; mouse: Cd19). None of the sorted NK cell populations analyzed included immune cells from other cell lineages.

Unsupervised hierarchical clustering

With spleen or blood samples alone: We identified the intersection of the genes displaying variable expression between individual samples for each cluster. The mean expression of the intersection of these genes between clusters was then calculated and hierarchical clustering was performed with these values.

With mixed spleen and blood samples: We identified the intersection of the genes with variable expression between individual samples for each organ. Mean expression was calculated for the union of these genes, between tissues, and hierarchical clustering was performed with the values obtained.

For all samples, a Euclidean distance was used for genes and clusters.

Principal component analysis

PC component clustering analysis was performed on the mean expression of all genes between clusters (even those removed from the clustering analysis). PCA gene loadings for PC1 and PC2 corresponded to the genes making the largest contributions, accounting for 40% of the total information for each PC.

Heatmap

Heatmaps were generated from the scaled expression (log TPM values) values for the discriminating gene sets defining each subset, with adjusted p-value < 0.05 with a non-parametric Wilcoxon rank-sum test. The color scale is based on the z-score distribution, from −2 (purple) to 2 (yellow).

Gene annotations

Cell membrane, secreted 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

Enrichment scores (p-values) for selected numbers of GO annotations were calculated with a hyper-geometrical statistical test with a threshold of 0.05. The data were plotted as -log10 of the p-value after Benjamini and Hochberg correction. The threshold of significance was set at -log10(0.05). For mouse data, enrichment was calculated for the union of the marker genes across samples. The background was the union of all genes expressed in mouse samples. For human data, enrichment was calculated for the marker genes from one gold-standard sample. The background was all the genes expressed in this sample. Samples were ranked according to the lowest p-value of the leader subsets for each GO annotation.

Comparison with known NK cell subsets

Module scores for ILC1 signatures (Robinette et al., 2015), CD27- CD11b+ and CD27+ CD11b- (mouse) or CD56bright and CD56dim (human) gene expression programs defined by Chiossone et al. 2009 or Hanna et al. 2004, respectively, were determined with AddModuleScore from Seurat for each of our NK cells at the single-cell level. Briefly, the mean expression for each gene in the defined expression programs was calculated for each cell and the aggregated expression of control gene sets was then subtracted. All analyzed genes were binned on the basis of mean expression, and the control genes were randomly selected from each bin. A violin plot representation was used to assess the distribution of module scores for each NK cell, grouped by subset. Statistical significance was calculated using a Wilcoxon rank sum test with continuity correction or using a Kruskal-Wallis with Dunn’s multiple-comparison 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.

Ortholog-based comparison of human and mouse transcriptomic signatures

We compared human and mouse organ-specific or subset-specific transcriptomic signatures by manually transforming mouse data into human orthologs according to the information in public databases (NCBI). Genes on the gene lists of both species were identified by pairwise comparison.

Assessment of the robustness of human NK cell subsets in a machine-learning protocol

A machine-learning approach proposed by Valentine Svensson was used (http://www.nxn.se/valent/2018/3/5/actionable-scrna-seq-clusters). We performed the analysis by logistic regression, using 50% of our cells for training and 50% for testing. Receiver operating characteristic (ROC) curves were constructed, using the false-positive rate (x axis) and the true positive rate (y axis) for all possible thresholds of probabilities given by the logistic regression. The area under the curve (AUC) was then calculated to assess the quality of cluster assignment, and, thus, the strength of subsets identification.

Data and Software Availability

All sequencing data have been deposited at NCBI GEO depository and are accessible with the accession number GEO: GSE119562 or using the link https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE119562.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Anti-human CD127 clone HIL-7R-M21, Pe-Cy7	BD Biosciences	cat # 560822; RRID: AB_2033938
Anti-human CD14 clone M5E2, BUV737	BD Biosciences	cat # 564444; RRID: NA
Anti-human CD16 clone 3G8, BUV496	BD Biosciences	cat # 564653; RRID: NA
Anti-human CD160 clone BY55, Alexa Fluor® 488	BD Biosciences	cat # 562351; RRID: AB_11153688
Anti-human CD161 clone DX12, Pe-Cy5	BD Biosciences	cat # 551138; RRID: AB_394068
Anti-human CD18 clone L130, BV395	BD Biosciences	cat # 744556; RRID: NA
Anti-human CD19 clone SJ25C1, BUV737	BD Biosciences	cat # 564304; RRID: AB_2716867
Anti-human CD3 clone UCHT1, BUV737	BD Biosciences	cat # 564308; RRID: NA
Anti-human CD45 clone 2D1, APC-H7	BD Biosciences	cat # 560178; RRID: AB_1645479
Anti-human CD52 clone 4C8, Alexa Fluor® 647	BD Biosciences	cat # 563610; RRID: NA
Anti-human CD82 clone 423524, Alexa Fluor® 647	BD Biosciences	cat # 564341; RRID: NA
Anti-human CD83 clone HB15e, BV421	BD Biosciences	cat # 566264; RRID: NA
Anti-human CD99 clone TÜ12, BV786	BD Biosciences	cat # 743045; RRID: NA
Anti-human NKp46 clone 9E 2, BV605	BD Biosciences	cat # 743710; RRID: NA
Anti-mouse CD107a clone 1D4B FITC	BD Biosciences	cat # 553793; RRID: AB_395057
Anti-mouse CD11b clone M1/70, PE-CF594	BD Biosciences	cat # 562287; RRID: AB_11154216
Anti-mouse CD19 clone 1D3, FITC	BD Biosciences	cat # 557398; RRID: AB_396681
Anti-mouse CD3 clone 145-2C11, FITC	BD Biosciences	cat # 553062; RRID: AB_394595
Anti-mouse CD3 clone 145-2C11, PE-Cy5	BD Biosciences	cat # 553065; RRID: AB_394598
Anti-mouse CD3 clone 500A2 V500	BD Biosciences	cat # 560771; RRID: AB_1937314
Anti-mouse CD45.2 clone 104, Alexa Fluor® 700	BD Biosciences	cat # 560693; RRID: AB_1727491
Anti-mouse CD45.2 clone 104, BV786	BD Biosciences	cat # 563686; RRID: NA
Anti-mouse IFNγclone B27 APC	BD Biosciences	cat # 554702; RRID: AB_398580
Anti-mouse IFNγclone XMG1.2 BV421	BD Biosciences	cat # 563376; RRID: NA
Anti-mouse KLRG1 clone 2F1, BUV395	BD Biosciences	cat # 740279; RRID: NA
Anti-mouse NK1.1 clone PK136 PerCP-Cy™5.5	BD Biosciences	cat # 561111; RRID: AB_10564092
Anti-mouse NK1.1 clone PK136, BV510	BD Biosciences	cat # 563096; RRID: NA
Anti-mouse NK1.1 clone PK136, PE-TR	BD Biosciences	cat # 562864; RRID: NA
Anti-mouse NKp46 clone 29A1.4, BV421	BD Biosciences	cat # 562850; RRID: NA
Anti-human CD14 clone RMO52, FITC	Beckman Coulter	cat # IM0645U; RRID: AB_130992
Anti-human CD19 clone J3-119, FITC	Beckman Coulter	cat # A07768; RRID: NA
Anti-human CD3 clone UCHT1, FITC	Beckman Coulter	cat # A07746; RRID: NA
Anti-human CD56 clone N901(NKH-1), PE	Beckman Coulter	cat # A07788; RRID: AB_2636814
Anti-human CD45 clone HI30, APC-Cy7	Biolegend	cat # 304014; RRID: AB_314402
Anti-human CX3CR1 clone 2A9-1, BV421	Biolegend	cat # 341620; RRID: AB_2687148
Anti-mouse CD27 clone LG.3A10, PE-Cy7	Biolegend	cat # 124216; RRID: AB_10639726
Anti-mouse CD28 clone E18, FITC	Biolegend	cat # 122007; RRID: AB_604062
Anti-mouse CD3 Alexa Fluor® 647 clone 17A2	Biolegend	cat # 100209; RRID: AB_389323
Anti-mouse CD90.2 clone 30-H12, BV785	Biolegend	cat # 105331; RRID: AB_2562900
Anti-mouse ITGB2 clone H155-78, APC	Biolegend	cat # 141009; RRID: AB_2564305
Anti-mouse NKp46 clone 29A1.4 purified	Biolegend	cat # 137601; RRID: AB_10551441
Anti-mouse CD19 clone 1D3, PE-Cy5	eBioscience	cat # 15-0193-83; RRID: AB_657673
Anti-mouse NKp46 clone 29A1.4, PE-Cy7	eBioscience	cat # 25-3351-82; RRID: AB_2573442
Anti-mouse NK1.1 clone PK136 purified	eBioscience	cat # 14-5941-81; RRID: AB_467735
Anti-mouse NKp46 clone 29A1.4 PerCP-eFluor 710	eBioscience	cat # 46-3351-80; RRID: AB_1834442
Anti-mouse CCR2 clone 475301, PE	RD Systems	cat # FAB5538P; RRID: AB_10718414
Chemicals, Peptides, and Recombinant Proteins
BD Golgi Plug (Brefeldin A solution)	BD Biosciences	cat # 555029
BD Golgi stop (Monensin solution)	BD Biosciences	cat # 554724
Mouse FcBlock	BD Biosciences	cat # 553142; RRID: AB_394657
Lympholyte mammals	Cedarline	cat # CL5120
DMSO	Eurobio	cat # CP-10
Lymphocyte separation medium	Eurobio	cat # CMSMSL01-0U
DPBS	GIBCO	cat # 14200-067
Fetal calf serum	GIBCO	cat # 10270-106
HEPES	GIBCO	cat # 15630056
L-Glutamine	GIBCO	cat # 25030081
Penicillin-Streptomycin	GIBCO	cat # 15140122
Sodium Pyruvate	GIBCO	cat # 11360039
RPMI	GIBCO	cat # 21875-034
BSA	ID Bio	cat # 1000-70
Dead cell marker, BV510	Invitrogen	cat # L34957
Dead cell marker, UV	Invitrogen	cat # L23105
EDTA	Invitrogen	cat # 15575-038
Red blood cell lysis	Invitrogen	cat # 00-4333-57
Proleukin® Human IL-2 (Aldesleukine)	Merck	cat # NA; RRID: NA
Normal mouse serum	Sigma-Aldrich	cat # M5905
Deposited Data
scRNA-seq data	This paper	GEO: GSE119562
Software and Algorithms
ade4	N/A	https://github.com/cran/ade4/
biomaRt	Durinck et al., 2009	https://github.com/grimbough/biomaRt
Cell Ranger	10x Genomics	https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger
Gostats	Falcon and Gentleman, 2007	https://www.bioconductor.org/packages/release/bioc/html/GOstats.html
R	R Development Core Team, 2008	https://www.r-project.org/
R Studio	N/A	https://www.rstudio.com/
Seurat	Satija et al., 2015	https://github.com/satijalab/seurat/