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

Abstract

Regulatory T cells (T_reg cells) perform two distinct functions: they maintain self-tolerance, and they support organ homeostasis by differentiating into specialized tissue T_reg cells. We found that epigenetic modifications defined the molecular characteristics of tissue T_reg cells. Tagmentation-based whole-genome bisulfite sequencing revealed more than 11,000 regions that were methylated differentially in pairwise comparisons of tissue T_reg cell populations and lymphoid T cells. Similarities in the epigenetic landscape led to the identification of a common tissue T_reg 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 T_H2 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 T_reg cells integrate multiple waves of epigenetic reprogramming that define their tissue-restricted specialization.

Figure 1. DNA methylome analysis of tissue-resident T_reg 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 T_reg cells and LN T_conv 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 T_reg cells and LN T_conv 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 T_reg (yellow), LN T_conv (green), liver T_reg (black), fat T_reg (purple) and skin T_reg (blue).

Figure 2. Transcriptome analysis of tissue-resident T_reg cells and correlation with epigenetic dataset.

(a) MA plots calculated from RNA sequencing data of fat vs. LN T_reg, skin T_reg vs. LN T_reg cells, fat vs. skin T_reg, and LN T_reg vs. T_conv 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.

Figure 3. Methylation changes of T_reg-specific epigenetic signature.

(a) Methylation mean difference (LN T_reg – LN T_conv) and corresponding log2 RNA expression for promoter and intragenic DMRs identified between LN T_reg and T_conv 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 T_reg (orange line), LN T_conv (green line), Fat T_reg (purple line), Skin T_reg (blue line) and Liver T_reg (grey line) for known T_reg 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 T_reg and T_reg precursor cells (d), LN T_reg, T_conv cells and in vitro induced T_reg cells (iT_reg) (e), tissue-isolated T_reg cells (f), and spleen-derived pT_reg cells, spleen-derived tT_reg cells and splenic Tconv cells (g). Yellow represents unmethylated and blue methylated CpG, while numbers depict quantity of analyzed reads.

Figure 4. Identification of epigenetic and transcriptional changes in tissue-resident T_reg cells.

(a) Methylation mean difference and corresponding log2 RNA expression plot as described in Fig. 3a for DMRs identified in the comparison fat T_reg vs. LN T_reg (upper panel) and skin T_reg vs. LN T_reg (lower panel). (b) Heatmaps of candidates (106 genes) for DMR methylation (left) and gene expression (middle) shown for skin, fat and LN T_reg cells. Right panel, gene expression of candidate genes in T_conv cells isolated from LN, skin and fat. T_reg and T_conv 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 T_reg (yellow) and LN T_conv (green).

Figure 5. Validation of the common tissue T_reg 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 T_reg cells (CD19^–MHCII^–CD3⁺CD8^–CD4⁺CD25⁺Foxp3⁺) and corresponding T_conv 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 T_reg 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 T_reg (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 T_reg (yellow) as in (c). Quantification of CD103 data are based on (e), where contour plots indicate CD103 expression on T_reg 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)

Figure 6. Fat and skin T_reg cells are T_H2-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 T_reg (purple), skin T_reg (blue), LN T_reg (yellow) and LN T_conv (green). (c) Diagnostic T_H2 up- (red, left panel) or downregulated (blue, right panel) signature genes are plotted for the comparison skin vs. LN T_reg (upper panel) and fat vs. LN T_reg (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 T_reg or T_conv 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 T_reg 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 T_reg, skin T_reg, liver T_reg, LN T_conv and LN T_reg 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.

Figure 7. Identification of tisT_regST2 cells.

(a) TisT_regST2 cells (CD8^–CD19^–MHCII^–CD3⁺CD4⁺CD25⁺Foxp3⁺ST2⁺KLRG1⁺GATA-3⁺) population (R1) and identification of this population in different tissues. KLRG1^–ST2^– T_reg cells served as controls (R2). Contour plots and histograms are based on concatenated files of four or more biological replicates. (b) Quantification of tisT_regST2 frequency in different tissues. (c) Frequency distribution of tisT_regST2 in pT_reg and tT_reg population in the colon. pT_reg cells were defined as CD8^–CD19^–MHCII^–CD4⁺CD25⁺Foxp3⁺HELIOS^–RORγt⁺, and tT_reg were defined as CD8^–CD19^–MHCII^–CD4⁺CD25⁺Foxp3⁺HELIOS⁺RORγt^–. tisT_regST2 were identified as KLRG1⁺ST2⁺ of the respective T_reg type. Statistical evaluation with two-tailed unpaired student’s t test (n=10). (d) Spleen cells were stimulated as described. TisT_regST2 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^– T_reg 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) TisT_regST2 from spleen (tisT_regST2 Spl) and ST2^–KLRG1^– T_reg cells (T_reg 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 tisT_regST2 cells derived from spleen. Gene expression of tisT_regST2-associated markers Helios, Gata3, Klrg1, Tigit, and skin associated marker Itgae, Ahr and Gpr55. t-SNE analysis of single tisT_regST2 cells with a skin T_reg signature. Expression of Ahr is indicated by dot size. Means are shown with SD.

Figure 8. Characterization of the tisT_regST2 population.

(a) Frequency of tisT_regST2 in mice. Data are derived from fat, skin, lymph nodes and spleen from mice of 5, 10, 15, 20, and 25 weeks of age. TisT_regST2 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^– T_reg cells were isolated from CD45.2⁺ Foxp3^Cre,YFP donor animals and injected into DT-treated CD45.1⁺ Foxp3^GFP,DTR recipients. After ten days, presence of tisT_regST2 was evaluated by flow cytometry. Contour plots are representative examples of four replicates. (c) Gene-expression analysis of in vitro cultivated and activated tisT_regST2 (fat) and KLRG1^–ST2^– (spleen) T_reg cells. Cultivated for six days with anti-CD3 and anti-CD28 microbeads and IL-2. Gene expression of tisT_regST2 genes was measured by qPCR. Colors indicate gene expression levels with red=high and blue=low relative expression. (d) Flow cytometry anaylsis of T_reg 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 T_reg 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.