The pioneer transcription factors Foxa1 and Foxa2 regulate alternative RNA splicing during thymocyte positive selection

During positive selection at the transition from CD4CD8 double-positive (DP) to single positive (SP) thymocyte, TCR signalling results in appropriate MHC-restriction and signals for survival and progression. We show that the pioneer transcription factors Foxa1 and Foxa2 are required to regulate RNA splicing during positive selection. Foxa1 and Foxa2 had overlapping/compensatory roles. Conditional deletion of both Foxa1 and Foxa2 from DP thymocytes reduced positive selection and development of CD4SP, CD8SP and peripheral naïve CD4T-cells. Foxa1 and Foxa2 regulated expression of many genes encoding splicingfactors and regulators, including Mbnl1, H1f0, Sf3b1, Hnrnpa1, Rnpc3, Prpf4b, Prf40b and Snrpd3. Within the positively selecting CD69DP cells, alternative RNA splicing was dysregulated in the double Foxa1/Foxa2 conditional knockout, leading to >850 differentially used exons (DEU). Many genes important for this stage of T-cell development (Ikzf1-3, Ptprc, Stat5a, Stat5b, Cd28, Tcf7) and splicing-factors (Hnrnpab, Hnrnpa2b1, Hnrnpu, Hnrnpul1, Prpf8) showed multiple DEU. Thus Foxa1 and Foxa2 are required during positive selection D ev el o pm en t • A cc ep te d m an us cr ip t to regulate alternative splicing of genes essential for T-cell development, and that by also regulating splicing of splicing-factors, they exert widespread control of alternative splicing.


Introduction
The production of T-cells in the thymus involves multiple stages of development during which haematopoietic precursors give rise to mature T-cells that can differentiate into functional effector T-cells. During this process, progenitors cells which do not express the coreceptor molecules CD4 and CD8 (CD4 − CD8 − double negative (DN) cells) differentiate to become CD4 + CD8 + double positive (DP) cells, which give rise to both CD4 single positive (SP) and CD8SP populations. Maturation from DP to SP follows successful rearrangement of the Tcr locus, and requires TCR signalling: positive selection results in appropriate MHC restriction of SP cells, and is followed by negative selection of potentially self-reactive clones and selection of regulatory T-cells (Treg) (1)(2)(3). The strength and duration of the TCR signal that a developing cell receives broadly determine its fate, with the strongest signals leading to negative selection or CD4 Treg differentiation, usually at the SP stage in the medulla, intermediate signals leading to positive selection usually in the cortex, and weaker signals or lack of TCR signalling leading to death by neglect (4). For DP thymocytes undergoing positive selection, TCR signal strength and duration also influence CD4 and CD8 lineage choice.
Those cells receiving stronger and longer TCR signals tend towards the CD4SP fate, while weaker/more transient signals favour CD8SP fate, and fate decisions are also influenced by Little is known about the function of Foxa1 and Foxa2 in the immune system. In mouse models of autoimmune inflammation, ectopic Foxa1 expression has been shown to drive differentiation and suppressive function of a novel subset of Treg (37). We recently found that expression of Foxa1 and Foxa2 in thymic epithelial cells (TEC) is required to maintain normal T-cell development and homeostasis of thymic and spleen Treg populations.
Conditional deletion of Foxa1 and Foxa2 from TEC increased the proportion of medullary TEC, but reduced cell-surface MHCII expression on TEC, leading to a smaller thymus with a reduction in conventional CD4 T-cell differentiation but an increase in the CD4 Treg population (25). In the thymus, Foxa1 and Foxa2 are also expressed in developing T-cells, and Foxa2 is a transcriptional target of Shh signalling after pre-TCR signal transduction (12,26).
Here we investigate the function of Foxa1 and Foxa2 in developing T-cells. We show that Foxa1 and Foxa2 are required at the transition from DP to SP T-cell, where they regulate RNA splicing. Conditional deletion of both Foxa1 and Foxa2 led to a reduction in positive selection, differentiation and maturation of the SP populations, and a reduction in the peripheral CD4 T-cell pool.

CD4SP and CD8SP development are impaired in the Foxa1/2cKO thymus
We examined Foxa1 and Foxa2 expression by QRT-PCR in FACS-sorted developing thymocytes from the DN3 stage onwards (Fig 1A). Foxa1 and Foxa2 were detected in all stages, and showed reciprocal patterns of expression, with Foxa1 expression upregulated in DN4, following by decline in DP cells. In contrast, Foxa2 was highly expressed in DN3 and DP stages, and expressed at lower levels in DN4 cells. Both Foxa1 and Foxa2 showed higher expression in CD4SP than in CD8SP cells.
To establish if Foxa1 and Foxa2 play a T-cell intrinsic role in thymocyte development, we conditionally deleted Foxa1 and/or Foxa2 from T-cells from the DP stage of development onwards, by crossing mice carrying a single or double loxP-flanked Foxa1 and/or Foxa2 allele to mice in which Cre is driven by the CD4-promotor. Foxa1 fl/fl CD4cre + and Foxa2 fl/fl CD4cre + Development • Accepted manuscript mutant mice are referred to as Foxa1cKO and Foxa2cKO, respectively. Foxa1 fl/fl Foxa2 fl/fl CD4cre + double mutant mice are referred to as Foxa1/2cKO. Foxa1 and Foxa2 were effectively deleted in Foxa1/2cKO CD4SP and CD8SP thymocytes, as expression of Foxa1 and Foxa2 were below detection by QRT-PCR, but detected in the control (Foxa1 fl/fl Foxa2 fl/fl CD4cre -) CD4SP and CD8SP populations (Fig 1B).
We compared thymocyte populations in Foxa1cKO, Foxa2cKO and Foxa1/2cKO thymus with their control littermates. The number of thymocytes were not significantly different between the conditional Foxa1 and/or Foxa2 mutants and their controls (Fig 1C). The Foxa1cKO thymus contained normal proportions of DP and SP populations (Fig 1D), but there was an increase in the proportion of DP cells in the Foxa2cKO compared to control ( Fig   1E,G). In the Foxa1/2cKO thymus, the proportion and number of the DP population were increased and the proportion and number of CD4SP and CD8SP populations were reduced ( Fig 1F-I). The fact that while the phenotype of the double Foxa1/2cKO thymus was more pronounced than that of the Foxa2cKO, the Foxa1cKO thymus appeared grossly normal, suggests that Foxa1 can compensate for Foxa2 at this developmental transition, and we therefore decided to use the double Foxa1/2cKO mice to investigate the impairment of differentiation from DP to SP cell.
The ratios of DP:SP, DP:CD4SP and DP:CD8SP were all increased in Foxa1/2cKO compared to control, while the ratio of CD4SP:CD8SP was decreased ( Fig 1J). Thus, deficiency of both Foxa1 and Foxa2 in thymocytes led to less efficient development of SP populations and bias towards the generation of CD8 lineage cells over CD4SP. During differentiation from DP to SP stage, cell surface TCR expression is upregulated. The proportion of TCR hi thymocytes was decreased in the absence of Foxa1 and Foxa2 (Fig 1K). When we gated on TCR hi cells, and compared subset distribution, we found an increase in the proportion of DP and CD4 + CD8 lo intermediate cells, and a decrease in the proportion of CD4SP population, indicating that fewer cells were completing positive selection and a partial arrest at the transition from DP and CD4 + CD8 lo to CD4SP (Fig 1L).

Foxa1/2 deficiency influences maturation of SP cells and peripheral T-cell populations
We then investigated the maturation status of SP cells in the Foxa1/2cKO thymus by expression of HSA, CD69, Qa2 and CD62L. After positive selection, SP thymocytes retain high expression of heat-stable antigen (HSA (CD24)) and CD69, and as they mature HSA and CD69 are downregulated, Qa2 is upregulated, and expression of CD62L indicates that thymocytes are mature and ready to egress from the thymus. The proportion of HSA + CD69 + cells in both CD4SP and CD8SP populations was significantly decreased in Foxa1/2cKO compared to control thymus (Fig 2A). Furthermore, there was a significant decrease in the percentage of Qa2 + and CD69 -CD62L + cells in the CD4SP compartment in Foxa1/2cKO compared to control (Fig 2B-C), whereas we did not detect significant differences in Qa2 and CD62L expression in the CD8SP population ( Supplementary Fig 2). Thus, in the absence of Foxa1 and Foxa2, not only were there fewer CD4SP and CD8SP cells, but their maturation was also affected, such that fewer mature CD62L + CD4SP were produced to egress from the thymus. In contrast, we found no significant difference in the proportion of CD4 + CD25 + Foxp3 + Tregs, NKT cells or NK cells in the Foxa1/2cKO thymus compared to control (Fig 2D,E).
In the spleen and lymph nodes changes in the CD4+ T-cell populations mirrored the thymus in the Foxa1/2cKO compared to control, with reductions in the proportion and number of conventional CD4 + T-cells overall and of naive CD44 -CD62L + CD4 + T-cells in spleen (Fig 2F-H).
In contrast, we did not detect significant differences in the number of peripheral CD8 + Tcells in the Foxa1/2cKO compared to control, suggesting that despite the reduction in CD8SP cells in the thymus, the peripheral CD8+ T-cell compartment is subject to homeostatic control and can expand to reach its normal size ( Fig. 2F-G).

Foxa1/2 deficiency reduces TCR signalling
TCR signal strength is one factor which determines positive selection in the thymus and tonic TCR signalling is again required for maintenance and homeostasis of peripheral T-cell populations after egress from the thymus. Therefore, as a proxy to measure TCR signal strength in thymocyte subsets, we compared expression of proteins whose levels of expression are directy determined by TCR signalling between Foxa1/2cKO and control.

Development • Accepted manuscript
Expression of intracellular Nr4a1 is induced as an early consequence of TCR signal transduction and high expression requires relatively strong TCR signalling (38,39). The proportion of intracellular Nr4a1 hi cells was significantly reduced in Foxa1/2cKO DP, CD4SP and CD8SP populations compared to control, suggesting that in the absence of Foxa1 and Foxa2 fewer cells had reached the threshold of TCR signal strength required to induce high levels of Nr4a1 (Fig 3A). Levels of cell-surface CD5 expression correlate with the strength of TCR signal transduction that a developing T-cell has received (40). As expected, mean fluorescence intensity (MFI) of anti-CD5 staining was lower in CD8SP and DP cells than in CD4SP cells (Fig 3B). MFI of CD5 was lower in Foxa1/2cKO DP, CD4SP and CD8SP populations than in their control counterparts, consistent with reduced TCR signal strength in the absence of Foxa1 and Foxa2.

Foxa1 and Foxa2 promote positive selection in DP thymocytes
Positive selection signals for survival and further maturation and requires interactions between self-peptide:MHC complexes in the thymic cortex, which may take place over a number of days, involving multiple or prolonged TCR-MHC interactions (41)(42)(43). TCR signalling for positive selection leads to cell-surface CD69 expression followed by upregulation of the cell-surface TCR complex. The proportions of DP cells that expressed cell surface TCR hi and CD69 were decreased in the Foxa1/2cKO thymus compared to control ( Fig 3C-D), consistent with a reduction in positive selection. To test this and to investigate mechanisms that might account for the reduced transition from DP to SP, we carried out transcriptome analysis on cells undergoing positive selection by RNA-sequencing FACSsorted CD69 + DP cells. RNA-sequencing identified only 176 differentially expressed genes (DEG) between Foxa1/2cKO and control datasets (adjusted (FDR) p<0.05), of which 109 (62%) were more highly expressed in Foxa1/2cKO than control ( Fig 3E and Table S1). Of the 176 DEG, 85 (~48%) had previously been identified as Foxa1/2 targets in genome-wide ChipSeq analysis of Foxa1/2 binding sites in dopaminergic neuronal progenitors and 56 of those genes with verified Foxa1/2 binding sites (approximately two thirds) showed higher expression in the Foxa1/2cKO than control (Fig 3F and Table S1) (34,44). DEG included genes involved in T-cell development and function, and approximately one quarter of DEG were genes whose transcription has been shown to be regulated during positive selection (45) (Fig   3G and Table S1). Among these, DEG that were more highly expressed in Foxa1/2cKO than control datasets included genes known to reduce TCR signal strength: Cbl, a ubiquitin ligase, which negatively regulates TCR signalling (46) and Themis, which can attenuate TCR signalling during repertoire selection (47); and also genes associated with the CD8 lineage: Cd8b1; and Lyst, a lysosomal trafficking regulator required for CTL lytic granules (48), consistent with the CD4 lineage being more severely affected in the conditional knockouts.
To test in an unbiased way, if Foxa1 and Foxa2 are required for the transcriptional response to TCR signalling for positive selection, we then used canonical correspondence analysis (CCA) to compare the overall pattern of gene expression in our datasets with transcriptome data from publically available datasets prepared from DP thymocytes that were receiving different strengths of TCR signals during selection (GSE38909) (49). The GSE38909 dataset contains TCR-transgenic AND DP thymocytes stimulated with a positively selecting peptide (gp250) or a non-selecting control peptide (49). We selected the 2000 most significantly differentially expressed genes from the GSE38909 dataset between DP thymocytes stimulated with the non-selecting control peptide and DP thymocytes stimulated with the positively selecting peptide and used these to generate a scale of unstimulated to TCRsignalling-for-positive-selection. We plotted our datasets against this scale. The CCA segregated the datasets by genotype: both control datasets were on the positive side of the axis, which corresponded to the transcriptional pattern induced by the positively selecting peptide, consistent with the fact that CD69 + DP thymocytes have initiated positive selection by TCR signal transduction ( Fig 3H). In contrast, both Foxa1/2cKO datasets fell on the negative side of the axis, thus showing an overall pattern of transcription that is closer to that of unstimulated DP cells than that of their control counterparts.
To confirm this impact on positive selection, we subdivided thymocytes by cell surface TCR and CD69 expression into four different stages: TCR lo/neg CD69 -(pre-selection thymocytes), TCR int CD69 + (intermediate transition, undergoing positive selection), TCR hi CD69 + (thymocytes after TCR signalling for initiation of positive selection) and TCR hi CD69 -(more mature population). The proportion of pre-selection thymocytes was higher in the Foxa1/2cKO than control thymus, whereas there was no significant difference in the proportion of the TCR int CD69 + population, and both TCR hi CD69 + and TCR hi CD69populations were reduced, indicating that Foxa1 and Foxa2 expression in thymocytes promotes initation of the process of positive selection, and also progression of cells during the differentiation process ( Fig 3I). The loss of cells from the TCR hi CD69 + population, but not its precursor TCR int CD69 + population, suggests impairment at a late stage of positive selection and a failure to progress, and consistent with this, the proportion of cells undergoing cell death (Annexin V + ) was increased in the CD69 + DP population in the Foxa1/2cKO thymus compared to control (Fig 3J). When we gated on the most mature TCR hi CD69thymocyte population, and compared CD4/CD8 subset distribution, we found that the proportion of DP cells increased by more than two-fold in the Foxa1/2cKO thymus compared to control, whereas the proportion of CD4SP cells was significantly reduced ( Fig   3K), demonstrating a clear requirement for Foxa1 and Foxa2 in normal differentiation to CD4SP and positive selection, and suggesting that differentiation is dysregulated in the Foxa1/2cKO, so that more cells that had upregulated cell surface TCR and downregulated CD69 were unable to progress beyond the DP stage.

Foxa1 and Foxa2 regulate exon usage
To investigate how Foxa1/2 might regulate positive selection we identified genes that encode verified transcription factors amongst the DEG between Foxa1/2cKO and control CD69 + DP datasets to look for known regulators of positive selection or differentiation at this developmental stage which might function downstream of Foxa1/2. Twenty-four transcription factors were found, of which 13 had previously been verified as Foxa1/2 targets in genome-wide ChipSeq analysis in neuronal progenitors ( Fig 4A and Table S1).
Amongst the transcription factors whose expression were downregulated by Foxa1/2 deletion there were no obvious down-stream candidates that might function to promote positive selection, although several showed changes in expression consistent with the developmental phenotype: for example Ptma, an anti-apoptotic gene (50); and Notch3 which is transcriptionally regulated during positive selection but not required for positive selection (45,51). Likewise, we found no obvious candidate genes which might function to inhibit positive selction or TCR signal strength amongst those genes which encode transcription factors and whose expression were upregulated in the absence of Foxa1/2.
Ikzf2 was more highly expressed in Foxa1/2cKO than control, but is not required for differentiation at this stage of development; Elk4 was also more highly expressed, but is required for positive selection and upregulated by TCR signal transduction, rather than functioning as a negative regulator (39,(52)(53)(54).
We then carried out Gene ontology (GO) term enrichment analysis of all DEG between Foxa1/2cKO and control CD69 + DP datasets. This revealed over-representation of genes associated with terms connected with immunity and T-cell mediated immunity, including lymphocyte migration, cytokine signalling, cell killing and apoptotic process (Fig 4B and   Table S2). Thus, as well as an overall enrichment in immunity related terms, as expected from the cell type, there was enrichment for terms related to processes required for differentiation to SP, such as lymphocyte migration, as differentiating cells must migrate to the medulla to complete their maturation and repertoire selection. Several terms related to cytokine signalling are also pertinent to the partial arrest at the DP and CD4 + CD8 lo stage in the Foxa1/2 cKO, as cytokine signalling rescues cells whose MHCI-restricted TCR signal has been interrupted by down-regulation of CD8, allowing differentiation to CD8SP, whereas cytokine signalling is inhibited for differentiation to the CD4SP compartment.
Additionally, there was enrichment for the term 'chromatin silencing', related to the known function of Foxa transcription factors as pioneer factors. The GO analysis also showed enrichment for several terms related to RNA splicing, and for terms connected to cell matrix adhesion (Fig 4B). Amongst the top 20 most statistically significant DEG were genes associated with these processes (Fig 4C): for example, Mmp14, a membrane bound matrix metalloproteinase involved in breakdown of the extracellular matrix(55); H1f0, an H1 linker histone necessary for the condensation of nucleosome chains into higher-order chromatin structures, which is involved in the regulation of mRNA splice site recognition (56,57); and Mbnl1, a regulator of alternative splicing, which functions in the control of T-cell development, and is over-expressed in mixed lineage leukemia (58)(59)(60)(61).

Development • Accepted manuscript
Several other RNA splicing genes were differentially expressed between Foxa1/2cKO and control datasets (Fig 4D). Given this, and the fact that the GO enrichment analysis identified several terms associated with RNA splicing, we hypothesized that Foxa1/2 might regulate the transition from DP to SP cell and positive selection by influencing RNA splicing in developing T-cells. To test this, we compared differential exon usage (DEU) between Foxa1/2cKO and control CD69 + DP RNA-sequencing datasets. This analysis identified 852 events of differential exon usage (FDR p<0.05) between the conditional knockout and control, which involved 628 different genes, 222 (~35%) of which had previously been shown to bind to Foxa1/2 in genome-wide ChipSeq analysis of Foxa1/2 binding sites in dopaminergic neuronal progenitors (Fig 4E and Table S3). Conditional deletion of Foxa1/2 therefore led to a greater number of changes in exon usage affecting more genes than the number of individual genes that were differentially expressed. Intersection between DEG and genes which contained differentially used exons revealed just 22 DEG (whose overall expression was differentially regulated and which also showed differential usage of individual exons) (Fig 4F), 18 of which (>80%) had previousy verified Foxa1/2 binding sites (Table S1). This intersection included several genes involved in the regulation of RNA splicing, with known splicing variants ( Fig 4G): Ranbp2, a nucleoporin protein, controls alternative splicing patterns during nuclear speckle formation(62); Sf3b1, a well known splicing factor, mutations of which lead to myelodysplasia and anemia by globally disrupted splicing (63,64); and Smg1 which is mutated in Acute Myeloid Leukemia, its depletion resulting in disruption of alternative splicing (65,66).
GO term enrichment analysis of genes which showed DEU revealed over-representation of genes associated with terms connected with the known functions of Foxa1/2 in metabolic processes (metabolic process, cellular glucose homeostasis, response to insulin) and with their known functions in other tissues as epigenetic regulators and pioneer transcription factors (e.g. chromosome organization, chromatin remodeling, regulation of histone methylation, DNA methylation, DNA conformational change) (Fig 5A and Table S4). Many terms associated with mRNA splicing were also over-represented (Fig 5A). Genes involved in mRNA splicing showed multiple changes in exon usage, indicating that mRNA splicing factors are themselves subject to alternative splicing in developing T-cells, and that the regulatory Development • Accepted manuscript effects on splicing of the Foxa transcription factors may be amplified by regulation of splicing of components of the splicing machinery (Fig 5B-D). Additionally, the enrichment analysis highlighted positive regulation of NFK signalling, a pathway which has been shown to regulate alternative splicing in T-cells (67).

Foxa1 and Foxa2 regulate alternative splicing of essential genes for T-cell development
GO terms associated with T-cell development and function showed over-representation of DEU genes (Fig 6A): for example, thymic T-cell selection; T-cell receptor signalling pathway; regulation of CD8-positive  differentiation; lymphocyte differentiation; CD4-positive alpha-beta T-cell differentiation. Thus, Foxa1/2 control exon usage for genes involved in processes essential for the transition from DP to SP cell. Multiple differentially used exons were found in key genes for this transition, such as Ikzf1, Ptprc, Stat5a, Stat5b, Cd28 and Klf13 ( Fig 6B) (9,(68)(69)(70)(71), and many other essential genes for T-cell development showed differential usage of a single exon (for example: Notch1, Rorc, Socs1, Pten, Orai1, Cbfb and Ets2) (Fig 6B and Table S3). Many of the T-lineage genes which contained significantly differential exon usage have been previously described to be alternatively spliced during lymphocyte development or activation, such as the Ikaros family members, Fyn, Ptprc (CD45), Stat5a, Stat5b and Cd28 (72,73) (Fig 6B-E). However, the differential exon usage caused by absence of Foxa1 and Foxa2 did not always correspond to well-described alternative spliced variants of these genes. For example, alternative splicing of Ptprc exons 4, 5 and 6 is functionally important in T-cell differentiation (73), but in our datasets despite variation in expression of exon 5, only exons 2, 9, 13 were significantly different between control and conditional knockout after adjustment for false discovery (p<0.05) (Fig 6B and   E). Likewise, for Ikzf1, Ikzf2 and Ikzf3 we identified differentially expressed exons that were distinct from the alternative splice variants described in mouse thymocytes (74), indicating that absence of Foxa1 and Foxa2 led to wide dysregulation in the RNA splicing of genes that would normally display splice variants (Fig 6B-C).

Development • Accepted manuscript
Given the many changes in exon usage observed in genes required for T-cell development and positive selection, it seemed likely that Foxa1/2's function during positive selection is to regulate RNA splicing of essential genes for developmental progression. Therefore, to investigate the impact of Foxa1/2 on the normal splicing of regulators of development, we identified transcription factors within the genes that showed DEU between Foxa1/2cKO and control. Of the 628 genes that showed DEU, 97 (15.5%), encoded transcription factors ( Fig   7A and Table S2), many of which are important during T-cell development, and are specifically required for the transition from DP to SP cell (for example, Yy1, Trim28, Tcf7, Tcf3, Tcf12, Stat5a, Stat5b, Sp3, Rorc, Klf13, Cbfb) (Fig 7B).

Discussion
Here we show that the pioneer transcription factors Foxa1 and Foxa2 regulate alternative RNA splicing during T-cell development at the transition from DP to SP thymocyte.
Conditional deletion of Foxa1 and Foxa2 from developing DP thymocytes led to reduced positive selection and a partial arrest at the transition from DP to SP thymocyte, with a reduction in SP cell maturation, and a reduced peripheral naïve CD4 T-cell population.
Conditional deletion of Foxa1 and Foxa2 also led to significant changes in the expression of genes that regulate RNA splicing in cells undergoing positive selection, and concomitantly to >850 significantly differentially used exons.
Alternative RNA splicing is a mechanism that enables cells to generate many different proteins from a limited number of genes and is important in the regulation of development processes, including T-cell development (60,67,72,75,76). Our study indicates that Foxa1 and Foxa2 regulate the mRNA splicing of many genes that are important for progression at this developmental transition, and we therefore propose that aberrant RNA splicing of multiple genes may account at least in part for the reduction in positive selection and differentiation to SP cell observed when Foxa1 and Foxa2 were conditionally deleted from DP cells. In support of this, many of the genes that showed DEU, but were not DEG, between Foxa1/2cKO and control datasets are well-known to be required at this transition (eg . Socs1, Ptprc, Rasgrp1,); and many encode transcriptional regulators of this stage of development (eg . Ikzf1, Stat5a, Stat5b, Cbfb, Tcf7, Tcf3, Tcf12, Klf13, Sp3, Rorc). GO term Development • Accepted manuscript enrichment analysis also highlighted terms associated with positive selection and differentiation to SP. In contrast, few DEG encoded known relevant transcription factors or regulators of TCR repertoire selection or differentiation from DP to SP, and Themis, a DEG with verified Foxa1/2 binding site, which was more highly expressed in the conditional knockout, and which can modulate TCR signalling in thymocytes, also showed differential exon usage. The way in which the transcriptional activity of Foxa1/2 relate to their regulation of splicing will require further investigation, as there is increasing evidence that splicing can occur co-transcriptionally as well as post-transcriptionally (77). Interestingly, approximately one third of genes which showed DEU but were not DEG had previously been shown to bind directly to Foxa1/2 in neuronal progenitors, suggesting that Foxa1/2 may act directly to regulate their splicing.
Overall, we identified only 176 DEG between Foxa1/2cKO and control CD69 + DP cells, and of these more than 60% were more highly expressed in the conditional knockout than control, suggesting either that they were directly repressed by Foxa1/2, most likely by association of Foxa1/2 with a co-repressor, or that Foxa1/2 activate the transcription of an intermediate transcriptional repressor. Approximately 48% of DEG have been previously been shown to bind Foxa1/2 in whole genome ChipSeq screen of neuronal progenitors, suggesting that they were also likely to be direct targets (directly bound by Foxa1/2) in developing T-cells.
These were distributed between DEG that were up-or down-regulated in the Foxa1/2cKO compared to control, and ~51% of up-regulated genes have previously been shown to bind Foxa1/2 by ChipSeq screening in other cell types (34,44). We therefore think it is most likely that absence of Foxa1/2 in CD69 + DP cells led to increased expression of these DEG because of a direct repressive impact of Foxa1/2 binding to these sites and recruitment of a corepressor in control cells, rather than indirectly by transcriptional activation of an unknown downstream repressor of transcription. Although in cell transfection assays Foxa factors behave as transcriptional activators, Foxa1/2 have previously been shown to have repressor activity in other developing tissues, and Foxa2 has been shown to interact with the transcriptional co-repressor Tle family of proteins (34,78,79). Further investigation of Foxa2/Tle interactions in the regulation of differentiation from DP to SP cell will therefore Development • Accepted manuscript be important, given that Tle1, Tle3 and Tle4 are together required for commitment to the CD8 lineage at the transition from DP to SP (80). Foxa1/2 conditional deletion led to increased transcription of several RNA splicing regulators, including Sf3b1, Smg1, and Mbnl1, mutation of which are associated with haematological malignancies and aberrant RNA splicing, but also significantly decreased transcription of splicing factors (Prpf40b, Hnrnpa1, Snrpd3) and the histone linker H1f0, which may be required for mRNA splice site recognition (56,(58)(59)(60)(61)(63)(64)(65)(66). Overall loss of Foxa1 and Foxa2 led to broad changes in splicing events, indicating that Foxa1 and Foxa2 are important regulators of mRNA processing during T-cell development. Foxa1/2 also regulated the mRNA slicing of many genes that encode splicing factors, suggesting that their influence on the splicing machinery is further amplified by regulation of alternative splicing of splicing components, as well as by a more direct effect on the transcription of splicing regulators, RNA binding and histone linkers proteins.
Mbnl1 is highly expressed in the thymus and a recent study demonstrated that its constitutive knock out led to a hyperplastic thymus with retention of thymocytes and many mis-splicing events (60). Mbnl1 was upregulated approximately two-fold in CD69 + DP cells in the absence of Foxa1 and Foxa2 and was one of the most significantly differentially expressed genes in our datasets. In Foxa1/2cKO DP thymocytes we observed differential exon usage in some of the same genes that were affected when Mbnl1 was consititutively deleted: for example, Tcf7 exon(E)10 was affected by absence of Mbnl1, whereas our datasets showed changes in Tcf7 E1 and E6; Map4k4 E20 was affected in Mbnl1-/-thymus, but Foxa1/2cKO led to changes in Map4k4 E15 and E33; Sptan1 E23 was affected in Mbnl1-/-thymus, but Foxa1/2cKO changed usage of Sptan1 E2, E11, E12 and E49(60).
Conditional deletion of Foxa1 alone from DP thymocytes did not grossly affect T-cell development, and the double Foxa1/2cKO had a stronger impact on differentiation from DP to SP than the single Foxa2cKO, suggesting overlapping or redundant functions, so that while Foxa1 may not be required at this developmental transition it can partially replace the requirement for Foxa2. Foxa1 and Foxa2 have compensatory and partially redundant roles in other tissues, including liver (22) and a recent study showed that conditional deletion of the third Foxa protein (Foxa3) in addition to Foxa1 and Foxa2 in adult liver abrogated liver gene regulatory networks, destroying liver tissue homeostasis and function in the adult (30).
It will therefore be important in the future to investigate the potential compensatory role of Foxa3 in T-cell development.
All mice were bred and maintained at University College London (UK) under UK Home Office regulations.
In this paper, we will refer these mice as Foxa1/2cKO. Similarly, the Foxa1 flox/flox CD4-cre + is termed Foxa1cKO and the Foxa2 flox/flox CD4-cre + is termed Foxa2cKO. For all experiments their control genotype is the CD4-crelittermate with the same floxed alleles as the experimental animals.
The PCR conditions for CD4-cre transgene were 1 min at 94°C, 1 min at 61°C, and 1 min at 72°C for 35 cycles. Foxa1 and Foxa2 WT and floxed gene primers and PCR conditions were as described (25).

RNA-sequencing and data analysis
CD69 + CD4 + CD8 + thymocytes from Foxa1/2cKO and control thymus were FACS sorted. Each sample was sorted independently from different mice, with 2 mice of each genotype per sort. RNA was prepared as described (84) and sequenced by UCL Genomics on the Illumina Next Seq 500. The sequenced data are publically available in the GEO database under accession number GSE169602. The Basespace Sequence Hub was used for both FASTQ generation and RNA-Seq Alignment to the mouse reference genome UCSC mm10 (RefSeq gene annotation). Aligned reads were counted by the HTSeq python package with "union" overlap resolution mode (85). The Bioconductor package DESeq2 (1.30.0) was used to test for differential expression (86). P-values were plotted in a histogram resulting in a hill-shape, indicating an overestimation of the variance in the null distribution ( Supplementary Fig 1).
Therefore, the z-scores returned by DESeq2 were used as input to the CRAN package fdrtool (1.2.16) to re-estimate the null variance and subsequently the p-values (87). Adjusted pvalues were then calculated by Benjamini-Hochberg false discovery correction (5%). Genes with adjusted p-values less than 0.05 were considered as DEGs.
The Bioconductor package DEXSeq (1.36.0) was used to test for differential exon usage (88,89). The Python script dexseq_prepare_annotation.py was used to prepare the genome annotation with the following parameter: "-r no". Subsequently, dexseq_count.py was used to generate counts of exons using the BAM alignment files as input with the following parameters "-p yes -s no -a 0". The exon count files were then inputted into DEXSeq.
To identify over-represented gene ontology terms in the set of DEG and genes containing DEU, we used the Bioconductor package GOseq (1.42.0) (90). For data visualisation, we used regularized logarithm (rlog) transformed counts generated by DESeq2 as input for heatmaps, which were generated using the CRAN package pheatmap (1.0.12): rows were centred; unit variance scaling was applied to rows; and both rows and columns were clustered using Pearson correlation distance and average linkage, where red represents higher expression and blue lower expression on a linear correlation scale. A value of 1 indicates a positive association, while a value of -1 indicates a negative association, and a value of 0 indicates no association. Venn diagrams were generated using the CRAN package VennDiagram (1.6.20).

Development • Accepted manuscript
To identify verified transcription factors amongst DEG and DEU genes we merged the most recent updated lists from the Riken mouse transcription factor database (91) and the mouse transcription factor list from The Animal Transcription Factor DataBase (92) and intersected our gene lists with this merged list.
Canonical correspondence analysis is a multivariate analysis that allows the comparison of experimental transcriptome data with publically available datasets from other laboratories (93). CCA was performed using the CCA function of CRAN package vegan, as previously described (83). The GSE38909 dataset was used as the environmental variable and our dataset was regressed onto it. The GSE38909 dataset contains TCR-transgenic AND DP thymocytes stimulated with a positively selecting peptide (gp250) or a non-selecting control peptide Hb (49). To represent environmental variables of interest, the 2000 most significant DEG (lowest p-values, calculated by moderated eBayes adjusted for false positives) between DP thymocytes stimulated with a non-selecting control peptide and DP thymocytes stimulated with the positively selecting peptide were used, to generate a scale of unstimulated to TCR-signalling-for-positive-selection, and we regressed our datasets onto this axis.  The y-axis represents the biological process GO terms and the x-axis represents the percentage of genes found in that GO term category. Dot size represents the number of genes and the color indicates the p-value.