Genome-wide DNA-methylation landscape defines specialization of regulatory T cells in tissues.

Regulatory T cells (Treg cells) perform two distinct functions: they maintain self-tolerance, and they support organ homeostasis by differentiating into specialized tissue Treg cells. We found that epigenetic modifications defined the molecular characteristics of tissue Treg cells. Tagmentation-based whole-genome bisulfite sequencing revealed more than 11,000 regions that were methylated differentially in pairwise comparisons of tissue Treg cell populations and lymphoid T cells. Similarities in the epigenetic landscape led to the identification of a common tissue Treg cell population that was present in many organs and was characterized by gain and loss of DNA methylation that included many gene sites associated with the TH2 subset of helper T cells, such as the gene encoding cytokine IL-33 receptor ST2, as well as the production of tissue-regenerative factors. Furthermore, the ST2-expressing population was dependent on the transcriptional regulator BATF and could be expanded by IL-33. Thus, tissue Treg cells integrate multiple waves of epigenetic reprogramming that define their tissue-restricted specialization.

Regulatory T cells (Treg) are critical to maintain self-tolerance. They modulate the functions of different immune cells, thereby affecting a variety of conditions, including autoimmunity, cancer, allergy and inflammation1, 2. In addition, it is becoming increasingly clear that specialized Treg cells in tissues are important to promote organ homeostasis, a function that was initially only attributed to tissue-resident macrophages3. In fat (visceral adipose tissue), Treg cells support metabolic functions and express PPAR-γ, a master-regulator of adipocyte differentiation3, 4, 5, and the IL-33R alpha chain (ST2)6. Other examples of tissue homeostasis promoted by specialized Treg cells include injured skeletal muscles and lungs after influenza A infection7, 8. In both cases, Treg cells present in damaged tissues produce amphiregulin (AREG), an epidermal growth factor receptor ligand important for tissue repair7, 8.

The molecular mechanisms by which tissue-resident Treg cells acquire and stabilize their ‘tissular’ program are poorly understood. Epigenetic modifications have been linked to establishing tissue-resident characteristics in macrophages9, 10. Similar mechanisms could be important to shape the tissue identity of Treg cells.

Our methylome analysis revealed 11,000 differential methylated regions (DMRs) associated with about 4,000 genes. Shared epigenetic profiles led to the identification of a common tissue Treg population, characterized by the epigenetic reprogramming of parts of the T-helper 2 (TH2) pattern and production of the tissue regenerative factor AREG. Our data suggest that epigenetic events shape the characteristics and function of tissue Treg cells.

Results

Identification of differentially methylated regions

To investigate the tissue-specific program of Treg cells, we performed low-input tagmentation-based whole-genome bisulfite sequencing (TWGBS) to decipher the DNA methylome of Treg cells isolated from different tissues. Utilizing Foxp3-reporter mice, we isolated Treg cells from abdominal fat depots, skin, liver and inguinal lymph nodes (LN), and included conventional CD4+ T cells (Tconv) from LN as a control population (Fig. 1a and Supplementary Fig. 1). Three independent replicates per sample were performed and robust data were obtained for all samples with reproducible replicates, with about 7x108 total reads per group and an average 20-fold coverage for each CpG per population (Supplementary Fig. 2a). In pairwise comparisons, a strict definition of at least 30% difference in DNA methylation was chosen and revealed about 11,000 unique differentially methylated regions (DMRs) (Fig. 1b, c and Supplementary Fig. 2b). The average length of a DMR was about 1 kilo base pairs (kb), and annotation with genomic features derived from Refseq illustrated that the majority of DMRs were located in intragenic regions (63%), whereas promoter and intergenic regions comprised only 9% and 28%, respectively (Fig. 1c and Supplementary Fig. 2c,d). We observed a peak of DMRs located immediately downstream of transcription start sites (TSS) (Fig. 1d). Principal component analysis (PCA) showed that the methylation patterns of fat and skin Treg were more similar to each other and rather distinct from LN Treg and Tconv (Fig. 1b, e). These results indicate that Treg cells in tissues have a distinguishable methylation pattern.

Figure 1

DNA methylome analysis of tissue-resident Treg cells.

(a) Circos plot illustrating tagmentation-based whole-genome bisulfite sequencing (TWGBS) methylation data for chromosomes 1-19 and sex chromosomes X and Y for fat, skin, liver, and LN Treg cells and LN Tconv cells. Cumulated methylation values are shown in brackets. Color codes indicate mean methylation and mean methylation difference, respectively. (b) Unsupervised hierarchical cluster of 11,744 differentially methylated regions (DMRs) identified via TWGBS of fat, skin, liver, and LN Treg cells and LN Tconv cells. DMRs were identified as at least 30% differentially methylated in pairwise comparisons. Three replicates per group are shown, with numbers indicating the respective replicate. Colors indicate low (yellow) or high (blue) DMR methylation. (c) Group-wise comparison of DMRs and stratification of DMRs into promoter-resident, intragenic or intergenic based on their genomic location, as shown in more detail in Supplementary Fig. 2d. Numbers are derived from pairwise comparisons of indicated groups. (d) Average distance of DMRs from the transcription start site (TSS) of the closest gene. (e) Principal Component Analysis (PCA) of the different groups based on DMR methylation. Three replicates per group are shown. Colors indicate cell type, with LN Treg (yellow), LN Tconv (green), liver Treg (black), fat Treg (purple) and skin Treg (blue).

RNA transcriptome analysis supports DNA methylation pattern

We performed RNA sequencing (RNAseq) of Treg cells and Tconv cells from tissues. The RNA transcriptome analysis revealed substantial gene expression differences between Treg cells from tissues and LN. 3,072 and 4,698 genes were differentially expressed between fat or skin and LN Treg cells, respectively. The comparison between fat and skin-resident Treg cells showed 552 genes to be differentially expressed (Fig. 2a). Unsupervised hierarchical clustering of the RNAseq data confirmed this notion and grouped fat and skin Treg cells close together, whereas liver Treg located closer to LN Treg cells (Fig. 2b). Hierarchical clustering of DMR methylation data resulted in a very similar grouping, suggesting that gene expression and DNA methylation patterns are interconnected (Fig. 2c). Indeed, integrated analysis of both datasets showed a negative correlation of gene expression and DMR methylation (median -0.303) (Fig. 2d), indicating that, in many cases, hypomethylation of a DMR correlated with the expression of the corresponding gene and vice versa. Such a correlation has been observed in similar studies11.

Figure 2

Transcriptome analysis of tissue-resident Treg cells and correlation with epigenetic dataset.

(a) MA plots calculated from RNA sequencing data of fat vs. LN Treg, skin Treg vs. LN Treg cells, fat vs. skin Treg, and LN Treg vs. Tconv cells. Numerators are plotted as mean of normalized counts. Fold changes are presented as log2. Significantly up- or downregulated genes (p<0.05) are highlighted in red, with the respective numbers shown above or below. Selected genes are labeled. (b) Unsupervised hierarchical clustering of gene expression data. Three replicates per group are shown, with numbers indicating the respective replicate. (c) Unsupervised hierarchical clustering of DMR methylation data. Three replicates per group are shown. (d) Pearson correlation between gene expression and DMR methylation. DMRs were intragenic or within 5kb of the nearest gene. Red line represents the median.

Treg cell-specific epigenetic signature

Methylation-based analyses of candidate regions have previously helped distinguish Treg from Tconv cells. The most prominent example is a Treg-specific demethylated region in the Foxp3 gene, located in the first intron and termed conserved non-coding sequence 2 (CNS2)1, 12. This analysis has been extended by using methylated DNA immunoprecipitation (MeDIP) to analyze differences between Treg and Tconv cells from lymphatic organs13. That study identified a Treg cell-specific CpG hypomethylation pattern that was established in the thymus and included, in addition to Foxp3, other Treg signature genes13. Since our data set included Treg and Tconv cells from LN, we first focused our analysis on this signature established in the thymus. Pairwise comparison between LN Treg and Tconv cells revealed 339 DMRs (Fig. 1c). When plotting the mean methylation difference (LN Treg – LN Tconv) of promoter and intragenically located DMRs against RNA expression data of the corresponding genes, we identified a clear anti-correlation of demethylation being associated with increased gene expression, and gain of methylation with gene repression (Fig. 3a). Our data confirmed Treg-specific hypomethylation at sites described in the earlier study13, e.g. at Ctla4, Ikzf2, Ikzf4, and Il2ra (Fig. 3b), while we also identified several novel hypomethylated sites linked to genes such as Ccr6, Cish and Ift80. Furthermore, we identified previously unappreciated hypermethylated regions in genes that were underrepresented in Treg cells, such as Itk, Satb1, Cox10, Fam78a and Tgfbr2 (Supplementary Fig. 3).

Figure 3

Methylation changes of Treg-specific epigenetic signature.

(a) Methylation mean difference (LN Treg – LN Tconv) and corresponding log2 RNA expression for promoter and intragenic DMRs identified between LN Treg and Tconv cells. Selected demethylated and upregulated genes are highlighted in red, hypermethylated and downregulated genes in blue. Linear regression line in grey. (b) Methylation profile of LN Treg (orange line), LN Tconv (green line), Fat Treg (purple line), Skin Treg (blue line) and Liver Treg (grey line) for known Treg function-related genes Foxp3, Ctla4, Ikzf2, Ikzf4, Il2ra, and Il2rb. Each line represents average methylation values derived from three individual replicates. Little ticks at the bottom of each plot represent location of individual CpGs. Arrows indicate gene direction, black bars gene body regions, red bars annotated promoter regions, blue bars exons. Methylation levels are beta values ranging from 0 (unmethylated) to 1 (methylated). (c) Detailed analysis of the Foxp3 gene with superimposed annotation of introns and exons as well as promoter region (PRO) and conserved non-coding regions 1-3 (CNS). Each circle represents one CpG and the color-code represents degree of methylation from yellow (low) to blue (high). Areas R1-R3 labeled in red represent regions for amplicon-based validation via bisulfite sequencing. (d-g) PCR amplicon sequencing of bisulfite-converted genomic DNA. Thymic Treg and Treg precursor cells (d), LN Treg, Tconv cells and in vitro induced Treg cells (iTreg) (e), tissue-isolated Treg cells (f), and spleen-derived pTreg cells, spleen-derived tTreg cells and splenic Tconv cells (g). Yellow represents unmethylated and blue methylated CpG, while numbers depict quantity of analyzed reads.

Since TWGBS allows resolution at a single CpG level, we used this to study the Foxp3 gene. The Treg-specific demethylation of Foxp3 went far beyond the CNS2 region initially described12, and spanned the entire first intron. Starting from exon 2, the remaining gene was methylated in Treg cells (Fig. 3c and Supplementary Fig. 4). To verify our whole-genome sequencing data with a complementary method, we selected several regions in the Foxp3 gene and performed PCR-based amplification of bisulfite converted DNA, followed by sequencing of these amplicons (Fig. 3d-f and Supplementary Table 1). The amplicon sequencing data confirmed and validated the whole-genome methylation data and, additionally, established that demethylation of the entire intron 1, as well as an upstream region 1, could occur during differentiation of Treg cells in the thymus (Fig. 3d and Supplementary Fig. 5). Furthermore, this Foxp3 pattern was shown to be present in naive Treg cells, but not in in vitro transforming growth factor-β (TGF-β)-induced (iTreg) cells (Fig. 3e and Supplementary Fig. 6). In addition, methylation at the Foxp3 gene was not responsible for differences in Foxp3 expression levels (Supplementary Fig. 7).

To further investigate the universality of this Treg-specific methylation signature, we isolated peripherally-induced (pTreg) and thymus-derived Treg cells (tTreg). pTreg cells are characterized by the expression of the transcription factor RORγt and do not express the transcription factor HELIOS (Ikzf2), whereas tTreg cells do not express RORγt, but express HELIOS14. We sorted both populations from the spleen and performed PCR-based amplification of bisulfite converted DNA for regions in Ctla4, Ikzf2, Il2ra and Foxp3. Although Ctla4, Il2ra and Foxp3 were uniformly hypomethylated in both Treg populations, Ikzf2 was completely methylated in pTreg cells, identifying differences in methylation of Ikzf2 as an epigenetic mark to differentiate between pTreg and tTreg cells (Fig. 3g and Supplementary Fig. 8).

In summary, TWGBS is a powerful method to study Treg cell-specific methylation differences at a single CpG level.

Epigenetic landscape of tissue Treg cells

While pairwise comparison between LN Treg and Tconv cells identified 339 DMRs, the number of DMRs between fat or skin Treg versus LN Treg was about 5-fold larger (1,593 and 1,645 DMRs, respectively) (Fig. 1c). Many of the DMRs were shared between fat versus LN and skin versus LN Treg cells, indicating either common effector/memory or specific tissue Treg programs (Fig. 4a). We extracted 106 genes that showed differential methylation and corresponding gene expression changes in both comparisons (Fig. 4b). To differentiate between common effector/memory and tissue Treg programs, we compared this signature with RNAseq data from Tconv cells isolated from the same peripheral tissue (Fig. 4b and Supplementary Fig. 9). Although individual genes of this list could be identified as common effector/memory related, e.g. Foxp1 or Lef1, the majority of the 106 genes were part of tissue Treg programs and not just common effector/memory related. Pairwise comparisons of Treg and CD4 non-Treg cells from the same tissue revealed that the signature of 106 genes was highly significantly biased towards tissue Treg cells (Supplementary Fig. 9b).

Figure 4

Identification of epigenetic and transcriptional changes in tissue-resident Treg cells.

(a) Methylation mean difference and corresponding log2 RNA expression plot as described in Fig. 3a for DMRs identified in the comparison fat Treg vs. LN Treg (upper panel) and skin Treg vs. LN Treg (lower panel). (b) Heatmaps of candidates (106 genes) for DMR methylation (left) and gene expression (middle) shown for skin, fat and LN Treg cells. Right panel, gene expression of candidate genes in Tconv cells isolated from LN, skin and fat. Treg and Tconv gene expression data were row-normalized together. Color codes indicate high (blue) or low (yellow) methylation (left panel) or high (red) and low (blue) gene expression (right panel). (c) Methylation profiles of indicated genes from (b) are plotted. Plots are as described in Fig. 3b. Each line represents average methylation values derived from three individual replicates. Plotted are skin (blue), fat (purple), and LN Treg (yellow) and LN Tconv (green).

In many cases, methylation at the promoter sites was similar between the groups, but started to be differential just after the TSS in the first intron, as observed in the Foxp3 gene (Fig. 4c and Supplementary Fig. 10). For selected genes, we confirmed differential expression at the protein level via flow cytometry (Fig. 5a). For example, Klrg1 displayed hypomethylation in tissue Treg cells, and was expressed by more than 80% of Treg cells in the skin and fat, but only by about 10% of Treg cells from LN. Hypomethylation and expression of TIGIT and ST2 behaved similarly. The strong changes in protein expression observed for TIGIT, KLRG1 and ST2, were not found in Tconv cells isolated from skin and fat (Fig. 5a).

Figure 5

Validation of the common tissue Treg signature and identification of tissue-specific patterns.

(a) Flow cytometry analysis of TIGIT, KLRG1, ST2, BCL-2, TCF7 and LEF1 in fat, skin and LN Treg cells (CD19–MHCII–CD3+CD8–CD4+CD25+Foxp3+) and corresponding Tconv cells (CD19–MHCII–CD3+CD8–CD4+CD25– Foxp3–). Contour plots are concatenated files representative of four replicates, quantification at bottom part. Individual mice are shown (n=4). Statistical evaluation based on one-way ANOVA with Bonferoni post-test (***=p<0.001). (b) Differences in DMR methylation and gene expression of the comparison Treg cells from fat vs. skin as described in Fig. 3a. (c) Methylation pattern of the gene Pparg and corresponding gene expression. Shown are skin (blue), fat (purple), and LN Treg (yellow). Each line represents one individual replicate (n=3). Gene expression is plotted as reads per kilobase per million mapped reads (Rpkm). Significance based on RNA sequencing calculations as described in methods section and indicated by asterisks. (d) Methylation of Itgae (CD103) and Gpr55 in skin (blue), fat (purple), and LN Treg (yellow) as in (c). Quantification of CD103 data are based on (e), where contour plots indicate CD103 expression on Treg cells from LN, fat and skin measured via flow cytometry. Statistical evaluation based on one-way ANOVA with Bonferoni post-test (n=9). Means are shown with standard deviation (SD)

In addition to shared characteristics, we were also interested in differences between tissue Treg cells (Fig. 5b). The comparison between fat and skin Treg cells revealed that the key transcription factor of visceral adipose tissue Treg differentiation, Pparg, was hypomethylated in fat Treg cells and, concomitantly, the gene was also highly expressed (Fig. 5c). On the other hand, skin Treg cells had several specifically hypomethylated gene loci, including Ahr, Icos, Itgae (CD103) and Gpr55, and all corresponding genes were overexpressed in skin Treg cells (Fig. 5b, d, e and Supplementary Fig. 11).

Overall, these results showed that the epigenetic landscape of fat and skin Treg cells shared many characteristics.

Epigenetic reprogramming of TH2-associated loci

To further delineate common principles of tissue Treg cells, a gene ontology (GO) analysis of the 106 genes with differential methylation and expression patterns (Fig. 4b) was performed. At the top of the list were GO terms describing T cell differentiation (Fig. 6a). Four transcription factors were particularly enriched in these lists, namely Gata3, Irf4, Rora, and Batf, all of which were hypomethylated and overexpressed in fat and skin Treg cells (Fig. 6b and Supplementary Fig. 12). GATA-3 and IRF4 are key transcription factors determining the TH2 fate decision of CD4+ T cells15. To probe whether fat and skin Treg cells show a generalized type-2 profile, a selected list of diagnostic TH2 probes16 was used to plot skin Treg versus LN Treg and fat Treg versus LN Treg cells. In both cases, TH2 overexpressed as well as TH2 underrepresented genes were significantly biased towards the respective tissue Treg site, indicating a type-2 polarization (Fig. 6c). For example, Il1rl1 (encoding for ST2) and Il10, both associated with type-2 conditions15, 17, were hypomethylated in tissue Treg cells. Il10 displayed two hypomethylation regions, one about 9 kb upstream of the promoter and a second intragenic region. The upstream region falls into a previously described DNaseI hypersensitive site of the Il10 locus in TH2 cells18. The Il1rl1 (ST2) gene was hypomethylated in the first intron, and about 90% of the fat Treg and 60% of skin Treg cells expressed ST2, as compared to less than 10% in the spleen counterpart (Fig. 6d and Supplementary Fig. 13). If tissue Treg cells from fat and skin were type-2 biased, it should be possible to recapitulate parts of their phenotype by treating lymphoid Treg cells with IL-4. Indeed, IL-4 treatment of lymphoid Treg cells strongly induced the expression of Gata3, Irf4, Il1rl1 and Il10, but repressed Ms4a4b, a TH1-associated gene, underrepresented in tissue Treg cells, in a dose-dependent manner (Fig. 6e and Supplementary Fig. 14).

Figure 6

Fat and skin Treg cells are TH2-like polarized.

(a) Gene Ontology (GO) term analysis of the gene list described in Fig. 4b. (b) Methylation profile and protein (GATA-3) or gene expression for Gata3, Irf4, Rora and Batf as described in Fig. 3b. Average methylation values derived from three individual replicates are shown for fat Treg (purple), skin Treg (blue), LN Treg (yellow) and LN Tconv (green). (c) Diagnostic TH2 up- (red, left panel) or downregulated (blue, right panel) signature genes are plotted for the comparison skin vs. LN Treg (upper panel) and fat vs. LN Treg (lower panel). Numbers indicate the number of genes per site. Individual genes are highlighted. Chi2-test for proportions. (d) Methylation profile and protein or gene expression of Il1rl1 and Il10 as described in (b). ST2 expression quantification based on Fig. 5a and statistically evaluated with one-way ANOVA (n=4). (e) Culture of Treg or Tconv cells with 25 ng/mL IL-4 (+) or control (-) followed by qPCR-based evaluation of target gene expression. Statistical evaluation with unpaired two-tailed student’s t test (n=4). (f) Tissue Treg methylation data were grouped into 4 specific methylation-expression clusters (top panel). Based on these clusters, binding motif enriched factors were calculated and illustrated in a heat map (middle panel). Gene expression values of cluster 2 binding motif enriched factors for fat Treg, skin Treg, liver Treg, LN Tconv and LN Treg are shown as heatmap (bottom panel), where colors indicate relative expression (high=red; low=blue). Individual replicates are shown (n=3). Means are shown with SD.

Demethylation at CG sites can allow transcription factors (TF) to modulate gene transcription19. We studied the TF binding motifs that were enriched in hypomethylated DMRs found in fat and skin Treg cells, and represented in cluster 2 (Fig. 6f; upper panel). Among the twelve TF with an enrichment of their binding motif in cluster 2, eleven were also overexpressed in fat and skin Treg cells, including Rora, Batf, and different members of the Jun and Fos family (Fig. 6f; lower panel). The BATF-JUN complex promotes DNA binding of the transcription factor IRF4 and the IRF4–JUN–BATF heterotrimeric complex was shown to be critical for IRF4-mediated transcription in T cells20, thereby, the complex may reinforce type-2 polarization.

Characterization of tissue ST2+ Treg cells

We called the here identified TH2-biased subset of Treg cells, expressing ST2 and dominating the Treg population in fat and skin tissue, ‘tissue Treg ST2 cells’ (tisTregST2). They are characterized by epigenetic and gene expression differences of 106 genes (Fig. 4b), such as Gata3, Irf4, Batf, Rora, Maf, Il1rl1, Il10, CD200r1, Tigit, and Klrg1. Fat and skin Treg cells showed demethylation and overexpression of Maf (c-maf), an originally described TH2 cell-associated TF able to bind to the Il10 promoter and induce Il10 transcription18. In addition, one of the fundamental characteristics of tisTregST2 cells is the expression of AREG, a TH2-associated epidermal growth factor receptor ligand21. The Areg gene locus harbored two hypomethylated regions upstream of the promoter in fat and skin Treg cells (Supplementary Fig. 15).

Since expression of ST2, KLRG1, TIGIT and GATA-3 characterized tisTregST2 cells in fat and skin, we used these markers to screen a number of organs for the presence of this cell type. While fat and skin had the highest fraction of tisTregST2 within the Foxp3-positive Treg compartment (about 80-90% and 50-60%, respectively), other peripheral organs such as lung, bone marrow and liver contained 10-20% of them, and lymphoid organs had the lowest fraction with less than 5%, (Fig. 7a, b). We could not detect a ST2+KLRG1+ population in the Tconv compartment in skin, liver, blood, BM, lung, and only a minor (< 5%) population in fat tissue (Supplementary Fig. 16a). While tisTregST2 cells from the different organs similarly expressed, for example, high levels of GATA-3, the corresponding Tconv population did not resemble tisTregST2 characteristics (Fig. 7a and Supplementary Fig. 9b, 16b, c).

Figure 7

Identification of tisTregST2 cells.

(a) TisTregST2 cells (CD8–CD19–MHCII–CD3+CD4+CD25+Foxp3+ST2+KLRG1+GATA-3+) population (R1) and identification of this population in different tissues. KLRG1–ST2– Treg cells served as controls (R2). Contour plots and histograms are based on concatenated files of four or more biological replicates. (b) Quantification of tisTregST2 frequency in different tissues. (c) Frequency distribution of tisTregST2 in pTreg and tTreg population in the colon. pTreg cells were defined as CD8–CD19–MHCII–CD4+CD25+Foxp3+HELIOS–RORγt+, and tTreg were defined as CD8–CD19–MHCII–CD4+CD25+Foxp3+HELIOS+RORγt–. tisTregST2 were identified as KLRG1+ST2+ of the respective Treg type. Statistical evaluation with two-tailed unpaired student’s t test (n=10). (d) Spleen cells were stimulated as described. TisTregST2 cells were identified as CD45+TCRβ+CD4+CD8–Foxp3+KLRG1+ST2+ cells (R2, red) and stained for intracellular expression of AREG and IL-10. KLRG1–ST2– Treg cells (R1, blue) were used as control. Right panel, quantification. Statistical evaluation with one-way ANOVA and Bonferoni post-test (n=4). Additional controls shown in Supplementary Fig. 19. (e) TisTregST2 from spleen (tisTregST2 Spl) and ST2–KLRG1– Treg cells (Treg Spl) from spleen were isolated using FACS and expression of different genes was analyzed via qPCR. Color code depicts gene expression value (red=high, blue=low). (f) Methylation of CG dinucleotides at the Gata3, Klrg1 and Lef1 DMR locus is shown for indicated groups, with blue indicating high and yellow indicating low methylation levels. (g) Expression of Grp55 and Pparg, based on the gene expression dataset derived from (e). (h) scRNASeq analysis of 101 tisTregST2 cells derived from spleen. Gene expression of tisTregST2-associated markers Helios, Gata3, Klrg1, Tigit, and skin associated marker Itgae, Ahr and Gpr55. t-SNE analysis of single tisTregST2 cells with a skin Treg signature. Expression of Ahr is indicated by dot size. Means are shown with SD.

To further study the influence of T cell receptor (TCR) signaling and activation on the tisTregST2 population within the Treg compartment in tissues, we separated Treg cells based on different expression levels of CD44, an effector/memory marker. The tisTregST2 population was almost exclusively located in the CD44-high effector/memory compartment (Supplementary Fig. 17a), which is in accordance with previous literature of tissue-resident Treg cells being of effector/memory phenotype22. Since high expression of CD44 in tisTregST2 suggested a previous activation event via TCR-signaling, we analyzed this population in Nr4a1 (which encodes Nur77)-GFP reporter mice, where TCR signal strength is measured by reporter activity23. Therefore, we subdivided the Treg pool in GFP-high, GFP-intermediate and GFP-negative fractions to study whether the presence of tisTregST2 cells depends on ongoing TCR signaling. Although the frequency of tisTregST2 was much lower in the lymphatic organs as compared to the fat tissue (1-2% vs. 90%, respectively), the per-organ fraction of tisTregST2 was not influenced whether current TCR-signaling was on or off (Supplementary Fig. 17b).

Next, we wanted to elaborate whether tisTregST2 cells could be part of the induced Treg population in tissues. Previous publications have described that the colon harbors two distinct Treg populations, pTreg and tTreg cells14. Unlike HELIOS+RORγt– tTreg cells, pTreg cells are induced by commensal bacteria in the colon. TisTregST2 cells were only present in the tTreg population of the colon, where they represented about 40% of all tTreg cells. The RORγt-positive pTreg compartments in the colon and spleen were completely devoid of tisTregST2 cells (Fig. 7c and Supplementary Fig. 18a). Importantly, tisTregST2 cells located in the colon expressed higher levels of GATA-3 compared to pTreg and ‘non-tisTregST2’ tTreg cells in the same tissue (Supplementary Fig. 18b). In summary, tisTregST2 cells can be identified in virtually every peripheral tissue.

TisTregST2 in the spleen

We investigated whether the small population of tisTregST2 found in lymphatic organs resembled the tisTregST2 pattern found in peripheral tissues. One of the important characteristics that we found in tisTregST2 cells in skin and fat was the production of AREG and IL-10. To this end, we stimulated fat and spleen tisTregST2 cells. About 80% of KLRG1+ST2+ tisTregST2 cells from spleen produced AREG and showed elevated levels of IL-10 when compared to spleen KLRG1–ST2– Treg cells, and thereby demonstrated a very similar effector profile comparable to tisTregST2 cells isolated from fat tissue (Fig. 7d and Supplementary Fig. 19a, b). CD8+ T cells and Tconv did not produce AREG and IL-10 under the same conditions (Supplementary Fig. 19a). To further validate the similarities, we sorted the tisTregST2 population from spleen and analyzed the characteristic tisTregST2 profile, including additional genes that showed epigenetic changes in fat and skin Treg populations such as lamin A (Lmna) or oxysterol binding protein like 3 (Osbpl3). All analyzed genes were differentially expressed in the tisTregST2 fraction from spleen compared to the global splenic Treg pool, including Gata3, Rora and Irf4, matching the profile of tisTregST2 cells (Fig. 7e). Especially the effector molecules Il10 and Areg were highly overrepresented with 65-fold and 16-fold difference, respectively, validating our functional protein expression data. To investigate whether the tisTregST2 characteristic methylation changes could also be observed in the spleen counterpart, we sequenced PCR-based amplicons of bisulfite converted DNA based on DMRs identified in our whole-genome approach. Genes specifically demethylated in fat and skin tisTregST2 cells, such as Gata3 or Klrg1, were also hypomethylated in splenic tisTregST2 cells, but not splenic Klrg1–ST2– Treg cells (Fig. 7f).

To further dissect the tisTregST2 cells found in lymphoid organs, we first looked at individual marker expression. Virtually all Treg cells in the skin expressed high levels of CD103, whereas less than 5% of fat-resident tisTregST2 cells expressed this marker Analysis of tisTregST2 cells from spleen showed that about 40% express CD103, indicating a heterogeneity of the tisTregST2 population in this organ (Supplementary Fig. 19c). Induction of PPAR-γ in the fat tisTregST2 population represented an additional functional tissue-based adaptation (Fig. 5c and Supplementary Fig. 20), comparable to demethylation and expression of the cannabinoid receptor gene Gpr55 in skin Treg cells (Fig. 5d and Supplementary Fig. 20). When looking at the markers Pparg and Gpr55, which indicate tissue-restricted adaption of the tisTregST2 pool in fat and skin, we found both to be highly over-represented in the spleen tisTregST2 population (>30-fold) (Fig. 7g). To extend this analysis, we performed single cell RNA sequencing (scRNAseq) of tisTregST2 cells isolated from spleen. As expected, all single cells expressed Helios, Gata3 and Klrg1. In contrast to this, genes that were biased towards the skin-resident tisTregST2 population, like Itgae (CD103), Ahr or Gpr55 were only expressed in a fraction of individual cells, indicating the presence of subgroups of tisTregST2 in the spleen, probably representing individual tissue characteristics (Fig. 7h). This suggests that fat and skin-resident Treg cells contain a recirculating fraction present in the lymphatic tisTregST2 pool.

TisTregST2 are distinct

To analyze whether tisTregST2 are a distinct differentiation state of Treg cells, we first analyzed the presence of these cells over time. Our data showed that the fraction of tisTregST2 cells among Treg cells in different tissues was stable over a time period of 5 to 25 weeks of age, indicating homeostasis of the tisTregST2 compartment within tissues (Fig. 8a and Supplementary Fig. 21a, b). To understand the origin of tisTregST2 cells, we studied whether lymphoid organ-derived, tisTregST2-depleted Treg cells could be the precursor of tisTregST2 cells found in tissues. To generate space in the Treg compartment, host Treg cells were depleted by injecting diphtheria toxin into Foxp3-DTR mice. Congenically marked KLRG1–ST2– Treg cells were injected and skin and fat tissue was analyzed for the presence of tisTregST2 cells after 10 days. These experiments showed that lymphoid organ Treg cells have the ability to seed the peripheral tissues and differentiate into tisTregST2 cells (Fig. 8b).

Figure 8

Characterization of the tisTregST2 population.

(a) Frequency of tisTregST2 in mice. Data are derived from fat, skin, lymph nodes and spleen from mice of 5, 10, 15, 20, and 25 weeks of age. TisTregST2 were gated as described in Fig. 7. Additional plots are shown in Supplementary Fig. 21. Contour plots are concatenated files representative of five replicates. (b) CD45.2+KLRG1–ST2– Treg cells were isolated from CD45.2+ Foxp3Cre,YFP donor animals and injected into DT-treated CD45.1+ Foxp3GFP,DTR recipients. After ten days, presence of tisTregST2 was evaluated by flow cytometry. Contour plots are representative examples of four replicates. (c) Gene-expression analysis of in vitro cultivated and activated tisTregST2 (fat) and KLRG1–ST2– (spleen) Treg cells. Cultivated for six days with anti-CD3 and anti-CD28 microbeads and IL-2. Gene expression of tisTregST2 genes was measured by qPCR. Colors indicate gene expression levels with red=high and blue=low relative expression. (d) Flow cytometry anaylsis of Treg cells in tissues of Batf-deficient mice (Batf–/–) versus Batf-sufficient wildtype mice (Batf+/+). Dot plots are concatenated files representative of six to ten replicates. Individual mice are shown. Statistical evaluation based on two-tailed unpaired student’s t test (n=6-10). (e) Flow cytometry analysis of Treg cells in tissues of wildtype mice treated with of IL-33 or PBS. Dot plots are concatenated files representative of four replicates. Additional plots are shown in Supplementary Fig. 21. Individual mice are shown. Statistical evaluation based on two-tailed unpaired student’s t test (n=4). Means are shown with SD.

To study the stability of already differentiated tisTregST2 cells, we cultured fat-derived tisTregST2 under well-defined conditions with anti-CD3+CD28 beads and IL-2 for six days in vitro. TisTregST2 cells showed a very stable expression pattern of characteristic markers such as Il1rl1, Pparg, Osbpl3, Batf, Tigit, Gata3, and Klrg1 comparable to freshly isolated tisTregST2 from fat. In parallel, cultivated ‘non-tisTregST2’ Treg cells did not upregulate these characteristic genes, indicating that expression of these genes is not a mere function of Treg activation and that the tisTregST2 program is not a temporary state, but a stable program (Fig. 8c).

In order to identify an essential transcription regulator for tisTregST2 cells, we focused on BATF. As described earlier, our data showed that Batf is over-expressed in tisTregST2 cells, its gene locus is heavily hypomethylated in tissue Treg cells from fat and skin, and we identified an enrichment of BATF DNA binding motifs in regions of genes that were specifically hypomethylated in tissue Treg cells. Indeed, BATF deficient mice showed severely reduced numbers of tisTregST2 in all organs analyzed, including skin, fat, lung, bone marrow, LN and spleen, while numbers of other Treg cells in the same tissues were not reduced. These data show that BATF is an essential transcriptional regulator of tisTregST2 cells (Fig. 8d).

As tisTregST2 cells are characterized by the expression of ST2, we wanted to understand if IL-33 could act as a growth factor amplifying the tisTregST2 pool in-vivo. Administration of recombinant IL-33 substantially expanded the tisTregST2 population in all organs tested. Numbers of tisTregST2 cells were increased by 10-fold in fat, 5-fold in skin, 13-fold in liver, and 60-fold in lung tissue (Fig. 8e). This expansion did not change the identity, as GATA-3 was still over-expressed in expanded cells as compared to ‘non-tisTregST2’ Treg cells isolated from the same organs (Supplementary Fig. 21c). Collectively, these data showed that tisTregST2 cells are a distinct state. They require the transcription factor BATF, and can be expanded via IL-33 signaling in situ.

Discussion

The present study provides evidence that tissue Treg cells undergo extensive epigenetic reprogramming. Changes in the methylome can be used to determine the underlying functional programs. The similarities in the epigenetic landscape between fat- and skin-resident Treg cells allowed us to identify a common tissue Treg population, characterized by the expression of KLRG1 and ST2, a TH2-like program, and the production of tissue regenerative factors.

Classically, Treg cells were viewed as regulators of other immune cells. With the characterization of fat Treg cells, this notion has been extended to a second critical function, supporting organ homeostasis3. In visceral adipose tissue, about 80-90% of Treg cells represent the here described tisTregST2 phenotype. In skin, this number was somewhat lower at about 50-60% of Treg cells, and in lung and liver between 10-20%. These findings explain why fat and skin-resident Treg cells shared a closer relation in the DNA methylome analysis as compared to liver Treg cells. It is intriguing that this TH2-biased tisTregST2 subset is present in virtually every organ. ST2, as well as the transcription factors BATF and IRF4, were shown to be required for fat Treg cell differentiation6. Based on our data, we could show that BATF and IL-33 are not just important for fat Treg cells, but for tisTregST2 in all tissues, extending the perspective from an adipose centered view to a global one.

Our data also show that tissue Treg cells integrate epigenetic changes from multiple differentiation steps. The first specific epigenetic reprogramming occurs during thymic differentiation and stabilizes the universal Treg identity13. The second line of epigenetic modifications solidifies the functional ‘tisTregST2’ specialization via selective hypomethylation of a signature that includes more than 100 genes. Thereby, the tisTregST2 population acquires a unique reprogramming landscape. On top of this tisTregST2 specialization program, we find tissue-specific epigenetic reprogramming. In the fat Treg population, we identified, among others, methylation differences in the Pparg gene. In addition, skin Treg cells revealed several interesting epigenetic differences, where Ahr and Gpr55 could be of specific relevance. Ahr signaling is important for immune cells and their function at barrier organs such as the skin24. Grp55, as a cannabinoid receptor, is associated with algesia linked to inflammatory and neuropathic pain, and could enable nociception by skin Treg cells in an organ with a strong pain perception25, 26.

Where does the epigenetic reprogramming take place? Are most organs independently able to induce the common tisTregST2 program, and, in addition, add organ specific flavors? Alternatively, one or very few organs could induce the tisTregST2 reprogramming and, by circulation, cells reach the individual tissues where they further specialize. In the latter, the fat tissue could be such a candidate as the vast majority of Treg cells in this tissue are tisTregST2 reprogrammed. However, parabiosis experiments detected only a low degree of chimerism in the fat tissue Treg compartment, indicating little exchange of tissue Treg cells via the circulation27. It is also possible that tisTregST2 cells found in spleen or LN represent a backup population that can quickly be recruited to support distressed organ homeostasis. In addition, it is very unlikely that local conversion from Tconv into Treg is responsible since our data show that tisTregST2 are only present in the tTreg, but not pTreg fraction, and are demethylated at the Ikzf2 locus.

TisTregST2 do not express the TH2-associated cytokines IL-4 and IL-13, but express IL-10 and AREG. Areg expression could presumably be induced by ST2 signaling, and was shown to be important for tissue repair in the lung8. The tisTregST2 population could therefore represent the prototype of tissue-repair prone Treg cells mediating tissue homeostasis by using the tissue regenerative factor AREG. This could have clinical implications in the adoptive transfer of Treg cells to treat autoimmune and graft versus host diseases28. Deliberate type-2 conditioning by IL-4 and IL-33 during the in vitro expansion of blood-derived Treg cells might amplify their therapeutic potential, especially in respect to supporting tissue repair functions. Indeed, ST2-dependent protective Treg functions could be demonstrated in the colon and adipose tissue6, 29.

As Foxp3 demethylation is used as a diagnostic marker to detect Treg cells12, 30, peripheral reprogramming events can be used to study the functional capacity of Treg cells. We identified thousands of DMRs with single CpG resolution, characterizing the universal Treg identity and the peripheral reprogramming. These differences will allow the design of diagnostic probes for amplicon-based sequencing to follow the origin and cell fate of Treg cells in various pathological conditions. Indeed, the here described difference in methylation of Ikzf2 between tTreg and RORγt-positive pTreg cells is such an example. Analyzing the epigenetic landscape is more than a complementary approach to describe Treg cells; it will help to define Treg cell identities and the permanent underlying molecular programs.

Data availability Statement

The source data are deposited as Fastq files from all TWGBS, RNA-Seq, and single-cell RNA-Seq in the European Nucleotide Archive (ENA) with the accession code PRJEB14591 (http://www.ebi.ac.uk/ena/data/view/PRJEB14591).

Online Methods

Mice

Wildtype C57BL/6, congenic B6.SJL-PtprcaPepcb/BoyCrl (CD45.1+), congenic B6.PL-Thy1a/CyJ (CD90.1+), and Nr4a1-GFP (C57BL/6-Tg(Nr4a1-eGFP/cre)820Khog/J; Jackson Stock Number: 016617)23 mice were obtained from Charles River Breeding Laboratories (Wilmington, MA, USA) or the Jackson Laboratory (Bar Harbor, ME, USA). B6N.129(Cg)-Foxp3tm3Ayr mice (Foxp3.IRES-DTR/GFP)31 were bred to CD45.1+ or CD90.1+ mice in the animal facility of the German Cancer Research Center (DKFZ). B6N.129(Cg)-Foxp3tm3Ayr and B6.129(Cg)-Foxp3tm4(YFP/cre)Ayr/J, Jackson (Foxp3.IRES-YFP/Cre)32 were used to sort YFP or GFP positive Treg cells.

All animals used in this study were male and between 15 and 30 weeks old, unless otherwise indicated. Animals were housed under specific pathogen-free conditions at the DKFZ animal care facility, and the governmental committee for animal experimentation (Regierungspräsidium Karlsruhe, Germany) approved all experiments involving animals.

Tissue digestion for flow sorting of cells

For cell isolation, we used Foxp3GFP or Foxp3YFP reporter mice (Foxp3.IRES-DTR/GFP or Foxp3.IRES-YFP/Cre). T cells were extracted from gonadal visceral adipose tissue (called fat in the study), skin, liver, and inguinal lymph nodes. For lymph nodes, single-cell suspensions were established and red blood cells were lysed. Epididymal visceral abdominal tissue was first mechanically dissected followed by digestion with a buffer containing collagenase II (1 mg/mL), BSA (20 mg/mL) and DNAse (20 µg/mL) for 45 minutes at 37°C in a slow shaking waterbath. Afterwards, the cell suspension was shortly incubated with 0.5 Mol EDTA-H2O, pelleted and further filtered. To isolate cells from skin tissue, back skin area of sacrificed skin was depilated by shaving and hair removal cream, followed by mechanically dissection and digest using a buffer containing collagenase IV (4 mg/mL), FCS (2% vol/vol), and DNAse (10 µg/mL) for 45 minutes at 37°C in a slow shaking waterbath, followed by filtration steps. Liver tissue cells were isolated from perfused animals, mechanically dissected and treated with digestion buffer containing collagenase II (1 mg/mL), BSA (5 mg/mL), and DNAse (20 µg/mL) for 45 minutes at 37°C in a slow-shaking waterbath followed by a Percoll gradient centrifugation step and filtrations. Cells from lung were isolated from perfused animals following a digestion step similar to fat cells (collagenase II (1 mg/mL), BSA (20 mg/mL) and DNAse (20 µg/mL)).

After purification, cells were filtered through a 70µM filter mesh and either stained for flow cytometry-based isolation of target cells or pre-enriched with anti-CD25 magnetic beads followed by column-based isolation. Cells were stained with antibodies as indicated in the next paragraph. Live/dead cell exclusion was performed with a fixable live/dead stain.

Flow sorting and DNA/RNA purification

Stained cells were twice sorted, first “pre-sorted” with moderate purity settings on a BD ARIA II or III cell-sorting machine (4-way sort, “enrich mode”) into flow buffer. Then, cells were re-acquired and sorted again with high purity settings (4-way sort, “4-way purity mode”) directly into either DNA or RNA lysis buffer. Aliquots were sorted into flow buffer for post-sort purity controls, as illustrated in Supplementary Figure 1. Treg cells from tissues were defined as CD3+CD4+CD8–CD45+CD25+Foxp3(GFP)+ population, while Tconv cells were sorted as CD3+CD4+CD8–CD45+Foxp3(GFP)– or as CD3+CD4+CD8–CD19–CD45+CD90.1+Foxp3(GFP)– population. Spleen pTreg were defined as CD4+CD8–CD19–MHCII–CD25+Foxp3(GFP)+HELIOS–RORγt +, while spleen tTreg were defined as CD4+CD8–CD19–MHCII–CD25+Foxp3(GFP)+HELIOS+RORγt–. TisTregST2 cells were defined as CD3+CD4+CD8–CD45+CD25+Foxp3(GFP)+KLRG1+ST2+ population. In the thymus, double negative 1 (DN1) thymocytes were sorted as CD4–CD8–CD25–CD44+, early thymic Treg precursors as CD4+CD8–TCRb+CD69+CD25–Foxp3(GFP)–, late thymic Treg precursors as CD4+CD8–TCRb+ CD69+CD25+Foxp3(GFP)–, and mature thymus Treg cells as CD4+CD8–TCRb+CD69–CD25+Foxp3(GFP)+. From spleen and lymph nodes, we isolated CD44-positive Treg cells (CD4+CD8–CD44+CD25+CD25+ Foxp3(GFP)+), CD44-positive Tconv cells (CD4+CD8–CD44+CD25– Foxp3(GFP)–), naive Treg cells (CD4+CD8–CD44–CD25+ Foxp3(GFP)+), naive Tconv cells (CD4+CD8–CD44–CD25–Foxp3(GFP)–), in vitro induced Treg cells (CD4+CD8–CD25+Foxp3(GFP)+) and in vitro activated Tconv cells (CD4+CD8–CD25+Foxp3(GFP)–). In addition, we sorted Foxp3high, Foxp3intermediate, and Foxp3low expressing Treg cells (CD3+CD4+CD8–CD19–CD25+Foxp3(GFP)high/int/low), CD25negative Foxp3-expressing cells (CD3+CD4+CD8–CD19–CD25–Foxp3(GFP)+) and Tconv cells (CD3+CD4+CD8–CD19–CD25–Foxp3(GFP)-) from spleen and lymph nodes. Genomic DNA was isolated using a gDNA Microprep Kit (Zymo Research) and concentrations were measured with a Qubit® fluorometer. RNA was isolated with the RNEasy mini kit (Quiagen) and concentration was determined with a 2100 Bioanalyzer instrument (Agilent technologies).

Tagmentation-based whole genome bisulfite sequencing (TWGBS)

We applied TWGBS for very low input DNA amounts according to 33 with some modifications. Double stranded pre-adapters consisted of oligo Tn5MErev ([phos]CTGTCTCTTATACACATCT) and either methylated oligo Tn5mC-Apt1 (TcGTcGGcAGcGTcAGATGTGTATAAGAGAcAG) or methylated oligo Tn5mC-Apt2 (GTcTcGTGGGcTcGGAGATGTGTATAAGAGAcAG); the lower case c base indicates 5-methyl cytosine. Pre-adapters were combined at 1:1 ratio to generate a 10 µM load adapter mixture. The transposome was assembled by mixing 12 µL load adapter and 10µL Ez-Tn5 transposase (Epicentre). About 10 ng genomic DNA and 5pg unmethylated λ DNA were used for tagmentation with 1 µl transposome. Tagmented DNA was purified with AmPure beads (Beckman Coulter) and repaired with Bst DNA polymerase (NEB) and 5mC-dNTP mix (Zymo). After a further bead purification, the DNA was bisulfite-converted with the EZ DNA methylation kit (Zymo) according to the manufacturer’s instructions. From each converted DNA sample, four differently barcoded sequencing libraries were generated by PCR (95°C, 3 min; 12 cycles 95°C, 20 sec; 62°C, 15 sec; 72°C, 40 sec) using Kapa 2G Robust HotStart ReadyMix (Kapa Biosystems), SYBRGreen reagent (Life Technologies), primer Tn5mCP1 (AATGATACGGCGACCACCGAGATCTACACTCGTCGGCAGCGTC) and barcoded Tn5mC reverse primers (CAAGCAGAAGACGGCATACGAGAT (- 8 bases barcode-) GTCTCGTGGGCTCGG). Library pools were 100bp paired-end sequenced using Illumina HiSeq 2000.

Ultra-low RNA Sequencing

cDNA was generated and amplified using 0.8 ng of total RNA and SMARTer Ultra Low Input RNA for Illumina Sequencing - HV (Clontech Laboratories, Inc.) according to the manufacturer’s protocol. Then, sequencing libraries were prepared using the NEXT ChIP-Seq Library Prep Master Mix Set for Illumina (New England Biolabs) according to the manufacturer's instructions with the following modifications: The adapter-ligated double-stranded cDNA (10µl) was amplified using NEBNext Multiplex Oligos for Illumina (New England Biolabs, 25 µM primers), NEBNext High-Fidelity 2x PCR Master Mix (New England Biolabs) and 15 cycles of PCR. Final libraries were validated using Agilent 2100 Bioanalyzer (Agilent Technologies) and Qubit flourometer (Invitrogen), normalized and pooled in equimolar ratios. 50bp single-read sequencing was performed on the Illumina HiSeq 2000 v4 according to the manufacturer’s protocol.

Read pair preprocessing of TWGBS data

To determine whether read pairs originated from the original strand in the genome (C converted to T in read 1 and G converted to A in read 2 during bisulfite treatment) or the reverse complementary (G converted to A in read 1 and C converted to T in read 2), we performed the following preprocessing before the alignment: Since the Cs in CpG sites are mainly methylated thus unconverted, they were excluded from the analysis, i.e., CpG, TpG, GpC and GpT dinucleotides were masked. For each read in every read pair, we than calculated the base ratios T/C and A/G denoted as R1T/C, R1A/G and R2T/C, R2A/G. We then compared these ratios between the first and second read in a read pair to determine if they are assigned in (R1-R2) or (R2-R1) order or cannot be assigned to any order. The following rules were applied: (1) Read pairs where the condition R1T/C > R2T/C and R2A/G > R1A/G held true, we assumed the read pair came from the original strand and assigned them the (R1-R2) read order. (2) Read pairs where the conditions R1T/C < R2T/C and R2A/G < R1A/G were met, we assumed the read pair came from the reverse complementary strand and read order (R2-R1) was assigned. (3) Read pairs that did not meet one of the conditions above were ambiguous thus disregarded from the analysis.

Mapping of whole-genome bisulfite sequencing data and methylation calling

The TWGBS data were processed as described in 34: The mm10 reference genome (GRCm38.73) was transformed in silico for both the top strand (C to T) and bottom strand (G to A) using MethylCtools35. Before alignment, adapter sequences were trimmed using SeqPrep (https://github.com/jstjohn/SeqPrep). The first read in each read pair was then C-to-T converted and the second read in the pair was G-to-A converted. The converted reads were aligned to a combined reference of the transformed top (C to T) and bottom (G to A) strands using BWA (bwa-0.6.2-tpx)36 with default parameters, yet disabling the quality threshold for read trimming (-q) of 20 and the Smith-Waterman for the unmapped mate (-s). After alignment, reads were converted back to the original states, and reads mapped to the antisense strand of the respective reference were removed. Duplicate reads were marked, and the complexity determined using Picard MarkDuplicates (http://picard.sourceforge.net/). Reads with alignment scores less than 1 were filtered before subsequent analysis. Total genome coverage was calculated using the total number of bases aligned from uniquely mapped reads over the total number of mappable bases in the genome. At each cytosine position, reads that maintain the cytosine status were considered methylated, and the reads that have cytosine converted to thymine were considered unmethylated. Only bases with Phred-scaled quality score of ≥20 were considered. In addition, the 5 bp at the two ends of the reads were excluded from methylation calling according to M-bias plot quality control. For the TWGBS libraries, the first 9 bp of the second read and the last 9 bp before the adaptor of the first read were excluded from methylation calling.

Calling of differentially-methylated regions (DMRs)

The raw counts of methylated and unmethylated reads for each CpG site from different libraries were merged for each replicate. BSmooth37 was used (default parameters, version 1.2.0) to call differentially methylated regions (DMRs) for all possible ten pairwise tissue-tissue comparisons. Each comparison contained three replicates in each group (3vs3 comparison). To account for possible false positive DMRs reported by BSmooth, due to smoothing in non-informative / low coverage regions, we applied additional filtering procedures. On autosomes we selected for DMRs which had a mean CpG-coverage per DMR greater or equal to five in all six replicates in a particular comparison based on raw read counts. We applied a paired Wilcoxon test on each DMR based on the beta-values of the six replicates. To correct for multiple testing, the resulting p-values have been Benjamini-Hochberg corrected and only DMRs with a p-value lower or equal 0.05 were selected. Due to the fact that all animals in this study were male mice, we expected lower coverage for the sex chromosomes. On allosomes, we pooled the raw counts from the three replicates in each group and selected DMRs with a coverage of >10. All DMRs from all chromosomes were filtered requiring a minimal mean methylation value difference of 0.3.

Principal Component Analysis (PCA) and hierarchical clustering on DMRs

The methylation values of called DMRs from all ten possible tissue-tissue comparisons were joined into a single matrix. Duplicate DMRs (exact same start and end bp position) were eliminated from the matrix. A PCA on all replicates using R (3.1.2) and prcomp was performed. A complete-linkage clustering was performed allowing for clustering of rows (DMRs) and columns (replicates) using R (3.1.2) and pheatmap.

Identification of promoter, intragenic, intergenic regions and TSS

Unique DMRs from all ten comparisons have been annotated (bedtools-2.24.0 closest38) using parameters -d to report distance and -t “first” to handle ties. We used RefSeq September 2013 version and defined promoters as 2000 bp upstream if the gene is located on the plus strand and 500 downstream if located on minus strand. The TSS was defined as the first bp downstream if on plus strand and the first bp upstream if on minus strand. Whole-genome data were visualized in a circos plot as described in39.

Calculation of average distance of DMRs from TSS

To calculate the distance of DMRs from the transcription start site (TSS), we first normalized all DMRs to the same DNA strand, subtracted the mean DMR genomic localization from the TSS coordinates, grouped the results into 500 bp clusters from -10000 to +10000 bp, and quantified the number of DMRs in each cluster. DMRs with more than 10,000 bp distance from the respective TSS were disregarded.

Pearson Correlation between gene expression and methylation

For this analysis, we included DMRs, which were intragenic or within 5kb up/downstream of the nearest gene according to RefSeq annotation. Methylation values and RPKM values were associated with each other using the RefSeq gene identifiers. For each gene we applied a correlation test using Pearson correlation with a cutoff of less or equal -0.7 or greater/equal 0.7 allowing for negative as well as positive correlation.

Hierarchical clustering of DMRs and RNA data

From the DMRs with an absolute Pearson correlation coefficient of bigger or equal 0.7, we created a new matrix of methylation values including all replicates. We applied complete-linkage hierarchical clustering to generate a heatmap using R (v3.1.2) and pheatmap. We further applied R (v3.1.2) and mclust(v5.1)[1] to estimate the number of clusters using parameters: mclust(…, G=1:20). The resulting number of clusters was four and we cut the resulting tree for the DMRs from the complete-linkage using cutree (…, k=4) into four different clusters.

Motif analysis

For each of the four identified clusters we merged overlapping DMRs (bedtools-2.24.0 merge) to avoid possible bias by overrepresentation and extracted the corresponding genomic DNA sequences (bedtools-2.24.0 fasta). All genomic regions in the four clusters where scored with the JASPAR motif library40 using the total binding affinity (TBA) score. Briefly, the TBA score is computed for each genomic region by summing for each position the maximum PWM score between the plus and minus strand41. Then, for each PWM, the regions were ranked according to their TBA score in decreasing order. For each PWM, we determined the recovery curve for the regions in a specific cluster, and the area under the curve (AUC) was computed over the first 300 regions. AUC values were converted into a z-score by computing the mean and standard deviation of the AUC over 1000 randomizations of the ranks. For each set of regions inside a cluster, the z-scores for the PWMs were represented as a heatmap. Only motifs which have an absolute z-score above 3 in one of the clusters are displayed.

Mapping of RNA sequencing data, statistical evaluation and plotting

For all samples, low quality bases were removed with Fastq_quality_filter from the FASTX Toolkit 0.0.13 (http://hannonlab.cshl.edu/fastx_toolkit/index.html) with 90 percent of the read needing a quality phred score > 20. Homertools 4.742 were used for PolyA-tail trimming, and reads with a length < 17 were removed. PicardTools 1.78 (https://broadinstitute.github.io/picard/) were used to compute the quality metrics with CollectRNASeqMetrics. With STAR 2.343, the filtered reads were mapped against mouse genome 38 using default parameters. Count data were generated using HTSeq44 for the genes. For the comparison with DESeq245, the input tables containing the replicates for groups to compare were created by a custom perl script. For DESeq2, DESeqDataSetFromMatrix was applied, followed by estimateSizeFactors, estimateDispersions, and nbinomWald testing. The result tables were annotated with gene information (gene symbol, gene type) derived from the gencode.vM8.gtf file. For the RPKM table of non-coding RNAs, a custom perl script separated non-coding RNA genes from protein coding genes and calculated their RPKM values from the HTSeq count values. Mapping filtered reads against protein coding transcripts using a custom pipeline generated the RPKM table of the protein coding genes. These genes (status “KNOWN”) were extracted from the Mouse EnsEMBL (Rel. 80) database. Mapping was carried out with bowtie2 version 2.2.446 against union mouse genes: every gene is represented by a union of all its transcripts (exons). The count values (RPKM and raw counts) were calculated by running CoverageBed from Bedtools v2.17.038 of the mapped reads together with a specific mouse annotation file for protein coding genes (based on Ensembl 80) in gtf format and parsing the output with custom perl scripts. MA plots were generated as described in47. For hierarchical clustering of RNA data a complete-linkage clustering was performed using R (3.1.2) and pheatmap.

Unsupervised clustering, computing of heatmaps, methylation plotter

Unsupervised hierarchical clustering was performed with Gene Pattern software/hierarchical clustering tool (http://genepattern.broadinstitute.org/). Column and row distance measure were set to Pearson correlation with pair-wise average linkage clustering. No normalizations or transformations have been performed. Heatmaps were generated with gene pattern software/heatmap viewer. Methylation plots (Figure 3C, Supplementary Figure 13c) were generated with methylation plotter48.

Isolation of RNA and reverse transcription followed by qPCR

Sorted cell populations were lysed and RNA was isolated using the RNeasy Mini Kit (Quiagen). Synthesis of cDNA was performed with SuperScript Reverse Transcriptase II and oligo(dT) primers (both Life Technologies) according to manufacturer’s instructions. Real-time PCR was performed using a ViiA7 instrument (Applied Biosystems) and Taqman master mix (Applied Biosystems). Gene expression values were normalized to housekeeping genes (Hprt).

Purification and bisulfite conversion of genomic DNA

Sorted cell populations were resuspended in PBS and genomic DNA was purified according to manufacturers guidelines using the DNEasy Blood and Tissue kit (Quiagen). For pTreg and tTreg isolation, cells were fixed and stained with the Foxp3 Fix/perm buffer set (eBiosciences) according to manufacturer’s instruction. Cells were then sorted into lysis buffer included in the QIAamp DNA micro kit (Quiagen). Reverse cross-linking of gDNA was performed for one hour at 56°C and one hour at 90°C, followed by gDNA isolation. DNA purity and concentration was measured with a NanoDrop or Qubit photometer. Bisulfite-conversion was performed using the EpiTect Bisulfite Conversion Kit (Quiagen) and converted DNA was used immediately after purification or aliquoted and stored at -20°C.

Computation and testing of bisulfite-DNA primers

Genomic DNA was in-silico bisulfite-converted using the Bisulfite Primer Seeker software (http://www.zymoresearch.com/tools/bisulfite-primer-seeker). Primer sequences were calculated based on manufacturers recommendations. Primer pairs were tested on bisulfite-converted genomic DNA to determine optimal annealing temperature range and cycle number for each specific reaction. Once parameters were optimized, adaptor sequences for 454 sequencing and barcodes to distinguish individual samples were attached to each primer pair sequence and synthesized. Alternatively, Illumina adapter sequences were attached to primer pairs for Illumina sequencing. Based on optimal annealing temperature and PCR cycle number, primers were used to generate PCR amplicons from bisulfite-converted DNA for each cell type tested. An overview of primers used for our sequencing experiments is listed in the key resources table. Once genomic DNA had been bisulfite-converted, PCR reactions with bisulfite-specific primers were performed. PCR amplicons were separated from primer dimers on 1%-2% agarose gels and visualized using ethidium bromide. Specific bands were excised under UV light exposure and DNA amplicons were purified using a Quick Gel Extraction Kit (Life Technologies). Equimolar amounts of amplicons were combined and processed on a GS Junior Sequencer (Roche). Sequence reads were aligned to the bisulfite-converted mouse genome and methylation levels were visualized in heat maps. Alternatively, Illumina-tagged amplicons were processed on an Illumina MiSeq V3 machine with Paired-End 300bp or Paired-End 250bp settings. Raw data were aligned to the mouse genome and CpG methylation was calculated. Detailed genomic positions of amplicon data are either labeled in the respective graph or depicted here: Foxp3 R3 (Figure 3g; CG#1: 7,583,950; CG#2: 7,583,986; CG#3: 7,584,002; CG#4: 7,584,036; CG#5: 7,584,063; CG#6: 7,583,950), Ikzf2 (Figure 3g; CG#1: 69,670,284; CG#2: 69,670,291; CG#3: 69,670,370; CG#4: 69,670,377; CG#5: 69,670,385), Il2ra (Figure 3g; CG#1: 11,645,653; CG#2: 11,645,705; CG#3: 11,645,718; CG#4: 11,645,738), Ctla4 (Figure 3g; CG#1: 60,912,472; CG#2: 60,912,521; CG#3: 60,912,536; CG#4: 60,912,573), Gata3 (Figure 7f; CG#1: 9,868,708; CG#2: 9,868,720; CG#3: 9,868,768; CG#4: 9,868,789; CG#5: 9,868,820; CG#6: 9,868,844; CG#7: 9,868,855; CG#8: 9,868,858; CG#9: 9,868,883; CG#10: 9,868,948; CG#12: 9,868,958), Lef1 (Figure 7f; CG#1: 131,116,109; CG#2: 131,116,113; CG#3: 131,116,152; CG#4: 131,116,173; CG#5: 131,116,187; CG#6: 131,116,191; CG#7: 131,116,196; CG#8: 131,116,257; CG#9: 131,116,288; CG#10: 131,116,302; CG#11: 131,116,327, CG#12: 131,116,356).

Flow cytometry

Tissues were isolated and digested as described previously. Single-cell suspension were stained with surface antibodies for 30 minutes at 4°C. If applicable, cells were fixed and permeabilized with the Foxp3 Fix/Perm Buffer Set (eBiosciences) for one hour at RT, followed by intracellular staining. Antibodies used for flow cytometry experiments were conjugated to either Brilliant Violet 411, eFluor 506, Brilliant Violet 605, Brilliant Violet 711, Brilliant UV 395, Brilliant UV 737, PE, AF488/FITC, PE-Cy7, AF647/APC, APC-Cy7, or PerCP-Cy5.5 and are listed in key resources table. Following life dead exclusion dyes were used: Fixable viability dye eFluor506 and eFluor780 (eBiosciences). If applicable, biotin-labeled primary antibodies were stained with fluorchrome-labeled secondary antibodies for 30 minutes at 4°C. Primary Lef1 and Tcf7 antibodies were stained with an anti-rabbit-IgG AF 647-conjugated secondary antibody. Areg and IL-10 stainings were performed as intracellular staining. Samples were acquired on a BD LSRII or a LSR Fortessa 5-laser flow cytometer. Biological replicates were measured separately. Data were analyzed with FlowJo software, and, in figures, concatenated files of biological replicates are shown. Concatenation was performed with FCS Concat (Cytobank).

IL-4 co-culture of Treg cells; iTreg cell induction

Spleen and lymph nodes from Foxp3.IRES-DTR/GFP animals were harvested. Treg cells were pre-enriched with CD25 bead-based positive selection, while Tconv cells were pre-enriched with CD4 bead-based positive selection (autoMACS Pro Separator, Miltenyi Biotec). Pre-enriched Treg and Tconv cells were sorted by FACS to isolate pure KLRG1–ST2– populations from both cell types. Cells were then seeded at 100,000 cells per well with anti-CD3/CD28 beads (Life technologies), IL-12-23p40 blocking mAb, IFNγ blocking mAb, murine IL-2 (Peprotech) and escalating doses of murine IL-4 (Peprotech). Cells were incubated for six days at 37°C and afterwards lysed for immediate RNA isolation, cDNA synthesis and real-time PCR, as described previously.

To generate in vitro induced Treg cells (iTreg cells), Tconv cells were pre-enriched by negative selection with CD8, CD11b, CD11c, CD19, CD25, and CD49b antibodies and magnetic beads. Afterwards, cells were treated with murine TGF-β (Peprotech) and incubated with CD3/28 microbeads for six days at 37°C.

Measurement of IL-10 and AREG production in spleen T cells

Spleens from Foxp3.IRES-DTR/GFP animals were harvested. Single-cell suspensions were established and red blood cells lysed. Afterwards, cells were treated with 1X PMA and Ionomycin cell stimulation cocktail plus transport inhibitor (eBiosciences) or 1X transport inhibitor only. In addition, all samples were subjected to metalloproteinase inhibitor treatment with 10 µM Marimastat (Sigma). Cells were stimulated for 4 hours at 37°C, followed by surface staining, fixation with the Foxp3 Fix/Perm Buffer Set (eBiosciences), and subsequent measurement via flow cytometry.

Single-cell RNA sequencing

To generate the single-cell RNA-sequencing (scRNA-seq) data, we isolated tisTregST2 from pooled spleens, followed by cell capture, cDNA synthesis and amplification on the C1 Single-Cell Auto Prep IFC (Fluidigm). In two separate runs, a total of 127 cells corresponded to single cells (as confirmed by visual inspection of the captured cells). Sequencing libraries were produced with the Illumina Nextera XT kit according to an adopted Fluidigm protocol. All single cells from one C1 run (about 65 cells on average) were pooled and sequenced 1x50 bp reads on an Illumina HiSeq 2000 machine. For each cell, reads were aligned to the murine genome (ERCC sequences concatenated to GRCm38.p4 version 84) with STAR43 version 2.5. On average 70 % of the reads were uniquely mapped. Raw counts were quantified from position-sorted alignment files with HTSeq-count44 using mode 'union' and default quality thresholds of 10. In order to remove bias introduced by low quality data we performed the quality control using the scater package as described in49. Cells were removed as low quality if one or more of the following conditions were met: low library size (0 cells), low number of captured transcripts (2 cells), exceeding mitochondrial content (11 cells) or exceeding ERCC sequences (13 cells). In total 101 cells remained for further analysis. The scran package49 was used to normalize raw counts with deconvolution and to identify HVGs, (with ERCC-fit, span=0.2, FDR=0.05).

Stability of tisTregST2 in vitro

KLRG1+ST2+ tisTregST2 from fat were isolated by FACS as described previously. In addition, KLRG1-ST2- Treg from spleen were isolated by FACS. Samples were frozen for 0h timepoint analysis or tisTregST2 and KLRG1-ST2- cells were incubated with CD3/CD28 microbeads (Dynabeads, Thermo Fisher Scientific) and 5000 U/mL IL-2 for six days at 37°C. Then, cells were lysed and RNA was isolated. Gene expression was determined by qPCR with Taqman probes as described previously. Four biological replicates have been performed.

Transfer of congenically-labeled KLRG1–ST2– Treg cells into Treg-depleted hosts

Treg cells were isolated from CD45.2+ Foxp3YFP, Cre animals with CD25 bead-based pre-enrichment and a reduced staining protocol (CD8–CD19–CD25+Foxp3-YFP+KLRG1–ST2–). Cells were then injected into CD45.1+ recipient animals intravenously. Control animals received PBS injections. In addition, all recipient animals were treated with peritoneal injections of Diptheria Toxin (DT) to eliminate all host-resident Treg cells. DT injections were repeated after 24 hours. Animals were analyzed ten days after Treg/PBS injections have been performed. Four recipient animals were treated with Treg injections from five donor animals. Four recipients were treated with PBS injections. Two animals were used as untreated controls to ensure proper gating.

Expansion of tisTregST2 with IL-33 in vitro

Foxp3GFP animals were injected with either 5μg recombinant murine IL-33 or PBS into the peritoneum on day 1 and day 3 (Biolegend). On day six, mice were sacrificed and frequency and percentage of tisTregST2 in tissues was evaluated by flow cytometry.

Bioinformatics and statistical analysis

Massive parallel sequencing data were statistically tested as described above, and these statistical values were used in Fig. 2a (MA plots), 5c (Pparg), 5d (Gpr55), 6b (Irf4, Rora, Batf), 6d (Il10), and Supplementary Fig. 12, 13a, 13b, 15a, 15b and 20 with n=3 for all comparisons. Data based on flow cytometry or real-time PCR were tested with unpaired two-tailed student’s t tests in Fig. 6e (n=4), 7c (n=10), 7d+g (n=4), 8d (n=6-10), 8e (n=4) and in Supplementary Fig. 18a (n=10), 19a (n=4), 19b (n=4), 21b (n=5), 21c (n=4). One-way ANOVA with Bonferoni post-test was used in Figures 5a (n=4-19), 5d (CD103, n=9), 6b (GATA-3, n=4-19), 6d (ST2, n=4), and Supplementary Fig. 13c (n=4-19), 14b (n=6-7), 16c (n=4-19), 17b (n=4), 19c (n=8). Dunnett posttest to compare all columns vs. control column was utilized in Supplementary Fig. 14a (n=4). All graphs represent the mean of at least three biological replicates ± Standard Deviation (SD). Statistical significance is indicated by asterisks: ***=p<0.001; **=p<0.01; *=p<0.05. To identify the TH2 bias of tissue Treg cells, we used a dataset derived from in vitro differentiated TH1, TH2, TH17, and iTreg cells16. Genes specifically up- or down-regulated in TH2 polarized cells were identified in the comparisons: TH2 vs. TH1, TH2 vs. TH17, and TH2 vs. naïve T cells. TH2-specific genes had to be differential in all three comparisons (>2-fold). We plotted both TH2-specific gene lists (up and down) on our gene expression dataset (Fat Treg vs. LN Treg and Skin Treg vs. LN Treg). The significance of bias was evaluated by chi square testing.

Data availability

Fastq files from TWGBS, RNA-Seq, and single-cell RNA-Seq that support the findings of this study have been deposited in European Nucleotide Archive (ENA) with the accession code PRJEB14591 (http://www.ebi.ac.uk/ena/data/view/PRJEB14591).

Antibodies
Pacific Blue anti-mouse CD3epsilon antibody	Biolegend	AB_2028475
Brilliant Violet 711 anti-mouse CD3 antibody	Biolegend	AB_2563945
APC anti-mouse CD4 antibody	Biolegend	AB_312719
APC/Cy7 anti-mouse CD4 antibody	Biolegend	AB_312699
Biotin anti-mouse CD4 antibody	Biolegend	AB_312711
Brilliant Violet 421 anti-mouse CD4 antibody	Biolegend	AB_11219790
FITC anti-mouse CD4 antibody	Biolegend	AB_312713
Brilliant Violet 711 anti-mouse CD4 antibody	Biolegend	AB_2562099
Brilliant Violet 605 anti-mouse CD4 antibody	Biolegend	AB_2563054
PE anti-mouse CD4 antibody	Biolegend	AB_312715
PE/Cy7 anti-mouse CD4 antibody	Biolegend	AB_312729
PerCP/Cy5.5 anti-mouse CD4 antibody	Biolegend	AB_893326
Brilliant UV 395 anti-mouse CD4 antibody	BD Biosciences	Cat# 563790
Brilliant UV 737 anti-mouse CD4 antibody	BD Biosciences	Cat# 564933
Biotin anti-mouse CD8a antibody	Biolegend	AB_312743
PE/Cy7 anti-mouse CD8a antibody	Biolegend	AB_312761
PerCP/Cy5.5 anti-mouse CD8a antibody	Biolegend	AB_2075238
Biotin anti-mouse/human CD11b antibody	Biolegend	AB_312787
Biotin anti-mouse CD11c antibody	Biolegend	AB_313773
APC/Cy7 anti-mouse CD19 antibody	Biolegend	AB_830707
Biotin anti-mouse CD19 antibody	Biolegend	AB_313639
APC anti-mouse CD25 antibody	Biolegend	AB_312861
Biotin anti-mouse CD25 antibody	Biolegend	AB_312853
PE anti-mouse CD25 antibody	Biolegend	AB_312857
PE/Cy7 anti-mouse CD25 antibody	Biolegend	AB_312865
Brilliant Violet 711 anti-mouse CD25 antibody	Biolegend	AB_2564130
Pacific Blue anti-mouse/human CD44 antibody	Biolegend	AB_493683
Brilliant Violet 421 anti-mouse/human CD44 antibody	Biolegend	AB_10895752
Brilliant Violet 605 anti-mouse/human CD44 antibody	Biolegend	AB_2562451
Brilliant Violet 421 anti-mouse CD45 antibody	Biolegend	AB_10899570
APC/Cy7 anti-mouse CD45 antibody	Biolegend	AB_312981
Pacific Blue anti-mouse CD45 antibody	Biolegend	AB_493535
Biotin anti-mouse CD49b (pan-NK cells) antibody	Biolegend	AB_313411
APC anti-mouse CD62L antibody	Biolegend	AB_313099
APC/Cy7 anti-mouse CD62L antibody	Biolegend	AB_830799
PerCP/Cy5.5 anti-mouse CD62L antibody	Biolegend	AB_2285839
Alexa Fluor 647 anti-mouse CD103 antibody	Biolegend	AB_535952
PE anti-mouse CD103 antibody	Biolegend	AB_1133989
Brilliant Violet 605 anti-mouse CD127 (IL-7Ralpha) antibody	Biolegend	AB_2562114
Brilliant Violet 421 anti-mouse CD127 (IL-7Ralpha) antibody	Biolegend	AB_11218800
PE anti-mouse CD200 R (OX2R) antibody	Biolegend	AB_2074080
PE/Cy7 anti-mouse I-A/I-E antibody	Biolegend	AB_2290801
Pacific Blue anti-mouse I-A/I-E antibody	Biolegend	AB_493527
APC/Cy7 anti-mouse I-A/I-E antibody	Biolegend	AB_1659252
PE anti-mouse/human KLRG1 (MAFA) antibody	Biolegend	AB_10574005
Brilliant Violet 421 anti-mouse/human KLRG1 antibody	Biolegend	AB_2565613
Brilliant Violet 605 anti-mouse/human KLRG1 antibody	Biolegend	AB_2563357
Anti-Mouse ST2 (IL-33R) biotinylated antibody	Affymetrix eBioscience	AB_2572809
Brilliant Violet 421 anti-mouse IL-33R (ST2) antibody	Biolegend	AB_2565634
PE anti-mouse IL-33R (ST2) antibody	Biolegend	AB_2561914
PE anti-mouse TIGIT (Vstm3) antibody	Biolegend	AB_10895760
PE/Cy7 anti-mouse TIGIT (Vstm3) antibody	Biolegend	AB_2565649
Mouse Amphiregulin biotinylated affinity Purified PAb antibody	R and D Systems	AB_2060662
Alexa Fluor 647 anti-Bcl-2 antibody	Biolegend	AB_2274702
Alexa Fluor 488 anti-Bcl-2 antibody	Biolegend	AB_2028390
Anti-Mouse/Rat Foxp3 Alexa Fluor 647 antibody	Biolegend	AB_763538
Anti-Mouse/Rat Foxp3 Biotin antibody	Biolegend	AB_763540
Anti-Mouse/Rat Foxp3 PE antibody	Biolegend	AB_465936
Alexa Fluor 647 anti-GATA3 antibody	Biolegend	AB_2563217
PE anti-GATA3 antibody	Biolegend	AB_2562723
PE anti-mouse IL-10 antibody	Biolegend	AB_315362
Anti-Human/Mouse T-bet PE antibody	Biolegend	AB_925762
LEF1 (C12A5) Rabbit mAb antibody	Biolegend	AB_823558
TCF1 (C63D9) Rabbit mAb antibody	Biolegend	AB_2199302
PE anti-mouse human HELIOS antibody	Biolegend	AB_10660749
Brilliant V 421 anti-mouse Rorγt antibody	BD Biosciences	Cat# 562894
Goat Anti-Rabbit IgG (H+L) antibody, Alexa Fluor 647 Conjugated	Thermo Fisher	AB_10562581
Fixable Viability Dye eFluor 506	Affymetrix eBioscience	Cat# 65-0866-18
Fixable Viability Dye eFluor 780	Affymetrix eBioscience	Cat# 65-0865-18
APC/Cy7 Streptavidin	Biolegend	Cat# 405208
eFluro450 Streptavidin	Affymetrix eBioscience	Cat# 48-4317-82
FITC Streptavidin	Biolegend	Cat# 405201
PE Streptavidin	Biolegend	Cat# 405204
PE/Cy7 Streptavidin	Biolegend	Cat# 405206
PerCP/Cy5.5 Streptavidin	Biolegend	Cat# 405214
Brilliant UV 395 Streptavidin	BD Biosciences	Cat# 564176
Brilliant UV 737 Streptavidin	BD Biosciences	Cat# 564293
APC Streptavidin	Biolegend	Cat# 405207
Brilliant Violet 421 Streptavidin	Biolegend	Cat# 504421
Brilliant Violet 421 Streptavidin	Biolegend	Cat# 504421
Primers & Probes
Taqman Probe for Hprt	Thermo Fisher	Mm01318746_g1
Taqman Probe for Il2ra	Thermo Fisher	Mm01340213_m1
Taqman Probe for Foxp3	Thermo Fisher	Mm00475162_m1
Taqman Probe for Il7r	Thermo Fisher	Mm00434295_m1
Taqman Probe for Pparg	Thermo Fisher	Mm01184322_m1
Taqman Probe for Tbx21	Thermo Fisher	Mm00450960_m1
Taqman Probe for Irf4	Thermo Fisher	Mm00516431_m1
Taqman Probe for Gata3	Thermo Fisher	Mm00484683_m1
Taqman Probe for Il10	Thermo Fisher	Mm01288386_m1
Taqman Probe for Il1rl1	Thermo Fisher	Mm00516117_m1
Taqman Probe for Rora	Thermo Fisher	Mm01173766_m1
Taqman Probe for Areg	Thermo Fisher	Mm00452865_m1
Taqman Probe for Osbpl3	Thermo Fisher	Mm00452865_m1
Taqman Probe for Tcf7	Thermo Fisher	Mm00493445_m1
Taqman Probe for Lef1	Thermo Fisher	Mm00550265_m1
Taqman Probe for Lmna	Thermo Fisher	Mm00497783_m1
Taqman Probe for Ms4a4b	Thermo Fisher	Mm00649916_m1
Taqman Probe for Klrg1	Thermo Fisher	Mm00516879_m1
Taqman Probe for Tigit	Thermo Fisher	Mm03807522_m1
Taqman Probe for Foxp1	Thermo Fisher	Mm00474848_m1
Taqman Probe for Batf	Thermo Fisher	Mm00479410_m1
Taqman Probe for Grp55	Thermo Fisher	Mm02621622_s1
Bisulfite-DNA primer for Foxp3 Upstream R1	Forward Primer	AGGATGTTAGGGTATTAAAGGTTGG
Bisulfite-DNA primer for Foxp3 Upstream R1	Reverse Primer	CCAATTTTCCTAAAACCAACAATAT
Bisulfite-DNA primer for Foxp3 R2 (R2.1)	Forward Primer	AGTGTTTAGTTTTTGTTTTTTTTTTAGGTTTTG
Bisulfite-DNA primer for Foxp3 R2 (R2.1)	Forward Primer	TCTTACRTAACACAAACAAAAAATCAATTAAATAC
Bisulfite-DNA primer for Foxp3 R2 (R2.2)	Reverse Primer	GTTGTGGTATTGTGTTTGGTATATG
Bisulfite-DNA primer for Foxp3 R2 (R2.2)	Forward Primer	ACACTTAATTCAAATAATCAAATTCATAAT
Bisulfite-DNA primer for Foxp3 R3 (TSDR)	Reverse Primer	TGGGTTTTTTTGGTATTTAAGAAAG
Bisulfite-DNA primer for Foxp3 R3 (TSDR)	Forward Primer	AAAAAACAAATAATCTACCCCACAA
Bisulfite-DNA primer for Foxp3 R4 (Exon 8)	Forward Primer	TGAAAGGTTATAATGAAATGATAAGTTTAA
Bisulfite-DNA primer for Foxp3 R4 (Exon 8)	Forward Primer	ATTACCATAACTTCCCCAAAAATAC
Bisulfite-DNA primer for Foxp3 R5 (Exon 11)	Forward Primer	TGATTGTTAATTTTGTTTTTGATTG
Bisulfite-DNA primer for Foxp3 R5 (Exon 11)	Reverse Primer	CAACCTCAATCTCATAATTTTAACC
Bisulfite-DNA primer for Gata3 DMR	Forward Primer	AGAAGGGAAGAAAAATAAAGGAGAGAAAATTATTAG
Bisulfite-DNA primer for Gata3 DMR	Reverse Primer	CCTCCCTCATCATTCTAAATATTACCTC
Bisulfite-DNA primer for Klrg1 DMR	Forward Primer	GAGTTGGGGTAGTAGAGAGTTTTATTTTG
Bisulfite-DNA primer for Klrg1 DMR	Reverse Primer	ACTATACCTCATAAATATATACAATCACTATCTCTCC
Bisulfite-DNA primer for Lef1 DMR	Forward Primer	TTGTTGTTAATTGGGGAAATATTTGAAGTTGTTG
Bisulfite-DNA primer for Lef1 DMR	Reverse Primer	AACRTAAAACTAAAAACAAACAACTAATTATAACC
Bisulfite-DNA primer for Ctla4 Exon 2 (derived from 13)	Forward Primer	TGGTGTTGGTTAGTAGTTATGGTGT
Bisulfite-DNA primer for Ctla4 Exon 2 (derived from 13)	Reverse Primer	AAATTCCACCTTACAAAAATACAATC
Bisulfite-DNA primer for Ikzf2 Intron 3a (derived from 13)	Forward Primer	AGGATGGTTTTTATTGAAGGTGAT
Bisulfite-DNA primer for Ikzf2 Intron 3a (derived from 13)	Reverse Primer	ATACACACCAAACAAACACTACACC
Bisulfite-DNA primer for Il2ra Intron 1a (derived from 13)	Forward Primer	TTTTAGAGTTAGAAGATAGAAGGTATGGAA
Bisulfite-DNA primer for Il2ra Intron 1a (derived from 13)	Reverse Primer	TCCCAATACTTAACAAAACCACATAT
Critical Chemicals
Collagenase II	Sigma Aldrich	Cat# C6885
Collagenase IV	Sigma Aldrich	Cat# C5138
Bovine Serum Albumin BSA	Sigma Aldrich	Cat# A4503
DNAse	Roche	Cat# 11284932001
Ez-Tn5 transposase	Epicentre	Cat# EZI011RK
Bst DNA polymerase	New England Biolabs	Cat# M0275S
5mC-dNTP mix	Zymo Research	Cat# D1030
Power SYBR Green Master Mix	Thermo Fisher	Cat# 4367659
Taqman Gene Expression Master Mix	Thermo Fisher	Cat# 4359016
NEBNext Multiplex Oligos for Illumina	New England Biolabs	Cat# E7335, E7500
NEBNext High Fidelity 2x PCR Master Mix	New England Biolabs	Cat# M0541
SuperScript II Reverse Transcriptase	Thermo Fisher Scientific	Cat# 18064071
Oligo d(T) 12-18 Primer	Thermo Fisher	Cat# 18418012
Anti-biotin Microbeads	Miltenyi Biotec	Cat# 130-090-385
Dynebeads Mouse T-Activator CD3/CD28	Thermo Fisher	Cat# 11456D
Anti-Mouse IL-12/IL-23 p40 Functional Grade Purified 500 ug Antibody	Affymetrix eBioscience	AB_469233
Purified anti-mouse IFN-gamma antibody	Biolegend	AB_315396
Recombinant murine IL-2	Peprotech	Cat# 212-12
Recombinant murine IL-4	Peprotech	Cat# 214-14
Recombinant murine IL-33	Biolegend	Cat# 580506
Recombinant human TGF-beta	Peprotech	Cat# 100-21
PMA/Ionomycin stim cocktail plus transport inhibitors	Affymetrix eBioscience	Cat# 00-4975-03
Marimastat	Sigma Aldrich	Cat# M2699
Critical Commercial Kits
Quick gDNA Microprep Kit	Zymo Research	Cat# D3020
RNEasy Mini Kit	Quiagen	Cat# 74104
Ampure Beads	Beckman Coulter	Cat# A63881
EZ DNA methylation kit	Zymo Research	Cat# D5002
Kapa 2G Robust HotStart ReadyMix	Kapa Biosystems	Cat# KK5701
Smarter Ultra Low Input RNA	Clontech Labs	Cat# 634888
NEXT ChIP-Seq Library Prep Master Mix	New England Biolabs	Cat# E6240
QIAamp DNA micro kit	Quiagen	Cat# 56304
DNEasy Blood and Tissue kit	Quiagen	Cat# 69504
EpiTect Bisulfite conversion kit	Zymo Research	Cat# 59104
PureLink Quick Gel Extraction Kit	Thermo Fisher	Cat# K210012
Foxp3 / Transcription Factor Staining Buffer Set	Affymetrix eBioscience	Cat# 00-5523-00