T-follicular helper cells are epigenetically poised to transdifferentiate into T-regulatory type 1 cells

Abstract

Chronic antigenic stimulation can trigger the formation of interleukin 10 (IL-10)-producing T-regulatory type 1 (TR1) cells in vivo. We have recently shown that murine T-follicular helper (TFH) cells are precursors of TR1 cells and that the TFH-to-TR1 cell transdifferentiation process is characterized by the progressive loss and acquisition of opposing transcription factor gene expression programs that evolve through at least one transitional cell stage. Here, we use a broad range of bulk and single-cell transcriptional and epigenetic tools to investigate the epigenetic underpinnings of this process. At the single-cell level, the TFH-to-TR1 cell transition is accompanied by both, downregulation of TFH cell-specific gene expression due to loss of chromatin accessibility, and upregulation of TR1 cell-specific genes linked to chromatin regions that remain accessible throughout the transdifferentiation process, with minimal generation of new open chromatin regions. By interrogating the epigenetic status of accessible TR1 genes on purified TFH and conventional T-cells, we find that most of these genes, including Il10, are already poised for expression at the TFH cell stage. Whereas these genes are closed and hypermethylated in Tconv cells, they are accessible, hypomethylated, and enriched for H3K27ac-marked and hypomethylated active enhancers in TFH cells. These enhancers are enriched for binding sites for the TFH and TR1-associated transcription factors TOX-2, IRF4, and c-MAF. Together, these data suggest that the TR1 gene expression program is genetically imprinted at the TFH cell stage.

Keywords: T-follicular helper cells; T-regulatory type 1 cells; TFH cells; TR1 cells; epigenetics; immunology; immunoregulation; inflammation; mouse; pMHCII-based nanomedicines; peptide-major histocompatibility complex class II - based nanomedicines; transdifferentiation.

Figure 1.. The TFH-to-TR1 conversion is associated with massive closure of previously open chromatin regions (OCRs) and gene silencing.

BDC2.5mi/I-A^g7-NP-induced Tet+ (TET); KLH-DNP-induced TFH (TFH), and Tconv cells (Tconv) from female nonobese diabetic (NOD) mice (n=4, 8, and 1 mice, respectively) were sorted and processed for 10X GEX+ATAC multiome. (A) Left: multidimensional analysis of scRNAseq and ATACseq data using weighted nearest neighbor (WNN). Cell-type prediction is based on K-means clusterization and predicted cell subtype classification of scRNAseq data. Middle and right: scRNAseq (middle) and scATACseq (assay for transposable-accessible chromatin using sequencing) (right) UMAP (Uniform Manifold Approximation and Projection) plots with K-means clusterization and predicted cell subtype classification based on scRNAseq data. The legends’ colors correspond to different T-cell types and the legends’ shape correspond to sample type. (B) Hierarchical clustering of multiome data based on Euclidean distance of all the clusters found in all sample types. (C) Bar plot with the number of differentially OCRs in TFH.1 (including Tet+ TFH.1 and KLH-DNP-induced TFH.1 cells), TR1/TR1-like cells, as compared to Tconv cells (adjusted p<0.05). Color depicts the proportion of OCRs that are shared with KLH-DNP-induced TFH.1 (vs. Tconv cells) (yellow) or that appear de novo in TR1/TR1-like cells (gray). (D) Volcano plot of Wilcox differential analysis of scRNAseq data from the multiome dataset comparing TFH.1 and TR1/TR1-like cells. Not overlapping dots were labeled. Genes mentioned in the text are boxed. (E–F) Jitter plots depicting log2FC in gene expression between TR1/TR1-like and TFH.1 cells as measured by scRNAseq (adjusted p<0.05) for genes associated with OCRs found in TFH.1 cells (left; closed in TR1/TR1-like cells) or TR1/TR1-like (right; closed in TFH.1 cells) (E), or in both TFH.1 and TR1/TR1-like cells (F) as measured by scATACseq (adjusted p<0.05). Dot (gene) colors define the -log10 (adjusted p) for RNA expression of Wilcox test; dot sizes are proportional to the number of differential OCRs associated with each gene. All genes are labeled except when overlapping. Chromatin closure by TR1 cells was significantly associated with gene downregulation as determined via Chi-square test. TFH, T-follicular helper; TR1, T-regulatory type 1; KLH, keyhole limpet hemocyanin.

Figure 1—figure supplement 1.. High transcriptional similarity between pMHCII-NP- and vaccine-induced TFH.1 cells.

Figure corresponds to an MA plot comparing the transcriptome of the TFH.1 cells from the KLH-induced TFH cell pool and the TFH-like cells contained within the BDC2.5mi/I-A^g7-NP-induced Tet + pool. There were only five differentially expressed genes (|log2FC|>0.5 and adjusted p<0.05; Actb, Ifi27l2a, Inpp4b, Nav2, and Tmsb10) between these two subsets. NP, nanoparticle; pMHCII, peptide-major histocompatibility complex class II; TFH, T-follicular helper; KLH, keyhole limpet hemocyanin.

Figure 2.. Chromatin accessibility in BDC2.5mi/I-A^g7-NP-induced Tet+ and KLH-DNP-induced T-follicular helper (TFH) cells vs. gene expression.

(A) Euler’s plot comparing the number and sharing of differential open chromatin regions (OCRs) between BDC2.5mi/I-A^g7-NP-induced Tet+ and KLH-induced TFH cells relative to their Tconv counterparts (n=4, 3, and 3, respectively) (adjusted p<0.01). (B) Associations between the number and status of differentially open or closed chromatin sites and gene expression as measured via bulk ATACseq (adjusted p<0.01) and RNAseq (adjusted p<0.01 and |log2FC|>2), respectively. Data correspond, from left to right, to BDC2.5mi/I-A^g7-NP-induced Tet+ vs. Tconv and KLH-induced TFH vs. Tconv. Each dot represents a gene, and its size is proportional to the number of associated OCRs. Dot color represents RNAseq differential analysis’ -log10 of false discovery rate (FDR). Only genes with the highest log2FC value in each condition were labeled. (C) Associations between the status of chromatin accessibility in BDC2.5mi/I-A^g7-NP-induced Tet+ vs. KLH-induced TFH cells and gene expression as measured via bulk ATACseq (adjusted p<0.01) and RNAseq (adjusted p<0.01 and. |log2FC|>2), respectively. Each dot represents a gene, and its size is proportional to the number of associated OCRs. Dot color represents RNAseq differential analysis’ -log10 of FDR. Only genes with the highest log2FC value in each condition were labeled. p-Values in (B) and (C) were calculated using Chi-square. KLH, keyhole limpet hemocyanin.

Figure 3.. Genome-wide distribution of H3K4me3, H3K27me3, and H3K27ac marks and active enhancers in BDC2.5mi/I-A^g7-NP-induced Tet+ cells and KLH-DNP-induced T-follicular helper (TFH) cells vs. differential gene expression.

(A) Clustering heatmaps of all the regions enriched for H3K4me3 (left), H3K27me3 (middle), and H3K27Ac (right) deposition in BDC2.5mi/I-A^g7-NP-induced Tet+ cells, KLH-DNP-induced TFH cells, and Tconv cells. The intensity of the red color is proportional to the magnitude of enrichment vs. the corresponding background. (B) Bar plot comparing the relative percentages (X-axis) and absolute numbers (number annotations in each bar) of differentially marked regions for each histone modification (false discovery rate [FDR]<0.01). Top, middle, and bottom rows correspond to KLH-DNP-induced TFH vs. Tconv, BDC2.5mi/I-A^g7-NP-induced Tet+ vs. Tconv, and BDC2.5mi/I-A^g7-NP-induced Tet+ vs. KLH-DNP-induced TFH comparisons, respectively. KLH, keyhole limpet hemocyanin; NP, nanoparticle.

Figure 3—figure supplement 1.. Differential gene deposition of H3K4me3, H3K27me3, and H3K27ac marks in BDC2.5mi/I-Ag7-NP-induced Tet+ cells, KLH-DNP-induced T-follicular helper (TFH) cells, and Tconv cells vs. differential gene expression.

(A) Pie charts showing the relative distribution of H3K4me3 (top), H3K27me3 (middle), and H3K27Ac (bottom) marks in different gene regions. (B) Jitter plots comparing differences in gene expression between BDC2.5mi/I-A^g7-NP-induced TR1/TR1-like and TFH.1 cells (log2 fold change between the two cell clusters in the single-cell Multiome [scMultiome] dataset) for genes linked to open chromatin regions (OCRs) shared between BDC2.5mi/I-A^g7-NP-induced Tet+ and KLH-DNP-induced TFH cells, as a function of shared or differential H3K4me3 (top), H3K27me3 (middle), and H3K27Ac (bottom) deposition status (false discovery rate [FDR]<0.01). Gene expression was not associated with differential histone deposition in any case (Wilcox test for differential gene expression: adjusted p<0.05): H3K4me3: Pearson’s Chi-square test p=0.53; H3K27me3: Fisher exact test p=1; H3K27Ac: Pearson’s Chi-square test p=0.23. Color depicts the -log10(adjusted p-value) of scRNAseq analysis. Bubble size corresponds to the number of marked regions annotated to specific genes. Most differentially expressed and TR1/TR1-like and TFH-relevant genes included in the 106 genes list from Supplementary file 1 are labeled. KLH, keyhole limpet hemocyanin; NP, nanoparticle.

Figure 4.. Gene methylation status in BDC2.5mi/I-A^g7-NP-induced Tet+, KLH-DNP-induced T-follicular helper (TFH), and Tconv cells.

(A) Venn diagram of differentially methylated regions (DMRs) in BDC2.5mi/I-A^g7-NP-induced Tet+ and KLH-DNP-induced TFH vs. Tconv cells. (B) Venn diagrams of differentially hypomethylated (left) or hypermethylated (right) regions shared by BDC2.5mi/I-A^g7-NP-induced Tet+ cells and KLH-DNP-induced TFH cells as compared to their Tconv counterparts. Gene names found in the 106 TFH/TR1/Treg gene list in Supplementary file 1, with shared methylation status, are indicated. (C) Bar plot showing the numbers of differentially hypo- or hypermethylated regions (separated by gene region: promoter, gene body, or distal intergenic) in BDC2.5mi/I-A^g7-NP-induced Tet+ and KLH-DNP-induced TFH vs. Tconv cells, respectively. DMRs are classified, from left to right, into those only found in BDC2.5mi/I-A^g7-NP-induced Tet+ cells (only Tet), KLH-DNP-induced TFH cells (only TFH), or found in both subsets (shared by tet+ and TFH). (D) Pie-donut chart showing the distribution of DMRs (hyper- and hypomethylated status) in genes associated with open chromatin regions (OCRs) (n=328) shared by TFH.1 and TR1/TR1-like cells (identified via single-cell Multiome [scMultiome]). (E) Jitter plot comparing differential gene expression between TR1/TR1-like and TFH.1 cells (as determined by scMultiome) and differential methylation status, focusing on the genes that remain open at the TR1/TR1-like cell stage, as determined by scMultiome. DMR-associated genes are classified based on their cell pool specificity, i.e., only found in BDC2.5mi/I-A^g7-NP-induced Tet+ cells (only Tet+), KLH-DNP-induced TFH cells (only TFH), or both (shared by Tet+ and TFH). Color depicts the methylation status (hypo- or hypermethylated) of the regions associated with these genes. No significant correlation between methylation status and differential gene expression was found. KLH, keyhole limpet hemocyanin; NP, nanoparticle.

Figure 4—figure supplement 1.. Differences in gene region distribution of differentially methylated regions (DMRs) among BDC2.5mi/I-A^g7-NP-induced Tet+, KLH-DNP-induced T-follicular helper (TFH), and Tconv cells.

Pie charts show the relative distribution of DMRs in BDC2.5mi/I-A^g7-NP-induced Tet+ vs. Tconv (left), KLH-DNP-induced TFH vs. Tconv (middle), and BDC2.5mi/I-A^g7-NP-induced Tet+ vs. KLH-DNP-induced TFH (right). KLH, keyhole limpet hemocyanin; NP, nanoparticle.

Figure 4—figure supplement 2.. Differential DNA methylation at or near the Il2, Il10, and Il21 loci between BDC2.5mi/I-A^g7-NP-induced Tet+ or KLH-DNP-induced T-follicular helper (TFH) and Tconv cells.

Lolliplots comparing the location of differentially methylated CpGs in Il2 (left), Il10 (middle), and Il21 (right) in BDC2.5mi/I-A^g7-NP-induced Tet+ vs. Tconv (top), KLH-DNP-induced TFH vs. Tconv (middle), and BDC2.5mi/I-A^g7-NP-induced Tet+ vs. KLH-DNP-induced TFH cells’ comparisons (bottom). KLH, keyhole limpet hemocyanin; NP, nanoparticle.

Figure 4—figure supplement 3.. Relative contribution of different epigenetic marks to changes in gene expression during the TFH-to-TR1 cell conversion (TFH vs. Tconv).

Heatmap depicting the presence of different epigenetic marks (from left to right: differential methylation, differential open chromatin regions [OCRs], differential H3K27ac deposition, differential H3K4me3 deposition, and differential H3K27me3 deposition) in antigen-induced TFH vs. Tconv cells. Data correspond to differentially expressed genes with shared OCRs between TR1 and TFH cells. Differential epigenetic data is scaled for each technique and when multiple genomic regions are associated to a gene, the average is provided. Genes are ordered by decreasing fold change of differential gene expression. TFH, T-follicular helper; TR1, T-regulatory type 1.

Figure 4—figure supplement 4.. Relative contribution of different epigenetic marks to changes in gene expression during the TFH-to-TR1 cell conversion (Tet+ vs. Tconv).

Heatmap depicting the presence of different epigenetic marks (from left to right: differential methylation, differential open chromatin regions [OCRs], differential H3K27ac deposition, differential H3K4me3 deposition, and differential H3K27me3 deposition) in BDC2.5mi/IA^g7-NP-induced Tet+ vs. Tconv cells. Data correspond to differentially expressed genes with shared OCRs between TR1 and TFH cells. Differential epigenetic data is scaled for each technique and when multiple genomic regions are associated to a gene, the average is provided. Genes are ordered by decreasing fold change of differential gene expression. TFH, T-follicular helper; TR1, T-regulatory type 1.

Figure 4—figure supplement 5.. Relative contribution of different epigenetic marks to changes in gene expression during the TFH-to-TR1 cell conversion (Tet+ vs. TFH).

Heatmap depicting the presence of different epigenetic marks (from left to right: differential methylation, differential open chromatin regions [OCRs], differential H3K27ac deposition, differential H3K4me3 deposition, and differential H3K27me3 deposition) in BDC2.5mi/IA^g7-NP-induced Tet+ vs. KLH-DNP-induced TFH cells. Data correspond to differentially expressed genes with shared OCRs between TR1 and TFH cells. Differential epigenetic data is scaled for each technique and when multiple genomic regions are associated to a gene, the average is provided. Genes are ordered by decreasing fold change of differential gene expression. KLH, keyhole limpet hemocyanin; TFH, T-follicular helper; TR1, T-regulatory type 1; NP, nanoparticle.

Figure 5.. Relative contribution of different epigenetic marks to changes in the expression of differentially expressed transcription factor (TF)-coding genes during the TFH-to-TR1 cell conversion.

Heatmap depicting the presence of different epigenetic marks (from left to right: differential methylation, differential open chromatin regions [OCRs], differential H3K27ac deposition, differential H3K4me3 deposition, and differential H3K27me3 deposition) in BDC2.5mi/I-A^g7-NP-induced Tet+ vs. KLH-DNP-induced TFH cells. Data correspond to differentially expressed TF-coding genes between TR1 and TFH.1 cells as determined by the single-cell Multiome (scMultiome) analyses. Differential epigenetic data is scaled for each technique and when multiple genomic regions are associated to a gene, the average is provided. Genes are arranged from most to least upregulated, followed by least to most downregulated (fold change). No differentially enriched sites for H3K27Ac histone deposition were associated with differential expression of these genes. KLH, keyhole limpet hemocyanin; TFH, T-follicular helper; TR1, T-regulatory type 1; NP, nanoparticle.

Figure 6.. T-regulatory type 1 (TR1) cells inherit hypomethylated active enhancers from their T-follicular helper (TFH) precursors.

(A) Euler’s plot showing total active enhancer sharing by BDC2.5mi/I-A^g7-NP-induced Tet+ cells, KLH-DNP-induced TFH cells, and Tconv cells. Overlapping of active enhancers in BDC2.5mi/I- A^g7-NP-induced Tet+ cells and KLH-DNP-induced TFH cells was significantly higher than overlapping between BDC2.5mi/I-A^g7-NP-induced Tet+ and Tconv cells (Pearson’s Chi-squared test with Yates' continuity correction p<2.2e-16). (B) Jitter plot comparing differential gene expression between TR1/TR1-like and TFH.1 cells (as determined by single-cell Multiome [scMultiome]), focusing on the genes that remain open at the TR1/TR1-like cell stage, and both distribution and number of active enhancers per gene as a function of their cell pool specificity (i.e. only found in BDC2.5mi/I- A^g7-NP-induced Tet+ cells (only Tet+), KLH-DNP-induced TFH cells (only TFH), or both (shared by Tet+ and TFH)). Color and size depict the region type (gene body or intergenic) and the number of active enhancers per gene, respectively. Gene names are displayed for all the genes except when more than 20 dots are in the same region of the plot. No significant correlation between active enhancer distribution and differential gene expression (Wilcox test for differential analysis: adjusted p<0.05) was found using Pearson’s Chi-square test, p=0.76. (C) Figure displays genome tracks for the various readouts examined herein in Tconv CD4+ T-cells (left), KLH-DNP-induced TFH cells (middle), and BDC2.5mi/I-A^g7-NP-induced Tet+ cells (right). Tracks correspond, from top to bottom, to reads for RNAseq (n=4), ATACseq (n=3), ChIPseq (n=1) for H3K4me3, H3K27me3, and H3K4me3 immunoprecipitation, respectively; ChiPseq (n=1) for Tox- 2 (TFH cells), Irf4 (Th17 cells), and cMaf (Th17 cells) deposition, respectively (see main text for references), active enhancers and differential methylation, respectively. Visualization reads were normalized to total sequencing depth per sample using BPM (bins per million) and, where replicates were available, height mean per bin was also computed. ChIPseq data was also normalized for input (non-immunoprecipitated) sequenced reads. Height (y-axis) is equivalent to the normalized number of mapped reads in each region. Active enhancers (ACT ENH) were predicted as overlapping regions of open chromatin region (OCR) (ATACseq) and H3K27Ac deposited peaks (ChIPseq H3K27Ac) and are depicted as absent or present in each region. DIFF METH shows differentially methylated regions (n=3) obtained comparing BDC2.5mi/I-A^g7-NP-induced Tet+ cells and KLH-DNP-induced TFH to Tconv cells. Height corresponds to the relative mean methylation value for each region. KLH, keyhole limpet hemocyanin; NP, nanoparticle.

Figure 7.. Extensive sharing of differentially methylated active enhancers carrying Tox2, IRF4, and/or c-MAF-binding sites in BDC2.5mi/I-A^g7-NP-induced Tet+ and KLH-DNP-induced T-follicular helper (TFH) cells vs. Tconv cells.

(A) Euler’s plots with differentially methylated active enhancers for each comparison: BDC2.5mi/I- A^g7-NP-induced Tet+ cells vs. Tconv (Tet+; red), KLH-DNP-induced TFH vs. Tconv (TFH; orange), and Tconv vs. BDC2.5mi/I-A^g7-NP-induced Tet+ cells (Tconv; green). Pie charts on the left show the methylation status (hypomethylated or hypermethylated) of each subset of active enhancers for the same subset-to-subset comparisons. The BDC2.5mi/I-A^g7-NP-induced Tet+ and KLH-DNP-induced TFH cell subsets share significantly more differentially methylated active enhancers among each other than with Tconv cells (Pearson’s Chi-squared test with Yates' continuity correction; p<2.2e-16). (B) Jitter plot of differentially methylated active enhancers in BDC2.5mi/I- A^g7-NP-induced Tet+ cells (Tet+) and KLH-DNP-induced TFH (TFH) vs. Tconv cells vs. differential gene expression along the TFH.1-TR1 axis as determined by single-cell Multiome (scMultiome). Plot corresponds to genes that remain accessible as TFH.1 cells transdifferentiate into TR1/TR1-like cells. The differentially methylated enhancers linked to these genes were classified as being specific for BDC2.5mi/I-A^g7-NP-induced Tet+ cells (only Tet+), exclusive for KLH-DNP-induced TFH cells (only TFH) or shared by both Tet+ and TFH cells. Gene bubble color depicts -log10 of adjusted p-value from scRNAseq differential analysis, and bubble size depicts the number of regions associated per gene. Gene names are displayed for all the genes except when more than 20 dots are in the same region of the plot. No statistically significant associations between the presence of differentially hypo- or hypermethylated active enhancers and gene expression differences were found (Chi-square test: p=0.32). (C) Overlap of active enhancers with Tox2, IRF4, and c-MAF binding sites. Histogram plot compares the relative proportion of the active enhancers linked to genes that are accessible in both TR1/TR1-like and TFH.1 cells (as defined via scMultiome) and that are upregulated in TR1/TR1-like vs. TFH.1 cells (as determined via scMultiome), which overlap with binding sites for Tox-2 in TFH cells, and IRF4 or c-MAF in Th17 cells. Bars correspond from left to right, to: (1) binding sites shared by both BDC2.5mi/I-A^g7-NP-induced Tet+ and KLH-DNP-induced TFH cells (shared Tet+ & TFH); (2) exclusively found in the BDC2.5mi/I- A^g7-NP-induced Tet+ pool (only Tet+); (3) shared by BDC2.5mi/I-A^g7-NP-induced Tet+, KLH-DNP-induced TFH and Tconv cells (shared Tet+ & TFH & Tconv); and (4) shared by BDC2.5mi/I-A^g7-NP-induced Tet+ and Tconv cells (shared Tet+ & Tconv). KLH, keyhole limpet hemocyanin; TR1, T-regulatory type 1; NP, nanoparticle.