Landscape of innate lymphoid cells in human head and neck cancer reveals divergent NK cell states in the tumor microenvironment

Natural killer (NK) cells comprise one subset of the innate lym- phoid cell (ILC) family. Despite reported antitumor functions of NK cells, their tangible contribution to tumor control in humans remains controversial. This is due to incomplete understanding of the NK cell states within the tumor microenvironment (TME). Here, we demonstrate that peripheral circulating NK cells differentiate down two divergent pathways within the TME, resulting in different end states. One resembles intraepithelial ILC1s (ieILC1) and possesses potent in vivo antitumor activity. The other ex- presses genes associated with immune hyporesponsiveness and has poor antitumor functional capacity. Interleukin-15 (IL-15) and direct contact between the tumor cells and NK cells are required for the differentiation into CD49a + CD103 + cells, resembling ieILC1s. These data explain the similarity between ieILC1s and tissue- resident NK cells, provide insight into the origin of ieILC1s, and identify the ieILC1-like cell state within the TME to be the NK cell phenotype with the greatest antitumor activity. Because the pro-portions of the different ILC states vary between tumors, these findings provide a resource for the clinical study of innate immune responses against tumors and the design of novel therapy. function These findings indicate that a better understanding of the different ILC states within human tumors will identify key immune cell subsets that have important influence on tumor growth and behavior. recently the role of ILCs in human begun studied, is The evidence, in mice, that peripheral NK cells can convert into ILC1-like cells within the tumor microenvironment 33) has not been reported in humans. In this study, we investigated the landscape of the ILC subsets in human HNSCC by single-cell RNA sequencing (scRNA-seq) of sorted ILCs from primary tu- mors. We describe a spectrum of ILC states that are present within the tumors and demonstrate that peripheral circulating NK cells differentiate into two divergent terminal states in the presence of tumor cells and IL-15. One state is an NK cell subset that resembles CD49a + CD103 + ieILC1s with potent antitumor functional capacity, and the other is an NR4A -expressing CD49a – NK cell subset with poor antitumor functional capacity. The understanding of these two NK cell states within human tumors may have important implications for prognostication and ther- apy, and the results from this study provide an important resource for the future investigation of NK cells in human cancer. alteredstates ofNK cells,in ILC1s ILC3s, found. These led to the that NK cells differentiated in the microenvironment, and in using available cell lines and from healthy analyzed for quality on a 5300 Fragment Analyzer with 96 capillary array (Agilent), using High Sensitivity NGS Fragment kit-DNF-474 (Agilent). Plates with a minimum DNA concen- tration of 0.2 ng/ μ L per well (average range 0.5 ng/ μ L to 1.5 ng/ μ L) were used for the final library preparation. First, the cDNA concentration was nor- malized with nuclease free water (IDT; Cat# 11-04-02-01) to a concentration less than or equal to 1 ng per well and transferred to a 384-well low volume serial dilution plate (TTP Labtech; Cat# 4150-05828) using the Mosquito × 1 (TTP Labtech). Then, library preparation steps as per Smart-seq2 protocol were performed using Nextera XT DNA Library Prep Kit (Illumina; FC-131-1096) in a 384-Well Rigi-Plate PCR Microplate (TTP Labtech; Cat# PCR- 384-RGD-C) using Mosquito HTS (TTP Labtech). One-microliter volumes of indexed libraries from each well of the 384-well plate were pooled together and purified using 0.8 × Agencourt AMPure XP Beads (Beckman Coulter; Cat# A63880). Concentration of the purified pooled library was determined using the Bioanalyzer (Agilent). The pooled library was run on the HiSeq 2500 high-output v4 sequencer (Illumina), with 2 × 125 bp read length. One 384-cell pooled library was run on one lane of the HiSeq 2500, with an average output of 1 × 10 6 reads per cell.

Natural killer (NK) cells comprise one subset of the innate lymphoid cell (ILC) family. Despite reported antitumor functions of NK cells, their tangible contribution to tumor control in humans remains controversial. This is due to incomplete understanding of the NK cell states within the tumor microenvironment (TME). Here, we demonstrate that peripheral circulating NK cells differentiate down two divergent pathways within the TME, resulting in different end states. One resembles intraepithelial ILC1s (ieILC1) and possesses potent in vivo antitumor activity. The other expresses genes associated with immune hyporesponsiveness and has poor antitumor functional capacity. Interleukin-15 (IL-15) and direct contact between the tumor cells and NK cells are required for the differentiation into CD49a + CD103 + cells, resembling ieILC1s. These data explain the similarity between ieILC1s and tissueresident NK cells, provide insight into the origin of ieILC1s, and identify the ieILC1-like cell state within the TME to be the NK cell phenotype with the greatest antitumor activity. Because the proportions of the different ILC states vary between tumors, these findings provide a resource for the clinical study of innate immune responses against tumors and the design of novel therapy. ILC | natural killer cells | intratumoral | ieILC1 | HNSCC T he role of NK cells in host protection against human tumor formation and growth has been unclear. Intratumoral NK cells, found in most solid tumors, are often dysfunctional and, under certain conditions, can promote tumor growth (1,2). Conversely, the detection of activated NK cells within tumors, particularly in head and neck squamous cell carcinoma (HNSCC), has also been associated with significantly better clinical outcomes (3)(4)(5). Previous studies examining NK cells within human tumors have focused on narrow definitions based on limited surface marker profiles. The emerging recognition of plasticity between NK cells and other ILC family members has suggested that a broader examination of ILCs within solid tumors may provide insight into the different intratumoral NK cell subsets and their relationships.
The ILC family of lymphocytes plays diverse roles in tissue homeostasis and host defense (6). Although ILCs do not express recombined antigen receptors found on T cells and do not undergo clonal expansion or selection when stimulated, they can be grouped into subsets that are akin to subsets of T cells. Based on the cytokines they produce and the transcription factors they express, ILC subsets are categorized into group 1 ILC, group 2 ILC, and group 3 ILC, which are considered analogous to the CD4 + T helper cell Th1, Th2, and Th17 cell subsets, respectively. Specifically, group I ILC includes ILC1s, intraepithelial ILC1s (ieILC1s) and NK cells, which produce IFNγ and TNFα and are dependent on the transcription factor T-BET, similar to Th1 cells (7)(8)(9)(10). ILC2s produce IL-4, IL-5, and IL-13 and express GATA-3, similar to Th2 cells (11)(12)(13)(14). ILC3s produce IL-17 and IL-22, similar to Th17 cells, and are dependent on the transcription factor RORγt (15)(16)(17).
NK cells are considered to be a member of the group 1 ILC (18), but, unlike ILC1s, NK cells express granzymes and perforin and have cytolytic activity, making them analogous to cytotoxic CD8 + T cells. The expression of the transcription factor EOMES by NK cells, in addition to T-BET, also distinguishes NK cells from ILC1s. Finally, unlike ILC1s which are primarily tissue Significance NK cells have been observed to be present within most solid tumors, but the antitumor activity of the intratumoral NK cells has been unclear. In this study, we examined the entire spectrum of innate lymphoid cells within primary human tumors and demonstrate that peripheral NK cells in the tumor microenvironment differentiate into heterogeneous cell states, resulting in either a hyporesponsive NK cell subset or a highly active NK cell phenotype that closely resembles intraepithelial ILC1s and that has potent antitumor properties. Importantly, this differentiation into ieILC1-like cells occurs when NK cells are cocultured with epithelial cells, providing important information for NK cell immunotherapy approaches. resident, NK cells both circulate in the peripheral blood and can reside within tissue compartments, including solid tumors. Within humans, the conventional NK (cNK) cells that circulate in the blood are primarily CD56 dim CD3 -CD19cells. These peripheral NK cells can be distinguished from tissue-resident NK (trNK) cells, which resemble CD56 bright cells (19) and express a combination of tissue residency markers, including the chemokine receptor CXCR6, integrin α1 (ITGA1, CD49a), and integrin αE (ITGAE, CD103), depending on particular tissue type and location (18,20). Another member of the group 1 ILC is ieILC1 cells, which are characterized by the expression of surface markers CD49a, CD103, and NKp44 (9). Because ieILC1s express both EOMES and T-BET, some have proposed that these cells are more closely related to NK cells than classic helper ILC1s (18,21). The ieILC1s, which reside in mucosal and tonsillar tissue, also express CXCR6, similar to trNK cells described in other tissues such as the bone marrow (BM), spleen, and liver (18,22).
There is evidence for plasticity between ILC subsets, which can be induced by changes in the tissue microenvironment (23)(24)(25)(26)(27)(28)(29). Within murine group 1 ILC, there is evidence for plasticity between NK cells and ILC1-like cells (30)(31)(32)(33), and TGF-β signaling is involved in the conversion of NK cells into ILC1-like cells. However, the significance of this plasticity in disease contexts is not well understood. For instance, in murine transplantable tumor models using B16 melanoma and MCA-induced sarcomas, the TGFβ−mediated conversion of NK cells to an intermediate ILC1 (intILC1) state resulted in poor control of tumor growth (32). In contrast, in the murine MMTV-PyMT mammary tumor model, a CD49a hi CD103 hi ILC1-like population was demonstrated to have a potent immunosurveillance function (31). These findings indicate that a better understanding of the different ILC states within human tumors will identify key immune cell subsets that have important influence on tumor growth and behavior.
It is only recently that the role of ILCs in human cancer has begun to be studied, and there is still very little that is understood (34). The evidence, in mice, that peripheral NK cells can convert into ILC1-like cells within the tumor microenvironment (32,33) has not been reported in humans. In this study, we investigated the landscape of the ILC subsets in human HNSCC by single-cell RNA sequencing (scRNA-seq) of sorted ILCs from primary tumors. We describe a spectrum of ILC states that are present within the tumors and demonstrate that peripheral circulating NK cells differentiate into two divergent terminal states in the presence of tumor cells and IL-15. One state is an NK cell subset that resembles CD49a + CD103 + ieILC1s with potent antitumor functional capacity, and the other is an NR4A-expressing CD49a -NK cell subset with poor antitumor functional capacity. The understanding of these two NK cell states within human tumors may have important implications for prognostication and therapy, and the results from this study provide an important resource for the future investigation of NK cells in human cancer.

Results
Human HNSCC Contains Heterogeneous ILC Subsets and States. Using scRNA-seq, we took an unbiased approach to evaluate the heterogeneity of ILC states within the tumor microenvironment of human HNSCC. We obtained histologically confirmed HNSCC tumor samples from patients undergoing tumor resection, as well as lymph node tissue and peripheral blood (Table 1). These HNSCC tumors are all of epithelial origin, and the tumor sites included the mucosa of the oral cavity or oropharynx. Some were associated with the human papilloma virus (HPV), as determined by either HPV in situ hybridization or by p16 staining, which is a surrogate marker of HPV infection in these tumors. There was no correlation of CD49a or CD103 expression with tumor site or HPV status. We also obtained peripheral blood from healthy volunteers. The tissue and blood were processed into single-cell suspensions (Fig. 1A) to sort comprehensively all ILCs, including NK cells, for scRNA-seq library preparation and sequencing. Human peripheral NK cells (also referred as conventional NK cells) are phenotypically defined as CD56 + lymphocytes lacking expression of CD3 and CD19. Although ILC1s, ILC2s, and ILC3s are characterized by expression of CD127 (IL-7Rγ), some subsets express CD56 (18). Likewise, some NK cells, particularly CD56 bright NK cells, express CD127 (18). Thus, to include all intratumoral ILCs and NK cells, we sorted lineagenegative lymphocytes expressing CD56 and/or CD127 for our analysis (SI Appendix, Fig. 1A).
Dimension reduction of the scRNA-seq data illustrated eight distinct clusters of ILCs (Fig. 1B). Cell indexing of the fluorescence-activated cell sorter (FACS)-sorted cells revealed that the CD127 + cells were primarily within clusters 5 and 6, whereas the CD56 + CD127cells were distributed among all the clusters (SI Appendix, Fig. 1B and Fig. 1B). Further, the cells sorted from the peripheral blood clustered separately (cluster 1, peripheral NK cluster) from cells sorted from the tumor and normal lymph node tissue (SI Appendix, Fig. 1C). In order to understand the different intratumoral clusters, we estimated the fractional compositions of the scRNA-seq profiles by CIBER-SORTx (35), using previously published gene signature matrices derived from bulk RNA-seq profiles (36,37) of sorted human ILC subsets ( Fig. 1 C and D). Using this approach, we were able to objectively assign identities to the seven different intratumoral clusters in an unbiased manner. We identified two NK cell clusters (clusters 2 and 3, designated NK-1 and NK-2), one ILC1 cluster (cluster 5), one ILC3 cluster (cluster 6), two ieILC1 clusters (clusters 7 and 8), and one cluster that appeared to represent a transition state between NK cell and ILC1 (cluster 4, designated NK-ILC1 intermediate) ( Fig. 1 B-D). While cells from all the clusters were present in all the patient samples, there was a significant heterogeneity in the relative proportion of each population between patient samples ( Fig. 1E and SI Appendix, Fig. 1D).
The presence of ILC subsets within primary HNSCC tumors was further confirmed by flow cytometry (Fig. 1F and SI Appendix, Fig. 1 E and F). Within the CD45 + Lin -CD127 + CD56compartment, we found CRTH2 -CD117cells consistent with ILC1s, and CRTH2 -CD117 + cells consistent with ILC3s. The latter were EOMES -T-BET -RORγt + (Fig. 1G). Similar to our scRNA-seq data, there was no evidence of the presence of CD127 + CRTH2 + ILC2s in the tumors. Within the NK cell gate (CD45 + Lin -CD56 + ), we observed cells expressing CD49a and CD103 that were EOMES hi T-BET int RORγt -, consistent with ieILC1s ( Fig. 1 F and G and SI Appendix, Fig. 1 E and F). Some of these cells expressed the proliferation marker Ki-67 ( Fig. 1 F and SI Appendix, Fig. 1G), suggesting that a subset of ieILC1s was actively cycling. Also, intratumoral CD56 + CD49a + ieILC1 cells expressed high levels of CD69, suggesting an activated and/ or tissue-resident state (SI Appendix, Fig. 1H). Similar to the scRNA-sequencing analysis ( Fig. 1E and SI Appendix, Fig. 1D), we observed significant heterogeneity in the proportion of ieILC1 cells between different patient samples by flow cytometry (Fig. 1H and SI Appendix, Fig. 1 I and J).

Intratumoral ILC Populations Display Distinctive Functional and
Effector States. In order to study the function and possible relationship between populations, we first analyzed the gene expression profiles among the different populations (Dataset S1). Peripheral NK cells were characterized by expression of FGFBP2, CXC3R1, and FCGR3A (encoding CD16) ( Fig. 2A). The NK-1 population had a signature of type I interferon (IFN) activation, whereas the NK-2 population had a signature notable for expression of the NR4A2 orphan nuclear receptor transcription factor gene (also known as NURR1) (Fig. 2 A and B), which was recently shown to negatively regulate T cell responsiveness (38)(39)(40). Both the ILC1 and ILC3 clusters expressed IL7R (CD127) and SELL (CD62L) (Fig. 2 A and B). The ILC3 cluster expressed the canonical genes RORC, KIT, IL23R, and AHR ( Fig. 2 A and B and SI Appendix, Fig. 2A). We observed an intermediate NK-ILC1 population with genes that were also expressed in the NK-2 and ILC1 clusters ( Fig. 2A), including AREG (SI Appendix, Fig. 2A), which has been observed to be up-regulated in in CD56 bright NK cells (41). CIBERSORTx fractional compositions showed that the gene expression profiles of cells in the NK-1, ILC1, and NK-ILC1 intermediate clusters were more similar to CD56 bright NK cells when compared to CD56 dim NK cells (SI Appendix, Fig. 2B).
Consistent with the previously proposed hypothesis that ieILC1s are related to NK cells (18), we found that all ieILC1, NK-1, NK-2, and peripheral NK cell clusters expressed NK-related genes GZMA (granzyme A), GZMB (granzyme B), GZMH (granzyme H), PRF1 (perforin), and GNLY (granulysin) (Fig. 2). However, the log fold change plots and the feature plots showed that the ieILC1s had the highest levels of GZMA, GZMB, PRF1, and GNLY ( Fig. 2 B and C and SI Appendix, Fig. 2C), suggesting that these cells may have high cytolytic capacity. In contrast, the NK-2 cluster had relatively lower expression of effector molecule genes (GZMA, GZMB, and PRF1) and was notable for high expression of NR4A1 and NR4A2 (Fig. 2). Thus, among the intratumoral NK cell-related subsets, the ieILC1 subset appeared to have the greatest cytolytic capacity, whereas the NK-2 subset expressed genes suggesting a dysfunctional phenotype.

NK Cells Have Two Distinct Intratumoral Differentiation Trajectories,
Leading to Either an ieILC1 Phenotype or an NK-2 Phenotype. The tumor microenvironment is known to affect the gene expression of immune cells, including their effector function phenotypes (42). In order to understand the relationship between the peripheral cNK cells and the intratumoral NK cell subsets, we used pseudotime ordering (43) of cells that had primarily NK cell signatures identified by CIBERSORTx deconvolution of the subset clusters ( Fig. 1 C and D). This enabled the inference of potential trajectories of NK cell differentiation among the different cell states seen in our scRNA-seq data. Relative to the cNK cells, two divergent trajectories, passing through an intermediate NK-1 population, were observed: one resulting in the ieILC1 state and another resulting in the NR4A2-expressing NK-2 state ( Fig. 3 and SI Appendix, Fig. 3).
Peripheral NK Cells Can Differentiate into CD49a + CD103 + Cells and NR4A-Expressing NK Cells In Vitro. The scRNA-seq data prompted us to investigate the plasticity of cNK cells to differentiate into ieILC1-like cells or NK-2−like cells in vitro. Our analysis of The Cancer Genome Atlas (TCGA) database showed that the expression of IL-15 in HNSCC tumors was highly correlated with top genes expressed by the ieILC1 cluster ( Fig. 4A) but not with the NK-2 gene signature (SI Appendix, Fig. 4A). Given that IL-15 is an important cytokine for NK cells and that it is involved in the generation of CD49a hi CD103 hi ILC1-like cells in the murine MMTV-PyMT mammary tumor model (31), we hypothesized that IL-15 may be involved in the differentiation of NK cells into ieILC1-like cells in human tumors. In vitro coculture of peripheral cNK cells, which lack expression of CD49a and CD103 (SI Appendix, Fig. 4B), with HNSCC cells resulted in the upregulation of CD49a and CD103 on a subset of NK cells, and this up-regulation was highly dependent on the presence of IL-15 ( Fig. 4 B and C). In addition to HNSCC cells, other tumor cell lines also induced an ieILC1-like phenotype with variable efficiencies (SI Appendix, Fig. 4C). Similar to our observation of primary intratumoral CD49a + ieILC1 cells (SI Appendix, Fig. 1 E, F, and J), CD103 tended to be coexpressed with CD49a following in vitro differentiation ( Fig. 4C and SI Appendix, Fig. 4C). Phenotypically, in vitro differentiated CD56 + CD49a + CD103 + cells were NKp44 + NKp80 -CD94 + CD16 -GzmB + and expressed T-BET and EOMES (Fig. 4D). The phenotype of CD49a + CD103 + cells was similar to that of CD49a + CD103cells, suggesting that, overall, in vitro differentiated CD49a + cells were phenotypically similar to primary intratumoral ieILC1s (SI Appendix, Figs. 1 F and G and 2 A and B). In contrast, the phenotype of CD49a -CD103 + cells was similar to that of CD49a -CD103cells (Fig. 4D), and CD49acells were characterized by high expression of NR4A genes, particularly NR4A2 (Fig. 4E), similar to the NK-2 cells we identified in primary tumors ( Fig. 2 A and B). Importantly, the differentiation of NK cells into CD49a + CD103 + cells required contact with the tumor cells (Fig. 4 F and G), suggesting the involvement of unknown receptor−ligand interactions.  Immature NK Cells Preferentially Differentiate into CD49a + CD103 + Cells. Peripheral NK cells are divided into two main subsets: CD56 dim and CD56 bright NK cells, which differ in their function and maturation status. CD56 bright NK cells have an immature phenotype compared to CD56 dim NK cells, and represent a minor subset in peripheral blood but are highly enriched in tissues, including malignant tissues (44). We assessed whether there was a difference in the plasticity of peripheral CD56 bright and CD56 dim NK cells to differentiate into the ieILC1-like cells. We observed that CD56 bright NK cells up-regulated CD49a to a greater extent than CD56 dim NK cells when cultured with HNSCC cells and IL-15 ( Fig. 5 A and B), and the percentage of CD49a + CD103 + cells was higher following culture with CD56 bright NK cells compared to CD56 dim NK cells (Fig. 5C).
NK cells develop through a multistage process in the BM, secondary lymphoid tissues, and peripheral blood (45). In the peripheral blood, there are three distinct latter stages of NK cell development that are observed: stage 4b (CD94 + NKp80 + CD16 − CD57 − ), stage 5 (CD94 +/− NKp80 + CD16 + CD57 − ), and stage 6 (CD94 − NKp80 + CD16 + CD57 + ), which represents fully differentiated cells. While most CD56 bright NK cells are stage 4b, some are stage 5. Therefore, in order to study more accurately the degree of plasticity of peripheral NK cells and their ability to take on characteristics of ieILC1s, we assessed differences between stage 4b and stage 5 NK cells. We sorted stage 4b cells and two distinct subsets of stage 5 cells, with the latter being distinguished by the expression of CD94 (CD94 + stage 5 cells being less mature than CD94stage 5 cells) and cultured them with HNSCC and IL-15 (Fig. 5D). We found that the capacity of stage 4b cells to differentiate into CD49a + CD103 + was much greater than that of stage 5 CD94 + , which, in turn, was much greater than the capacity of stage 5 CD94 − cells. Overall, this indicates that circulating immature NK cells preferentially differentiate into ieILC1-like cells, compared to mature NK cells.
Since CD103 has been shown to be a TGF-β−imprinted gene (46) and because TGF-β is involved in NK-ILC1 differentiation (32,33,47), we assessed whether TGF-β played a role in the induction of the ieILC1 phenotype in vitro. In support of this, we found that TGFβ could be detected in cultured supernatants of HNSCC tumor cells and that inhibition of TGF-β signaling reduced the proportion of ieILC1-like cells that were seen after NK cells are cocultured with the tumor cells (SI Appendix, Fig.  5). Although culture of NK cells with IL-15 and soluble TGF-β (in the absence of tumor cells) was able to induce the expression of CD49a and CD103 (SI Appendix, Fig. 5C), in accordance with previous observations (48), the appearance of ieILC1-like cells, when NK cells are cocultured with HNSCC cells, required contact between the cells (Fig. 4F). Membrane-bound TGF-β could not be detected on the PCI-13 cells; thus, another unknown mechanism related to TGF-β signaling appears to be involved.
Conventional NK Cells Differentiate into CD49a + CD103 + Cells in the Tumor Microenvironment and Have Potent Antitumor Activity. Considering the in vitro differentiation of CD49a + CD103 + cells (Fig. 4 B and C), we hypothesized that IL-15 would promote the differentiation of cNK cells into ieILC1-like cells in vivo. Human HNSCC tumor cells were used to establish tumors in the subcutaneous compartment of immunodeficient NOD-scid IL2Rgamma null (NSG) mice, which lack T cells, B cells, and NK cells. Peripheral CD56 + NK cells (SI Appendix, Fig. 6A) or CD56 bright NK cells (Fig. 7A), sorted from healthy volunteers, were directly injected into the tumors, and the mice were treated with or without IL-15. After 6 d, tumors were harvested, and we assessed the intratumoral NK cells by flow cytometry. We observed that a subset of NK cells had up-regulated the ieILC1 markers CD49a and CD103 and that this was dependent on exogenous administration of IL-15 to the mice (Fig. 7A and SI Appendix, Fig. 6A). Similar to our in vitro findings (Fig. 6 A-C), we observed that, in the presence of IL-15, the in vivo-differentiated CD49a + cells produced higher amounts of IFNγ compared to the CD49acells following stimulation (SI Appendix, Fig. 6B).
Because CD49a + ieILC1-like cells produce higher levels of IFNγ, TNFα, and granzymes than CD49acells and are able to degranulate more readily upon stimulation (Fig. 6 A-C), we hypothesized that CD49a + cells would be able to control tumor growth in vivo better than CD49a -NK cells. To assess this, in vitro-differentiated CD49a + or CD49a -NK cells were FACS sorted and mixed separately with HNSCC tumor cells. Equal numbers of total cells of these separate mixes were injected into the subcutaneous compartment of NSG mice (Fig. 7B). We observed that CD49a + ieILC1-like cells were able to control tumor growth to a significantly greater degree than the CD49a -NK cells derived from the same in vitro culture. Together, our data indicate that the tumor microenvironment can shape the plasticity of NK cells by polarizing them either into an ieILC1-like phenotype or an NK-2−like phenotype, which may have opposing effects on the tumor and clinical behavior (see SI Appendix, Fig. 7).

Discussion
In this study, we comprehensively assessed the panorama of intratumoral ILCs, including NK cells, in human HNSCC using an unbiased scRNA-seq approach and found seven cell states, distinct from peripheral blood NK cells (Fig. 1B). Previous tumor scRNA-seq studies included NK cells among many different cell types (other immune cells, tumor cells, and stromal cells); however, because the NK cell/ILC compartment in tumors is relatively small (∼1 to 5% of the immune cells), insights into intratumoral NK cell/ILC biology have been limited. A recent scRNA-seq study in melanoma included only CD56 + NK cells, limiting the spectrum of ILC subsets, and relied heavily on integration of gene expression profiles between samples in their analysis (49). Here, despite having a larger number of source patients and a greater diversity of samples, we did not require any batch correction or rely on any integrative approaches in our analysis of the entire ILC population.  Our study revealed that, within the tumor microenvironment, peripheral cNK cells transition through an intermediate NK-1 state and can differentiate into two distinct and divergent terminal states. One is a CD49a + CD103 + ieILC1-like subset with high granzyme and perforin expression and high capacity to control tumor growth, and the other is an NR4A2-expressing CD49asubset with low capacity to control tumor growth.
Mechanistically, we found that the differentiation of cNK cells into ieILC1-like cells required IL-15 and was dependent on contact between NK cells and HNSCC cells. Also, we found that immature NK cells (i.e., CD56 bright stage 4b cells) have the greatest potential to differentiate into ieILC1-like cells among the circulating NK cells. In accordance with previous studies showing that TGF-β is involved in the imprinting of CD103, and the differentiation of cNK into ILC1 cells (32,33,48), we found that TGF-β signaling in the NK cells plays a role in the acquisition of the ieILC1-like phenotype. This observation is somewhat surprising, given that the ieILC1-like cells obtained under our conditions display much high effector functions and that TGF-β has been shown to have suppressive effects on NK cell function (50). It is possible that the combination of TGF-β with IL-15 has quite different effects on NK cells than TGF-β alone, similar to the role of TGF-β in the differentiation of both Treg cells and Th17 cells (51).
An ieILC1-like NK cell subset has been observed by mass cytometry studies of human lung and colorectal tumors previously; however, in that study, the ieILC1-like cells had lower expression of granzyme B and perforin compared to a CD103 -NK cell population (52). It is possible that the CD103cells in the lung are predominantly peripheral blood NK cells, given the highly vascularized nature of the lung tissue. It is also possible that different tumor types (e.g., squamous cell carcinoma vs. adenocarcinoma) induce different ieILC1 functional states. Interestingly, the CD49a + CD103 + ieILC1-like population that we observed in HNSCC resembled the murine CD49a hi CD103 hi GzmB hi ILC1-like tissue-resident lymphocyte population observed in the tumor microenvironment of MMTV-PyMT tumors (31). Similar to the human ieILC1-like cells that we identified here, the murine ILC1-like population also exhibited potent cytotoxicity and capacity to control tumor growth. Among the two ieILC1 clusters that we observed, one had a distinctive cycling phenotype, which may be reflective of greater IL-15 signaling or other stimuli present within the tumor microenvironment.
NR4A2 was recently demonstrated to confer a hyporesponsive phenotype in cytotoxic T cells and to promote Treg differentiation in CD4 + T cells (38-40, 53, 54). Our finding that NR4A2 + CD49a -NK cells show poor tumor control in vivo is consistent with the inhibitory function of NR4A2 on T cells. While the underlying mechanism remains to be determined, it is interesting to note that NR4A2 + CD49a -NK cells have significantly lower IFNγ and TNFα expression. In future studies, it will be important to understand the role of NR4A2 in NK cell function and the tumorassociated factors governing the differentiation of peripheral NK cells into divergent NK cell states.

Materials and Methods
Study Design. We used a total of eight primary tumor samples and four paired nodal metastases from HNSCC patients, as well as peripheral blood mononuclear cells (PBMC) from four HNSCC patients and two healthy volunteers for scRNA-seq analysis. In vitro experiments used a sample size of at least three different donors per experiment. In vitro and in vivo experiments were replicated at least three times. No data were excluded from analysis.
The goal of the study was to examine the spectrum of ILCs present within the tumor microenvironment of human HNSCC. For this, primary tissue samples were obtained from patients consented at Stanford University. Our prespecified hypothesis was that all known ILCs would be observed; however,  Harvested xenograft tumors growing in NSG mice were minced into small fragments, which were mechanically dissociated using a gentleMACS Dissociator (Miltenyi Biotec) followed by enzymatic digestion using Collagenase/ Hyaluronidase (Stemcell Technologies; Cat# 07912) for 30 min at RT with constant rotation. Then, the mixture was filtered with a 70-μm cell strainer (Falcon; Cat# 352350) to obtain single-cell suspension. Mice and In Vivo NK Cell Differentiation. NSG mice were housed and bred in an animal facility under standard temperature, humidity, and timed lighting conditions and were provided with mouse chow and water ad libitum. All animal experimentation protocols were approved by the Stanford University Animal Care and Use Committee. Mice of ages 7 wk to 10 wk were used in our experiments. One million PCI-13 cells were injected subcutaneously in the flank of the mouse, two injections per mouse. After 10 d to 15 d, when the tumors reached a group average of 200 mm 3 to 300 mm 3 , mice were sublethally irradiated 300cG using a Faxitron irradiator (Faxitron). Then, 2 × 10 6 to 5 × 10 6 sorted NK cells were injected intratumorally, per tumor. Where indicated, human recombinant IL-15 (BioLegend; Cat# 570306) was injected intraperitoneally on the day of NK cell intratumoral injection and every 2 d thereafter. After 6 d, tumors were harvested and processed (as previously described), and the function and phenotype of the intratumor NK cells was studied.
Functional Characterization of NK Cells. NK cells were stimulated in vitro for 4 h with either PCI-13 (E:T = 1:1); a mixture of cytokines containing IL-12 (10 ng/mL), IL15 (10 ng/mL), and IL-18 (50 ng/mL); PMA (at 50 ng/mL) and ionomycin (500 ng/mL); or control condition. At 20 min into the incubation, protein transport inhibitors (BD; Cat# 555029 and #554724) and anti-CD107 antibody were added to the culture. After the stimulation, NK cells were harvested, surface stained, and fixed, and intracellular content on IFNγ was determined by intracellular staining.
In vivo antitumor potential of NK cells was determined by mixing PCI-13 with sorted NK cells (NK cells:PCI-13 cells ratio of 1:2.5), and the mixture was injected subcutaneously in the flanks of NSG mice; tumor growth was measured over time.
To compare the killing capacity of NK cells, xCELLigence Real-time Cellular Analysis (Agilent Technologies) was performed. On day 0, the xCELLigence microtiter plate was loaded with 10,000 target cells; on day 1, 10,000 sorted effector cells (at E:T ratio of 1:1) were loaded into wells in the presence of IL-15 (10 ng/mL). The cytolytic activity was measured based on the electric impedance (cell index) between the biosensors, and the obtained cell index values were normalized.
Cytokine Secretion. Sorted NK cells were stimulated with PMA (50 ng/mL) and ionomycin (500 ng/mL) for 12 h, in 96-well round-bottom plates, at 50,000 NK cells per well in duplicate. Then, supernatants were harvested, frozen at −80°C, and analyzed by 62-plex Luminex Immuno-Assay by the Human Immune Monitoring Center facility at Stanford University.
PCI-13 cells were cultured for 3 d, in 96-well flat-bottom plates, at 10 4 cells per well, in duplicate. Then, culture supernatants were harvested, frozen at −80°C, and analyzed by 62-plex Luminex Immuno-Assay by the Human Immune Monitoring Center facility at Stanford University.
Antibody Staining for Phenotypic Analysis and Cell Sorting. The list of antibodies used to sort and phenotype NK cells and ILCs is available in SI Appendix, Then, the cells were resuspended in either FACS buffer for single-cell sorting for scRNA-seq or NK cell culture media for bulk sorting for in vivo and in vitro experiments, or were fixed and permeabilized using the Foxp3/Transcription Factor Staining Buffer Set (eBioscience; Cat# 00-5523-00) followed by intracellular staining according to manufacturer's protocol for FACS analysis.
qRT-PCR. Sorted cells were frozen on TriZol Reagent (Life Technologies; Cat# 15596018). RNA was extracted with the RNeasy mini kit (QIAGEN; Cat# 74104), and complementary DNA (cDNA) was made using the Maxima First Strand cDNA Synthesis kit (Thermo Fisher Scientific; Cat# K1642). Quantitative gene expression was performed using the Taqman Gene Expression Assay with the recommended best coverage primers for glyceraldehyde-3-phosphate dehydrogenase (GAPDH), NR4A1, NR4A2, and NR4A3 (Thermo Fisher Scientific). Each gene expression assessment was measured in triplicate. Gene expression was determined relative to GAPDH messenger RNA (mRNA) expression using the formula 2 (−(NR4A-GAPDH)) . In some cases, NR4A gene expression was normalized to that of CD49a-NK cells (which was set to one), to allow easy comparison with the mRNA expression on other NK cell population.
Preparation of Single-Cell Libraries and Sequencing. Sorted single cells were processed for scRNA-seq library preparation using the previously described Smart-seq2 protocol (55). Briefly, the RNA from each cell was reverse transcribed and PCR amplified (25 cycles). Amplified cDNA product was purified using 0.55× Agencourt AMPure XP Beads (Beckman Coulter; Cat# A63880) on Biomek FX Automated Workstation (Beckman Coulter). A 2-μL aliquot of purified cDNA sample from each well was analyzed for quality on a 5300 Fragment Analyzer with 96 capillary array (Agilent), using High Sensitivity NGS Fragment kit-DNF-474 (Agilent). Plates with a minimum DNA concentration of 0.2 ng/μL per well (average range 0.5 ng/μL to 1.5 ng/μL) were used for the final library preparation. First, the cDNA concentration was normalized with nuclease free water (IDT; Cat# 11-04-02-01) to a concentration less than or equal to 1 ng per well and transferred to a 384-well low volume serial dilution plate (TTP Labtech; Cat# 4150-05828) using the Mosquito ×1 (TTP Labtech). Then, library preparation steps as per Smart-seq2 protocol were performed using Nextera XT DNA Library Prep Kit (Illumina; FC-131-1096) in a 384-Well Rigi-Plate PCR Microplate (TTP Labtech; Cat# PCR-384-RGD-C) using Mosquito HTS (TTP Labtech). One-microliter volumes of indexed libraries from each well of the 384-well plate were pooled together and purified using 0.8× Agencourt AMPure XP Beads (Beckman Coulter; Cat# A63880). Concentration of the purified pooled library was determined using the Bioanalyzer (Agilent). The pooled library was run on the HiSeq 2500 high-output v4 sequencer (Illumina), with 2 × 125 bp read length. One 384-cell pooled library was run on one lane of the HiSeq 2500, with an average output of 1 × 10 6 reads per cell.
Library Size Normalization, Dimensional Reduction, and Visualization. Further downstream processing was performed using the R package Seurat (59), starting from the gene expression matrix comprising raw counts. Specific parameters included the following: 1) LogNormalization was performed using a scale factor of 10,000. 2) Variable features were identified using the "vst" selection method.
3) The statistical significance of principal components was determined by the JackStraw function. 4) Both JackStraw and Elbow plots were used to determine the significant principal components for further evaluation.
An initial round of clustering was performed using the top 12 principal components at a resolution of 0.9. Visualization was performed using the Uniform Manifold Approximation and Projection (UMAP) dimensional reduction technique. Ten clusters were identified, including three clusters comprising CD3 + /CD4 + T cells and CD86 + antigen-presenting cells (SI Appendix, Fig. 10). Cells from these three contaminating clusters were removed. The remaining cells, comprising NK cells and ILCs, were reanalyzed starting from the gene expression matrix comprising raw counts. Clustering was performed using the top 12 principal components at a resolution of 0.8, resulting in eight distinct UMAP clusters. Each cell was annotated by its tissue of origin, patient source, and CD56 and CD127 protein expression based on flow cytometry data.
Gene expression markers (differentially expressed genes) representative of each cluster were identified by the Seurat "FindAllMarkers" function using the Wilcoxon Rank Sum test and fulfilling the following thresholds: 1) minimum detection fraction of 0.25 and 2) log fold change of at least 0.25. The top 20 differentially expressed genes for each cluster were represented in a heat map. For fold change comparisons between two clusters, the "FindMarkers" function was used. For the log fold change plots, genes with a fold change of three or greater and an adjusted P value of < 0.05 were considered as significant.
Deconvolution of Gene Expression Profiles. In order to assign phenotypic identities to each single-cell cluster, CIBERSORTx, a gene signature deconvolution approach, was used (35,60).
RNA-seq gene expression profiles from Yudanin et al. (36) comprising raw counts of bulk-sorted human NK, ILC1, ILC2, and ILC3 cells from the ileum, jejunum, and spleen were downloaded from the Gene Expression Omnibus (GSE126107). Reads per kilobase million (RPKM) values were obtained using the "rpkm" function of the edgeR package (61).
The log2-normalized RPKM gene expression profiles of bulk-sorted peripheral human NK cells (CD56 bright, CD56dim/CD57pos, CD56dim/ CD57neg), as well as ieILCs and ILC3s from human tonsil (GSE112813), were obtained from the supplemental information accompanying Collins et al. (37) The log2-normalized RPKM gene expression matrix was converted to RPKM for consistency.
A reference signature matrix was created for each of the two datasets using the CIBERSORTx web interface (https://cibersortx.stanford.edu), using the default parameters (kappa = 999, q value = 0.01, number of genes 300 to 500). The fractional composition of each single cell was imputed using this signature matrix, with "B-mode" batch correction enabled to account for the different gene expression profiling techniques, and quantile normalization disabled as recommended by the authors.
Pseudotime Analysis. Pseudotime analysis was performed using the R package Monocle version 3.0 (43). The log-normalized gene expression matrix for the top 50 differentially expressed genes of each cluster (total of 400 genes) was used as input. Dimension reduction was performed using the reduc-eDimension function and the DDRtree reduction method with a maximum of two components. Pseudotime was determined by the orderCells function selecting the peripheral cluster as the root state.
The expression of marker genes along the ieILC1 and NK-2 trajectories was plotted using the plot_genes_branched_pseudotime function, with the curves representing Loess-smoothed expression of marker genes along each trajectory.
Correlation of IL-15 Gene Expression and ieILC1 Signature. The expression values of IL15 were correlated with the values of an ieILC1 signature, across head and neck (HNSC) primary tumor samples from TCGA (n = 522). The transcripts-per-million (TPM) counts of IL15 and each of the top 20 ieILC1 marker genes used for generating the ieILC1 signature were log-transformed using formula log 2 (TPM + 1), and scaled to mean zero and unit variance across samples. The values of the ieILC1 signature were obtained by averaging the log-transformed and scaled values of the top 20 differentially expressed genes in the ieILC1 cluster.
Statistical Analysis for In Vitro and In Vivo Experiments. Results are expressed as mean ± SE. For statistical analysis, matched-pairs t test, two-tailed with a 95% confidence level, was performed using Prism software version 8.3.1 (GraphPad Software, LLC). P values are shown when statistically difference was reached.
Data Availability. The mapped BAM files from the blood and tumor specimens of the HNSCC patients have been uploaded to the database of Genotypes and Phenotypes (dbGaP Study Accession: phs002002.v1.p1, "Landscape of Intratumoral NK Cell and ILC in Head and Neck Squamous Cell Carcinoma", Study ID: 38031).