PD-1 marks dysfunctional regulatory T cells in malignant gliomas

The introduction of immunotherapies that target immune checkpoint receptors expressed on T cells represents a new paradigm for treating cancer (1). Developing a clear immunological understanding of how these therapies work — particularly with respect to the phenotype and function of the cells they target — should enable further improvements of these therapies in the clinic. Two approved antibodies target coinhibitory receptors expressed on T cells — cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4) (2, 3) and programmed cell death protein 1 (PD-1) (4, 5). These receptors impart negative feedback through signaling pathways that modify the thresholds for activation and functional responses (6). The dramatic clinical responses following treatment with anti-PD-1 monoclonal blocking antibodies in certain malignant tumors along with the significant expression of PD-1 on tumor-infiltrating CD4 T cells suggest a role for these cells and this pathway in modulating tumor immune responses (7).

CD4⁺CD25^hiFoxP3⁺ Tregs express specific coinhibitory and costimulatory receptors involved in signaling pathways that shape effector functions, including CTLA-4 (8), T cell immunoglobulin and mucin domain-containing 3 (TIM-3) (9), and T cell immunoreceptor with immunoglobulin and ITIM domains (TIGIT) (10, 11). Tregs play a central role in the regulation of immune responses, becoming dysfunctional in autoimmunity (12, 13) and arguably hyperfunctional in the context of immune suppression in tumors (14–17). PD-1 is expressed upon T cell activation, providing negative feedback on the effector functions of T cells during inflammation, and therefore suppresses autoimmune responses (18). The role of PD-1 expression on CD8 T cells in chronic infection, including HIV and hepatitis C infection, has been extensively studied in both mice and humans (19–22). Following chronic antigenic stimulation during viral infections, increased expression of PD-1 on CD8⁺ T effector cells has been associated with an “exhausted” phenotype, distinguished by progressive loss of effector functions, including proliferation and cytokine production (23). The molecular signature of CD4 cell exhaustion in mice was also recently described and shows similarities to CD8 T cell exhaustion (24). However, the key molecular and functional characteristics of circulating and tumor-infiltrating Tregs expressing PD-1 is not known.

PD-1 is expressed by a significant percentage of tumor-infiltrating lymphocytes, including Tregs in the context of many cancers, such as lung (25) and brain (26). The function of PD-1^hi Tregs in tumors, however, is also not clear. Signaling through the PD-1 pathway favors the induction of Tregs in the periphery in mice (27). In contrast, liver-infiltrating Tregs from patients with hepatitis C displayed impaired in vivo expansion and suppression activity through the interaction between PD-1 and PD-L1 (28). Given the importance of PD-1 in regulating immune responses in tissues and its paradoxical role as both an activation and an exhaustion marker, understanding the function of PD-1-expressing Tregs and their role in regulating tumor immune responses remains of critical clinical importance.

Here, we examine the functional and molecular features of Tregs from both healthy subjects and from tissue of patients with malignant gliomas, using a combination of ex vivo functional assays, RNA sequencing, and cytometry by time of flight (CyTOF). In healthy subjects, circulating Tregs expressing high amounts of PD-1 were impaired in their ability to suppress CD4⁺ effector T cells, produced high quantities of IFN-γ, showed partial demethylation at the Treg-specific demethylated region (TSDR) locus, and had eroded telomeres. Interrogation of Tregs from malignant gliomas demonstrated high PD-1 expression, with loss of suppressor function associated with secretion of IFN-γ, reflecting the phenotype of circulating PD-1^hi Tregs in healthy subjects. Transcriptional profiling of tumor-resident Tregs revealed several key genes that are coexpressed with PD-1, including ICOS, TIM-3, lymphocyte-activation gene 3 (LAG-3), and STAT4. Treatment of patients with anti-PD-1–blocking mAb further increased the secretion of IFN-γ by Tregs, perhaps redirecting less differentiated Treg populations toward a more dysfunctional state and further enhancing tumor clearance. Thus, high PD-1 expression on human Tregs identifies a population of dysfunctional, exhausted Tregs secreting IFN-γ that exist in healthy individuals and are enriched in malignant glioma infiltrates as they attempt to modulate the inflammatory tumor immune responses.

To assess the functional state of circulating PD-1^hi Tregs of healthy individuals, CD25^hiCD127^lo Tregs were sorted into PD-1^hi and PD-1^– populations. High expression of PD-1 was found on approximately 5% of both CD4 effector T cells and Tregs (Supplemental Figure 1, A and B; supplemental material available online with this article; doi:10.1172/jci.insight.85935DS1). PD-1^hi Tregs expressed similarly high levels of FoxP3 as PD-1^– Tregs, suggesting their purity as Tregs (Supplemental Figure 1, C and D). The suppressive function of PD-1^hi Tregs was assessed using standard ex vivo suppression assays. Sorted PD-1^– and PD-1^hi Tregs (CD4⁺CD25^hiCD127^–) and T effector cells (CD4⁺CD25^loCD127⁺) were cocultured and stimulated with αCD3/αCD28/αCD2 for 4 days. The percentage of T effector cells entering the cell cycle with or without Tregs was then assessed by dilution of CSFE (Figure 1, A and B). Strikingly, PD-1^hi Tregs exhibited a significant reduction in suppressor function.

Figure 1

Programmed cell death protein 1 (PD-1)^hi Tregs do not suppress and produce IFN-γ. (A) Representative suppression of T effector proliferation by Tregs sorted on programmed cell death protein 1 (PD-1)^hi expression. Tregs and T effector cells were stimulated for 4 days with αCD3/αCD28/αCD2. (B) Percentage (mean ± SD) of suppression of 6 individuals (Student’s paired t test). (C) Representative production of IFN-γ of sorted Tregs stimulated with aCD3/aCD28/aCD2 in presence of IL-2 for 4 days and restimulated with PMA and ionomycin in the presence of brefeldin A for 4 hours. (D) Percentage (mean ± SD) of cytokine-producing Tregs of 6 individuals (Student’s paired t test). (E) Quantification of phosphorylated FoxO-1 (Ser319). Representative histogram from one donor and median fluorescence intensity (MFI, mean ± SD) for 6 donors (Student’s paired t test). (F) Proliferation of sorted Tregs stimulated with aCD3/aCD28/aCD2 in presence of IL-2 for 4 days. One donor representative of 4.

Dysfunctional human Tregs in subjects with autoimmune disease and those induced with the high concentrations of sodium chloride found in tissue secrete the proinflammatory cytokine IFN-γ (29–31). Therefore, we assessed whether PD-1^hi Tregs also produced IFN-γ by culturing sorted Treg populations for 4 days, followed by stimulation with PMA and ionomycin in the presence of brefeldin A. PD-1^hi Tregs produced significantly more IFN-γ than PD-1^– Tregs (Figure 1, C and D). In Tregs, resistance to production of IFN-γ upon T cell receptor stimulation is achieved through inhibition of the PI3K/Akt pathway, resulting in nuclear localization of unphosphorylated FoxO-1 and transcriptional repression of the ifng gene (32, 33). To examine whether IFN-γ secretion by PD-1^hi Tregs is a hallmark of increased PI3K/Akt activity, we assessed the ex vivo phosphorylation of FoxO-1 in sorted PD-1^– and PD-1^hi Tregs (Figure 1E). Indeed, we observed increased phosphorylation of FoxO-1 in PD-1^hi Tregs, indicative of its exclusion from the nucleus, suggesting that the PI3K/Akt pathway is hyperactive in PD-1^hi Tregs. Blocking IFN-γ has been shown to partially restore suppressive activity of autoimmune IFN-γ⁺ Tregs, but this perturbation did not improve the suppressive function of PD-1^hi Tregs (Supplemental Figure 2). These results suggest that dysfunctional PD-1^hi Tregs produce increased IFN-γ, potentially through hyperactivity of the PI3K/Akt pathway, though the production of IFN-γ does not fully explain loss of functional suppression.

Tregs replicate extensively in vivo, and a reduction in this activity has been associated with dysfunctional Tregs in autoimmune disease (34). In order to assess how impaired suppressor function of PD-1^hi Tregs related to entry into cell cycle, we stimulated sorted PD-1^– and PD-1^hi Tregs with αCD3/αCD28/αCD2 for 4 days. PD-1^– Tregs displayed limited proliferative activity, consistent with the well-documented in vitro partial anergy of Tregs, while the PD-1^hi population displayed a trend to a reduction in the frequency of cells undergoing divisions (Figure 1F and Supplemental Figure 3). Moreover, while the PD-1^– compartment proliferation correlated with increase of FoxP3 expression, PD-1^hi Tregs failed to upregulate, or even lost, FoxP3 expression upon entering cell cycle, despite comparable levels of FoxP3 prior to stimulation (Supplemental Figure 1D).

PD-1 exerts its coinhibitory activity by preventing CD28-mediated activation of PI3K, resulting in reduced activity of Akt (35). Inhibition of the Akt/mTOR axis is critical for Treg development (36) and for the subsequent stability of the suppressive function of Tregs, including repression of IFN-γ (32, 37). Therefore, the dysfunction of PD-1^hi Tregs may not result from direct PD-1–mediated signaling. Indeed, antibodies blocking PD-1, its ligand PD-L1, or both together during suppression assays failed to rescue the regulatory activity of PD-1^hi Tregs (Supplemental Figure 4). Taken together, these data indicate that high expression of PD-1 delineates dysfunctional Tregs.

To further characterize the molecular phenotypes of circulating PD-1^hi Tregs, we performed transcriptional profiling of sorted PD-1^hi, PD-1^–, and total Tregs isolated from healthy individuals by RNA sequencing. T cells were sequenced to a depth of approximately 10 million reads per sample (Supplemental Table 1), and all samples from healthy controls were included for subsequent analysis based on their high-quality sequencing metrics (Supplemental Figure 5A). Biological replicates showed consistent expression of several housekeeping genes and CD4 T cell–specific genes, including CD4, ANK3, and FoxP3 (Supplemental Figure 5B). A total of 791 differentially expressed genes, including TIM3 and LAG3, distinguished PD-1^hi Tregs from PD-1^– Tregs (Figure 2A and Supplemental Table 2). Functional classification of differentially expressed genes that were either upregulated or downregulated identified 18 gene modules in PD-1^hi Tregs (Supplemental Table 3). One gene module included the coinhibitory molecules LAG3, PD-1, and TIM3 as well as several surface receptors, such as CXCR3, CCR6, CD226, and CCR5 (Figure 2B and Supplemental Figure 6). A total of 6 gene modules were downregulated in PD-1^hi Tregs, including a large group of DNA-binding proteins (Supplemental Table 4). We confirmed differential expression of several of these surface receptors, including CCR6, CXCR3, CD226, and LAG-3, by flow cytometry (Supplemental Figure 7). Additionally, a recent report identified CD226 as the marker of IFN-γ production by human Tregs (11). Taken together, these data show that PD-1 expression marks a distinct population of human Tregs that coexpress several surface receptors, including the coinhibitory receptors LAG-3 and TIM-3.

Figure 2

Programmed cell death protein 1 (PD-1)^hi Tregs are transcriptionally distinct from PD-1^– Tregs. (A) Volcano plot showing differentially expressed genes identified through RNA sequencing (n = 3 PD-1^hi and 3 PD-1^–). (B) DAVID functional classification was performed on differentially expressed genes to identify functional gene modules. A summary of one functional gene module enriched for surface markers identified in PD-1^hi Tregs is shown as a heatmap. Full heatmaps are shown in Supplemental Figure 6. (C) Enrichment of gene sets for each indicated category. Gene sets were broadly classified by function, and representative enriched gene sets are shown. Full gene lists are shown in Supplemental Table 5.

We next performed gene set enrichment analysis (GSEA) to better understand the biological functions enriched in circulating PD-1^hi Tregs. Several gene sets related to cell cycle and proliferation, metabolic function, and the immune system showed significant enrichment in PD-1^hi Tregs (FDR < 0.05) (Figure 2C and Supplemental Table 5). The top-scoring gene sets (FDR < 0.05) included genes associated with CTLA-4, PD-1, and IFN-γ signaling gene sets (Figure 2C). Foxp3 has been shown to directly interact with the CTLA-4 gene, which is significantly upregulated in activated human Tregs (38). We also found significant enrichment of the PI3Cl AKT pathway, mTOR pathway, and PI3K pathway (Supplemental Table 5), supporting our functional data that suggests hyperactive Akt/mTOR signaling in PD-1^hi Tregs. In contrast, PD-1^– Tregs only showed 4 enriched sets related to defensins and general transcription.

Interestingly, in circulating PD-1^hi Tregs, we found enrichment of several gene sets related to cell cycle and proliferation, such as G1 to S phase signaling, cell cycle checkpoints, and apoptosis (Supplemental Figure 8). Additionally, a small increase in annexin V staining in PD-1^hi Tregs was confirmed by flow cytometry (Supplemental Figure 7). A closer examination of the leading edge of these gene sets revealed several cyclin-dependent kinases and cyclin proteins, such as CDK1 and cyclin B1/A, as well as many proteasome subunits. Our transcriptional data, coupled with lowered proliferation, suggest that circulating PD-1^hi Tregs are functionally impaired, with limited proliferative capacity, but also have increased expression of genes related to proliferation, cell cycle, apoptosis, and inhibitory signaling pathways. This potential paradox led us to examine whether or not these cells exist in a state of replication-induced exhaustion. We performed GSEA with known exhaustion, Th-specific, memory, and activation gene signatures curated from the literature (Supplemental Table 6). We found that PD-1^hi Tregs were enriched in Th17, Th1, and Th2-like signatures and in several exhaustion gene signatures derived from both CD8 and CD4 exhausted T cells (Figure 3A). These include a coinhibitory signature derived from chronic LCMV infection in mice (24) and a signature derived from analysis of transcription factors in acute and chronically exhausted PD-1^hi CD8 T cells (39). Not surprisingly, PD-1^hi Tregs were also enriched in gene signatures specific to activated CD4⁺ Tregs. In contrast, PD-1^– Tregs were enriched in one gene set that was downregulated during chronic infection in mice. These data, coupled with the enrichment of cell cycle and proliferation gene sets, suggest that PD-1^hi Tregs are in the process of transitioning to a terminally differentiated and exhausted state.

Figure 3

Programmed cell death protein 1 (PD-1)^hi Tregs display shorter telomeres and may be clonally exhausted. (A) Gene set enrichment analysis of RNA-sequencing data with exhaustion signatures derived from the literature. The FDR for PD-1^hi and PD-1^– Treg populations (n = 3 each) is shown. (B) Relative telomere length (mean ± SD) assessed by flow cytometry ex vivo in sorted Tregs in 5 individuals (Student’s paired t test). (C) Percentage of telomere-specific demethylated region (TSDR) demethylation (mean ± SD) assessed ex vivo in sorted Tregs and T effector cells in 3 individuals (1-way ANOVA with Tukey’s multiple comparisons test). (D) FoxP3 expression intensity (mean ± SD) in sorted PD-1^hi and PD-1^– Tregs ex vivo, relative to the total Treg population from 8 individuals (Student’s t test).

We then sought to determine if PD-1^hi Tregs were indeed terminally differentiated. Telomeres protect the integrity of chromosome ends as the cell divides, and their length can indicate previous cell divisions. It has been reported that highly proliferative human memory Tregs have shortened telomeres and low telomerase activity (40). Moreover, this population displays close continuity with memory CD4 effector T cells, as established by analysis of T cell receptor repertoires, indicating a peripheral acquisition of regulatory properties (40). We measured the telomeres of sorted PD-1^hi, PD-1^–, and total Tregs and observed a significantly shorter telomere length in PD-1^hi Tregs, indicating replication-induced erosion (Figure 3B). In order to assess whether PD-1^hi Tregs showed hallmarks of a peripherally induced population, we assessed the degree of demethylation of the TSDR in the FOXP3 locus. Only thymic-induced Tregs present a fully demethylated TSDR, ensuring stable FoxP3 expression (41). We observed less TSDR demethylation in PD-1^hi Tregs compared with that in PD-1^– or the total Treg populations, although the locus was still highly demethylated compared with the level of effector CD4 T cells (Figure 3C). Additionally, we observe no decrease in FoxP3 expression in PD-1^hi Tregs (Figure 3D). These observations suggest that circulating PD-1^hi Tregs may represent a peripherally committed, replicatively exhausted population with lower TSDR demethylation.

Increased numbers of tumor-infiltrating Tregs in the brain tissue of patients with malignant gliomas is prognostic of poor survival, suggesting that there are increases in suppressor T cell function within the tumor environment (42, 43). Since tumor-infiltrating CD4 cells can express PD-1, we investigated the phenotype of PD-1^hi Tregs infiltrating gliomas as compared with that of those isolated from the blood to determine whether these Tregs were similarly dysfunctional as those characterized from healthy subjects. Tregs (CD4⁺CD25^hiCD127^lo) and T effector cells (CD4⁺CD25^loCD127⁺) were isolated from fresh tumor resections and matched peripheral blood of patients with gliomas of different WHO tumor grade (44) and were assessed for the presence of PD-1^hi Tregs (Supplemental Table 7). In agreement with previous studies (43), we observed an enrichment of Tregs within the CD4⁺ T cells infiltrating glioblastoma multiforme (GBM) tumors and not within gliomas of lower grade (Figure 4A). GBM-infiltrating Tregs expressed higher levels of CD25 and FoxP3, indicative of enhanced activation (Figure 4B). There was no difference in the frequency of PD-1^hi Tregs between the peripheral blood of healthy individuals and GBM patients, but there was a marked enrichment in PD-1^hi Tregs among tumor-infiltrating Tregs in malignant gliomas (Figure 4C). Interestingly, we also observed significantly increased numbers of TIM-3⁺ and PD-1^hiTIM-3⁺ Tregs (Figure 4, D and E, and Supplemental Figure 9) in tumors, which may suggest a role for TIM-3 in conjunction with PD-1 in decreased suppressive ability of tumor-derived Tregs (45). Since the circulating PD-1^hi Tregs in healthy individuals produced IFN-γ, we next examined whether malignant glioma tumor-infiltrating Tregs display an enhanced production of this proinflammatory cytokine. Indeed, after ex vivo restimulation with PMA and ionomycin, we observed a significant increase in IFN-γ, but not IL-10, secreted by tumor-infiltrating Tregs (Figure 4F).

Figure 4

Programmed cell death protein 1 (PD-1)^hi Tregs are enriched in tumors from glioblastoma multiforme patients, and they produce IFN-γ. (A) Frequency of Tregs among total CD4⁺ T cells from the blood (white) or tumor infiltrate (black) in 20 glioblastoma multiforme (GBM) patients (mean ± SD, 2-way ANOVA with Sidak’s multiple comparisons test; n = 20; G1–7, G10–18, G30–33). (B) Mean fluorescence intensity (MFI) of CD25 and Foxp3 within a gated infiltrating Treg population normalized to mean fluorescence intensity of patient peripheral blood–derived Tregs (mean ± SD, n = 16; G1–7, G10–12, G14–18, G30). Dotted line represents normalized patient peripheral blood (n = 7). (C) Frequency of PD-1^hi Tregs from the blood of 16 GBM patients and 16 healthy controls (mean ± SD, 2-way ANOVA with Sidak’s multiple comparisons test; n = 16; G1–9). (D) Frequency of Tim-3⁺ and (E) Tim-3⁺PD-1^hi Tregs in peripheral blood (white) and tumor populations of GBM patients (black) (mean ± SD, n = 18; G1–18). Dotted line represents average of healthy controls (n = 7; 2-way ANOVA with Sidak’s multiple comparisons test). (F) GBM patient Tregs were sorted and stimulated with PMA and ionomycin for 4 hours. IFN-γ and IL-10 release was measured by Luminex analysis (mean ± SD, Student’s t test; n = 6; G24–29).

We then assessed the functionality of Tregs isolated directly from malignant GBM compared with those from blood. Blood-derived Tregs displayed an intact suppressive activity that was directly proportional to the ratio between Tregs and T effector cell, but tumor-derived Tregs, which we have shown to be predominantly PD-1^hi, had a markedly reduced suppressive activity (Figure 5A; we note that the paucity of Tregs isolated from the brain limited analysis to a Treg/T effector cell ratio of 1:4 only). We also compared the functionality of Tregs from the blood of patients with GBM with respect to expression of PD-1 by sorting blood-derived PD-1^hi and PD-1^– Tregs and testing them in suppression assays at a Treg/T effector cell ratio of 1:2. Consistent with observations for circulating PD-1^hi Tregs from healthy individuals, circulating Tregs expressing PD-1 in GBM patients had reduced suppressive activity. Moreover, coexpression of PD-1 and Tim-3 on Tregs impaired their functionality even further (Supplemental Figure 10). Future studies should explore more deeply the functional state of Tregs expressing PD-1 alone or in conjunction with and TIM-3 in the circulation and the tumors of GBM patients.

Figure 5

Programmed cell death protein 1 (PD-1)^hi Tregs from glioblastoma multiforme patients display reduced suppression, enrichment for exhaustion signatures, and a network of coexpressed genes. (A) Peripheral blood or tumor-infiltrating Tregs from glioblastoma multiforme (GBM) patients cocultured with peripheral blood–derived CD4 effectors at a 1:4 ratio (mean ± SD, n = 8; G19–23, G34–36). Percent suppression calculated upon the proliferation of peripheral blood CD4 effectors without Tregs for each patient after stimulation with αCD3/αCD28/αCD2 for 5 days and ³H-thymidine addition over the last 18 hours (Student’s paired t test). (B) Gene set enrichment analysis (GSEA) of exhaustion signatures in blood and tumor-derived PD-1^hi versus PD-1^– Tregs from GBM patients (FDR is shown, n = 4 donors). (C) Network representation of correlation analysis (Pearson) identifying genes that strongly correlate with PD-1 in tumor-derived PD-1^hi Tregs (n = 4 donors). A full list is shown in Supplemental Table 10.

Next, we sought to determine if tumor-infiltrating PD-1^hi Tregs differed transcriptionally from their counterparts in the blood of patients. We isolated PD-1^hi and PD-1^– Tregs from the blood and tumors of GBM patients and performed RNA sequencing. High-quality biological replicates were included for subsequent analysis on the basis of sequencing metrics and expression of housekeeping genes, and overall these showed high correlation across biological replicates (Supplemental Figures 5 and 11, and Supplemental Table 1). Overall, when comparing PD-1^hi or PD-1^– Tregs in tumors versus those in blood, we found 412 and 732 differentially expressed genes, respectively. Filtering for the top differentially expressed genes (top 150) and comparing PD-1^hi to PD-1^– Tregs, we found 41 genes that were commonly up in tumor-derived Tregs, including IFN-γ (Supplemental Figure 12), as well as gene sets unique to PD-1^hi and PD-1^– Tregs (Supplemental Figure 13). Comparing PD-1^hi or PD-1^– Tregs from the blood of healthy control and GBM patients, we found only a handful of differentially expressed genes (Supplemental Table 8). Additionally, the PD-1^hi Tregs isolated from tumors were enriched in exhaustion signatures compared with their counterparts in the blood (Figure 5B and Supplemental Table 9). We performed co-correlation analysis to determine which genes correlated with PD-1 expression (Supplemental Table 10). Several exhaustion markers (PD-L1, LAG3, and ICOS), transcription factors (STAT4, SUMO3), and signaling molecules (PYCARD, TNFSF4) were significantly correlated with PD-1^hi expression in tumors but not in blood. Top-scoring molecules correlated with PD-1^hi expression in tumor-derived Tregs (P > 0.9) were then analyzed for functional connectivity to determine if these molecules are known to interact (Figure 5C). We found that these molecules are highly coexpressed, with several known to physically interact. Overlaying gene ontology functions onto this network revealed that this PD-1^hi Treg gene set includes regulation of proliferation and T cell costimulation. Taken together with the functional suppression data, transcriptional profiling reveals that tumor-derived PD-1^hi Tregs, like those derived from healthy donors, are exhausted. However, PD-1^hi Tregs derived from tumors have a unique signature of genes co-correlated with PD-1^hi expression compared with those in circulation. These results suggest that the majority of Tregs in circulation may not fully reflect the extent of dysfunction of the cells residing at the site of disease in GBM.

PD-1 signaling ultimately represses the transcription of IFN-γ by inducing nuclear translocation of FoxO-1; however, we found that PD-1^hi Tregs display higher levels of FoxO-1 phosphorylation. Thus, we sought to determine whether blocking PD-1 in vivo affected the production of IFN-γ by Tregs. We isolated Tregs from the circulation of 8 subjects with recurrent GBMs in a phase I trial of an anti-PD-1 monoclonal antibody and assessed for IFN-γ production after 3 infusions separated by 14 days of treatment. Following in vivo treatment, there were significant increases in IFN-γ production by circulating Tregs (Figure 6A). These data suggest that anti-PD-1 treatments in vivo may drive less-differentiated PD-1^intermediate Tregs, in the middle of the continuum to exhaustion, toward a dysfunctional state (Figure 6B), and ultimately loosen tight regulation on CD8 populations. Thus, PD-1^hi Tregs both in the circulation of healthy subjects and in tumors are dysfunctional and exhibit enhanced secretion of IFN-γ.

Figure 6

Anti–programmed cell death protein 1 treatment shifts the profile of Tregs in vivo to exhaustion and IFN-γ production. (A) Production of IFN-γ relative to baseline (mean ± SEM) by sorted circulating Tregs from 8 glioblastoma multiforme (GBM) patients undergoing nivolumab treatment (3 mg/kg every 2 weeks) over the first 5 weeks. Tregs were restimulated for 6 hours with PMA and ionomycin in the presence of brefeldin A. Significance (1-way ANOVA) is P < 0.05 for week 2 and P < 0.001 for week 5. (B) Summary of the potential effect of anti–programmed cell death protein 1 (anti-PD-1). Tregs were gated as CD127^lo and Foxp3^hi for A and B. Gating strategy is shown in Supplemental Figure 12. (C) Conditional density resampled estimate of mutual information (DREMI) scores computed with a noise threshold of 0.85 and an arch-transformed range at 0.3 to maximize protein abundance. Heatmap of DREMI scores of the strength of the statistical dependency between each indicated protein and PD-1 for blood week 1, blood week 7, and tumors. (D) Conditional density rescaled visualization (DREVI) plots of each indicated protein against PD-1 for Tregs from blood week 1, blood week 7, and tumors. Red = high dependency, blue = low dependency.

Finally, we sought to determine how protein expression profiles of single Tregs changed over the course of treatment with α-PD-1-blocking antibody (nivolumab) and how the tumor and blood compartments in patients differed. We used CyTOF (46) to profile blood obtained at week 1 and 7 after treatment as well as matched tumors from 4 patients (Supplemental Tables 7 and 11). CyTOF data were gated in silico for Tregs as CD4⁺CD8^–CD45⁺CD3⁺FoxP3⁺CD127^– (Supplemental Figure 14). As with our transcriptional data, we found that tumor and blood profiles were strikingly different and displayed marked heterogeneity (Supplemental Figure 15). Visualizing the data using t-distributed stochastic neighbor embedding, we found distinct clusters of cells, such as a strong CTLA-4 cluster in week 7 blood and IFN-γ in blood and tumors (Supplemental Figure 16). However, these analyses make it challenging to identify relationships between each molecule measured. To better quantify the strength of the relationship between PD-1 and other proteins measured, we computed density resampled estimate of mutual information (DREMI) scores (47) and displayed these as a heatmap ordered by score within the tumor (Figure 6C). Calculated DREMI scores are robust and have an average standard deviation of 0.02 in week 1 blood, 0.03 in week 7 blood, and 0.03 in brain (Supplemental Figure 17). We found that CXCR3, granzyme B, CCR6, ICOS, and LAG3 had high scores within the tumor, indicating that these molecules have a strong relationship with PD-1 within tumor tissue. We also found several molecules, such as CD27, LAG3, and IL-10, whose strength of dependency on PD-1 increased over the course of treatment from week 1 to week 7, suggesting that α-PD-1 treatment shifted the phenotype of Tregs in circulation. This, along with our RNA-sequencing data, suggests distinct differences with respect to correlations with PD-1 expression between blood- and tumor-derived Tregs, potentially driven by treatment. We next sought to determine the nature of the relationship between top-ranking markers identified by the DREMI scores and PD-1 in tumor-derived Tregs. We then used conditional density rescaled visualization (DREVI) plots (47) to visualize the relationship between PD-1 and the top 4 proteins identified by their DREMI score in order to determine the nature of the relationship between these molecules and PD-1. DREVI plots are shown in Figure 6D for CXCR3, granzyme B, IL-10, CCR6, and LAG3 in week 1 blood, week 7 blood, and tumors. Overall, through the course of treatment, profiles were shifted in blood to reflect those seen in tumors. This is especially true for CXCR3, which shifted from a discontinuous dependency on PD-1, indicative of a bistable state of two distinct populations, to a monotonic linear response with PD-1, suggesting that this protein may be continuously expressed with PD-1 in a graded fashion that is driven by treatment. Taken together, these data indicate that there are distinct differences between blood- and tumor-derived Tregs, where coexpression of LAG3 and CXCR3 may be predictive of highly dysfunctional Tregs in the context of GBM, and treatment with α-PD-1 may drive Tregs toward this dysfunctional state.

We describe a decrease in the suppressive activity of tumor-infiltrating Tregs and circulating Tregs from GBM patients associated with PD-1^hi expression. We found that PD-1^hi Tregs in both the circulation of healthy individuals and that of patients with malignant gliomas share impaired suppressive activity and similar molecular signatures of T cell exhaustion. Circulating CD4⁺ PD-1^hi Tregs lose their proliferative potential but maintain the capacity to produce cytokines with partial demethylation at the TSDR locus. While it should be noted that fate-mapping studies cannot be performed in humans and these PD-1^hi Tregs may represent induced or highly activated FoxP3 T effector cells as opposed to thymically derived Tregs, the PD-1^hi Tregs have the same or increased levels of FoxP3 expression as compared with PD-1^– Tregs. These findings are consistent with those of recent studies demonstrating recruitment of methyltransferases at the TSDR under proinflammatory conditions (48) and loss of TSDR demethylation and eventual FoxP3 expression in Tregs where the PI3K pathway was hyperactive (49), leading us to speculate that in PD-1^hi Tregs loss of demethylation may be related to chronic stimulation and loss of control of a well-described Treg-stabilizing signaling pathway. Given that PD-1 is a repressor of the PI3K/Akt axis, the finding that PD-1^hi Tregs produce IFN-γ and display enhanced phosphorylation of FoxO-1 is surprising. Nevertheless, PD-1^hi Tregs from healthy donors and GBM patients also overexpressed the ICOS receptor, which has been linked both to phosphorylation of FoxO-1 (50) and to IFN-γ production by Tregs (51). It is thus of interest that ICOS was among the top-ranking genes coexpressed with PD-1 in PD-1^hi Tregs from tumors in GBM patients. These data could suggest that expression of PD-1 in human Tregs may be induced as a compensatory mechanism to stabilize the PI3K/Akt pathway and repress IFN-γ. Alternatively, a recent assessment of the level of activity of FoxO-1 in resting versus activated Tregs in mice reports that inactivation of FoxO-1 is a prerequisite for migration of Tregs to peripheral sites. As such, enhanced phosphorylation of FoxO-1 in PD-1^hi Tregs might indicate that they are indeed recirculating from the tissue (52).

We identified a molecular state of Treg exhaustion found in infiltrating malignant GBM characterized by the high expression of PD-1 and IFN-γ secretion. A similarly exhausted RNA signature with IFN-γ secretion was observed in PD-1^hi Tregs in the circulation of healthy subjects. The observation that Treg function can be restored by blocking IFN-γ in Tregs from patients with autoimmune disease but not in PD-1^hi Tregs suggests that Tregs in autoimmunity are in a “pre-exhausted” transition state. In contrast, we speculate that there are multiple functional populations of tumor-infiltrating Tregs that can differentiate to an exhausted state identified by PD-1^hi, LAG-3, and CXCR3 expression as well as IFN-γ secretion that can be further disabled by anti-PD-1 monoclonal antibody treatment. This finding is consistent with a previous report on the role of the PD-1/PD-L1 axis in promoting Treg induction through inhibition of the Akt/mTOR pathway (27) as well as with the observation that PD-1^hi Tregs may represent a peripherally induced population as related to their reduced TSDR demethylation. As such, anti-PD-1 therapy may work by destabilizing a crucial pathway for the function and de novo induction of Tregs. The observation of dysfunctional Tregs in cancer is in contrast with the current view that Tregs are hyperfunctional and suppress antitumor immune responses. Nevertheless, a role for Tregs in promoting tissue repair independent from T effector suppressor functions through secretion of growth factors has been recently reported (53). It can be speculated that Tregs could directly promote tumor proliferation through production of oncogenic growth factors rather than by suppression of the antitumoral response. These considerations, coupled with our CyTOF analysis of tumor-infiltrating T cells in patients treated with PD-1-blocking antibody, suggest further single-cell molecular characterization of tumor-infiltrating Tregs that may allow more targeted immunomodulatory treatments in cancer. Finally, more broad comparisons between Tregs in autoimmunity and cancer should identify key biomarkers of dysfunction in these contrasting immune states.

Human patients

All human tissue and peripheral blood were collected from patients with brain tumors, undergoing surgical resection, or healthy individuals (Supplemental Table 7). Tumors were histologically characterized by the Pathology Department at the Yale–New Haven Hospital according to the 2007 WHO classification 59, and immunohistochemical Ki-67 staining was performed to measure the proliferative index of each tumor. Immune monitoring samples were taken from patients as part of the CheckMate 143 clinical trial (ClinicalTrials.gov NCT02017717) of nivolumab for recurrent GBM. In this trial, nivolumab is given intravenously at a dose of 3 mg/kg every 2 weeks, and blood was drawn as indicated after treatment.

Flow cytometry

For immunophenotyping analyses, flow cytometry was performed with the following anti-human surface and intracellular antibodies and dyes: anti-CD4 (RPA-T4), anti-CD8 (RPA-T8), anti-CD25 (M-A251), anti-CD127 (HIL-7R-M21), anti-PD-1 (EH12.1), anti-PDL1 (MIH1) anti-IFN-γ (B27), anti-KI67 (B56), anti-CCR6 (11A9), anti-CXCR3 (IC6), anti-CCR7 (3D12), anti-CXCR5 (RF8B2), and V500 Annexin-V stain (all from BD Biosciences); anti-TIM-3 (F38-2E2), anti-CTLA-4 (L3D10), anti-CD226 (TX25), and anti-IL-10 (JES3-9D7) (all from Biolegend); anti-Foxp3 (PCH101), anti-IL-2 (MQ1-17H12), and anti-ICOS (ISA3) (all from eBioscience); anti-LAG-3 (polyclonal goat IgG, R&D Systems); and anti-phosphoSer319 FoxO-1 (polyclonal rabbit IgG, Bioss). Cell viability was assessed using Live/Dead Cell Viability Assays (Life technologies). Samples were run using a BD LSRFortessa for all analysis, and FACS sorting was carried out on a BD FACS Aria II. FlowJo software (TreeStar Inc.) was used for data collection and analysis.

Tumor-infiltrating lymphocyte isolation

For extraction of infiltrating lymphocyte populations, freshly resected tumor specimens were manually disrupted, followed by digestion with collagenase IV (2.5 mg/ml) and DNase I (0.2 mg/ml) (Worthington Biochemical Corporation) for 1 hour. Tumor homogenates were passed through a 70-μm cell strainer prior to separation on discontinuous 70-30% Percoll (Sigma-Aldrich) gradients. Matched fresh PBMCs were isolated by Ficoll-Hypaque gradient centrifugation.

Luminex

Multiplex assays for cytokine secretion were performed on cell culture supernatants after stimulation with PMA (50 ng/ml) and ionomycin (250 ng/ml) using the Milliplex Human Cytokine/Chemokine Magnetic Bead Panel (Millipore).

Suppression assays (CFSE and thymidine)

Tregs (CD4⁺CD25⁺CD127^–) were sorted alongside CD4 T effector cells (CD4⁺CD25^–CD127⁺) from either tumor-infiltrating lymphocytes or PBMCs. To assess functional potential of tumor-infiltrating lymphocytes, both infiltrating and peripheral blood Tregs were cultured alongside peripheral blood CD4 effector T cells at ratios ranging from 1:1 to 1:8 with no fewer than 2,000 effector T cells in the presence of 2 Treg Inspector Beads (anti-CD2/anti-CD3/anti-CD28) (Miltenyi) per cell. After 72 hours, cells were pulsed with 0.5 μCi 3H-thymidine and harvested 18 hours later. 3H-thymidine incorporation was measured using a microbeta counter. To assess suppressive potential of PD-1^hi or PD-1^– Tregs, both CD4 effector T cells and Tregs from peripheral blood were further sorted based upon PD-1 surface expression into PD-1^hi and PD-1^– populations, alongside a total population that was not sorted by PD-1 expression. CD4 effector T cells were stained with CellTrace CFSE (1 μM) and Tregs were stained with CellTrace Violet (1 μM) (both Life technologies) to assess proliferation through dilution of the dye. 3 × 10³ Tregs were cultured with 6 × 10³ CD4 effector T cells for 4 days before viability staining, and samples were run on a BD LSRFortessa (BD Biosciences) and analyzed using FlowJo software (TreeStar Inc.).

Blocking antibodies

Human anti-PD-1 (M3, mouse IgG1) and anti-PD-L1 (5H1, mouse IgG1) were gifts from Leiping Chen (Yale School of Medicine). Human anti-IFN-γ (AF-285-NA, R&D Systems) were used at a concentration of 5 μg/ml, with appropriate isotype control antibodies whenever used.

Telomere length

The telomere length of sorted lymphocyte populations was assessed using the TelomerePeptide Nuclear Antigen Kit for flow cytometry (Dako). Relative telomere length was calculated as the ratio between the telomere signal from each sorted lymphocyte population and the 1301 control cell line (Sigma-Aldrich).

FOXP3 methylation analysis

Genomic DNA was isolated from sorted Treg (CD4⁺CD25⁺CD127^–) and CD4 effector (CD4⁺CD25^–CD127⁺) populations, and DNA methylation analysis of the FOXP3 locus was carried out by Epiontis.

RNA sequencing

T cells from healthy controls and GBM patients were sorted into RNA lysis buffer for RNA sequencing as described previously (54). Briefly, RNA was extracted using the NucleoSpin RNA XS Kit (Macherey-Nagel) according to the manufacturer’s instructions. cDNA synthesis and amplification were performed using SMARTer Ultra Low Input RNA for the Illumina Sequencing High Volume Kit (Clontech) according to the manufacturer’s instructions, with the following for LD amplification: 98°C for 10 seconds, 65°C for 30 seconds, and 68°C for 3 minutes 18 times and 72°C for 10 minutes. cDNA was normalized to an input of 0.5 ng total after quantification using Picogreen (Invitrogen) for preparing sequencing libraries. Paired-end sequencing libraries were prepared using the Nextera XT DNA sample Prep Kit (Illumina) according to the manufacturer’s instructions. Libraries were pooled in an equimolar ratio and sequenced on a HiSeq 2500 sequencer with 200 cycles per lane (Illumina).

Reads were mapped against UCSC’s known gene annotation of the hg19 human genome assembly using RSEM (55) v 1.2.15 and bowtie (56) 1.0.1 (parameters -p 4 --output-genome-bam --calc-ci --paired-end --bowtie-chunkmbs 1024). Posterior estimates of counts per genes were retrieved for processing in the R statistical environment v. 3.1.1 using the DESeq package (57). After normalization (estimateSizeFactors) and variance estimates across all samples using the estimateDispersions function in a method similar to that of this program’s original implementation (sample-blind estimate and fit-only sharing mode), differential expression was performed either pairwise between every combinations of sample labels (nbinomTest) or by building a GLM model to test every condition independently (fitNbinomGLMs). P values were adjusted using the Benjamin-Hochberg procedure (p.adjust).

Quality control was performed on forward and reverse fastq files from each sample using custom python scripts. Quality control metrics used were number and percentage of mapped reads and reads mapping to exonic, intronic, CDS, 5′ and 3′ untranslated regions, intergenic regions, and different classes of noncoding RNAs. Additional metrics reflecting mRNA enrichment were derived from exonic/intronic read density ratios, and biases in RNA integrity were detected using the 5′/3′ ratios. Samples kept for downstream analysis are shown in Supplemental Table 1. The log2FPKM values for mapped transcripts in each sample were used to create the input.gct file for GSEA (v2.0) as previously described (58). GSEA (classic scoring scheme, 15 set minimum, 1,000 permutations) was performed using C2CP gene sets and custom gene sets curated from the literature. Data were further analyzed using DAVID for GO enrichment and the GeneE data visualization tool (http://www.broadinstitute.org/cancer/software/GENE-E/). Correlation analysis with PD-1 expression was performed by computing a vector of log-transformed PD-1 expression levels across tumor-infiltrating and blood-derived Tregs in patients. Pearson’s correlation coefficients were computed between these PD-1 expression vectors and similar vectors for all the other genes under consideration. P values and correlation coefficients were calculated using cor.test in R. Genes with correlation values greater than 0.9 were visualized using GeneMania (59). All RNA sequencing data are available in the database of genotypes and phenotypes (accession number phs001079.v1.p1).

Mass cytometry analysis

PBMCs and tumor-infiltrating lymphocytes were isolated as described above. The cells were counted, adjusted to 10⁶/ml, and incubated at 37°C in complete RPMI medium with or without PMA (250 ng/ml) and ionomycin (50 ng/ml). After 4 hours, the cells were fixed with 2% PFA for 10 minutes at room temperature and stored at –80°C until use. Samples were thawed in batches of 20 at 37°C, barcoded using the Cell-ID 20-Plex Barcoding Kit (Fluidigm). The 20-barcoded samples were washed 3 times with staining buffer (Fluidigm) and then combined into a single sample, incubated with Fc block (BD Biosciences), followed by 30 minutes of staining at room temperature with a cocktail of metal-conjugated antibodies (Supplemental Table 10) in staining buffer. The combined sample was washed once with 2% PFA and incubated for 12 hours with an iridium-containing DNA intercalator (Fluidigm). The cell sample was washed twice with PBS and once with water. The cell concentration adjusted to 10⁶/ml in water supplemented with 10% EQ 4-element calibration beads (Fluidigm), filtered, and injected successively to a CyTOF II mass cytometer in a series of 500 μl aliquots. At the end on the acquisition, the files resulting from the different injections were concatenated into one fcs file, and the calibration beads were removed using the CyTOF II built-in software. Detector sensitivity was performed at the beginning and at the end of each use of CyTOF II using polystyrene normalization beads containing lanthanum-139, praseodymium-141, terbium-159, thulium-169, and lutetium-175.

CyTOF analysis

DREVI for 2-dimensional relationship visualization. In order to visualize 2-dimensional relationships in single-cell data, we use a recently developed method known as DREVI (47). DREVI is designed to highlight the shape and spread of 2-variable relationships in single-cell data. It depicts the Y-variable as a function of the X-variable despite noise in measurements and uneven coverage along the dynamic range of the X-variable. DREVI computes a 2-dimensional kernel density estimate p(X,Y) on a fine-grained grid spanning the dynamic range of the 2 variables using a heat-diffusion method (60). Then, the density estimate is column-normalized to obtain the conditional-density estimate p(Y|X). The result is visualized as a rescaled heatmap where the densest regions along the y axis are shown in red for each X-slice, and the sparser regions are shown in blue. This way, the X-Y relationship is revealed even in regions of sparse data, because each x-slice is renormalized individually. As compared with an ordinary density estimate, the visualization shows how the y axis changes with the x axis rather than simply highlighting where the majority of the data is concentrated. CyTOF data were further visualized using t-distributed stochastic neighbor embedding (61, 62).

DREMI metric for relationship strength. DREMI is a measure of functional dependency between two variables and was developed coordinately with DREVI (47). The DREMI metric is computed by starting with the conditional-density estimate for DREMI (as described in the previous section). The conditional-density estimate is then denoised using a threshold that eliminates noisy regions of low density in each X-slice. Then, data are resampled from the denoised conditional-density estimate, generating an even number of samples for each slice of the x axis. Finally, the x and y axes are repartitioned into large partitions (using 8 partitions), and mutual information is computed on the resampled points. This method of computing relationship strength has the advantage that it automatically filters out sparse, unreplicated regions of the 2-dimensional plane, while providing a score on the functional X-Y relationship through the full dynamic range of the x axis, regardless of where the majority of cells are concentrated.

Louvain method for unsupervised clustering

The Louvain method (63) is a method that finds “communities” or unsupervised partitions of a graph. This method was also recently used to cluster AML populations (64). For single-cell analysis, a cell-cell similarity graph is drawn by including a vertex in the graph for each cell and connecting 2 vertices x and y when y is within the 30 nearest neighbors of x. The next step in the Louvain method is agglomerating cells into a cluster as it optimizes a cluster quality known as a modularity. Modularity represents the fraction of the edges that fall within a given group minus the expected fraction if edges were distributed at random. The advantage of the Louvain method is that it scales to high dimensions and large numbers of cells and does not require the prior knowledge of the number of clusters or groups within a population.

Statistics

Statistical analysis was carried out using GraphPad Prism (GraphPad Inc.). Paired and unpaired 2-tailed Student’s t tests were used when appropriate. P < 0.05 was considered statistically significant.

Study approval

All patients with brain tumors and healthy donors provided samples after providing informed consent and approval by the Institutional Review Board at the Yale School of Medicine.