Single-cell Genomics Unveils Critical Regulators of Th17 cell Pathogenicity (original) (raw)

Cell. Author manuscript; available in PMC 2016 Dec 3.

Published in final edited form as:

PMCID: PMC4671824

NIHMSID: NIHMS736435

Jellert T. Gaublomme,1,2,3,* Nir Yosef,4,* Youjin Lee,1,5,* Rona S. Gertner,2,3 Li V. Yang,6 Chuan Wu,5 Pier Paolo Pandolfi,7 Tak Mak,8 Rahul Satija,1,9,10 Alex K. Shalek,1,11,12 Vijay K. Kuchroo,1,5,** Hongkun Park,1,2,3,#** and Aviv Regev1,13,#**

Jellert T. Gaublomme

1Broad Institute of MIT and Harvard, 415 Main Street, Cambridge, Massachusetts 02142, USA

2Department of Chemistry & Chemical Biology, Harvard University, 12 Oxford Street, Cambridge, Massachusetts 02138, USA

3Department of Physics, Harvard University, 17 Oxford Street, Cambridge, Massachusetts 02138, USA

Nir Yosef

4Department of Electrical Engineering and Computer Science and Center for Computational Biology, University of California, Berkeley 94720, USA

Youjin Lee

1Broad Institute of MIT and Harvard, 415 Main Street, Cambridge, Massachusetts 02142, USA

5Evergrande Center for Immunologic Diseases, Harvard Medical School and Brigham and Women’s Hospital, Harvard Institutes of Medicine, 77 Avenue Louis Pasteur, Boston, Massachusetts 02115, USA

Rona S. Gertner

2Department of Chemistry & Chemical Biology, Harvard University, 12 Oxford Street, Cambridge, Massachusetts 02138, USA

3Department of Physics, Harvard University, 17 Oxford Street, Cambridge, Massachusetts 02138, USA

Li V. Yang

6Department of Internal Medicine, Department of Anatomy and Cell Biology, Brody School of Medicine, East Carolina University, Greenville, North Carolina 27834, USA

Chuan Wu

5Evergrande Center for Immunologic Diseases, Harvard Medical School and Brigham and Women’s Hospital, Harvard Institutes of Medicine, 77 Avenue Louis Pasteur, Boston, Massachusetts 02115, USA

Pier Paolo Pandolfi

7Cancer Research Institute, Beth Israel Deaconess Cancer Center, Department of Medicine and Pathology, Beth Israel Deaconess Medical Center, Harvard Medical School, Boston, MA 02215, USA

Tak Mak

8Ontario Cancer Center, Princess Margaret Hospital, Toronto, ON M5G 2M9, Canada

Rahul Satija

1Broad Institute of MIT and Harvard, 415 Main Street, Cambridge, Massachusetts 02142, USA

9NewYork Genome Center, New York 10013, USA

10Center for genomics and systems biology, New York University, New York 10012, USA

Alex K. Shalek

1Broad Institute of MIT and Harvard, 415 Main Street, Cambridge, Massachusetts 02142, USA

11Institute for Medical Engineering & Science and Department of Chemistry, Massachusetts Institute of Technology, Cambridge, Massachusetts 02142, USA

12Ragon Institute of Massachusetts General Hospital, Massachusetts Institute of Technology, and Harvard University, Boston, Massachusetts 02139, USA

Vijay K. Kuchroo

1Broad Institute of MIT and Harvard, 415 Main Street, Cambridge, Massachusetts 02142, USA

5Evergrande Center for Immunologic Diseases, Harvard Medical School and Brigham and Women’s Hospital, Harvard Institutes of Medicine, 77 Avenue Louis Pasteur, Boston, Massachusetts 02115, USA

Hongkun Park

1Broad Institute of MIT and Harvard, 415 Main Street, Cambridge, Massachusetts 02142, USA

2Department of Chemistry & Chemical Biology, Harvard University, 12 Oxford Street, Cambridge, Massachusetts 02138, USA

3Department of Physics, Harvard University, 17 Oxford Street, Cambridge, Massachusetts 02138, USA

Aviv Regev

1Broad Institute of MIT and Harvard, 415 Main Street, Cambridge, Massachusetts 02142, USA

13Howard Hughes Medical Institute, Department of Biology, Massachusetts Institute of Technology, Cambridge, Massachusetts 02142, USA

1Broad Institute of MIT and Harvard, 415 Main Street, Cambridge, Massachusetts 02142, USA

2Department of Chemistry & Chemical Biology, Harvard University, 12 Oxford Street, Cambridge, Massachusetts 02138, USA

3Department of Physics, Harvard University, 17 Oxford Street, Cambridge, Massachusetts 02138, USA

4Department of Electrical Engineering and Computer Science and Center for Computational Biology, University of California, Berkeley 94720, USA

5Evergrande Center for Immunologic Diseases, Harvard Medical School and Brigham and Women’s Hospital, Harvard Institutes of Medicine, 77 Avenue Louis Pasteur, Boston, Massachusetts 02115, USA

6Department of Internal Medicine, Department of Anatomy and Cell Biology, Brody School of Medicine, East Carolina University, Greenville, North Carolina 27834, USA

7Cancer Research Institute, Beth Israel Deaconess Cancer Center, Department of Medicine and Pathology, Beth Israel Deaconess Medical Center, Harvard Medical School, Boston, MA 02215, USA

8Ontario Cancer Center, Princess Margaret Hospital, Toronto, ON M5G 2M9, Canada

9NewYork Genome Center, New York 10013, USA

10Center for genomics and systems biology, New York University, New York 10012, USA

11Institute for Medical Engineering & Science and Department of Chemistry, Massachusetts Institute of Technology, Cambridge, Massachusetts 02142, USA

12Ragon Institute of Massachusetts General Hospital, Massachusetts Institute of Technology, and Harvard University, Boston, Massachusetts 02139, USA

13Howard Hughes Medical Institute, Department of Biology, Massachusetts Institute of Technology, Cambridge, Massachusetts 02142, USA

*Co-first authors

**Co-senior authors

Supplementary Materials

1.

GUID: 5CEE2A56-D5F4-45BF-BCD1-EA23DCDE92B7

3.

GUID: 771F326F-E2B9-470D-AA11-EE1B57ABFA84

4.

GUID: B537AC62-F76A-4426-B503-C8873DBCCEAE

5.

GUID: F05F3B00-F66B-41AA-A059-116DEF3A6008

6.

GUID: 02B79845-C3E1-4BEB-9154-B48DA5927F1B

7.

GUID: 90079E00-876B-4DE8-86DF-BFE8A9ABEAFE

8.

GUID: DA1D0FAD-DB86-4DC2-A0E5-1ABC5BA660BF

9.

GUID: 610B3F12-89E8-49DB-9F71-2ED9D399BF4A

10.

GUID: D09A3A90-EB08-4DF3-93F4-9E5803CF885A

11.

GUID: D4C94355-8CC0-42AA-A5F4-AD87C8757939

12.

GUID: A259D7F3-6A39-4EDD-B3A4-6F6D23A0B48E

13.

GUID: BD3A5164-0228-4D61-BCC9-28166CCD59A3

14.

GUID: 3F3B7191-B67B-4174-8179-A3E4B3C4D5A2

15.

GUID: 41918A6B-99A0-4F58-993A-86868C21D404

16.

GUID: 5419A3E7-4565-4E4A-B96C-1C0197657FBC

2.

GUID: 9D71334D-29F9-47E2-B795-9865E76125DF

SUMMARY

Extensive cellular heterogeneity exists within specific immune-cell subtypes classified as a single lineage, but its molecular underpinnings are rarely characterized at a genomic scale. Here, we use single-cell RNA-seq to investigate the molecular mechanisms governing heterogeneity and pathogenicity of Th17 cells isolated from the central nervous system (CNS) and lymph nodes (LN) at the peak of autoimmune encephalomyelitis (EAE) or differentiated in vitro under either pathogenic or non-pathogenic polarization conditions. Computational analysis relates a spectrum of cellular states in vivo to in vitro differentiated Th17 cells, and unveils genes governing pathogenicity and disease susceptibility. Using knockout mice, we validate four new genes: Gpr65, Plzp, Toso and Cd5l (in a companion paper). Cellular heterogeneity thus informs Th17 function in autoimmunity, and can identify targets for selective suppression of pathogenic Th17 cells while potentially sparing non-pathogenic tissue-protective ones.

INTRODUCTION

The immune system strikes a balance between mounting proper responses to pathogens and avoiding autoimmune reactions. In particular, as part of the adaptive immune system pro-inflammatory IL-17-producing Th17 cells mediate clearance of fungal infections and other pathogens (Hernandez-Santos and Gaffen, 2012) and maintain mucosal barrier functions (Blaschitz and Raffatellu, 2010), but are also implicated in pathogenesis of autoimmunity (Korn et al., 2009).

Mirroring this functional diversity, in vitro polarized Th17 cells can either cause severe autoimmune responses upon adoptive transfer (‘pathogenic’, polarized with IL-1β+IL-6+IL-23) or have little or no effect in inducing autoimmune disease (‘non-pathogenic’, polarized with TGF-β1+IL-6) (Ghoreschi et al., 2010; Lee et al., 2012).

Analysis of these states has been limited however, by relying either on genomic profiling of cell populations, which cannot distinguish distinct states within them, or on tracking a few known markers by flow cytometry (Perfetto et al., 2004). Single-cell RNA-seq (Shalek et al., 2013; Shalek et al., 2014; Trapnell et al., 2014) opens the way for a more unbiased interrogation of cell states, including in limited in vivo samples.

Here, we use single-cell RNA-seq to show that cells isolated from the draining LNs and CNS at the peak of EAE exhibit diverse functional states, and relate them to a spectrum spanning from more regulatory to more pathogenic cells observed in Th17 cells polarized in vitro. Genes associated with these opposing states include both canonical known regulators and novel candidates. We validated four high-ranking candidates – Gpr65, Plzp, Toso and Cd5l (the latter in a companion study, Wang et al.) – with knockout mice, uncovering substantial effects on in vitro differentiation and in vivo EAE development.

RESULTS

RNA-seq profiling of single Th17 cells isolated in vivo and in vitro

We profiled the transcriptome of 976 Th17 cells (ultimately retaining 722 cells, below), either harvested in vivo or differentiated in vitro (Figure 1A and Table S1, Experimental Procedures). In vivo, we induced EAE by myelin oligodendrocyte glycoprotein (MOG35–55) immunization, harvested CD3+CD4+IL-17A/GFP+ cells from the draining LNs and CNS at the peak of disease and profiled them promptly. In vitro, we profiled CD4+ naïve T cells at 48h of activation under TGF-β1+IL-6 or IL-1β+IL-6+IL-23. We prepared mRNA SMART-Seq libraries using microfluidic chips (Fluidigm C1), followed by transposon-based library construction.

An external file that holds a picture, illustration, etc. Object name is nihms736435f1.jpg

Single-cell RNA-seq of Th17 cells in vivo and in vitro

(A) Experimental setup. (B–E) Quality of single-cell RNA-seq. Scatter plots (B–D) compare transcript expression (FPKM+1, log10) from the in vitro TGF-β1+IL-6 48hr condition, between two bulk population replicates (B), the ‘average’ of single-cell profile and a matched bulk population control (C), or two single cells (D). Histograms (E) depict the distributions of Pearson correlation coefficients (X axis) between single cells and their matched population control and between pairs of single cells. (F,G) Comparison to RNA Flow-FISH. (F) Expression distributions by RNA-seq and RNA Flow-FISH at 48h under the TGF-β1+IL-6 in vitro condition. Negative control: bacterial DapB gene. (G) Bright-field and fluorescence channel images of RNA Flow-FISH in negative (left) and positive (right) cells.

See also Figure S1, Table S1, related to Figure 1.

We removed 254 cells (~26%) by quality metrics (Supplemental Experimental Procedures) and we controlled for quantitative confounders and batch effects (Experimental Procedures, Figure S1A,B). We retained ~7,000 appreciably expressed genes (fragments per kilobase of exon per million (FPKM) > 10 in at least 20% of cells in each sample) for in vitro experiments and ~4,000 for in vivo ones. To account for expressed transcripts that are not detected (false negatives) due to the limitations of single-cell RNA-seq (Deng et al., 2014; Shalek et al., 2014), we down-weighted the contribution of less reliably measured transcripts (Figure S1C, Experimental Procedures). Following these filters, expression profiles tightly correlated between population replicates (Figure 1B), and between the average expression across single cells and the matching population profile (r ~ 0.65–0.93; Figure 1C, S1D, S2, and Table S1). However, we found substantial differences in expression between individual cells in the same condition (r ~ 0.45–0.75 Figure 1D, 1E, S1D), comparable to previous observations in other immune cells (Shalek et al., 2014). We validated the observed expression patterns for eight representative genes with flow RNA-fluorescence in situ hybridization (Supplemental Experimental Procedures) (Figure 1F, 1G, S1E). While most transcripts (e.g., Irf4, Batf, Actb) are expressed unimodally or nearly so (e.g., Rorc), some key transcripts (e.g., Il17a) exhibit a bimodal distribution, suggesting functional variation.

A functional annotation of single-cell heterogeneity shows that Th17 cells span a spectrum of states in vivo

To study the main sources of cellular variation in vivo, we used a principal component analysis (PCA, Figure 2A) followed by a novel analysis for functional annotation of the principal component (PC) space based on the single-cell expression of gene signatures of known T cell states (Figure 2B, Experimental Procedures), such that each signature is scored for its association with each PC. To identify transcription factors that may orchestrate this heterogeneity, we identified factors whose targets are strongly enriched (Fisher exact test, p<10−5) in genes that correlated with each PC (Pearson correlation, FDR < 0.05; Figures 2E,F, Table S3).

An external file that holds a picture, illustration, etc. Object name is nihms736435f2.jpg

Th17 cells span a progressive trajectory of states from the LN to the CNS

(A) PCA separates CNS-derived cells from LN-derived cells. Shown are 302 cells in the space of the first two PCs. Numbered circles: signatures that significantly correlate with PC1 or PC2 (p < 10−6, Table S2, Supplemental Experimental Procedures). (B) Functional annotation. Top to bottom: Gene signatures are defined from literature, and a signature score is calculated for each single cell. The Pearson correlation coefficient is calculated between the signature score and the PC loading, for each cell and PC, and plotted on the PCA plot (circled numbers in (A)) (Supplemental Experimental Procedures). (C) Five progressive Th17-cell states from the LN to the CNS. Plot as in (A), but with Voronoi cells, defined by signatures (colored circles, Table S2) characterizing the cells populating the extremities of PCA space: Five signature-specific subpopulations are marked. The self-renewing state was observed in two technical replicates of one of the two in vivo biological replicates, potentially due to differences in disease induction or progression. (D) Example genes that distinguish each sub-population. Cumulative distribution function (CDF) plots of expression for key selected genes. Dotted/solid line corresponds to CNS/LN cells respectively, where appropriate. (E,F) Transcription factors (nodes) whose targets are significantly enriched in PC2 (E) or PC1 (F). Nodes are sized proportionally to fold enrichment (Table S3) and colored according to the loading of the encoding gene in the respective PC (loadings were normalized to have zero mean and standard deviation of 1).

See also Figure S2–S4, Table S2–4, related to Figure 2.

The first PC (PC1) positively correlates with a recently defined effector vs. memory signature following viral infection (Crawford et al., 2014), and negatively correlates with a signature characterizing memory T cells (Wherry et al., 2007) (Figure 2A, Table S2). This suggests that cells span from an effector (high positive PC1 scores) to a memory (low negative PC1 scores) phenotypes. PC2 separates cells by their source of origin (CNS or LN) and correlates with a transition from a naïve-like state (negatively correlated with PC2; p<10−42, Figure 2A, Table S2) with low cell-cycle activity (FDR<5%) to a Th1-like effector or memory effector state (positively correlated with PC2, Figure 2A, p<10−21 and p<10−57, respectively).

A trajectory of progressing cell states from the LN to the CNS

To further explore the diversity of LN and CNS cells, we used five of the key signatures discovered by our functional annotation to define Voronoi regions that divide the PCA space into subpopulations of cells (Supplemental Experimental Procedures, Figure 2C, Table S2). We identified genes characterizing each group by differential expression (KS test, FDR<0.05; Table S4). For brevity, we assign new labels to these subpopulations (Th17 self-renewing, Th17/Th1-like effector, etc.), based on strong correlation with previous signatures and known genes (below); any such label may inevitably fall short of capturing the complex underlying biology.

Overall, the cells gradually progress through from a self-renewing-like state in the LN to a pre-Th1 effector-like phenotype in the LN and CNS, to a Th1-like effector state and a Th1-like memory state in the CNS, and finally a less functional state in the CNS. First, self-renewing-like Th17 cells in the LN (Figure 2C) are characterized by: (1) a signature of Wnt signaling (p<10−4, KS test comparing the signature score to all other sub-populations; Figure 2A) (Reya et al., 2003) and high expression of Tcf7 (Figure 2D), a key Wnt target and transcription factor regulating the stem cell-like state of Th17 cells (Muranski et al., 2011); (2) high Cd62l expression (Figure 2D),, a known naïve state marker (De Rosa et al., 2001); and (3) up-regulation (Figure 2D) of Cd27, a pro-survival gene lacking in short-lived T cells (Dolfi et al., 2008; Hendriks et al., 2000). Our analysis (Figure 2E) suggests that Etv6, Med12 and Zfx drive this self-renewing population (Discussion). Next, cells from the LN and CNS adopt similar (overlapping) cell states in the central region of our PCA plot (Figure 2C), reflecting effector Th17-like cells with a pre-Th1-like phenotype, characterized by induction of receptors for IFN (Ifngr2) and IL-18 (Il18r1, Figure 2D) (Holzer et al., 2013), and of chemokine receptors Cxcr6 (Figure 2D) (Aust et al., 2005) and Ccr2 (Figure 2D) (Mahad and Ransohoff, 2003), which may all poise the cells for recruitment to the CNS. In turn, IL-17A/GFP+ sorted cells acquire a Th17/Th1-like effector phenotype in the CNS (Figure 2C), with up-regulation (p<10−3, KS test, Table S4) of: _Ifn_-γ (Figure 2D), Rankl, (Tnfsf11) (Nakae et al., 2007) (Komatsu et al., 2014), and cell cycle genes (e.g., Geminin, Figure 2D), a strong correlation with a salt-induced pathogenic Th17 cell signature (Wu et al., 2013) (Figure 2A), and association with both canonical Th17 TFs (Stat3 and Hif1a) and Th1-associated factors, including Rel and Stat4 (Figure 2E), which are associated with EAE (Hilliard et al., 2002; Mo et al., 2008) or with human autoimmune disease (Gilmore and Gerondakis, 2011). These cells could either be stable ‘double producers’ or reflect Th17 plasticity into the Th1 lineage (Discussion). Next, Th1-like memory cells detected in the CNS (Figure 2C) correlate highly with both a memory phenotype (negative PC1) and a Th1-like phenotype (positive PC2), upregulating (p<10−3, KS test, Table S4) memory signature genes (e.g., Nur77, Samsn1, Il2ra, Il2rb, Tigit, Ifngr1 and Il1r2) and pro-inflammatory genes (Csf2 and Gpr65; Figure 2D), and associated with Hif1a regulation (Figure 2F), crucial for controlling human Th17 cells to become long-lived effector memory cells (Kryczek et al., 2011). While IL-2 is a growth factor for Th1 cells, IL-2 also affects Th17 differentiation and stability (Quintana et al., 2012). Finally, a few Th17 cells may even acquire a somewhat senescent-like phenotype in the CNS (negative PC1 and PC2 scores; Figure 2C), correlating with a signature comparing CD4 T cells at day 30 during a chronic infection to those in during acute infection (Table S2), and down regulation (p~10−2, KS test, Table S4) of some T-cell activation genes (Figure 2D, Table S4).

Further supporting the interpretation of gradual progression, most in vivo cells are maximally correlated with bulk profiles at 48–72h during Th17 cell differentiation in vitro (Yosef et al., 2013) (Figure S4B), but some cells, especially from the LN, correlate most strongly with earlier time points. Indeed, timepoint annotations positively correlate with PC2 (_r_~0.5; p<10−27, Figure 2A and Figure S4D). Finally, a population based signature comparing profiles of EAE Th17 cells to lamina propria lymphocyte (LPL) Th17 cells (Supplemental Experimental Procedures, Figure S3), which are known to assume a regulatory phenotype (Esplugues et al., 2011; O'Connor et al., 2009), correlates with PC1 in vivo (p < 10−26, Figure 2A, Table S2), indicating that EAE-derived Th17 cells adopt a stronger effector phenotype than gut-derived Th17 cells.

In vitro derived cells span a spectrum of pathogenicity states with similarities and distinctions from in vivo isolated cells

Given the limited cell availability from in vivo samples, and the fact that cells are obtained as a mixed “snapshot” of an asynchronous process, it is difficult to further characterize their distinct pathogenic potential. A complementary strategy is to profile Th17 cells differentiated in vitro, and compare in vivo and in vitro profiles.

We analyzed single-cell RNA-seq profiles of 420 Th17 cells derived under non-pathogenic (TGF-β1+IL-6, unsorted: 130 cells from 2 biological replicates and TGF-β1+IL-6, sorted for IL-17A/GFP+: 151 cells from 3 biological replicates) and pathogenic conditions (Il-1β+IL-6+IL-23, sorted for IL-17A/GFP+: 139 cells from 2 biological replicates) (Figure 3A).

An external file that holds a picture, illustration, etc. Object name is nihms736435f3.jpg

A spectrum of pathogenicity states in vitro

(A) PCA plot of Th17 cells differentiated in vitro. PC1 separates cells from most (left) to least (right) pathogenic, as indicated both by the differentiation condition (color code), and by the correlated signatures (numbered circles). (B–D) Key signatures related to pathogenicity. CDFs of the single-cell scores for key signatures for the three in vitro populations (colored as in A): (B) a signature distinguishing the in vivo Th17/Th1-like memory sub-population (blue in Figure 2C); (C) a signature distinguishing the in vivo Th17 self-renewing sub-population (green in Figure 2C); and (D) a signature of pathogenic Th17 cells (Lee et al., 2012). (E) CDFs of expression level (FPKM+1, log10) of Il10 for the three in vitro populations.

See also Table S2, related to Figure 3.

Using our functional annotation approach (Figure 2B), we find that i_n vitro_ differentiated Th17 cells vary strongly in a key pathogenicity signature (Lee et al., 2012), reflecting their respective conditions (Figure 3A,D). High pathogenicity scores were associated with IL-17A/GFP+ sorted cells polarized under the pathogenic condition (Figure 3A,D), whereas IL-17A/GFP+ sorted cells from non-pathogenic conditions correlate highly with expression of regulatory cytokines (e.g., IL-10), and their targets, which are barely detected in pathogenic cells (Figure 3E). There is a zone of overlap in cell states between the pathogenic and non-pathogenic conditions (Figure 3A) with cells polarized under the non-pathogenic (TGF-β1+IL-6) condition that were not sorted to be IL-17A/GFP+ spanning the broadest pathogenicity spectrum (Figure 3A,D). A signature from IL-23R−/− cells differentiated with IL-1β+IL-6+IL-23 (Y.L. and V.K.K, unpublished data) correlates highly with the more regulatory cells, confirming the role of the IL-23 pathway in pathogenicity (Figure 3A).

As expected, the most cells’ profiles correlate with bulk profiles at 48–50h (Yosef et al., 2013) (Figure S4A). The correlated time points match with variation along PC2 (Figure 3A), but not PC1, suggesting that pathogenicity is not a reflection of the cell’s position along the differentiation trajectory, but is an orthogonal aspect of cell state.

To relate the in vitro differentiated cells to their in vivo counterparts, we scored the in vitro cells for signatures of immune-related genes that characterize the in vivo identified subpopulations (Figure 2C; Figure 3B,C; Supplemental Experimental Procedures). Cells derived in the non-pathogenic conditions score higher for the Th17 self-renewing-like signature (p<10−10, KS test; Table S2 and Figure 3A,C), whereas those derived in pathogenic conditions resemble the Th17/Th-1 like memory phenotype identified in the CNS (p<10−3, KS test; Figure 3A,B and Table S2).

Co-variation with pro-inflammatory and regulatory modules in Th17 cells highlights novel candidate regulators

We analyzed each gene’s variation in expression across the unsorted cells from the TGF-β1+IL-6 differentiation condition. About 35% (2,252) of the detected genes are expressed in >90% of the cells (Figure 4A) with a unimodal distribution: these include housekeeping genes (p<10−10, hypergeometric test, Figure S5A&B), the Th17 signature cytokine _Il17f_, and transcription factors that are essential for Th17 differentiation (_e.g., Batf, Stat3, Rorc and Hif1a_). Bimodally expressed genes– with high expression in at least 20% of the cells and much lower (often undetectable) levels in the rest – include inflammatory and regulatory cytokines and their receptors (_e.g., Il17a, Il10, Il21, Ccl20, Il24, Il27ra_, Figure 4A). Expression variation may be more strongly related to pathogenicity than differentiation. Most (>75%) cells express pioneer and master transcription factors for the Th17 lineage (Rorc, Irf4, and Batf), but some also express transcripts encoding key genes from other T-cell lineages (e.g., Stat4 for Th1 cells, Ccr4 for Th2 cells), suggesting the presence of previously reported “hybrid” double-positive cells (Antebi et al., 2013), and/or reflecting our model of duality in the Th17 transcriptional network (Yosef et al., 2013). The expression of many key immune genes varies more than that of other transcripts with the same mean expression level (Figure S5C), even when only considering the expressing cells (Figure S5D), implying a greater degree of diversity in immune-gene regulation. Such patterns must be interpreted with caution, because some (e.g., Il17a, Il24 and Ccl20), but not all (e.g., Il9), of the transcripts with bi-modal patterns are lowly expressed and thus may be less reliably detected, and because transcription bursts coupled with transcript instability may lead to ‘random’ fluctuations.

An external file that holds a picture, illustration, etc. Object name is nihms736435f4.jpg

Modules of genes that co-vary with pro-inflammatory and regulatory genes across single cells

(A) Single-cell expression distribution of genes. The heat map shows for each gene (row) its expression distribution across single cells differentiated under the TGF-β1+IL-6 condition for 48h (without further IL-17A-based sorting). Color scale: proportion of cells expressing in each of the 17 expression bins (columns). Genes are sorted from more unimodal (top) to bimodal (bottom). (B) Modules co-varying with pro-inflammatory and regulatory genes. Heat map of the Spearman correlation coefficients between the single-cell expression levels of signature genes of pathogenic T cells (Lee et al., 2012) or of other CD4+ lineages (columns) and the single-cell expression of any other bimodally expressed gene (rows) in cells differentiated under the TGF-β1+IL-6 condition at 48h. Genes are clustered. (C) The modules co-varying with pro-inflammatory and regulatory genes distinguish key variation. Each cell (TGF-β1+IL-6, 48h) is colored by a signature score comparing the two co-variation modules. (D) Expression of key module genes. Each panel shows the PCA plot of (C) where cells are colored by an expression ranking score of a key gene, denoted on top. (E) A ranking of the top 100 candidate genes co-varying with pro-inflammatory or regulatory genes (out of 184; Table S5), sorting from high (left) to lower (right) ranking scores (bar chart, Supplemental Experimental Procedures).

See also Figure S5 and S6, Table S2&S5 related to Figure 4.

To overcome these challenges, we analyzed co-variation between transcripts across cells (Figure 4B), reasoning that if variation reflects distinct cell states, entire gene modules should robustly co-vary across cells. Focusing on significant co-variation (Spearman correlation; FDR<0.05) between bimodally expressed transcripts (expressed by less than 90% of cells; Figure 4B, rows, Supplemental Experimental Procedures) and a curated set of bimodally expressed immune response genes (Figure 4B, columns), we find two key transcript modules: one that co-varies with known pro-inflammatory Th17 cytokines, such as Il17a and Ccl20, and another that co-varies with known regulatory genes such as Il10, Il24, and Il9.

Using these modules as signatures to annotate the original in vitro cell states (Figure 3A & 4C), we find that a signature comparing the module co-varying with pro-inflammatory genes to the module co-varying with regulatory genes strongly correlates with the most pathogenic cells (Figure 4C,D). We find further support from additional signatures and analyses: (1) a negative correlation between PC1 scores and a curated pathogenicity signature (Lee et al., 2012) and a positive correlation between PC1 and a Tr1exTh17 cells vs. Th17 cell signature (Gagliani et al., 2015)(Figure 4C, Table S2); (2) correlation with the inflammatory CNS derived Th17 cells in vivo (Figure 2A); (3) enrichment of genes in the co-variation modules (rows of Figure 4B) for immune response genes (using MSigDB (Liberzon et al., 2011); Table S5), for genes generically associated with inflammatory bowel disease (Jostins et al., 2012) and rheumatoid arthritis (Okada et al., 2014) (p<10−6; hypergeometric test), and for genes up-regulated in cortical lesions derived from patients with progressive multiple sclerosis (Fischer et al., 2013) (P<0.02, hypergeometric test).

These co-variation modules highlight novel putative regulators, many not detected or prioritized by previous population-level approaches (Ciofani et al., 2012; Yosef et al., 2013). We prioritized follow-up candidates with a computational ranking scheme (Experimental Procedures). While the genes from our co-variation matrix (rows, Figure 4B) tend to be highly ranked compared to all genes also in bulk-population data (p<10−10, Wilcoxon Ranksum test) or rankings (Ciofani et al., 2012) (Table S7, Supplemental Experimental Procedures), they do not necessarily stand out in bulk population rankings (Figure S6), highlighting the distinct signal from single-cell profiles. Based on our ranking and knockout mice availability, we chose four genes, novel to Th17 function, for functional follow up: Gpr65, Toso, Plzp and Cd5l, (the latter presented in Wang et al.).

GPR65 promotes Th17-cell pathogenicity and is essential for EAE

GPR65, a glycosphingolipid receptor, is a member of the module co-varying with pro-inflammatory genes (Figure 4B), and is highly expressed in our Th1-like effector/memory cells (Figure 2D). Genetic variants in the GPR65 locus are associated with multiple sclerosis (International Multiple Sclerosis Genetics et al., 2011), ankylosing spondylitis (International Genetics of Ankylosing Spondylitis et al., 2013), inflammatory bowel disease (Jostins et al., 2012), and Crohn’s disease (Franke et al., 2010).

Naïve _Gpr65_−/− T cells differentiated with TGF-β1+IL-6 or IL-1β+IL-6+IL-23 for 96 hours show a ~40% reduction of IL-17A+ cells compared to WT controls (Figure 5A), as measured by intracellular cytokine staining (ICC, Supplemental Experimental Procedures). Memory cells from _Gpr65_−/− mice also showed a ~45% reduction in IL-17A+ cells after reactivation with IL-23 (Figure S7A). There was a reduced secretion of IL-17A (p<0.01) and IL-17F (p<10−4) (Figure 5B) and increased IL-10 secretion (p<0.01, Figure S7C) under pathogenic differentiation conditions in the knockout cells (Supplemental Experimental Procedures).

An external file that holds a picture, illustration, etc. Object name is nihms736435f5.jpg

GPR65, TOSO and PLZP are validated as T-cell pathogenicity regulators

(A,B) Reduction in IL17A-producing cells in GPR65−/− T-cells differentiated in vitro. (A) Intracellular cytokine staining for IFN-γ and IL-17A of CD4+ WT or GPR65−/− cells differentiated for 96h. (B) Quantification of secreted IL-17A and Il-17F by cytometric bead assays (CBA) in corresponding samples. * p < 0.05, ** p < 0.01, *** p < 0.001. (C) Reduced IL-17A and IFN-γ production by GPR65−/− memory (CD62L−CD44+CD4+) T cells in a recall assay (Supplemental Experimental Procedures). (D) Loss of GPR65 reduces tissue inflammation and autoimmune disease in vivo. RAG-1−/− mice (n = 10 per category) were reconstituted with 2×106 naïve WT or GPR65−/− CD4+ T-cells, and induced with EAE one week post transfer. Error bars: standard deviation. (E) Transcriptional impact of a loss of GPR65, TOSO and PLZP. Shown is the significance of enrichment (−log10 (P-value); hypergeometric test, Y axis) of genes that are dysregulated compared to WT during the TGF-β1+IL-6 differentiation of GPR65−/− (96h), PLZP−/− (48h) and TOSO−/− (96h) cells. (F,G) Reduction in IL17A-producing cells in TOSO−/− T cells differentiated in vitro. (F) Intracellular cytokine staining as in (A) but for WT or TOSO−/− CD4+ T-cells, activated in vitro for 96h. (G) Quantification of secreted IL-17A and IL-17F for WT or TOSO−/− CD4+ T cells, as in (B). (H) Reduced IL-17A production by TOSO−/− LN memory T cells in a recall assay as in (C). (I) Hampered IL-17A production by PLZP−/− CD4+ T cells in an in vitro recall assay (Supplemental Experimental Procedures). Intracellular cytokine staining for IFN-γ (Y axis) and IL-17A (X axis). (J) Quantification of secreted IL-17A and IL-17F of a MOG35–55 recall assay for littermate controls and PLZP−/− mice at 96h post ex vivo. All experiments are a representative of at least three independent experiments with at least three experimental replicates per group.

See also Figure S7, Table S6 related to Figure 5.

RNA-seq profiles (Supplemental Experimental Procedures) of populations of _Gpr65_−/− Th17 cells, differentiated under both non-pathogenic and pathogenic conditions for 96 hours, show that genes up-regulated in _Gpr65_−/− cells are most strongly enriched (p<10−28, hypergeometric test, Figure 5E) for the genes characterizing the more regulatory Th17 cells (positive PC1, Figure 4C, Table S6, Supplemental Experimental Procedures). Moreover, ChIP-Seq analysis (Ciofani et al., 2012) indicates that Rorγt, the Th17 cell master transcription factor, binds the promoter region of Gpr65.

To determine the effect of loss of GPR65 on autoimmune disease, we reconstituted RAG-1−/− mice with naïve WT or _Gpr65_−/− CD4+ T cells, and induced EAE (Supplemental Experimental Procedures). In the absence of GPR65-expressing T cells, mice are protected from EAE (Figure 5D, Figure S7D) and far fewer IL-17A and IFN-γ positive cells are recovered from the LN and spleen compared to WT controls (Figure S7B). Furthermore, in vitro restimulation with MOG35–55 of the spleen and LN cells from immunized mice showed that loss of GPR65 resulted in dramatic reduction of MOG35–55-specific IL-17A or IFN-γ positive cells (Figure 5C), suggesting that GPR65 regulates encephalitogenic T cells generation in vivo.

TOSO is implicated in Th17 pathogenicity

TOSO (FAIM3) is an immune-cell specific surface molecule that negatively regulates _Fas_-mediated apoptosis (Hitoshi et al., 1998) and a member of the module co-varying with regulatory genes (Figure 4B). While this may naively suggest that TOSO would enhance regulatory mechanisms, TOSO-knockout mice are resistant to EAE (Lang et al., 2013). Toso could therefore be a negative regulator of the non-pathogenic state.

Supporting this hypothesis, _Toso_−/− cells showed a defect in the production of the pro-inflammatory cytokine IL-17A for both differentiation conditions (Figure 5F,G) and memory _Toso_−/− cells stimulated with IL-23 lacked IL-17A production (Figure S7E). In a MOG35–55 recall assay, CD3+CD4+ _Toso_−/− T cells showed no IL-17A production (Figure 5H). This supports a role for TOSO as a promoter of pathogenicity.

Population RNA-seq analysis shows that loss of TOSO results in suppression of the key regulatory genes (e.g., Il24, Il9, Procr; Table S6), consistent with an IL-10 reduction measured by ELISA (Figure S7G), and a reduced FOXP3+ cell count during Treg differentiation (TGF-β1, Figure S7F). On the other hand, in the pathogenic condition, Il17a is down regulated in the absence of TOSO. Enrichment analysis with respect to PC1 of the non-pathogenic condition suggests that TOSO knockout cells, rather than up-regulating regulatory genes, down-regulate genes associated with a more pro-inflammatory phenotype (Figure 5E). Toso is also bound by Rorγt (Ciofani et al., 2012), providing an additional Th17 specific mechanism of action.

MOG35–55 -stimulated _Plzp_−/− cells have a defect in generating pathogenic Th17 cells

The transcription factor PLZP (ROG, ZBTB32), is a known repressor of the Th2 master regulator GATA3 and regulates cytokine expression (Miaw et al., 2000) in T-helper cells. We hypothesized that Plzp regulates pathogenicity in Th17 cells, but we could not undertake an EAE experiment since PLZP−/− mice were not available on an EAE-susceptible background.

While in vitro differentiated _Plzp_−/− cells produced similar IL-17A levels as WT controls (Figure S7H), a MOG35–55 recall assay revealed a defect in IL-17A production with increasing MOG35–55 concentration during restimulation (Figure 5I). When reactivated in the presence of IL-23, which expands in vivo generated Th17 cells, _Plzp_−/− cells also produced less IL-17A (Figure S7I). Plzp appears to influence the expression of a wider range of inflammatory cytokines, as _Plzp_−/− T cells secreted less IL-17A, IL-17F (Figure 5J), IFN-γ, IL-13 and CSF2 (Figure S7J).

Based on RNA-seq profiles at 48 hours of non-pathogenic differentiation of _Plzp_−/− cells, Irf1, Il9 and other transcripts of the module co-varying with regulatory genes are up regulated (Table S6), whereas transcripts from the module co-varying with pro-inflammatory genes (e.g., Ccl20, Tnf, Il17a) are repressed, and genes characterizing the more pro-inflammatory cells (PC1, Figure 4C) are strongly enriched among the down-regulated genes (Figure 5E).

DISCUSSION

Here, we show how variation and co-variation in single cell profiles can be leveraged to identify key regulatory modules and the factors that may control them, to dissect Th17 cell pathogenicity, beyond differentiation.

In vivo, we used variation to infer the life cycle of Th17 cells. Processes such as self-renewal, observed in the LN, may provide a pool of cells that are precursors for differentiating Th17 cells to effector/ memory formation in the CNS. The Th1-like phenotype we observe in the CNS may be the most pathogenic (Bending et al., 2009; Lee et al., 2009; Muranski et al., 2011) and might facilitate memory cell formation, as the entry of Th1 cells into the memory pool is well established (Sallusto et al., 1999). It is unclear if cells that adopt a Th1 phenotype are stable ‘double producers’ or if they show plasticity towards a Th1 fate.

We used transcription factor target enrichment analysis in vivo to nominate key regulators of each state. For example, we predict that Med12, Etv6 and Zfx drive the Th17 self-renewing-like subpopulation in the LN. While neither has been linked to Th17 self-renewal, each has been associated with self-renewal and related functions in other cells (Rocha et al., 2010) (Hock et al., 2004) (Tsuzuki and Seto, 2013) (Galan-Caridad et al., 2007). For the pathogenic effector and memory cells observed in the CNS during EAE, we assign a prominent role to known Th17/Th1 transcription factors such as Hif1a, Fosl2, Stat4 and Rel.

In vitro, we used strong co-variation, most pronounced under the least pathogenic and most variable conditions to rank candidate genes, such as Cd5l and Gpr65, based on their association with known regulatory and pro-inflammatory genes. Consistently, a lack of both Cd5l and Gpr65 significantly alters EAE disease progression. Genes similarly associated with pro-inflammatory functions, which we have not yet followed up on include: Gem, Tmem109, and cd226. Conversely, Foxp1, a member of the module co-varying with regulatory genes, was highly expressed in the LN-derived Th17 self-renewing subpopulation and the gut-derived Th17 cells (Figure S3). Foxp1 negatively regulates IL-21, a driver of Th17 generation (Korn et al., 2007), and dampens T-cell activation (Wang et al., 2014). Co-variation of a gene with a particular module does not, however, necessarily indicate similar function of this gene with other genes in the module, as we have seen for TOSO. Another example, Lag-3, is up-regulated during T cell activation, but suppresses it (Grosso et al., 2007). This is consistent with a model where regulators with opposite, antagonistic functions are co-regulated.

Whereas population-based expression profiling has identified genes that govern the differentiation states of Th17 cells, single-cell RNA-seq provides new granularity to unveil potent candidates for manipulation of pathogenicity of Th17 cells without affecting nonpathogenic Th17 cells that may be critical for tissue homeostasis and for maintaining barrier functions.

EXPERIMENTAL PROCEDURES

Mice, EAE induction and cell isolation

C57BL/6 WT and CD4−/−(2663) mice were obtained from Jackson Laboratory. IL-17A–GFP+ mice were obtained from Biocytogen. GPR65−/−, PLZP−/− and TOSO−/− mice were provided by Li Yang, Pier Paolo Pandolfi and John Coligan, respectively. All animals, unless noted otherwise, were housed and maintained in a conventional pathogen-free facility at the Harvard Institute of Medicine in Boston (Supplemental Experimental Procedures). EAE induction and disease analysis, isolation of T cells from EAE mice at the peak of disease, isolation of T cells and in vitro differentiation, isolation of memory cells and recall assays, and isolation of T cells from lamina propria was performed as described in the Supplemental Experimental Procedures.

RNA-Seq

Whole transcriptome amplification of cell lysates was performed by SMART-Seq (Ramskold et al., 2012) using the Fluidigm C1 Single-Cell Auto Prep System, followed by Nextera XT DNA Sample preparation (Illumina), as described in the Supplemental Experimental Procedures. We collected at least two independent biological replicates for each in vivo and in vitro condition, and two technical replicates for two in vivo conditions.

RNA-seq preprocessing

RNA-seq reads alignment, transcript quantification were performed as described in the Supplemental Experimental Procedures. We used log transform and quantile normalization to further normalize the expression values (FPKM) within each batch of samples (i.e., all single cells in a given run), and accounted for low (or zero) expression by adding a value of 1 prior to log transform. For each library we computed quality scores using Fastqc, Picard tools, and in-house scripts, excluding poor libraries from further analysis and further adjusting for the quality scores (Supplemental Experimental Procedures).

Batch correction

We performed batch correction separately for in vivo and in vitro samples. A filtered gene set consists of the genes that have an expression level exceeding 10 FPKM in at least 20% of the cells of a given sample (Supplemental Experimental Procedures).

Taking into account false negatives using a weighted analysis

To account for the effect of each gene’s expression and each cell’s quality on the probability of false negatives with zero transcript abundance we construct for each cell a false-negative curve (FNC) representing the false-negative rate as a function of transcript abundance in the bulk population, and use this to weight subsequent analyses (Supplemental Experimental Procedures).

Signature scores and gene set enrichment analysis (GSEA)

To interpret the functional implications of the variation between cells, we assembled a set of gene signatures that are indicative of various cell states. A typical signature is comprised of a “plus” subset and a “minus” subset. A strong match will have extreme, and opposite values for the expression of genes in the two sets (Supplemental Experimental Procedures).

Gene ranking

We rank genes in the co-variation modules that significantly correlate (Spearman correlation with FDR<0.05, using the Benjamini-Hochberg scheme) with at least one of the genes in the curated set of bimodally expressed immune response genes (columns of Figure 4B) by five criteria (Supplemental Experimental Procedures): (1) correlation with the first PC in the in vitro derived Th17 cells (using TGF-β1+IL6) (2) correlation with the first and (3) second PCs in the in vivo derived Th17 cells; (4) correlation with immune-related genes that are specified in the columns of Figure 4B. (5) a similar analysis using a curated pathogenicity signature (genes that are positively or negatively associated with pathogenic Th17 cells based on population-level experiments (Lee et al., 2012)).

Additional analyses and details are in the Supplemental Experimental Procedures.

Supplementary Material

1

3

4

5

6

7

8

9

10

11

12

13

14

15

16

2

Acknowledgements

We thank J. E. Coligan for providing TOSO−/−mice; J. Shuga for scientific discussions; D. Kozoriz, T. Rogers and M. Tam for assistance with cell sorting; T.J. Diefenbach for assistance with the Imagestream system; C. Wu for help with EAE experiments; O. Rozenblatt-Rosen, E. Shefler, C. Guiducci for project management; the Broad Genomics Platform for all sequencing work; and L. Gaffney for help with artwork. Work was supported by an NIH grant P50 HG006193 (A.R. and H.P.), by the Koch Institute Support (core) Grant P30-CA14051 from the National Cancer Institute (A.R.) and the Klarman Cell Observatory at the Broad Institute (A.R., H.P. and V.K.K.). V.K.K. was supported by NIH grants NS030843, AI039671, NS076410, AI056299, National MS Society, New York and Crohn’s and Colitis Foundation of America. C.W. was supported by a Career Transition Award from the National Multiple Sclerosis Society (TA3059-A-2), and NIH grant 4R00AI110649. R.S. was supported by NIH NSRA grant F32 HD075541. A.R. is an Investigator of the Howard Hughes Medical Institute.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Contributions

J.T.G., N.Y., Y.L., A.K.S., V.K.K., H.P. and A.R. conceived the study and designed experiments. A.R. and N.Y. devised analyses and N.Y. developed computational methods. N.Y., J.T.G., Y.L. and R.S analyzed the data. J.T.G., Y.L. and R.S.G. conducted the experiments. L.V.Y., P.P.P. and T.M. provided knockout mice. J.T.G, N.Y., V.K.K., H.P. and A.R. wrote the paper with input from all the authors.

Accession Numbers

All RNA-seq data is submitted to GEO and accession numbers will be provided during proofs.

References