FoxP3 forms a head-to-head dimer in vivo and stabilizes its multimerization on adjacent microsatellites

Authors: Fangwei Leng, Ryan Clark, Wenxiang Zhang, Thibault Viennet, Cuidie Wang, Haribabu Arthanari, Xi Wang, Sun Hur

Affiliations: Howard Hughes Medical Institute and Program in Cellular and Molecular Medicine, Boston Children’s Hospital, Boston, MA 02115, USA; Department of Biological Chemistry and Molecular Pharmacology, Harvard Medical School, Boston, MA 02115, USA; Institute of Immunology, Chinese Institutes for Medical Research, Capital Medical University, Beijing 100069, China; These authors contributed equally; Shanghai Institute of Immunology, Shanghai Jiaotong University, Shanghai 200025, China; Department of Cancer Biology, Dana-Farber Cancer Institute, Boston, MA 02215, USA; Department of Chemistry and iNANO, Aarhus University, Aarhus, Denmark; Lead contact

License: CC BY 4.0 This is an open access article under the CC BY license (https://creativecommons.org/licenses/by/4.0/).

Article links: DOI: 10.1016/j.celrep.2025.116633  |  PubMed: 41313681  |  PMC: PMC12820563

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (2.3 MB)

INTRODUCTION

Transcription factors (TFs) are sequence-specific DNA-binding proteins central to virtually all gene regulatory networks. Of the ~1,600 TFs characterized to date, about half contain a single DNA-binding domain (DBD)^1,2 and have long been thought to adopt a single, conserved DNA-binding mode. However, emerging evidence reveals a more complex landscape. One example is FoxP3, which was recently shown to adopt multiple conformations depending on the DNA sequence.^3–5 The relationships and functional distinctions between these distinct binding modes, however, remain unclear.

FoxP3 belongs to the forkhead family, whose members adopt the conformation of a winged-helix forkhead domain (FKH) to bind DNA.⁶ Although early reports suggested otherwise,^7,8 FoxP3’s forkhead domain also assumes this canonical fold.³ Unlike most forkhead TFs, however, FoxP3 binds poorly to isolated consensus forkhead motifs (FKHMs).^3–5 Instead, it requires composite recognition sites–either an inverted repeat with a 4-nt gap (IR-FKHM)³ or microsatellites harboring simple tandem repeats of T2G, T3G, T4G, or T5G (collectively referred to as TnG repeats).^4,5,9 In both cases, a loop adjacent to FKH, namely Runx1-binding region (RBR, Figure 1A), plays a critical role as it mediates FoxP3-FoxP3 interactions.^3,10 On IR-FKHM, FoxP3 forms a head-to-head (H-H) dimer through the RBR-RBR interactions.³ On tandem repeats of TnG sequence (n = 2–5), FoxP3 forms complex multimers in various stoichiometries and configurations, using both RBR and FKH.^4,5 The multimerization bridges 2–4 DNA segments^4,5 and stabilizes chromatin loops in vivo.^4,11,12 Although FoxP3 is constitutively dimeric through its coiled coil (CC) domain¹³ (Figure 1A), these additional RBR-mediated interactions are essential for cooperative binding to IR-FKHM and TnG repeats.^3–5

A central question is whether these in vitro observed conformations occur in vivo. Genomic analyses of FoxP3-occupied loci––via cleavage under targets and release using nuclease-sequencing (CNR-seq),^12,14 chromatin immunoprecipitation sequencing (ChIP-seq)^14,15 or ChIP-seq with exonuclease treatment (ChIP-exo)⁹––have consistently shown strong enrichment of TnG repeats,^4,5,9 supporting a model where FoxP3 multimerizes to stabilize chromatin loops in Tregs. By contrast, IR-FKHM motifs associated with H-H dimerization were previously reported to show only minimal enrichment within a relatively limited but highly reproducible set of FoxP3-bound regions.³ Functionally disentangling H-H dimerization from multimerization using mutagenesis has also proven difficult. This is because their interaction interfaces overlap, while some residues involved in multimerization lie outside the RBR, and those critical for H-H dimerization, which are mainly within the RBR, are also essential for multimerization.^4,5 Notably, while multimerization on TnG is common to all four members of the FoxP family (FoxP1–4),⁴ H-H dimerization is unique to FoxP3,³ raising the question of whether it reflects an evolutionarily adopted function specific to FoxP3.

We hypothesized that H-H dimerization occurs in vivo, but was undetected because it takes place on DNA sequences beyond IR-FKHM, and/or because previous analyses of in vivo occupancy data were too restrictive. To test these possibilities, we conducted a more systematic in vitro pull-down sequencing together with broader genomic data analyses. These results revealed that FoxP3 recognizes a wider range of motifs as an H-H dimer than previously appreciated, and that these H-H motifs, including IR-FKHM, are indeed enriched within the FoxP3 genomic targets in vivo. These H-H sites often lie immediately adjacent to TnG repeats, where H-H dimerization assists FoxP3 multimerization on short, suboptimal TnG repeats. Together, these findings uncover a previously unrecognized role for H-H dimerization in expanding FoxP3’s architectural functions and provide previously uncharacterized mechanistic insight into its mode of action.

RESULTS

FoxP3 recognizes diverse DNA motifs paired with FKHM through H-H dimerization

To capture the full repertoire of FoxP3’s sequence preferences, we performed an unbiased pull-down assay followed by deep sequencing (PD-seq) using a library of degenerate DNA oligonucleotides. Unlike our previous genome-based PD-seq⁴–which predominantly identified TnG repeats due to their genomic abundance¹⁶ and the large fragment sizes that favor multimerization⁴–we here specifically aimed to isolate the dimeric binding mode. For this reason, we designed oligos containing a randomized 7-nt segment (site 1) followed by a 4-nt spacer and a fixed reverse-complement FKHM (rcFKHM, site 2), all flanked by 27 nt upstream (NBS1) and 35 nt downstream (NBS2) of a non-binding sequence (Figure 1A; Table S1). Note that if site 1 were an FKHM, this design would recreate an IR-FKHM that facilitates FoxP3 H-H dimerization.³ For clarity, all sequences are reported on the forward strand.

Recombinant FoxP3 protein fused with a maltose-binding protein (MBP) tag was used to pull down via amylose resin, as in the previous study.⁴ Since full-length FoxP3 could not be purified due to solubility issues and non-specific proteolytic cleavage, we used a construct lacking the N-terminal disordered region (residues 188–423, FoxP3^ΔN). Both the MBP tag and the truncation of the N-terminal portion have little impact on DNA sequence specificity.^3,4 MBP-FoxP1^ΔN (residues 326–579) construct with the equivalent boundary was used for comparison (Figure 1A).

We first ranked individual sequences based on read counts normalized to input (Figure 1B, see the full sequence list in Table S2). The two biological replicates of FoxP3-PD-seq were highly concordant (Figure S1A), whereas FoxP3 and FoxP1 PD-seq displayed a striking divergence (Figures 1B and S1B). While FoxP3 binding strongly depended on site 1 sequence, FoxP1 exhibited markedly lower reliance on this region. Testing selected sequences from the top-, middle-, and bottom-ranked groups based on the FoxP3 PD-seq further confirmed this conclusion (Figure 1C). These findings suggest that FoxP3 recognizes both site 1 and site 2, whereas FoxP1 engages only the rcFKHM in site 2. Moreover, minimal enrichment of FKHM in site 1 for FoxP1––despite its known ability to bind isolated FKHM––implies that site 1 occupancy is either incompatible with site 2 co-occupancy or not needed for FoxP1.

To identify motif preference, we initially employed conventional de novo motif discovery tools, such as MEME.¹⁷ However, whether applied to the entire 18-nt region encompassing site 1, a 4-nt gap, and rcFKHM or just site 1, MEME identified few discernible sequence features (Figures S1C and S1D). Instead, we systematically analyzed read per million (RPM) for each unique sequence of site 1. Allowing degeneracy at 2 out of 7 bases (2N motifs) markedly simplified the representation of highly enriched sequences, enabling 10 motifs to collectively account for the sequences with top 10% of enrichment values (Figure S1E, see Table S3 for the top 25 motifs). Pull-down assays of representative sequences confirmed that all 10 motifs are comparable to IR-FKHM in FoxP3 binding (Figure S1F). These 10 motifs were classified into two groups: (1) those containing TGTTT or TGTTG, closely resembling FKHM (TGTTTAC), and (2) those containing GCAT flanked by either T or C on each side, akin to a non-canonical FoxM1 motif¹⁸ but previously unrecognized for FoxP3 (Figure 1D).

Notably, all top 10 motifs paired with FKHM displayed a strong orientational bias (Figure 1E), suggesting that FoxP3’s binding to site 1 depends on its binding to site 2 sequence rcFKHM. We hypothesized that this bias reflects FoxP3’s H-H dimerization. Supporting this hypothesis, group 1 (G1) sequences in site 1 align with FKHM, forming sequences similar to IR-FKHM with a gap of 3–4 nt (Figure 1D); that is, motifs #2 and #4 in site 1 recreate IR-FKHM with a 4-nt gap, whereas motifs #5 and #8 generate 3-nt gapped variants (Figure 1D). This is consistent with previous observations that FoxP3 prefers IR-FKHM with a 4 nt-gap but can tolerate a 3-nt gap.³

To test whether FoxP3 also binds group 2 (G2) sequences paired with rcFKHM in an H-H orientation, we employed a FoxP3-NFAT cooperativity assay,³ which exploits NFAT’s ability to dock onto FoxP3’s “tail” only when their DNA sites are correctly phased^7,8 (Figure 1F). In this assay, a 7-nt test sequence placed 3 nt downstream of the NFAT site (GGAAA) will restore cooperativity if FoxP3 engages that sequence in the same register and orientation as FKHM. When we tested GCGCATC (with the conserved core sequence CGAT underscored), we observed NFAT-dependent enhancement of FoxP3-DNA binding exclusively at a 3-nt offset and only in the forward orientation (Figure 1G, left). When the core GCAT was shifted by 1 nt in another G2 sequence TGTGCAT, optimal cooperativity was observed at a 2-nt offset, again only in the forward orientation (Figure 1G, right). These results suggest that FoxP3 binds the NNGCATN motif in the same manner and register as FKHM, and the enriched sequences are consistent with H-H dimerization with a 4-nt gap.

Altogether, these findings suggest that FoxP3 forms an H-H dimer on a broader set of DNA sequences, including TGTTNNN or NNGCATN paired with rcFKHM at a 4 or 3 nt-gap, and that this property is specific to FoxP3 and not shared with FoxP1.

FoxP3 recognizes pairs of group 1 and 2 sequences in H-H configuration

Next, we examined how FoxP3 binds DNA when both sites 1 and 2 were allowed to vary, rather than fixing site 2 to rcFKHM. Specifically, we asked whether FoxP3 could bind the 10 motifs in Figure 1 when paired with each other (without rcFKHM) and whether the orientational bias for the H-H configuration is preserved. To test this, we performed FoxP3 PD-seq using a DNA oligo library based on the construct in Figure 1A, but with degenerate sequences in both sites 1 and 2 (random-random, Figure 2A). Because this expanded library contained 4¹⁴ (268,435,456) possible sequences, which is too large for comprehensive input sequencing, we sequenced only the PD-enriched samples to a depth of ~280 million reads each. For FoxP3, sharp enrichment was seen for fewer than 10,000 unique sequences, which differed vastly from that of FoxP1 PD-seq (Figure 2B), consistent with their difference in DNA specificity.

To examine the orientational bias of FoxP3, we generated heatmaps showing relative read counts of each pair of the 10 motifs in sites 1 and 2 in one of the four orientations (H-H, head-to-tail [H-T], tail-to-head [T-H], or tail-to-tail [T-T]). The comparison showed that H-H pairing is indeed preferred, compared to the other configurations (Figure 2C).

Since the binding positions of FoxP3 in the random-random library were no longer fixed (unlike in the random-rcFKHM library), we also allowed for −2, −1, 0, +1, and +2 bp shifts when scanning for motifs (Figure 2A). This, however, inevitably introduces the influence of both flanking and gap sequences. To correct this, we computed the baseline probability of observing each motif pair under random shifts and normalized the observed motif frequency in the PD samples by these baseline probabilities (see STAR Methods). Even after this positional relaxation, H-H dimerization remained markedly favored over other configurations (Figures 2C and 2D). This bias for H-H orientation was more pronounced in FoxP3 than in FoxP1 (Figures 2C and 2D), in agreement with the random-rcFKHM library result in Figure 1.

Note that random-random PD-seq results were not used for de novo motif discovery due to the ambiguity in the FoxP3 register, which can artificially influence motif enrichment. For example, although the inverted repeat of GCGYGCH (Y indicates T/C, and H indicates A/C/T) appeared enriched in the top 25 sequences from this library (Table S4), this does not represent a bona fide FoxP3 motif because GCGYGCH (or its reverse complement) was not enriched in the random-rcFKHM library (green lines in Figure S1E). Moreover, mutating the spacer or flanking sequences in IR-GCGTGCA abolished FoxP3 binding, whereas similar mutations had minimal impact on IR-FKHM binding (Figure S1G). This suggests that the apparent enrichment of IR-GCGYGCH reflects FoxP3 recognition of flanking or spacer elements rather than the core motif itself. Consistent with this, closer inspection revealed that IR-GCGYGCH within our random-random construct contains two inverted pairs of high-ranking motifs from Figure 1 in ±2 shifted positions (Figure S1H), likely explaining its strong enrichment. Given these complexities and the high risk of over-interpretation, we relied primarily on the random-rcFKHM library for de novo motif discovery, while using the random-random library mainly to validate pairing orientations.

Altogether, these data confirmed H-H dimerization of FoxP3 using pairs of both G1 and G2 sequences and the lack of such dimerization for FoxP1.

FoxP3 recognizes both H-H motifs and TnG repeats in Tregs

Using the identified FoxP3 motifs, we systematically assessed whether H-H motifs are enriched within FoxP3-occupied genomic regions in vivo. We first compared the frequency of individual motifs in FoxP3-bound regions–25,724 peaks defined using four individual FoxP3 Treg ChIP-seq data from two groups^14,15–with their frequency in open chromatin regions (OCRs) from Treg ATAC-seq¹⁹ (Figure 3A). This definition of FoxP3-bound regions is intentionally more relaxed than in the previous study,³ as discussed later. For FoxP3-bound regions, we also examined CUT&RUN-seq (CNR-seq, 18,552 peaks) datasets from Tregs^12,14 and recent ChIP-exo data⁹ (4,307 peaks), which provide precise FoxP3 footprints, though derived from ectopically expressed FoxP3 in mouse embryonic stem cells (mESCs) (Figure S2A). For significantly enriched motif pairs (p < 0.05; black circles in Figure 3A), fold enrichment values were visualized in heatmaps (Figures 3B and S2A). Motif frequencies were calculated only on the positive reference strand; analyzing the negative strand would convert T-H into H-T and invert the rows/columns for H-H and T-T, explaining the similarity between H-T and T-H maps and the diagonal symmetry observed in H-H and T-T maps.

The results reveal a strong enrichment of H-H pairs of G1-G1 motifs in FoxP3-bound regions across all three in vivo datasets (ChIP-seq, CNR-seq, and ChIP-exo; Figures 3B and S2A). The H-H motif pairs were found in 2,375 sites (2,270 unique peaks) within the total 25,724 FoxP3-ChIP peaks (184.96 per Mb), compared to 10,140 occurrences in 502,437 ATAC-peaks (80.34 per Mb). Among these 2,375 FoxP3-occupied H-H sites, the canonical IR-FKHM (inverted repeat of TGTTTNN, Figure 3B) appeared in 568 sites (24%), while the majority (76%) contained non-canonical motifs (i.e., non-IR-FKHM). FoxP3 preferentially bound H-H pairs of relaxed FKHM-like sequences (G1-G1) over G1-G2 or G2-G1, with G2-G2 pairings showing no significant enrichment (Figures 3B and 3C). This preference mirrored intrinsic affinity patterns observed in vitro using purified FoxP3 and DNA by electromobility shift assay (EMSA) (Figure 3D). Altogether, while IR-FKHM represents the most favored motif for FoxP3 in vivo and in vitro, non-canonical motifs collectively make substantial contributions to FoxP3’s H-H dimerization in cells.

In contrast to FoxP3, FoxP1 showed little sign of bias for H-H pairs both in vivo and in vitro. In CD8 T cells (chosen because they lack FoxP3, which can influence its genomic targeting), FoxP1 ChIP-seq data²⁰ did not show any H-H motif enrichment (Figures 3B–S2B), consistent with our in vitro finding that FoxP1 does not display intrinsic preference for H-H motifs (Figures 1 and 2).

In addition to the strong H-H enrichment, the second most prominent feature in the FoxP3 ChIP-seq heatmaps was the presence of H-T and T-H pairs of G1-G1 motifs (Figure 3B). We noted that TnG repeats––long tandem arrays of TnG or AnC sequences that are occupied by FoxP3 in vivo and in vitro^4,5,9––were included in the H-T/T-H motif pair analysis (Figure 3E). Specifically, 86.31% and 84.41% of significant H-T and T-H G1-G1 pairs, respectively, overlapped with TnG repeats within FoxP3 ChIP peaks (Figure 3F). By contrast, H-T/T-H arrangements of G2-G2 pairs were less enriched, likely because tandem repeats of GCAT (G2) and GCGT (G2-like), while capable of supporting FoxP3 multimerization (Figure S2C) and DNA bridging (Figures S2D and S2E), were less efficient than G1-G1 tandem repeats.

Note that these H-T/T-H pairs of G1-G1 motifs were also enriched in FoxP3 ChIP-exo and CNR-seq (Figure S2A), but were not enriched in oligo PD-seq (Figure 2B), presumably because the random-random oligo library only allowed 7 + 7 bp variation with fixed flanking and spacer sequences, which precluded TnG repeats. H-T/T-H G1-G1 pairs were also enriched in FoxP1 ChIP peaks (Figure 3B), where they too map primarily to TnG repeats (Figure 3F), consistent with the common ability of FoxP TFs to multimerize on these microsatellites.⁴

Overall, the paired motif analyses using genomic occupancy data indicate that FoxP3 employs both G1 and G2 motifs to form H-H dimers in vivo. In addition, FoxP3 forms H-T/T-H multimers on G1 motifs (TnG repeats). Although tandem repeats of G2-like sequences can also facilitate multimerization, they do so with markedly lower efficiency and enrichment. By contrast, T-T motifs were not enriched, suggesting that this configuration is not favored for FoxP3 dimerization or multimerization.

FoxP3-bound H-H motifs are enriched for a 4-nt gap and are preferentially associated with FoxP3-regulated genes

The observed enrichment of H-H motifs contrasts with a previous report³ that found no such enrichment. Given this contradiction, we further investigated the significance of the observed enrichment of H-H motifs. We first compared two analysis methodologies that led to the divergent outcomes. That earlier report³ was restricted to 5,047 ChIP peaks that are highly reproducible in both peak width and locations (defined by 50% reciprocal overlap across four ChIP-seq datasets^14,15). In contrast, our analysis in Figure 3 used peaks directly called from the same four datasets with a threshold of p < 0.001, yielding 25,724 peaks (see STAR Methods for details). The H-H motif-containing loci showed concordant ChIP signal, regardless of the datasets (Figure S3A).

We first asked whether the broader peak definition might have artificially contributed to H-H enrichment; we compared gap-size distributions of significant motif pairs in FoxP3-bound loci (25,724 ChIP peaks) vs. Treg OCRs (502K ATAC peaks). If enrichment were an artifact of peak broadening, no particular gap size should be favored. Indeed, OCRs showed no preference in gap size (Figure 4A), whereas FoxP3-bound regions were strongly enriched for a 4-nt gap across all significant motif pairs (Figure 4A)–the optimal spacing for H-H dimerization.³ Similar enrichment was observed in FoxP3-bound loci defined by 18,552 CNR peaks (Figure S3B). These results indicate that H-H motif enrichment reflects genuine FoxP3 binding properties, not artifacts of peak definition.

Why, then, was H-H motif enrichment missed in the earlier analysis³? Although most of the earlier ChIP peaks (4,738 out of 5,047) overlapped with our current 25,724 peak set, the smaller dataset may have lacked sufficient statistical power. Additionally, the narrow focus on the strict IR-FKHM rather than broader H-H motifs could have contributed further. Consistent with both of these ideas, re-analysis of the original 5,047 peaks with our expanded H-H definitions revealed that, while a modest 4-nt gap preference was observed for some non-canonical motifs (e.g., TGTTKNN_GATGCNN and NNGCATY_NNAAACA) (Figure S3C), the baseline distribution was uneven, diminishing statistical robustness. Together, these findings suggest that restrictive motif definitions combined with limited peak sets may have contributed to the lack of detectable H-H motif enrichment.³

To evaluate the functional significance of H-H motifs and compare them with that of TnG repeats, we analyzed their genomic distribution. Within the 25,724 FoxP3 ChIP peaks, both H-H motifs and TnG repeats were broadly distributed across various genomic features, from promoter-proximal regions to gene-distal sites (Figure 4B). Histone mark analyses further supported this widespread distribution (Figure S3D). To link motif binding to gene regulation, we focused on genes affected by acute FoxP3 degradation,²¹ as FoxP3 knockout models are known to alter hundreds of genes indirectly.^14,22 Strikingly, genes directly regulated by FoxP3 contained significantly higher densities of both H-H motifs and TnG repeats than unaffected genes (Figure 4C).

We further examined the fraction of FoxP3-regulated genes that harbor H-H motifs or TnG repeats within 2 kb. Based on differential gene expression analysis, genes most significantly affected by FoxP3 depletion (lower p values) were far more likely to contain these motifs, whereas this enrichment declined with the significance of the FoxP3′ effect (Figure 4D). Importantly, most motif-associated, FoxP3-regulated genes were upregulated upon FoxP3 loss, consistent with its predominant role as a suppressor.²¹ Many of these FoxP3-regulated genes encode factors involved in T cell biology (e.g., Rora, Pde7a, and Dpp4) (Figure 4E). Thus, these results support a model in which FoxP3 binds H-H motifs and TnG repeats to fine-tune the Treg transcriptional program.

H-H dimerization seeds FoxP3 multimerization on adjacent TnG repeats

Our analyses above suggest that although TnG repeats are more frequent than H-H motifs, the two often occur in close proximity and share similar genomic features, raising the question of how FoxP3 dimerization on H-H motifs relates to its multimerization on TnG repeats. A systematic analysis of FoxP3 ChIP peaks revealed that 68.9% of H-H–containing peaks (1,563 of 2,270) also harbored TnG repeats (Figure 5A). By contrast, only 26.9% of H-H-containing open chromatin regions (2,697 of 10,027 ATAC peaks) overlapped with TnG repeats (Figure 5A).

To better define the spatial relationship between H-H motifs and TnG repeats, we measured the distances between these motifs (border-to-border distance), designating positive values when a TnG repeat occurred downstream and negative values when upstream (Figure 5B, top). Overlaps were counted as a 0-nt gap. This analysis revealed a pronounced peak at 0 nt in FoxP3-bound loci (Figure 5B, bottom; examples in 5C), with 29.6% of H-H sites (702 of 2,375) lying immediately adjacent to or overlapping with a TnG repeat, compared with only 9.4% in FoxP3-free OCRs. Strikingly, these abutting motifs displayed a strong orientational bias, such that the proximal half-site of the H-H motif could seamlessly extend into the neighboring TnG repeat (“aligned”; Figures 5B and 5C).

Prompted by this close positioning, we asked whether H-H proximity enhances FoxP3 occupancy on TnG repeats. We compared ChIP signals at repeats with adjacent H-H motifs (0-nt gap) to those located >1 kb from any H-H site. We stratified TnG repeats into short (≤20 bp) and long (≥40 bp) groups (Figure 5D), because longer repeats inherently bind FoxP3 more efficiently.⁴ Consistent with this, FoxP3 occupancy increased with TnG repeat length, when there were no H-H motifs nearby (Figure 5E). H-H adjacency boosted FoxP3 occupancy for both short and long TnG repeats, but the effect was markedly greater for short repeats (Figure 5E). As a result, TnG repeats adjacent to H-H motifs exhibited little dependence on the repeat length, unlike TnG repeats without H-H. These findings indicate that H-H dimerization seeds or stabilizes FoxP3 multimerization on otherwise suboptimal, short TnG repeats.

In support of this notion, 5 repeats of TTTG [(TTTG)₅]–which alone is too short to drive robust FoxP3 multimerization,⁴ even with a single adjacent FKHM (lanes 13–14, Figure 5F)–greatly facilitated multimer assembly when an IR-FKHM was placed immediately adjacent (lanes 5–6, Figure 5F). Moreover, this positive effect was seen only when IR-FKHM was in the “aligned” orientation as (TTTG)₅ and with a 0- to −5-nt gap, but not in the “divergent” configuration (lanes 3–4) or with a 10-nt gap (lanes 5–12). Likewise, DNA-bridging on (TTTG)₅ was substantially enhanced by a neighboring H-H motif only in the aligned orientation with a 0- to −5-nt gap (Figure 5G).

Together, these findings demonstrate that H-H dimerization expands FoxP3’s TnG sequence specificity, enabling efficient multimerization and DNA bridging even on otherwise suboptimal, short TnG repeats.

RBR is the key determinant for H-H dimerization of FoxP3

Our finding that only FoxP3, and not FoxP1, can form H-H dimer, while both can form H-T multimers, raises the question regarding the molecular basis of this specificity. It is particularly intriguing given that FoxP3 and FoxP1 share conserved domain architectures and similar FKH DBDs (Figures 1A and 6A). We first investigated whether a small number of divergent residues within the FKH DBD (e.g., R358, N376, H377, P378, S400, E401, and R414; Figures S4A and S4B) underlie FoxP3’s unique H-H dimerization. However, mutating these residues minimally affected FoxP3’s preference for IR-FKHM over single isolated FKHM (sFKHM) (Figure S4C), suggesting that they are not the primary determinants of FoxP3’s capacity for H-H dimer formation.

Previous studies showed that the hydrophobic loop (namely RBR) immediately preceding FKH DBD is important for both H-H dimerization and H-T multimerization of FoxP3.^3–5 Notably, the sequence of FoxP3’s RBR differs considerably from that of FoxP1, FoxP2, and FoxP4 (Figure 6A). We therefore investigated whether RBR alone dictates the divergent ability to form the H-H dimer by FoxP3 and FoxP1. Note that the removal of RBR leads to misfolding of FKH of all four FoxP TFs, from the winged-helix conformation to domain-swap dimer,³ precluding the RBR deletion test. We instead swapped the RBR between FoxP3 and FoxP1 and examined the impact on their DNA binding behaviors (Figure 6B). We used RBR-FKH constructs (FoxP^RBR–FKH)–the minimal constructs required for proper folding and maintenance of the full-length protein’s DNA specificity.³ Unlike the FoxP^ΔN constructs, which form a constitutive dimer through a CC domain,¹³ FoxP^RBR–FKH is monomeric,³ allowing us to distinguish monomeric from dimeric binding based on distinct EMSA band patterns. Since monomeric binding is generally weak,³ NFAT was used to stabilize FoxP3-DNA interactions, enabling reliable capture of both monomeric and dimeric complexes.

Consistent with a previous report,³ FoxP3^RBR–FKH bound IR-FKHM as a dimer, but sFKHM as a monomer (Figures 6C and 6D). In contrast, FoxP1^RBR–FKH bound both IR-FKHM and sFKHM as a monomer (Figures 6C and 6D), even when the protein was in excess over DNA. This indicates that FoxP1 is inherently incompatible with double occupancy on IR-FKHM. Remarkably, swapping the RBRs interchanged their binding modes. That is, a chimera with FoxP1^RBR and FoxP3^FKH bound DNA in the same manner as FoxP1, not FoxP3 (Figure 6C). Conversely, a chimera comprising the FoxP3^RBR and FoxP1^FKH bound DNA in the same manner as FoxP3, not FoxP1 (Figure 6D). Furthermore, a similar shift in binding mode was observed upon exchanging RBR loops between FoxP3 and either FoxP2 or FoxP4 (Figures 6E and 6F). When swapping the FoxP3 RBR loop with that of FoxP3 from different species, we preserved the dimeric binding modes (Figures S4D and S4E), although with a modest reduction in DNA binding affinity for some RBR species. Collectively, these findings indicate that RBR is a modular element and the sole determinant for H-H dimerization in FoxP3.

To further delineate the specific region within the RBR loop responsible for FoxP3’s distinct H-H dimerization property, we divided the loop into three segments (A, B, and C; Figure 6A), each comprising 20 amino acids. Through targeted sequence-swapping, we identified fragment C (amino acids 317–336; RBR-C) as critical for FoxP3’s H-H dimerization. Substituting FoxP3’s RBR-C with the corresponding region from FoxP1 abolished FoxP3’s ability to form dimers on the IR-FKHM element, whereas introducing FoxP3’s RBR-C into FoxP1 conferred dimerization capacity (Figures 6C and 6D).

A comparison of RBR-C sequences revealed that FoxP3 contains a higher density of aromatic residues than that of FoxP1/2/4: six out of 20 residues in FoxP3’s RBR-C are aromatic, compared with only two in FoxP1/2/4 (Figure 6A). Mutating these aromatic residues (Y330A, F331A, and Y333A) significantly reduced FoxP3’s ability to bind IR-FKHM as a dimer (Figure 6G). These aromatic residues also contribute to FoxP3’s multimerization on TnG repeat DNA and FoxP3’s transcriptional activity in CD4 T cells,⁵ underscoring the importance and the dual role of these aromatic RBRs in mediating both H-H dimerization and multimerization.

In addition to differences in aromatic content, secondary structure predictions indicated an α-helix in the RBR-C region of FoxP1, whereas FoxP3 was predicted to contain only a short, four-amino acid helical stretch (Figure S4F). NMR chemical shifts confirmed a strong α-helical propensity for FoxP1 residues 477–488 within RBR-C (Figure S4G). In case of FoxP3, NMR signals were weaker and displayed line broadening––indicative of conformational exchange occurring on the NMR timescale––which precluded direct comparisons between FoxP1 and FoxP3 RBR conformations. Nonetheless, the previously reported structures of FoxP3 H-H dimer³ and H-T multimers^4,5 revealed both α-helical and non-helical, disordered RBC-C conformations, consistent with its intrinsically lower helical propensity.

Taken together, these results pinpoint RBR-C as the primary driver of H-H dimerization in FoxP3, possibly due to differences in hydrophobicity and helical propensity.

DISCUSSION

Previous studies have shown that FoxP3 binds DNA in two distinct modes: a multimeric conformation when binding long (>30 bp) TnG repeats (n = 2–5)^4,5 and a H-H dimeric conformation on shorter motifs such as IR-FKHM.³ FoxP3 multimerization occurs at numerous sites, facilitating chromatin looping through DNA bridging and global transcriptional regulation^4,5 (Figure 7A). In contrast, in vivo evidence and the functions of FoxP3’s H-H dimerization have remained elusive. Using a degenerate oligonucleotide pull-down approach, we identified a broader repertoire of sequences that support H-H binding (Figure 7B). Both FKHM-like (group 1) and GCAT-containing motifs (group 2) can support FoxP3 binding in H-H orientation with a preferred 4-nt gap. These motifs were enriched not only in vitro but also at FoxP3-occupied sites in vivo, as revealed by ChIP-seq, ChIP-exo, and CNR-seq. Their enrichment was overlooked in our earlier analysis³ largely due to the use of a limited number of ChIP peaks, which reduced statistical power. Importantly, both H-H motifs and TnG repeats were more frequently associated with FoxP3-regulated genes than FoxP3-independent genes in Tregs, underscoring their functional relevance.

What, then, are the molecular functions of FoxP3’s two binding modes? Our analyses showed that TnG repeats are more common (~50% of ChIP peaks) compared to the H-H motifs (~10%). Despite their low frequency, H-H motifs often reside adjacent to TnG repeats, where they can “seed” or stabilize FoxP3 multimerization on the neighboring sequences, particularly on shorter or suboptimal repeats that would otherwise fail to form stable multimers (Figure 7B). Once established, FoxP3 multimerization promotes DNA bridging and the stabilization of chromatin loops.^4,5 Because FoxP3 primarily binds pre-accessible regions marked by H3K27ac, individual looping events may not necessarily drive large transcriptional shifts. Instead, its binding and bridging at thousands of loci across the genome may have subtle, yet global modulation of gene expression dynamics. This model is in keeping with the finding that FoxP3 typically exerts only modest effects (<2- to 3-fold) across hundreds of genes, with its regulatory impact highly context-dependent.²¹ We therefore propose that H-H motifs serve a complementary role to TnG repeats, extending FoxP3’s architectural influence to a wider genomic repertoire and thereby broadening its ability to maintain the Treg transcriptional landscape.

Limitations of the study

Our study demonstrates that FoxP3 forms H-H dimers on genomic DNA in vivo and can promote multimerization on adjacent TnG repeats. However, whether this represents the sole function of H-H dimerization remains unknown. Future work should examine differential cofactor recruitment and locus-specific transcriptional effects of H-H dimers versus TnG multimers.

STAR★METHODS

RESOURCE AVAILABILITY

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Sun Hur (sun.hur@crystal.harvard.edu).

Materials availability

All plasmids generated in this study are available from the lead contact with a completed materials transfer agreement.

Data and code availability

• FoxP3 and FoxP1 pull-down-seq data have been deposited to the GEO database with the accession code of GEO: GSE294472. The NMR data have been deposited in the Biological Magnetic Resonance Bank (BMRB) with accession code 53015.
• The custom codes used in this manuscript are deposited at Zenodo (https://doi.org/10.5281/zenodo.15306890).
• Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS

All data analyzed in this study were obtained from publicly available databases and previously published studies. No new human or animal samples were collected. All in vitro experiments were performed using commercially available reagents and kits according to the manufacturers’ instructions. No cell lines or animal models were used in this study.

METHOD DETAILS

Material preparation

Plasmids

For bacterial expression, the genes encoding mouse FoxP3^ΔN (residues 188–423) and mouse FoxP1 ^ΔN (residues 342–579) were inserted into pET50b between Xmal and HindIII sites, and into a pMAL-c2 vector between BamHI and Xbal sites respectively. The pMAL-c2 vector was modified to contain an HRV 3C protease target sequence between EcoRI and BamHI sites. Other constructs including mouse FoxP3^RBR-forkhead (residues 284–423), mouse FoxP1^RBR-forkhead (residues 427–423), mouse FoxP2^RBR-forkhead (residues 447–588), mouse FoxP4^RBR-forkhead (residues 418–557), FoxP3^RBR-forkhead from different species, and chimeric constructs were generated by inserting genes into pET50b between Xmal and HindIII sites. All FoxP3 mutations including Y330A/F331A/Y333A, N376R/H377N/P378A, S400N/E401V and R414Q were generated by site-directed mutagenesis using Phusion High Fidelity (New England Biolabs) DNA polymerases.

DNA oligos

Single-stranded DNA (ssDNA) oligonucleotides were synthesized by IDT (Integrated DNA Technologies). To generate the degenerate DNA library for FoxP3 and FoxP1 PD-seq, 30 μL of ssDNA (200 μM) was annealed with 60 μL of a complementary primer (Klenow-1, 200 μM) in 10 μL 10× NEB Buffer 2 (final volume: 100 μL). Annealing was performed by heating to 95°C for 5 min followed by slow cooling to room temperature. The reaction was supplemented with 10 μL of 10 mM dNTP mix and 2 μL of Klenow fragment (DNA Polymerase I, large fragment; NEB), and incubated at 25°C for 30–60 min to synthesize the complementary strand. The resulting dsDNA was purified using the QIAquick Nucleotide Removal Kit (Qiagen). DNA purity and integrity assessed by TBE gel electrophoresis and Agilent TapeStation analysis. To prepare biotin-labeled dsDNA for DNA bridging assays, long forward ssDNA was mixed with a short biotinylated complementary ssDNA at equimolar concentrations. After annealing, the Klenow fragment was used to extend the complementary strand, generating full-length biotin-labeled dsDNA. Double-stranded DNA (dsDNA) oligos for EMSA assay, pulldown assay and DNA bridging assay were annealed from single-stranded, complementary oligos. Sequence of all the DNA oligos used are shown in Table S1.

Protein expression and purification

All recombinant proteins in this paper were expressed in BL21(DE3) at 18°C for 16–20 h following induction with 0.2 mM IPTG. Cells were lysed by high-pressure homogenization using an Emulsiflex C3 (Avestin). All proteins are from the Mus. musculus sequence, unless mentioned otherwise. His₆-MBP-tagged FoxP3^ΔN (residues 188–423) and FoxP1 ^ΔN (residues 342–579) used for PD-seq and DNA PD were purified using Ni-NTA affinity column and Superdex 200 Increase 10/300 (GE Healthcare) SEC column in 20 mM Tris-HCl pH 7.5, 500 mM NaCl, 2 mM DTT. FoxP3^ΔN (residues 188–423) was expressed as a fusion protein with an N-terminal His₆-NusA tag. After purification using Ni-NTA agarose, the protein was treated with HRV3C protease to cleave the His₆-NusA-tag and were further purified by a series of chromatography purification using HiTrap Heparin (GE Healthcare), Hitrip SP (GE Healthcare) and Superdex 200 Increase 10/300 (GE Healthcare) columns. The final size-exclusion chromatography (SEC) was done in 20 mM Tris-HCl pH 7.5, 500 mM NaCl, 2 mM DTT. NFAT1 protein (residues 394–680) was also expressed as a fusion protein with an N-terminal His₆-NusA tag. After purification using Ni-NTA agarose, the His₆-NusA-tag was removed using the HRV3C protease and was further purified by SEC on Superdex 75 Increase 10/300 (GE Healthcare) column in 20 mM Tris-HCl pH 7.5, 500 mM NaCl, 5% Glycerol, 2 mM DTT. His₆-NusA tagged proteins are all purified using Ni-NTA affinity column and Superdex 200 Increase 10/300 (GE Healthcare) SEC column in 20 mM Tris-HCl pH 7.5, 500 mM NaCl, 2 mM DTT. All RBR-forkhead proteins (including all chimeric proteins) were expressed as a fusion protein with an N-terminal His₆-NusA tag. After purification using Ni-NTA agarose, the protein was treated with HRV3C protease to cleave the His₆-NusA-tag and were further purified by Superdex 75 Increase 10/300 (GE Healthcare) column in 20 mM Tris-HCl pH 7.5, 500 mM NaCl, 2 mM DTT.

DNA pull-down assay using MBP tagged proteins

0.4 μM purified MBP-mFoxP3^ΔN or MBP-mFoxP1^ΔN protein was incubated with 0.1 μM DNA in the incubation buffer (20 mM Tris-HCl pH 7.5, 100 mM NaCl, 1.5 mM MgCl₂) for 20 min. The protein-DNA mixture was then incubated with 25μL Amylose Resin (New England Biolabs) for 30 min with rotation at room temperature. The bound DNA was recovered using proteinase K (New England Biolabs), purified using QIAquick Nucleotide Removal kit (QIAGEN) and analyzed on 10% Novex TBE Gels (Invitrogen). DNA was visualized by SYBR Gold staining. For NFAT–FoxP3 cooperativity assay, 0.2 μM NFAT protein (residues 394–680) was included along with 0.4 μM MBP-mFoxP3^ΔN and 0.1 μM DNA. All other steps were performed as described above.

DNA bridging assay

Biotin-DNA (bait, 0.1μM) was incubated with Streptavidin Agarose (25μL, Thermo Fisher) in buffer B (20mM Tris-HCl pH 7.5, 100mM NaCl, 1.5mM MgCl₂, 5mM DTT) for 30 min by rotating the mixture at room temperature. Agarose beads were washed three times with buffer B and incubated with non-biotinylated DNA (prey, 0.1 μM) and 0.4 μM purified mFoxP3^ΔN protein. After incubation for 30 min with rotation, bead-bound DNA was recovered using proteinase K (New England Biolabs), purified by QIAquick Nucleotide Removal kit (QIAGEN) and analyzed on 10% Novex TBE Gels (Invitrogen). DNA was visualized by Sybr Gold staining.

Electrophoretic mobility shift assay (EMSA)

DNA (0.2μM) was mixed with the indicated amount of FoxP3 and NFAT in buffer (20mM HEPES pH 7.5, 150mM NaCl, 1.5mM MgCl₂ and 2mM DTT), incubated for 10 min at 4°C and analyzed on 3–12% gradient Bis-Tris native gels (Life Technologies) at 4°C. After staining with SYBR Gold stain (Life Technologies), Sybr Gold fluorescence was recorded using iBright FL1000 (Invitrogen) and analyzed with iBright Analysis Software.

Negative stain EM

FoxP3^ΔN (0.4μM) was incubated with DNA (0.05μM) in buffer (20mM HEPES pH 7.5, 150mM NaCl, 1.5mM MgCl₂ and 2mM DTT) at RT for 10 min. The samples were diluted 10-fold with buffer, immediately adsorbed to freshly glow-discharged carbon-coated grids (Ted Pella) and stained with 0.75% uranyl formate as described before.³⁸sImages were collected using a JEM-1400 transmission electron microscope (JEOL) at 50,000× magnification.

NMR experiments

NMR experiments were measured on a 700 mM sample of ¹³C-¹⁵N murine FoxP1 RBR-FKH (aa 433–579) at 25°C. The buffer was 20 mM MES pH = 6.0, NaCl 150 mM, TCEP 2 mM, ²H₂O 10% v/v. Spectra were collected on a Bruker spectrometer operating at 700 MHz ¹H Larmor frequency, equipped with a triple-channel ¹H, ¹³C, ¹⁵N cryogenically cooled probe and z-shielded gradients. All data were processed using NmrPipe³⁹ and analyzed using CCPNmr Analysis.⁴⁰ Backbone ¹H, ¹⁵N and ¹³C resonances were assigned using a standard set of triple resonance experiments: HSQC, HNCA, HN(CO)CA, HNCO, HN(CA)CO, HN(CA)CB, CBCA(CO)NH. 3D experiments were recorded using non-uniform sampling (NUS), where Poisson-Gap sampling was used to select 10% of the Nyquist grid⁴¹ and the NUS data was reconstructed using the hmsIST protocol.⁴² 80% of expected resonances could be assigned, including 98% of RBR resonances. H_N, N, CO, Cα and Cβ chemical shifts were used to derive dihedral angle values and secondary structure propensities using the software TALOS-N.⁴³ Notably, two conformations of the RBR-C regions were assigned: one with high helical propensity and one with highly disordered character. ¹H-¹⁵N HSQC spectra were compared at 700 and 50 mM and it was found that the disordered conformation does not exist significantly at lower concentrations. Chemical shift assignments have been deposited in the BMRB under accession number 53015.

FoxP3/FoxP1 PD-seq

Binding reactions were carried out by incubating 0.4 μM purified MBP-FoxP3^ΔN or MBP-FoxP1^ΔN protein with 0.1 μM dsDNA library in binding buffer (20 mM Tris-HCl pH 7.5, 100 mM NaCl, 1.5 mM MgCl₂) for 30 min at room temperature. DNA-protein complexes were captured using amylose resin (New England Biolabs), followed by three washes with binding buffer to remove non-specifically bound DNA. Elution of bound DNA was performed by using proteinase K (NEB) and DNA cleanup using the QIAquick Nucleotide Removal Kit (Qiagen). Eluted DNA was amplified using PCR with barcoded primers compatible with Illumina sequencing (primer sequences listed in Table S1), libraries were pooled and subjected to deep sequencing.

Random-rcFKHM library PD-seq analysis

PD library and input samples were sequenced on an Illumina NovaSeq 6000 platform, generating paired-end 150 bp reads to a depth of 20 million reads each. Raw reads were merged with BBMerge²⁴ (v 38.84), combining R1 and R2 reads into single reads at average merge rates of 98.67% (input), 99.47% (FoxP1 PD), and 99.16% (FoxP3 PD). Biopython²⁵ (v1.85) was used to extract trimmed sequences matching the structure NNNNNNN–TCGA–GTAAACA (site 1–gap–rcFKHM). Trimmed sequences were combined across replicates and were organized into a Pandas²⁶ (v 2.2.3) DataFrame, with each row corresponding to a unique 7-nucleotide sequence (excluding any sequences containing ‘N’), yielding exactly 4⁷ sequences. FoxP3/Input fold enrichment values, calculated from the combined replicates, were used to sort sequences for consistent ordering across all plots.

To identify enriched sequence patterns, we identified 7-nucleotide motifs containing two variable positions, denoted as ‘N’. For each of the 4⁷ selected sequences, we systematically generated all possible combinations of two variable positions, resulting in 21 degenerate “2N motifs” per sequence. Duplicate motifs were then removed across all sequences to avoid redundancy. For each resulting 2N motif, we enumerated all 16 possible nucleotide combinations at the two ‘N’ positions (4 possibilities at each position). The counts for these 16 fully specified sequences were summed to calculate the total count for the 2N motif. Enrichment ratios (FoxP3/Input) were computed for each 2N motif by comparing their aggregate counts in the FoxP3 and input datasets. Plots related to PD-seq motif counts and enrichment were visualized using Matplotlib²⁹ (v 3.9.4).

Random-random library PD-seq analysis

FoxP3 PD and FoxP1 PD samples were sequenced on an Illumina NextSeq 2000 platform, generating single-end 150 bp reads to a depth of 280 million reads each. Reads with Phred quality score below 20 in the variable regions were discarded, resulting in ~200 million reads each. Biopython²⁵ (v1.85) was employed to extract trimmed sequences matching the structure: NNNNNNN-TCGA-NNNNNNN (site 1–gap–site 2) and were organized into a Pandas²⁶ (v 2.2.3) DataFrame, with each row corresponding to a unique 18-nucleotide sequence, resulting in 104,367,238 unique FoxP1 sequences and 91,730,585 unique FoxP3 sequences. The DataFrame included columns for the normalized PD-seq counts (in RPM) of each sequence in the FoxP1 and FoxP3 samples.

Using the top 10 sequences from the random-rcFKHM PD-seq (Figure 1E) as a reference, we generated all pairwise motif combinations, resulting in 100 combinations (10X10). For each pair, we examined four orientations: Head-to-Head (H-H), Tail-to-Tail (T-T), Head-to-Tail (H-T) and Tail-to-Head (T-H). This resulted in 400 unique 18-nt long paired motifs, each with a fixed 4-nucleotide spacer “TCGA” between them.

Since the binding positions of FoxP3 in the random+random library were no longer fixed (unlike in the FKHM+random library), we scanned for motifs with −2, −1, 0, +1, and +2 nt offsets around the central gap. Each 18-nt sequence was transformed into five shifted variants that simulate 0, ±1, and ±2 base pair offsets around the central spacer (see Figure 2A). Each transformed sequence in FoxP3 and FoxP1 was checked against the precomputed dictionary containing all 400 paired motifs. If any shifted version matched a motif pair, the original read was counted as supporting that motif interaction. Since fixed sequences in the 4-nt gap and flanking regions introduced bias when motifs were counted during shifting, we computed the baseline likelihood of observing each motif pair under random offsets. We then normalized the observed motif frequency in PD samples by dividing by these baseline probabilities. Normalized frequencies were visualized as heatmaps and bar graphs using Matplotlib²⁷ (v 3.9.4).

Motif analysis using in vivo dataset (ChIP-seq, CNR-seq, ATAC-seq, and ChIP-exo-seq)

As in our previous work,⁴ FoxP3 CNR-seq^12,14 and FoxP3 Treg ChIP-seq^15,44 datasets were trimmed by trimmomatic (v 0.3.6),³⁴ aligned to mm10 using Bowtie2²⁹ (v2.3.4.3) and sorted by read position using SAMtools³⁰ (v1.9). Peaks were called using MACS3³¹ (v3.0.3) with a q-value cutoff of 0.001. Treg ATAC-seq peaks were downloaded from Immgen project,²⁸ and merged using Bedtools³² (v 2.31.0), which were referred to as open chromatin regions (OCRs). FoxP1 ChIP-seq²⁰ peaks were called using the same pipeline, but with a q-value cutoff of 0.01, which were chosen to yield a peak count comparable to FoxP3. Finally, FoxP3 ChIP-exo peaks were defined as ±10 bp around reported ChIP-exo summits.⁴⁵ All peaks that overlap with blacklist regions⁴⁶ were removed. This led to final sets of FoxP3 ChIP peaks (n = 25,724), CNR peaks (n = 18,552), ChIP-exo peaks (n = 4,307), FoxP1 ChIP peaks (n = 35,767) and ATAC peaks (n = 502,437).

To evaluate the presence of motif pairs in ChIP peaks, we extracted the genomic sequence from the GRCm38.p6 mouse reference genome⁴⁷ for each peak and placed them in Pandas²⁶ (v 2.2.3) DataFrames. Each sequence was scanned for exact matches to the 400 motifs identified in Figure 2C (100 motif pairs in H-H, H-T, T-H and T-T configurations). Note that only the positive strand of the genomic sequence was used to identify sequence matches. Matches within a given motif pattern were non-overlapping, as the search advanced past each match before continuing, to prevent inflated counts from repetitive or sliding matches of the same motif. However, overlaps between different motif patterns were allowed to avoid search-order bias and ensure all valid matches were captured. To control for baseline motif frequency, we compared the frequency of each motif pair in FoxP3 ChIP peaks (or CNR or ChIP-exo peaks) to that in OCRs. We first normalized motif counts in ChIP peaks and OCRs with per-million base pair in their respective total peak lengths, and used exact binomial test⁴⁸ to identify motifs that have significantly higher normalized counts in ChIP peaks than in OCRs (p-value <0.05). For each of these motifs, the adjusted fold change (FC) of normalized counts in ChIP peaks over OCRs were visualized in the form of diameter scale in Figure 3A and heatmap in Figure 3B using Matplotlib²⁷ (v 3.9.4).

Motif gap size analysis was performed using similar approaches, where normalized counts (per million bp of their respective total peak lengths) of motif pairs with a gap size of 1–10 nt were measured in ChIP peaks and OCRs, and were plotted using Matplotlib²⁹ (v 3.9.4).

Analysis of relationships between TnG repeats and paired motifs

Defining genomic loci with TnG repeats

FIMO from MEME suite (v5.5.8)³⁵ was used to identify TnG-repeat-like elements. The TnG-repeat motif in our previous report¹⁷ was used as a query motif,⁴ and a search was performed against the mouse (GRCm38) genome, FoxP3 ChIP peaks or OCRs. A p-value cutoff of 8e-5 was used in FIMO results, and FIMO regions with repetitive 6TGs, 6ACs, 6TAs, 6TCs, 6AGs, 6CGs, 12Ts, 12Gs, 12Cs or 12As were removed to be consistent with FoxP3’s sequence specificity.¹⁷ The filtered TnG regions were exported using the python pyranges³⁶ (v 0.1.4) package.

Relationships between TnG repeats and H-T/T-H motifs

To examine the relationship between TnG repeats and H-T/T-H pairs of group 1-group 1 (G1-G1) motifs, we first compiled the genomic coordinates of H-T and T-H motif pairs within FoxP3 ChIP peaks using pyranges³⁶ (v 0.1.4) package. Only the G1-G1 motifs that are significantly enriched (p < 0.05) in FoxP3 ChIP peaks over OCRs were considered. We then calculated the percentage of H–T and T–H motifs that overlapped with TnG repeats defined above and plotted the values using Matplotlib²⁷ (v 3.9.4).

Relationships between TnG/H-H and differential gene expression

RNA-seq datasets^12,14,21 were acquired by mapping the reads to mm10 reference using splice-aware STAR mapper,⁴⁹ and exon reads counted using gencode mm10 annotation.⁵⁰ Deseq2 was used to calculate adjusted p-value and Log2FC values.⁵¹ Genes were stratified into significant (adjusted p-value <0.05) and non-significant groups. To quantify the enrichment of TnG or H–H features near these genes, we measured feature density (bp per Mb) within distance windows ranging from 1 kb to 1 Mb from each gene. Densities were computed using pyranges (v0.1.4). Separately, to evaluate whether proximity to TnG or H–H features explained differential expression, the nearest distance from each gene to TnG repeats or H–H motifs was measured. Genes were grouped by adjusted p-value bins and log2 fold change direction, and the proportion of genes in each bin was calculated.

Relationships between TnG repeats and H-H motifs

To assess how frequently H–H motifs and TnG repeats co-localize within the same FoxP3-bound regions or OCRs, each FoxP3 ChIP peaks or ATAC peaks examined for the presence of TnG repeats and H-H motifs using the pyranges³⁶ (v 0.1.4) package.

To investigate the spatial and orientational relationships between H-H motifs and nearest TnG repeats, the genomic center and borders of each H-H motif were calculated, and the nearest TnG repeat was identified based on border-to-border distance. H-H motifs that had overlapping coordinates due to similar motif structure were merged to avoid duplicates. If the H-H motif coordinates overlapped with the coordinates of a TnG repeat, the border-to-border distance was recorded as zero. Otherwise, we measured the border-to-border distance as positive when the TnG occurred downstream (from the end of the H-H to the start of the TnG), and negative when it occurred upstream (from the end of the TnG to the start of the H-H). To assess the orientation of the nearest TnG repeat, we used strand information associated with the overlapping TnG region identified via the FIMO³⁵ search. In cases where the TnG (positive strand) appeared upstream of the H-H, or CAn (negative strand) appeared downstream, the orientation was considered “aligned”. For comparison, the same analysis was performed in FoxP3 ChIP peaks and FoxP3-free OCRs (n = 328,369). The latters were defined as Treg ATAC-seq¹⁹ peaks located at least 10 kb from any FoxP3 ChIP-seq peak and with FoxP3 ChIP-seq signal below the minimum signal in FoxP3 ChIP peaks, as measured by area under curve (AUC) for a 200 bp window around the center of each TnG repeat.

To compare FoxP3 binding to TnG repeats with and without adjacent H-H motifs, we calcualted FoxP3 ChIP-seq AUC for a 200 bp window around the center of each TnG repeat. TnG repeats with 0 nt gap with nearest H-H motif were compared against those located at least 1 kb away from any H–H motif.

Genomic feature and histone ChIP-seq analysis

Genomic feature analysis was performed using ChIPseeker.⁵² To compare H3K4me3, H3K27ac and ATAC signal intensity, H3K4me3 and H3K27ac ChIP–seq and ATAC–seq data¹⁵ were trimmed by trimmomatic (v0.3.6)³⁴ and mapped to the mm10 genome using Bowtie2²⁹ (v2.3.4.3). The intensity was calculated within 5 kb upstream and downstream of centers of TnG repeats and centers of H–H motifs using deeptools³⁷ (v 3.5.6) bamCoverage and computeMatrix.

QUANTIFICATION AND STATISTICAL ANALYSIS

All quantification and statistical analyses were performed either with GraphPad Prism or R. The statistical tests used, the value of n, what n represents, the graph metrics, and the significance values are indicated in the figure legends and/or on the figures directly.

References

1. SA Lambert, A Jolma, LF Campitelli, PK Das, Y Yin, M Albu, X Chen, J Taipale, TR Hughes, MT Weirauch. The Human Transcription Factors.. Cell, 2018. [DOI | PubMed]
2. KR Nitta, A Jolma, Y Yin, E Morgunova, T Kivioja, J Akhtar, K Hens, J Toivonen, B Deplancke, EE Furlong, J Taipale. Conservation of transcription factor binding specificities across 600 million years of bilateria evolution.. eLife, 2015. [DOI | PubMed]
3. F Leng, W Zhang, RN Ramirez, J Leon, Y Zhong, L Hou, K Yuki, J van der Veeken, AY Rudensky, C Benoist, S Hur. The transcription factor FoxP3 can fold into two dimerization states with divergent implications for regulatory T cell function and immune homeostasis.. Immunity, 2022. [DOI | PubMed]
4. W Zhang, F Leng, X Wang, RN Ramirez, J Park, C Benoist, S Hur. FOXP3 recognizes microsatellites and bridges DNA through multimerization.. Nature, 2023. [DOI | PubMed]
5. F Leng, R Merino-Urteaga, X Wang, W Zhang, T Ha, S Hur. Ultrastable and versatile multimeric ensembles of FoxP3 on microsatellites.. Mol Cell, 2025. [DOI | PubMed]
6. S Dai, L Qu, J Li, Y Chen. Toward a mechanistic understanding of DNA binding by forkhead transcription factors and its perturbation by pathogenic mutations.. Nucleic Acids Res., 2021. [DOI | PubMed]
7. Y Wu, M Borde, V Heissmeyer, M Feuerer, AD Lapan, JC Stroud, DL Bates, L Guo, A Han, SF Ziegler. FOXP3 controls regulatory T cell function through cooperation with NFAT.. Cell, 2006. [DOI | PubMed]
8. HS Bandukwala, Y Wu, M Feuerer, Y Chen, B Barboza, S Ghosh, JC Stroud, C Benoist, D Mathis, A Rao, L Chen. Structure of a domain-swapped FOXP3 dimer on DNA and its function in regulatory T cells.. Immunity, 2011. [DOI | PubMed]
9. S Arora, J Yang, T Akiyama, DQ James, A Morrissey, TR Blanda, N Badjatia, WKM Lai, MSH Ko, BF Pugh, S Mahony. Joint sequence & chromatin neural networks characterize the differential abilities of Forkhead transcription factors to engage inaccessible chromatin.. Preprint at bioRxiv., 2023. [DOI]
10. M Ono, H Yaguchi, N Ohkura, I Kitabayashi, Y Nagamura, T Nomura, Y Miyachi, T Tsukada, S Sakaguchi. Foxp3 controls regulatory T-cell function by interacting with AML1/Runx1.. Nature, 2007. [DOI | PubMed]
11. RN Ramirez, K Chowdhary, J Leon, D Mathis, C Benoist. FoxP3 associates with enhancer-promoter loops to regulate T(reg)-specific gene expression.. Sci. Immunol., 2022. [DOI | PubMed]
12. Z Liu, DS Lee, Y Liang, Y Zheng, JR Dixon. Foxp3 orchestrates reorganization of chromatin architecture to establish regulatory T cell identity.. Nat. Commun., 2023. [DOI | PubMed]
13. X Song, B Li, Y Xiao, C Chen, Q Wang, Y Liu, A Berezov, C Xu, Y Gao, Z Li. Structural and biological features of FOXP3 dimerization relevant to regulatory T cell function.. Cell Rep., 2012. [DOI | PubMed]
14. J van der Veeken, A Glasner, Y Zhong, W Hu, ZM Wang, R Bou-Puerto, LM Charbonnier, TA Chatila, CS Leslie, AY Rudensky. The Transcription Factor Foxp3 Shapes Regulatory T Cell Identity by Tuning the Activity of trans-Acting Intermediaries.. Immunity, 2020. [DOI | PubMed]
15. Y Kitagawa, N Ohkura, Y Kidani, A Vandenbon, K Hirota, R Kawakami, K Yasuda, D Motooka, S Nakamura, M Kondo. Guidance of regulatory T cell development by Satb1-dependent superenhancer establishment.. Nat. Immunol., 2017. [DOI | PubMed]
16. CG de Boer, J Taipale. Hold out the genome: a roadmap to solving the cis-regulatory code.. Nature, 2024. [DOI | PubMed]
17. TL Bailey, J Johnson, CE Grant, WS Noble. The MEME Suite.. Nucleic Acids Res., 2015. [DOI | PubMed]
18. S Nakagawa, SS Gisselbrecht, JM Rogers, DL Hartl, ML Bulyk. DNA-binding specificity changes in the evolution of forkhead transcription factors.. Proc. Natl. Acad. Sci. USA, 2013. [DOI | PubMed]
19. TS Heng, MW Painter, C Immunological Genome Project. The Immunological Genome Project: networks of gene expression in immune cells.. Nat. Immunol., 2008. [DOI | PubMed]
20. Z Zhu, G Lou, XL Teng, H Wang, Y Luo, W Shi, K Yihunie, S Hao, K DeGolier, C Liao. FOXP1 and KLF2 reciprocally regulate checkpoints of stem-like to effector transition in CAR T cells.. Nat. Immunol., 2024. [DOI | PubMed]
21. C Jager, P Dimitrova, Q Sun, J Tennebroek, E Marchiori, M Jaritz, R Rauschmeier, G Estivill, A Obenauf, M Busslinger, J van der Veeken. Inducible protein degradation reveals inflammation-dependent function of the T(reg) cell lineage-defining transcription factor Foxp3.. Sci. Immunol., 2025. [DOI | PubMed]
22. D Zemmour, LM Charbonnier, J Leon, E Six, S Keles, M Delville, M Benamar, S Baris, J Zuber, K Chen. Single-cell analysis of FOXP3 deficiencies in humans and mice unmasks intrinsic and extrinsic CD4(+) T cell perturbations.. Nat. Immunol., 2021. [DOI | PubMed]
23. G Tang, L Peng, PR Baldwin, DS Mann, W Jiang, I Rees, SJ Ludtke. EMAN2: an extensible image processing suite for electron microscopy.. J. Struct. Biol., 2007. [DOI | PubMed]
24. B Bushnell, J Rood, E Singer. BBMerge – Accurate paired shotgun read merging via overlap.. PLoS One, 2017. [DOI | PubMed]
25. PJ Cock, T Antao, JT Chang, BA Chapman, CJ Cox, A Dalke, I Friedberg, T Hamelryck, F Kauff, B Wilczynski, MJ de Hoon. Biopython: freely available Python tools for computational molecular biology and bioinformatics.. Bioinformatics, 2009. [DOI | PubMed]
26. W McKinney, SM van der Walt. Proceedings of the 9th Python in Science Conference)., 2010
27. J Hunter. Matplotlib: A 2D graphics environment.. Comput. Sci. Eng., 2007. [DOI]
28. CR Harris, KJ Millman, SJ van der Walt, R Gommers, P Virtanen, D Cournapeau, E Wieser, J Taylor, S Berg, NJ Smith. Array programming with NumPy.. Nature, 2020. [DOI | PubMed]
29. B Langmead, SL Salzberg. Fast gapped-read alignment with Bowtie 2.. Nat. Methods, 2012. [DOI | PubMed]
30. P Danecek, JK Bonfield, J Liddle, J Marshall, V Ohan, MO Pollard, A Whitwham, T Keane, SA McCarthy, RM Davies, H Li. Twelve years of SAMtools and BCFtools.. GigaScience, 2021. [DOI | PubMed]
31. Y Zhang, T Liu, CA Meyer, J Eeckhoute, DS Johnson, BE Bernstein, C Nusbaum, RM Myers, M Brown, W Li, XS Liu. Model-based analysis of ChIP-Seq (MACS).. Genome Biol., 2008. [DOI | PubMed]
32. AR Quinlan, IM Hall. BEDTools: a flexible suite of utilities for comparing genomic features.. Bioinformatics, 2010. [DOI | PubMed]
33. R Leinonen, H Sugawara, M Shumway. The sequence read archive.. Nucleic Acids Res., 2011. [DOI | PubMed]
34. AM Bolger, M Lohse, B Usadel. Trimmomatic: a flexible trimmer for Illumina sequence data.. Bioinformatics, 2014. [DOI | PubMed]
35. CE Grant, TL Bailey, WS Noble. FIMO: scanning for occurrences of a given motif.. Bioinformatics, 2011. [DOI | PubMed]
36. EB Stovner, P Saetrom. PyRanges: efficient comparison of genomic intervals in Python.. Bioinformatics, 2020. [DOI | PubMed]
37. F Ramirez, DP Ryan, B Gruning, V Bhardwaj, F Kilpert, AS Richter, S Heyne, F Dundar, T Manke. deepTools2: a next generation web server for deep-sequencing data analysis.. Nucleic Acids Res., 2016. [DOI | PubMed]
38. M Ohi, Y Li, Y Cheng, T Walz. Negative Staining and Image Classification – Powerful Tools in Modern Electron Microscopy.. Biol. Proced. Online, 2004. [DOI | PubMed]
39. F Delaglio, S Grzesiek, GW Vuister, G Zhu, J Pfeifer, A Bax. NMRPipe: a multidimensional spectral processing system based on UNIX pipes.. J. Biomol. NMR, 1995. [DOI | PubMed]
40. WF Vranken, W Boucher, TJ Stevens, RH Fogh, A Pajon, M Llinas, EL Ulrich, JL Markley, J Ionides, ED Laue. The CCPN data model for NMR spectroscopy: development of a software pipeline.. Proteins, 2005. [DOI | PubMed]
41. SG Hyberts, K Takeuchi, G Wagner. Poisson-gap sampling and forward maximum entropy reconstruction for enhancing the resolution and sensitivity of protein NMR data.. J. Am. Chem. Soc., 2010. [DOI | PubMed]
42. SG Hyberts, AG Milbradt, AB Wagner, H Arthanari, G Wagner. Application of iterative soft thresholding for fast reconstruction of NMR data non-uniformly sampled with multidimensional Poisson Gap scheduling.. J. Biomol. NMR, 2012. [DOI | PubMed]
43. Y Shen, A Bax. Protein backbone and sidechain torsion angles predicted from NMR chemical shifts using artificial neural networks.. J. Biomol. NMR, 2013. [DOI | PubMed]
44. RM Samstein, A Arvey, SZ Josefowicz, X Peng, A Reynolds, R Sandstrom, S Neph, P Sabo, JM Kim, W Liao. Foxp3 exploits a pre-existent enhancer landscape for regulatory T cell lineage specification.. Cell, 2012. [DOI | PubMed]
45. N Yamada, MJ Rossi, N Farrell, BF Pugh, S Mahony. Alignment and quantification of ChIP-exo crosslinking patterns reveal the spatial organization of protein-DNA complexes.. Nucleic Acids Res., 2020. [DOI | PubMed]
46. HM Amemiya, A Kundaje, AP Boyle. The ENCODE Blacklist: Identification of Problematic Regions of the Genome.. Sci. Rep., 2019. [DOI | PubMed]
47. AD Yates, P Achuthan, W Akanni, J Allen, J Allen, J Alvarez-Jarreta, MR Amode, IM Armean, AG Azov, R Bennett. Ensembl 2020.. Nucleic Acids Res., 2020. [DOI | PubMed]
48. P Virtanen, R Gommers, TE Oliphant, M Haberland, T Reddy, D Cournapeau, E Burovski, P Peterson, W Weckesser, J Bright. SciPy 1.0: fundamental algorithms for scientific computing in Python.. Nat. Methods, 2020. [DOI | PubMed]
49. A Dobin, CA Davis, F Schlesinger, J Drenkow, C Zaleski, S Jha, P Batut, M Chaisson, TR Gingeras. STAR: ultrafast universal RNA-seq aligner.. Bioinformatics, 2013. [DOI | PubMed]
50. A Frankish, M Diekhans, AM Ferreira, R Johnson, I Jungreis, J Loveland, JM Mudge, C Sisu, J Wright, J Armstrong. GENCODE reference annotation for the human and mouse genomes.. Nucleic Acids Res., 2019. [DOI | PubMed]
51. MI Love, W Huber, S Anders. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.. Genome Biol., 2014. [DOI | PubMed]
52. G Yu, LG Wang, QY He. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization.. Bioinformatics, 2015. [DOI | PubMed]