Induction and transcriptional regulation of the co-inhibitory gene module in T cells

Abstract

The expression of co-inhibitory receptors, such as CTLA-4 and PD-1, on effector T cells is a key mechanism for ensuring immune homeostasis. Dysregulated expression of co-inhibitory receptors on CD4⁺ T cells promotes autoimmunity, whereas sustained overexpression on CD8⁺ T cells promotes T cell dysfunction or exhaustion, leading to impaired ability to clear chronic viral infections and diseases such as cancer1,2. Here, using RNA and protein expression profiling at single-cell resolution in mouse cells, we identify a module of co-inhibitory receptors that includes not only several known co-inhibitory receptors (PD-1, TIM-3, LAG-3 and TIGIT) but also many new surface receptors. We functionally validated two new co-inhibitory receptors, activated protein C receptor (PROCR) and podoplanin (PDPN). The module of coinhibitory receptors is co-expressed in both CD4⁺ and CD8⁺ T cells and is part of a larger co-inhibitory gene program that is shared by non-responsive T cells in several physiological contexts and is driven by the immunoregulatory cytokine IL-27. Computational analysis identified the transcription factors PRDM1 and c-MAF as cooperative regulators of the co-inhibitory module, and this was validated experimentally. This molecular circuit underlies the co-expression of co-inhibitory receptors in T cells and identifies regulators of T cell function with the potential to control autoimmunity and tumour immunity.

Headnote

The expression of co-inhibitory receptors, such as CTLA-4 and PD-1, on effector T cells is a key mechanism for ensuring immune homeostasis. Dysregulated expression of co-inhibitory receptors on CD4⁺ T cells promotes autoimmunity, whereas sustained overexpression on CD8⁺ T cells promotes T cell dysfunction or exhaustion, leading to impaired ability to clear chronic viral infections and diseases such as cancer1,2. Here, using RNA and protein expression profiling at single-cell resolution in mouse cells, we identify a module of co-inhibitory receptors that includes not only several known co-inhibitory receptors (PD-1, TIM-3, LAG-3 and TIGIT) but also many new surface receptors. We functionally validated two new co-inhibitory receptors, activated protein C receptor (PROCR) and podoplanin (PDPN). The module of coinhibitory receptors is co-expressed in both CD4⁺ and CD8⁺ T cells and is part of a larger co-inhibitory gene program that is shared by non-responsive T cells in several physiological contexts and is driven by the immunoregulatory cytokine IL-27. Computational analysis identified the transcription factors PRDM1 and c-MAF as cooperative regulators of the co-inhibitory module, and this was validated experimentally. This molecular circuit underlies the co-expression of co-inhibitory receptors in T cells and identifies regulators of T cell function with the potential to control autoimmunity and tumour immunity.

We used single-cell RNA sequencing (scRNA-seq) to analyse coinhibitory and co-stimulatory receptor expression in 588 CD8+ and 316 CD4+ tumour-infiltrating lymphocytes (TILs) from B16F10 mouse melanoma3. We found that the expression of Pdcd1 (also known as PD-1), Tim3 (Havcr2), Lag3, Ctla4, 4-1BB (Tnfrsf9) and Tigit strongly co-vary in CD8+ TILs. CD4+ TILs showed a similar pattern with the additional co-expression of Icos, Gitr (also known as Tnfrsf18) and Ox40 (Tnfrsf4) (Fig. 1a, top). Single-cell mass cytometry (cytometry by time of flight, CyTOF) confirmed the surface co-expression of these receptors (Fig. 1a, bottom, Supplementary Information 1). The expression of PD-1, LAG-3, TIM-3 and TIGIT was tightly correlated on both CD8+ and CD4+ TILs (Fig. 1a, bottom). Clustering analysis (t-stochastic neighbourhood embedding (t-SNE)4, Methods) showed two groups of CD8+ TILs (clusters 1 and 2) (Fig. 1b, Extended Data Fig. 1a, c), with PD-1, LAG-3, TIM-3 and TIGIT mainly expressed in cluster 1 cells (Fig. 1b, Extended Data Fig. 1c), in addition to LILRB4 (Extended Data Fig. 1a) and co-stimulatory receptors of the TNF receptor family, 4-1BB, OX40 and GITR. By contrast, ICOS and CD226 were less restricted to cluster 1 (Extended Data Fig. 1a). We further observed two discrete clusters of CD4+ TILs (clusters 3 and 4), with co-expression of PD-1, TIM-3, LAG-3 and TIGIT restricted to cluster 3 (Fig. 1b, Extended Data Fig. 1c).

The co-expression of co-inhibitory receptors on CD8+ and CD4+ T cells suggests a common trigger. One candidate is IL-27, a heterodimeric member of the IL-12 cytokine family that suppresses autoimmunity5, induces IL-10-secreting type 1 regulatory T (Treg) cells6,7 and induces expression of TIM-3 and PD-L1 on CD4+ and CD8+ T cells8,9. Activation of CD4+ and CD8+ T cells in the presence of IL-27 induced the expression of TIM-3, LAG-3 and TIGIT at both the mRNA (Fig. 1c) and protein levels (Extended Data Fig. 2a). mRNA expression of Tim3 (Havcr2), Lag3 and Tigit was reduced in IL-27RAdeficient T cells, whereas Pdcd1 expression was unaffected by IL-27 in vitro (Fig. 1c, Extended Data Fig. 2a).

CyTOF analysis showed that the loss of IL-27RA resulted in the loss of cells in cluster 1 of CD8+ TILs and cluster 3 of CD4+ TILs (Fig. 1d, P = 5 × 10-23 and 6.8 × 10-7 for CD8+ and CD4+, respectively, hypergeometric test; Extended Data Fig. 1b-d), indicating a key role for IL-27 in driving co-inhibitory receptor co-expression in both CD4+ and CD8+ T cells in vivo. Although PD-1 expression was not dependent on IL-27 in vitro, it was dependent on IL-27RA signalling in vivo. Consistent with the induction of IL-10 by IL-275-7, we observed reduced IL-10 in IL-27RA-knockout CD8+ TILs (Extended Data Fig. 2b).

scRNA-seq of CD8+ and CD4+ TILs from wild-type and IL-27RAknockout mice (Fig. 1e, Extended Data Fig. 3a, b, Methods) revealed distinct clusters of CD8+ (cluster 5) and CD4+ (cluster 4) TILs that highly expressed the co-inhibitory receptors Pdcd1, Tim3, Lag3 and Tigit. The expression of these genes was decreased in CD8+ TILs from IL-27RA-knockout mice, whereas the expression of only Tim3 and Lag3 was decreased in CD4+ TILs from IL-27RA-knockout mice (Fig. 1e). Thus, IL-27 drives a module of co-inhibitory receptors that are strongly co-expressed in vivo together with IL-10.

The co-inhibitory receptor module could be part of a larger IL-27- driven inhibitory gene program. We analysed the mRNA profiles of CD4+ and CD8+ T cells stimulated in the presence or absence of IL-27. IL-27 induced similar expression programs in CD4+ and CD8+ T cells (Extended Data Fig. 4a, b). We identified 1,201 genes with IL-27- dependent expression (Methods). We compared the IL-27-driven gene program to the gene signatures for four different states of T cell nonresponsiveness: CD8+ T cell exhaustion in both cancer3 and chronic viral infection10, and antigen-specific11 and non-specific (anti-CD3 antibody12) CD4+ T cell tolerance. We found a significant overlap with all of these signatures (Methods, Extended Data Fig. 4c-f).

Projection of the IL-27 and CD8+ cancer T cell exhaustion overlap signature onto the single-cell profiles of CD8+ TILs marked a distinct subset of cells (Fig. 2a, panel I). This subset scored highly for the overlap signatures between the IL-27-driven gene program and each of the other three states of T cell non-responsiveness (Fig. 2a, panels II-IV). The transcriptional program induced in IL-27RA-knockout TILs was active in a complimentary subset of TILs (Fig. 2a, panel V, Methods). The control signature from cells stimulated with IL-27 in vitro showed bimodal distribution and by itself did not detect the same population of cells (Fig. 2a, panel VI). From these analyses, we identified a co-inhibitory gene module (272 genes) that is shared across several states of T cell non-responsiveness (Supplementary Table 2). Within this module, we identified a set of 57 genes that encode cell-surface receptors and cytokines, including TIM-3, LAG-3, TIGIT and IL-10 (Fig. 2b), which we further stratified by their expression in cancer and chronic viral infections (Fig. 2c). Two surface molecules, PROCR and PDPN, were highly expressed in the setting of cancer (Fig. 2c). Activation of naive CD4+ and CD8+ T cells in vitro in the presence of IL-27 induced the expression of PROCR and PDPN (Extended Data Fig. 5a). In vivo, PROCR and PDPN exhibited IL-27-dependent coexpression with PD-1 and TIM-3 on CD8+ TILs (Extended Data Fig. 5b).

PROCR+CD8+ TILs exhibited an exhausted phenotype, producing less TNF and IL-2 and more IL-10 than PROCR-CD8+ TILs (Extended Data Fig. 5c). The growth of B16F10 melanoma was inhibited in PROCR hypomorph (Procrdelta/delta, hereafter Procrd/d)13 mice (Fig. 2d), and Procrd/d CD8+ TILs mice exhibited enhanced production of TNF, but no difference in the production of IL-2, IFN-γ or IL-10 (Fig. 2e). Procrd/d TILs exhibited a decreased frequency of TIM-3high and PD-1high CD8+ T cells, suggesting that PROCR signalling promotes a severely exhausted phenotype in CD8+ T cells14 (Fig. 2f). Adoptive transfer of CD8+ T cells that lack PROCR revealed a T cell-specific role for PROCR in constraining tumour growth (Extended Data Fig. 5d).

Although PDPN can limit CD4+ T cell survival in inflamed tissues15, its role in T cell exhaustion is unknown. We observed a significant delay in B16F10 tumour growth in mice with PDPN deficiency in T cells (PDPN conditional knockout (cKO)) (Fig. 2g). PDPN-deficient CD8+ TILs exhibited enhanced TNF production but no significant difference in IL-2, IFN-γ or IL-10 (Fig. 2h). The frequency of TIM-3high and PD-1high CD8+ TILs was decreased, indicating a reduced accumulation of T cells with a severely exhausted phenotype in PDPN cKO mice14 (Fig. 2i). Consistent with previous data15, PDPN-deficient PD-1+TIM-3+ CD8+ TILs had higher expression of IL-7RA, indicating that PDPN may limit the survival of CD8+ TILs in the tumour microenvironment (Extended Data Fig. 5e, f).

We identified the transcription factor PRDM1 as a candidate regulator of the co-inhibitory module. PRDM1 is induced in vitro by IL-27 in CD4+ and CD8+ T cells (Extended Data Fig. 6a), is enriched in TILs with high expression of the IL-27 co-inhibitory module (Extended Data Figs. 3c-f, 6b, c, Methods), and is overexpressed in exhausted CD8+ TILs (P = 0.0004, t-test, Extended Data Fig. 6d). Network analysis based on profiling of naive CD8+ T cells from mice with a T cell-specific deletion of PRDM1 (PRDM1 cKO) stimulated with IL-27, showed that PRDM1 regulates several genes in the IL-27 co-inhibitory module (Extended Data Fig. 6e, P = 2.32 × 10-12; hypergeometric test; Methods). This was further supported by PRDM1 chromatin immunoprecipitation followed by sequencing (ChIP-seq) data16 (P = 2.9 × 10-8, hypergeometric test; Extended Data Fig. 6e, Methods).

CD8+ TILs from B16F10 tumour-bearing PRDM1 cKO mice expressed lower levels of TIM-3, PD-1 and PROCR (Fig. 3a); however, there was no difference in tumour growth compared to wildtype controls (Fig. 3b), indicating that the reduction of co-inhibitory receptor expression in PRDM1 cKO mice was insufficient to promote effective anti-tumour immunity. We therefore examined whether other transcription factors may regulate the co-inhibitory module and compensate for the absence of PRDM1. We analysed CD8+ TILs from PRDM1 cKO mice for the expression of genes from the IL-27-driven gene signature and the signature for exhausted CD8+ TILs (Methods, Supplementary Table 3). We found that only a few genes were upregulated in PRDM1 cKO CD8+ T cells, including one transcription factor, c-MAF (P < 0.05; Fig. 3c). Indeed, c-MAF is induced by IL-27, is coexpressed with PRDM1 in T cells after IL-27 stimulation (Extended Data Fig. 6a), and can regulate IL-10 expression17 and T cell exhaustion18. In addition, many genes (226 genes, P = 5.34 × 10-5, hypergeometric test) in the co-inhibitory gene module have a binding motif and a reported binding event for c-MAF within their promoter regions19.

CD8+ TILs from c-MAF cKO mice exhibited decreased expression of several co-inhibitory receptors (Fig. 3d). PRDM1 and c-MAF each affected co-inhibitory receptor expression only partially (Fig. 3e). As in the PRDM1 cKO mice, c-MAF cKO mice did not show any differences in tumour growth relative to controls (Fig. 3f). Notably, PRDM1 expression in c-MAF cKO TILs was similar to that in wild-type TILs, indicating that PRDM1 might drive the expression of the co-inhibitory gene module in the absence of c-MAF.

We addressed whether PRDM1 and c-MAF could act cooperatively to regulate co-inhibitory receptor expression. We found no evidence for a physical interaction between PRDM1 and c-MAF (data not shown); we therefore examined whether they shared targets. We combined the network analysis for PRDM1 (Extended Data Fig. 6e) with c-MAF ChIP-seq data19 and c-MAF targets (Methods). We observed 121 genes in the co-inhibitory module that are affected (RNA-seq) or have a direct binding event (ChIP-seq) for both PRDM1 and c-MAF (Fig. 4a), but that are not affected in either individual knockout. This is consistent, among other possibilities, with compensatory (for example, 'OR' gate logic) regulation20. Examination of ATAC-seq (assay for transposase-accessible chromatin using sequencing)21,22 and ChIP- seq data for PD-1, TIM-3, LAG-3 and TIGIT shows that PRDM1 and c-MAF can bind both overlapping and non-overlapping sites in the loci of these receptors and can synergistically trans-activate TIM-3 expression (Extended Data Fig. 7).

Mice with a T cell-specific deletion in both PRDM1 and c-MAF (PRDM1/c-MAF conditional double-knockout (cDKO)) showed normal development of CD4+ and CD8+ T cells in terms of frequency and expression of memory or activation markers, although the frequency of FOXP3+ Treg cells was increased (Extended Data Fig. 8a). CD4+ and CD8+ TILs from cDKO mice bearing B16F10 melanomas exhibited a near absence of PD-1, TIM-3, LAG-3, TIGIT, PDPN and PROCR expression (Fig. 4b, Extended Data Fig. 8b). Moreover, cDKO CD8+ TILs exhibited enhanced IL-2 and TNF production (Extended Data Fig. 8c). In contrast to singly deficient mice, cDKO mice showed significant control of B16F10 tumour growth despite the increased frequency of Treg cells (Fig. 4c). We addressed whether PRDM1 and c-MAF have a cell-intrinsic role in CD8+ and CD4+ T cells in controlling tumour growth by using an adoptive transfer model. Although CD8+ T cells from cDKO were able to inhibit tumour growth with decreased expression of co-inhibitory molecules, these effects were stronger when PRDM1 and c-MAF were lacking in both CD4+ and CD8+ T cells (Fig. 4d, Extended Data Fig. 8d). We examined the roles of PRDM-1 and c-MAF in tumour antigen-specific T cell responses using the MC38-OVA tumour model. We observed a significant reduction in tumour growth in mice receiving cDKO T cells as compared to mice receiving wild-type T cells (Extended Data Fig. 8e). We also observed an increase in ovalbumin (OVA)-specific T cells in the tumour draining lymph nodes and in OVA-specific IFN-γ- and TNF-producing CD8+ T cells in both the tumour infiltrate and the periphery in mice receiving double-knockout T cells (Fig. 4e, f, Extended Data Fig. 8f). Lastly, we observed an increase in CD8+Ki67+ T cells in the periphery of mice receiving double-knockout T cells (Fig. 4f).

We tested for non-additive effects between PRDM1 and c-MAF by using a binomial generalized linear model to compare the effect of single knockouts to the cDKO, and found that 149 out of 940 differentially expressed genes (adjusted P < 0.05, likelihood ratio test and false discovery rate (FDR) correction) between wild-type and cDKO CD8+ TILs have non-additive (that is, synergistic) effects (Extended Data Fig. 9, Methods).

Examination of the transcriptional signatures of cDKO CD8+ TILs showed significant overlap with those of CD8+TIM-3-PD-1- TILs (Fig. 4g, P = 2.8 × 10-7, one-sample Kolmogorov-Smirnov test; Extended Data Fig. 10a-c), suggesting that the loss of both c-MAF and PRDM1 increases the proportion of non-exhausted CD8+ effectors that exist normally in tumours. We scored the individual scRNA-seq profiles of CD8+ TILs for the cDKO 940 gene signature and found that the expression of the cDKO gene signature and the co-inhibitory gene module signature mark mutually exclusive populations of TILs (Extended Data Fig. 10e). The cDKO signature showed significant overlap with PD-1+CXCR5+CD8+ T cells, which may represent precursors for functional effectors in chronic lymphocytic choriomeningitis virus (LCMV) infection23 (Extended Data Fig. 10d, e, P = 1 × 10-13, one-sample Kolmogorov-Smirnov test). Furthermore, the IL-27RA-knockout TIL signature also showed significant overlap with this PD-1+CXCR5+CD8+ T cell signature (P < 2.2 × 10-16, one-sample Kolmogorov-Smirnov test; Fig. 2a, Extended Data Fig. 10e). Collectively, our data indicate that the loss of c-MAF and PRDM1 preferentially results in loss of the co-inhibitory gene module expression and acquisition of a more responsive effector T cell state.

In conclusion, we identified a co-inhibitory gene module, which is expressed in several settings of both CD4+ and CD8+ T cell nonresponsiveness, along with its transcriptional regulators. The discovery of this module provides a basis for the identification of novel co-inhibitory and co-stimulatory receptors that may have an important role in T cell regulation.

Online content

Any Methods, including any statements of data availability and Nature Research reporting summaries, along with any additional references and Source Data files, are available in the online version of the paper at https://doi.org/10.1038/s41586- 018-0206-z.

Received: 16 August 2016; Accepted: 27 April 2018;

Published online 13 June 2018.

Sidebar

Acknowledgements We thank M. Collins for discussions, D. Kozoriz, J. Xia and Z. Chen for technical advice, S. Riesenfeld for computational advice, N. Paul and J. Keegan for CyTOF, and L. Gaffney for artwork. This work was supported by grants from the National Institutes of Health, the American Cancer Society, the Melanoma Research Alliance, the Klarman Cell Observatory at the Broad Institute, and the Howard Hughes Medical Institute.

Author contributions N.C., A.M., P.R.B., A.C.A., O.R.-R., A.R. and V.K.K. designed the experiment; N.C., A.M., S.K., J.N., C.D.B., P.R.B., J.D.B. and A.R. developed analytical tools; N.C., A.M., T.K., N.A., J.N., N.D.M., M.S.K., C.W., H.Z., T.L., Y.E. and P.R.B. performed experiments; A.M. and M.S. performed computational analysis. N.C. and A.M. wrote the original draftof the paper and P.R.B., A.C.A., A.R. and V.K.K. reviewed and edited the paper; A.C.A., A.R. and V.K.K. supervised the project.

Competing interests A.C.A. is a member of the SAB for Potenza Therapeutics and Tizona Therapeutics. V.K.K. has an ownership interest and is a member of the SAB for Potenza Therapeutics and Tizona Therapeutics. A.C.A.'s and V.K.K.'s interests were reviewed and managed by the Brigham and Women's Hospital and Partners Healthcare in accordance with their conflict of interest policies. A.R. is an SAB member for Thermo Fisher and Syros Pharmaceuticals and is a consultant for Driver Group.

Additional information

Extended data is available for this paper at https://doi.org/10.1038/s41586- 018-0206-z.

Supplementary information is available for this paper at https://doi. org/10.1038/s41586-018-0206-z.

Reprints and permissions information is available at http://www.nature.com/ reprints.

Correspondence and requests for materials should be addressed to A.C.A. or A.R. or V.K.K.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Methods

Mice. C57BL/6 wild-type, IL-27RA-knockout, and Prdm1fl/flmice were obtained from the Jackson Laboratory. Maffl/fl, Pdpnfl/flmice and Procrd/d mice were previously described13,15,25. Pdpnfl/flmice were initially obtained from C. Buckley and crossed to CD4-Cre mice to obtain conditional deletion in T cells. CD4-Cre mice were purchased from Taconic. Prdm1fl/fland Maffl/flmice were crossed to CD4-Cre mice to generate doubly deficient T cell conditional knockout mice. All experiments were performed in accordance with the guidelines outlined by the Harvard Medical Area Standing Committee on Animals.

Tumour experiments. B16F10 melanoma cells (ATCC) (5 × 105) were implanted into the right flank of C57BL/6 mice. Tumour size was measured in two dimensions using a caliper. TILs were isolated by dissociating tumour tissue in the presence of 2.5 mg ml-1 collagenase D for 20 min before centrifugation on a discontinuous Percoll gradient (GE Healthcare). Isolated cells were then used in various assays of T cell function. For antigen specific analysis, we applied adoptive transfer tumour experiments using T cells from PRDM1/c-MAF cDKO mice, CD4+ or CD8+ T cells sorted from cDKO mice or littermate controls were transferred into RAG1- knockout mice at a 2:1 ratio (CD4: 1 million per mouse and CD8: 0.5 million per mouse) 2 days before subcutaneous injection of B16-OVA or MC38-OVA tumour. B16-OVA was a giftfrom K. Wucherpfennig, and MC38-OVA was a giftfrom M. Smyth. For adoptive transfer tumour experiments using T cells from Procrd/d mice, CD4+ T cells from wild-type and CD8+ T cells from wild-type or Procrd/d mice were isolated by cell sorting (BD FACS Aria) and transferred into RAGdeficient recipient mice at a 2:1 ratio (WT CD4+: 1 million per mouse and WT or Procrd/d CD8+: 0.5 million per mouse) 2 days before tumour implant. Although we did not use blinding or randomization experimental approaches, at least five animals of target gene knockout and control mice were used to adequately power biological validation experiments throughout the article. All mice used are C57BL/6 background, both male and female, 6-12 weeks of age, 15-25 g. Each experiment was performed using age- and sex-matched controls (Supplementary Table 5).

CyTOF. Antibodies were labelled using MaxPar Metal Labelling Kits (DVS) by The Longwood Medical Area CyTOF Antibody Resource and Core. In some experiments, TILs were enriched using Dynabeads FlowComp Mouse Pan T (CD90.2) Kit (Invitrogen). Cells were washed and resuspended in CyTOF PBS (PBS plus 0.05% sodium azide and 0.5% BSA) and stained viability marker Rhodium (DVS) following the cocktail of antibodies against cell-surface molecules for 30 min. Cells were washed again and resuspended in CyTOF PBS with 4% paraformaldehyde. After 10 min fixation, cells were washed and barcoded with Cell-ID intercalators (DVS). Before analysis, cells were resuspended in water with beads and loaded to the CyTOF Mass Cytometer (DVS). CyTOF data were recorded in dual-count according to Fluidigm's recommended settings that calibrated on the fly, combining pulse-count and intensity information. Data obtained as mass peaks for the channels are processed according to cell event selection criteria. These criteria include cell viability selection (Pt195), single-cell selection (Intercalator-Ir), and barcoding selection (Pt194 and Pt198) to identify single-cell events from wild-type and knockout TILs for further analysis.

To obtain clusters of cells similar in their protein expression patterns, cells were clustered using k-means algorithm. Optimal cluster number was estimated using the within groups sum of squared error (SSE) plot followed by gap statistics with bootstrapping and first standard error max method. These methods suggested 9 clusters as optimal in the multidimensional space. Applying k-means clustering with k = 9 on our CyTOF data, resulted in clear distinction between cluster 1 and 2 of the CD8+ TILs and cluster 3 and 4 of the CD4+ TILs. This separation could be further visualized by two-dimensional nonlinear embedding of the protein expression profiles using t-SNE4. The t-SNE plot can then be overlaid by k-means clustering results to reflect a non-biased approach to the clusters or with intensity of the different markers.

Flow cytometry. Single-cell suspensions were stained with antibodies against CD4 (RM4-5), CD8 (53-6.7), PD-1 (RMP1-30), LAG-3 (C9B7W), TIGIT (GIGD7), TIM-3 (5D12), PROCR (eBio1560) and PDPN (8.1.1.) obtained from BioLegend. Fixable viability dye eF506 (eBioscience) was used to exclude dead cells. For intracytoplasmic cytokine (ICC) staining, cells were stimulated with phorbol myristate acetate (50 ng ml-1) and ionomycin (1 µg ml-1) or with OVA 323-339 peptide for antigen specific experiments. Permeabilized cells were then stained with antibodies against IL-2, TNF, IFN-γ or IL-10. All data were collected on a BD LSR II (BD Biosciences) and analysed with FlowJo software (Tree Star). In brief, data were analysed by FSC/SSC gates of starting cell population (the gating strategy is exemplified in Supplementary Fig. 1). Positive gates were set based on fluorescence minus one (FMO) controls in each setting for cell surface molecules and based on unstimulated sample for ICC staining.

In vitro T cell differentiation. CD4+ and CD8+ T cells were purified from spleen and lymph nodes using anti-CD4 microbeads and anti-CD8a microbeads (Miltenyi Biotech) then stained in PBS with 0.5% BSA for 15 min on ice with anti-CD4, anti-CD8, anti-CD62L, and anti-CD44 antibodies (all from Biolegend). Naive CD4+ or CD8+ CD62LhighCD44low T cells were sorted using the BD FACSAria cell sorter. Sorted cells were activated with plate-bound anti-CD3 (2 µg ml-1 for CD4 and 1 µg ml-1 for CD8) and anti-CD28 (2 µg ml-1) in the presence of recombinant mouse IL-27 (25 ng ml-1) (eBioscience). Cells were collected at various time points for RNA, intracellular cytokine staining and flow cytometry.

qPCR. Total RNA was extracted using RNeasy columns (Qiagen). Reverse transcription of mRNA was performed in a thermal cycler (Bio-Rad) using iScript cDNA Synthesis Kit (Bio-Rad). qPCR was performed in the Vii7 Real-Time PCR system (Applied Biosystems) using the primers for Taqman gene expression (Applied Biosystems). Data were normalized to the expression of Actb.

Nanostring RNA analysis, expression profiling of TILs. We analysed gene expression in CD8+ TILs from PRDM1 or c-MAF cKO mice bearing B16F10 melanoma collected on day 14 after tumour implantation, using a custom nanostring codeset of 397 genes representing both the IL-27-driven gene signature (245 genes) and the dysfunctional CD8+ TIL gene signature (245 genes) (Supplementary Table 3). Expression values were normalized by first adjusting each sample based on its relative value to all samples. This was followed by subtracting the calculated background (mean.2sd) from each sample with additional normalization by housekeeping geometric mean, in which housekeeping genes were defined as: Hprt, Gapdh, Actb and Tubb5. Differentially expressed genes were defined using the function that fits multiple linear models from the Bioconductor package limma in R26 with P < 0.05.

Microarray processing and analysis. Naive CD4+ and CD8+ T cells were isolated from wild-type or IL-27RA-knockout mice, and differentiated in vitro with or without IL-27. Cells were collected at 72 h for CD8+ and 96 h for CD4+, and Affymetrix GeneChip Mouse Genome 430 2.0 Arrays were used to measure the resulting mRNA levels at these time points. Individual CEL files were RMA normalized and merged to an expression matrix using the ExpressionFileCreator of GenePattern with default parameters27. Gene-specific intensities were then computed by taking for each gene j and sample i the maximal probe value observed for that gene. Samples were then transferred to log-space by taking log2(intensity).

Differentially expressed genes were annotated as genes with FDR-corrected ANOVA < 0.05 computed between the CD4 with or without IL-27 stimulation (CD4+ IL-27 and T helper type 0 (TH0)) subpopulations (1,202 genes). A total of 468 genes were differentially expressed between wild-type CD8+ T cells stimulated in the presence or absence of IL-27 (P < 0.05). Two hundred and thirtyfour genes were shared between these two differentially expressed gene lists (P = 2.25 × 10-157, hypergeometric test, background = 16,618 (union of genes expressed)). A list of 972 cell surface/cytokines genes of interest that include: cytokines, adhesion, aggregation, chemotaxis and other cell-surface molecules (Supplementary Table 4) composed using Gene Ontology (GO) annotation in Biomart was used to generate the gene subset in Fig. 2b and c.

RNA-seq gene expression profiling of tumour infiltrating cells. Tumourinfiltrating CD8+ T cells were isolated from wild-type, IL-27RA-knockout, PRDM1 cKO, c-MAF cKO, and PRDM1/c-MAF cDKO tumour-bearing mice via FACS sorting on a FACSAria (BD Biosciences). Tumour-infiltrating CD8+ T cells were processed using an adaptation of the SMART-Seq 2 protocol28, using 5 µl of lysate from bulk CD8+ T cells as the input for each sample during RNA cleanup via SPRI beads (approximately 2,000 cells lysed on average in RLT).

RNA-seq reads were aligned using Tophat29 (mm9) and RSEM-based quantification30 using known transcripts (mm9), followed by further processing using the Bioconductor package DESeq in R31. The data were normalized using TMM normalization. The TMM method estimates scale factors between samples that can be incorporated into currently used statistical methods for DE analysis. Postprocessing and statistical analysis was carried out in R30. Differentially expressed genes were defined using the differential expression pipeline on the raw counts with a single call to the function DESeq (adjusted P < 0.1). Heat map figures were generated using pheatmap package (https://CRAN.R-project.org/package= pheatmap).

scRNA-seq. CD4+ and CD8+ TILs from wild-type or IL-27RA-knockout mice bearing B16 melanomas were sorted into 96-well plates with 5 µl lysis buffer comprised of buffer TCL (Qiagen) plus 1% 2-mercaptoethanol (Sigma). Plates were then spun down for 1 min at 3,000 r.p.m. and immediately frozen at -80 °C. Cells were thawed and RNA was isolated with 2.2× RNAClean SPRI beads (Beckman Coulter Genomics) without final elution32. The beads were then air-dried and processed immediately for cDNA synthesis. Samples were then processed using the Smart-seq2 protocol33, with minor modifications applied to the reverse transcription step (M.S.K. and A.R., in preparation). This was followed by making a 25 µl reaction mix for each PCR and performing 21 cycles for cDNA amplification. Then 0.25 ng cDNA from each cell and 0.25 of the standard Illumina NexteraXT reaction volume were used in both the tagmentation and final PCR amplification steps. Finally, libraries were pooled and sequenced (50 × 25 paired-end reads) using a single kit on the NextSeq500 5 instrument. All CD4+ TIL (wild-type and IL-27RA-knockout) scRNA-seq data were generated as part of this study. CD8+ TIL single-cell data include wild-type CD8+ TIL data from Singer et al.3 and wild-type and IL-27RA-knockout CD8+ single-cell data generated as part of this study. scRNA-seq data preprocessing and expression. Initial preprocessing was performed as previously described3. In brief, paired reads were mapped to mouse annotation mm10 using Bowtie34 (allowing a maximum of one mismatch in seed alignment, and suppressing reads that had more than 10 valid alignments) and TPMs were computed using RSEM30, and log2(TPM + 1) values were used for subsequent analyses.

Next, we filtered out low-quality cells and cell doublets, maintaining for subsequent analysis the cells that had (1) 1,000-4,000 detected genes (defined by at least one mapped read), (2) at least 200,000 reads mapped to the transcriptome, and (3) at least 50% of the reads mapped to the transcriptome, ending with a total of 707 CD4+ and 825 CD8+ wild-type TILs and 376 CD4+ and 394 CD8+ IL-27RAknockout TILs. We restricted the genes considered in subsequent analyses to be the genes expressed at log2(TPM + 1) ≥ 2 in at least 20% of the cells.

After removal of low-quality cells, the data were normalized using quantile normalization followed by principal component analysis (PCA). Principal components 1-10 were chosen for subsequent analysis owing to a drop in the proportion of variance explained following principal component 10. We used t-SNE4 to visualize single-cells in a two-dimensional nonlinear embedding.

scRNA-seq clustering and differential expression analysis. For the coupled dataset of wild-type and IL-27RA-knockout TILs, we followed the analysis previously described35. We performed batch correction using ComBat36 and the batchcorrected expression matrix was then reduced using PCA; principal components 1-13 were chosen for subsequent analysis owing to a drop in the proportion of variance explained following principal component 13. Next, we cluster the cells based on their principal component scores using the Louvain-Jaccard method using 40 nearest neighbours, and the 13 principal components37,38; 11 clusters were detected. We then compared the composition of each cluster in terms of total number and percentage of wild-type and IL-27RA-knockout cells and found cluster 5 to be enriched for wild-type CD8 TILs cells (P = 0.0357, one sample t-test, Extended Data Fig. 3c, d). Projecting the IL-27 co-inhibitory gene module onto the scRNA-seq data highlighted clusters 4 and 5 (CD4 and CD8, respectively) (Extended Data Fig. 3e), further showing that in addition to the decrease in the expression of the co-inhibitory receptors: PD-1, TIM-3, LAG-3 and TIGIT (Fig. 1e), a significant decrease in the total IL-27 co-inhibitory gene module signature score is observed with lack of IL-27 signalling (P = 0.01, t-test, Extended Data Fig. 3f). Last, we searched for differentially expressed genes between clusters 4 and 5 and the rest of the clusters using a nonparametric binomial test35.

Signature analysis of other states of T cell non-responsiveness. Given that orthogonal approaches were used to generate the various signatures, we first addressed the robustness of each signature before the comparative analysis. First, to address some of the concerns regarding the definition of these signatures, we sub-sampled the genes in each of the signatures and observed the resulting changes by projection on the single-cell data. These changes were quantified by randomly selecting decreasing subsets of genes from each signature (100%, 90%...30%) and calculating the average silhouette width of the cells that scored high for the different generated signatures, based on Euclidian distance between the principal component values used to generate the t-SNE plot. This analysis shows that the signatures are relatively resilient to this procedure up to 60% of the original signature (Extended Data Fig. 4e).

Second, we calculated a signature P value per cell. The P value is calculated by generating random sets of signatures that are composed of genes with a similar average and variance expression levels as the original signature. This was followed by comparing the generated scores to the score obtained from the original signature. Cells that had a statistically significant score (adjusted P < 0.05) were marked by a plus symbol '+' (Extended Data Fig. 4f).

For viral exhaustion, a microarray dataset10 was downloaded, followed by RMA. A signature of viral exhaustion was defined as the genes that are differentially expressed between chronic and acute viral infection on day 15 and day 30. Genes were ranked based on a t-test statistic and fold change, each gene rank was then adjusted for multiple hypotheses testing using FDR. A threshold of fold change > 1.1 and FDR < 0.2 was applied.

For antigen-specific tolerance, data11 were downloaded. Two groups were defined, group 1 that includes the PBS and 0.008 µg treated samples (treatment number 1) versus group 2-80 µg (treatment number 5 and 6). After log2 transformation and quantile normalization, the Limma package was used to estimate the fold changes and standard errors by fitting a linear model for each gene for the assessment of differential expression. Genes with P < 0.05 were selected: 1,845 genes were upregulated of which 88 were defined as cytokine and cell surface molecules26,39,40.

For antigen non-specific tolerance, data12 were downloaded. Robust multi-array average (RMA) and quantile normalization were applied for background correction and normalization using the ExpressionFileCreator module of GenePatterns. Differentially expressed genes were defined using signal-to-noise ratio, following FDR correction. Differentially expressed genes were identified as genes having a FDR < 0.2 between mRNA expression profiles of naive CD4+ or CD4+GFP/IL-10+ T cells isolated from the spleen or central lymph nodes of B6NODF1IL-10-GFP mice following nasal treatment with anti-CD3, which attenuates the of progressive phase of experimental autoimmune encephalomyelitis.

For cancer, data3 were obtained. In brief, mRNA samples from CD8+TIM- 3-PD-1- (double negative) TILs, CD8+TIM-3-PD-1+(single positive), and CD8+TIM-3+PD-1+ (double positive) TILs were measured using Affimetrix GeneChip Mouse Genome 430 2.0 Arrays, expression values were RMA normalized, corrected for batch effects using ComBat36 and gene-specific intensities were then computed by using the maximal prob intensity per gene, values were transferred to log-space by taking log2(intensity). Differentially expressed genes were defined as genes with either an FDR-corrected t-test P value smaller or equal to 0.2 computed between the double negative and double positive subpopulations and a fold-change of at least 1.5 between the two subpopulations.

The IL-27 co-inhibitory gene module was defined as a union of the overlap between the IL-27-driven gene program (1,201 genes; see 'Microarray processing and analysis' section) and each of the four different states of T cell non-responsiveness mentioned above (272 genes, Supplementary Table 2).

For the IL-27RA-knockout signature, mRNA samples from FACS sorted CD8+ TILs from wild-type and IL-27RA-knockout mice bearing B16 melanomas were measured an adaptation of the SMART-Seq 2 protocol28 (see 'RNA expression profiling of tumour infiltrating cells' section). Differentially expressed genes were defined as genes with either an FDR-corrected t-test P value smaller or equal to 0.2 computed between the wild-type and IL-27RA-knockout, or a fold-change of at least 1.5 between the two subpopulations. IL-27RA-knockout signature was defined as 929 differentially expressed genes in IL-27RA-knockout CD8+ TILs compared to wild-type CD8+ TILs.

Single-cell gene signature computation. As an initial step, the data were scaled (z-score across each gene) to remove bias towards highly expressed genes. Given a gene signature (list of genes), a cell-specific signature score was computed by first sorting the normalized scaled gene expression values for each cell followed by summing up the indices (ranks) of the signature genes. For gene-signatures consisting of an upregulated and downregulated set of genes, two ranking scores were obtained separately, and the downregulated associated signature score was subtracted from the upregulated generated signature score. A contour plot was added on top of the t-SNE space, which takes into account only those cells that have a signature score above the mean to further emphasis the region of highly scored cells.

Network construction. Networks were generated using Cytoscape version 3.2.141. The network model is based on coupling in vitro RNA-seq gene expression data of naive CD8+ T cells from knockout (PRDM1 or c-MAF) and wild-type controls stimulated in the presence of IL-27 and previously published ChIP-seq data for c-MAF and predicted PRDM1-binding sites by motif scan. More specifically, differentially expressed genes between wild-type control and knockout were defined using the function that fits multiple linear models from the Bioconductor package limma in R26 with FDR <0.05. We used published c-MAF ChIP-seq data19 and PRDM1 ChIP-seq data16. In addition, potential PRDM1-binding sites were detected using FIMO (MEME suite; http://meme-suite.org/doc/fimo. html). Association to gene promoters was based on the following thresholds (upstream = 5,000, downstream = 500 of transcription start site) and the overlap with the co-inhibitory module was found to be significant (P = 0.009 hypergeometric, background of 20,000 genes). In the network presentation, we visualize all the genes that are part of the IL-27 inhibitory module (Fig. 4a, Extended Data Fig. 6e).

Reporting summary. Further information on experimental design is available in the Nature Research Reporting Summary linked to this paper.

Data availability. Sequence data that support the findings of this study have been deposited in the Gene Expression Omnibus (GEO) with the accession code GSE113968. All other data are available from the corresponding authors upon reasonable request.

References

1. Wherry, E. J. & Kurachi, M. Molecular and cellular insights into T cell exhaustion. Nat. Rev. Immunol. 15, 486-499 (2015).

2. Anderson, A. C., Joller, N. & Kuchroo, V. K. Lag-3, Tim-3, and TIGIT: co-inhibitory receptors with specialized functions in immune regulation. Immunity 44, 989-1004 (2016).

3. Singer, M. et al. A distinct gene module for dysfunction uncoupled from activation in tumor-infiltrating T cells. Cell 166, 1500-1511 (2016).

4. Maaten, L. H. G. Visualizing data using t-SNE. J. Mach. Learn. Res. 9, 2579-2605 (2008).

5. Fitzgerald, D. C. et al. Suppression of autoimmune inflammation of the central nervous system by interleukin 10 secreted by interleukin 27-stimulated T cells. Nat. Immunol. 8, 1372-1379 (2007).

6. Awasthi, A. et al. A dominant function for interleukin 27 in generating interleukin 10-producing anti-inflammatory T cells. Nat. Immunol. 8, 1380-1389 (2007).

7. Stumhofer, J. S. et al. Interleukins 27 and 6 induce STAT3-mediated T cell production of interleukin 10. Nat. Immunol. 8, 1363-1371 (2007).

8. Zhu, C. et al. An IL-27/NFIL3 signalling axis drives Tim-3 and IL-10 expression and T-cell dysfunction. Nat. Commun. 6, 6072 (2015).

9. Hirahara, K. et al. Interleukin-27 priming of T cells controls IL-17 production in trans via induction of the ligand PD-L1. Immunity 36, 1017-1030 (2012).

10. Doering, T. A. et al. Network analysis reveals centrally connected genes and pathways involved in CD8+ T cell exhaustion versus memory. Immunity 37, 1130-1144 (2012).

11. Burton, B. R. et al. Sequential transcriptional changes dictate safe and effective antigen-specific immunotherapy. Nat. Commun. 5, 4741 (2014).

12. Mayo, L. et al. IL-10-dependent Tr1 cells attenuate astrocyte activation and ameliorate chronic central nervous system inflammation. Brain 139, 1939-1957 (2016).

13. Castellino, F. J. et al. Mice with a severe deficiency of the endothelial protein C receptor gene develop, survive, and reproduce normally, and do not present with enhanced arterial thrombosis after challenge. Thromb. Haemost. 88, 462-472 (2002).

14. Sakuishi, K. et al. Targeting Tim-3 and PD-1 pathways to reverse T cell exhaustion and restore anti-tumor immunity. J. Exp. Med. 207, 2187-2194 (2010).

15. Peters, A. et al. Podoplanin negatively regulates CD4+ effector T cell responses. J. Clin. Invest. 125, 129-140 (2015).

16. Mackay, L. K. et al. Hobit and Blimp1 instruct a universal transcriptional program of tissue residency in lymphocytes. Science 352, 459-463 (2016).

17. Apetoh, L. et al. The aryl hydrocarbon receptor interacts with c-Maf to promote the differentiation of type 1 regulatory T cells induced by IL-27. Nat. Immunol. 11, 854-861 (2010).

18. Giordano, M. et al. Molecular profiling of CD8 T cells in autochthonous melanoma identifies Maf as driver of exhaustion. EMBO J. 34, 2042-2058 (2015).

19. Ciofani, M. et al. A validated regulatory network for Th17 cell specification. Cell 151, 289-303 (2012).

20. Capaldi, A. P. et al. Structure and function of a transcriptional network activated by the MAPK Hog1. Nat. Genet. 40, 1300-1306 (2008).

21. Karwacz, K. et al. Critical role of IRF1 and BATF in forming chromatin landscape during type 1 regulatory cell differentiation. Nat. Immunol. 18, 412-421 (2017).

22. Sen, D. R. et al. The epigenetic landscape of T cell exhaustion. Science 354, 1165-1169 (2016).

23. Im, S. J. et al. Defining CD8+ T cells that provide the proliferative burst after PD-1 therapy. Nature 537, 417-421 (2016).

24. Chen, L. & Flies, D. B. Molecular mechanisms of T cell co-stimulation and co-inhibition. Nat. Rev. Immunol. 13, 227-242 (2013).

25. Wende, H. et al. The transcription factor c-Maf controls touch receptor development and function. Science 335, 1373-1376 (2012).

26. Smyth, G. K. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 3, Article3 (2004).

27. Reich, M. et al. GenePattern 2.0. Nat. Genet. 38, 500-501 (2006).

28. Tirosh, I. et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science 352, 189-196 (2016).

29. Trapnell, C., Pachter, L. & Salzberg, S. L. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics 25, 1105-1111 (2009).

30. Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323 (2011).

31. Anders, S. & Huber, W. Differential expression analysis for sequence count data. Genome Biol. 11, R106 (2010).

32. Shalek, A. K. et al. Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells. Nature 498, 236-240 (2013).

33. Picelli, S. et al. Full-length RNA-seq from single cells using Smart-seq2. Nat. Protocols 9, 171-181 (2014).

34. Langmead, B., Trapnell, C., Pop, M. & Salzberg, S. L. Ultrafast and memoryefficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25 (2009).

35. Shekhar, K. et al. Comprehensive classification of retinal bipolar neurons by single-cell transcriptomics. Cell 166, 1308-1323 (2016).

36. Johnson, W. E., Li, C. & Rabinovic, A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8, 118-127 (2007).

37. Blondel, V. D., Guillaume, J. L., Lambiotte, R. & Lefebvre, E. Fast unfolding of communities in large networks. J. Stat. Mech. https://doi.org/10.1088/1742- 5468/2008/10/P10008 (2008).

38. Levine, J. H. et al. Data-driven phenotypic dissection of AML reveals progenitorlike cells that correlate with prognosis. Cell 162, 184-197 (2015).

39. Smyth, G. K. in Bioinformatics and Computational Biology Solutions using R and Bioconductor. Statistics for Biology and Health (eds. Gentleman, R. et al.) 397-420 (Springer, New York, 2005).

40. Davis, S. & Meltzer, P. S. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and Bi°Conductor. Bioinformatics 23, 1846-1847 (2007).

41. Lopes, C. T. et al. Cytoscape Web: an interactive web-based network browser. Bioinformatics 26, 2347-2348 (2010).