A high OXPHOS CD8 T cell subset is predictive of immunotherapy resistance in melanoma patients

The development of immune checkpoint inhibitor (ICI) therapies, including anti-programmed cell death 1 (PD-1)/programmed death-ligand 1 (PD-L1) and anti-cytotoxic T lymphocyte antigen 4 (CTLA4), represents a major advance in treating patients with melanoma. However, not all patients respond to these therapies (Garon et al., 2019; Garon et al., 2015; Larkin et al., 2019; Ribas and Wolchok, 2018), indicating that additional strategies are needed to improve responses. Among the immune repertoire, CD8 T cells play critical roles in ICI therapies and thus are an important cell type to focus on to improve immunotherapy (Wei et al., 2017).

We and others showed that tumor-infiltrating lymphocytes (TILs), including CD8 T cells, express high levels of immune checkpoint receptors that include PD-1, which can suppress cytotoxic T cell effector function (Maybruck et al., 2017; Pfannenstiel et al., 2019; Wherry, 2011; Wherry and Kurachi, 2015; Zhang et al., 2011). These TILs often display so-called exhausted behaviors characterized by antigen unresponsiveness and reduced cytotoxic effector function (Sharma et al., 2017; Wherry, 2011; Wherry and Kurachi, 2015). While therapeutic blockade of PD-1/PD-L1 and/or CTLA4 can revive effector T cell function, this is not always complete and may explain the lack of anti-tumor therapeutic efficacy (Larkin et al., 2019; Wei et al., 2018; Zimmer et al., 2017).

Importantly, CD8 T cell tumor infiltration remains one of the most correlated factors for anti-PD-1/PD-L1 immunotherapy response across cancer types, as well as tumor mutational burden and high PD-1/PD-L1 expression levels (Gros et al., 2014; Lee and Ruppin, 2019; Samstein et al., 2019). Today, no predictive model of ICI response exists that is robust enough to implement in the treatment algorithm for melanoma patients (Agur et al., 2016; Kogan et al., 2012; Liu et al., 2019). Other studies have focused on gene expression and transcriptomic data of T cells to predict patient response to anti-PD-1 therapy (Auslander et al., 2018; Gide et al., 2019; Hugo et al., 2016; Sade-Feldman et al., 2018a); these strategies required invasive biopsies and resulted in moderate prognostic value (Auslander et al., 2018; Hugo et al., 2016; Topalian et al., 2016).

Data from the National Cancer Institute (NCI) surgery branch indicate that there are overlapping CD8 subpopulations both in the TIL and in the periphery, which was recently corroborated at the transcriptomic level (Fairfax et al., 2020; Gros et al., 2016; Wu et al., 2020). Comparing our studies of tumor-associated dysfunctional CD8 T cells with suppressor cell function and those from the NCI surgery branch (both studies identified PD-1 expression as critical) led us to hypothesize that this population exists in both the tumor and peripheral blood (Gros et al., 2016; Gros et al., 2019; Pfannenstiel et al., 2019; Zhang et al., 2011). We further hypothesized that identification of such a subpopulation would elucidate new genetic mechanisms that could be used to not only monitor treatment outcome but also identify future therapeutic targets. Indeed, our single-cell transcriptome comparisons between purified CD8 T cells from TILs and peripheral blood lymphocytes (PBLs) from melanoma patients identified several CD8 subpopulations and underlying genetic programs. Specifically, we identified three overlapping populations in TILs and PBLs. One of the populations displayed unexpectedly high levels of cytotoxic and exhausted markers (e.g., PD-1), as well as increased levels of metabolic activities, specifically oxidative phosphorylation (OXPHOS). These high-OXPHOS CD8 T cells also have elevated levels of CD38/CD39 ectonucleotidases (we coin these CD8⁺ T_OXPHOS cells). To further validate our findings, we isolated PBL and TIL CD8⁺ T_OXPHOS cells from melanoma patients and confirmed that they are more pronounced in ICI-resistant patients with elevated bioenergetics, including increased glucose metabolism, ATP production, and mitochondria activity. We then generated a predictive ICI therapy response model using these CD8⁺ T_OXPHOS cells that can be assessed by either TILs or PBLs. Ultimately, we generated a gene expression profile (GEP), which was validated using both published clinically annotated transcriptomic data of CD8 TILs and a new cohort of our patients, including interrogating their CD8 PBLs. Thus, through comprehensive and careful gene expression analysis of individual CD8 T cells, we establish an ICI response GEP that can be employed via a blood-based approach. Our work establishes an ICI-predictive platform, and the CD8 subpopulation that forms the basis of this assay illustrates new targetable pathways to potentially enhance immunotherapy and improve outcomes of melanoma patients.

To develop a noninvasive method to predict effective immunotherapeutic responses, we carefully characterized individual CD8 T cells from tumors and peripheral blood of melanoma patients. We performed single-cell RNA sequencing (scRNaseq) on isolated CD8 T cells from both melanoma PBLs (CD8-mPBLs) and TILs (CD8-mTILs) from eight advanced melanoma patients (Fig. 1 and Table S1). Compared with recent scRNaseq studies of melanoma (Sade-Feldman et al., 2018b), our dataset is uniquely tailored to capture data from CD8 T cells with paired PBLs and TILs from the same patient. Through barcoding and combined analyses, we ensured matched processing and efficient downstream deconvolution of both samples. CD8-mPBLs and CD8-mTILs from the same patient were set as pairs and analyzed with flow cytometry for their proportion correlation analyses, and consistent with previous reports, the proportions of total CD8 T cells in PBLs were moderately correlated with those in TILs (R = 0.29) among the examined patients (Fig. 2 A and Fig. S1; Shao et al., 2014).

Next, we combined scRNaseq results of 173,061 CD8-mPBL/mTILs and performed transcriptomic analysis. Using t-distributed stochastic neighbor embedding (t-SNE) visualization, we uncovered a high level of heterogeneity of CD8 T cells in CD8-mPBL/mTIL, with a total of 20 clusters (clusters with <1% of cells in all samples were removed for further analysis; Fig. 2, B–D; Fig. S1; Fig. S2; and Table S2). Interestingly, both CD8-mTILs and CD8-mPBLs appeared activated in TCR signaling, components of the phosphatidylinositol 3-kinase pathway, CD28 and OX40 signaling, MHC II and NF-κB signaling, and the NRF2-mediated oxidative stress response (Fig. 2 E, Table S3, and Table S4), based on the significant differential expression of genes (false discovery rate [FDR]–adjusted t test, P < 0.05) between CD8-mTILs and CD8-mTILs. In addition to enhanced pathways for T cell activation, CD8-mTILs exhibited traits of T cell exhaustion, stress responses, apoptosis, and suppressed immune checkpoint and PPAR signaling pathways (Fig. 2 E and Table S3). These observations are a testament to scRNaseq, showing that tumor-bearing conditions impart a systemic impact on circulating immune cells in addition to tumor-infiltrating cells (Bai et al., 2015; Lee et al., 1999; Manjarrez-Orduño et al., 2018).

CD8-mPBLs and CD8-mTILs were quite distinct, with only three clusters shared between them (Fig. 2 D and Fig. 3 A). These three CD8-mPBL/mTIL shared clusters (clusters 2, 6, and 15) may permit tracking of CD8 T cells and help decipher their intracellular programming when exposed to the tumor environment, thus facilitating identification of important factors that dictate patient responsiveness to immunotherapy. To test this, we first analyzed the transcriptomic correlation of clusters 2, 6, and 15 with all other clusters and found that cluster 2 showed similarity with most other clusters, with the exception of clusters 6 and 15 (Fig. 3 B). Cluster 6 showed similarity with clusters 5, 12, 14, 18, and 20, whereas cluster 15 was similar to clusters 4, 9, 16, and 19. This nonoverlapping correlation of these three clusters supports their unique cellular programs and potential functions. Indeed, the 3D plot of all cells clearly depicts their distinct transcriptional profiles (Fig. S2, A–C). Next, we extracted signature genes (FDR-adjusted t test, P < 0.05) of the clusters (Fig. 3, C and D; and Table S5) that were found in all eight patients (Fig. S2, E and F) and performed a pathway analysis (Fig. 3 E, Table S6, and Table S7).

Cluster 2 expressed genes consistent with naive or resting T cell behavior, such as similar activation levels of DNA damage checkpoint regulation. Housekeeping RhoGDI (Rho GDP–dissociation inhibitor) pathways that are necessary to inhibit T cell activation were also highly active in cluster 2 cells (Allenspach et al., 2001; Burkhardt et al., 2008; Lin et al., 2003; Fig. 3 E). Interestingly, HIPPO signaling was activated in cluster 2 (Fig. 3 E), which includes the key inhibitor MOB1A, subunits of protein phosphatase, and 14–3-3, confirming the housekeeping function of these genes is consistent with that of naive or resting CD8 T cells. In addition, cluster 2 also contained suppressed Cdc42 signaling and T cell exhaustion signaling pathways that would ensure normal cell cycle and prevention of apoptosis (Fig. S2 F and Table S4).

Cluster 15 displayed pathway signatures mostly opposite to clusters 2 and 6, with prominent activation of OXPHOS, as well as activated glycolysis and NER (nucleotide excision repair) pathway (Fig. 3 E). Contrary to clusters 2 and 6, inhibition of HIPPO and RhoGDI pathways were also observed in cluster 15 (Fig. 3 E). Not surprisingly, T cell exhaustion markers were highly up-regulated in cluster 15, including the immune checkpoint molecules PD-1, TIM-3, LAG3 (Fig. 3, D and F), LAYN, and CXCL13 (Fig. 3 F). Remarkably, ectonucleotidases CD38 and CD39, which are involved in NAD⁺ and ATP regulation, were also significantly up-regulated in cluster 15 (Fig. 3 F). These characteristics of having not only expression of exhaustion markers but also clear metabolic activities add to our previous work of functionally suppressive CD8 T cells with high levels of checkpoints (Pfannenstiel et al., 2019).

Compared with clusters 2 and 15, cells in cluster 6 were transcriptionally suppressed, with fewer than 100 genes (out of 11,000 total detected gene IDs in the whole dataset) up-regulated compared with all other clusters. Furthermore, most signaling pathways were inhibited or had low activation, and only a few pathways remained relatively active, such as the enhanced Sirtuin, RhoGDI, and sumoylation pathways. These analyses suggested that some CD8 T cells in mPBLs resemble CD8 T cells in mTILs, with three similar populations that are transcriptionally distinct from each other. To examine the potential bias effect of different numbers of T cells from different patients, we also analyzed our dataset with a t-SNE process using an equal number of cells from each patient. Consistently, we not only observed comparable cluster distributions but also were able to identify the three PBL/TIL shared clusters with similar phenotypes and properties (Fig. S3). Thus, we show that that these three overlapping clusters are not due to an artifact of sample selection, a multidimensional reduction process, or a unique effect of a limited number of patients but are a phenomenon of the whole melanoma cohort in this study.

To further evaluate the cell state transition, we randomly selected equal number of cells (1,000 cells from each sample) to perform the trajectory analysis using the Monocle package (Qiu et al., 2017) to reveal cell fate differentiation along the branch point trajectories. The overall distribution of CD8-mPBLs and CD8-mTILs formed a three-branched plot, with cluster 2, 6, or 15 at the tip of each branch (Fig. 4, A and B). Interestingly, CD8-mPBL continuously span from cluster 2 to cluster 15, with few isolated cells extending toward cluster 6, suggesting that cluster 2→15 is a tumor-specific pattern (Fig. 4, B and C). Further, CD8-mTILs are located along all three branches (cluster 2→15, cluster 2→6, and cluster 15→6), suggesting that CD8-mPBLs in cluster 6 are likely those migrating to and from tumors and the periphery. Since cell trajectory is calculated based on overall transcriptome profiles of each cell and assigned their relative position based on relative simulation, the location of clusters 2, 6, and 15 on the plot further confirms that these clusters represent the three most distinct transcriptomic patterns, suggesting three distinct function types. This adds to the importance of clusters 2, 6, and 15 as the only clusters that have overlapping populations between TILs and PBLs. Setting the branch of cluster 2 as root in Monocle, we tracked whole-transcriptome gene expression changes along the trajectory branches and identified gene clusters that were altered along the pseudotime progression from cluster 2→6 or cluster 2→15 (Fig. 4 D). Interestingly, immune checkpoint genes, including CTLA4, PD-1, and LAG3, were highly expressed at the cluster 15 branch, whereas the PIK3IP1 gene, which is important for suppressing T cell activation, displayed an increased expression pattern in the cluster 2→6 direction (Uche et al., 2018). In contrast, genes contributing to apoptosis and T cell exhaustion pathways (CALM1, HLA-DQA1, HLA-DRA, CTLA4, PD-1, LAG3, TIGIT, and TIM-3) were highly enriched in cluster 15, but not in cluster 2 or 6 (Fig. 4, D and E; and Fig. S4).

By pseudotemporal trajectory analysis, we discovered three distinct programs that influence the gradual transition between the PBL/TIL shared clusters 2, 6, and 15 (Fig. 5 A). In CD8-mTILs, program 1 progressively transforms naive T cells (cluster 2) into a more “exhausted” CD8 T cell state (clusters 13, 20, and 19) along the trajectory and eventually transition to cluster 15 at the most distal of the trajectory branch (Fig. 5 A). While program 2 drives the transition of naive CD8 T cells to cluster 6, program 3 transforms cluster 15 to cluster 6 (Fig. 5 A). However, in CD8-mPBLs, only one program induces CD8 T cell transition, where program 1 is shown to trigger the transition of naive CD8 T cells (cluster 2) to cluster 15 (Fig. 5 A). Of note, only program 1 significantly fosters the transition of naive CD8 T cells (cluster 2) into an exhausted CD8 T cell state (cluster 15) in both PBLs and TILs, with other clusters distributed along the program transitional states, indicating that CD8-mPBLs display certain molecular and gene profiles resembling CD8-mTILs (Fig. 5 A).

Subsequently, we further analyzed the expression patterns of cell surface proteins and cytokines, which are vital in dictating CD8 T cell characteristics and functions in all the three programs (1, 2, and 3). Cytokines or secreted polypeptides that are known to maintain a resting T cell state, such as NOG and LTB, are enriched in cluster 2, whereas secreted molecules associated with T cell activation, such as CD40LG, CCL28, FLT3LG, TNFSF8, and NRG2, displayed an extended pattern along program 1 (cluster 2→15). Cytokines important for a strong inflammatory response in tumors are enriched at the cluster 15 end of program 1 (Fig. 5 B). Further, dynamic changes of cell surface receptors also displayed similar patterns (Fig. 5 B). Overall, cluster 15 CD8 T cells expressed high levels of genes for MHC II (HLA-DRB1, HLA-DKA, HLA-DPA1, and CD74), receptors for cell-to-cell interaction (CD2, CD47, CD27, IL2RG, and CXCR3), and dysfunctional markers and immune checkpoints such as PD-1, CTLA4, LAG3, and CXCL13. A recent review reported that CXCL13 was overexpressed in terminal dysfunctional CD8 T cells, suggesting a positive correlation between CXCL13 and dysfunctional anti-tumor efficacy (van der Leun et al., 2020). Interestingly, IL6R, IL7R, and CCR7 were down-regulated from resting/naive cells (cluster 2) to cluster 15 CD8 T cells (Fig. 5 B). More importantly, we identified key immune checkpoint genes and transcription factors along the three program directions in CD8-mPBLs and CD8-mTILs (Fig. 5 C). For example, in these programs, we identified TOX, which is a critical transcriptional factor that reprograms and drives CD8 T cells into an exhausted state during cancer progression and chronic infection (Mann and Kaech, 2019). We found that TOX increases along program 1, along with multiple checkpoints, and is highest in cluster 15 CD8 T cells (Fig. 5 C).

By transcriptomic analysis, we next investigated the metabolic signaling in all clusters. Strikingly, we observed a clear pattern of cell metabolic shifting during cell transition from naive or effector state to a dysfunctional exhaustion state (Fig. 6). Among all clusters, cluster 15 displayed elevated OXPHOS, glycolysis, glucose, and lipid transportation activity (Fig. 6 A). Genes involved in metabolic pathways, including OXPHOS, glycolysis, glucose, and lipid metabolism, were up-regulated in multiple clusters leading to cluster 15, such as PGK1, ATP5F1A, SLC2A1, and LDLR (Fig. 6 A). In fact, cluster 15 had the highest cell metabolic rate, with the most pronounced activation of glycolysis and OXPHOS pathways (previously shown in Fig. 3 E; see also Fig. 6, B–D). Remarkably, the OXPHOS pathway underwent a more significant change relative to glycolysis that enabled cells to transition from the state represented by cluster 2 to the state represented by cluster 15 (Fig. 6 E). There is visible segregation of a high OXPHOS and glycolysis activation state in clusters 4 (TILs), 15 (shared PBLs/TILs), and 16 (PBLs) from the other clusters, with cluster 15 having the maximal activation of both OXPHOS and glycolysis (Fig. S4 A).

Cluster 15, enriched with exhaustion markers (shown in Fig. 4 D), has prominent levels of both OXPHOS and exhausted markers (Fig. 6 F). Similar to cluster 15, we also observed high OXPHOS and exhaustion levels in both clusters 4 and 16 (Fig. 6 F). On the other hand, clusters 14 and 9 in TILs have high exhaustion levels but low OXPHOS activity (Fig. 6 F), thus indicating that they are classical exhausted TILs (van der Leun et al., 2020; Wherry and Kurachi, 2015). Commonly, exhausted T cells have low or mild expression of cytotoxic genes, which is very similar to clusters 4 and 9 (Fig. S4 B). However, we unexpectedly observed that clusters 4, 15, and 16 are enriched with cytotoxic genes (e.g., PRF1 and GZMB); in particular, IFN-γ is elevated in cluster 15 (Fig. 3 D), suggesting that these three clusters are distinctly different from classic exhausted T cells. Overall, our data demonstrate that clusters 4, 15, and 16 conglomerate into a unique T cell subset characterized by high exhaustion and OXPHOS states as well as increased cytotoxic gene expression (Fig. 6 F; and Fig. S4, A–C).

We next examined the biological progression of the various clusters based on their metabolic and transcriptional programs using pseudotime trajectory analysis. Consistent with the data shown in Fig. 3, we noticed a trifurcation of CD8 T cells into three distinct branches: resting/naive (cluster 2 branch, including clusters 1, 2, 12, 17, and 18), high OXPHOS with exhaustion markers (cluster 15 branch, including distal and related clusters 4, 15, and 16 [circled in Fig. 6 G]), or inactive dormant (cluster 6) cells. Other clusters are distributed along the program transitional states (Fig. 6 H). These results suggest that transcriptional programs in CD8-mPBLs have features similar to those found in CD8-mTILs, specifically in the formation of cluster 15, with clusters 4 and 16 having the closest proximity to cluster 15 (Fig. 6 H).

To functionally identify these unique subset of T cells, we looked for potential cell surface markers that could differentiate the transcriptionally related clusters 4, 15, and 16 from other clusters. By trajectory plot, we observed a gradual enrichment of PD-1, CD38, and CD39 expression leading toward cluster 15, inferring that these markers could be used for stratification (Fig. S4 D). We indeed found that both CD38 and CD39 are able to distinguish clusters 4 and 15 from other clusters, but not in cluster 16, which is found in PBLs (Fig. 6 H). Consistently, clusters 4 (TILs) and 15 (shared PBLs and TILs) also have high levels of combined expression of PD-1 and CD38/CD39, but not cluster 16, which may have to do with not being under the direct influence of the tumor microenvironment (Fig. S4, E and F). Together, clusters 4, 15 and 16 form a unique spectrum of CD8 T cells defined transcriptionally by their distinct features of high OXPHOS, exhaustion (e.g., PD-1), and cytotoxic gene expression, which we coin CD8⁺ T_OXPHOS. However, for cellular-level evaluation, to validate these transcriptional findings, the ability to identify PD-1, CD38, and CD39 expression will at least allow for studies of CD8⁺ T_OXPHOS clusters 4 and 15.

We reported that many melanoma patients have a significant number of CD8 TILs expressing high levels of immune checkpoints and that these same cells have active immune suppression of autologous healthy lymphocytes (Pfannenstiel et al., 2019). In our single-cell transcriptomic analysis, we found that CD8⁺ T_OXPHOS cells have the strongest metabolic signal within CD8 TILs. These findings led us to explore whether purified CD8 TILs from patients enriched with CD8 T cell subsets expressing high levels of immune checkpoints have a metabolic rate. Due to the limitations of the Seahorse glucose oxidation assay, which requires a large number of cells, we did a simple comparison as a starting point with healthy T cells.

We first performed the Seahorse assay to validate the involvement of glucose oxidation during glycolysis by determining the oxygen consumption rate (OCR) of healthy CD8 PBLs and CD8⁺PD-1⁺ TILs (Fig. 6 A shows that CD8⁺ T_OXPHOS cells account for the largest signal of glucose transportation in melanoma patients’ CD8 T cells). Melanoma CD8 TILs were indeed fueled by glucose oxidation as validated by higher OCR compared with healthy PBLs (Fig. 7 A). Of note, these CD8 TILs not only depended on glucose oxidation but also had high initial OCRs (Fig. 7 A). Interestingly, CD8 TILs had strong dependency on oxidation of glucose, as demonstrated by increased glucose dependency rate (Fig. 7 B). Glycolysis, in concert with OXPHOS, generates energy through ATP production. Thus, we measured ATP production in TILs and PBLs of melanoma patients. Using an ATP assay, we validated that peripheral CD8 T cells of ICI-nonresponder melanoma patients produced significantly higher ATP as compared with responders (Fig. 7 C), although we do note that the changes are not as striking in PBLs (P < 0.01) compared with TILs (P < 0.0001). Our biobank of surgical specimens contains TILs mainly from two sources: (1) ICI nonresponders and (2) patients who have not yet undergone systemic therapy (naive). CD8 TILs from nonresponders produced significantly higher ATP compared with naive patients (Fig. 7 C). The reduction of differential membrane potential is equivalent to a reduced capacity for ATP production and other related mitochondria functions (Little et al., 2018). Hence, to answer whether these high-bioenergy-state CD8 T cells were generated by the CD8⁺ T_OXPHOS cell subpopulation within them, we performed flow cytometry combined with metabolic analyses using tetramethylrhodamine methyl ester (TMRM), a fluorescent dye that accumulates in the functional mitochondria caused by differential membrane potential, as a marker of OXPHOS (Davis et al., 2020; Scaduto and Grotyohann, 1999). PBL and TIL CD8 T cells were costained with enriched-surface markers identified in our scRNaseq analysis, including PD-1, CD38, and CD39 (Fig. S5 D), and MitoTracker to determine the mitochondrial membrane potential relative to mitochondrial total mass.

Although expression of CD38 and CD39 in CD8⁺PD-1⁺ T cells does help characterize the CD8⁺ T_OXPHOS cells, additional transcriptional assays are required for exact determination of this immune population. We compared CD8⁺ T_OXPHOS cells to other CD8 subpopulations; they expressed a unique transcriptomic signature, even when compared with effector and central memory cells as well as exhausted CD8 T cells (Fig. S5 A). Still, within this PD-1⁺ subgroup (that have high CD38/39), we observed increased levels of OXPHOS (Fig. 7 D), consistent with our CD8⁺ T_OXPHOS cells designation, which is distinct from the classically defined PD-1⁺ exhausted T cells, which have been shown to have impaired mitochondrial function (Schurich et al., 2016; Vardhana et al., 2020). Within this population of CD8⁺PD-1⁺CD38⁺CD39⁺ T cells, we illustrate two levels of TMRM: (1) high TMRM (TMRM^hi) and (2) low TMRM (TMRM^lo; for a schematic representation of flow cytometry in TILs, see Fig. 7 D). Nonresponder patients mainly have TMRM^hi and negligible TMRM^lo levels, while naive patients have higher TMRM^lo but lower TMRM^hi in TILs (Fig. 7 D). Although nonresponders have higher TMRM^total in PBLs, no obvious difference was observed in CD8⁺PD-1⁺CD39⁺ and CD8⁺PD-1⁺CD38⁺CD39⁺ TILs (Fig. S5 B). TMRM^lo in CD8⁺PD-1⁺CD39⁺ and CD8⁺PD-1⁺CD38⁺CD39⁺ T cells in both PBLs and TILs do not stratify for nonresponders (Fig. S5 C). Most importantly, significantly higher mitochondrial membrane potential was observed in PBL CD8⁺PD-1⁺CD39⁺ and CD8⁺PD-1⁺CD38⁺CD39⁺ T cells of nonresponders of melanoma patients (Fig. 7 E). However, CD38 as a marker did not correlate with high mitochondrial membrane potential in PBLs (data not shown). In TILs, nonresponders have elevated mitochondrial membrane potential in both CD8⁺PD-1⁺CD39⁺ and CD8⁺PD-1⁺CD38⁺CD39⁺ cells, but not in CD8⁺PD-1⁺CD38⁺ T cells, compared with naive melanoma patients (Fig. 7 F).

Importantly, our data suggest the high possibility of stratifying responders and nonresponders using PBL-based metabolic activity assays (ATP and TMRM). However, we did not observe significant stratification between patients before therapy (naive) and after failed ICI therapy (nonresponders). The likely reason is that the naive group contains both responders and nonresponders, masking the stratifying effect. With longer follow-up, we will compare naive patients who did or did not respond to ICI; this is part of an ongoing and distinct effort. Further, these findings highlight the possibility that these metabolic pathways are not only phenotypic but also functional. Therefore, targeting them to modulate CD8⁺ T_OXPHOS cells could be a novel way to improve ICI therapy.

Based on the association between CD8⁺ T_OXPHOS cells and ICI resistance, we further developed a predictive model of immunotherapy responses using signature genes enriched in cluster 15–like cells, which include not only clusters 4 and 15 but also PBL cluster 16. To screen for a predictive immunotherapy response gene signature, we first performed whole-transcriptome gene coexpression analyses in the three unique PBL/TIL shared populations (clusters 2, 6, and 15). Genes that displayed significant positive correlation with PD-1 expression were highly expressed in clusters 4 and 16 but were most pronounced in cluster 15 (Fig. 8 A and Fig. S5 D). By gene ontology analyses, genes coexpressed with PD-1 in CD8⁺ T_OXPHOS cell clusters were enhanced in the mitochondrial dysfunction pathway, Cdc42 signaling, and CTLA4 signaling (Fig. 8 B). Interestingly, OXPHOS was the most activated pathway in cluster 15–like cells, further substantiating the importance of immune metabolism in checkpoint-based immunotherapy.

To develop a GEP predictive of immunotherapy response, we selected the most significantly elevated top 20 coexpressed with PD-1 (Fig. 8 A). One of these is the immune checkpoint LAG3, which binds to MHC II and is associated with exhausted CD8 T cells in human tumors (Chen and Mellman, 2017). Another, GAPDH, a key enzyme of glucose metabolism, was overexpressed in CD8⁺ T_OXPHOS cells. We also identified several genes less known for their roles in cancer ICI therapy (including serglycin, a small proteoglycan important for cytotoxic T cell secretory function, and cystatin F [CST7], a cysteine peptidase inhibitor) that may regulate immune cell function in tumors (Grujic et al., 2005; Perišić Nanut et al., 2017). Besides CST7, CD74 and HLA genes (HLA-A and HLA-C) known to be crucial for CD8 T cell function were also found to be up-regulated in both CD8 PBLs and TILs with high PD-1 coexpression.

Using the 21-gene GEP (including PD-1), we built a logistic regression model using a training cohort from a published study with scRNaseq dataset of TILs from melanoma patients undergoing immunotherapy (GSE120575; Sade-Feldman et al., 2018b). Our model is designed to predict the status of each T cell in a patient. Then, the proportion of nonresponding cells in each patient was rescaled to 0∼10, depicting the nonresponse score (NRS; see Materials and methods for details). A median NRS of 5 was selected as the threshold score, and patients with an NRS >5.0 were considered nonresponders. Using the training dataset, we predicted the ICI response in 11 of the 12 patients, yielding an accuracy of 92%, with nonresponders having significantly higher NRSs than responders (t test, P = 0.02; Fig. 8 C).

To further validate our model, we applied it to four additional datasets, including one published dataset and three independent validation patient cohorts from our institution. Due to the limited availability of clinical immunotherapy response data from public repositories, we validated our model using a published scRNaseq dataset of TILs from nonmelanoma skin cancer patients receiving anti-PD-1 immunotherapy with only one prediction error (GSE123813; dataset 1; Fig. 8 C; Yost et al., 2019). Next, we went on to further examine the accuracy of our model using patient cohorts from our institution. We first employed scRNaseq of CD8⁺PD-1⁺ T cells from five nonresponder patients with matched PBLs–TILs (dataset 2: five PBLs and five TILs; note three of these patients [six samples, paired TILs–PBLs] were from our original discovery set [GSE138720]), and then we performed scRNaseq on two more refractory patients’ TIL and PBL CD8 T cells (GSE153098). In addition, we tested the model using a bulk RNA sequencing (RNaseq) platform. Based on our power calculation, 14 blood samples are required to achieve a power of 0.9 (see Materials and methods for details). To ensure we achieved this, we collected and performed bulk RNaseq on a total of 32 blood samples (15 responders and 17 nonresponders) from melanoma patients (dataset 3; Fig. 8 C) and four lung cancer patients (dataset 4; Fig. S5 E) treated with immunotherapy. Remarkably, combining analysis of validation sets 1–3, our model achieves 88% accuracy (note, in the melanoma, it achieved accuracy of 89% in set 2 and 88% in set 3, Fig. 8 C). Surprisingly, it achieves 100% accuracy in a small cohort of lung cancer patients (not included in the overall accuracy calculation of melanoma; see Fig. S5 E). These results validate our GEP model as an efficient predictor of immunotherapy responses; this algorithm is effective whether the source material is from the tumor or obtained less invasively using blood-based approaches. In all of the validation cohorts, our model assigned significantly higher NRSs to nonresponders than responders (t test, P = 1.56 × 10⁻⁸; Fig. 8 C) and effectively discriminated them with an area under the receiver-operating characteristic curve (ROC-AUC) of 0.90 (Fig. 8 D). Given that blood-based approaches have greater potential for utility in oncology, as they do not require invasive acquisition of CD8 TILs, we termed our prediction platform, for simplicity, the noninvasive circulating T cells model (NiCir). Collectively, our validation data demonstrate that NiCir has high predictive accuracy in TILs or PBLs and can be applied to multiplatform datasets. Of note, given NiCir’s ability to predict ICI responses from peripheral blood CD8 T cells, it can now be used as a noninvasive tool to predict immunotherapy responses in the clinical setting.

For a number of years, our group has defined tumor-induced dysfunctional changes in infiltrating T cells (Maybruck et al., 2017; Montes et al., 2008; Pfannenstiel et al., 2019; Pfannenstiel et al., 2018; Zhang et al., 2011). Our recent work showed that in cancers like melanoma, a large percentage of tumor-associated CD8 T cells express high levels of checkpoints like PD-1, and these same cells can actively suppress other immune effectors (Pfannenstiel et al., 2019). Thus, expressing classical exhaustion markers (e.g., checkpoints) does not necessarily define a CD8 T cell as anergic and can be associated with a form of active dysfunction; little is known about this phenomenon outside of nonmalignant pathologies (Flippe et al., 2019; Van Kaer, 2010). Instead, targeting suppressor cells like CD4 T regulatory cells, myeloid-derived suppressor cells, and M2-macrophages remains the main focus in the majority of cancer studies (Klug et al., 2013; Liu et al., 2018; Ugel et al., 2015). CD8 T cells, a critical player in immunotherapy effectiveness, have been extensively studied, but mainly in terms of their exhaustiveness or inefficiency of tumor infiltration (Huang et al., 2017; Jansen et al., 2019; Mariathasan et al., 2018). Thus, if additional TIL CD8 T cell heterogeneity exists within the population targeted by checkpoint inhibitors, then this should prompt additional investigation to improve this line of immunotherapeutics.

Interestingly, CD8⁺PD-1⁺ T cells in the periphery of cancer patients were found to contain overlapping TCRs with CD8 T cells in the tumor and are also the most likely to yield (after expansion and adoptive transfer) an anti-tumor response (Gros et al., 2016). The NCI surgery branch that made this discovery surmised that the peripheral and intratumoral TCR sharing CD8⁺PD-1⁺ T cells could be interrelated (Gros et al., 2016). If this peripheral blood–intratumoral relationship exists, then it may provide insight into the same tumor-derived CD8 T cells expressing checkpoints that we have been studying and were found to be affected by ICIs (Pfannenstiel et al., 2019). We hypothesized that we can leverage this PBL–TIL relationship to gauge the status of immunotherapy resistance and learn from the related biochemical pathways within CD8 T cell subpopulations.

Using scRNaseq, we discovered that only three unique transcriptionally shared/overlapping clusters in PBL-CD8 and TIL-CD8 T cells exist, two of which we determined spawn from naive CD8 T cells through specific reprogramming, as illustrated in our trajectory analysis; we defined them above as inert dormant and CD8⁺ T_OXPHOS cells. Unlike classic exhausted CD8 T cells, CD8⁺ T_OXPHOS cells in particular appear to contain many of the hallmarks of exhausted CD8 T cells that we previously studied, including checkpoint expression (e.g., PD-1), but they also expressed cytotoxic markers (e.g., IFN-γ and GZMB), corroborating recent reports illustrating that exhausted T cells are not homogeneous but a gradual transition (van der Leun et al., 2020; Yost et al., 2019). This is not completely surprising, as we and others had already identified that suppressor CD8 T cells (which we later found express multiple checkpoints) regulate other cells via expression of IFN-γ (Hidalgo and Halloran, 2002; Mele et al., 2003; Montes et al., 2008; Robb and Hill, 2012; Robb et al., 2011). In fact, studies show that IFN-γ, depending on the conditions, can act as an effector or suppressor immune mediator, and the ramifications of our findings are under investigation (Lee and Ashkar, 2018). These CD8⁺ T_OXPHOS cells also have up-regulated metabolic pathways, consistent with the fact that dysfunctional CD8 T cells are functionally diverse (van der Leun et al., 2020). Our discovery of metabolically active CD8⁺ T_OXPHOS cells supports recent work that CD8 T cells in TILs with PD-1 expression have two distinct states, low and high metabolism (Hartmann et al., 2021). Interestingly, in their study, which used CYTOF (cytometry by time of flight), Hartmann et al. (2021) found that cells that, based on our work, would likely be CD8⁺ T_OXPHOS cells interface directly with the tumor, whereas the more classical exhausted CD8 T cells were found in the peripheral areas of the tumor, thus supporting a unique role for CD8⁺ T_OXPHOS cells in tumor immunology.

Recent scRNAseq studies have shown that TIL CD8 T cells with hallmarks of exhaustion (e.g., expression of PD-1) contain subpopulations that are actively cycling and proliferating (Li et al., 2019). However, few studies have focused on identifying the mechanistic drivers of these active dysfunctional cells. Most reports study lymphocytes in general and not solely lymphocytes on CD8 T cells (Li et al., 2019; Sade-Feldman et al., 2018b; Tirosh et al., 2016; van der Leun et al., 2020). Despite recent attempts by using scRNaseq to investigate and identify CD8 T cells in the peripheral blood and TILs, only Wu et al. and Fairfax et al. have reported the identification of shared populations (and only within cytotoxic CD8 T cells; they did not identify the CD8⁺ T_OXPHOS cells found in our study; Fairfax et al., 2020; Wu et al., 2020). By focusing on enriched CD8 T cells in both TILs and PBLs, with sufficient cell numbers from individualized paired patient samples, we had greater sensitivity and resolution compared with most previous studies, allowing us to better analyze these underexplored CD8 subpopulations.

Recent studies have explored the role of TIL T cell metabolism in regulating immunotherapy, where metabolic properties of T cells exert an essential role in anti-tumor immunity (Guerra et al., 2020; Hartmann et al., 2021; Leone and Emens, 2018; Leone et al., 2019; Scharping et al., 2021). The metabolic state of peripheral blood CD8 T cells is an even more immature field (Hatae et al., 2020; Varanasi et al., 2020; Voss et al., 2021). Remarkably, our data indicate that higher metabolic activity, including up-regulation of OXPHOS and glycolysis within CD8⁺ T_OXPHOS cells in TILs and PBLs, is associated with ICI resistance. To translate this work for both potential immunotherapeutic targeting and cellular validation, we identified uniquely high expression of the cell surface molecules CD38 and CD39 on CD8⁺ T_OXPHOS cells, whose inhibitors are being assessed in current ongoing clinical trials targeting solid cancers like melanoma (Allard et al., 2017; Lokhorst et al., 2015; Perrot et al., 2019; van der Leun et al., 2020). Both CD38 and CD39 are ectonucleotidases that play a role in NAD⁺/ATP regulation, which is connected to several metabolic pathways (Kishton et al., 2017; Sukumar et al., 2017; Thommen and Schumacher, 2018). Based on these findings, it is possible that these ectonucleotidases are associated with the high metabolic state that we observed in CD8⁺ T_OXPHOS cells.

Another group demonstrated that activation of PBL Ki-67⁺PD-1⁺CD8⁺ T cells (where Ki-67 represents increased proliferation capacity and assumedly increased metabolism) before anti-PD-1 therapy correlates with poor ICI response (Huang et al., 2017). The same group later reported that patients with sustained high levels of Ki-67 in PD-1⁺CD8⁺ T cells in blood have better outcomes after a single dose of neoadjuvant anti-PD-1 therapy (Huang et al., 2019). To understand this discrepancy, they hypothesized that these proliferating T cells likely convert to a more reinvigorated and functional state after treatment with anti-PD-1 therapy (Huang et al., 2019). Whether these cells in PBLs represent CD8⁺ T_OXPHOS cells or other populations that can contain these makers is unknown. However, if these Ki-67⁺PD-1⁺CD8⁺ exhausted T cells are CD8⁺ T_OXPHOS cells, then it would imply that patients whose CD8⁺ T_OXPHOS cells become reinvigorated by anti-PD-1 therapy have better outcomes, and identifying why this does not occur ubiquitously could create strategies to improve survival in all melanoma patients. Our study, which illuminates many of the mechanisms found to be activated in CD8⁺ T_OXPHOS cells, offers many targets that may become new therapeutic avenues.

Anti-PD-1–based ICI has revolutionized the care of melanoma; unfortunately, a correlative assay that reliably aids clinicians in predicting who will respond continues to be an unmet need (Halim, 2015; Masucci et al., 2016). Currently, only three predictive biomarkers were approved by the US Food and Drug Administration for ICI therapies for any cancer: protein expression of PD-L1 by immunohistochemistry, microsatellite instability status by PCR phenotyping, and tumor mutation burden (Boyiadzis et al., 2018; Davis and Patel, 2019; Prasad and Addeo, 2020). In recent years, most biomarker exploration has been focused on searching for methods that use minimally invasive approaches (Hofman et al., 2019; Nixon et al., 2019). To date, despite a large amount of effort, whether it be with tumor tissue or PBLs, and using various “omic” technologies or immunohistochemistry-based studies, there is no standard correlative assay used for ICI therapy in the treatment of melanoma (Auslander et al., 2018; Fattore et al., 2021; Havel et al., 2019).

Based on our discovery set, we developed a GEP-coined NiCir that we validated with a published dataset (Yost et al., 2019) and four independent datasets, including from TILs and PBLs from our own patients. NiCir had an average accuracy of 88% in the validation datasets, indicating good prediction value across a spectrum of different sequencing platforms. NiCir showed excellent predictive power, with an area under the curve (AUC) of 0.90, as compared with several previously reported immunotherapy response predictions; for example, IMPRES (immunopredictive score) achieved an AUC of 0.83 for overall accuracy (Auslander et al., 2018); a radiographical characteristic in non–small cell lung cancer patients using artificial intelligence with an AUC of up to 0.76 (Trebeschi et al., 2019); and tumor mutational burden in non–small cell lung cancer predicting ICI response with an AUC range of 0.554–0.755 (Wang et al., 2019). We plan to translate this work into an assay that can be easily applied in a clinical setting for more timely decision making.

In summary, our results increase the field’s knowledge of the heterogeneity of dysfunctional CD8 T cells in cancer patients. Our work implicates the potential of manipulating the immunometabolism of T cell populations like CD8⁺ T_OXPHOS cells to regulate and improve anti-tumor immunity. Moreover, the presence of CD8⁺ T_OXPHOS cells in particular can be exploited even with blood-based approaches for a rapid predictive ICI efficacy assay to improve clinical outcomes in melanoma patients.

All human tissue was obtained at the Cleveland Clinic under a protocol approved by Cleveland Clinic’s institutional review board, and written informed consent was obtained from each patient. All patient samples used for scRNaseq were immunotherapy nonresponders. Peripheral blood samples of immunotherapy responders were used as additional validation cohort for the predictive model (NiCir).

Peripheral blood mononuclear cells were purified from buffy coats by centrifugation over a Ficoll-Hypaque gradient according to the manufacturer’s protocol (GE Healthcare). Matched PBLs and tumor specimens were obtained from patients with melanoma. PBLs were obtained by venipuncture and isolated by centrifugation over a Ficoll-Hypaque gradient. After surgical resection, tumor specimens were rinsed with antibiotic-containing media and minced with crossed scalpels under sterile conditions. Enzymatic digestion was then used to dissociate tumor tissue using 1,500 U/ml collagenase IV (Gibco/Life Technologies), 1,000 U/ml hyaluronidase (Sigma), and 0.05 mU/ml DNase IV (Gibco) in RPMI for 1 h at 37°C followed by mechanical agitation. The resulting single-cell suspensions were separated from debris by centrifugation over a Ficoll-Hypaque gradient followed by being immediately gated and sorted for LIVE > CD3⁺ > CD8⁺ using LIVE/DEAD (Invitrogen; #L10119) CD3 (BioLegend; #300406) and CD8 (BioLegend #300912). The negative population was dumped using dump channel PE-Cy7 (YG) for CD4 (BioLegend; #300454), CD19 (BioLegend; #302215), and CD56 (BioLegend; #302215). On average, sorted cell counts were between 18,605 and 297,000, with an average of 91% efficiency for scRNaseq. Sorted cells were analyzed with flow cytometry to ensure over 90% purity (data not shown).

The oxygen consumption and glucose dependency of healthy donor PBLs and melanoma TILs were determined using the Seahorse XF Mito Fuel Flex Test (Agilent; 103260–100) with the XFe96 Bioanalyzer (Agilent) according to the manufacturer’s instructions. Briefly, enriched CD8 T cells of PBLs and sorted CD8⁺PD-1⁺ T cell TILs were plated into Agilent Seahorse XF96 microplate in triplicate (75,000 cells/well) overnight, followed by 1-h incubation at 37°C in a non-CO₂ incubator on the day of the assay. The pyruvate inhibitor UK5099 (2 mM) was injected into the microplates to determine the glucose oxidation dependency. Baseline OCR was monitored at 24 min, followed by sequential pyruvate inhibitor UK5099 (2 mM) injections, and the OCR readings were recorded for a duration of 2 h at 8-min intervals. Glucose OCR and dependency were calculated using the precent dependency equation based on the instructions, and graphs were plotted with GraphPad Prism software.

The ATP bioluminescence assay was used to quantify ATP levels. TIL and PBL samples were labeled with CD8 microbeads (Miltenyi Biotec; #130-045-201, MACS) and then positively selected for CD8 T cells by passing them through the LS column (Miltenyi Biotec; #130-042-401, MACS) according to the manufacturer’s instructions. ATP levels of the positive-selection CD8 T cells were determined using a Luminescent ATP Detection Assay kit (Abcam; ab113849) following the manufacturer’s instructions. Briefly, CD8 T cells were incubated with detergent for 5 min, followed by incubation with substrate solution for 10 min at room temperature in the dark with agitation. Then, luminescence was quantified on a microplate reader (Molecular Devices; SpectraMax M2). The luminescence values of ATP standards were determined in a similar manner. ATP concentration was calculated according to an ATP-standard curve.

Human-specific antibodies and key reagents used are tabulated in Table S7. Cell viability was determined using LIVE/DEAD Fixable DEAD Cells Stain Kit (Thermo Fisher Scientific). Human PBLs (naive, n = 4; responders, n = 5; nonresponders, n = 7) and TILs (naive, n = 5; nonresponders, n = 6) from melanoma patients were stained with CD3, CD8, PD-1, CD38, CD39 (all from BioLegend), TMRM (Thermo Fisher Scientific), and MitoTracker Green (Thermo Fisher Scientific; for details, refer to Table S8). Cells were incubated for 20 min at 37°C with TMRM and MitoTracker Green according to the manufacturer’s instructions. Next, the cells were stained with antibodies for 13 min at 4°C in FACS buffer, followed by flow cytometry analysis. Compensation controls and fluorescence minus one controls were used to determined cell populations. Flow cytometry was performed using a BD LSRFortessa and analyzed with FlowJo software, and normalized mean fluorescence intensity (MFI) graphs were plotted with GraphPad Prism software.

All cells were resuspended in DPBS with 0.04% BSA, and immediately processed for scRNaseq as follows. Cell count and viability were determined using trypan blue on a Countess FL II, and ∼12,000 cells were loaded for capture onto the Chromium system using the v2 single-cell reagent kit according to the manufacturer’s protocol (10x Genomics). Following capture and lysis, cDNA was synthesized and amplified (12 cycles) as per the manufacturer’s protocol (10x Genomics). The amplified cDNA from each channel of the Chromium system was used to construct an Illumina sequencing library and sequenced on an Illumina NovaSeq with 150-cycle sequencing (asymmetric reads per 10x Genomics). Illumina basecall files (*.bcl) were converted to FASTQs using CellRanger v3.0, which uses bcl2fastq for FASTQ file generation. FASTQ files were then aligned to GRCh38 human reference genome and transcriptome using the CellRanger v3.0 software pipeline with default parameters as reported previously (Larkin et al., 2019); this demultiplexes the samples and generates a gene versus cell expression matrix based on the barcodes and assigns unique molecular identifiers (UMIs) that enable determination of the individual cell from which the RNA molecule originated.

PBLs were sorted for CD8 T cells using a selective gating strategy based on staining with the following monoclonal antibodies: negative gating for staining with LIVE/DEAD Aqua (Invitrogen; #L34957), CD14 (BD Biosciences; #564444), CD19 (BioLegend; #302242), CD56 (BD Biosciences; #564058), and CD4 (Invitrogen; #Q10008) to exclude all dead cells and monocytes, B cells, natural killer cells, and CD4 T cells, followed by positive gating on CD3 (BD Biosciences; #557943) and CD8 (Invitrogen; #MHCD0817) to selectively sort 15,000 cells from a pure CD8 population. Purity tests for sorting were performed routinely. RNA was purified from CD8 T cells using RNeasy Micro Kits (Qiagen). Libraries of T cell RNA samples were prepared using Illumina’s Nextera XT DNA Library Prep Kit and sequenced on Illumina NovaSeq 6000. The libraries were generated following the Library Prep Kit manual. The libraries were pooled and prepared for sequencing following the “Standard Normalization Method” within the Illumina Novaseq6000 System Denature and Dilute Libraries Guide. The BCL files were demultiplexed using Illumina’s bcl2fastq v2.20.0.422. Reads in FASTQ files were aligned to GRCh38 human reference genome and transcriptome with Hisat2. Tags per million reads were calculated from aligned reads using StringTie.

t-SNE dimensionality reduction of cells based on whole transcriptomes was generated by 10x Genomics Cell Ranger pipeline (version 3.0), as reported previously (Larkin et al., 2019). Dimensionality of gene-barcode matrices was first reduced to 10 principal components using principal-component analysis (PCA). PCA-reduced data were further reduced to 2D space using the t-SNE method and visualized in the Loupe Cell Browser (10× Genomics) and/or R. Graph-based clustering of cells was conducted in the PCA space; a sparse nearest-neighbor graph of the cells was built first, and Louvain modularity optimization was then applied. The number of nearest neighbors was logarithmically in accordance with the number of cells. In the last step, repeated cycles of hierarchical clustering and merging of cluster pairs that had no significant differential expression were performed until no more cluster pairs could merge.

3D t-SNEs of cells were calculated using the “cellranger reanalyze” command with modified parameters and visualized using R plotly package (https://plot.ly).

Gene differential expression analyses of each cluster were conducted by cellranger (10x Genomics). The log₂ fold change of a certain gene’s expression (UMIs) in one cluster compared with all other clusters, and the corresponding adjusted P values, were calculated for each cluster.

To compare similarity between the three shared clusters (2, 6, and 15) and other clusters, Pearson correlations were calculated between each pair using the clusters’ average gene expressions (mean UMIs of genes). The results were visualized using Circos Table Viewer (http://mkweb.bcgsc.ca/tableviewer/), with ribbons connecting two clusters representing significant similarity (50th percentile of all Pearson distances) between them (Krzywinski et al., 2009).

P values of Pearson correlation were probabilities of T > t or < −t, where T follows a t distribution with n − 2 degrees of freedom.

The P values were adjusted with FDR correction, and the genes of significant adjusted P values (< 0.05) were then ranked by their Pearson correlation coefficient r calculated in the three shared clusters from the most significant to the least significant correlation.

Signaling pathway enrichment and upstream regulator prediction were analyzed using Ingenuity Pathway Analysis (IPA; Qiagen), unless indicated specifically.

Single-cell trajectories of cells were built using the Monocle package (version 2.8.0) as previously reported (Larkin et al., 2019). Briefly, whole-transcriptome trajectories of T cells were built using 9,080 genes that were expressed (UMI >1) in more than 1% of the total T cells. Dimensionality reduction was conducted using the DDRTree method, and the minimum spanning tree was plotted using the plot_cell_trajectory function.

Expression levels of the signature genes along trajectories (i.e., from the cluster 2 end to the cluster 6 or cluster 15 end) were visualized using the plot_genes_branched_heatmap function; on the generated heatmap, gene expression levels were smoothened using the VGAM package, rescaled to a −3 to 3 range, and hierarchically clustered.

The prognosis model to predict response of patients to ICI therapy was built using PD-1 and the top 20 genes of the PD-1 coexpression genes (21 predictor genes in total) that were mostly correlated (most significant P values) with PD-1. Responders/nonresponders were defined following the RECIST criteria (i.e., complete response and partial response for responders or stable disease and progressive disease for nonresponders; Eisenhauer et al., 2009).

The Gene Expression Omnibus (GEO) database was surveyed for all available data that were provided with both T cell transcriptomics and matching patient treatment response records to immune checkpoint therapy. A total of two published datasets, GSE120575 (Sade-Feldman et al., 2018b) and GSE123813 (Yost et al., 2019), were returned from this survey, and both were used for this study. We randomly assigned GSE120575 as training data, and GSE123813, along with our two additional datasets, were used as independent validation sets.

The model with β₀and β_i determined from the regression test were applied to calculate P (probability) for each individual CD8 T cells from the dataset and returned with a predicted status (1, responder; 0, nonresponder). Subsequently, the proportions of all calculated CD8 T cells in one sample were calculated and presented for each patient with the scored output using the abovementioned formula.

The model was then further validated in two additional datasets. Dataset 1 was from a previously published study (GSE123813; Yost et al., 2019), in which CD8 tumor-infiltrating T cells were collected before ICI therapy from three and six patients who responded or did not respond to anti-PD-1 therapy, respectively. The response status of each T cell was predicted, and the percentage of predicted responding/nonresponding T cells in each patient was calculated. NRSs of patients were calculated by scaling the percentage of predicted nonresponding T cells to a 0∼10 range. Patients that had NRSs higher than 5.0 (more than 50% of T cells are predicted to respond) were considered responders to ICI therapy; otherwise, patients were predicted as nonresponders. Dataset 2 was collected in the present study, including peripheral CD8 T cells of 32 melanoma patients (15 responders and 17 nonresponders) profiled by RNaseq. Patients’ gene expression levels were input into out prediction model. For CD8 T cell scRNaseq data, the model predicted status of each T cells, and the ratio of nonresponding T cells in a patient was rescaled to 0∼10 as the patient’s NRS. For bulk RNaseq data, the output values of our logistic regression model (ranging 0∼1) were rescaled to 0∼10 as patients’ NRSs. Patients that had NRSs higher than 5.0 were considered responders to ICI therapy. The P values comparing the NRSs between true responders and true nonresponders were calculated by t test. A P value < 0.05 is considered significant.

Unless otherwise stated, P values of gene differential expression were determined by two-tailed Welch’s t test. P values of enrichment of pathways, upstream regulators, and gene ontology terms were generated by the corresponding bioinformatics tools. All statistics calculations were conducted using R unless otherwise stated.

To determine the size of samples for validating our NiCir model, we conducted power calculation to determine the required sample size for NiCir using the rocr package in R. Power calculation suggested that in a balanced feature model (kappa = 1), with a power of 0.9 and a significance level of 0.05 for the logistic regression prediction model, a proper sample size should be at least 14 samples (seven responders and seven nonresponders), as shown in this case of NiCir for a AUC of 0.9 (as indicated in our training set, Fig. 8 B), a proper sample size should be at least 14 samples (seven responders and seven nonresponders). We validated the NiCir model using nine previously published samples (Yost et al., 2019) and 32 samples collected in the present study (Fig. 8 C). We collected samples from both responders and nonresponders to ensure a proper kappa = 1 for the model evaluation.

Fig. S1 shows scRNaseq data of melanoma patients’ PBLs and TILs. Fig. S2 shows clusters of CD8 T cell profiles from melanoma patients. Fig. S3 shows the characteristics of three clusters containing equal numbers of CD8 PBLs and TILs from each patient. Fig. S4 shows the characterization of CD8 T cells. Fig. S5 shows a heatmap of CD8 T cell state depicted by signature genes, a heatmap of the top 1,000 PD-1–coexpressing genes, MFI of TMRM^total and TMRM^low in both CD8 PBLs and TILs, and the performances of NiCir’s prediction in the validation dataset (GSE152590) of PBL samples from lung cancer patients collected at Cleveland Clinic core. Table S1 shows the patient characteristics. Table S2 shows scRNaseq data of CD8 T cells from melanoma patients. Table S3 lists the pathways in CD8 PBLs and TILs from melanoma patients. Table S4 lists genes enriched in pathways in CD8 PBLs and TILs. Table S5 lists the signature genes for each cluster. Table S6 lists the enriched pathways in clusters 2, 6, and 15. Table S7 lists the genes enriched in pathways in clusters 2, 6, and 15. Table S8 lists the antibodies and reagents used for flow cytometry.

T cell scRNaseq and RNaseq data generated during this study are available through National Center for Biotechnology Information GEO accession numbers GSE138720, GSE153098, and GSE171256. The other T cell scRNaseq data used in this study are from previously published studies and are available through GEO accession numbers GSE120575 (Sade-Feldman et al., 2018b) and GSE123813 (Yost et al., 2019).

We are grateful to Drs. Andrew Henry Lichtman (Brigham and Women’s Hospital), Pramod K. Srivastava (UConn Health), Stephen Safe (Texas A&M University), and Xiaoxia Li (Lerner Research Institute, Cleveland Clinic) for valuable advice and comments on this project. We are very thankful for Dr. Tomas Radivoyevitch for his contribution in power calculation. We thank the Genomics Core at the Lerner Research Institute of Cleveland Clinic for 10× Genomics single-cell library preparation.

This research is supported by Cleveland Clinic internal funding.

Author contributions: B.R. Gastman and B. Zhou conceived and designed the project; B.R. Gastman, P. Funchain, J.M. Song, C.M. Diaz-Montero, Y.P. Phoon, Y.F. Tian, C. Cameron, and S. Thapaliya collected patient samples and prepared or designed cell sorting for scRNAseq; C. Li performed computational processing of scRNAseq and RNAseq data. C. Li and B. Zhou developed analytic strategies and the prediction model; Y.P. Phoon performed flow cytometry; A. Thongkum performed the ATP assay; L. Qu, A.J. Matz, B. Tamilselvan, J.B. Golden, and M. Cartwright prepared library and conducted RNAseq of the validation dataset; C. Li, B. Zhou, B.R. Gastman, A. Vella, A. Rodriguez, K. Karlinsey, L. Qu, S. Thapaliya, Y.P. Phoon, M. Cameron, and C. Cameron interpreted data and conducted omic or functional analyses; B. Zhou, C. Li, K. Karlinsey, B.R. Gastman, and Y.P. Phoon wrote the manuscript; and C. Bonin, A. Vella, A. Rodriguez, A. Menoret, A. Rodriguez, K. Karlinsey, L. Qu, S. Thapaliya, M. Cameron, C. Cameron, Y.P. Phoon, B. Zhou, and B.R. Gastman revised the manuscript.