Regulatory T (T reg) cells, a specialized subset of CD4+ T cells, are essential to prevent fatal autoimmunity. Expression of the T reg lineage-defining transcription factor Foxp3, and therefore their differentiation in the thymus, is dependent upon T cell receptor (TCR) and interleukin-2 (IL-2) signaling. Here, we report that the majority of IL-2–producing cells in the thymus are mature CD4 single-positive (CD4SP) thymocytes and that continuous IL-2 production sustained thymic T reg cell generation and control of systemic immune activation. Furthermore, single-cell RNA sequencing analysis of CD4 thymocyte subsets revealed that IL-2 was expressed in self-reactive CD4SP thymocytes, which also contain T reg precursor cells. Thus, our results suggest that the thymic T reg cell pool size is scaled by a key niche factor, IL-2, produced by self-reactive CD4SP thymocytes. This IL-2–dependent scaling of thymic T reg cell generation by overall self-reactivity of a mature post-selection thymic precursor pool may likely ensure adequate control of autoimmunity.

Regulatory T (T reg) cells, a specialized subset of CD4+ T cells defined by expression of Foxp3, are critical for protection against life-threatening autoimmunity (Sakaguchi et al., 1995; Kim et al., 2007; Lahl et al., 2007). Neonatal thymectomy studies in mice demonstrated that the thymus serves as the main source of suppressive T reg cells (Asano et al., 1996). T reg cell differentiation in the thymus proceeds in a two-step fashion. Upon strong TCR stimulation by self-peptide presented on MHC class II complexes, CD4 single-positive (CD4SP) thymocytes undergo changes in gene expression, among which up-regulation of CD25, the high-affinity subunit of the IL-2 receptor, is particularly important. Induction of CD25 expression endows immediate Foxp3 T reg cell precursors with the ability to effectively compete for the limiting amount of IL-2 in the thymic microenvironment (Burchill et al., 2008; Lio and Hsieh, 2008). The resulting signaling through the IL-2 receptor leading to STAT5 phosphorylation in CD25+Foxp3CD4+ thymocytes facilitates induction and maintenance of Foxp3, the X chromosome–linked lineage-defining transcription factor, whose continuous expression is required for T reg cell differentiation and function (Burchill et al., 2007).

The size of the T reg cell population needs to be tightly controlled, to ensure suppression of autoimmunity and inflammation mediated by self-reactive T cells, while allowing for protection against microbial and abiotic challenges (Lee et al., 2012). One factor contributing to T reg cell differentiation in the thymus is the random generation of self-reactive TCRs. However, according to the two-step model, cell-intrinsic TCR signaling in T reg precursors confers potential to give rise to T reg cells, but by itself is insufficient for their differentiation. Thus, the mechanisms of scaling of the T reg cell population to the overall pool of self-reactive T cells in the thymus remains unknown.

On a general level, the size of differentiating cell populations is defined by “niche” cells, which serve as a source of a limiting key factor(s) and are frequently of a developmental origin distinct from their “client” cells. The former, by acting in a cell-extrinsic manner, affect the pool size of differentiating cells or their precursors. Considering the rather unique function of T reg cells as a dominant negative feedback to control self-reactivity, we hypothesized that contrary to a customary distinct ontogeny of niche cells, a cell-extrinsic mechanism controlling T reg population size is the limiting IL-2 produced by self-reactive thymocytes themselves. This ensures a newly generated T reg cell population size proportional to the size of the developing self-reactive T cell pool. Here we report that a population of self-reactive CD4SP thymocytes (with a signature of recent TCR engagement) are the major source of IL-2 in the thymus that is required for efficient T reg cell differentiation and maintenance of immune homeostasis.

Continuous dependence of thymic T reg cell generation on IL-2

Germline deficiency in IL-2 or a T reg cell–specific deletion of the IL-2 receptor subunits results in an early onset fatal autoimmunity and wasting disease, characterized by a severe reduction of T reg cells in the thymus and periphery (D’Cruz and Klein, 2005; Fontenot et al., 2005a; Chinen et al., 2016). Even though T reg cells express the high affinity IL-2 receptor (IL2Rα/CD25), they themselves do not produce IL-2 and are dependent on paracrine sources of IL-2 (Thornton and Shevach, 1998; Pan et al., 2009). Since early induction of systemic autoimmune disease including severe anemia in mice congenitally lacking IL-2 can obfuscate assessment of a role of IL-2 in T reg cell generation and does not allow its assessment in adulthood, we induced ablation of a conditional Il2 allele (Il2fl) in adult mice and observed a profound decrease in the frequency of T reg cells in the thymus and effector CD4 T cell expansion and activation in the periphery (Fig. 1 A; data not shown). This observation highlighted a continuous dependence of thymic T reg cell generation on IL-2, consistent with its proposed role as a niche factor.

IL-2–producing cells in the thymus

The cellular sources of IL-2 that promote T reg cell differentiation in the thymus and maintenance in secondary lymphoid organs or in nonlymphoid tissues have been the subject of considerable debate, and specific features of these cells remained unknown. Among candidate thymic IL-2 producers have been thymocytes, thymic dendritic cells (DCs), and stromal cells (Weist et al., 2015; Owen et al., 2018). To identify the specific cellular sources of IL-2, we generated Il2cre/+Rosa26lox-STOP-lox-Tom mice, in which cells with a developmental or "life-history" of IL-2 expression were permanently marked by a fluorescent reporter protein (tdTomato; tdTom), whose expression is induced by Cre recombinase driven by the endogenous Il2 locus. A small percentage of cells in the thymus were tagged with tdTom (Il2cre–fate-mapped [FM] cells), and these cells were mostly TCRβ-expressing CD4+ thymocytes (Fig. 1, B and C; and Fig. S1 A and B). Labeled (“IL-2 fate-mapped”) CD4SP thymocytes could already be detected in neonatal mice as early as postnatal day 3 (P3) and continued to increase into adulthood (Fig. 1 D, left panel). Of note, in mice at P3, the thymus was the only tissue in which labeled cells were detectable, and P3 coincided with the first T reg cells emerging from the thymus, supporting the idea that thymocyte production of IL-2 is a limiting factor for T reg cell generation. We also consistently detected FM T reg cells in the thymus, spleen, and lung of neonatal and adult mice, suggesting a life-history of IL-2 expression in precursors of some T reg cells, even though their proportion was markedly lower than that of FM non–T reg cells (Fig. 1 D, right panel). The frequencies of labeled cells in the thymus were lower than in spleen or lung, which might be a reflection of thymic emigration during Cre-mediated recombination and expression of the fluorescent reporter or a preferential expansion of those cells in the periphery. In contrast, we failed to detect tdTom-expressing DCs or CD45 stromal cells (Fig. S1 C). Consistent with these observations, ablation of a conditional Il2 allele in DCs using a highly specific Cre driver (Il2Zbtb46cre/+) did not affect T reg cell development in the thymus, and the peripheral maintenance of immune homeostasis was unimpeded, ruling out thymic DCs as a source of IL-2 biologically meaningful for T reg cell generation (Fig. S1, D–F).

To assess the overall potential for stimulation-induced IL-2 protein production by different thymocyte subsets, we assayed bulk thymocytes ex vivo and detected IL-2 production by intracellular cytokine staining in CD4SP thymocytes, but to a much lesser degree in CD8SP, TCRγδ+ T cells, and natural killer T cells (Fig. 1 E). Since CD4SP thymocytes represented the major potential source of IL-2 in the thymus, we sought to more precisely define the IL-2 producers among them, specifically, by dividing them into four major subsets: thymic T reg cells (CD4+Foxp3+), CD4+CD25+Foxp3 cells, and CD4+CD25Foxp3 cells that were further subdivided into immature and mature CD4SP cells based on CD24 and CD62L cell surface expression (for gating strategy, see Fig. S2). These four subpopulations were assessed for chromatin accessibility at the Il2 locus measured by ATAC-seq (Assay for Transposase-Accessible Chromatin using sequencing) analysis. Strikingly, only the CD4+CD25+Foxp3 subset exhibited an active Il2 locus marked by open chromatin, suggesting that this population is enriched in cells that express IL-2, have the potential to express IL-2, or have a recent history of its expression (Fig. 1 F). We surmised that these cells may have recently received strong TCR stimulation upon self-peptide:MHC class II recognition that led to IL-2 expression and CD25 up-regulation. Indeed, in line with previously published data, we found that CD4+CD25+Foxp3 cells expressed high levels of Nur77 when compared with the mature CD4SP subset, supporting the idea that this population is enriched in self-reactive T cells (Fig. 1 G; Moran et al., 2011). To further delineate the ability of these various subsets to secrete IL-2, we sorted the populations described above and stimulated them ex vivo. Immature CD4SP and thymic T reg cells expressed and secreted little IL-2. The majority of IL-2 was produced by mature CD4SP and to a lesser degree by CD4+CD25+ thymocytes (Fig. 1, H and I). These data represent the potential of all cells within each subpopulation to produce IL-2 after encountering strong TCR stimulation through PMA and ionomycin (Fig. 1 H) or TCR cross-linking (Fig. 1 I), irrespective of their TCR affinity for self-antigens, the predominant source of TCR ligands encountered by thymocytes in an unperturbed animal.

Thymocytes, but not peripheral T cells, provide IL-2 for T reg cell generation

To confirm the role of thymocyte-derived IL-2 in thymic T reg cell generation, we employed various genetic approaches. Deleting an Il2fl allele in all T cells early during thymocyte differentiation using proximal Lck promoter-Cre (Il2pLckcre) or CD4-Cre transgenes (data not shown) revealed a profound defect in T reg cell differentiation in the thymus (Fig. 2 A). The frequency of CD25+Foxp3 CD4SP was unchanged, but the CD25+Foxp3+ T reg cell compartment was reduced by ∼50%. Of note, the overall frequency of Foxp3+ thymocytes was only slightly decreased, due to an increase in CD25Foxp3+ thymocytes, a recently described putative T reg cell differentiation intermediate (Owen et al., 2019). This reduction of T reg cells was also evident in secondary lymphoid organs and was accompanied by increases in activated and proliferating CD4 and CD8 T cells and elevated IFN-γ production observable already in young mice (Fig. 2, B–F). These findings raised the question whether a small number of IL-2–producing activated peripheral CD4+CD25+Foxp3 T cells re-enter the thymus and serve as a source of IL-2 for T reg cell differentiation. To exclude the latter possibility, we ablated IL-2 very late in thymocyte differentiation and in peripheral T cells using a distal Lck (dLck) promoter–driven Cre transgene (Il2dLckcre). We found no discernible impact on T reg cell differentiation in the thymus (Fig. 2 G); yet peripheral T cell activation was observed, albeit less pronounced than in Il2pLckcre mice (Fig. S3). In addition, our IL-2 fate mapping in neonatal mice described above also strongly supports a thymic origin of the IL-2 producers at the time of T reg cell emergence (Fig. 1 D). These data combined with the ATAC-seq results suggest that self-reactive CD25+Foxp3CD4SP thymocytes, a subset that also contains T reg precursor cells, are the relevant source of IL-2 that promote T reg cell generation in the thymus.

TCR-induced IL-2 production by CD4SP thymocytes promotes T reg cell generation

To explore the relationship between producers and consumers of IL-2, we set up an in vitro model of thymic T reg cell differentiation based on a two-step process requiring strong TCR stimulation followed by exposure to IL-2, leading to induction of Foxp3 expression. Both immature and mature CD4SP cells, upon stepwise TCR stimulation and exposure to IL-2, could up-regulate Foxp3 protein expression. The efficiency of Foxp3 induction scaled with the amount of IL-2 provided, and was higher with mature CD4SP as a starting population (Fig. 3 A). TCR stimulation was absolutely required for Foxp3 induction, as IL-2 alone was not sufficient to induce Foxp3 in either subset (Fig. 3 B). If starting with CD4+CD25+ CD4SP cells, which have presumably experienced strong TCR signaling, exposure to IL-2 alone was sufficient to induce Foxp3 expression, in agreement with a previous report (Fig. 3 C; Lio and Hsieh, 2008). Exposure to conditioned media after withdrawal of TCR stimulation was sufficient to induce Foxp3 expression in mature CD4SP, and the inducing activity was reduced by ≥50% upon antibody-mediated blockade of IL-2 (Fig. 3 B). To isolate the Foxp3 expression–promoting effect of IL-2 signaling from TCR stimulation and consequent IL-2 production, we next employed stimulation with the bacterial super-antigen staphylococcal enterotoxin B (SEB), which strongly stimulates T cells expressing TCR Vβ8 chains (White et al., 1989). Expectedly, SEB preferentially induced Foxp3 expression in Vβ8-bearing T cells in a dose-dependent manner, which was strongly impaired in the presence of blocking antibodies against IL-2 (Fig. 3 D, gray bars). Up-regulation of CD25, a surrogate for TCR stimulation, was also preferentially observed in TCR Vβ8-expressing cells, but was largely unaffected when IL-2 was neutralized (data not shown). Interestingly, we consistently observed differentiation of Vβ8-negative self-reactive T cells, even if markedly less prominent (Fig. 3 D, black bars). These results suggest that IL-2 produced by Vβ8-expressing thymocytes in response to SEB stimulated Foxp3 expression in “bystander” T reg cell precursors. Similarly, treating mice with SEB led to an overall expansion of T reg cells and IL-2 FM cells in the thymus at 48 h after injection. While the majority of the increase was contributed by thymocytes receiving strong TCR stimulation through Vβ8+ TCR chains, we observed consistent increases in Vβ5- and Vβ6-bearing T reg cells, albeit to a smaller degree (Fig. 3, E–G). Importantly, the expansion of T reg cells was blunted in Il2pLckcre mice injected with SEB (Fig. 3, H–J). These results suggest that T reg precursor cells, upon receiving a strong TCR signal, can produce IL-2, which in turn promotes their own differentiation into T reg cells, or differentiation of thymocytes, which have received a TCR signal before paracrine IL-2 signaling.

Single cell RNA-sequencing (scRNA-seq) analysis of IL-2 production by CD4SP thymocytes

To gain insight into characteristic features of thymocytes producing IL-2 in situ, we performed scRNA-seq analysis of mature CD25Foxp3CD4SP, CD25+Foxp3, and Foxp3+ CD4SP thymocytes. After initial processing, a total of 11,874 individual cellular transcriptomes were analyzed using unsupervised graphical clustering and visualized using t-distributed stochastic neighbor embedding (t-SNE). The three analyzed subsets had largely distinct distribution patterns in the resulting t-SNE plots (Fig. 4 A). The cells yielded 12 distinct clusters (clusters 0–11; Fig. 4, C and D). Performing unsupervised clustering of the top 30 differentially expressed genes from each cluster revealed 7 distinct groups of similarly regulated genes (Fig. 4 E). Clusters 1, 5, and 6 contained predominantly thymic T reg cells, enriched for expression of Foxp3, Il2rb, Izumo1r, and Ctla4 strongly associated with a T reg cell signature (Zemmour et al., 2018; Fig. 4, D and E, group V). Cluster 6 in particular seemed to be enriched for genes that are associated with T reg cell activation (Ccl5, Maf, Icos, and Nrp1; Fig. 4, B and E). Cluster 11 cells were defined by a strong IFN signaling signature (Stat1, Ifit1, and Irf7), potentially encompassing cells undergoing late stages of T cell maturation (Xing et al., 2016; Fig. 4 E, group VI). Cluster 9 cells were enriched for gene expression associated with recent TCR signaling (Nr4a1, CD69, Rel, and Orai1). These same genes were also enriched, albeit to a somewhat lesser degree, in clusters 4, 7, and 8, which contained CD25+Foxp3 and Foxp3+ cells, probably a reflection of recent TCR engagement (Fig. 4, D and E, group II).

To analyze the most significant sources of phenotypic variation in our dataset, we performed diffusion map analysis (Fig. 4 F). The first three components correlated with signatures for IL2-STAT5 signaling/T reg cell activation, apoptosis, and TCR signaling (Fig. 4, F and G). Diffusion components 1 and 3 recapitulated T reg cell differentiation in the thymus, including the need for TCR signaling to progress from the mature CD4SP stage to facilitate up-regulation of CD25 (Il2ra) to compete for IL-2 and subsequent STAT5-dependent induction of Foxp3 (Fig. 4, B, F, and G). The second diffusion component highlighted apoptosis, specifically in the autoreactive CD25+Foxp3 thymocyte population, as reflected by the heightened expression of the pro-apoptotic Bcl2l11 gene, whereas most other cells expressed increased levels of the anti-apoptotic Bcl2 gene (Fig. 4, B, F, and G). These data suggest at least three cell fates for mature CD4SP thymocytes: death by apoptosis largely restricted to a proportion of the CD4+CD25+ thymocytes or S1PR1-mediated emigration as a CD4+Foxp3 or a CD4+Foxp3+ T cell (Fig. 4 E, group VII; Matloubian et al., 2004). Within our dataset, we were able to identify 428 cells with detectable Il2 transcript (Fig. 4 B). Their individual transcriptomes, visualized by t-SNE projection (Fig. 5 A), were enriched for genes associated with activation of the NF-κB pathway (Nfkb1, Nfkbia) likely downstream of recent TCR engagement evidenced by increased CD69 expression. This TCR signaling signature is in strong support of our assumption that the IL-2 producers are enriched in the self-reactive pool of CD4SP thymocytes, since self-antigens are predominantly presented by MHC class II in the thymus of an unperturbed mouse. Consistent with a known role for NF-κB signaling in promoting survival, we observed increased Bcl2 expression, which offers protection from TNF-mediated death, and the type I IFN pathway (Isg15, Stat1) also required for late-stage maturation of T cells in the thymus (Fig. 5, B–D; Xing et al., 2016). About a third of the cells co-expressed Foxp3 and Il2, consistent with the idea that at least a fraction of IL-2–producing self-reactive precursor cells can differentiate into thymic T reg cells and in agreement with the presence of Il2cre-FM cells within the Foxp3+ thymocyte and peripheral population (Fig. 5 Band1 D). Of note, a small proportion of IL-2–expressing cells also expressed pro-apoptotic Bim (Bcl2l11) and the co-stimulatory receptor 4-1BB (Tnfrsf9), and those cells were mostly found in the CD25+Foxp3 thymocyte population (Fig. 5 B). This observation suggests that some self-reactive thymocytes, upon TCR stimulation, produce IL-2 but then undergo apoptosis (Fig. 5, C and D). Interestingly, 4-1BB signaling is thought to promote IL-2 production in peripheral conventional T cells by overriding T reg cell–mediated suppression, a mechanism potentially also at work in the thymus (Barsoumian et al., 2016).

In aggregate, our studies suggest that in the thymus, IL-2 was expressed predominantly by self-reactive CD4SP thymocytes, which also harbor CD25+Foxp3 T reg precursor cells. IL-2–expressing thymocytes exhibited a gene expression pattern characteristic of TCR engagement, NF-κB signaling, and tonic type I IFN signaling and served as a key source of IL-2 for developing thymic T reg cells (Fig. 5 E). Thus, self-reactive thymocytes regulate the size of the T reg cell population generated in the thymus through production of a niche factor: IL-2.

Mice

Il2fl, Il2cre/+, Foxp3GFP, and Foxp3Thy1.1 mice have been described previously (Fontenot et al., 2005b; Liston et al., 2008; Popmihajlov et al., 2012; Yamamoto et al., 2013). Il2fl mice were provided by Vandana Kalia (University of Washington School of Medicine, Seattle, WA), Foxp3GFP and Foxp3Thy1.1 were maintained in house, and the rest of the mice were obtained from The Jackson Laboratory. Nur77RFP mice were generated using the bacterial artificial chromosome (BAC) clone RP24-366J14, carrying a GFP cassette at the 5′ end of the Nr4a1 (Nur77) gene, a gift from Kristin Hogquist (University of Minnesota, Minneapolis, MN; Moran et al., 2011). The GFP was replaced with the turboRFP coding sequence using standard recombineering protocols, and the Nur77RFP BAC clone was validated by PCR, Sanger sequencing, and SpeI digestion. Purified BAC DNA was then injected into C57BL/6J embryos, and offspring were screened for the presence of the transgene by PCR and bred to Foxp3GFP reporter mice for use in this study. Il2fl mice were used to generate Il2pLckcre mice by breeding with B6.Cg-Tg(Lck-cre)548Jxm/J, Il2dLckcre mice by breeding with B6.Cg-Tg(Lck-icre)3779Nik/J, Il2Zbtb46cre/+ mice by breeding with B6.Cg-Zbtb46tm3.1(cre)Mnz/J, and Il2UbcreERt2/+ mice by breeding with B6.Cg-Ndor1Tg(UBC-cre/ERT2)1EjB/J. To induce deletion in the latter strain, mice were placed on a tamoxifen citrate–containing diet (TD.130860; Envigo) at 6–8 wk of age for 4 wk. For genetic fate mapping of IL-2–expressing cells, Il2cre/+ mice were crossed with ROSA26lsltom/+ mice (B6.Cg-Gt(ROSA)26Sortm14(CAG-tdTom)Hze/J) and either Foxp3GFP or Foxp3Thy1.1 reporter mice. All animals were maintained in the Memorial Sloan Kettering Cancer Center animal facility under specific pathogen–free conditions, and the experiments were performed according to the institutional guidelines (Institutional Animal Care and Use Committee 08–10-23). Mice aged 5–8 wk were used for all experiments unless otherwise specified in the figure legend. Mice were age matched and litter matched. Mice were sex matched when possible.

Tissue preparation and cell isolation

For analysis of single-cell suspensions, thymus, spleen, LNs, and Peyer’s patches were mechanically dissociated between frosted glass slides or with the back of a syringe plunger and filtered through a 100-µm nylon mesh. For isolation and analysis of DCs, spleen and thymus were enzymatically digested with 300 µg/ml Liberase TL (Roche) and 50 µg/ml DNaseI (Sigma-Aldrich) in RPMI-1640 for 20 min at 37°C while shaking. Splenic DCs were further enriched using magnetic CD11c MicroBeads UltraPure and LS columns (Miltenyi Biotec). For analysis of immune cells from liver, tissues were isolated after phosphate-buffered saline (DPBS) perfusion, and liver tissue was enzymatically digested with 1 mg/ml collagenase A from Clostridium histolyticum (Roche) and 1 U/ml DNaseI in RPMI-1640 with 2% FBS, 100 µg/ml penicillin/streptomycin, 20 mM L-glutamine, and 10 mM Hepes for 30–45 min at 37°C while shaking in the presence of three 0.25-inch ceramic spheres (MP Biomedicals). Digested tissue was filtered through 100-µm cell strainers and further enriched by 40% Percoll (GE Healthcare) centrifugation and resuspended in RPMI-10 (RPMI-1640, 10% FBS, 20 mM L-glutamine, 100 µg/ml penicillin/streptomycin, 10 mM Hepes, 1 mM sodium pyruvate, 1× MEM nonessential amino acids, and 10−5 M 2-mercaptoethanol).

To assess cytokine production after ex vivo re-stimulation, single-cell suspensions were incubated for 3–4 h at 37°C with 5% CO2 in the presence of 50 ng/ml PMA and 500 ng/ml ionomycin with 1 µg/ml brefeldin A (all Sigma-Aldrich) and 2 µM monensin (BioLegend).

Flow cytometry and antibodies

All flow cytometry data were collected on a BD LSRII (Becton Dickinson) flow cytometer or the Aurora spectral analyzer (Cytek) and analyzed using FlowJo 10.5 software (TreeStar Software). Single-cell suspensions were re-suspended in DPBS and stained with Ghost Dye Violet 510 or Ghost Dye Red 780 (Tonbo Biosciences) for dead-cell exclusion in the presence of anti-CD16/32 (Tonbo Biosciences) to block binding to Fc receptors for 10 min at 4°C. Cell-surface antigens were stained for 20 min at 4°C in FACS buffer (DPBS, 0.5% BSA, and 2 mM EDTA). Cells were analyzed unfixed, fixed with 2% paraformaldehyde (Thermo Fisher Scientific) for later analysis, or fixed and permeabilized with the BD Cytofix/Cytoperm kit (for intracellular cytokine staining) or the eBioscience Foxp3/Transcription Factor staining buffer set. Fixation, permeabilization, and intracellular staining were performed as recommended by the manufacturer. All cells were re-suspended in FACS buffer and filtered through 100-µm nylon mesh before analysis on the flow cytometer. 123count eBeads (eBioscience) were added at 5,000 beads per sample to quantify absolute cell numbers.

Antibodies used in this study were purchased from BioLegend, BD Biosciences, eBioscience, Thermo Fisher Scientific, or Tonbo Biosciences. The following clones were used: CD45 (30-F11), CD45.1 (A20), CD45.2 (104), NK1.1 (PK136), Ter-119 (Ter-119), TCRγδ (GL3), TCRβ (H57-597), TCR Vβ5.1, 5.2 (MR9-4), TCR Vβ6 (RR4-7), TCR Vβ8.1, 8.2 (KJ16-133.18), CD3ɛ (145-2C11), CD4 (RM4-5, GK1.5), CD8a (53–6.7), CD8b (YTS156.7.7), CD62L (MEL-14), CD44 (IM7), CD24 (M1/69), CD25 (PC61), CD90.1 (HIS51), CD90.2 (30-H12, 53–2.1), Ki67 (SolA15), B220 (RA3-6B2), CD19 (6D5, 1D3), CD11c (HL3), MHCII (I-A/I-E; M5/114.15.2), Foxp3 (FJK-16s), IFNγ (XMG1.2), and IL-2 (JES6-5H4).

Cell sorting and in vitro T reg cell induction assays

Thymi from Foxp3GFP or Foxp3Thy1.1 reporter mice were dissected and disassociated into a single-cell suspension. Cells were re-suspended in depletion buffer (DPBS, 2% FBS, and 2 mM EDTA) and enriched for CD4SP thymocytes by depleting CD8SP and double-positive cells using the Dynabeads Flowcomp mouse CD8 kit (Life Technologies). The enriched cells were stained with fluorochrome-conjugated anti-CD45, anti-TCRβ, anti-CD4, anti-CD24, anti-CD62L, and anti-CD25 and biotinylated anti-CD8, anti-TCRγδ, and anti-NK1.1 (dump gate), followed by fluorochrome-conjugated streptavidin before sorting on a BD FACSAria II sorter. The populations were sorted using gates shown in Fig. S2.

A typical in vitro T reg cell induction assay was set up with 2 × 104 to 5 × 104 cells per well of a flat-bottom 96-well tissue culture plate (in most cases immature and mature CD4SP thymocytes). The 96-well plates had been either coated with anti-CD3 and anti-CD28 at 1 µg/ml in DPBS for 12 h at 4°C or left untreated. Cells were cultured in RPMI-10 for 10–12 h in those plates and then transferred to V-bottom 96-well plates and re-suspended in fresh RPMI-10 supplemented with recombinant human IL-2 (Roche) at the indicated concentrations or re-suspended in their own media (conditioned media) for an additional 12 h. To block the effects of IL-2, a cocktail of IL-2–blocking antibodies was added (S4B6-1, JES6-1A12; 5 µg/ml; BioXcell) throughout.

For the SEB-mediated induction of Foxp3 expression, mature CD4SP thymocytes were co-cultured with CD11c bead–enriched splenic DCs that had been further FACS purified as CD11c+MHCIIhi cells at a 2:1 ratio. SEB (Millipore) was added at a range of concentrations as indicated. After 12 h of culture, TCR stimulation was blocked by addition of I-Ab specific antibody Y3P (1 µg/ml), and the cells were cultured for an additional 12 h before Foxp3 protein induction was assessed.

ELISA

Quantification of IL-2 protein levels in cell culture supernatants was performed using the Ready-SET-Go! Mouse IL-2 kit (eBioscience) according to the manufacturer’s instructions. ELISAs were read at OD450 on a Synergy HTX instrument (BioTek). IL-2 levels were calculated through extrapolation from a standard curve using Prism 7 software.

In vivo mouse treatments

For SEB treatment, mice were injected with 20 µg SEB (Millipore) or DPBS i.v. under isoflurane anesthesia and analyzed 48 h after injection.

Statistical analysis

Graph Pad Prism (v. 7.0d) was used for all statistical analysis for biological experiments. We calculated P values by unpaired Student’s t test, and differences between groups were considered statistically significant with a P value cut-off of 0.05 or less. Data are represented as means ± SEM. In all cases, ****, P ≤ 0.0001; ***, P ≤ 0.001; **, P ≤ 0.01; *, P ≤ 0.05; ns, P > 0.05; ND, not detectable/below limit of detection.

ATAC-seq and computational analysis

Thymi from Foxp3GFP mice were dissociated, enriched for CD4SP cells, and stained as described above. Cells were then sorted into DNA Lo-bind tubes as follows: Foxp3-GFP+ CD4+ CD8 TCRβ+ (Foxp3+), CD25+ Foxp3-GFP CD4+ CD8 TCRβ+ (CD25+Foxp3), CD25 Foxp3-GFP CD4+ CD8 TCRβhi CD24lo (mature CD4SP), and CD25 Foxp3-GFP CD4+ CD8 TCRβlo CD24hi (immature CD4SP). There were two biological replicates per cell population with 40–50,000 cells.

Preparation of the ATAC-seq library

Immediately after sorting, cells were processed to generate ATAC libraries as previously described with minor modifications (Buenrostro et al., 2015). Briefly, cells were pelleted in a fixed-rotor mini-centrifuge at 500 ×g for 5 min at 4°C. Cells were then washed in 1 ml cold PBS and pelleted again. Supernatant was aspirated, and cells were re-suspended in 50 μl ice-cold cell lysis buffer (10 mM Tris-HCl, pH 7.4; 10 mM NaCl; 3 mM MgCl2; and 0.1% NP-40) to disrupt plasma membranes. Nuclei were pelleted at 1,000 ×g for 10 min. Supernatant was aspirated, and nuclei were re-suspended in 50 μl transposition reaction mixture (25 μl TD buffer; 2.5 μl TDE1; and 22.5 μl double-distilled H2O; Nextera DNA sample preparation kit). Samples were incubated in a ThermoMixer at 1,100 rpm for 45 min at 42°C. DNA was then purified using a QIAGEN MinElute PCR Purification Kit, per manufacturer instructions. DNA was eluted in 10 μl buffer EB. Libraries were then barcoded and amplified with NEBNext High-Fidelity Master Mix and primers previously described (Buenrostro et al., 2013; 50 μl reaction with 10 μl DNA and 2.5 μl of 25 μM primers; one cycle of 5 min at 72°C and 30 s at 98°C; five cycles of 10 s at 98°C, 20 s at 63°C, and 1 min at 72°C). A quantitative PCR on the product determined that an additional 7 cycles (10 s at 98°C, 20 s at 63°C, and 1 min at 72°C) were required. The library was purified and size selected with Agencourt AMPure beads: 45 μl of PCR product was incubated with 18 μl of beads, and supernatant was collected (beads bound larger >2,000-bp fragments). Supernatant (63 μl) was then incubated with an additional 63 μl of beads for 5 min, supernatant was removed, beads were washed twice with 75% ethanol, and DNA was eluted into 50 μl H2O by incubating for 2 min. Samples underwent quality control and quantification on an Agilent BioAnalyzer by the Integrated Genomics Operation core at Memorial Sloan Kettering Cancer Center before pooling and sequencing on an Illumina HiSeq.

Computational analysis

Approximately 30–40 million paired-end ATAC-seq reads from different samples were aligned to the mouse reference genome assembly GRCm38 using STAR (Dobin et al., 2013). The genome was first indexed using the mode as: $STAR --runMode genomeGenerate. To align the reads, the default splicing alignment behavior was disabled. Only uniquely mapped read pairs were reserved in the output BAM files. The command for aligning the reads was as follows: $STAR --runMode alignReads --genomeLoad NoSharedMemory --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --outBAMcompression 6 --outFilterMultimapNmax 1 --outFilterMismatchNoverLmax 0.06 --outFilterMatchNminOverLread 0.55 --outFilterMatchNmin 30 --alignIntronMax 1 --alignEndsType Local --clip3pAdapterSeq CTGTCTCTTA. With these parameters, ~23–32 million reads were uniquely mapped to the accessible chromatin regions. After removal of duplicates by Picard (http://broadinstitute.github.io/picard/), a chromatin accessibility atlas was created using all replicates from all cell types. Briefly, ATAC-seq peaks for individual replicates were called using MACS2 (Zhang et al., 2008) using the following parameters: $macs2 callpeak -g mm -p 1e-2 --nomodel --shift 75 --extsize 200 --keep-dup all --call-summits. Reproducibility of the peaks among replicates of a cell type was evaluated as previously described (Li et al., 2011). After removing the irreproducible peaks (confidence score > 0.05), the chromatin accessibility atlas was constructed by aggregating all the reproducible peaks from the four cell types.

scRNA-seq and computational analysis

Mature CD4+, CD4+CD25+Foxp3, and Foxp3+ thymocytes were sorted from CD4SP-enriched thymocytes pooled from three Foxp3GFP mice as described above. Cells were re-suspended at 106/ml in DPBS with 0.04% BSA.

scRNA-seq and preprocessing

Individual samples were loaded on 10X Genomics Chromium System. Libraries were prepared following 10X Genomics protocols (Chromium Single Cell 3′ Reagent Kits v2 Chemistry) and sequenced on NovaSeq 6000 (Illumina S2 flow cell, paired-end). FASTQ files were processed using the Sequence Quality Control pipeline (Azizi et al., 2018), and reads were aligned to the mouse genome mm10 from ENSEMBL GRCm38. Cells that contained a percentage of mitochondrial transcripts >5% were filtered out, resulting in 4,623 CD4+CD25+ cells, 4,151 Foxp3+ cells, and 4,548 mature CD4+ cells that passed QC metrics, with a median of 3,523 molecules/cell. Cells with total molecule counts of <1,000, determined by lower mode of distribution of total molecules per cell, were additionally filtered out to remove putative empty droplets. The resulting count matrices from all samples were then combined, resulting in a final set of 11,874 cells, and normalized to median library size, where library size is defined as total molecule counts per cell. The normalized data are then log transformed as log (0.1 + counts) for downstream analysis.

t-SNE visualization

For dimensionality reduction and visualization of data, we further excluded genes with very low or very high expression of transcripts (log average expression <0.012 or >3 and a dispersion >0.5), and principle component analysis was performed on log-transformed normalized data. By using 30 principal components, where the number of principal components was determined by the knee-point in eigenvalues, single-cell gene expression was visualized in a two-dimensional projection with t-SNE (van der Maaten and Hinton, 2008), as shown in Fig. 4 A.

Clustering

Clustering of cells was performed using PhenoGraph (Levine et al., 2015) and setting k = 30 nearest neighbors (Fig. 4 C). One cluster of cells was subsequently removed since it was uniformly comprised of cells with a relatively lower number of total molecules compared with other populations. Another cluster was removed because of its disparity with the rest of the data (t-SNE projected this cluster as a separate component that was comprised of cells from all three populations). Significant differentially expressed genes in each cluster were identified using a t test implemented in Single-Cell Analysis in Python (Wolf et al., 2018) with default parameters. The top 30 differentially expressed genes are shown in Fig. 4 E; ribosomal genes were excluded from the heat map due to their higher capture rates in scRNA-seq.

Diffusion component analysis

Due to the continuous structure of the data, we used diffusion maps (Coifman et al., 2005) as a nonlinear dimensionality reduction technique to find the major nonlinear components of variation across cells representing possible trajectories and gradients. We used a Gaussian kernel of perplexity = 30, symmetric Markov normalization, and t = 8 diffusion steps implemented in the tool MAGIC (van Dijk et al., 2018). Ribosomal genes were excluded in computing diffusion components. Fig. 4 F displays the top three diffusion components.

Annotation of diffusion components

To identify processes associated with the nonlinear trajectories across cells, we computed the Pearson correlation between each diffusion component and the expression of predefined gene signatures as described elsewhere (Azizi et al., 2018). The expression of a signature is defined as the mean normalized expression across all genes listed in the signature. In addition to predefined signatures, we also included Hallmark signatures that showed significant enrichment (P < 0.05; false discovery rate < 0.1) in any of the diffusion components. This was achieved by first computing the correlation between each gene and each diffusion component, followed by ranking genes according to the correlation values and performing gene set enrichment analysis (Subramanian et al., 2005) on the preranked gene list per component. For signatures comprised of human genes, the genes were converted to their mouse orthologues. The full set of signatures used for annotation is provided in Table S1.

Analysis of IL-2–expressing cells

For a focused analysis of IL-2–expressing cells, we selected cells with nonzero expression of Il2 (in raw un-imputed data), achieving a subset of 428 cells that are projected using t-SNE in Fig. 5 A. Diffusion component analysis and annotation were performed independently on this subset using the same parameters as above (Fig. 5, C and D).

Data availability

CD4SP thymocyte subset ATAC-seq and scRNA-seq data have been deposited to the Gene Expression Omnibus under accession nos. GSE134366 and GSE134902, respectively.

Online supplemental material

Fig. S1 shows that limiting IL-2 deletion to the DC compartment has no effect on T reg cell development. Fig. S2 shows gating strategy to subset CD4SP thymocytes. Fig. S3 compares immune activation upon IL-2 deletion early or late during T cell development. Table S1 lists gene signatures used in scRNA-seq analysis.

We thank A. Bravo and J. Verter for help with animal husbandry. We thank L. Mazutis and the staff at the Memorial Sloan Kettering Single Cell Sequencing Core facility for performing sample processing and library preparation. We thank all the members of the Rudensky laboratory for technical input and discussion.

This work was supported by the National Institutes of Health/National Cancer Institute Cancer Center Support Grant P30 CA008748, National Institutes of Health grant R37 AI034206, and the Ludwig Center at the Memorial Sloan Kettering Cancer Center. A.Y. Rudensky is an investigator with the Howard Hughes Medical Institute.

The authors declare no competing financial interests.

Author contributions: S. Hemmers and A.Y. Rudensky conceptualized the study, designed experiments, and interpreted all the data. S. Hemmers performed all in vitro and in vivo experiments, with the exception of the ATAC-seq experiment performed by S. Dikiy. M. Schizas and E. Azizi performed bioinformatic analyses and data visualization of the scRNA-seq data. Y. Zhong and S. Dikiy performed bioinformatic analyses and data visualization of the ATAC-seq data. Y. Feng and G. Altan-Bonnet provided unique reagents. S. Hemmers and A.Y. Rudensky wrote the manuscript with input from co-authors. A.Y. Rudensky supervised and acquired funding for the study.

Asano
,
M.
,
M.
Toda
,
N.
Sakaguchi
, and
S.
Sakaguchi
.
1996
.
Autoimmune disease as a consequence of developmental abnormality of a T cell subpopulation
.
J. Exp. Med.
184
:
387
396
.
Azizi
,
E.
,
A.J.
Carr
,
G.
Plitas
,
A.E.
Cornish
,
C.
Konopacki
,
S.
Prabhakaran
,
J.
Nainys
,
K.
Wu
,
V.
Kiseliovas
, and
M.
Setty
.
2018
.
Single-cell map of diverse immune phenotypes driven by the tumor microenvironment
.
bioRxiv.
(Preprint posted April 2, 2018)
Barsoumian
,
H.B.
,
E.S.
Yolcu
, and
H.
Shirwan
.
2016
.
4-1BB signaling in conventional T cells drives IL-2 production that overcomes CD4+CD25+FoxP3+ T regulatory cell suppression
.
PLoS One.
11
:e0153088.
Buenrostro
,
J.D.
,
P.G.
Giresi
,
L.C.
Zaba
,
H.Y.
Chang
, and
W.J.
Greenleaf
.
2013
.
Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position
.
Nat. Methods.
10
:
1213
1218
.
Buenrostro
,
J.D.
,
B.
Wu
,
H.Y.
Chang
, and
W.J.
Greenleaf
.
2015
.
ATAC-seq: A method for assaying chromatin accessibility genome-wide
.
Curr. Protoc. Mol. Biol.
109
:
1
9
.
Burchill
,
M.A.
,
J.
Yang
,
C.
Vogtenhuber
,
B.R.
Blazar
, and
M.A.
Farrar
.
2007
.
IL-2 receptor β-dependent STAT5 activation is required for the development of Foxp3+ regulatory T cells
.
J. Immunol.
178
:
280
290
.
Burchill
,
M.A.
,
J.
Yang
,
K.B.
Vang
,
J.J.
Moon
,
H.H.
Chu
,
C.W.
Lio
,
A.L.
Vegoe
,
C.S.
Hsieh
,
M.K.
Jenkins
, and
M.A.
Farrar
.
2008
.
Linked T cell receptor and cytokine signaling govern the development of the regulatory T cell repertoire
.
Immunity.
28
:
112
121
.
Chinen
,
T.
,
A.K.
Kannan
,
A.G.
Levine
,
X.
Fan
,
U.
Klein
,
Y.
Zheng
,
G.
Gasteiger
,
Y.
Feng
,
J.D.
Fontenot
, and
A.Y.
Rudensky
.
2016
.
An essential role for the IL-2 receptor in Treg cell function
.
Nat. Immunol.
17
:
1322
1333
.
Coifman
,
R.R.
,
S.
Lafon
,
A.B.
Lee
,
M.
Maggioni
,
B.
Nadler
,
F.
Warner
, and
S.W.
Zucker
.
2005
.
Geometric diffusions as a tool for harmonic analysis and structure definition of data: diffusion maps
.
Proc. Natl. Acad. Sci. USA.
102
:
7426
7431
.
D’Cruz
,
L.M.
, and
L.
Klein
.
2005
.
Development and function of agonist-induced CD25+Foxp3+ regulatory T cells in the absence of interleukin 2 signaling
.
Nat. Immunol.
6
:
1152
1159
.
Dobin
,
A.
,
C.A.
Davis
,
F.
Schlesinger
,
J.
Drenkow
,
C.
Zaleski
,
S.
Jha
,
P.
Batut
,
M.
Chaisson
, and
T.R.
Gingeras
.
2013
.
STAR: ultrafast universal RNA-seq aligner
.
Bioinformatics.
29
:
15
21
.
Fontenot
,
J.D.
,
J.P.
Rasmussen
,
M.A.
Gavin
, and
A.Y.
Rudensky
.
2005
a
.
A function for interleukin 2 in Foxp3-expressing regulatory T cells
.
Nat. Immunol.
6
:
1142
1151
.
Fontenot
,
J.D.
,
J.P.
Rasmussen
,
L.M.
Williams
,
J.L.
Dooley
,
A.G.
Farr
, and
A.Y.
Rudensky
.
2005
b
.
Regulatory T cell lineage specification by the forkhead transcription factor foxp3
.
Immunity.
22
:
329
341
.
Kim
,
J.M.
,
J.P.
Rasmussen
, and
A.Y.
Rudensky
.
2007
.
Regulatory T cells prevent catastrophic autoimmunity throughout the lifespan of mice
.
Nat. Immunol.
8
:
191
197
.
Lahl
,
K.
,
C.
Loddenkemper
,
C.
Drouin
,
J.
Freyer
,
J.
Arnason
,
G.
Eberl
,
A.
Hamann
,
H.
Wagner
,
J.
Huehn
, and
T.
Sparwasser
.
2007
.
Selective depletion of Foxp3+ regulatory T cells induces a scurfy-like disease
.
J. Exp. Med.
204
:
57
63
.
Lee
,
H.M.
,
J.L.
Bautista
,
J.
Scott-Browne
,
J.F.
Mohan
, and
C.S.
Hsieh
.
2012
.
A broad range of self-reactivity drives thymic regulatory T cell selection to limit responses to self
.
Immunity.
37
:
475
486
.
Levine
,
J.H.
,
E.F.
Simonds
,
S.C.
Bendall
,
K.L.
Davis
,
A.D.
Amir
,
M.D.
Tadmor
,
O.
Litvin
,
H.G.
Fienberg
,
A.
Jager
,
E.R.
Zunder
, et al
.
2015
.
Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis
.
Cell.
162
:
184
197
.
Li
,
Q.H.
,
J.B.
Brown
,
H.Y.
Huang
, and
P.J.
Bickel
.
2011
.
Measuring reproducibility of high-throughput experiments
.
Ann. Appl. Stat.
5
:
1752
1779
.
Lio
,
C.W.
, and
C.S.
Hsieh
.
2008
.
A two-step process for thymic regulatory T cell development
.
Immunity.
28
:
100
111
.
Liston
,
A.
,
K.M.
Nutsch
,
A.G.
Farr
,
J.M.
Lund
,
J.P.
Rasmussen
,
P.A.
Koni
, and
A.Y.
Rudensky
.
2008
.
Differentiation of regulatory Foxp3+ T cells in the thymic cortex
.
Proc. Natl. Acad. Sci. USA.
105
:
11903
11908
.
Matloubian
,
M.
,
C.G.
Lo
,
G.
Cinamon
,
M.J.
Lesneski
,
Y.
Xu
,
V.
Brinkmann
,
M.L.
Allende
,
R.L.
Proia
, and
J.G.
Cyster
.
2004
.
Lymphocyte egress from thymus and peripheral lymphoid organs is dependent on S1P receptor 1
.
Nature.
427
:
355
360
.
Moran
,
A.E.
,
K.L.
Holzapfel
,
Y.
Xing
,
N.R.
Cunningham
,
J.S.
Maltzman
,
J.
Punt
, and
K.A.
Hogquist
.
2011
.
T cell receptor signal strength in Treg and iNKT cell development demonstrated by a novel fluorescent reporter mouse
.
J. Exp. Med.
208
:
1279
1289
.
Owen
,
D.L.
,
S.A.
Mahmud
,
K.B.
Vang
,
R.M.
Kelly
,
B.R.
Blazar
,
K.A.
Smith
, and
M.A.
Farrar
.
2018
.
Identification of cellular sources of IL-2 needed for regulatory T cell development and homeostasis
.
J. Immunol.
200
:
3926
3933
.
Owen
,
D.L.
,
S.A.
Mahmud
,
L.E.
Sjaastad
,
J.B.
Williams
,
J.A.
Spanier
,
D.R.
Simeonov
,
R.
Ruscher
,
W.
Huang
,
I.
Proekt
,
C.N.
Miller
, et al
.
2019
.
Thymic regulatory T cells arise via two distinct developmental programs
.
Nat. Immunol.
20
:
195
205
.
Pan
,
F.
,
H.
Yu
,
E.V.
Dang
,
J.
Barbi
,
X.
Pan
,
J.F.
Grosso
,
D.
Jinasena
,
S.M.
Sharma
,
E.M.
McCadden
,
D.
Getnet
, et al
.
2009
.
Eos mediates Foxp3-dependent gene silencing in CD4+ regulatory T cells
.
Science.
325
:
1142
1146
.
Popmihajlov
,
Z.
,
D.
Xu
,
H.
Morgan
,
Z.
Milligan
, and
K.A.
Smith
.
2012
.
Conditional IL-2 gene deletion: consequences for T cell proliferation
.
Front. Immunol.
3
:
102
.
Sakaguchi
,
S.
,
N.
Sakaguchi
,
M.
Asano
,
M.
Itoh
, and
M.
Toda
.
1995
.
Immunologic self-tolerance maintained by activated T cells expressing IL-2 receptor alpha-chains (CD25). Breakdown of a single mechanism of self-tolerance causes various autoimmune diseases
.
J. Immunol.
155
:
1151
1164
.
Subramanian
,
A.
,
P.
Tamayo
,
V.K.
Mootha
,
S.
Mukherjee
,
B.L.
Ebert
,
M.A.
Gillette
,
A.
Paulovich
,
S.L.
Pomeroy
,
T.R.
Golub
,
E.S.
Lander
, and
J.P.
Mesirov
.
2005
.
Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles
.
Proc. Natl. Acad. Sci. USA.
102
:
15545
15550
.
Thornton
,
A.M.
, and
E.M.
Shevach
.
1998
.
CD4+CD25+ immunoregulatory T cells suppress polyclonal T cell activation in vitro by inhibiting interleukin 2 production
.
J. Exp. Med.
188
:
287
296
.
van der Maaten
,
L.
, and
G.
Hinton
.
2008
.
Visualizing data using t-SNE
.
J. Mach. Learn. Res.
9
:
2579
2605
.
van Dijk
,
D.
,
R.
Sharma
,
J.
Nainys
,
K.
Yim
,
P.
Kathail
,
A.J.
Carr
,
C.
Burdziak
,
K.R.
Moon
,
C.L.
Chaffer
,
D.
Pattabiraman
, et al
.
2018
.
Recovering gene interactions from single-cell data using data diffusion
.
Cell.
174
:
716
729.e27
.
Weist
,
B.M.
,
N.
Kurd
,
J.
Boussier
,
S.W.
Chan
, and
E.A.
Robey
.
2015
.
Thymic regulatory T cell niche size is dictated by limiting IL-2 from antigen-bearing dendritic cells and feedback competition
.
Nat. Immunol.
16
:
635
641
.
White
,
J.
,
A.
Herman
,
A.M.
Pullen
,
R.
Kubo
,
J.W.
Kappler
, and
P.
Marrack
.
1989
.
The V β-specific superantigen staphylococcal enterotoxin B: stimulation of mature T cells and clonal deletion in neonatal mice
.
Cell.
56
:
27
35
.
Wolf
,
F.A.
,
P.
Angerer
, and
F.J.
Theis
.
2018
.
SCANPY: large-scale single-cell gene expression data analysis
.
Genome Biol.
19
:
15
.
Xing
,
Y.
,
X.
Wang
,
S.C.
Jameson
, and
K.A.
Hogquist
.
2016
.
Late stages of T cell maturation in the thymus involve NF-κB and tonic type I interferon signaling
.
Nat. Immunol.
17
:
565
573
.
Yamamoto
,
M.
,
Y.
Seki
,
K.
Iwai
,
I.
Ko
,
A.
Martin
,
N.
Tsuji
,
S.
Miyagawa
,
R.B.
Love
, and
M.
Iwashima
.
2013
.
Ontogeny and localization of the cells produce IL-2 in healthy animals
.
Cytokine.
61
:
831
841
.
Zemmour
,
D.
,
R.
Zilionis
,
E.
Kiner
,
A.M.
Klein
,
D.
Mathis
, and
C.
Benoist
.
2018
.
Single-cell gene expression reveals a landscape of regulatory T cell phenotypes shaped by the TCR
.
Nat. Immunol.
19
:
291
301
.
Zhang
,
Y.
,
T.
Liu
,
C.A.
Meyer
,
J.
Eeckhoute
,
D.S.
Johnson
,
B.E.
Bernstein
,
C.
Nusbaum
,
R.M.
Myers
,
M.
Brown
,
W.
Li
, and
X.S.
Liu
.
2008
.
Model-based analysis of ChIP-Seq (MACS)
.
Genome Biol.
9
:
R137
.

Author notes

Y. Feng’s present address is Department of Immunology, St. Jude Children’s Research Hospital, Memphis, TN.

This article is distributed under the terms of an Attribution–Noncommercial–Share Alike–No Mirror Sites license for the first six months after the publication date (see http://www.rupress.org/terms/). After six months it is available under a Creative Commons License (Attribution–Noncommercial–Share Alike 4.0 International license, as described at https://creativecommons.org/licenses/by-nc-sa/4.0/).

Supplementary data