A multiple redundant genetic switch locks in the transcriptional signature of T regulatory cells

Bioinformatic prediction of Treg transcriptional control

Transcriptional regulation is governed by extensive and interconnected networks of genes. This complexity can be tackled by computational methods that start from a large number of gene expression datasets and reconstruct plausible regulatory models, and infer and rank potential connections between target genes and a set of putative regulators ^{33, 34}. These algorithms, typically based on multiple regression or related approaches, analyze the pairwise variation between regulator(s) and target across a large number of related datasets, in response to a range of perturbations (differentiation, genetic or chemical perturbations). Here, we used 129 gene-expression profiles previously generated on the M430.v2 microarray platform from various CD4⁺ T cells: primary Treg and Tconv cells from various anatomical locations and of different surface phenotypes, TGF-β-induced FoxP3⁺ cells, cells from mutant mice (Rara, Foxp3 deficiencies), and Akt- or various TF-transfectants (Supplementary Table 1). We selected as potential regulators 2021 transcription-control factors from GeneOntology annotation (conventional TFs as well as chromatin modifiers), and 603 target genes that compose the canonical Treg signature (407 and 196 over- or under-expressed in Tregs, respectively; Fig. 1a) ¹¹. The Context Likelihood of Relatedness (CLR) algorithm ³⁵ was used, a relevance network reconstruction method that operates by combining the relative strength of coexpression between a regulator and potential targets. The results are listed in Supplementary Table 2, the top regulators shown in Table 1, and Fig. 1b. Reassuringly, top predicted regulators included FoxP3 and other factors previously associated with Treg function such as Eos (Ikzf4) and Helios (Ikzf2) ^{23, 36}, but also some novel TFs not previously associated with Tregs such as Lef1 or Gata1. Some of the predicted regulators were themselves differentially expressed in Treg versus Tconv cells (Ikzf2, Ikzf4, Lef1), while others were only modestly so (Gata1).

Many target genes were predicted to be influenced by several TFs (Fig. 1b), making it difficult to infer the regulators of the Treg signature. In keeping with our hypothesis of additive transcriptional control by a panel of TFs, we started from these predicted regulators and formulated an optimization process on the ILOG Cplex package (IBM) to determine a set of TFs that would, in combination, account for the greatest fraction of the Treg signature. Under this model, 10 TFs could explain over half of the Treg signature (Fig. 1c). FoxP3 led the list, covering 10.8 % of the signature, lower but in the same range as estimates from transduction experiments ¹¹. Most of the TFs were predicted to be both stimulatory or repressive, depending on the target, although some seemed to be only activating (Gata1, Hdac9).

Loss of function validation of computational predictions

A set of complementary gain- and loss- of function experiments was undertaken to validate whether the computational predictions had actual biological relevance. We chose a subset of TFs based on availability of knockout mice and/or enforced-expression vectors. First, FoxP3’s predicted targets were analyzed in a comparison of TGF-β-induced cultures of CD4⁺ cells from wildtype or FoxP3-deficient scurfy mice, which brings forth those transcripts strictly dependent on Foxp3 ¹¹. Predicted FoxP3 targets were skewed to the extremes of the distribution, more so than the Treg signature as a whole (Fig. 2a; p=6.9×10^-6) suggesting the validity of computational prediction of FoxP3 targets.

We then analyzed the transcriptomes of Treg cells in a set of knockout mice available for several of these predicted cofactors: i) a complete knockout of Eos; these mice are viable and fertile, with no noted autoimmune phenotype, and have normal Treg numbers and phenotypes (Supplementary Fig.1, and R.G. and D.R., unpublished), perhaps contrary to expectations ²³; ii) a promoter deletion of Gata1 ³⁷, in which Treg and other T cells appear normal (Supplementary Fig. 1, and J.H., unpublished); other known Gata1 target genes were affected in those mice. iii) a conditional knockout of the Xbp1 gene (Xbp1^fl/fl × Mx1-cre) ³⁸; Treg populations in lymphoid organs are again normal in these mice (Supplementary Fig. 1, and S.A. and L.G., unpublished); iv) a knockout of Helios in which Tregs appear normal ³⁹. Gene-expression profiles were generated from purified splenic CD4⁺CD25^hi Tregs of these mice and their wildtype littermates (Fig. 2b). No bias was detected in any of the mutant Tregs, whether of the Treg signature as a whole (Fig. 2b) or of the computationally predicted targets of each TFs (Supplementary Fig.2) Thus, while each of these TFs may impact on the Treg signature when varying naturally within the diverse cell types used for the computational analysis, the Treg signature was robust to the complete elimination of any one of them.

Gain-of-function validation of computational predictions

We then performed gain-of-function experiments by retrovirally transducing cDNAs encoding a number of candidate TFs, alone or together with FOXP3, into CD4⁺ Tconv cells activated with anti-CD3+CD28 beads (human FOXP3, which has a very comparable transcriptional signature relative to that of mouse FoxP3 ¹¹ was used to allow distinction from the endogenous transcript). These manipulations resulted in expression levels of FoxP3 and other TFs in the same general range as found in normal ex vivo T cells (Supplementary Fig. 3). Cells expressing each TF, alone or together with FOXP3, were sorted 3 days after transduction for gene expression profiling (Fig. 3a). Consistent with previous findings ^{11, 20}, FOXP3 alone could influence only a limited number of Treg signature genes (FOXP3, Fig. 3b). Enforced expression of cofactors alone had even less effect (EOS, Fig. 3b), but robust induction of Treg Up signature genes and repression of Treg Down signature genes were observed when FOXP3 and cofactors were both present (FOXP3+EOS, Fig. 3b). This was validated by RT-PCR for representative genes in independent samples (Fig. 3c). Such a synergistic outcome was seen with each of the 7 candidate TFs tested (Supplementary Fig. 4), but five of them (EOS, IRF4, GATA1, LEF1 and SATB1, hereafter “quintet”) had a striking effect, cooperating with FOXP3 to similarly shift most of the Treg signature (Fig. 3d and Supplementary Table 3). Indeed, each of this quintet of TFs, together with FOXP3, regulated the same set of genes, as shown by the direct comparison of Fig. 3e. This synergy was not an artifact of the dual transductions, as cells transduced with FOXP3 and a control TF (Pbx1) were similar to those expressing FOXP3 alone (Fig. 3b, 3f). Rates of cell division were also identical in singly- and doubly-transduced cells, as measured by CFSE dilution (Supplementary Fig. 5). This response was different from that prompted by Helios or XBP1, although the latter also synergized with FOXP3, as shown by the integrated Treg signature index (Fig. 3f).

We then asked whether the combination of two quintet TFs could induce the Treg signature, without FoxP3. Indeed, the combination of EOS+LEF1, or of GATA1+SATB1) had a partial effect, including a modest induction of Foxp3 (Supplementary Fig. 6).

Affymetrix 1.0 ST microarrays contain features that span the length of the transcripts, allowing us to parse signals from the transfected genes vs. their endogenous homolog (Fig. 4a and Supplementary Table 4). FOXP3 plus any of the quintet factors modified the expression of endogenous TF transcripts: induction of Foxp3, Eos, Irf4 and Xbp1 and repression of Lef1 and Satb1. Thus, the introduction of any of the quintet factors synergized with FOXP3 to induce a widespread reassortment of the cell’s regulatory TF balance, in an autoassembly of the Treg profile. This involved the induction of TFs that are normally over-expressed (Eos, Irf4) and the repression of those that are under-expressed (Lef1, Satb1) in Treg cells (Supplementary Fig. 7). Thus, these results indicate a synergistic effect between FoxP3 and cofactors that propagates to other TFs and locks-in the Treg signature. Accordingly, the genes affected here include the Treg signature genes found to be FoxP3-independent in previous studies ¹¹.

How the cofactors operate

We then asked how the quintet cofactors might elicit this transition. It was not through stabilization of the FoxP3 protein, whose abundance, measured by intracellular staining, was identical whether or not the cells were co-transduced with a quintet factor, across a range of expression values of the co-transcribed IRES-Thy1.1 reporter (Fig. 4b). It was possible that quintet factors have a quantitative influence on FoxP3’s activity, simply displacing a threshold of activity below which FoxP3 would be ineffective. To test this hypothesis, we sorted and profiled matched bins of FOXP3-transduced cells bearing various levels of FOXP3, alone or co-transduced with GATA1, chosen as a representative of the quintet factors (Fig. 4c). As might be expected, increasing FOXP3 did have a more substantial transcriptional impact. But the highest levels of FOXP3, when alone, were unable to match the induction of Treg signature genes together with GATA1. The cooperating effect of GATA1 was apparent at all levels of FOXP3. Thus, the quintet factors were not merely providing a quantitative boost to FoxP3, but instead enhanced its transcriptional activity. The nuclear/cytoplasmic distribution of FOXP3 was unchanged, being almost exclusively nuclear whether transduced alone or together with a quintet factor (Fig. 4d).

These effects also suggested that FoxP3 interacts molecularly with the quintet factors within nuclear complexes. Such interactions have already been demonstrated for Irf4 and Eos ^{23, 24}, so we tested the other three. Indeed, reciprocal co-immunoprecipitation in transduced cells showed an interaction between FoxP3 and GATA1, SATB1 and LEF1, but not with the control TF Pbx1 (Supplementary Fig.8).

The synergizing activity of the cofactors, most dramatically observed with quintet factors, and not explained by quantitative effects on FoxP3 or its global cellular localization, could have two interpretations: first, through cooperative binding, that the cofactors recruit FoxP3 to genomic locations close to Treg signature genes; second, that the cofactors enhance the activity of FoxP3 molecules already bound to DNA. To distinguish between these scenarios, we used chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) to assess how quintet factors affect the genome-wide localization of FoxP3. Chromatin was prepared from primary CD4⁺ Tconv cells transfected with FLAG-FoxP3 alone or together with GATA1 (the preparation of double-transduced cells in the numbers needed for ChIP-seq was technically very demanding, so we chose GATA1 as a representative), immunoprecipitated with anti-FLAG, and the bound DNA determined by Illumina deep sequencing (see Supplementary Table 5 for ChIP-seq statistics). Immunoprecipitation with anti-PolII, or whole cell extract, provided genome-wide controls for transcriptional start sites (TSS) or for sequencing non-homogeneity, respectively. Overall, summing the genome-wide signals relative to TSS locations showed that FoxP3 predominantly localized in the vicinity of known TSS, as expected (Fig. 5a; in other experiments, ChIP-seq with irrelevant control antibody showed a paucity of signal around the TSS ⁴⁰, substantiating the significance of signals observed here). The data allowed statistically robust detection of more than 5,000 FoxP3-binding sites (MACS p<10^-7; Supplementary Table 6). Many of these sites were corroborated by comparison with similar data from ex vivo Tregs, kindly provided by R. Samstein and A. Rudensky). To further validate these data we computed, across the range of genes showing significant peaks of FoxP3 binding, the distribution of genes whose expression affected by transduction of FoxP3 and cofactor. As might be expected, the genes with the highest FoxP3 binding were enriched for those activated by FoxP3 in the microarray data (Fig. 5b; of the 57 genes with FoxP3 binding peak height >75, 12.2% had a FoldChange after transduction >1.6, vs 4.7% in the whole dataset, p=0.008). This was not an absolute, however, and many sites of strong FoxP3 binding did not correspond to significant transcriptional change in expression, as often observed in ChIP-seq data. In addition, transcripts repressed by FoxP3 were not over-represented among those with highest FoxP3 binding.

In cells doubly-transduced with FoxP3+GATA1, we did not observe novel sites of significant FoxP3 binding. Rather, there was a quantitatively enhanced occupancy by FoxP3 at the same locations, as illustrated for the FoxP3-binding site in the first intron of Icos (Fig. 5c), or when genome-wide binding was quantitated in parallel (Fig. 5d). Thus, the quintet factors do not spread the recruitment of FoxP3 to different genomic locations, but seem to functionally stabilize, and enhance the activity of, FoxP3 independently bound to its target sites.

Signature lock by feedback loops: computational modeling

The results depicted in Fig. 3d were puzzling: how could five distinct TFs, of widely different structure, DNA sequence specificity, and functional activity, elicit the same transcriptional outcome? This was particularly paradoxical for Lef1 and Satb1, which are repressed in Treg cells ²⁶ (Supplementary Fig. 7); Of note, however, is that the retroviral vectors used only contain the SATB1 or LEF1 coding regions, and lack the 3’- and 5’-UTR regions, which have been shown to be involved in the regulation of endogenous Satb1 expression ²⁶; the transduced version lacks these controls, and likely leads to ‘constitutive’ expression of LEF1 and SATB1 during the culture period). A plausible interpretation was suggested by the effect on endogenous TF expression (Fig. 4a): the Treg signature, with the regulatory factors it includes, is organized with regulatory feedback loops, both positive and negative, such that it has the capacity to ‘auto-assemble’ and lock-in once the expression of FoxP3 and some cofactors is pushed beyond that of Tconv cells. Intuitively, such locking-in could be achieved by positive feedback, but also by double-negative inhibition of repressive factors.

To assess the plausibility of this intuition, we used computational simulation to ask whether such a self-reinforcing system that incorporated repressed regulators could actually function. A mathematical model was developed to simulate the dynamics of such a system (Fig. 6). The model consists of the main regulator F, with its active conformation F* (F to F* transition can mean either quantitative induction, post-translational modifications (e.g. acetylation) or complex formation, as suggested by co-IP assays, that potentiate or stabilize F), and a set of downstream regulatory factors of the Treg signature, either up- or down-regulated by F (U_i, D_i; Fig. 6a) (see details in Supplementary Information). These signature genes themselves control smaller sub-networks, some of which are pure effectors (U₄, D₃), while others can regulate F* (e.g. U_1-3, D_1-2). Output of these subnetworks, which themselves can operate with AND or OR logic, then control a larger set of signature genes, but most importantly influence the F to F* conversion.

To leave the model computationally manageable, subnetwork calculations were bypassed, and cross-regulatory influences between cofactors, which are likely to occur, were omitted. Differential equations paired up with Hill functions described the biochemical kinetics engaged in the model. After fixing a reasonable parameter set, this minimal model successfully reproduced the bistable features of the Treg program, and the outcome of the experimental perturbations (Figs. 2, 3). After in silico transduction (Fig. 6b), Treg signature genes remained off when FoxP3 or any of the cofactors were expressed singly, but all signature genes transitioned to the Treg state, and were locked in, when FoxP3 was overexpressed together with either of the cofactors, including those repressed in Tregs (e.g. D₁). The model showed no effect of the single in silico knockout of any of the cofactors (Fig. 6c), consistent with experimental results, but the Treg signature shut off with extinction of FoxP3. The latter deviated somewhat from experimental results, as thymi of mice with an inactivated FoxP3 protein do contain cells with some Treg features, including partial activation of the Treg signature and a transcriptionally active Foxp3 locus ^{17, 18}. The discrepancy could be resolved in the model by postulating that the differentiation of Treg cells triggers, directly or indirectly, the transient expression of FoxP3 and one of its cofactors (for instance, U₃ in Fig. 6d, top). In the simulation, transient activation of F and U₃ resulted in the activation of the whole network (Fig. 6d, bottom left), but only if the external inducing conditions were modeled to activate both F and U3. This activation was then stable upon removal of the inductive stimulus. With inactive FoxP3, however, cells showed only partial Treg features, which reverted to the Tconv state some time after withdrawal of differentiating stimulus (Fig. 6d, bottom right), a scenario consistent with the unstable “Treg wannabes” mentioned above.

Thus, the simulation arrived at a model of Treg cell differentiation compatible with most experimental outcomes and with several aspects of Treg physiology.

Mice

C57BL/6 (B6) mice were purchased from the Jackson Laboratory. GATA1 promoter mutant mice on the Balb/c background ³⁷ carrying a deletion within the double GATA-site 21 bp upstream of the first hematopoietic exon, were obtained from The Jackson Laboratory. Xbp1^f/f × Mx1-Cre conditional knockout mice on the B6 background have been described ³⁸; 5-6 weeks old mice were intraperitoneally injected 3 times with 250 μg of poly(I:C) (invivoGen) with 2 days intervals to induce Cre expression. Mice were used for experiments 6 weeks after the final poly(I:C) injection; Helios (Ikzf2) KO mice have been described ³⁹; Eos (Ikzf4) knockout mice were generated by first inserting LoxP sites flanking exons 1 through 4, and then crossing with a germ-line Cre to generate constitutively deficient mice. Targeting of the genomic locus was validated by Southern blot analysis using 5’ and 3’ external probes, deletion of exon 1-4 was confirmed using quantitative RT-PCR. Homozygous Eos-deficient mice are viable and fertile with no apparent abnormality. All mice cared for in accordance with the ethical guidelines of the Center for Animal Resources and Comparative Medicine at Harvard Medical School under Institutional Animal Care and Use Committee (IACUC) approved procedures (protocol 02954).

CLR Network Regulatory Prediction

For this analysis we compiled 129 previously published gene expression datasets obtained from purified CD4⁺ T cell populations in several experimental contexts: ex vivo conventional T (Tconv) or regulatory T (Treg) cells from anatomical locations, cultured Treg cells, TGFβ-induced FoxP3⁺ cells, retinoic acid treated cultures ^{11, 12, 48}. The Affymetrix M430v2 microarray raw data were preprocessed with the RMA algorithm in GenePattern ⁴⁹, and averaged expression values were used for analysis.

For a robust definition of the transcriptional signature of mature Treg cells, results from several independent experiments had been brought together. The consensus peripheral Treg cell signature had been defined by calculating the Treg to Tconv fold change (FC) ratios and retaining only those genes that showed a consistent 1.5-fold overexpression or underexpression in Treg cells in all four datasets, using Affymetrix M430v2 arrays. This resulted in a total of 603 genes, 407 overexpressed and 196 underexpressed in Treg cells respectively ¹¹.

For prediction of regulatory connections, we used the CLR (Context Likelihood of Relatedness) algorithm ³⁵ that operates by combining the relative strength of coexpression between a regulator and potential targets. The CLR algorithm builds upon the relevance network strategies by applying a background correction step. First we reconstructed a relevance network between 2021 transcriptional regulators and 603 Treg signature genes ¹¹. After computing mutual information (MI) values between all pairwise TF-Treg signature gene pairs, the algorithm compares the MI between a TF/gene pair to the background distribution of MI scores of all genes associated with the TF or all TFs associated with the gene of interest. After this background correction, the most probable interactions are those whose combined scores stand significantly above the background distribution of MI scores. The background corrected CLR scores were filtered at a false discovery rate (FDR) =0.005, which was computed with Bonferroni correction, to generate the Treg CLR network. CLR computations were performed in Matlab (MathWorks).

In the second phase, starting from the CLR scores we formulated an optimization problem whose objective was to identify the TFs that together influence the Treg signature the most and account for the most number of genes. This is a mixed integer optimization problem which we solved using the CPLEX9.0 optimization package (ILOG) for AMPL.

Retroviral transduction

For enforced expression of FOXP3 and other TFs, the retroviral expression vector MSCV-IRES-Thy1.1/GFP ¹¹ was used throughout. The human FOXP3, GATA1, EOS, IRF4, LEF1, SATB1, XBP1 cDNAs were obtained from human ORFeome. Mouse Helios cDNA was kindly provided by Stephen Smale. Plat-E cells were plated one day before and transfected with these plasmids, together with a packaging construct pCL-Eco using Lipofectamine 2000 (Invitrogen), according to the manufacturer’s instructions.

For primary CD4⁺ T cell transduction, cell suspensions were prepared from spleens and lymph nodes of 6-8 week-old B6 mice by physical dissociation and red blood cells were lysed in 0.8% Ammonium Chloride lysis buffer. CD4⁺ T cells were negatively purified by magnetic selection (labeling with phycoerythrin-conjugated anti-CD11b (M1/70), -CD11c (N418), -CD19 (6D5), -CD8α (53-6.7), -CD25 (PC61), and -NK1.1 (PK136). After washing, anti-PE beads (Miltenyi Biotec, #130-048-801) were added to the cell suspension and subsequently, and CD4⁺ Tconv cells were purified using MACS LD columns (Miltenyi) to purity > 95%. Cells were then activated with anti-CD3/CD28 beads (Invitrogen) at a ratio of one bead per cell, with addition of 20 U/ml of recombinant human IL-2 (Proleukin; Chiron) in complete culture medium (RPMI 1640 supplemented with 10% fetal calf serum, 2 mM L-Glutamine, 100 U/ml penicillin/streptomycin and 50 μM 2-Mercaptoethanol). T cells were cultured for 24 hours and then were spin-infected (2000 rpm, 32°C, 2 hours) with retrovirus supernatants. Cells were then cultured for an additional 72 hours. Infected cells were sorted by flow cytometry, as CD4⁺ cells further gated on Thy1.1 and GFP that denote expression of FOXP3 and of the other cofactors, respectively. For the experiments shown in Fig. 4b,c, infected T cells were sorted into different fractions based on the intensities of Thy1.1 expression (high, intermediate and low, respectively) for microarray profiling and FOXP3 protein analysis.

To verify and quantitate FOXP3 expression in these transfectants, sorted cells were fixed and permeabilized for intracellular staining with anti-FoxP3 (eBioscience) according to manufacturer’s instructions, and were either analyzed by flow cytometry for quantitative FOXP3 protein expression, or by confocal microscopy for FOXP3 protein localization.

Gene expression profiling

For analysis of gene expression in knockout mice or after retroviral transduction, sorted cell populations were lysed in TRIzol reagent, and RNA was prepared according to the manufacturer‘s instructions (Invitrogen). RNA amplification was conducted for two rounds using the MessageAmp aRNA kit (Ambion), followed by biotin labeling using the BioArray high yield RNA transcription labeling kit (Enzo Life Sciences, Inc.), and purified using the RNeasy mini kit (QIAGEN). The resulting cRNAs were hybridized to Mouse Gene 1.0 ST arrays (Affymetrix). These steps followed the ImmGen pipeline and were performed at Expression Analysis, Inc (Durham, NC). Data were normalized with the RMA algorithm implemented in Affymetrix Power Tools after first pre-filtering to remove unannotated probes and visualized on GenePattern Multiplot module.

We developed a Treg signature index to estimate the global expression of Treg signature genes in the transfectants (Fig. 3f). First, we calculated the fraction of signature genes up-regulated under various conditions (F); then we calculated the median value of the fold-change relative to control transfectant for all Up signature genes (M), and the Treg Up signature index was established as be I_Up=F*M*2. As expected, the I_Up equals 1 in controls. Similar calculation was done for Treg Down signature genes, and a composite Treg signature index was calculated as I=[I_Up+I_Down]/2.

To distinguish the expression of transduced TFs from that of their endogenous counterparts, we used feature-level analysis of the 1.0 ST microarray data. The Affymetrix Mouse Gene 1.0 ST Array offers whole-transcript coverage, as each of the 28,853 genes is represented on the array by approximately 27 oligonucleotides (“features”) spread across the full length of the gene. This characteristic allowed us to distinguish the expression of mouse endogenous TFs from the transgene TFs, which mostly were human origin. First, nucleotide sequences of all the features (25-mer oligonucleotides) for one particular gene (e.g. Foxp3) were retrieved and the sequence similarities between mouse and human were analyzed using NCBI Blast tools. Features with strong dissimilarity (more than 5 mismatches among 25 nucleotides) between mouse and human were considered as mouse-specific probes, and their expression values averaged and normalized to arrive at the values shown in Fig. 4a and Supplementary Table 4.

Immunoprecipitation (IP) and Immunoblotting (IB)

293 cells double-transfected with vectors for Flag-FoxP3 plus any of the following TFs: SATB1, LEF1, GATA1 or Pbx1, lysed on ice with hypotonic solution (10 mm HEPES, 1.5 mM MgCl2, 10 mM KCl, and 0.05% NP-40 like/IgePal Ca-630) supplemented with EDTA-free complete protease inhibitors (Roche). Nuclear pellets were subsequently treated with nuclear lysis buffer (20 mM HEPES, 300 mM NaCl, 20 mM KCl, EDTA-free complete protease inhibitor cocktail) and MNase (Nuclease S7; Roche). Chromatin digestion was stopped by adding EDTA to 5 mM, and post-nuclear supernatants were incubated with Protein-G Sepharose beads coupled to antibodies for IP (Flag, M2, Sigma; FoxP3, FJK-16s, eBioscience; GATA1, Ab28839, Abcam; LEF1, Ab124271, Abcam; SATB1, 611182, BD; Control IgG) overnight at 4°C with constant rotation. Bound proteins were eluted by boiling, separated by SDS-PAGE, and electro-transferred to PVDF. After blocking (2 hr in 5% milk/1x PBS 0.02% Tween20), blots were probed 1 hour at room temperature with antibodies for IB.

Chromatin immunoprecipitation and sequencing

Mouse primary CD4+ T cells, transduced and sorted as above, were used in this assay. Chromatin immunoprecipitation (ChIP) was done as described ⁵⁰. Briefly, ~107 cells were cross-linked with formaldehyde (11%). Cell lysates were sonicated (8 cycles of 30” at 60” intervals, on ice; Misonix), incubated with 10 μg antibodies (anti-Pol-II (total) (sc-899, Santa Cruz); anti-FLAG-FoxP3 (M2, Sigma); anti-GATA1 (ab28839, Abcam), which were pre-bound with protein G-conjugated Dynal beads (Invitrogen). Immunoprecipitated DNA was purified and used for library construction using ChIP-Seq DNA Sample Prep Kit for Illumina sequencing ⁵⁰. Sequences were aligned to the genome using Bowtie software (ver. 0.12.7) to NCBI Build 36 (UCSC mm9) of the mouse genome. Peaks of binding were called with MACS software (1.4.0rc2). The number of reads in each tag pileup were first normalized relative to the total number of reads in the sample. To accurately compare the local tag densities in peak regions of the different samples (particularly for FoxP3 binding in the samples transduced with either FoxP3, or FoxP3+GATA1), values were rescaled by a constant, which was calculated from the integrated values of the noise in regions devoid of any FoxP3-binding peaks (7 regions ranging from 60 to 650 Kb). This correction stemmed from the assumption that the experimental noise should be constant even when true signal (and hence the total number of reads) might be expected to vary between parallel samples, and that a normalization factor calculated from the genome background level allows appropriate compensation for variability in amplification during the construction of sequencing libraries.

Mathematical Modeling

Based on the experimental findings, we developed a mathematical model to describe the kinetic and dynamic behaviors of the regulation networks governing Treg cells, and to test which modes of regulatory connectivity might account for the peculiar behavior of this network.

The Treg signatures defined in this manuscript consist of the expression profiles of 603 genes that are differentially expressed: Treg relative to Tconv. Among them, 407 genes are up-regulated and the rest 196 genes are down-regulated. Regardless of the diversity of underlying regulatory architectures, the essential features of the FoxP3-regulated T cell network can be mimicked by a simple model as illustrated in Figure 6a in the main text. This model consists of two master species F and F*, representing the native FoxP3 and its modified functional form^#1 (See the text box for definition), and a set of downstream genes U_i’s and D_i’s. These genes mimic correspondingly the up- and down-regulated signature genes^#2 in Treg cell’s regulation program. They are further arranged in parallel based on our experimental evidence that Treg signature is robust to the full elimination of any one of signature factors. In this model, FoxP3 can transit from its native (F) to functional (F*) state with the mediation from its cofactors^#3. Conversely, the functional form promotes the production of its native. Although both forms of FoxP3 regulate signature genes, the modified one has a much stronger regulation effect than the native form. In addition, among all of the signature genes, some subsets are straight downstream genes with no feedbacks (pure effectors^#4, such as U4 and D3), while others are capable of regulating the transition of FoxP3 from its native form to modified form (cofactors^#3, such as U1 and D1).

To simplify our model while still capturing the major picture, we omitted many intermediate details of the Treg program. In a real Treg cell, the regulation network is likely much more complicated than the diagram we drew here. For instance, the regulation of FoxP3 onto the signature genes can be indirect and implemented through intermediate subnetworks and other signature genes. Similarly, the cofactors likely regulate the transition of FoxP3 via intermediate species. Nevertheless, this simplified model remains the key features of the regulation program, enabling us to investigating its core function and kinetic and dynamic behaviors for a mechanistic understanding of the system.

1. Biochemical reactions in Treg’s regulatory program

Biomolecular events engaged in the regulation of Treg consist mainly of two classes of reactions: (1) Bindings and dissociations of Promoter-FoxP3-cofactor complex; (2) Productions and degradations of reactant molecules, including FoxP3 and signature genes. Molecular species participating in the program include FoxP3, signature genes (proteins), promoter, FoxP3-protein complex, promoter-FoxP3 complex, and promoter-FoxP3-protein complex, as listed in Supporting Table M1. Detailed biochemical reactions are listed in Supporting Table M2.

It is important to point that only cofactors, not pure effectors, are able to interact with FoxP3 to form functional forms (F*) and consequently promote the expression of downstream genes. Although we generalized the reactions for all signature genes in the Supporting Table M2 without distinguishing differences between cofactors and effectors, it can be simply addressed by choosing appropriate rate constants. For instance, if a signature gene (X_i) is a pure effector, we can assign an infinity (practically, a very large number) to the corresponding dissociation constant W_i to implement this situation where the formation of FoxP3-protein complex is prohibited.

2. Mathematical model

With the network diagram shown in Fig. 6a and the corresponding reactions listed above, we can derive our mathematical model, using the fast reaction arguments for the binding and dissociation of operators with corresponding transcription factors, as follows

where U_j and D_j are the concentrations of representative up- and down-regulated signature genes, and F_t is total concentration of the native FoxP3 (F) and its modified complex (F*). X is a vector referring to the total signature genes, i.e., X = (U₁, U₂, …, U_n1, D₁, D₂, …, D_n2) where n₁ and n₂ are the numbers of up- and down-regulated genes respectively. α_s0, α_s and γ_s (s = u_j, d_j, f) are the rate constants for the basal production, regulated production, and degradation of the gene s. The function H_s (F_t, X) (s = u_j, d_j, f) is the hybrid production rate of the gene s that is co-regulated by F_t and X as

where n is the total number of signature genes (up- and down-regulated genes), i.e., n = n₁ + n₂.

For the simplest case where there is only one up-regulated and one down-regulated genes, the above model can be simplified as:

with the corresponding hybrid protein production function as

where C_0,s, C_1,s, C_2,s, C_3,s are the folds of change for the expression of species s (s=u, d, f) upon activation when promoters are bound by nothing, original FoxP3, FoxP3-U complex, and FoxP3-D complex respectively.

3. In silico overexpression and knockout experiments

Although the above mathematical model is capable of describing the complete signatures (603 genes) of the Treg program, we here, for simplicity, use a seven-species network as a representative of the whole network to explore the program’s network feature through in silico experiments. In this simplified circuit, four of the signature genes are up-regulated (U₁, U₂, U₃, U₄) and the other three are down-regulated (D₁, D₂, D₃), as indicated in Fig. 6a. Among the seven species, five of them (U₁, U₂, U₃, D₁ and D₂) are cofactors that facilitate the transition of FoxP3 from its native to functional forms while the other two (U₄, D₃) are pure effectors with no regulatory input.

To begin with, we chose a set of parameters from commonly used and biologically reasonable parameter space for this model (See Supporting Table M3 for details). The parameters include the rate constants of the basal production, activated production, and degradation for the four up-regulated genes (α_u0,i, α_u,i, γ_u,i, where i=1-4 corresponds to the gene U₁, U₂, U₃, and U₄), the three down-regulated genes (α_d0,i, α_d,i, γ_d,i, where i=1-3 corresponds to the gene D₁, D₂, D₃), and total FoxP3 (α_f0, α_f, γ_f). The parameter set also contains C_0-8,i, the folds of change for the expression of species i (i=1-8, corresponding to U₁, U₂, U₃, U₄, D₁, D₂, D₃, F_t), when promoters are bound respectively by nothing, original FoxP3, FoxP3-U₁, FoxP3-U₂, FoxP3-U₃, FoxP3-U₄, FoxP3-D₁, FoxP3-D₂, and FoxP3-D₃. Here C_5,1-8 and C_8,1-8 are set all zeros because the 5th and 8th genes (U₄ and D₃) are pure effectors and do not regulate the expression of any gene. Additionally, the model contains M₁₋₈, the dissociation constants of FoxP3-protein complexes ([O-F, O-FU₁, O-FU₂, O-FU₃, O-FU₄, O-FD₁, O-FD₂, O-FD₃]), and W₁₋₇, the dissociation constants of FoxP3- signature proteins complex (F-U₁, F-U₂, F-U₃, F-U₄, F-D₁, F-D₂, F-D₃). Here W₄ and W₇ are set as 10²⁰ to match the facts that both U4 and D3 are pure effectors and do not regulate the transition of FoxP3.

With the above parameter set, our model successfully mimics the bistability feature of T cell’s regulatory program: All up-regulated genes remain in their low expressions while those down-regulated remains high when the program starts with the conventional state (Tconv) (first panel, top row of Fig. 6b); However, when the program is in the regulatory state (Treg) (first panel, Fig. 6c), the whole expression profiles are opposite and remain opposite. Both two states are stable with the same parameter set and their final state is primarily determined by their initial conditions.

4. Activation of the Treg regulatory program

The above overexpression and knockout experiments successfully illustrated the robustness of Treg’s regulation network, and all of our computational results are consistent with the wet lab data from both this and previous studies except the last panel in Fig. 6c: In our numerical experiment, knocking out FoxP3 fully shut off the Treg program while previous studies showed that a small population of T cells still developed partial signatures of Treg cells in transgenic FoxP3 mutant mice that expressed an inactive FoxP3 protein ^17,18. This intriguing discrepancy could be resolved by postulating an upstream “triggering factor” that promotes the expression of FoxP3 as well as other intermediate sub-networks of the Treg program (note that “triggering factor” may be a single factor or a combination of factors activated in concert).

We here numerically investigated the Treg differentiation process (left panel of Fig. 6d) by simulating the induction of FoxP3 together with a cofactor (here U₃) by the external “triggering factor” (how the external “triggering factor” is actually induced, or whether it corresponds to a single factor or a coordinated combination thereof are not relevant here). We revised our mathematical model (Eq. 1) to incorporate the reactions as follows

where I₁ and I₂ are two intermediate species and X is the external upstream factor (X is different from the signature gene vector X). The last two equations in our original 7-species model are now revised to include the production induction from the intermediate species. To simulate the process, we chose the following common parameters for the Hill function: α_I1 0 = 0.0, α_I1 = 0.3, γ_I1 = 0.3, θ_x1 = 2, W_x1 = 10.0, α_I20 = 0.0, α_I2 = 0.25, γ_I2 = 0.25, θ_x2 = 2, W_x2 = 10.0, α_xu3 = 60.0, θ_I1 = 2, W_I1 = 10.0, α_xf = 50.0, θ_I2 = 2, and W_I2 = 100.0.

To mimic the induction process, we started by simulating FoxP3-negative wild type cells for a sufficient time so that it achieved a steady state. This was followed by a step-like increase of the external triggering factor(s), which lasts for 15 time units before removal. The result (the middle panel of Fig. 6d) shows that the transient induction from the upstream factor turned on the whole program (FoxP3 and seven signature genes) and the whole expression profiles remained stable even after the removal of external triggering factor(s).

We also conducted a numerical assay for the case where FoxP3 is expressed but is functionally inactive (as in the inactivating knock-in insertions from the Rudensky and Chatila laboratories). Depicted in the right panel of Fig. 6d, a subset of the signature genes (U₃ and FoxP3) showed expression profiles similar to those of Treg cells in a transient time window after the induction. Although these genes eventually returned to the conventional state, the transient period showed partial features of a Treg cell.

5. Parameter exploration

In addition to reproducing the experimental findings, we can further use our model to uncover the system’s behaviors that are not explored in our wet-lab experiments. One particularly interesting feature of this system is the bifurcation of the dynamics. A similar question is when the bistability of the program will disappear. To address that, we evaluated the system’s behavior by proportionally scaling all of the dissociation constants (M₁, M₂, …, M₈) and simulating the system for each scaled constant set using two different initial conditions [102, 86, 91, 80, 5.8, 3.1, 8.3, 96] and [6, 4, 2, 7.5, 90, 85, 90, 5] (for the variables [U₁, U₂, U₃, U₄, D₁, D₂, D₃, F_t]). These two initial conditions are chosen from the two distinct final steady states obtained in Fig. 6c and Fig. 6b. The figure below shows the final steady state of the system respecting to the scaling factor: The left panel corresponds to the first initial condition above, the middle panel is from the second initial condition, and the right panel is the overlap of FoxP3 in the left and middle panels. The results show that the system has two stable states, corresponding to the conventional and regulatory states, when varying the scaling factor from 10⁻² to 10^0.65. In other words, the systems settle in different final states when starting with different initial conditions. However, the system always arrives at the same set of final solutions regardless of its initial condition when the scaling factor is greater than 10^0.65, indicating that the system becomes monostable (conventional state). It seems intriguing but makes biological sense: The dissociation constants indicate how easily FoxP3 and its modified forms escape from the promoters of signature genes, which anti-correlates with the production rates of downstream genes. Thereby, higher dissociation constants result in lower expression levels of signature genes including those serving as cofactors, which subsequently lowers the transition rate for FoxP3 from its native to functional forms and hence decreases the productions of all signature genes. Once the constants go below a threshold, the Treg state becomes unstable and the system is incapable of sustaining in it in regardless of its initial state and, as a consequence, the system always settles in the monostable conventional state.

Statistical analysis

Statistical significance was determined with Kolmogorov-Smirnov test.