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

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 co-inhibitory 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-seq (scRNA-Seq) to analyze co-inhibitory and co-stimulatory receptor expression in 588 CD8+ and 316 CD4+ tumor-infiltrating lymphocytes (TILs) from B16F10 melanoma[3]. We found that PD-1, Tim-3, Lag-3, CTLA-4, 4–1BB, and TIGIT strongly co-vary in CD8+ TILs. CD4+ TILs showed a similar pattern with the additional co-expression of ICOS, GITR, and OX40 (). Single-cell mass cytometry (CyTOF) confirmed the surface co-expression of these receptors (). Expression of PD-1, Lag-3, Tim-3, and TIGIT was tightly correlated on both CD8+ and CD4+ TILs (). Clustering analysis (t-SNE4, Methods) showed two groups of CD8+ TILs (clusters 1 and 2) () where PD-1, Lag-3, Tim-3, and TIGIT were mainly expressed in cluster 1 cells () as were LILRB4 (), and co-stimulatory receptors of the TNF-receptor family, 4–1BB, OX-40, and GITR. In contrast, ICOS and CD226 were less restricted to cluster 1 (). We further observed two discrete clusters of CD4+ TILs (clusters 3 and 4) wherein PD-1, Tim-3, Lag-3, and TIGIT co-expression was restricted to cluster 3 ().

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 autoimmunity[5], induces IL-10-secreting Type 1 regulatory (Tr1) cells[6,7], and induces expression of Tim-3 and PD-L1 on CD4+ and CD8+ T cells[8,9]. Activation of CD4+ and CD8+ T cells in the presence of IL-27 induced Tim-3 (Havcr2), Lag-3, and TIGIT at mRNA () and protein levels (). Expression of Tim-3, Lag-3, and TIGIT was reduced in IL-27R-deficient T cells, whereas PD-1 (Pdcd1) expression was unaffected by IL-27 in vitro ().

CyTOF analysis showed that loss of IL-27ra resulted in loss of cells in cluster 1 of CD8+ TILs and cluster 3 of CD4+ TILs ( p-value= 5×10−23 and 6.8×10−7 for CD8+ and CD4+ respectively, hypergeometric test, ), 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 wasn’t dependent on IL-27 in vitro, it was dependent on IL-27R signaling in vivo. In line with the induction of IL-10 by IL-27[5−7], we observed reduced IL-10 in IL27ra KO CD8+ TILs ().

scRNA-seq of CD8+ and CD4+ TILs from WT and IL27ra KO mice () revealed distinct clusters of CD8+ (cluster 5) and CD4+ (cluster 4) TILs that highly expressed the co-inhibitory receptors PD-1, Tim-3, Lag-3, and TIGIT. Expression of these genes was decreased in CD8+ TILs from IL27ra KO mice, while only Tim-3 and Lag-3 were decreased in CD4+ TILs from IL27ra KO mice (). 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 analyzed 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 (). 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 non-responsiveness: CD8+ T cell exhaustion in both cancer[3] and chronic viral infection[10] and antigen-specific[11] or non-specific (anti-CD3 antibody[12]) CD4+ T cell tolerance. We found significant overlap with all of these signatures (Methods, ).

Projection of the IL-27/CD8+ cancer T cell exhaustion overlap signature onto the single-cell profiles of CD8+ TILs marked a distinct subset of cells (, 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 ( panels II-IV). The transcriptional program induced in IL27ra KO TILs was active in a complimentary subset of TILs (Methods, panel V). 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 ( panel VI). From these analyses, we identified a co-inhibitory gene module (272 genes) that is shared across multiple states of T cell non-responsiveness (). Within this module, we identified a set of 57 genes encoding cell surface receptors and cytokines, including Tim-3, Lag-3, TIGIT, and IL10 (), which we further stratified by their expression in cancer and chronic viral infections (). Two surface molecules, Procr (protein C receptor) and Pdpn (podoplanin) were highly expressed in the setting of cancer (). Activation of naïve CD4+ and CD8+ T cells in vitro in the presence of IL-27 induced the expression of Procr and Pdpn (). In vivo, Procr and Pdpn exhibited IL-27 dependent co-expression with PD-1 and Tim-3 on CD8+ TILs ().

Procr+ CD8+ TILs exhibited an exhausted phenotype, producing less TNFα and IL-2 and more IL-10 than Procr- CD8+ TILs (). Growth of B16F10 melanoma was inhibited in Procr hypomorph (Procrd/d)[13] mice (), and Procrd/d CD8+ TILs mice exhibited enhanced TNFα production, but no difference in IL-2, IFN-γ, or IL-10 (). Procrd/d TILs exhibited a decreased frequency of Tim-3hiPD-1hi CD8+ T cells suggesting that Procr signaling promotes a severely exhausted phenotype in CD8+ T cells[14] (). Adoptive transfer of CD8+ T cells lacking Procr revealed a T cell specific role for Procr in constraining tumor growth ().

Although Pdpn can limit CD4+ T cell survival in inflamed tissues[15], its role in T cell exhaustion is unknown. We observed a significant delay in B16F10 tumor growth in mice with Pdpn deficiency in T cells (Pdpn cKO) (). Pdpn-deficient CD8+ TILs exhibited enhanced TNFα production but no significant difference in IL-2, IFN-γ, or IL-10 (). The frequency of Tim-3hiPD-1hi CD8+ TILs was decreased, indicating a reduced accumulation of T cells with a severely exhausted phenotype in Pdpn cKO[14] (). Consistent with previous data[15], 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 tumor microenvironment ().

We identified the transcription factor (TF) Prdm1 as a candidate regulator of the co-inhibitory module. Prdm1 is induced in vitro by IL-27 in CD4+ and CD8+ T cells (), is enriched in TILs with high expression of the IL-27 co-inhibitory module ( and and Methods), and is overexpressed in exhausted CD8+ TILs (p-value= 0.0004, t-test, ). Network analysis based on profiling of naïve CD8+ T cells from mice with a T cell specific deletion of Prdm1 (Prdm1 cKO) stimulated with IL-27, showed that Prdm1 regulates multiple genes in the IL-27 co-inhibitory module (, p-value= 2.32×10−12; hypergeometric test; Methods). This was further supported by Prdm1 Chip-seq data[16] (p-value= 2.9×10−8 respectively, hypergeometric test; ).

CD8+ TILs from B16F10 tumor-bearing Prdm1 cKO mice expressed lower levels of Tim-3, PD-1, and Procr (); however, there was no difference in tumor growth compared to wild type (WT) controls (), indicating that the reduction of co-inhibitory receptor expression in Prdm1 cKO mice was insufficient to promote effective anti-tumor immunity. We therefore examined whether other TFs may regulate the co-inhibitory module and compensate for the absence of Prdm1. We analyzed 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; ). We found that only a few genes were upregulated in Prdm-1 cKO CD8+ T cells, including one TF, c-Maf (p-value < 0.05) (). Indeed, c-Maf is induced by IL-27, is co-expressed with Prdm1 in T cells upon IL-27 stimulation (), and can regulate IL-10 expression[17] and T cell exhaustion[18]. Additionally, many genes (226 genes, p-value 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 regions[19].

CD8+ TILs from c-Maf cKO mice exhibited decreased expression of multiple co-inhibitory receptors (). Interestingly, Prdm1 and c-Maf each impacted co-inhibitory receptor expression only partially (). As in the Prdm1 cKO mice, c-Maf cKO mice did not show any differences in tumor growth relative to controls (). Notably, Prdm1 expression in c-Maf cKO TILs was similar to that in WT TILs, indicating that Prdm1 might drive 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); therefore we examined whether they shared targets. We combined the network analysis for Prdm1 () with c-Maf ChIP-seq data[19] and c-Maf targets (Methods). We observed 121 genes in the co-inhibitory module that are affected (RNAseq) or have a direct binding event (ChIP-Seq) for both Prdm1 and c-Maf (), but that are not affected in either individual knockout. This is consistent, among other possibilities, with compensatory (e.g., “OR”) regulation[20]. Examination of ATACseq[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 ().

Mice with a T cell specific deletion in both Prdm1 and c-Maf (Prdm1/c-Maf cDKO) showed normal development of CD4+ and CD8+ T cells in terms of frequency and expression of memory/activation markers, although the frequency of Foxp3+ Treg was increased (). 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 (). Moreover, cDKO CD8+ TILs exhibited enhanced IL-2 and TNFα production (). In contrast to singly deficient mice, cDKO mice showed significant control of B16F10 tumor growth despite the increased frequency of Treg (). We addressed whether Prdm1 and c-Maf play a cell-intrinsic role in CD8+ and CD4+ T cells in controlling tumor growth by using an adoptive transfer model. Although CD8+ T cells from cDKO were able to inhibit tumor 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 (). We examined the roles of Prdm-1 and c-Maf in tumor antigen-specific T cell responses using the MC38-OVA tumor model. We observed a significant reduction in tumor growth in mice receiving cDKO T cells as compared to mice receiving WT T cells (). We also observed an increase in Ova-specific T cells in the tumor draining lymph nodes and in OVA-specific IFN-γ and TNF-α producing CD8+ T cells in both the tumor infiltrate and in the periphery in mice receiving DKO T cells (). Lastly, we observed an increase in CD8+ Ki67+ T cells in the periphery of mice receiving DKO T cells ().

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 (adj. p-value<0.05, likelihood ratio test and FDR correction) between WT and cDKO CD8+ TILs have non-additive (i.e. synergistic) effects ().

Examination of the transcriptional signatures of cDKO CD8+ TILs showed significant overlap with those of CD8+ Tim-3-PD-1- TILs ( p-value = 2.8×10−7 one-sample Kolmogorov-Smirnov test, , p-value=0.008), suggesting that loss of both c-Maf and Prdm1 increases the proportion of non-exhausted CD8+ effectors that exist normally in tumors. We scored the individual scRNA-seq profiles of CD8+ TILs for the cDKO 940 gene signature and found that expression of the cDKO gene signature and the co-inhibitory gene module signature mark mutually exclusive populations of TILs (). The cDKO signature showed significant overlap with PD-1+CXCR5+CD8+ T cells, which may represent precursors for functional effectors in chronic LCMV infection[23] (, p-value = 1×10−13 one-sample Kolmogorov-Smirnov test). Furthermore, the IL27ra KO TILs signature also showed significant overlap with this PD-1+CXCR5+CD8+ T cell signature (p-value < 2.2×10−16 one-sample Kolmogorov-Smirnov test, ). Collectively, our data indicate that 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 multiple settings of both CD4+ and CD8+ T cell non-responsiveness, 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 play an important role in T cell regulation.

Methods

Mice

C57BL/6 wild-type (WT), IL27ra KO, and Prdm1 fl/fl mice were obtained from the Jackson Laboratory (Bar Harbor, ME). c-Maf fl/fl, Pdpn fl/fl mice and Procr delta/delta mice were previously described[13,15,26]. Pdpn fl/fl mice were initially obtained from Christopher Buckley (University of Birmingham, Birmingham, UK) and crossed to CD4Cre mice to obtain conditional deletion in T cell. CD4Cre mice were purchased from Taconic (Hudson, NY). Prdm1 fl/fl and c-Maf fl/fl mice were crossed to CD4Cre mice to generate doubly deficient T cell conditional knockout mice. All experiments were performed in accordance to the guidelines outlined by the Harvard Medical Area Standing Committee on Animals (Boston, MA).

Tumor Experiments

5 × 105 B16F10 melanoma cells (ATCC) were implanted into the right flank of C57BL/6 mice. Tumor size was measured in two dimensions using a caliper. TILs were isolated by dissociating tumor tissue in the presence of 2.5 mg/ml 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 tumor 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 KO mice at a 2:1 ratio (CD4: 1 million/mouse and CD8: 0.5 million/mouse) 2 days before subcutaneous injection of B16-OVA or MC38-OVA tumor. B16-Ova was kind gift from Kai Wucherpfennig (Dana-Farber Cancer Institute, Boston, MA) and MC38-Ova was kind gift from Mark Smyth (QIMR Berghofer, Queensland Institute of Medical Research, Brisbane Australia). For adoptive transfer tumor experiments using T cells from Procrd/d mice, CD4+ T cells from WT and CD8+ T cells from WT or Procrd/d mice were isolated by cell sorting (BD FACS Aria) and transferred into Rag deficient recipient mice at a 2:1 ratio (WT CD4+: 1 million/mouse and WT or Procrd/d CD8+: 0.5 million/mouse) 2 days before tumor implant. Although we did not blinding or randomization, at least 5 animals of target gene knock out 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–25g. Each experiment was performed using age, sex matched controls ().

CyTOF

Antibodies were labeled using MaxPar® Metal Labeling 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 + 0.05% sodium azide + 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 WT TILs and KO 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 SE 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 non-linear embedding of the protein expression profiles using t-stochastic neighborhood embedding (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), and Tim-3 (5D12), Procr (eBio1560), and Pdpn (8.1.1.) were obtained from BioLegend (San Diego, CA). Fixable viability dye eF506 (eBioscience) was used to exclude dead cells. For intra-cytoplasmic cytokine staining, cells were stimulated with (PMA) (50ng/ml, Sigma-Aldrich, MO), ionomycin (1µg/ml, Sigma-Aldrich, MO). 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 analyzed with FlowJo software (Tree Star).

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, CA). Naïve 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 for CD4 and 1μg/ml for CD8) and anti-CD28 (2μg/ml) in the presence of rmIL-27 (25ng/ml) (eBioscience). Cells were harvested at various time points for RNA, intracellular cytokine staining, and flow cytometry.

Real-time PCR

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). Real-time PCR was performed in the Vii7™ Real-Time PCR system (Applied Biosystems) using the primers for Taqman gene expression (Applied Biosystems). Data was normalized to the expression of ACTB.

Nanostring RNA analysis

Expression profiling of TILs.

We analyzed gene expression in CD8+ TILs from Prdm1 or c-Maf cKO mice bearing B16F10 melanoma collected on day 14 after tumor implantation, using a custom nanostring code-set of 397 genes representing both the IL-27-driven gene signature (245 genes) and the dysfunctional CD8+ TIL gene signature (245 genes) (). 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, where housekeeping genes were defined as: Hprt, Gapdh, Actin and Tubb5. Differentially expressed genes were defined using the function that fits multiple linear models from the Bioconductor package limma in R[27] with p-value<0.05.

Microarray processing and analysis

Naïve CD4+ and CD8+ T cells were isolated from WT or IL27ra KO mice, and differentiated in vitro with or without IL-27. Cells were collected at 72 hours for CD8+ and 96 hours 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 parameters[28]. 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+ IL27 and Th0) subpopulations (1,202 genes). 468 genes were differentially expressed between WT CD8+ T cells stimulated in the presence or absence of IL-27 (p-value<0.05). 234 genes were shared between these two differentially expressed gene lists (p-value = 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 () composed using GO annotation in Biomart was used to generate the gene subset in and .

RNAseq gene expression profiling of tumor infiltrating cells

Tumor infiltrating CD8+ T cells were isolated from WT, IL27ra KO, Prdm1 cKO, c-Maf cKO, and Prdm1/c-Maf cDKO tumor bearing mice via FACS sorting on a FACSAria (BD Biosciences). Tumor infiltrating CD8+ T cells were processed using an adaptation of the SMART-Seq 2 protocol[29], using 5uL of lysate from bulk CD8+ T cells as the input for each sample during RNA cleanup via SPRI beads (~2,000 cells lysed on average in RLT).

RNA-seq reads were aligned using Tophat[30] (mm9) and RSEM-based quantification[31] using known transcripts (mm9), followed by further processing using the Bioconductor package DESeq in R[32]. The data was normalized using TMM normalization. The TMM method estimates scale factors between samples that can be incorporated into currently used statistical methods for DE analysis. Post-processing and statistical analysis was carried out in R[31]. Differentially expressed genes were defined using the differential expression pipeline on the raw counts with a single call to the function DESeq (adjusted p value<0.1). Heatmap figures were generated using pheatmap package[33].

Single-cell RNA-seq

CD4+ and CD8+ TILs from WT or IL27ra KO 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 one minute at 3000rpm and immediately frozen at −80°C. Cells were thawed and RNA was isolated with 2.2x RNAClean SPRI beads (Beckman Coulter Genomics) without final elution[34]. The beads were then air-dried and processed immediately for cDNA synthesis. Samples were then processed using the Smart-seq2 protocol[35], with minor modifications applied to the reverse transcription (RT) step (MSK and AR., 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 ¼ 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+ TILs (WT and IL27ra KO) single-cell RNA-seq data was generated as part of this study. CD8+ TILs single-cell data includes WT CD8+ TILs data from[3] and WT and IL27ra KO CD8+ single-cell data generated as part of this study.

Single-cell RNA-seq data preprocessing and expression

Initial preprocessing was performed as described in[3]. Briefly, paired reads were mapped to mouse annotation mm10 using Bowtie[36] (allowing a maximum of one mismatch in seed alignment, and suppressing reads that had more than 10 valid alignments) and TPMs were computed using RSEM[31], 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+ WT TILs and 376 CD4+ and 394 CD8+ IL27ra KO TILs. We restricted the genes considered in subsequent analyses to be the genes expressed at log2(TPM+1)≥2 in at least twenty percent of the cells.

After removal of low quality cells the data was normalized using quantile normalization followed by PCA analysis. PCs 1–10 were chosen for subsequent analysis due to a drop in the proportion of variance explained following PC10. We used tSNE[4] to visualize single-cells in a two-dimensional non-linear embedding.

Single-cell RNA-seq clustering and differential expression analysis

For the coupled dataset of WT and IL27ra KO TILs we followed the analysis described in [37]. We performed batch correction using ComBat[38] and the batch-corrected expression matrix was then reduced using PCA, PCs 1–13 were chosen for subsequent analysis due to a drop in the proportion of variance explained following PC13. Next, we cluster the cells based on their PC scores using the Louvain-Jaccard method using 40 nearest neighbors, and the 13 PCs[25,39]; 11 clusters were detected. We then compared the composition of each cluster in terms of total number and percentage of WT and IL27ra KO cells and found cluster 5 to be enriched for WT CD8 TILs cells (p-value= 0.0357, one sample t-test, ). Projecting the IL-27 co-inhibitory gene module onto the single-cell RNA-seq data highlighted clusters 4 and 5 (CD4 and CD8 respectively) (), further showing that in addition to the decrease in the expression of the co-inhibitory receptors: PD-1, Tim-3, Lag-3 and TIGIT (), a significant decrease in the total IL-27 co-inhibitory gene module signature score is observed with lack of IL-27 signaling (p-value=0.01, t-test, ). Last, we searched for differentially expressed genes between clusters 4/5 and the rest of the clusters using a nonparametric binomial test[37].

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 prior to 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 tSNE plot. This analysis shows that the signatures are relatively resilient to this procedure up to 60% of the original signature ().

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-value<0.05) were marked by ‘+’ ().

For viral exhaustion: Microarray dataset[10] 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 false discovery rate (FDR). A threshold of fold change>1.1 and FDR<0.2 was applied.

For antigen-specific tolerance: Data[11] 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 value < 0.05 were selected: 1,845 genes were upregulated of which 88 were defined as cytokine and cell surface molecules[27,40,41].

For antigen non-specific tolerance: Data[12] was 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 (SNR), following FDR correction. Differentially expressed genes were identified as genes having a FDR<0.2 between mRNA expression profiles of naïve CD4+ or CD4+ GFP/IL-10+ T cells isolated from the spleen or cLNs of B6NODF1IL10:GFP mice following nasal treatment with anti-CD3 which attenuates the of progressive phase of EAE.

For cancer: Data[3] was obtained. Briefly, mRNA samples from CD8+Tim-3-PD-1- (DN) TILs, CD8+Tim-3-PD-1+(SP), and CD8+Tim-3+PD-1+ (DP) TILs were measured using Affimetrix GeneChip Mouse Genome 430 2.0 Arrays, expression values were RMA normalized, corrected for batch effects using ComBat[38] 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 DN and DP 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 Methods: Microarray processing and analysis) and each of the four different states of T cell non-responsiveness mentioned above (272 genes, ).

For IL27ra KO signature: mRNA samples from FACS sorted CD8+ TILs from WT and IL27ra KO mice bearing B16 melanomas were measured an adaptation of the SMART-Seq 2 protocol[29] (see Method: RNA expression profiling of tumor infiltrating cells). 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 WT and IL27ra KO and a fold-change of at least 1.5 between the two subpopulations. IL27ra KO signature was defined as 929 differentially expressed genes in IL27ra KO CD8+ TILs compared to WT CD8+ TILs.

Single-cell gene signature computation

As an initial step, the data was 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 down-regulated associated signature score was subtracted from the up-regulated generated signature score. A contour plot was added on top of the tSNE 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.1[42]. The network model is based on coupling in vitro RNAseq gene expression data of naïve CD8+ T cells from KO (Prdm1 or c-Maf) and WT 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 WT control and KO were defined using the function that fits multiple linear models from the Bioconductor package limma in R[27] with FDR<0.05. We used published c-Maf ChIP-seq data[19] and and Prdm1 ChiP-seq data[16]. 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=5000, downstream=500 of TSS) and the overlap with the co-inhibitory module was found to be significant (p-value= 0.009 hyper geometric, background of 20,000 genes). In the network presentation, we visualize all the genes that are part of the IL-27 inhibitory module ().

CyTOF analysis of co-inhibitory and co-stimulatory receptor co-expression in TILs.

a) TILs were harvested from B16F10 melanoma tumor-bearing WT and IL27ra KO mice from Fig. 1b and analyzed using CyTOF (5000 cells from each). CyTOF data were analyzed using vi-SNE. Applying k-means clustering with (k=9) on the CyTOF data resulted in clear distinction between clusters 1, 2, 3 and 4. Polygons indicating clusters 1, 2 (in CD8+ T cells), 3 and 4 (in CD4+ T cells) are shown. Individual panels show expression of the indicated markers. b) Pie charts show the distribution of WT or IL27ra KO CD8+ and CD4+ TILs in clusters 1 and 2 (C1 and C2) of CD8+ TILs and clusters 3 and 4 (C3 and C4) of CD4+ TILs as defined in . c) Independent data of WT and IL27ra KO TILs samples from that shown in (5000 cells from each). Applying k-means clustering with (k=7) on the CyTOF data resulted in clear distinction between clusters 1, 2, 3 and 4. Polygons indicating clusters 1, 2 (in CD8+ T cells), 3 and 4 (in CD4+ T cells) are shown. d) vi-SNE plot highlighting the distribution of cells from WT (blue) and IL27ra KO (red) in CD8+ TILs clusters 1 and 2 and CD4+ TILs clusters 3 and 4. Pie charts show the distribution of WT or IL27ra KO CD8+ and CD4+ TILs in each cluster.

Figure 1.

Multiple co-inhibitory receptors are expressed as a module on CD4+ and CD8+ T cells

a) CD4+ and CD8+ tumor-infiltrating lymphocytes (TILs) were harvested from WT mice bearing B16F10 melanoma tumors. Top panels, co-expression analysis of co-inhibitory and co-stimulatory receptor mRNA expression as determined by single-cell RNA-seq for 316 CD4+ and 588 CD8+ TILs. Bottom panels, protein expression by CyTOF for 23,656 CD4+ and 36,486 CD8+ TILs. Spearman correlation, followed by dendrogram ordering of the matrix using Euclidian distance is shown. Data are from biologically independent experiments. b) TILs from WT mice bearing B16F10 melanoma were analyzed using CyTOF with a custom panel of antibodies against co-inhibitory and co-stimulatory cell surface receptors[2,24] (). Data were analyzed using vi-SNE. Polygons indicating clusters 1, 2 (in CD8+ T cells), 3 and 4 (in CD4+ T cells) are shown. Individual panels show expression of the indicated markers. c) Naïve T cells from either wild type (WT) or IL-27ra deficient (IL27ra KO) mice were stimulated with anti-CD3/CD28 in the presence or absence of IL-27. Indicated co-inhibitory receptors expression was examined by real-time PCR (qPCR) at 96hr (CD4) and 72hr (CD8). Data are from biologically independent animals. mean + s.e.m is shown. d) vi-SNE plot showing WT (red) and IL27ra KO (blue) cells. e) ScRNA-seq of TILs from mice bearing B16F10 melanoma. Data were analyzed using t-SNE. Polygons indicating cluster 4 (in CD4+ T cells, orange) and cluster 5 (in CD8+ T cells, blue) are shown. Individual panels show expression of the indicated markers. Bar graphs show the mean signal intensity for indicated co-inhibitory receptors from WT (CD4+ (n=849); CD8+ (n=1752)) and IL27ra KO (CD4+ (n=628); CD8+ (n=541)) TILs for CyTOF (d) or WT (CD4+ (n=707); CD8+ (n=825)) and IL27ra KO (CD4+ (n=376); CD8+ (n=394)) TILs for ScRNA-seq (e). Error bars indicate s.e.m. and *p < 0.05, **p < 0.01, ***p < 0.001; two-sided t-test.

IL-27 induces multiple co-inhibitory receptors on CD4+ and CD8+ T cells.

a) Naïve T cells from WT or IL27ra KO mice were stimulated in vitro with anti-CD3/CD28 in the presence or absence of IL-27. Expression of co-inhibitory receptors was determined by flow cytometry. Representative data of 3 biologically independent experiments are shown. b) Expression of PD-1, Tim-3, Lag-3, TIGIT, and IL-10 on CD8+ TILs obtained from WT and IL27ra KO mice bearing B16F10 melanoma was determined by flow cytometry. Thy1.1-IL-10 reporter mice crossed with WT and IL27ra KO mice were used for IL-10 expression analysis. Representative data of 3 biologically independent experiments are shown.

Single-cell RNA-seq expression analysis of WT and IL27ra KO TILs.

a) TILs were harvested from B16F10 melanoma tumor-bearing WT (707 and 825 for CD4+ and CD8+ respectively) and IL27ra KO (376 and 394 for CD4+ and CD8+ respectively) mice as in . t-SNE plot shows the presence of WT and IL27ra KO CD4+ and CD8+ TILs as indicated. b) Clustering using the Louvain-Jaccard method (40 nearest neighbors and 13 principal components[25]). c) The composition of each cluster in terms of total number (c) and percentage (d) of WT (red) and IL27ra KO (blue) cells. P-values (*p-value<0.05, **p-value<0.01, ***p-value<0.001) were calculated using one sample t-test. e) Projection of the IL-27 co-inhibitory module signature on the single-cell RNA-seq data. The contour plot marks the region of highly expressing cells by taking into account only those cells that have an expression value above the mean. f) Violin and box plots displaying the distribution of the IL-27 co-inhibitory module signature score compared between WT (72 and 98 for CD4+ and CD8+ respectively) and IL27ra KO (85 and 77 for CD4+ and CD8+ respectively) cells in clusters 4 and 5 (CD4+ and CD8+ respectively, *p-value=0.01, one-sided t-test. The lower and upper hinges in the boxplot correspond to the first and third quartiles and the horizontal line corresponds to the median).

Overlap of the IL-27-induced gene program with signatures from four states of T cell impairment/tolerance/dysfunction.

a) Pearson correlation between WT CD4+ and CD8+ T cells for the 1,201 genes that were differentially expressed between WT CD4+ T cells stimulated in the presence or absence of IL-27 (FDR<0.05). b) Expression profile of 118 differentially expressed genes (from (a)) encoding cell surface receptors and cytokines are shown as a heatmap. c) The IL-27-induced gene program (1,201 genes) was compared to T cell signatures obtained from four states of T cell non-responsiveness. Number of overlapping genes between the IL-27 gene program and each signature is depicted. P values (***p < 0.001) were determined by hypergeometric test: Nasal anti-CD3 – 4.7×10−21, Cancer – 1.2 ×10−33, antigen-specific tolerance – 4 ×10−14 and Viral exhaustion – 1.7 ×10-26. d) p-value statistics for the significance of the overlap between the IL-27-induced gene program (1,201) and genes induced in other states of T cell non-responsiveness using wilcoxGST and camera. e) Gene signatures from (c) were sub-sampled and projected onto the CD8+ single-cell TILs data. Changes were quantified by randomly selecting decreasing subsets of genes from each signature and calculating the average silhouette width of cells that scored high for the different generated signatures based on Euclidian distance between the principal component values used to generate the tSNE plot. The lower and upper hinges in the boxplot correspond to the first and third quartiles and the horizontal line corresponds to the median (Methods). f) Panels I-V, tSNE plots of the 588 CD8+ single-cell TILs (dots) harvested from WT mice bearing B16F10 melanoma tumor. Cells are colored in each panel by their signature score. The score reflects the relative average expression of the genes in the overlap of the IL-27 gene signature with the signatures for each of the indicated states of T cell non-responsiveness. Panel VI is a projection of a signature of the differentially expressed genes between CD8+ TILs from WT and IL27ra KO mice bearing B16 melanomas (Methods). The contour plot marks the region of highly scored cells by taking into account only those cells that have a signature score above the mean score. Cells that had a statistically significant score (adjusted p-value<0.05) were marked by ‘+’ (Methods).

Characterization of the role of Pdpn and Procr in CD8+ TILs

a) Pdpn and Procr protein and mRNA expression was determined in T cells from WT and IL27ra KO stimulated with anti-CD3/CD28 in the presence or absence of IL-27. CD4+ cells were analyzed at 96hr and CD8+ cells at 72hr. Data are representative flow cytometry and qPCR data from biologically independent animals. mean + s.e.m is shown. b) Representative flow cytometry data of 3 independent experiments showing Pdpn and Procr expression in PD-1+Tim-3+ CD8+ and PD-1-Tim-3- CD8+ TILs obtained from WT and IL27ra KO mice bearing B16F10 melanoma. c) TILs from WT mice bearing B16F10 melanoma were stimulated with PMA and Ionomycin. Cytokine production in Procr+ or Procr- CD8+ TILs is shown. Thy1.1-IL-10 reporter mice were used for IL-10 expression analysis. Data are from biologically independent animals. mean + s.e.m is shown. *p < 0.05; **p < 0.01, paired t-test. d) 5×105 CD8+ T cells from wild type or Procrd/d mice were transferred along with 1×106 wild type CD4+ T cells to Rag1 KO mice (N=5). On day 2, 5×105 B16F10 cells were implanted. Mean tumor size + s.e.m is shown. *P<0.05, repeated measures ANOVA, Sidak’s multiple comparisons test. e) TILs were obtained from WT and Pdpn cKO mice bearing B16F10 melanoma and stained for the expression of IL-7Ra. Representative flow cytometry data from 3 independent animals. f) Summary data of IL-7Ra expression are from biologically independent animals. mean + s.e.m is shown. *p < 0.05, one-sided t-test.

Prdm1 is a candidate regulator of the co-inhibitory module.

a) Log2 fold change RNA levels between naïve CD4+ or CD8+ T cells simulated with or without IL-27. Data are from two independent experiments. Shown are transcription factors that are part of the IL-27 co-inhibitory module (Differentially expressed transcription factors were annotated as genes with FDR-corrected ANOVA <0.05). b) Transcription factors that are both in the IL-27 co-inhibitory module and are also overexpressed in clusters 4 and 5 in the single-cell data (clusters that were enriched for the IL-27 signature, ). Differentially expressed genes between clusters 4/5 and the rest of the clusters were determined using binomcount.test (binomial distribution, Methods). Log effect corresponds to log proportion of expressing cells and p-value is calculated by the probability of finding n or more cells positive for the gene in clusters 4/5 given the fraction in the rest of the clusters. c) tSNE plot of Fig. 1e. showing the expression of Prdm1 in WT (707 and 825 for CD4+ and CD8+, respectively) and IL27ra KO (376 and 394 for CD4+ and CD8+, respectively) cells. d) Normalized RNA expression levels of Prdm1 in PD-1-Tim-3- (n=3) and PD-1+Tim-3+ (n=3) CD8+ TILs (mean + s.e. is shown, ***p = 0.0004, two-sided t-test). e) Network model based on RNAseq gene expression data of naïve CD8+ T cells from Prdm1fl/fl (WT) or CD4crePrdm1fl/fl (Prdm1 cKO) mice stimulated in the presence of IL-27 and actual binding events (ChIPseq) data for Prdm1[19]. Green arrows designate genes up-regulated by Prdm1, red arrows designate genes down-regulated by Prdm1, and dashed gray arrows mark binding events.

Genomic tracks surrounding the co-inhibitory molecules

Lag3 (a), Pd-1 (b), Tigit (c) and Tim-3 (d) with overlay of Chipseq data of Prdm1[16] and c-Maf[19] and ATACseq data of naïve CD4+ cells induced with IL27 for 72h and ATACseq data of CD8+ T cells 27 days following chronic viral infection[22]. Regions of binding sites common to both Prdm1 and c-Maf are indicated by the dotted rectangles. e) Luciferase activity in 293T cells transfected with pGL4.23 luciferase reporters for depicted enhancers of Tim-3 together with empty vector (control), constructs encoding Prdm1, c-Maf, or both. Firefly luciferase activity was measured 48h after transfection and is presented relative to constitutive Renilla luciferase activity.

Immune characterization of Prdm1 cKO, cMaf cKO, and Prdm1/c-Maf cDKO before and after tumor challenge.

a) Analysis of steady-state immune system in WT, c-Maf cKO, Prdm1 cKO, and Prdm1/c-Maf cDKO. Summary data for CD4, CD8, Foxp3, CD44, CD62L, and CD69 expression in spleen from WT, c-Maf cKO, Prdm1 cKO andPrdm1/c-Maf cDKO mice. Data are from biologically independent animals. mean + s.e.m is shown. *p < 0.05; **p <0.01; ****p < 0.0001, one-way ANOVA and Tukey’s multiple comparisons test. b) co-inhibitory receptor expression in CD4+ TILs from Prdm1/c-Maf cDKO mice. Top panels, representative flow cytometry data from 3 independent experiments for TILs from WT and Prdm1/c-Maf cDKO stained for PD-1, Tim-3, TIGIT, Pdpn, and Procr expression. Bottom panels show summary data. Data are from biologically independent animals. mean + s.e.m is shown *p < 0.05, two-sided t-test. c) Top panels, representative flow cytometry data from 3 independent experiments showing cytokine production from CD8+ TILs from WT and cDKO bearing B16F10 melanoma. Bottom panels, summary data. Data are from biologically independent animals. mean + s.e.m is shown. *p < 0.05, two-sided t-test. d) Co-inhibitory receptor expression on CD8+ TILs sorted from B16-OVA-bearing Rag1 KO mice that were transferred with Prdm1/c-Maf cDKO (n=4) or wild type (n=4) CD4+ and CD8+ T cells as indicated. Data are from biologically independent animals. mean + s.e.m is shown. *P<0.05, one-way ANOVA and Tukey’s multiple comparisons test. e) Rag1 KO mice were transferred with either wildtype or cDKO CD4+ and CD8+ (2:1 CD4:CD8 ratio) followed by subcutaneous injection of MC38-OVA. Mean tumor size + s.e.m is shown. ****P<0.0001, repeated measures ANOVA, Sidak’s multiple comparisons test. On Day 14 post tumor implantation mice were sacrificed and TILs, spleen and draining Lymph nodes were harvested. f) The frequency of antigen specific CD8+ T cells in the dLN of mice in (e).

Examination of additive and non-additive (synergistic) effects of Prdm1 and c-Maf.

a) A Heatmap showing all 940 DE genes between WT (n=5) and cDKO (Prdm1/c-Maf, n=4) and their expression in single KO (Prdm1 control n=7, Prdm1 KO n=3, cMaf control n=4 and cMaf KO n=3) mice. The red markings on the top of the heatmap indicate genes on whose expression the two knockouts have a statistically significant (p-value<0.05) non-additive effect in the cDKO (149 out of 940 DE genes). b) Volcano plot of the same analysis as in (a) for global gene expression. Genes whose expression in the two single knockouts have a statistically significant (p-value<0.05) non-additive effect in the cDKO (1144 out of 12,906 genes) and had abs (coefficient)>1 (779 out of 1144) are shown in orange.

Comparison of gene expression between Prdm1/c-Maf cDKO TILs and CD8+ TILs populations from wild type mice.

a) Barcode enrichment plot displaying two gene sets in a ranked gene list. The ranked gene list was defined as fold change in gene expression between Prdm1/c-Maf cDKO and WT CD8+ TILs. The three gene sets consist of differentially expressed genes between: PD-1+Tim-3+ CD8+ (DP, n=3) and PD-1-Tim-3- CD8+ (DN, n=3) TILs, PD-1+Tim-3+ CD8+ (DP) TILs and Memory CD8+ (n=3), and PD-1+Tim-3- CD8+ (SP, n=3) and PD-1-Tim-3- CD8+ (DN) TILs. b) This analysis was followed by four statistical tests (one-sample Kolmogorov-Smirnov test, mean-rank gene set test (wilcoxGST), hypergeometric, and competitive gene set test accounting for inter-gene correlation) for enrichment of these signatures in the cDKO expression profile. c) WT versus cDKO volcano plot. Green indicates genes that were up-regulated in the PD-1-Tim-3- CD8+ (DN) TILs and red indicates genes that were up-regulated in the PD-1+Tim-3+ CD8+ (DP) TILs. d) WT versus cDKO volcano plot. Red indicates genes that were up-regulated in PD-1+CXCR5+CD8+ T cells and green indicates genes that were up-regulated in PD-1+CXCR5-CD8+ T cells in chronic LCMV infection[23]. e). A tSNE plot of the 588 CD8+ TILs harvested from WT mice bearing B16F10 melanoma tumors, colored by the relative signature score for the co-inhibitory module (272 genes, ), the cDKO signature (shown in (g)), and the PD-1+CXCR5+CD8+ T cell signature from chronic virus infection[23]. The contour plot marks the region of highly scored cells by taking into account only those cells that have a signature score above the mean.