Mechanisms underlying the direct programming of mouse embryonic fibroblasts to thymic epithelial cells by FOXN1 Free

The thymus is a primary lymphoid organ and is the major source of self-restricted, self-tolerant naïve T cells, which are required for a robust adaptive immune system. However, through aging, the thymus is one of the earliest organs that starts losing its function, also known as thymic involution. Thymic involution leads to decreased immune function, which significantly increases the risk of cancer and auto-immune diseases (Chinn et al., 2012; Hale et al., 2006; Palmer et al., 2018). Thus, finding an effective method with which to rescue the loss of thymus function caused by thymus involution and thymus abnormalities could have significant translational impact.

Thymic epithelial cells (TECs) are the main functional resident cell types in the thymus. TEC-thymocyte interactions are required for T lineage commitment and for all stages of thymocyte differentiation, proliferation, selection and survival (Takahama, 2006). TECs are divided into cortical (cTEC) and medullary (mTEC) lineages, based on their location and function. There is evidence that both lineages derive from a common progenitor, and that mTEC lineages derive from cells that have initiated cTEC differentiation (the ‘cTEC first’ model) (Baik et al., 2013; Bornstein et al., 2018; Rossi et al., 2006; Takahama et al., 2017). The mechanisms inducing cTEC differentiation are not understood. Recent evidence indicates that Notch signaling is required for mTEC lineage commitment but must be suppressed for mTEC differentiation to proceed beyond the progenitor stage (Li et al., 2020; Liu et al., 2020).

Forkhead transcription factor FOXN1 is both necessary and sufficient for TEC differentiation and maturation, and is the key transcription factor required for fetal TEC differentiation. FOXN1 is expressed by all cells specified for TEC identity within the 3rd pharyngeal pouch by mouse embryonic day E11.5 (Vaidya et al., 2016). Null mutation of Foxn1 (nude) in mice, rats and humans disrupts both normal hair growth and thymus development, resulting in immunodeficiency (Blackburn et al., 1996; Nehls et al., 1994). Foxn1 also plays a crucial role in postnatal thymus function and maintenance. Foxn1 is expressed in both cTECs and mTECs during postnatal stages but significantly decreases during aging, which is considered the major cause of thymic involution (Bredenkamp et al., 2014a; Chen et al., 2009; O'Neill et al., 2016; Ortman et al., 2002; Rode et al., 2015; White et al., 2022).

Transplanting neonatal thymus or cytokine induction can partially rescue thymic function after involution, but these methods are all highly limited by resources and efficacy. Another approach is to use cellular programming to generate key cell types in vitro for functional studies or transplantation. Direct lineage programming is characterized by a dynamic conversion of cellular morphology and transcriptomes. We and our collaborators have shown that overexpression of Foxn1 is sufficient to reprogram mouse embryonic fibroblasts (MEFs) into functional induced thymic epithelial cells (iTECs) (Bredenkamp et al., 2014b). These iTECs express FOXN1 downstream targets Dll4, Ccl25 and Kitl, and exhibit epithelial cell-like morphology. However, iTECs cultured alone do not express MHC Class II (MHCII), a crucial functional component of TECs; although, upon co-culture with early thymic progenitors (ETPs), iTECs can support their maturation into single positive (SP) T cells. Finally, iTECs grafted with supporting mesenchyme under the kidney capsule of nude mice can generate an organ with distinct cortical and medullary regions that could generate and export T cells. These results suggest that iTECs placed in a crosstalk environment (upon co-culture and/or transplant) are competent to upregulate MHCII and to mature into both cortical and medullary lineages. They thus suggest the potential for both the generation of an in vitro system to study TEC function and TEC-T cell interactions, and for future therapeutic application of iTECs; however, there are significant challenges to overcome before the in vitro system can have broader utility.

In this study, we used both bulk and single-cell RNA-seq to dissect the Foxn1-induced direct reprogramming process (Treutlein et al., 2016). We sequenced iTEC transcriptomes at different stages of programming to characterize iTECs at the cellular and molecular level. We showed that iTEC generation consisted of early and late programming stages characterized by discrete gene expressions. In the early stage, iTECs were characterized by activation of FOXN1 downstream targets, cell-matrix re-organization and cell morphology changes, and inhibition of cell cycle and cell proliferation. We also demonstrated that releasing the proliferation block with overexpression of Myc improved the programming efficiency. In the late programming stages, our data showed that iTECs not only upregulate additional epithelial markers but also expressed both cTEC and mTEC markers, although only mTEC progenitor and not mature mTEC markers were observed. Our current data provide new evidence demonstrating that these mTEC progenitors are pre-existing in iTEC cultures. We further showed that Notch inhibition using DAPT in iTEC cultures induced mTEC terminal differentiation (albeit at a low level), as measured by the marker Aire. Finally, we found that co-culturing iTECs with CD3^Lo double-positive thymocytes increased MHCII expression in iTECs to similar levels found in fetal TECs. Together, this study provides new insight into the mechanisms behind iTEC programming, which is crucial for developing iTECs as a useful experimental tool both for understanding de novo TEC biology and for potentially translational pre-clinical applications.

To generate samples representing different iTEC programming time-points, we first modified the originally published method of iTEC generation using nucleofection-based direct transfection of a Cre-expressing plasmid, instead of 4-OHT treatment, to activate an inducible CreER transgene (Bredenkamp et al., 2014b). In brief, R26-CAG-STOP-Foxn1-IRES-GFP/+ (Fig. 1A,B; R26-iFoxn1) heterozygous MEFs from E13.5 embryos were transfected with a PGK-Cre plasmid on day 0. After 48 h of culture, about one-third of the MEFs removed the STOP cassette by Cre recombinase and activated the R26-iFoxn1 locus, which was identified by GFP expression (Fig. 1C). We termed these Foxn1- and GFP-expressing cells iFoxn1 MEFs. This method avoids the low level of ‘leaky’ Cre expression from CreER, allowing us to precisely control the onset of iTEC programming and providing a clean starting point for downstream analyses.

As an initial characterization of the programming process, we analyzed the iFoxn1 MEFs for the epithelium-specific marker EpCam and the TEC maturation marker MHCII. Flow cytometry showed that more than 10% of iFoxn1 MEFs express EpCam 13 days after the initial activation of transgenic Foxn1 expression; however, as in the previous report, MHCII protein was not detected (Bredenkamp et al., 2014b) (Fig. 1D). Also consistent with the previous report, many of the iFoxn1 MEFs were keratin 8 positive 10 days after transfection (Fig. 1F). No keratin 8-positive cells were observed in the mock transfected MEF controls (Fig. 1E).

We collected iFoxn1-MEFs at different time points after the initial iFoxn1 activation and used RT-qPCR to analyze the expression of several known FOXN1 target genes (Žuklys et al., 2016) in TECs (Fig. 2A,B). FOXN1 target genes were detected early after induction of Foxn1 expression. Dll4 and Kitl significantly increased as early as 24 h after transfection, and Ccl25 was significantly upregulated 2 days after transfection (Fig. 2A). The increase in these FOXN1 target genes over the course of the 12-day iTEC induction period did not follow a linear pattern (Fig. 2B). Instead, these three genes underwent rapid expression increases, starting at 7 days after transfection, accelerating further after 10 days (Fig. 2B). This pattern suggests that programming has two different stages: an early stage with initial low-level expression of target genes, followed by a later stage in which FOXN1 target gene expression is dramatically accelerated.

To further understand the dynamic process of iTEC programming, we designed a time point-based bulk-RNA-seq experiment, choosing our time points based on the induction of FOXN1 target genes (Fig. 2C). We used mock-transfected R26-iFoxn1/+ MEFs as our control, representing the initial MEFs transcriptome status (no Cre MEFs). We induced R26-iFoxn1/+ MEFs to express Foxn1 using Cre transfection, cultured for 2 days and sorted cells that had successfully induced iFoxn1 expression as GFP⁺ using flow cytometry (as in Fig. 1C). These GFP⁺ iFoxn1 MEFs represented initial programming and are referred to as iTEC2d. We also cultured iTECs for a total of 10 days after transfection to represent the late stage of programming, and refer to these iTECs as iTEC10d. For the final stage, we cultured iTECs for a total of 12 days and used EpCam antibody staining to separate them into two distinct groups: EpCam⁻ and EpCam⁺. Each RNA-seq sample consisted of cells pooled from multiple independent iTEC cultures. RNA-seq was performed in duplicate (MEFs, iTEC2d, EpCam⁻ and EpCam⁺) or triplicate (iTEC10d).

To achieve an initial understanding of the progression of iTEC programming, analysis of MEFs and variance stabilizing transformation (VSD) of the iTEC whole-transcriptome profiles were carried out, followed by principle component analysis (PCA) (Fig. 2D). We found that, even 48 h after the initial activation of ectopic Foxn1, iTEC gene expression was changed significantly when compared with MEFs, as indicated by the clear separation of iTEC2d and ‘no Cre’ MEF samples. The iTEC10d sample was further separated from the iTEC2d sample but did not show a linear increase in variance from ‘no Cre’ to iTEC10d through iTEC2d. This indicated that the iTEC programming process is not continuous at the gene expression profile level and, instead, has distinct gene expression patterns at the early and late stages. Finally, the EpCam⁺ samples were more distinct than EpCam⁻ samples compared with iTEC10d samples. The whole transcriptome analysis showed a linear progression of differential gene expression in the order of 10d, EpCam⁻ and EpCam⁺, which suggests that EpCam⁺ cells represent a more advanced programming state (Fig. 2D). The sample distance count matrix also showed this progression from MEFs to EpCam⁺ iTECs with a significant difference between each step of programming (Fig. 2E). Furthermore, each programming stage is generally reproducible, although two EpCam⁺ samples differed more from each other than did the replicates of earlier stages.

The gene count matrix of differentially expressed genes (DEGs) (adjusted P value <0.01) results further showed that each stage of programming had a distinct signature with sets of genes turning off and on across the timeline (Fig. 2F). Although many genes were gradually upregulated throughout the process, some genes were dynamically expressed. Other genes are downregulated in iTEC2d and iTEC10d samples compared with MEFs but then upregulated further in EpCam⁺ and EpCam⁻ samples. Similar to the global DEG analysis, our assessment of transcription factor changes revealed that specific sets of transcription factors were expressed dynamically at the early and late programming stages. Although some transcription factors start to be expressed in the early stage, many are not expressed until the late stage of programming (Fig. 2G), suggesting a time-dependent gene expression regulation.

Finally, we analyzed upregulated and downregulated genes at different stages by gene enrichment analysis (Fig. S1). We found that in early programming stages from MEFs to iTEC2d, iTECs were characterized by downregulation of cell shape, extracellular matrix (ECM) and other fibroblast-related genes, and upregulation of cytokine signaling and other immune-related categories (Fig. S1A,B). Between 2 days and 10 days, the cell cycle and DNA replication were inhibited, and genes related to epithelial characteristics were upregulated, indicating that the iTECs are transitioning into epithelial cells during this second stage of programming (Fig. S1C, D). The distinctions between EpCam⁻ and EpCam⁺ 12-day iTECs were also made more evident when compared with iTEC10d. EpCam⁻ iTECs contained more programmed cell death and cell cycle inhibition categories, whereas epithelial-related genes were downregulated, suggesting that they constitute cells failing to maintain the transdifferentiation process and may undergo cell death (Fig. S1E,F). In contrast, EpCam⁺ cells upregulated differentiation and metabolism-related genes, whereas genes related to mesenchymal properties were downregulated (Fig. S1G,H).

Together, our results show that Foxn1-induced iTEC generation proceeds through intermediate steps that are significantly different from both MEF controls and the late programming stage. These data indicate that each programming stage is characterized by a unique gene expression profile, rather than a gradual accumulation of gene expression changes that build through the process.

The initial iTEC report showed formation of a well-organized organoid and differentiation of iTECs into both cTECs and mTECs upon transplantation under the kidney capsule (Bredenkamp et al., 2014b). However, it did not characterize cTEC and mTEC lineage markers in detail in iTECs in culture. We performed MA plot analysis to investigate TEC lineage specification and differentiation-related gene expression patterns during the programming process (Fig. 3A-E).

Comparison of MEFs to iTEC 2d (Fig. 3A) showed upregulation of FOXN1 target genes Dll4, Kitl and Ccl25. Strikingly, cTEC markers appeared very early, with Krt18 (and to a lesser extent Krt8) already significantly upregulated in iTEC2d samples. Psmb11 (β5t), a cTEC marker and direct Foxn1 target, was also upregulated but expressed at a lower level than Krt8 and Krt18. In contrast, mTEC markers such as Krt5, Krt14 and Cldn4 were not detected at early stages (Fig. 3A), indicating that cTEC markers are expressed before mTEC markers during iTEC generation.

Between 2 days and 10 days, FOXN1 target genes were further upregulated. More remarkable is very strong upregulation of both cTEC (Krt8 and Psmb11) and mTEC markers (Fig. 3B). mTEC markers first appeared in iTEC10d samples, where they were highly differentially expressed compared with iTEC2d samples but still at much lower expression levels than the cTEC markers. Interestingly, Foxn1 target genes and cTEC markers remained stable in iTEC10d samples and in both EpCam⁻ (Fig. 3C) and EpCam⁺ (Fig. 3D) iTECs , but mTEC-related genes Krt14 and Krt5 were differentially upregulated specifically in EpCam⁺ iTECs at this latest stage of programming (Fig. 3E). Notably, only mTEC markers were upregulated in EpCam⁺ iTECs compared with EpCam⁻ iTECs. This progressive upregulation of TEC differentiation genes with delayed mTEC marker upregulation was further evident in a heat-map analysis (Fig. 3F). However, EpCam⁺ iTECs had much lower expression of these markers than fetal TECs, which could be due to lower expression of Foxn1 in the iTEC system, or because this is a population analysis so far, and the proportions are very low (Fig. 3G). We also detected cTEC markers β5T (translated from Psmb11), Krt8 and DLL4 at the protein level in iTECs (Fig. 3H).

To understand the iTEC programming process while understanding TEC lineage specification, we generated single-cell RNA-sequencing data of MEFs, and iTEC3d and iTEC12d samples. An unsupervised clustering of a batch-corrected merged dataset with Seurat (Hao et al., 2021) yielded nine high-quality cell clusters (Fig. 4A). We detected a Foxn1-high cell cluster located at the middle of the UMAP space that also expressed high levels of TEC markers, and a distinct set of proliferative cells. Despite having sorted iTECs based on their GFP expression (Fig. 1C), many of the iTECs retained MEF-like transcriptomes across these time points, suggesting Foxn1-susceptible and Foxn1-recalcitrant MEFs (Fig. 4B). The Foxn1-high cluster was detected in iTEC3d samples, but only appeared in large numbers in iTEC12d samples. We also detected a loss of the proliferative cell clusters 12 days post-programming, which we will explore later (Fig. 6).

Our RT-qPCR and RNA-sequencing data revealed an exponential increase in the relative expression of Ccl25, Dll4 and Kit-l (Figs 2A,B and 3F). Violin plots showed that this increase correlated with high levels of Foxn1 (Fig. 4C). Atop previously surveyed cTEC genes (Fig. S2), we also detected cells expressing a high level of Prss16, a thymus-specific serine protease (Fig. 4D-F). A small number of cells expressed high Cldn4 within the Foxn1^Hi iTEC cluster in iTEC12d samples (Fig. 4G-I), suggesting that mTEC progenitor (mTEPC) generation also requires a high level of Foxn1. When we surveyed the relationship between a cTEC marker, Krt8, and an mTEPC marker, Cldn4, it exhibited a Pearson correlation coefficient of 0.45, indicating a moderate positive correlation between the two genes. Furthermore, although there are many Krt8⁺Cld4⁻ cells, the majority of Cldn4⁺ cells are also Krt8⁺, suggesting the developmental trajectory of Cldn4⁺ mTECs is consistent with the ‘cTEC first’ model of TEC differentiation (Baik et al., 2013; Bornstein et al., 2018; Rossi et al., 2006; Takahama et al., 2017) (Fig. 4J). Of note, a small number of the iTEC12d samples expressed high levels of Avil, which is a tuft-like mimetic TEC marker (Fig. 4K-M).

Overall, our merged single-cell transcriptomic data validated our bulk-RNA-seq result, showing different stages of iTEC programming, and revealed that iTEC programming did not occur evenly across Cre-induced MEFs. There was a subset of MEFs that were predisposed to undergo transdifferentiation in response to Foxn1 overexpression, whereas other MEFs were recalcitrant. Successfully programmed MEFs yielded iTECs with high levels of archetypal cTEC, mTEC and tuft-like mimetic TEC markers. Only a small number of iTECs expressed high levels of canonical mTEC markers, which our results suggest could be further improved by sorting Foxn1-receptive MEFs before inducing Foxn1 expression.

The differential lineage-specific gene expression analysis showed that the mTEC progenitor marker Cldn4 was one of the highest expressed mTEC markers at later iTEC stages, but the terminal differentiation marker Aire was not detected in our bulk-RNA-seq data (Fig. 3). In our single-cell RNA-seq data of 12-day-old iTECs, only one out of 8290 cells expressed Aire (Fig. S2E). Given these data, we hypothesized mTEC differentiation is induced at the progenitor stage during programming, but further mTEC differentiation is inhibited or incomplete. Recently, we and others showed that Notch signaling plays an essential role in mTEC proliferation and differentiation (Li et al., 2020; Liu et al., 2020). Most crucially, too much Notch signaling after initial mTEC specification blocks mTEC terminal differentiation and mTECs arrest at Cldn4⁺ progenitor states (Goldfarb et al., 2016). These in vivo mechanisms are similar to what we observe in iTEC generation and suggest the possibility that persistent Notch signaling may suppress mTEC differentiation during iTEC programming in vitro.

To test this hypothesis, we first analyzed gene expression in the Notch signaling pathway (Fig. 5A,B). We found that multiple Notch receptors, ligands and downstream target genes were upregulated during early stages and maintained throughout late programming. To test whether this maintenance of Notch signaling blocked mTEC differentiation at the Cldn4⁺ stage, we designed a temporally controlled Notch inhibition experiment (Fig. 5C). We first allowed early programming to progress to establish initial mTEC identity. At day 5 after the initial transfection, we added the Notch signaling inhibitor DAPT daily for 7 days. DAPT-treated cultures reduced but did not eliminate Notch signaling, as indicated by a 50% reduction in Hes1 expression (Fig. 5D). Cldn4 was also downregulated to a similar degree, indicating that iTECs had reduced mTEC progenitor characteristics under DAPT treatment. mTEC differentiation was analyzed by measuring Aire expression, which is a crucial marker for mTEC maturation and terminal differentiation. Aire expression was detected in six out of 11 experiments. In five out of those 6 experiments, Aire expression was significantly upregulated in the DAPT cultures compared with DMSO controls (Fig. 5D). To test whether this Aire expression was sufficient to induce AIRE-dependent promiscuous gene expression (PGE), we tested Ins2, Chrna1 and Tpo expression by qPCR. However, these assays were negative, suggesting that the levels of Aire are too low to efficiently induce PGE, or that additional conditions are required other than Aire expression alone that are not present in the iTEC cultures. As an independent assessment of mTEC differentiation, we also measured Fezf2, an Aire-independent marker of mTEC differentiation. Fezf2 was expressed, but was not consistently different between DMSO- and DAPT-treated cultures (Fig. 5D).

Immunocytochemistry of iTECs12d treated with DMSO (Fig. 5E) or DAPT (Fig. 5F) shows AIRE only in Notch-inhibited iTECs12d samples. AIRE expression in these cells resembled punctate in vivo AIRE expression (Fig. 5G). Inhibiting Notch signaling did not affect iTEC programming efficiency, based on the percentage of the EpCam⁺ population, and it was also insufficient to induce MHCII protein expression (Fig. 5H). Taken together, these data suggest that the upregulation of Aire is a specific response to inhibiting Notch signaling during iTEC programming and that iTEC generation may follow the same mTEC differentiation mechanisms in vivo, which raises the possibility of establishing iTEC as a useful system for investigating mTEC differentiation and function.

FOXN1 promotes TEC proliferation in vivo (Chen et al., 2009; Chojnowski et al., 2014; Nowell et al., 2011); however, we observed the disappearance of proliferative cells in 12-day-old iTECs in vitro (Fig. 4A,B). To investigate this question at the molecular level, we interrogated our RNA-seq data for changes in cell cycle-associated genes during the early and late stages of iTEC programming. Our expression data showed that cell cycle inhibitors, including Rb, Cdkn1a (p21) and Gadd45b were upregulated beginning at the early iTEC2d stage, while cell cycle-promoting genes, such as Myc, Mcm7 and p53 (Trp53), were downregulated in the late stage (Fig. 6A,B).

Cell cycle arrest is common, and often required, in cell programming systems (Treutlein et al., 2016). This proliferation block prevents expansion of iTECs in culture and limits their utility as an experimental system. To test whether we could rescue cell cycle arrest without inhibiting the programming process, we used an inducible Rosa26 transgene for Myc that can be activated by Cre plasmid transfection. This transgene also expresses human CD2 antigen, which allows the detection of MYC-positive cells by huCD2 staining (Calado et al., 2012). We generated double transgenic MEFs (iFoxn1/+;iMYC/+) in which Cre expression during initial iTEC induction should activate both Foxn1 and Myc gene expression (Fig. 6C). Cells that have activated both alleles are GFP⁺ and CD2⁺. Using propidium iodide staining and cell trace analysis, we found that Myc overexpression partially rescued cell cycle arrest in iTECs (Fig. 6D,E). iFoxn1⁺iMyc⁺ iTEC cultures had increased S-phase frequency (Fig. 6D) and a dramatic shift to the left in cell trace analysis (Fig. 6E). This shift indicated dilution of the cell trace reagent, consistent with cells undergoing proliferation. Activation of Myc alone was unable to induce iTEC generation or EpCam expression, as detected by flow cytometry (Fig. 6F). Surprisingly, activation of the cell cycle did not block iTEC programming, but instead it further increased its efficiency, as measured by EpCam expression (Fig. 6G). These results suggest that cell cycle arrest is a side-effect of the programming process, rather than a necessary step, which would enable expansion of iTECs in culture and increase their utility.

We isolated and sequenced in vivo-derived total fetal TECs at E14.5 and compared this data with that of EpCam⁺ iTECs to determine their similarity at the transcriptome level. PCA analysis clearly demonstrated that all iTEC samples are still very different from fetal TECs (Fig. 7A). The sample distance count matrix showed a similar result (Fig. 7B). GO term analysis indicated that fetal TEC overall had higher expression of cell-adhesion molecules relevant to epithelial identity and TEC differentiation (including claudins, CD40 and CD28), and a variety of immune function-related gene categories (Fig. 7C). EpCam⁺ iTECs had higher expression of categories indicating retention of mesenchymal and/or fibroblast-like characteristics (Fig. 7D). MA plot analysis of EpCam⁺ iTEC versus fetal TEC also showed that both mTEC and cTEC markers were higher in fetal TEC, but the degree of difference was less in iTEC10d samples versus iTEC2d samples (Fig. 7E,F).

To investigate the gene expression differences between iTEC and fetal TEC in more detail, we first used a count matrix heatmap (Fig. 7G). This analysis clearly showed that one of the two EpCam⁺ samples was clearly more similar to fetal TECs, indicating a more advanced stage of programming. There were four genes that were not significantly upregulated in any iTEC samples compared with fetal TECs (blue/white in all iTEC samples, red/orange in fetal TECs). Three of these were MHCII genes (annotated as H2-A or H2-E), and the other was Tbata, which regulates TEC proliferation in vivo (Flomerfelt et al., 2010). The significant downregulation of Col6a3 (white in fetal TECs, orange in iTECs) in fetal TECs compared with all iTEC samples is consistent with persistent retention of fibroblast characteristics in iTECs. We also performed heatmap analysis of 18 of the most highly DEGs (Fig. 7H). Ten of these genes associated with MHCII, including Ciita, which encodes the major transcriptional regulator of MHCII expression (León Machado and Steimle, 2021). This analysis further highlighted the lack of MHCII-related genes as a defining feature in in vitro iTECs when compared with fetal TECs. In addition, three genes are associated with TEC function and/or proliferation (Psme1, Prss16 and Tbata), and one with Notch signal modulation (Lfng). This heatmap analysis also demonstrated that the more differentiated EpCam⁺ sample shared some characteristics with fetal TEC, including upregulation of MHCII-related genes (Psme1 and Rfxap) and Ifngr1, and downregulation of Creb1, Rfx1, Nfya and Ctsb. Similar patterns were seen with cell cycle-related genes (Fig. 7I) and Notch pathway-related genes (Fig. 7J). In both cases, one of the two EpCam⁺ samples was more similar to fetal TECs, while still retaining significant differences. These data are consistent with our results above and suggest that these pathways may be potential target candidates for improving iTEC programming to generate cells more comparable with fetal TEC.

Expression of MHCII genes H2-Aa, H2-Ab1 and H2-Eb1 was not detectable with our RNA-sequencing. However, we could detect a low, overall increase in their expression when we surveyed these genes with our single-cell RNA-sequencing dataset. About 2% of iTEC3d expressed low levels of H2-Ab1, which increased to 2.9% at the iTEC12d stage (Fig. 8A-C, Table 1). We calculated the average expression of H2-Aa, H2-Ab1 and H2-Eb1 in each sample, and made a graph with relative values per gene (Fig. 8D). The average expression of H2-Ab1 (2.4-fold), H2-Aa (2-fold) and H2-Eb1 (2-fold) increased in iTEC12d samples compared with iTEC3d samples (Fig. 8D), suggesting an overall increase in MHC class II gene expression in iTEC12d samples.

Although we did not detect MHCII protein on the surface of iTECs after DAPT treatment (Fig. 5H), the increase in mTEC differentiation seen with DAPT could impact MHCII at the gene expression level. To assess the effects of Notch signaling on MHCII gene expression, we performed single-cell RNA-sequencing after treating iTEC12d samples with DAPT and measured expression changes in H2-Aa, H2-Ab1 and H2-Eb1 (Fig. 8E-K). H2-Ab1 expression, which exhibited the largest increase in iTEC12d samples compared with iTEC3d samples, did not change when Notch signaling was inhibited with DAPT (Fig. 8E-G,K). However, Notch-inhibited iTEC12d samples compared to iTEC12d control samples showed 2.3- and 2.6-fold increases in H2-Aa and H2-Eb1 expression, respectively (Fig. 8H-K).

iTECs can support early T lineage progenitor (ETP) differentiation, which requires MHCII expression (Bredenkamp et al., 2014b). However, MHCII gene expression was not evaluated directly in those experiments. We designed a co-culture experiment to determine whether TEC-thymocyte crosstalk can successfully induce expression of MHCII. At the beginning of late-stage iTEC programming (5 days after transfection), we co-cultured 100,000 iTECs with 80,000 E16.5 double-positive thymocytes for 5 consecutive days before collecting them for single-cell RNA-seq. Co-culturing with fetal thymocytes increased MHCII gene expression: around two to three times more cells are expressing MHCII genes at approximately three times higher levels compared with DAPT controls (Fig. 8L-N). Next, we compared co-cultured and DAPT-treated iTECs to MEF and iTEC10d samples by RT-qPCR (Fig. 8O-Q). CD80 expression did not differ between any iTEC samples (Fig. 8O). However, we found a significant increase in the relative H2-Aa expression in Notch-inhibited iTEC10d samples compared with all controls (Fig. 8P), supporting our previous observation with the single-cell RNA-seq data (Fig. 8H-K). Co-culturing with E16.5 DP thymocyte significantly increased the H2-Aa expression level even further when compared with all samples (Fig. 8Q), showing that crosstalk with DP thymocytes is sufficient to induce MHCII gene expression. Combined, these data provide a stepping stone towards further understanding and dissecting the molecular mechanism behind MHCII gene expression, which will improve the iTEC in vitro platform.

The previous iTEC study demonstrated that expression of Foxn1 in MEFs induced a TEC-like differentiation state, with upregulation of FOXN1 direct targets and induction of epithelial characteristics (Bredenkamp et al., 2014b). That study also showed that co-culture with fetal ETP or transplant under the kidney capsule could induce further differentiation and support the T-cell development. These results indicate that the iTEC system could represent a novel platform for studying FOXN1 function and TEC-thymocyte interactions, and could present a pathway to therapeutic thymus transplant. The current study builds on the initial report, and provides a significant advancement of our understanding of iTEC programming and opens avenues for future improvement of this technique. We showed that iTEC programming is not linear, and defined characteristics of early and late stages. Although we cannot fully rule out the possibility that other factors are necessary for the transdifferentiation, our data indicate that programming begins with immediate upregulation of FOXN1 direct target genes that are not normally expressed in MEFs, consistent with the role of Forkhead transcription factors as pioneer factors (Zaret and Mango, 2016). During later stages, we observed strong upregulation of known FOXN1 targets, as well as gene expression and morphological changes not directly attributable to known FOXN1 in vivo functions. Our data are also generally consistent with iTEC programming following a path consistent with the ‘cTEC first’ model, including arrest at the mTEC progenitor stage (Baik et al., 2013; Bornstein et al., 2018; Rossi et al., 2006; Takahama et al., 2017). We rescued this arrest by inhibiting the Notch signaling pathway, which suggests that mTEC differentiation in iTECs follows similar pathways to those of endogenous fetal mTECs. Additionally, we improved the efficiency of iTEC programming by rescuing cell cycle arrest. Finally, we identified crucial gene expression differences between iTECs and fetal TECs. Notably, expression of MHC-related genes was severely downregulated in iTECs, but we demonstrated this can be significantly improved by co-culturing iTECs with fetal DP thymocytes.

One persistent question is why iTEC differentiation is not uniform across iTEC cultures. In our current iTEC direct programming protocol, the EpCam⁺ population is always a subset of the cultures, even though all iTECs are expressing the same Foxn1 transgene for the same amount of time. One possibility is that the cell cycle plays a partial role in this process, as rescuing the cell cycle block increased the percentage of Epcam⁺ iTECs (although it still was not 100%). Moreover, MEFs are a heterogeneous population (Shakiba et al., 2019), and our single-cell RNA-seq suggested that some MEFs are more susceptible to Foxn1 programming than others. Another variable is endogenous Foxn1 expression. We consistently detect activation at the endogenous Foxn1 locus at a low level during the late stage of programming [transgenic Foxn1 levels were 2037-fold higher than endogenous Foxn1; ±193 (s.e.m.), n=6], although it is unclear whether this occurs in all cells. Although this additional low level of Foxn1 expression is unlikely to have a significant impact on programming, it could indicate the degree to which programming is occurring at the epigenetic level. Additional analyses of heterogeneity during this process needs to be carried out to test these hypotheses.

We detected both cTEC and mTEC markers during iTEC programming, but mTECs were limited to the Cldn4⁺ progenitor stage. Although we demonstrated that blocking Notch signaling increased mTEC differentiation, it is unclear how cTEC and mTEC differentiation occurs at the individual cell level. Following the ‘cTEC first’ model, it is possible that all iTECs express cTEC markers before subsets of more differentiated cells start to express mTEC markers later in programming. Alternatively, it is possible that different subsets of cells initiate cTEC and mTEC differentiation directly. Our single-cell RNA-seq data suggest that mTECs were rising from cTECs due to the co-expression of cTEC markers with Cldn4, but further studies are needed to test this model in the iTEC system.

Several previous studies show that cell cycle arrest is necessary for direct lineage reprogramming of some cell types, such as neuron direct reprograming (Treutlein et al., 2016). Although cell cycle arrest occurs during iTEC direct programming, our data show that this arrest is not necessary. In contrast, removing this cell cycle block by forced Myc expression results in higher efficiency of EpCam⁺ iTEC generation. It is still unclear what role the cell cycle has during iTEC programming and how proliferating iTECs compare with cell cycle-arrested iTECs. It is possible that Myc overexpression has a larger impact on iTEC programming than on cell cycle regulation, as Myc expression both alone and as a component of the retinoblastoma pathway is implicated in TEC differentiation in vivo (Cowan et al., 2019; Garfin et al., 2013). However, we found that Myc expression alone had no detected impact in the absence of Foxn1 expression. Alternatively, this could be due to a technical artifact whereby Cre-induced gene damage causes reduced cell proliferation (Loonstra et al., 2001), which could be addressed with a control in which Cre is expressed without Foxn1 induction. Further experiments are needed to analyze the impact of rescuing cell cycle arrest and to further understand the mechanisms underlying the cell cycle-mediated effect on iTEC programming.

iTEC cultures alone do not significantly express MHC Class II, which is a crucial functional hallmark of TEC differentiation. In the original iTEC report, co-culture or transplant of iTEC supported the development of single positive CD4⁺ T cells (Bredenkamp et al., 2014b). These data suggest that iTEC-thymocyte crosstalk induced MHCII gene expression, although the impact of crosstalk on iTEC phenotypes was not directly assessed. Our single-cell RNA-seq datasets suggest that Foxn1 overexpression is sufficient to induce MHCII gene expression at low efficiency. Intriguingly, Notch-inhibited iTECs had increased expression of H2-Aa and H2-Eb1, suggesting a potential role of the Notch signaling pathway in MHCII expression. Co-culturing with fetal DP thymocytes was sufficient to increase the relative expression of the MHCII genes H2-Aa, H2-Ab1 and H2-Eb1. Thus, future experiments to identify the crosstalk elements that increase MHCII expression will be essential to design a fully viable iTEC culture system that can be used as a novel experimental platform for understanding TEC biology and function, and as an approach to generate human iTECs for in vitro or in vivo therapeutic uses.

Rosa26^{CAG−STOP−Foxn1−IRES−GFP/+} mice (Bredenkamp et al., 2014b) were backcrossed onto the C57BL/6J genetic background for at least five generations then maintained by intercrossing. Each batch of mouse embryonic fibroblasts (iFoxn1+ MEFs) was pooled from embryos in a single litter of E13.5 embryos generated from a Rosa26^{CAG−STOP−Foxn1−IRES−GFP/+} female crossed to a C57Bl/6J male. For timed matings, noon of the day of the vaginal plug was taken as day 0.5. All experiments with animals were approved by the UGA Institutional Animal Care and Use Committee.

MEFs were prepared from decapitated and eviscerated E13.5 embryos. The embryos were then finely minced and trypsinized into a single-cell suspension. Cells were plated in DMEM containing 10% fetal calf serum, 2 mM sodium pyruvate, 4 mM glutamine, 50 µg/ml streptomycin and 50 U/ml penicillin (DMEM/FCS). Each embryo was genotyped using the following primers to detect the iFoxn1 allele: iFoxn1 forward primer, 5′-TGG AGT AGG CGG GGA GAA GG-3′; iFoxn1 forward primer mid, 5′-TCG CCC TTC CCA ACA GTT GC-3′; iFoxn1 reverse primer, 5′-GCC CAC ACA CCA GGT TAG CC-3′. MEFs were freshly isolated and initially cultured in 60 mm plates until reaching 80% confluency and transferred to T-75 flasks until reaching ∼80-90% confluency. MEFs were then trypsinized and transferred into cryovials containing freezing medium (MEF medium with 10% DMSO). The cryovials in the isopropanol freezing container (to achieve −1°C/minute) were placed at −20°C for 1 h, then at −80°C overnight. Cryovials were then transferred into liquid nitrogen tanks until use. Once MEFs were thawed, they were immediately transferred to pre-warmed (37°C) MEF medium. MEF medium was replaced with fresh MEF medium the following day to remove residual DMSO. MEFs were expanded once or twice more for downstream experiments.

Primary iFoxn1/+ MEFs were thawed then cultured in DMEM containing 15% fetal calf serum, 2 mM sodium pyruvate, 4 mM glutamine, 50 µg/ml streptomycin and 50 U/ml penicillin (DMEM/FCS). Cells were passaged once at 70%∼80% confluency, then collected at 70% confluency for nucleofection. Primary iFoxn1/+ MEFs for the nucleofection experiment were passaged no more than five times for all analyses in this research. Nucleofection experiments were performed following the manufacturer's protocol for mouse embryonic fibroblasts: free program N-024 (Amaxa Nucleofector, Lonza). To induce iFoxn1 expression, pPGK-Cre-bpA (Addgene 11543) was used in the transfection according to the manufacturer's protocol. After nucleofection, cells were further cultured for 48 h until sorted by flow cytometry.

Adherent MEFs or iTEC2d samples were collected by trypsinization using 0.05% trypsin for 2 min at 37°C. iTEC10d and iTEC12d samples were collected using 0.25% trypsin for 15 min at 37°C. For fetal TEC collection, E15.5 fetal stage thymi were dissected and digested in 1 mg/ml collagenase/dispase (Roche), and passed through a 100 μm mesh to remove debris. Thymi were processed individually before genotyping, then pooled. PE-Cy7 conjugated anti-CD45 (BioLegend, 30-F11, 1:150) and APC-conjugated anti-EpCam (BioLegend, G8.8, 1:150) were used to isolate TEC populations. Anti-MHCII (BioLegend, 1:150) was also used for TEC analysis. Cell sorting was performed using a BioRad-S3 cell sorter.

Exactly 60 h (5 days) after transfection, 120,000 iTECs were treated with no treatment, DMSO control or 10 µM of the γ-secretase inhibitor DAPT (GSI-IX) (Apex Bio, A8200). The treatments were continued once every 24 h for 7 consecutive days. Two or three technical replicates and at least five biological replicates were carried out for analyzing the effect of Notch signal inhibition on iTECs. Cells were then collected as described in the RNA extraction section.

Exactly 60 h (5 days) after transfection, iTECs were treated with no treatment control, DMSO control or 10 µM of the γ-secretase inhibitor DAPT (GSI-IX) (Apex Bio, A8200), or co-cultured with E16.5 double-positive thymocytes. MEF medium, DMSO vehicle medium and DAPT medium were replenished every 24 h for 5 consecutive days before downstream analyses. Owing to the non-adhering nature of DP thymocytes, small quantities of MEF medium were added every day without removing the existing MEF medium, restricting the duration of the experiment.

RNA was prepared using the RNeasy mini kit (Qiagen) according to the manufacturer's instructions. All samples were DNase treated during the preparation. cDNA was synthesized from concentration-matched RNA using the iScript Reverse Transciption Supermix for RT-qPCR (Bio-Rad) according to the manufacturer's instructions. Thermal cycling conditions for the cDNA synthesis were as follows: 25°C for 5 min, 46°C for 20 min and 95°C for 1 min. Quantitative PCR was performed with Applied Biosystems 7500 Real-Time PCR System (ThermoFisher Scientific). Thermal cycling conditions were as follows: 50°C for 2 min, 95°C for 10 min, 40 cycles of 95°C for 15 s, then 60°C for 1 min. Relative expression was determined after normalization to the geometric mean of two reference genes (Gapdh and Tbp) and after normalization to Gapdh whenever Tbp was undetected. All samples were run with two or three technical triplicates and with no reverse transcribed controls. The primers used for RT-qPCR were purchased from ThermoFisher Scientific (TaqMan Gene Expression Assays). All probes were designed with FAM-MGB Reporter-Quencher.

Post-sorting, iTECs were placed on the gelatin-coated glass coverslips (Electron Microscopy Science, 72293-01) placed in wells of a 24-well plate. On the day of immunocytochemistry, iTECs were washed in sterile Ca²⁺- and Mg²⁺-free 1×PBS. iTECs were then fixed in 4% PFA (EM grade) (Electron Microscopy Sciences, 157-4) for 15 min at room-temperature. Post-incubation, cells were washed in sterileCa²⁺- and Mg²⁺-free 1×PBS for 5 min, three times at room (washing condition). To permeabilize, iTECs were incubated in 1×PBS (with 0.1% Triton X-100) for 15 min at room temperature. Post-washing, iTECs were incubated in blocking solution (10% donkey serum) for 45 min at room temperature. iTECs were then incubated in primary antibodies at 4°C overnight. Post-washing, iTECs were incubated in secondary antibodies conjugated with fluorophores for 1 h at room temperature, then cells were washed and incubated in DAPI for 5 min before being placed on the slides with fluorogel with TRIS buffer (Electron Microscopy Sciences, 17985-10). Primary antibodies used were as follows: AIRE monoclonal antibody (5H12), eBioscience (Invitrogen, 14-5934-82), anti-mDLL4 (R&D Systems, AF1389) and MHC class II (ER-TR3) (Santa Cruz Biotechnology, sc-59318). Secondary antibodies used were as follows: donkey anti-goat IgG (H+L), Alexa Fluor 488 (Invitrogen, A-11055), donkey anti-rat IgG (H+L), Alexa Fluor 594 (Invitrogen, A-21209), donkey anti-rabbit IgG (H+L) and Alexa Fluor 647 (Invitrogen, A-31573).

Bulk RNA-sequencing was performed in the Georgia Genomics and Bioinformatics Core (GGBC) using an Illumina NextSeq 500 sequencing platform to generate single end sequencing data. Raw reads were pre-processed with sequencing grooming tools FASTQC and Trimmomatic, and then further assembled by Hisat2, SAMtools and StringTie. Differential gene expression analysis between samples and graphical display (including count matrix plot, PCA, heatmap and other plots) were performed using DEseq2 and R package ggplot2. Gene ontology enrichment analysis was performed using online EnrichR program (Chen et al., 2013; Kuleshov et al., 2016).

Viabilities of the iTECs (three wells/24-well plate) were confirmed with a hemocytometer with cells stained with Trypan Blue. Libraries were generated with cells with >97% viability. 10,000 iTECs with Single Cell 3′ Reagent Kit (v3) were loaded into a 10X Chromium instrument. The quantity and quality of cDNA were examined before the final DNA library was sequenced on Illumina NextSeq2000 P3 with 200 cycles.

As per the 10X guideline, we used CellRanger (v7) to generate and pre-process FASTQ files. CellRanger generated filtered_feature_bc_matrix files underwent QC steps with Seurat (v.4.3.0) (Hao et al., 2021) in R (v.4.2.1). Individual datasets were merged and batch corrected with top 3000 features. The merged datasets were then analyzed with standard Seurat workflow.

We thank Clare Blackburn for helpful advice and discussions during the progress of this project. We thank J. Nelson in the Center for Tropical and Emerging Global Diseases Flow Cytometry Facility at the University of Georgia for flow cytometry and cell sorting technical support. This work was supported by the Biomedical Microscopy Core and the Georgia Genomics and Bioinformatics Core facilities.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.202730.reviewer-comments.pdf