High-dimensional single cell analysis identifies stem-like cytotoxic CD8⁺ T cells infiltrating human tumors

Despite the abundance of tumor-infiltrating lymphocytes (TILs) that has been associated with a better prognosis in several human cancers, the anti-tumor response mediated by immune cells is insufficient to induce tumor regression, mainly due to the potent immunosuppression found in the tumor microenvironment (TME; Croci et al., 2007; Gajewski et al., 2013). Among TILs, the number, localization, and quality of cytotoxic CD8⁺ T cells have been shown to play a major role in this regard (Galon et al., 2013). Specifically, increased infiltration of CD8⁺ T cells displaying an effector phenotype and expressing killer molecules such as granzyme B at the tumor invasive margin has been associated with favorable prognosis in colorectal (Pagès et al., 2005; Galon et al., 2006; Bindea et al., 2013) and other types of cancer (Schumacher et al., 2001; Hamanishi et al., 2007; Ganesan et al., 2017). Immunosuppression within the TME comprises several mechanisms that include, but are not limited to, the presence of suppressive populations such as tumor-associated macrophages/myeloid cells and CD4⁺ regulatory T cells, chronic antigenic stimulation, and inhibitory metabolites, cytokines, and ligands. All together, they generate a state of dysfunction in T cells also known as T cell exhaustion, characterized by poor proliferative capacity, diminished cytokine production and killing function, and increased expression of several inhibitory receptors on the cell surface (Wherry, 2011). Current immunotherapeutic strategies targeting inhibitory receptors such as Programmed Death-1 (PD-1) and Cytotoxic T Lymphocyte Antigen-4 (CTLA-4) by blocking antibodies reinvigorate T cells and further boost T cell infiltration, expansion, and effector functions (Leach et al., 1996; Freeman et al., 2000; Sharma and Allison, 2015; Wu et al., 2016; Wei et al., 2017), thus resulting in long-term disease stabilization or, in some cases, tumor regression. However, responses are only seen in a subset of patients (Hodi et al., 2010; Topalian et al., 2012; Jacquelot et al., 2017; Krieg et al., 2018), hence highlighting the need for improved strategies.

The CD8⁺ T cell compartment in peripheral blood and tissues is largely diverse, comprising several subsets with different degrees of specialization in phenotype, function, and gene expression (Mahnke et al., 2013; Farber et al., 2014). Recent application of high-content single cell technologies at the level of the whole transcriptome such as single cell RNA sequencing (scRNA-seq; Tirosh et al., 2016; Zheng et al., 2017) and cytometry by time-of-flight (CyTOF; Chevrier et al., 2017; Lavin et al., 2017) suggested that also CD8⁺ TILs are functionally heterogeneous, displaying different levels of T cell activation and exhaustion (Tirosh et al., 2016; Chevrier et al., 2017). Differential analysis of these states, which may coexist in the same T cell or be mutually exclusive, led to the identification of new molecular regulators of exhaustion, such as the transcription factor GATA-3, which inhibits T cell effector function (Singer et al., 2017). Along the same lines, recent data obtained from chronically infected mice showed that the exhausted PD-1⁺ T cell compartment is organized in a hierarchy of differentiation, comprising less-differentiated cells with intermediate levels of PD-1 (PD-1^int) that retain self-renewal capacity and differentiation potential upon PD-1 blockade and PD-1 high (PD-1^hi) cells that are terminally differentiated and functionally inefficient (He et al., 2016; Im et al., 2016; Leong et al., 2016; Utzschneider et al., 2016). This precursor–progeny relationship is reminiscent of that seen in the blood and lymphoid tissues in physiological situations, where CD8⁺ T memory stem cells (T_SCM) are at the apex of the differentiation program and are currently considered a major reservoir of long-term immunity (Gattinoni et al., 2011; Lugli et al., 2013a; Fuertes Marraco et al., 2015; Oliveira et al., 2015; Akondy et al., 2017). Whether a stepwise differentiation program is also present in the context of T cell exhaustion in human tumors is still not clear.

In this study, 27-parameter flow cytometry applied to primary non-small cell lung cancer (NSCLC) samples, the adjacent cancer-free tissue and the blood enabled us to investigate millions of single CD8⁺ T cells and thus to identify rare immunophenotypes that are exclusively present within the tumor. In combination with scRNA-seq, we report the occurrence of stem-like CD8⁺ T cells that lack the canonical T_SCM phenotype of circulating T cells. Rather, they simultaneously display a partial exhausted phenotype along with potent effector potential and enhanced self-renewal capacity. We further show that these CXCR5⁺ TILs are lost with disease progression.

To better characterize the diversity occurring at the level of the CD8⁺ T cell compartment specifically within the tumor, we developed a polychromatic flow cytometry panel capable of simultaneously investigating 27 parameters on millions of single cells from resected tumor (n = 53), paired adjacent cancer-free lung tissue (hereafter referred to as “normal”; n = 45), and peripheral blood (n = 22) samples of treatment-naive, NSCLC patients (Fig. 1 A and Table 1). A list of the markers included in the panel, providing information on differentiation, activation, proliferation, and exhaustion status of T cells (see Mahnke et al., 2013), is reported in Table S1. Initial data representation with t-distributed stochastic neighboring embedding (tSNE), used to reduce the dimensionality of the 27-parameter dataset (Amir et al., 2013), revealed that single cells from the three districts are substantially diverse in terms of their CD8⁺ T cell profile (Fig. 1 B, top row). To gain more insight on the immunophenotypes that preferentially characterize CD8⁺ T cells from the different specimens, we concatenated CD8⁺ T cell fluorescent data from all samples and analyzed them with Phenograph (Fig. 1 A), an unsupervised clustering algorithm suitable for high-dimensional single cell data (Levine et al., 2015). Phenograph identified 24 different CD8⁺ putative subpopulations (“clusters;” visualized on top of the tSNE plot; Fig. 1 B, bottom row), whose frequencies are subsequently determined in each sample type. Hierarchical metaclustering of these data clearly separated blood and tumor specimens, indicating that they are substantially different at the high-dimensional single cell level (Fig. 1 C). Instead, normal tissues were split into two major groups that were more similar to either the blood or the tumor samples.

To define the phenotypic identity of immunophenotypes resulting from Phenograph clustering, we calculated the integrated median fluorescence intensity (iMFI) values (Darrah et al., 2007) of each marker in each Phenograph cluster and visualized values in a heat map, along with the median frequency of the same clusters in the blood, normal tissues, and tumor samples (Fig. 1 D). Hierarchical metaclustering of Phenograph clusters (rows) and markers (columns) grouped subpopulations with similar immunophenotypes. Two of these clusters (1 and 3) were mostly present within the blood and displayed a phenotypic identity coincident with naive T cells, i.e., characterized by high levels of CD45RA, CCR7, and CD27 and by the lack of multiple memory markers, such as CD95, PD-1, and T-bet. Other clusters, including 24, 14, 5, and 22, were virtually absent from blood while present in comparable abundance in both the normal tissue and tumor samples, indicating that they are lung-specific subsets. Despite displaying a common CD95⁺ CD27^lo T-bet^– memory phenotype, these were mainly characterized by variable levels of CD69 and CXCR3 expression, among others, thereby indicating tissue-residency (Wong et al., 2016). The abundance of these phenotypic clusters was counterbalanced by the absence of CD27^int T-bet^int CD69^– transitional memory cells (cluster 2) that instead are highly frequent in the blood. Importantly, the tumor lacked multiple subsets of CD45RA^+/−CCR7^– terminal effector cells (clusters 18, 23, 4, and 8), characterized by variable levels of CD57 and T-bet expression, that are known to be highly cytolytic, but replicative senescent (Takata and Takiguchi, 2006; Chattopadhyay et al., 2009; Mahnke et al., 2013), leading us to conclude that the tumor is largely devoid of anti-tumor effectors. On the contrary, multiple clusters were specifically enriched in the tumor samples. These included both exhausted CD8⁺ T cell clusters 21, 13, and, to a lesser extent, 20, characterized by high expression of activation markers HLA-DR, CD25, and IRF4, and different combinations in expression of the inhibitory receptors PD-1, TIM-3, and TIGIT. Activated CD8⁺ T cell cluster 7, expressing only high levels of the activation markers, but not the proliferation marker Ki-67, was also enriched in lung tumors. Remarkably, Phenograph identified two memory clusters expressing CXCR5, i.e., 15 and 17, the latter featuring the highest expression. In addition, only the frequency of cluster 17, characterized by a CD69⁺ PD-1^int CD45RA^– CCR7^– effector-like phenotype, but lacking additional activation and exhaustion markers, was most different between the tumor and the adjacent lung tissue or the blood, suggesting this is a tumor-specific subpopulation (Fig. 1 D and Fig. S1).

It has been recently reported that chronically infected mice harbor, besides CXCR5^– TIM-3⁺ CD8⁺ T cells that are exhausted and thus incapable of persisting in the long term, CXCR5⁺ TIM-3^– CD8⁺ T cell progenitors that are PD-1^int and, upon PD-1 blockade, are responsible for controlling chronic viral infections (He et al., 2016; Im et al., 2016; Leong et al., 2016; Utzschneider et al., 2016; Wu et al., 2016). Therefore, these cells share phenotypic features with human CXCR5⁺ CD8⁺ TILs identified in Phenograph clusters 15 and 17. We initially confirmed that results obtained with Phenograph clustering could be obtained by standard approaches and indeed found that CXCR5⁺ TIM-3^– CD8⁺ T cells isolated by manual gating are enriched within tumors (Fig. 2, A and B), indicating that unsupervised clustering of multidimensional single cell data are a reliable tool for the identification of phenotypically different subpopulations. In addition to the CXCR5⁺ TIM-3^– and CXCR5^– TIM-3⁺ CD8⁺ T cells, our manual analysis identified a third subset of CD8⁺ TILs lacking both CXCR5 and TIM-3 expression that is infrequent in chronically infected mice, indicating that human TILs are characterized by additional phenotypic complexity. As the major difference between clusters 15 and 17 is not only CXCR5 abundance, but also expression of the tissue residency marker CD69, we next sought to determine to which extent CXCR5⁺ TIM-3^– CD8⁺ T cells in tumors express CD69. The vast majority of these cells were CD69⁺ compared with those from normal tissue and blood samples (Figs. 2, C–E). Instead, CD69^– CXCR5⁺ TIM-3^– CD8⁺ T cells were relatively infrequent in both normal lung tissues and tumors, and their abundance decreased only slightly compared with blood (Fig. 2 F), thereby indicating that Phenograph overestimates the frequency of this phenotype. As a whole, lung cancer-infiltrating CXCR5⁺ CD8⁺ T cells were largely different compared with those from blood and the normal tissue at the phenotypic level, as they displayed a uniform CD45RA^–CCR7^–, effector memory-like phenotype, further characterized by increased expression of CD27 and PD-1, high levels of CXCR3, and lacking CD57 and T-bet (Fig. S1). In addition, CXCR5⁺ TIM-3^– CD8⁺ T cells could be detected in primary colorectal carcinomas (1° colorectal carcinoma [CRC]; Fig. 2 G) and colorectal liver metastases (Fig. 2 H) at comparable frequencies, and, similarly to NSCLC-infiltrating TILs, they were for the vast majority CD69⁺.

We next determined how the CXCR5⁺ CD8⁺ TILs relate to other subpopulations within lung tumors. Pearson correlation analysis between the frequency of the CD69⁺ CXCR5⁺ cluster 17 and the additional clusters identified by Phenograph revealed a significant negative correlation with cluster 7 (HLA-DR⁺ CD98^hi activated cells), cluster 9 (Ki-67⁺ proliferating cells), and cluster 13 (exhausted cells; Fig. 2 I). Confirming what was previously revealed by bioinformatic analysis, Ki-67 was poorly expressed by CXCR5⁺ CD8⁺ T cells as determined by manual gating, while it increased progressively in CXCR5^– TIM-3^– and CXCR5^– TIM-3⁺ CD8⁺ TILs (Fig. 2, J and K). As a whole, tumor-infiltrating CXCR5⁺ CD8⁺ T cells are a quiescent population with hallmarks of tissue residency, but not of exhaustion and activation.

Due to their enhanced self-renewal capacity and multipotency in the context of chronic infection, CXCR5⁺ CD8⁺ T cells resemble stem-like memory T cells (He et al., 2016; Im et al., 2016; Leong et al., 2016; Utzschneider et al., 2016). The presence of a population with stem-like characteristics in human tumors accessible to PD-1 or, more in general, checkpoint blockade, would be extremely beneficial for the local anti-tumor immune response. Canonical human T_SCM cells found in lymphoid tissues and circulation bear a CD45RA⁺ CCR7⁺ CD27⁺ CD95⁺ phenotype (Zhang et al., 2005; Gattinoni et al., 2011; Lugli et al., 2013a) that is largely different from that of the identified CXCR5⁺ CD8⁺ TILs (Fig. 1 D). These T_SCM cells were virtually absent in NSCLC TILs and in lung tissues (Fig. S2), indicating that the lung does not support their development or infiltration and subsequent maintenance. As mouse CXCR5⁺ CD8⁺ T cells highly respond to PD-1 blockade, we sought to determine whether human tumor-infiltrating CXCR5⁺ CD8⁺ T cells have shared features. The gene expression profile of CXCR5⁺ CD8⁺ T cells from chronically infected mice is similar to that of CD4⁺ T follicular helper (Tfh) cells and of CD8⁺ T cell memory precursors (He et al., 2016; Im et al., 2016; Leong et al., 2016; Utzschneider et al., 2016), suggesting they interact with B cells in lymphoid tissues and retain a long-lived memory potential, respectively. Despite different tumor types having largely different biology, infiltrating CD8⁺ T cells share overlapping phenotypes, leading us to reason that a publicly available scRNA-seq dataset of CD45⁺ leukocytes isolated from melanoma patients could better define the identity of CXCR5⁺ CD8⁺ T cells at the global transcriptomic level (Tirosh et al., 2016). In this dataset, a total of 912 cells were labeled as CD8⁺, of which 58 expressed CXCR5, thus comprising 6.4% of the melanoma-infiltrating CD8⁺ T cells, a frequency similar to that of NSCLC-infiltrating CXCR5⁺ CD8⁺ T cells. Initial exploratory analysis by Pearson correlation between the expression of CXCR5 and of all the genes in the matrix revealed that transcripts typically present in long-lived memory T cells such as TCF7 (encoding TCF-1), WNT2, MYB, CD28, CCR7, BCL6, and SELL (encoding CD62L) are correlated with CXCR5 expression, while those present in exhausted or more differentiated effector CD8⁺ T cells such as HAVCR2 (encoding TIM-3), GZMB, PRF1, IFNG, and PRDM1 (encoding Blimp-1) are anti-correlated (Fig. 3 A). To confirm these results, we further subdivided CD8⁺ TILs from this dataset in three subsets based on CXCR5 and HAVCR2 expression levels. In particular, based on the combination of these two markers, we were able to identify cells that are CXCR5⁺ and express low levels of HAVCR2 (encoding TIM-3; see Materials and methods), CXCR5^– TIM-3^–, and CXCR5^– TIM-3⁺ T cells (reflecting those that were identified using flow cytometry). As the vast majority of CXCR5⁺ CD8⁺ T cells from melanoma expressed CD69 mRNA (Fig. S3 A), in line with protein data of lung tumors, we did not further subdivide between CD69^– and CD69⁺ for the analysis of gene expression (Fig. 3 A). In this way, we confirmed the preferential expression of TCF7, EOMES, BCL6, FOXO1, and CD28 by CXCR5⁺ T cells (Fig. 3, B and C) and revealing additional Tfh cell–associated genes such as SH2D1A (encoding SAP) and Ly9 (Fig. S3 B). Notably, a large proportion of these cells also expressed killer molecules such as GZMA, GZMK, GZMM, and PRF1 (the latter at slightly lower levels compared with CXCR5^– TIM-3⁺). Importantly, the same cells expressed inhibitory receptors such as PDCD1 (encoding PD-1) and TIGIT, albeit at lower levels compared with CXCR5^– TIM-3⁺ CD8⁺ T cells, whereas they lacked LAG3 and CTLA4 (Fig. 3 C). Besides CD28, also costimulatory receptors CD27 and TNFRSF9 (encoding 4-1BB) are expressed by CXCR5⁺ CD8⁺ T cells, yet the latter at a lower level compared with CXCR5^– TIM-3⁺ cells. Overall, CXCR5^– TIM-3^– cells had a pattern of gene expression that is intermediate between these two subsets (Fig. 3, B and C).

Gene set enrichment analysis (GSEA) revealed significant overrepresentation of germinal-center CD4⁺ Tfh and CD8⁺ naive, but not CD8⁺ effector T cell gene signatures within the CXCR5⁺ compared with CXCR5^– TIM-3⁺ CD8⁺ T cells (Fig. 3 D). Instead, naive and effector gene signatures where simultaneously more enriched when compared with CXCR5^– TIM-3^– (Fig. S3 C), further indicating that CXCR5⁺ CD8⁺ T cells have features of long-lived memory progenitors that are preferentially armed for effector functions. Most importantly, the genetic signature of murine CXCR5⁺ TIM-3^– cells isolated from chronically infected mice as from Im et al. (2016) was highly enriched in human CXCR5⁺ TILs compared with CXCR5^– TIM-3⁺ (Fig. 3 D) and CXCR5^– TIM-3^– CD8⁺ T cells (Fig. S3 C), further highlighting the similarity of human and murine subsets. In line with their follicular identity, the abundance of CXCR5⁺ TIM-3^– CD8⁺ T cells positively correlated with that of CD4⁺ CXCR5⁺ (comprising Tfh cells) and, at a lesser extent, of CD19⁺ B cells within the tumor (Fig. S3, D–G). As expected, the frequency of CD4⁺ Tfh highly correlated with that of CD19⁺ B cells, confirming their mutual relationship also in lung cancer (Fig. S3 H). Several of the transcripts found differentially expressed in melanoma-infiltrating CD8⁺ T cell subsets by scRNA-seq were then confirmed at the protein level in the NSCLC samples by using flow cytometry (TCF-1, Eomes, CD27, CD28, and PD-1; Fig. 3, E and F), thus further indicating that CXCR5⁺ CD8⁺ TILs share similar identity among different tumors. This flow cytometric analysis revealed additional differentially expressed markers that could not be measured by scRNA-seq due to low levels of transcripts in all populations (unpublished data), including T-bet (lower in CXCR5⁺) and CXCR3 (higher in CXCR5⁺). Most notably, CXCR3 is a signature marker of T_SCM cells (Gattinoni et al., 2011). Fig. S3 (I–K) additionally shows the phenotype of CD8⁺ TILs compared with naive (T_N) and effector memory T (T_EM) cells from the blood of a healthy donor, confirming that CXCR5⁺ TIM-3^– cells coexpress markers of T cell precursors along with intermediate expression of PD-1 and TIGIT. We conclude that CXCR5⁺ CD8⁺ TILs are armed for effector functions while retaining features of memory cells and express molecules that are candidate for immunotherapeutic intervention, such as PD-1, TIGIT, and CD27, but not LAG-3, CTLA-4, or TIM-3.

We next sought to determine whether the lack of IFNG, GZMH, and GZMB would result in decreased effector functions by CXCR5⁺ CD8⁺ TILs compared with other subsets. All sort-purified T cell subsets rapidly produced IFNγ, IL-2, TNF cytokines, and the degranulation marker CD107a upon PMA/ionomycin stimulation ex vivo (Fig. 4, A and B). In line with their exhausted phenotype, CXCR5^– TIM-3⁺ were the less functional in this regard, while CXCR5⁺ TIM-3^– expressed the highest level of CD107a (Fig. 4 B). Combinatorial analysis revealed the preferential coproduction of the four effector molecules by CXCR5⁺ TIM-3^– compared with both CXCR5^– TIM-3^– and CXCR5^– TIM-3⁺ cells (Fig. 4, C and D). These data indicate that CXCR5⁺ TIM-3^– CD8⁺ TILs retain rapid effector functions ex vivo and display enhanced polyfunctionality compared with other TIL subsets.

We next tested whether NSCLC-infiltrating CXCR5⁺ CD8⁺ T cells retain enhanced stem cell–like properties as suggested by their molecular signature. A hallmark of human memory stem cells is their capability to generate more effector progeny while self-renewing. Therefore, we sort-purified CXCR5⁺ TIM-3^–, CXCR5^– TIM-3^–_, and CXCR5^– TIM-3⁺ CD8⁺ T cells from tumor samples by flow cytometry (Fig. 5 A) and stimulated them with αCD3/28 antibody-conjugated beads in combination with IL-2 and IL-15 or IL-15 alone to determine their multipotency and self-renewal capabilities, respectively (Gattinoni et al., 2011). Fig. 5 (B and D) shows that CXCR5⁺ TIM-3^– and CXCR5^– TIM-3^– cells proliferated more vigorously than CXCR5^– TIM-3⁺ cells in response to TCR stimulation. These proliferating subsets displayed high levels of HLA-DR and CD25 activation markers and low levels of CD27 and CD28 (Fig. 5 K and Fig. S4), suggesting further differentiation. In line with murine data of adoptively transferred T cells in hosts bearing a chronic infection, only CXCR5⁺ TIM-3^– cells were able to generate all subsets of T cells identified by combinations of CXCR5 and TIM-3 expression (Fig. 5, E and H), thereby suggesting multipotency. In addition, these cells displayed lower levels of TIM-3 on the cell surface compared with other subsets (Fig. 5 G). CXCR5⁺ TIM-3^– proliferated more than CXCR5^– TIM-3⁺ cells also in response to IL-15 (Fig. 5 D). A similar trend could be seen for CXCR5^– TIM-3^–, although not reaching statistical significance. Importantly, the three subsets generated progenies with dramatic differences at the phenotypic level. While CXCR5⁺ TIM-3^– largely retained their original phenotype (Fig. 5, F and I), characterized by high levels of CD27 and CD28 costimulatory molecules and low levels of CD45RA (Fig. 5 L and Fig. S4), proliferating CXCR5^– TIM-3^– cells, the vast majority of which was CD45RA^– upon ex vivo purification, displayed high levels of CD45RA and, at a lesser extent, TIM-3, indicating differentiation toward the terminal effector phenotype (Fig. 5, J and L; and Fig. S4). Most notably, CXCR5 could not be detected in these cells, as well as in CXCR5^– TIM-3⁺ exhausted cells in response to IL-15, hence confirming the unidirectional CXCR5⁺ to CXCR5^– conversion. Thus, the tumor-infiltrating CXCR5⁺ CD8⁺ T cells have features of multipotency and stemness, as previously described for circulating T_SCM cells, in response to TCR-dependent and homeostatic stimulation, respectively.

The enhanced stemness of the tumor-infiltrating CXCR5⁺ CD8⁺ T cells suggest that these cells might be beneficial for the anti-tumor CD8⁺ T cell response, but their role in disease progression is still not defined. The maximum standardized uptake value (SUVmax), as obtained from positron emission tomography (PET) before surgery, is an indicator of the glycolytic capacity of the tumor and is considered a negative prognostic parameter for NSCLC (Takeuchi et al., 2014; Liu et al., 2016; Lopci et al., 2016). Indeed, patients with a pathological stage (pStage) I tumor, defined according to the international TNM classification, displayed a lower SUVmax compared with patients with pStage II or III in our cohort (Fig. 6 A). Therefore, we divided patients in two groups on the basis of the median value for this parameter across the entire cohort (high and low groups) for which both PET values and high-dimensional single cell profiling data were available (n = 33). We next determined the frequency of cluster 17 (i.e., CD8⁺ T cells with a preferential CD69⁺ CXCR5⁺ phenotype) in the tumor samples and revealed that this cluster is more prevalent in low glycolytic tumors (Fig. 6 B). Accordingly, we found that the same subpopulation was more prevalent within NSCLC samples of patients with early versus late disease stage (Fig. 6 C). In conclusion, the abundance of CXCR5⁺ CD8⁺ T cells in the tumor negatively correlates with disease progression.

We applied high-dimensional single cell approaches to analyze the tumor-infiltrating CD8⁺ T cells and define those phenotypic subpopulations that are preferentially enriched within the tumor compared with the blood and the adjacent cancer-free tissue. We adopted a novel flow cytometry technology capable to investigate 27 parameters simultaneously in millions of single T cells, a number still difficult to reach by CyTOF or scRNA-seq due to throughput and cost, respectively, and thus were able to identify also those rare subsets that are present in the TME. Specifically, we focused on CD8⁺ T cells that are amenable to immune checkpoint blockade (i.e., those that have traits of exhaustion) while retaining functional capacity, so to exploit their enhanced anti-tumor response. In this way, we revealed the presence of a subset of CD8⁺ T cells featuring a CXCR5⁺ TIM-3^– phenotype, long-lived memory characteristics, CD4⁺ Tfh-shared gene expression, and enhanced stem cell–like properties when compared with more differentiated subsets at the same site. The presence of these T cells exclusively in the TME is accompanied by increased frequencies of activated and exhausted phenotypes at the expense of CD27⁺ memory T cells and highly cytolytic, CD45RA⁺ T cells, which instead are normally abundant within the adjacent cancer-free lung tissue from the patient (Fig. 1 D). Our high-dimensional analysis identified a subset of nonneoplastic lung tissues that was more similar to blood samples rather than the tumors, possibly suggesting contamination from the circulating blood. Despite not being formally excluded, this is considered to be minor given the virtual absence of true naive cells (clusters 1 and 3) along with the presence of CD45RA^int/+ CCR7^– CD57⁺ phenotypes (clusters 18, 23, 4, and 8) that are normally abundant in lung tissues from healthy individuals (Thome et al., 2014) and of additional subsets that are rare in the blood (clusters 24, 14, 5, and 22). Overall, our data further confirm that the tumor is largely devoid of subpopulations capable to kill in an effective way, while supporting features of T cell exhaustion since the early disease stage. These results are in line with the recent characterization of a cohort of 73 renal cell carcinoma patients and a smaller cohort of 18 NSCLC patients by CyTOF, where patients at different stages of the disease were characterized by similar immune signatures at the level of T cells, including that of T cell exhaustion (Chevrier et al., 2017; Lavin et al., 2017).

The CXCR5⁺ population described here is the human equivalent of the partially exhausted, CXCR5⁺ PD-1^int CD8⁺ T cells expressing the transcription factor TCF-1, that, in chronically infected mice, are preferentially responsible for viral control following PD-1 blockade (He et al., 2016; Im et al., 2016; Leong et al., 2016; Utzschneider et al., 2016; Wu et al., 2016). Although, we could not directly address this enhanced responsiveness to PD-1 blockade by the tumor-infiltrating CXCR5⁺ CD8⁺ T cells due to practical reasons (i.e., low number of cells obtained after sorting), the flow cytometry and scRNA-seq data support this concept. By scRNA-seq of melanoma TILs and subsequent validation by high-dimensional flow cytometry and functional assays in lung TILs, we confirmed the follicular identity, a differentiation pattern similar to that of long-lived memory T cells, and partial exhaustion of single CXCR5⁺ CD8⁺ T cells and additionally revealed the simultaneous coexistence of stem cell–like behavior, polyfunctional effector molecule production, and cytotoxic potential. This is important because immune activation will result in immediate effector functions that can persist in the long term. Stem-like properties of T cells were shown to be beneficial in multiple instances, such as in the context of adoptive cell transfer where long-lived memory T cells are generally associated with superior anti-tumor responses (Gattinoni et al., 2005, 2011; Chapuis et al., 2013). Although CXCR5⁺ CD8⁺ T cells display a maturation phenotype that is largely different from that of canonical T_SCM cells, i.e., CD45RA^– CCR7^–, effector memory-like with tissue-resident characteristics, they also share some T_SCM properties, including the expression of CD27, CD28 and CXCR3 surface markers, and Eomes, but not T-bet transcription factors, self-renewal, and multipotency (Gattinoni et al., 2011; Lugli et al., 2013a; Fuertes Marraco et al., 2015; Oliveira et al., 2015; Pilipow et al., 2015; Akondy et al., 2017). Moreover, only CXCR5⁺ CD8⁺ TILs express TCF-1, a master regulator of persistence that is critical for the formation and maintenance of these cells in chronically infected mice. Despite acknowledging the limitations of in vitro assays to define a precursors–progeny relationship in humans, our data support the concept that a hierarchy of differentiation is also present at the tumor level in the context of exhaustion, albeit featuring different players, i.e., less differentiated CXCR5⁺ TIM-3^– stem-like CD8⁺ T cells giving rise to CXCR5^– and TIM-3⁺ effector subsets upon triggering of the TCR.

Fine mapping of the expression of inhibitory and costimulatory receptors with single cell resolution by both flow cytometry and scRNA-seq indicated that PD-1 and TIGIT could be preferential targets of immunotherapy with checkpoint inhibitors, while CD27 with immunostimulating antibodies favor effector differentiation (Ahrends et al., 2016). Differently, blockade directed to LAG-3 and CTLA-4, which are present on CXCR5⁺ CD8⁺ T cells from chronically infected mice, but are poorly expressed on human CXCR5⁺ CD8⁺ TILs, would be relatively inefficient. Despite not being formally demonstrated due to the difficulty in identifying tumor antigens, the presence of traits of exhaustion and the accumulation exclusively within tumors suggest that CXCR5⁺ CD8⁺ TILs are tumor-specific. It is at the moment unclear why CXCR5⁺ CD8⁺ T cells mainly form/accumulate in the TME compared with the cancer-free lung. Their follicular cell identity and their positive correlation with the abundance of Tfh and B cells within the tumor suggest they localize at the border of B cell areas in tertiary lymphoid structures, which have been shown in multiple instances to be associated with enhanced effector responses and increased survival in lung cancer (Dieu-Nosjean et al., 2008).

In conclusion, CXCR5⁺ CD8⁺ TILs are the ideal cells to be relieved from suppression as they serve as precursors of more differentiated effectors with enhanced functionality and persistence. However, they tend to be lost with disease progression. More differentiated CXCR5^– TIM-3^– CD8⁺ T cells retain some degree of functionality following stimulation despite their decreased self-renewal and degranulation potential and may serve as precursors of effector cells upon PD-1 blockade in advanced disease. However, this phenotype is poorly represented in the murine antigen-specific CD8⁺ T cell pool generated following chronic infection (He et al., 2016; Im et al., 2016; Leong et al., 2016; Utzschneider et al., 2016; Wu et al., 2016), thus raising the question whether tumor-specific CD8⁺ T cells are contained in the CXCR5^– TIM-3^– subset and can be primed intratumorally. Future studies should address whether tumor antigen-specific stem-like cells also exist in lymphoid tissues, and whether they are an additional reservoir of anti-tumor T cell responses.

All human experiments were approved by the Humanitas Clinical and Research Center Internal Review Board (approval nos. 1501, 1021, and 168/18). All patients provided written informed consent. In this research, we included 53 patients with diagnosis of NSCLC with clinical stage I or II who underwent a lung major resection (anatomical segmentectomy, lobectomy, or pneumonectomy) plus radical lymph node dissection with radical intent at the Humanitas Cancer Center. Tumor (n = 53) and adjacent cancer-free lung tissue samples from the same surgical specimen (n = 45) were obtained in the operating room in a sterile field by the surgeon while blood sample collection (10 ml; n = 22) were obtained by venipuncture at induction of anesthesia. Human intestinal specimens were obtained from patients with CRC undergoing resection surgery and liver tumor tissues from patients undergoing hepatectomy for hepatic metastatic disease from CRC. None of patients received preoperative chemotherapy or radiotherapy. Details on the patients’ characteristics are indicated in Table 1. Information on the pathological stage (pStage), determined by an institutional pathologist (some patients were restaged as III following examination of the tumor), was available for all patients, while results of the preoperative fluorodeoxyglucose (FDG) PET–computed tomography (CT) scan were available for 33 patients. PET scans were acquired ∼60 min after FDG administration in patients fasting for at least 6 h. Whole body images were obtained from the base of the skull to mid-thigh by means of integrated PET-CT tomographs. Reconstructed images were displayed on a GE ADW4.6 workstation (GE Healthcare). SUVmax was determined as the highest pixel value within the tumor mass that was identified with three-dimensional volumes of interest, as previously described (Lopci et al., 2016). SUVmax high and low groups (Fig. 6) were defined as those patients whose tumors have values that are above or below the median of the distribution, respectively.

Blood samples were collected in Vacutainer EDTA Tubes (BD), while the tumor and normal tissue samples were collected in RPMI-1640 medium supplemented with 10% FBS (Sigma-Aldrich), 1% penicillin-streptomycin, and 1% Ultra-glutamine (both from Lonza), hereafter referred to as R10. Peripheral blood mononuclear cells were isolated from the blood samples or from healthy donors (buffy coats) by Lympholyte-H (Cedarlane) density separation and frozen in liquid nitrogen according to standard procedures, as described previously (Lugli et al., 2013b). The tumor and normal tissue samples were processed into single cell suspensions by mechanical disaggregation. The samples were cut into small fragments in a size of ∼1 mm³, after which further disaggregation was performed using gentleMACS dissociator (Miltenyi Biotec). Subsequently, samples were passed through a 70-µm cell strainer (Falcon) and washed with physiological saline solution (NaCl 0.9%; Baxter). After centrifugation (300 g, 15 min), erythrocytes lysis was performed by incubating the cell suspension for 5 min in ACK lysis buffer (made in house; 155 mM NH₄Cl, 10 mM KHCO₃, and 0.1 mM EDTA-2Na [all from Sigma-Aldrich] diluted in H₂O), and washed again with saline. CRC samples were mechanically minced in small pieces, washed in HBSS without Ca²⁺and Mg² (Lonza) supplemented with 1% penicillin/streptomycin/amphotericin (Invitrogen) and followed enzymatic digestion with 1.5 µg/ml of collagenase IV (Life Technologies), 20 µg/ml of hyaluronidase, and 50 µg/ml of DNase (both from Sigma-Aldrich) for 2 h at 37°C/5% CO₂. Lymphocytes were separated by 30%/70% Percoll gradient centrifugation. Liver metastasis were cut in small pieces followed by enzymatic digestion in gentleMACS dissociator with 2 mg/ml collagenase D (Roche Diagnostic) for 45 min at 37°C/5% CO₂. Lymphocytes were isolated by 30%/70% Percoll density gradient. The cells were either frozen in liquid nitrogen according to standard procedures or, in case of sorting by flow cytometry, NSCLC leukocytes were further purified by centrifugation (790 g for 30 min without brake) on a 40%/70% Percoll (Sigma-Aldrich) gradient. After isolation of the interphase between the two layers, the cells were washed with HBSS without Ca²⁺ or Mg²⁺ (Lonza) and stained with fluorescently conjugated monoclonal antibodies (mAbs).

For the ex vivo immunophenotyping experiments, frozen samples were used and were thawed in R10 supplemented with 20 µg/ml DNase I from bovine pancreas (Sigma-Aldrich). After extensive washing with PBS (Sigma-Aldrich), the cells were stained immediately with the Zombie Aqua Fixable Viability kit (BioLegend) for 15 min at room temperature. Then, the cells were washed and stained with the combination of mAbs purchased from either BD Biosciences, BioLegend, or eBioscience, as listed in Table S1. mAbs were previously titrated to define the optimal concentration, as described (Lugli et al., 2013b). Chemokine receptors were stained for 20 min at 37°C, while all other surface markers (except CD3) were stained for 20 min at room temperature. Intracellular detection of CD3, Ki-67, and transcription factors was performed following fixation of cells with the FoxP3/transcription factor staining buffer set (eBioscience) according to manufacturer’s instructions and by incubating with specific mAbs for 30 min at 4°C.

In parallel, all samples were stained with a panel of mAbs directed to reveal the major immune populations, such as T cells (CD3⁺), B cells (CD19⁺), myeloid (CD11b⁺), and NK subsets (CD16, CD56, and NKp46). Additional panels used for further characterizing the CXCR5⁺ CD8⁺ T cell subset are listed in Table S1. Stainings were performed as described above, except for intracellular staining of TCF-1, for which cells were fixed with Transcription Factor Buffer set (BD Biosciences) according to the manufacturer’s instructions. For detection of cytokines, cells were fixed with BD Cytofix/Cytoperm Fixation/Permeabilization Solution kit (BD Biosciences) according to the manufacturer’s instructions. CD107a-specific antibody was added in the culture medium during the stimulation.

All data were acquired on a FACS Symphony A5 flow cytometer (BD Biosciences) equipped with five lasers (UV, 350 nm; violet, 405 nm; blue, 488; yellow/green, 561 nm; red, 640 nm; all tuned at 100 mW, except UV tuned at 60 mW) and capable to detect 30 parameters. Flow cytometry data were compensated in FlowJo by using single stained controls (BD Compbeads incubated with fluorescently conjugated antibodies), as described (Lugli et al., 2017).

Flow Cytometry Standard (FCS) 3.0 files were imported into FlowJo software version 9, and analyzed by standard gating to remove aggregates and dead cells, and identify CD3⁺ CD8⁺ T cells. 3,000 CD8⁺ T cells per sample were subsequently imported in FlowJo version 10, biexponentially transformed, and exported for further analysis in R (version 3.3.3) by a custom-made script that makes use of Bioconductor libraries and R statistical packages. Blood, normal lung tissue, and tumor samples were labeled with a unique computational barcode for further identification and concatenated in a single matrix by using the “cytof_exprsMerge function” (Chen et al., 2016). Data were analyzed using the Phenograph algorithm coded in the cytofkit package (version 1.6.5; Chen et al., 2016). Parameter K, indicating the number of nearest neighbors identified in the first iteration of the algorithm, was set at 60. Phenograph clusters were visualized using tSNE (cytofkit package, with default parameters). Clusters representing <0.5% were disregarded in subsequent analysis. The data were then reorganized and saved as new FCS files, one for each cluster, that were further analyzed in FlowJo to determine the frequency of positive cells for each marker and the corresponding MFI. These values were multiplied to derive the integrated MFI (iMFI, rescaled to values from 0 to 100; Fig. 1 D). The heat map, showing the iMFI of each marker per cluster, and the subsequent metaclustering was performed using the gplots R package. Hierarchical metaclustering of all samples, based on the frequency of Phenograph clusters (Fig. 1 C), was performed in R based on the Euclidean distance and Ward-linkage.

Normalized scRNA-seq counts were downloaded from Gene Expression Omnibus dataset (GSE72056). Analysis was restricted to the cells labeled as “T cells,” as previously defined by Tirosh et al. (2016). T cells were separated into CD4⁺ and CD8⁺ based on the normalized expression levels (E) of CD4 (E > 4), and CD8 (average of CD8A and CD8B, E > 4). The correlation between the expression profile of CXCR5 in the CD8⁺ TILs (n = 912), and that of all other genes present in the matrix was calculated using the “cor.test” function in R with Pearson as correlation methods. The CD8⁺ TILs (n = 912) were further divided in three groups based on CXCR5 and HAVCR2 (encoding TIM-3) expression levels. Specifically, we identified T cells that were CXCR5⁺ (E > 0; including all cells irrespectively of HAVCR2 expression, which was negligible in this subset; unpublished data), CXCR5^– TIM-3⁺ (E = 0 and E ≥ 5.930667, respectively), and CXCR5^– TIM-3^– (E = 0 and E ≤ 2.1527, respectively). Thresholds of HAVCR2 (TIM-3) expression were determined by calculating the 33.3 and 66.7 percentile of the HAVCR2 mRNA distribution. Therefore, if HAVCR2 mRNA levels were <33.3 percentile, cells were labeled as CXCR5^– TIM-3^–; otherwise, if >66.7 percentile, as CXCR5^– TIM-3⁺. Differentially expressed genes in each pairwise comparison of CD8⁺ T cell subsets were determined by the “FindMarkers” function coded in the Seurat R package (version 2.1.0; Satija et al., 2015) with default parameters except for the “logfc.threshold” that was set equal to 0. Subset-specific genes were defined as up-regulated in the pairwise comparisons between a cell subset and the other two. In this way we obtained 555, 46, and 1 specific gene for the CXCR5⁺, CXCR5^– TIM-3⁺, and CXCR5^– TIM-3^– populations, respectively. The heat map shown in Fig. 3 B was generated by using the “DoHeatmap” function coded in the Seurat R package (for the CXCR5⁺ subset, only the first 100 differentially expressed genes with the highest fold change resulting from the comparison with the CXCR5^– TIM-3⁺ or CXCR5^– TIM-3^– CD8⁺ T cells are shown). Violin plots, depicting the gene expression probability distributions across subsets, were obtained with the Seurat R package.

Microarray probe fluorescence signals of the publicly available dataset GSE84105 (sample IDs used for this analysis: GSM2227309, GSM2227310, GSM2227311, GSM2227312, GSM2227313, and GSM2227314) were converted to expression values using robust multiarray average procedure RMA (Irizarry et al., 2003) of Bioconductor affy package. Specifically, fluorescence intensities were background adjusted, normalized using quantile normalization and log₂ expression values for a total of 18,075 custom probe sets calculated using median polish summarization and custom chip definition files for Affymetrix Mouse Genome 430 2.0 arrays based on Entrez genes (Mouse4302_Mm_ENTREZG version 21.0.0). To identify genes that are differentially expressed, we compared the expression profiles of CXCR5⁺ TIM-3^– versus CXCR5^– TIM-3⁺ samples using limma algorithm coded in the same R package (Ritchie et al., 2015). By setting the false discovery rate and fold change thresholds at 0.05 and 2, respectively, we obtained 246 differentially expressed genes. Mouse Entrez IDs were converted in the corresponding human homologous genes using the HUGO Gene Nomenclature Committee database, and they were subsequently used for GSEA.

Over-representation analysis was performed using the GSEA tool described in Subramanian et al. (2005). GSEA was applied on the entire list of genes that compose the single cell RNaseq expression matrix (GSE72056). Genes were ranked based on their fold change calculated between groups of each pairwise comparison and analyzed by GSEA in preranked mode. We adopted the “classic” enrichment statistic, the recommended approach for RNA-sequencing data. Gene sets of interest (“GSE10239 NAIVE VS MEMORY CD8 TCELL UP,” “GSE21380 NON TFH VS GERMINAL CENTER TFH CD4 TCELL DN,” and “GSE10239 NAIVE VS DAY4.5 EFF CD8 TCELL DN”) were retrieved from the immunological signatures collection (c7.all.v6.1) of the Molecular Signatures Database.

CXCR5⁺ TIM-3^–, CXCR5^– TIM-3^–, and CXCR5^– TIM-3⁺, pregated as CD8⁺ Aqua^–, were isolated from freshly minced tumors with a FACSAria cell sorter (BD Biosciences) and stained with CellTrace CFSE kit (final concentration: 2 µM; Thermo Fisher) according to the manufacturer's protocol. Subsequently, cells were plated in R10 in U-bottom 96-well plates (10,000 cells/well) and stimulated with either αCD3/CD28 Dynabeads (Life Technologies; 2 beads/cell) in combination with IL-2 and IL-15 (Peprotech; 10 ng/ml each) or IL-15 alone (50 ng/ml) for 6 and 11 d at 37°C, respectively. T cells kept in low dose of IL-15 (1 ng/ml) served as non-proliferating, negative control. Cells were then stained with mAbs (Table S1) as described above. Proliferation index depicted in Fig. 5 D was calculated as follows: Proliferation index = MFI non-proliferating fraction / MFI proliferating fraction × % cells with diluted CFSE.

Freshly flow cytometry–sorted CXCR5⁺ TIM-3^–, CXCR5^– TIM-3^–, and CXCR5^– TIM-3⁺ were plated in R10 in U-bottom 96-well plates (10,000 cells/well) and stimulated with 10 ng/ml PMA and 0.5 µg/ml Ionomycin (both from Sigma-Aldrich) for 3 h at 37°C in the presence of Golgiplug (BD Biosciences) according to the manufacturer’s instructions. The antibody for the detection of degranulation marker CD107a was added to the culture medium. Cells were then stained with mAbs (Table S1) as described above.

Statistical analyses were performed using GraphPad Prism version 7, unless specified otherwise. Significance of differences for the frequency of single Phenograph clusters among blood, normal tissue, and tumor samples was determined using two-way ANOVA with Bonferroni post-hoc test. To compare distributions of manually gated subsets significance was determined by paired Wilcoxon t test, unless otherwise specified in the figure legends. Correlation of data were determined by calculating the Pearson correlation coefficient. P values <0.05 were considered statistically significant.

Fig. S1 shows the phenotypic analysis of CXCR5⁺ CD8⁺ T cells. Fig. S2 shows that bona fide CD8⁺ T_SCM cells are absent in lung tumors and adjacent tissues. Fig. S3 shows that the gene signature of CXCR5⁺ CD8⁺ T cells resembles that of both CD4⁺ Tfh cells and memory CD8⁺ T cells. Fig. S4 shows that proliferating CXCR5⁺ TIM-3^– CD8⁺ T cells are self-renewing and multipotent. Table S1 lists all of the antibodies for flow cytometry used in this study.

We wish to thank the laboratory of Massimiliano Pagani (INGM, Italy) for helping with single cell RNA-sequencing data analysis, Gianluca Rotta, and Jens Fleischer (BD Biosciences) for helping with FACSSymphony A5 instrument setup; Matteo Donadon and Guido Torzilli (Liver Unit, Humanitas) for colorectal liver metastasis samples; Antonino Spinelli (Colon and Rectal Unit, Humanitas) for colorectal carcinoma samples; Paola Spaggiari and Massimo Roncalli (Pathology Unit, Humanitas) for histological analysis; and Prof. Alberto Mantovani (Humanitas University) and the members of the Laboratory of Translational Immunology for critical discussion.

This work has been supported by grants from the Associazione Italiana per la Ricerca sul Cancro (AIRC; grants IG20607 to E. Lugli, IG14687 to D. Mavilio, and MFAG189323 to E. Lopci), Bristol Myers Squibb (preclinical grant no. CA209-9YM to E. Lugli), the Italian Ministry of Health (grant PE-2016-02363915 to D. Mavilio), and the Humanitas Clinical and Research Center grant to E. Lugli and D. Mavilio. J. Brummelman is a recipient of the “Fondo di beneficenza Intesa San Paolo” fellowship from AIRC. E.M.C. Mazza and P. Novellis are recipients of the 2017 Fondazione Umberto Veronesi postdoctoral fellowship. E. Lugli is an International Society for the Advancement of Cytometry Marylou Ingram scholar. Purchase of the BD FACSSymphony A5 has been in part defrayed by a grant from Italian Ministry of Health (agreement no. 82/2015).

The Laboratory of Translational Immunology receives reagents in kind from BD Biosciences as part of a collaborative research agreement. The authors have no additional financial interests.

Author contributions: E. Lugli and G. Veronesi wrote the protocol for ethical evaluation. M. Alloisio, P. Novellis, and G. Veronesi performed the surgeries. E. Lopci performed and generated the FDG PET-CT scans. J. Brummelman and G. Alvisi processed the patients’ samples. J. Brummelman and E. Lugli designed the experiments. J. Brummelman and G. Alvisi carried out the experiments. F.S. Colombo performed the quality check of the FACSSymphony A5 and assisted with sort experiments. J. Brummelman and G. Alvisi analyzed the flow cytometry data. E.M.C. Mazza performed the tSNE and Phenograph analysis and processed the single cell RNA-sequencing data with assistance of F. Ferrari and A. Grilli. D. Mavilio and J. Mikulak performed analysis on colorectal and liver metastasis tumors. J. Brummelman and E. Lugli interpreted the data and wrote the manuscript with input of all authors.