Research ArticleIMMUNOLOGY

Retracing the evolutionary emergence of thymopoiesis

See allHide authors and affiliations

Science Advances  27 Nov 2020:
Vol. 6, no. 48, eabd9585
DOI: 10.1126/sciadv.abd9585


The onset of lymphocyte development in the vertebrate primordial thymus, about 500 million years ago, represents one of the foundational events of the emerging adaptive immune system. Here, we retrace the evolutionary trajectory of thymopoiesis, from early vertebrates to mammals, guided by members of the Foxn1/4 transcription factor gene family, which direct the differentiation of the thymic microenvironment. Molecular engineering in transgenic mice recapitulated a gene duplication event, exon replacements, and altered expression patterns. These changes predictably modified the lymphopoietic characteristics of the thymus, identifying molecular features contributing to conversion of a primordial bipotent lymphoid organ to a tissue specializing in T cell development. The phylogenetic reconstruction associates increasing efficiency of T cell generation with diminishing B cell–generating capacity of the thymus during jawed vertebrate evolution.


The adaptive immune system arguably represents one of the major evolutionary novelties that distinguish vertebrates from their nonvertebrate ancestors. About half a billion years ago, specialized cell types, such as lymphocytes, and new organs, such as the thymus, emerged, contributing to the radical redesign of animal immunity (13). These morphological and cytological innovations likely occurred in the ancestor common to all vertebrates, since they are present in both of the two extant groups of vertebrates, the jawless and jawed vertebrates (3). Additional novelties, such as the chemokine and chemokine receptor systems, enabled new types of cellular interactions, for instance, between hematopoietic progenitor cells and the emerging thymopoietic tissue (4). The thymic microenvironment provides distinct molecular cues such as Notch ligands, chemokines, and cytokines that interact in a synergistic, context-dependent, and hierarchical manner and, in this way, determine the outcome of hematopoietic precursor cell differentiation (5). The key importance of the thymic microenvironment for T cell development was first revealed through studies of rodents homozygous for mutations at the nude locus (68), which disrupt thymus development, causing animals to be immune deficient. The nude locus encodes a transcription factor of the forkhead class (9), now designated as Foxn1. Subsequent studies in mice indicated that the Foxn1 gene is dispensable for the formation of the thymic anlage during embryonic development (8) but is required for the subsequent steps of differentiation of primitive precursor cells into the typical cortical and medullary subsets of the thymic epithelium (10). Foxn1-expressing thymic epithelial precursors (11) are bipotent (12, 13), with each cell able to give rise to a self-organizing thymopoietic unit containing cortical and medullary compartments capable of supporting T cell development (12). Thus, in the absence of Foxn1, the master regulator of the thymic epithelium in mammals, thymus differentiation is aborted, and T cell development fails (9). Attempts to identify the key functional elements of the thymic epithelium through transgenic expression of candidate genes has highlighted the roles of the chemokines Cxcl12 and Ccl25, the cytokine/stem cell factor Scf, and the Notch ligand Dll4 (5), which direct the early steps of T cell development in the thymic microenvironment. However, the individual components of the genetic network downstream of Foxn1 that regulates the emergence of distinct subsets of thymic epithelial cells (TECs) are unknown.

So far, it has not been possible to reconstruct the nature and evolutionary sequence of the steps that gave rise to a functional thymus, as no extant species representing intermediate phylogenetic stages along the transition from innate to adaptive immunity are available for analysis. To begin to address this problem, we have developed an in vivo reconstruction strategy that rests on the exchange of the mammalian version of the Foxn1 transcription factor with its evolutionarily distant relatives of nonvertebrate and vertebrate origin. Since the decisions controlling alternative differentiation pathways in the hematopoietic system are strongly influenced by extrinsic cues (5, 1417), we hypothesized that it should be possible to directly examine and compare the thymopoietic properties of distinct types of epithelial microenvironments in the thymus as a result of the activities of different Foxn1/4 family members.

Foxn1 and its paralog Foxn4 comprise a distinct two-member family of vertebrate wing-helix transcription factors that recognize a unique (G+C)-rich DNA target sequence, distinguished by the core tetranucleotide sequence, 5′-ACGC (18, 19). This sequence is also recognized by the DNA binding domain (fig. S1A) encoded by the Foxn4 gene of amphioxus (19), indicating that the evolutionary conservation of the protein sequence of the DNA binding domains encoded by this gene family (4) is mirrored in identical target recognition sequences. In addition to the centrally located DNA binding domain, Foxn1/4 transcription factors exhibit N-terminal domains of variable lengths and acidic transcriptional activation domains (20) in their C-terminal region (fig. S1A); of note, the activation domains are functionally interchangeable, as indicated by the fact that the C terminus of the amphioxus (lancelet) Foxn4 protein can replace that of the mouse Foxn1 protein in in vitro assays of transcriptional activation (21). These observations suggest that the N-terminal regions of Foxn1/4 proteins have been important targets for evolutionary modifications, possibly related to changes in function.

To determine the thymopoietic properties of different members of the Foxn1/4 transcription factor family in the context of the mouse hematopoietic system, we have developed a generic transgenic replacement strategy. It begins with Foxn1-deficient mice, in which the thymic epithelium is functionally inactive (10); it is important to note that, because immature progenitor cells persist in the organ anlage, a functional thymus can be formed upon reactivation of the Foxn1 function (12). We then introduce a transgenic construct that contains all regulatory sequences of the mouse Foxn1 gene (fig. S1B) that are necessary to direct tissue-specific and orthotopic expression of any complementary DNA (cDNA) to the endogenous Foxn1 expression domains into this Foxn1-deficient background. A previous study using this system demonstrated that the expression of the cognate mouse Foxn1 cDNA rescued the pleiotropic nude phenotype in the thymus and the skin (22), thus functionally validating the replacement strategy in vivo. In the first application of this method, we showed that the mouse Foxn4 transcription factor gene [which is paralogous to mouse Foxn1 (4), although not expressed in the mouse thymic epithelium (22)] is nonetheless capable of supporting a degree of lymphopoiesis in the reconstructed thymi (22). These results encouraged us to extend our studies to Foxn1/4 family members (fig. S1C) identifiable in extant representatives of evolutionarily more ancient chordates, using the phenotypes of mouse Foxn1 and mouse Foxn4 replacements as references.


Phylogeny of the Foxn1/4 gene family

For centuries, biologists have debated the evolutionary origin of vertebrates (23); the current scenario of chordate taxa suggests that lancelets are the sister group to tunicates and vertebrates (24, 25). This phylogenetic relationship is mirrored in our analysis of the chordate Foxn1/4 gene family, which suggests that an ancestral metazoan Foxn4 gene gave rise to the Foxn1 and Foxn4 genes of vertebrates (Fig. 1 and fig. S2) (4). Although tunicates are considered to be phylogenetically closest to vertebrates (24, 25), we chose the Foxn4 gene of amphioxus for our functional analyses, because tunicate development is substantially secondarily modified (26). Hence, to understand the role of Foxn4 in the emergence of specific aspects of the vertebrate body plan, lancelets appeared to be better proxies than tunicates. Moreover, we had previously shown that the amphioxus Foxn4 gene is expressed in the pharyngeal endoderm, the future site of the thymic epithelium in vertebrates (4). For the present experiments, we focused on one species of lancelets (Branchiostoma lanceolatum). To study vertebrate-specific aspects of the Foxn1/4 gene family, we turned to an extant representative of the most ancient group of jawed vertebrates, cartilaginous fishes. To this end, we identified and isolated the Foxn1 and Foxn4 genes of elephant shark (Callorhinchus milii) for our replacement studies.

Fig. 1 Phylogeny of representative Foxn1 and Foxn4 proteins.

Vertebrate Foxn1 and Foxn4 clades recapitulate the known phylogenetic relationships of vertebrates (representatives of mammals, monotremes/marsupials, birds, reptiles, amphibians, and bony and cartilaginous fishes are depicted). The Foxn1 proteins in vertebrates form a monophyletic clade. The Foxn4 proteins are paraphyletic in vertebrates and in tunicates. The Foxn4 proteins in lancelets form the base of the tree. The support from 1000 bootstrap replicates is shown as color-coded branches. The vertebrate Foxn1 and Foxn4 clades are boxed in different colors. The scale is in units of average amino acid substitutions per site.

Analytical steps in the characterization of thymopoietic function

To reconstruct the functional changes that occurred along the evolutionary trajectory of the thymic microenvironment in vertebrates, we examined the developmental fate of mouse hematopoietic progenitors in the different types of thymic microenvironments. Since the Foxn1-deficient thymic epithelium fails to differentiate and does not support lymphoid development (fig. S3) (69), any lymphopoietic activity in the reconstituted thymic microenvironment must be driven by the expression of the respective Foxn1/4 family member under study. The transgenic thymi were examined both for their lymphopoietic capacity and the characteristics of the epithelial microenvironment. Whereas the hematopoietic compartment was analyzed by cell surface markers and in situ analysis of tissue sections, TECs were additionally characterized by RNA sequencing (RNA-seq). Guided by our previous work (22), we focused our attention on three major phenotypic aspects of the thymi that were reconstituted by the Foxn1/4 gene family members of amphioxus and shark. With regards to the composition of the thymic microenvironment, we specifically addressed the question of whether they are capable of supporting the formation of a distinct medullary area, a key compartment associated with the selection of a self-tolerant T cell repertoire (27). With respect to the lymphopoietic properties of the reconstructed thymi, we paid particular attention to their capacity to support the differentiation of thymocytes throughout the known developmental trajectory (17); moreover, we scored the presence and localization of immature and mature B cells, since our previous studies indicated a propensity of the mouse Foxn4 gene to support B cell poiesis when expressed in the thymic epithelium (22), a capacity that the wild-type mouse thymus lacks.

Lymphopoietic properties of amphioxus Foxn4

We first examined the thymopoietic capacity of lancelet Foxn4 (Bl_Foxn4), the sole family member of the Foxn1/4 gene family in the genome of the cephalochordate B. lanceolatum. Lancelets lack an adaptive immune system, although cytological evidence for the presence of lymphocyte-like cells has been reported in some species (28, 29). Bl_Foxn4 exhibits some thymopoietic activity. After replacement of mouse Foxn1 (Mus musculus, hereafter Mm_Foxn1) with Bl_Foxn4, the transgenic thymic microenvironment predominantly harbors CD4+CD8+ double-positive (DP) immature thymocytes (Fig. 2A), although their number amounts to only about 1% of that in wild-type thymi (Fig. 2B). Two other features distinguish the hematopoietic compartment of Bl_Foxn4 thymi from the corresponding wild-type situation. First, very few, if any, single-positive T cells are detectable (Fig. 2A), and second, the CD4CD8 double-negative compartment is much larger than in wild-type mice (Fig. 2A). Notably, the absolute number of CD45+CD4CD8CD19+B220+IgM+CD93 mature B cells present in Bl_Foxn4 thymi (Fig. 2C) is only 10-fold lower than that in wild-type thymi (Fig. 2D), indicating that, relative to a mouse wild-type thymus, the transgenic microenvironment tends to be more favorable for mature B cells than for immature T cells.

Fig. 2 Thymopoietic capacity of the invertebrate lancelet Foxn4 gene.

(A) Flow cytometric analysis of CD45+ thymocytes, stained for T cell markers; wt (n = 6) and Bl_Foxn4 (n = 7). (B) Absolute numbers of total thymocytes and CD4+CD8+ DP T cells in wt (n = 7) and Bl_Foxn4 transgenic thymi (n = 6); ***P < 0.001; two-tailed t test for both groups. (C) Flow cytometric analysis of CD45+ thymocytes and B cell markers; wt (n = 6) and Bl_Foxn4 (n = 7). (D) Absolute numbers of immature CD45+CD4CD8CD19+IgMCD93+ B cells; P = 0.1058, two-tailed t test; wt (n = 6) and Bl_Foxn4 (n = 7). (E) Ly51 expression and UEA1 binding on EpCAM+CD45 TECs; wt (n = 5) and Bl_Foxn4 (n = 5). (F) Heatmap of differential gene expression patterns of TECs. Genes whose expression is associated with particular TEC subsets are indicated (green, cTEC-like; red, mTEC-like; blue, mature cTEC); analysis based on 18,808 protein-coding genes. Scale refers to the percentage of maximum and minimum values of transcript counts of individual genes. (G) Total numbers of CD45EpCAM+ TECs; wt (n = 5) and Bl_Foxn4 (n = 8). (H and I) T and B cell poietic indices calculated from ratios of means (±SDs correspond to propagated errors); ***P < 0.001 two tailed t test.

The unique lymphoid signature in the Bl_Foxn4-reconstituted thymi is the result of a microenvironment characterized by epithelial cysts (fig. S4, A and B), dominated by cells reminiscent of cortical TECs (cTECs), as determined by cell surface phenotype (EpCAM+Ly51+UEA1) (Fig. 2E) and the keratin signature (K8+K5) in tissue sections (fig. S4A); we note a conspicuous lack of medullary TECs (mTECs) both in flow cytometric (EpCAM+Ly51UEA1+) (Fig. 2E) and immunohistological (K8K5+) (fig. S4A) analyses. Accordingly, the epithelium has a global transcriptional profile heavily skewed toward a cTEC signature (high levels of Bmp4, Psmb11, Cxcl12, and Dll4 genes), with low levels of genes typically associated with mTECs (Trpm5 and Aire) (30); of note, the near-complete absence of mature TECs is also reflected in low levels of Mhc2 gene expression (Fig. 2F). In the wild-type thymus, the Aire gene is expressed in the CD80+MHCIIhi subset of mature mTEC and involved in the regulation of promiscuous gene expression that is required for proper negative selection of autoreactive T cells (31).

Overall, this constellation gives rise to a low thymopoietic index (calculated as the ratio of hematopoietic cells to TECs) for immature CD4+CD8+ T cells; the capacity of Bl_Foxn4 thymi to generate these immature T cells is two orders of magnitude lower than that of wild-type thymi, compatible with the immature status of the transgenic TECs. Because immature CD45+CD4CD8CD19+B220+IgMCD93+ (henceforth IgMCD93+) B cells are barely detectable in the transgenic thymi (Fig. 2, G to I), the reduction of B poietic capacity can be less confidently determined but appears to be at least 10-fold. Collectively, mice expressing Bl_Foxn4 instead of Mm_Foxn1 in TECs fail to support proper T cell development in the thymus, lack appreciable B cell–generative capacity, and are most likely incapable of executing the necessary positive and negative selection steps that normally give rise to a self-tolerant repertoire of mature T cells. As a consequence, transgenic mice rapidly (fig. S4C) succumb to an inflammatory syndrome, chiefly affecting the intestine (fig. S4D), occasionally accompanied by vitiligo (fig. S4E). Nonetheless, our data indicate that Bl_Foxn4 is preadapted to support the early stages of T lymphopoiesis when expressed in the context of a vertebrate immune system, attesting to the strong evolutionary conservation of the Notch signaling pathway in regulating the differentiation of hematopoietic precursors [reviewed in (14)].

Lymphopoietic properties of shark Foxn4

To assess the functional changes associated with the emergence of a vertebrate form of the Foxn4 gene, we next examined the lymphopoietic capacity of the Foxn4 gene of an extant representative of an evolutionarily ancient group of jawed vertebrates (Fig. 1 and figs. S1 and S2). When the thymic epithelium expresses the Foxn4 gene of the cartilaginous fish C. milii (Cm_Foxn4) instead of Mm_Foxn1, the lymphoid compartment (Fig. 3, A to E) features a large fraction of IgMCD93+ immature B cells, which amounts to about 10% of all hematopoietic cells. This represents an increase of three orders of magnitude when compared to the fraction of immature B cells among the thymocytes of wild-type mice, which averages about 0.01% of all thymocytes. In absolute terms, the Cm_Foxn4 thymi harbor only about 1% of the hematopoietic cell numbers in wild-type mice. Despite the comparatively low number of T cells in the transgenic thymi, their differentiation appears to proceed normally, as reflected by the presence of DP and single-positive thymocytes and the absence of systemic inflammation. These results suggest that the Cm_Foxn4 thymic microenvironment is conducive to both T and B cell differentiation (Fig. 3, F to H); this type of unusual lymphoid bipotency was previously observed in Mm_Foxn4 transgenic thymi (22).

Fig. 3 Lymphopoietic activities of C. milii Foxn4 and Foxn1 genes.

(A to C) Total numbers of CD45+ hematopoietic cells {wt (n = 20), Cm_Foxn4 (n = 7), Cm_Foxn1 (n = 11), and Cm_Foxn1/Cm_Foxn4 double transgenics [Cm_dtg (n = 5)]}, DP T cells [wt (n = 20), Cm_Foxn4 (n = 7), Cm_Foxn1 (n = 11), and Cm_dtg (n = 5)], and CD93+ immature B cells [wt (n = 13), Cm_Foxn4 (n = 7), Cm_Foxn1 (n = 7), and Cm_dtg (n = 5)] in thymi of mice with the indicated genotypes. (D and E) Representative flow cytometric profiles for mice summarized in (A to C). (D) wt (n = 16), Cm_Foxn4 (n = 7), Cm_Foxn1 (n = 9), Cm_dtg (n = 5). (E) wt (n = 13), Cm_Foxn4 (n = 7), Cm_Foxn1 (n = 7), and Cm_dtg (n = 5). (F) Total numbers of CD45EpCAM+ TECs; wt (n = 17), Cm_Foxn4 (n = 6), Cm_Foxn1 (n = 11), and Cm_dtg (n = 5). (G) T cell poietic indices; wt (n = 20), Cm_Foxn4 (n = 6), Cm_Foxn1 (n = 11), and Cm_dtg (n = 5). (H) B cell poietic indices; wt (n = 13), Cm_Foxn4 (n = 6), Cm_Foxn1 (n = 7), and Cm_dtg (n = 5). In (A) to (C) and (F) to (H), each data point represents one mouse. ***P < 0.001; one-way ANOVA with Tukey’s multiple comparison test in (A) to (C) and (F) to (H).

Compared to the situation in Bl_Foxn4 transgenic mice, the epithelial microenvironment of Cm_Foxn4 thymi exhibits the cell surface phenotype (Fig. 4A) and the keratin (Fig. 4B) and gene (Fig. 4C) expression patterns of a more mature cTEC compartment; this is particularly evident from the much higher expression levels of Prss16 and Ccl25 genes that are elevated to wild-type levels (Fig. 4C).

Fig. 4 Lymphopoietic activities of C. milii Foxn4 and Foxn1 genes.

(A) Representative flow cytometric analyses of Ly51 expression and UEA1 binding on TECs for mice; wt (n = 22), Cm_Foxn4 (n = 6), Cm_Foxn1 (n = 11), and Cm_dtg (n = 5). (B) Epithelial microenvironment of reconstructed thymi resolved by keratin 5 (K5) (in green) and K8 (in red) staining; 4-week-old mice. m, medulla; c, cortex. Note the small size of the Cm_Foxn4-reconstructed thymus. (C) Differential gene expression patterns in TEC transcriptomes relative to wild-type mice; Bl_Foxn4 (n = 4), Cm_Foxn4 (n = 4), and wt (n = 3). (D) Localization of Aire+ cells relative to cortical (K8) and medullary (K5) areas. (E) Differential gene expression patterns in TEC transcriptomes of transgenic relative to wild-type mice; wt (n = 3), Mm_Foxn4 (n = 6), and Cm_Foxn4 (n = 4). (F) Localization of B220+ B cells (in green) adjacent to ER-TR7+ mesenchyme (in red). PVS, perivascular space. (G) Differential gene expression patterns in TEC transcriptomes of transgenic relative to wild-type mice; wt (n = 3), Cm_Foxn4 (n = 4), Cm_Foxn1 (n = 3), and Cm_dtg (n = 5). In (C), (E), and (G), each data point represents the average value of at least three mice; all values are significantly different from the wild-type genotype (adjusted P values of <0.05), except those data points falling on “0” Arrows indicate the directions of changes in expression levels between transgenic mice.

Although no histologically distinct medullary region is detectable (Fig. 4B), the expression levels of genes indicative of different subsets of mTECs (Fig. 4C) are higher than in the Bl_Foxn4-driven epithelium. Although the morphology of the transgenic thymus lacks the characteristic zonation of cTECs and mTECs that is typical of the wild-type epithelium, low levels of Aire gene expression are detectable (Fig. 4C). Immunohistological analysis indicates the presence of only few Aire+ TECs, which are observed in the edges of epithelial cell–free regions in the reconstituted thymus, in contrast to the obvious Aire+ mTEC clusters in the wild-type thymus (Fig. 4D).

In line with the remarkable sequence conservation of vertebrate Foxn4 proteins (figs. S1, S2, and S5), the structures of the thymic microenvironments driven by Cm_Foxn4 (Fig. 4B) and Mm_Foxn4 (32) genes are essentially identical, despite 500 million years (Ma) of independent evolution. However, compared to the situation in Cm_Foxn4 transgenic mice, Mm_Foxn4-driven epithelia express lower levels of Prss16 and Aire, suggesting that with evolutionary progression, Foxn4 genes have gradually lost the capacity to support maturation of both cTECs and mTECs (Fig. 4E).

Previously, we have shown that the B cells generated in the Mm_Foxn4 thymus are preferentially located in close proximity to the vasculature, i.e., in the mesenchymal perivascular space (22). This finding highlighted the presence of anatomically distinct niches supporting the development of the two principal lymphocyte lineages in a primordial thymus; T cells differentiate in an epithelial environment, whereas B cells differentiate in a mesenchymal niche, similar to the situation in the sinusoidal environment of the bone marrow (33). In notable similarity, B cells in the Cm_Foxn4-expressing thymi are again found in the perivascular space (Fig. 4F), indicating that Foxn4 proteins favor the formation of anatomically separated domains specialized in either T or B cell differentiation.

Lymphopoietic properties of shark Foxn1

Next, we tested the thymopoietic capacity of a Foxn1 gene that is suggested by phylogenetic analysis to have emerged from the ancestral vertebrate Foxn1/4-like gene, following a gene duplication event at the base of vertebrate evolution (Fig. 1 and figs. S1 and S2). With respect to mammalian Foxn1 genes, the elephant shark Cm_Foxn1 gene is the evolutionarily most distant form of a Foxn1-like gene of jawed vertebrates examined here (Fig. 1 and figs. S1, S2, and S6). Despite more than 400 Ma of independent evolution, the lymphopoietic capacity of shark Cm_Foxn1 is remarkably similar to that of mouse Mm_Foxn1. The total numbers of CD4+CD8+ immature thymocytes in Cm_Foxn1-expressing thymi approach those observed in mouse wild-type thymi; moreover, differentiation into CD4+ and CD8+ single-positive T cells occurs; as a result, the relative proportions of immature and mature thymocytes are identical to wild-type thymi (Fig. 3, B and D). With respect to intrathymic B cell development, we find that the absolute numbers of immature B cells are about 10-fold lower in the Cm_Foxn1 transgenic thymi than in the corresponding Cm_Foxn4 reconstitution (Fig. 3, C and E); however, they are slightly increased when compared to the wild-type mouse thymus.

As expected from the near-normal hematopoietic compartment in Cm_Foxn1 thymi, their microenvironment also resembles a wild-type mouse thymus. This is evident from the cytometric profiles of Ly51 and UEA1 markers (Fig. 4A) and the anatomical segregation of K5+K8 (medullary) and K5K8+ (cortical) areas in tissue sections (Fig. 4B); except for a somewhat smaller overall size, the histology of the Cm_Foxn1 thymus is essentially identical to the wild-type form. As expected, Aire+ TECs are present in the distinct medullary areas, sparing the thymic cortex (Fig. 4E). The expression patterns of signature genes in Cm_Foxn1 TECs indicate the remarkable similarity to wild-type TECs in the mouse thymus (Fig. 4G). However, subtle differences in the gene expression profiles exist; the increased levels of Prss16 suggest the presence of a bias toward mature cTECs at the expense of certain mTEC populations, for instance, those expressing Trpm5 (Fig. 4G). Nonetheless, our observations suggest that the T cell–biased lymphopoietic properties of Foxn1-like genes, as exemplified by the mammalian thymus, had already emerged in the ancestor of jawed vertebrates and were maintained throughout subsequent evolution.

Reconstruction of a shark thymic microenvironment in the mouse

With respect to the reconstructions described above, it is important to note that they represent an artificial disentanglement of Cm_Foxn1 and Cm_Foxn4 functions, since in the thymus of cartilaginous fish, Foxn1 and Foxn4 paralogs are coexpressed (fig. S7) (22). However, the expression patterns of the two genes are not completely identical; whereas they are both expressed in the shark retina, the expression of Foxn4 in distinct cell clusters in the spinal cord is unique (fig. S7). To mimic the physiological coexpression of both genes in TECs of cartilaginous fish, we generated Cm_Foxn1/Cm_Foxn4 double-transgenic mice (hereafter Cm_dtg). When compared to the number of thymocytes in Cm_Foxn4 and Cm_Foxn1 single-transgenic mice, this constellation of coexpression results in an intermediate number of T cells, although the cellularity is closer to the situation in Cm_Foxn1 thymi (Fig. 3A). Unexpectedly, coexpression of Cm_Foxn1 and Cm_Foxn4 leads to markedly increased numbers of immature B cells, almost two orders of magnitude more than in the Cm_Foxn1 single transgenic (Fig. 3, C and E). The higher capacity for B cell development in the double-transgenic thymus does not affect their anatomical localization; B cells still reside in close proximity to the vasculature, as seen in thymi driven by the expression of Cm_Foxn4 alone (Fig. 4F). The TECs of the bipotent lymphoid organ in double-transgenic mice exhibit a phenotype intermediate between the single transgenics, as reflected in their cell surface phenotype (Fig. 4A); cTECs expressing high levels of Ly51 still predominate, although UEA1+ mTECs make up a quarter of the TEC compartment. The keratin expression pattern indicates that the medullary and cortical regions are not as precisely demarcated as in the Cm_Foxn1 single transgenic (Fig. 4B); however, the medulla contains a large number of Aire+ cells (Fig. 4E), in line with the gene expression pattern (Fig. 4G). In terms of lymphopoietic capacity, the double-transgenic TECs are distinguished by a markedly increased capacity for B cell development (about 10 times higher than that of the Cm_Foxn4 thymus), whereas the capacity for immature T cells is only slightly reduced compared to the Cm_Foxn1 single-transgenic thymus (Fig. 3, G and H). The unique morphological and functional characteristics of the double-transgenic thymus are reflected in an intermediate gene expression pattern, demonstrating that the combination of signature genes selected here faithfully report the thymopoietic characteristics of TECs. Collectively, our results indicate that the coexpression of Cm_Foxn1 and Cm_Foxn4 in the same TECs recapitulates the observed bipotent nature of lymphopoiesis in the thymi of extant cartilaginous fishes (34), providing an important validation of the biological relevance of the present reconstruction strategy.

Molecular basis of lymphopoietic specialization

Whereas Foxn4 proteins are well conserved during the course of vertebrate evolution (fig. S5), Foxn1 proteins are much more variable (Fig. 1 and fig. S6) in particular, with respect to the lengths and amino acid sequence compositions of their N-terminal domains. When viewed through the lens of the mouse proteins, the sequences encoded in coding exons 2 and 3 of mouse Foxn4 are replaced by unrelated sequences in a single exon (coding exon 2) of mouse Foxn1 (fig. S8). This observation suggested that, after gene duplication, the evolution of protein domain(s) in Foxn1 proteins was associated with increasing lymphoid selectivity, from bipotency of Foxn4 to unipotency of Foxn1. We set out to test this hypothesis by creating a chimeric protein, Foxn1*4, in which the sequence encoded by coding exon 2 of the mouse Foxn1 gene was replaced by coding exons 2 and 3 of mouse Foxn4 (Fig. 5A and fig. S8). The overall cellularity in the thymus of Mm_Foxn1*4 transgenic mice was reduced to about 20% of wild-type numbers, but the differentiation of T cells occurred normally, as revealed by the presence of normal percentages of CD4+ and CD8+ single-positive cells (Fig. 5, B and C). However, the CD4CD8 double-negative compartment was increased in the Mm_Foxn1*4 transgenic thymi, likely caused by an increase in the number of B cells (Fig. 5D). Expression of the Mm_Foxn1*4 chimeric protein led to a 10-fold increase of B cell poietic capacity, without appreciable loss of T cell poietic capacity, when compared to Mm_Foxn1 (Fig. 5, E and F).

Fig. 5 A chimeric Foxn1/Foxn4 protein induces B cell development in the thymic microenvironment.

(A) Schematic representation of the N-terminal domains of mouse Foxn4, Foxn1, and the Foxn1*4 chimera; boxes correspond to exons, and colored lines correspond to conserved amino acid residues in Foxn4 and Foxn1 clades (see figs. S2 and S8 for details). The > sign denotes the DNA binding and activation domains, which are not shown here. (B) Intermediate cellularity of Foxn1*4 thymi (***P < 0.001; two-tailed t test); wt (n = 5), Mm_Foxn4 (n = 6), and Foxn1*4 (n = 13). (C) Enlarged CD4CD8 DN compartment in Foxn1*4 thymi (P < 0.001; two-tailed t test); wt (n = 5) and Foxn1*4 (n = 13). (D) Moderately increased numbers of IgMCD93+ immature B cells (P = 0.293; two-tailed t test, compared to wt); wt (n = 3) and Foxn1*4 (n = 7). (E) Foxn1*4 supports intrathymic T cell development (P = 0.6654; two-tailed t test, compared to wt); wt (n = 5) and Foxn1*4 (n = 13). (F) Increased B cell development (*P = 0.0335; two-tailed t test, compared to wt); wt (n = 8), Mm_Foxn4 (n = 7), and Foxn1*4 (n = 13). In (E) and (F), data for Mm_Foxn4 transgenics (shaded area) are taken from (22). (G) Flow cytometric analyses of Ly51 expression and UEA1 binding of EpCAM+CD45 TECs; wt (n = 5) and Foxn1*4 (n = 13). (H) Epithelial microenvironment of reconstructed thymi resolved by keratin 5 (K5) (in green) and K8 (in red) staining; 4-week-old mice. (I) Localization of B220+ B cells (in green) adjacent to ER-TR7+ mesenchyme (in red); inset shows a higher magnification of the indicated region highlighting the perivascular space. (J and K) Differential gene expression patterns in TECs; see legend in Fig. 4 for details; wt (n = 3), Mm_Foxn4 (n = 6), Foxn1*4 (n = 3), and Cm_dtg (n = 5). In (B), (E), and (F), each data point represents one mouse.

Cell surface markers (Fig. 5G) of TECs in the Mm_Foxn1*4 thymus are reminiscent of the reconstructed thymus in Cm_dtg mice (Fig. 4A). This is evident from a 10-fold greater fraction of Ly51+UEA1 cTECs in the transgenic thymus, with a corresponding reduction of Ly51UEA1+ mTECs. This conclusion is supported by the pattern of K5 and K8 keratin expression and the histological features of cortical and medullary areas (Fig. 5H). Notably, as is the case in Cm_dtg transgenic thymi (Fig. 4F), B cells are situated in the perivascular space of the thymus of Mm_Foxn1*4-transgenic mice (Fig. 5I). Collectively, this finding demonstrates that the sequences in the evolutionarily dynamic N-terminal domains of Foxn1 and Foxn4 proteins have important roles in controlling the extent of B cell development in the thymic microenvironment. In further support of this conclusion, we find that the lancelet Foxn4 protein [which does not support B cell development (Fig. 2, C, D, and I)] lacks most of the conserved amino acid sequence signature of the vertebrate Foxn4 family in this domain (fig. S2). Last, as expected from flow cytometric profiles of TECs and the distinct histological features seen in the tissue sections, the gene expression signature of Foxn1*4-expressing TECs suggests a bias toward the cTEC lineage at the expense of mature mTECs (Fig. 5J), much like the situation in Cm_dtg mice (Fig. 5K).

Features associated with the transition from invertebrate to vertebrate Foxn1/4 functions

Our results provide a previously unattainable possibility to compare the transcriptional profiles of thymic microenvironments established by the activity of the cephalochordate Bl_Foxn4, the combinatorial activities of shark Cm_Foxn1 and Cm_Foxn4, and the activity of the mammalian Mm_Foxn1 genes. On the basis of the expression levels of selected TEC signature genes, a clear evolutionary trend becomes apparent: A gradual decrease in expression of cTEC-like genes is accompanied by an increase in expression of the genes that are characteristic of the mTEC compartment (Fig. 6A). This genetic constellation is associated with the more than 100-fold increase in overall T lymphopoietic capacity when comparing the thymic microenvironment established by the activities of the amphioxus Bl_Foxn4 and the mammalian Mm_Foxn1 genes. At the same time, appreciable B lymphopoietic capacity appears to have been a transient phenomenon, absent in Bl_Foxn4-driven and Mm_Foxn1-driven thymic microenvironments but present in the unique Foxn1 and Foxn4 coexpressing microenvironments characteristic of cartilaginous and teleost fishes (22). Hence, we became interested to determine the potential mechanistic basis of the bipotent microenvironment. Our previous experiments (22) suggested that the ratio of Dll4 and Il7 expression levels in the thymic microenvironment is an important determinant of its lymphopoietic properties, with respect to the balance between T and B cell poiesis. The Notch1/Dll4 signaling pathway is essential for the initiation of T cell development in the thymus (5, 15, 16, 35), whereas Il7 functions as a general lymphopoietic growth factor (36). The Dll4 gene is known to be a direct target of the Foxn1 (5, 37) and Foxn4 (4, 38) transcription factors, whereas the expression of Il7 is independent of Foxn1 activity (39). In support of the latter conclusion, we observed similar expression levels of the Il7 gene in all reconstituted TEC compartments examined here. Hence, since the denominator in the Dll4/Il7 ratio is a constant, it follows that the Foxn1/4-dependent expression levels of Dll4 determine the particular type of lymphopoietic activity in the thymus. In the reconstructions with Cm_Foxn1 and Cm_Foxn4 genes, we observed that, compared to the single transgenics, the Cm_dtg thymus had the highest capacity for B cell development (Fig. 6B), in line with the natural bipotency of the thymus of cartilaginous fish (34, 40). The ratio of T cell poiesis and B cell poiesis (here referred to as the T/B index) positively correlated with the ratio of expression levels of Dll4 and Il7 genes (Fig. 6C). In the comparison of Mm_Foxn4, Mm_Foxn1, and Mm_Foxn1*4 transgenics, the same positive relationship between T/B index and Dll4/Il7 ratio holds (Fig. 6, D and E). As compared to the situation in wild-type mice, Dll4 expression is several fold higher in Bl_Foxn4 transgenic microenvironments than it is in Cm_Foxn4-driven TECs, reflecting the T cell bias in the former and the presence of B cell development in the latter. Overall, the Dll4 expression levels vary by about one order of magnitude (Figs. 4, C, E, and G, and 5, J and K). Although it is likely that other factors affect the lymphopoietic properties of the thymic microenvironment, the modulation of Dll4 expression through Foxn1/4 transcription factors emerges as an evolutionarily conserved and functionally relevant mechanism by which the lymphopoietic capacity and the bias for or against B cell development of the thymus could be modulated.

Fig. 6 Transcriptional profiles of transgenic microenvironments highlight the formation of the thymic medulla during vertebrate evolution accompanied by a high unipotent capacity for intrathymic T cell development.

(A) Differential gene expression patterns in TECs; see legend in Fig. 2 for details. (B and D) T and B cell poietic indices; arrows indicate the altered balance between T and B cell generation. (C and E) Ratios of T and B cell indices as a function of the ratios of Dll4 and Il7 expression levels; Cm_Foxn4 (n = 6), Cm_Foxn1 (n = 7), Cm_dtg (n = 5), wt (n = 3), Mm_Foxn4 (n = 6), and Foxn1*4 (n = 7). ***P < 0.001 and *P < 0.05; one-way ANOVA with Tukey’s multiple comparison test; ±SDs correspond to propagated errors. (F) Schematic summarizing gene content and expression characteristics and associated lymphopoietic properties of thymi during vertebrate evolution.


The in vivo reconstitution experiments described here suggest a sequence of events during vertebrate evolution that culminated in the emergence of the T cell–biased thymus. We hypothesize that the primordial vertebrate Foxn4-like gene was expressed in the pharyngeal endoderm in the ancestor common to all vertebrates, much like its ortholog in lancelets (4). After the emergence of lymphocytes, which may have had their evolutionary origin in lymphocyte-like cells of tunicates (41), the Foxn4-expressing patch of pharyngeal endoderm cells may have supported the development of T-like cells; we further propose that this primordial type of lymphopoietic activity initially supported their development only up to the stage where their germline-encoded antigen receptors were expressed (Fig. 6F). Our previous studies demonstrated that the expression of the Notch ligand Dll4 and the chemokine Cxcl12 in Foxn1-deficient TECs suffices to support T cell differentiation up to the CD4+CD8+ DP stage (5), albeit at a much lower efficiency than we observe here in Bl_Foxn4 reconstitutions. Nonetheless, these results collectively suggest that a small number of effector molecules suffice to jump-start the formation of a lymphopoietic environment. We consider it likely that at this particular stage of immune system evolution, the facility of somatic recombination of antigen receptor genes and an associated quality control mechanism(s) mitigating any potential autoreactivity was not yet established (42). Nonetheless, lymphocytes in the vertebrate ancestor’s immune system may have expressed different kinds of germline-encoded antigen-specific receptors, analogous to pattern recognition receptors; variegated expression of these sensory modules would have afforded early vertebrates with the capability of immune responses and memory functions through clonal proliferation, akin to natural killer cells in mammals (43). Collectively, the phenotype of the Bl_Foxn4 epithelium defines a previously unidentified checkpoint during TEC differentiation, which marks the support of T cell development up to the CD4+CD8+ DP stage.

The emergence of a typical vertebrate-like Foxn4 gene (here exemplified by Cm_Foxn4) heralds another critical transition point in the immune systems of early vertebrates, as it established an environment supporting the development of the two principal lymphocyte lineages. Within the epithelial TEC compartment, it fostered the further development of T cells to reach the single-positive stage, albeit at low efficiency, and at the interface between epithelial and mesenchymal components of the microenvironment, it established conditions conducive to B cell development, as indicated by the presence of substantial numbers of immature B cells in and around the perivascular space. A comparison of protein sequences suggests that changes in the N-terminal domain of Bl_Foxn4 facilitated this process. Since the predicted Foxn4 protein of tunicates assumes an intermediate position in the phylogeny of chordate Foxn1/4 proteins (Fig. 1), it will be of interest to examine its capacity to support the development of B cells.

After the emergence of the Foxn1 gene in a primordial vertebrate, as a result of a gene duplication event, coexpression in the pharyngeal epithelium of the paralogous Foxn4 and Foxn1 genes was maintained; this coexpression signature still persists in extant cartilaginous (fig. S7) (4) and bony fishes (22, 44). Protein sequence comparisons indicate that the emergence of the Foxn1 paralog was accompanied by a radical modification of the amino acid sequence composition of the N terminus. This finding strongly supports the notion that exon replacement event(s) accompanied the emergence of the first Foxn1 genes that exchanged two exons of the Foxn4 gene by a single exon of substantially different sequence. In the Cm_Foxn1 protein, the N terminus is much shorter than it is in the living representatives of evolutionarily more recent taxa, such as the mouse Foxn1 protein, suggesting that structural features of this domain rather than primary amino acid sequence similarities underlie the equivalent functionalities of the Cm_Foxn1 and Mm_Foxn1 proteins.

Because the overall number of TECs present in the thymic lobes was essentially invariant in all reconstructed thymi and similar in magnitude to that in the wild-type mouse thymus, the main difference between Bl_Foxn4 and the vertebrate versions of Foxn1/4 proteins is the much lower lymphopoietic capacity of the former. This observation indicates that the key functions of the Foxn1/4 proteins are to qualitatively alter epithelial cell phenotypes rather than acting to simply quantitatively expand the epithelial compartment via proliferation. Our results indicate that the increase in lymphopoietic capacity occurred in a stepwise fashion, from Bl_Foxn4 to Cm_Foxn4 to CmFoxn1. In addition, the transition from Foxn4 to Foxn1 was accompanied by an increased size of the mTEC compartment, associated with a larger proportion of mature single-positive T cells among thymocytes. In the embryonic mouse thymus, cTECs develop earlier than mTECs (45, 46); this ontogenetic sequence closely resembles the phenotypes we observe in the phylogenetic sequence of Bl_Foxn4 to Cm_Foxn4 to Cm_Foxn1.

Another critical transition in the evolutionary trajectory of Foxn1/4 genes and vertebrate thymopoiesis occurred when cis-regulatory changes led to the loss of Foxn4 expression in TECs. This reorganization of genetic networks is exemplified by the reciprocal expression patterns of Foxn1 and Foxn4 in chicken tissues (fig. S9). As a result, T cell development became entirely dependent on Foxn1 (9, 47). At present, we do not know why a degenerate network structure [note that Foxn4 and Foxn1 are partially redundant in the teleost thymus (22)] was replaced by a nonredundant design. However, since the contribution of the thymus to B cell development was abolished in this process, loss of Foxn4 expression in TECs helped establish the strict anatomical segregation of developing lymphocyte lineages (48). As a result, the thymus was eventually transformed into an organ highly specialized for efficient T cell development (Fig. 6F).



C57BL/6J mice are maintained in the Max Planck Institute of Immunobiology and Epigenetics. Foxn1:Bl_Foxn4, Foxn1:Cm_Foxn4, and Foxn1:Cm_Foxn1 transgenic mice were constructed according to standard protocols by cloning a 27,970–base pair mouse Foxn1 promoter fragment (GenBank accession number Y12488; nucleotides 5680 to 33,650) upstream of either B. lanceolatum Foxn4 cDNA (GenBank accession number AJ252025.1; nucleotides 1 to 1590), C. milii Foxn4 cDNA (GenBank accession number FJ176202.1, nucleotides 56 to 1615), or a C. milii Foxn1 cDNA (GenBank accession number FJ176201.1; nucleotides 76 to 1584), followed by the bovine growth hormone polyadenylation sequence (32); cDNAs of Cm_Foxn4 and Cm_Foxn1 were synthesized by Eurofins Genomics and cloned into peX-K248 vector; in the Foxn1:Foxn1*4 construct, nucleotides 267 to 728 of M. musculus Foxn1 (GenBank accession number NM_008238.2) were replaced by nucleotides 170 to 445 of M. musculus Foxn4 (GenBank accession number AF323488.1). To generate transgenic mice, constructs were linearized and injected into FVB pronuclei according to standard protocols.

Transgenic mice were subsequently backcrossed to Foxn1-deficient mice (10) on a C57BL/6J background. Mice were kept in the animal facility of the Max Planck Institute of Immunobiology and Epigenetics under specific pathogen–free conditions. All animal experiments were performed in accordance with the relevant guidelines and regulations, approved by the review committee of the Max Planck Institute of Immunobiology and Epigenetics and the Regierungspräsidium Freiburg, Germany (license AZ 35-9185.81/G-14/57). Transgene expression levels were determined by RNA-seq, comparing transgene-derived BGH_PolyA transcripts to lacZ transcripts, the latter representing the activity of the targeted endogenous Foxn1 locus (10); the ratios of lacZ transcript counts to those emanating from the 5′ part of the Foxn1 gene [note that the 5′ end of the Mm-Foxn1 sequence is still detectable in the transcriptome of Foxn1−/− mice (10)] served as normalization. A negative feedback loop suppresses the endogenous Foxn1 locus activity for transgenes encoding proteins functionally equivalent to Foxn1; despite a 10-fold difference in expression levels, the transgenes are expressed in the range of the endogenous Foxn1 gene in wild-type mice (fig. S10). For identification of transgenes, the following primers were used for genotyping: Foxn1:Mm_Foxn1, SS40 (wild-type allele) + SS35 (wild-type and knockout allele) + JBS003 (knockout allele), Foxn1:Mm_Foxn1*4, XAH388 + SS6; Foxn1:Cm_Foxn1, XAH163 + RM19; Foxn1:Cm_Foxn4, XAH163 + RM22; Foxn1:Bl_Foxn4, JBS465 + JBS466; SS40, 5′-CTGTGAACTCAGCCATACTC; SS35, 5′-TGCACCAAGCCTCTGCTGGGA; JBS003, 5′-TCGCCTTCTTGACGAGTTCT; XAH388, 5′-CAGCAACTGATAAGGTCACC; SS6, 5′-ACAGAATTCTTCCAGCCATCA; XAH163, 5′-GTCCCTAATCCGATGGCTAGCTC; RM19, 5′-TATCGCGTGCACGAGTTGTA; RM22, 5′-GGTTAAAGTTCATGCGGCCG; JBS465, 5′-CCAGCTCCGAAACAGCCTAA; JBS466, 5′-GTCCTTTGTCGTCTGGTCGT.


Thymus organs were fixed for 120 min in 4% paraformaldehyde, washed with phosphate-buffered saline (PBS), and incubated in 20% sucrose overnight before mounting and snap-freezing in optimal cutting temperature (OCT) embedding compound. Tissue sections (8 μm) were cut using a cryostat and mounted onto precoated slides (Superfrost plus, Thermo Fisher Scientific). Slides were dried, followed by a 30-min blocking step using mouse immunoglobulin G (IgG) at 1:50 diluted in PBS + 0.5% bovine serum albumin (BSA) + 0.2% Tween. K5 K8 staining was performed with rabbit anti-K5 antibody (Ab) (PRB-160P, Covance) at 1:500 and rat anti-K8 Ab (Troma1, in-house) at 1:200. As secondary Ab, goat anti-rabbit Alexa Fluor 488 (A11008, Thermo Fisher Scientific) at 1:500 and donkey anti-rat Cy3 (AB_2340668, Jackson ImmunoResearch) at 1:500 were used. For ER-TR7 B220 staining, the rat anti–ER-TR7 Alexa Fluor 647 Ab (sc-73355 AF647, Santa Cruz Biotechnology) at 1:50 and rat anti-B220 Alexa Fluor 488 Ab (RA3-6B2, eBioscience) at 1:200 were used. Sections were mounted with Fluoromount G before analysis (Apotome, Zeiss). For combined K5/K8/Aire staining, sections were dried, blocked, and stained with unlabeled primary Abs as above and then stained with the secondary Abs donkey anti-rabbit IgG Cy3 (Jackson ImmunoResearch 711-165-152) and mouse anti-rat κ light chain Alexa Fluor 647 (MAR 18.5.28, purified and labeled in-house). Sections were then blocked with rat IgG and subsequently stained with rat anti-mouse Aire Alexa Fluor 488 Ab (5H12, eBioscience 14-5934-82) at 1:200.

RNA in situ hybridization

RNA in situ hybridization was carried out essentially as described (22) using the following probes: Gallus gallus Foxn1, nucleotides 1371 to 2557 in GenBank accession number XM_415816.6; G. gallus Foxn4, nucleotides 559 to 1779 in GenBank accession number NM_001083359.1; Scyliorhinus canicula Foxn1, nucleotides 37 to 346 in GenBank accession number FJ187748; and S. canicula Foxn4, nucleotides 1 to 422 combined from GenBank accession numbers Y11538, Y11539, and Y11540, respectively.

Flow cytometry

To generate single-cell suspensions for analytical and preparative flow cytometry of TECs, the procedures in (49) were followed. Note that the enzymatic cocktail required to liberate TECs destroys the extracellular domains of CD4 and CD8 surface markers (but not that of the CD45 molecule); hence, when analysis of thymocyte subsets was desired, thymocyte suspensions were prepared in parallel by mechanical liberation, best achieved by gently pressing thymic lobes through 40-μm sieves. Cell surface staining [anti-CD45 (30-F11), conjugated with phycoerythrin (PE) Cy7 (BioLegend); anti-EpCAM (G8.8), conjugated with allophycocyanin (APC; BioLegend); anti-Ly51 (BP-1) (6C3), conjugated with PE (eBioscience); UEA1, conjugated with fluorescein isothiocyanate (FITC; Vector Biosciences); anti-CD4 (GK1.5), conjugated with FITC (BioLegend); anti-CD8a (53-6.7), conjugated with APC (eBioscience); anti-CD19 (eBio1D3), conjugated with PerCPCy5.5 or PeCy7 (eBioscence); anti-B220 (CD45R) (RA3-6B2), conjugated with biotin (eBioscience); anti-IgM (II/4.1), conjugated with PE (eBioscience), anti-CD93 (C1qRp) (AA4.1), conjugated with APC (eBioscience); streptavidin conjugated with eFluor 450 or FITC (eBioscience)] was performed at 4°C in PBS supplemented with 0.5% BSA and 0.02% NaN3. Because of their small size in Bl_Foxn4 mice, the numbers of TECs and hematopoietic cells in thymi were determined independently; hence, B and T poietic indices were calculated from mean values with error propagation.

RNA-seq of TECs

Single-cell suspensions were prepared by TEC digest as described above. CD45EpCAM+ cells [negative enrichment using anti-CD45 magnetic-activated cell sorting (MACS) beads and anti–Ter-119 MACS beads, Miltenyi Biotec] were sorted directly into TRI reagent (T9424, Sigma-Aldrich). RNA isolation was performed according to standard protocols. Libraries were prepared using the Ultra RNA Library Prep Kit (Illumina). Samples were run on HiSeq2500 and sequenced to a depth of > 60 × 106 to 100 × 106 reads per sample.

Phylogenetic analysis

The relevant Foxn1 and Foxn4 amino acid sequences can be found under the following GenBank accession numbers: Foxn1: Rhincodon typus (XM_020525471); C. milii (XM_007896499); Amblyraja radiata (XM_033046380); Erpetoichthys calabaricus (XM_028808185); Acipenser ruthenus (XM_034049625); Lepisosteus oculatus (XM_015367325); Danio rerio (XM_009291615); Microcaecilia unicolor (XM_030186024.1); Nanorana parkeri (XM_018555008.1); Xenopus laevis (XM_018248776.1); Xenopus tropicalis (XM_018091796); Podarcis muralis (XM_028708697); Geotrypetes seraphini (XM_033921305); Anolis carolinensis (XM_016997340); G. gallus (XM_415816); Corvus cornix cornix (XM_019287554); Monodelphis domestica (XM_001375795); Ornithorhynchus anatinus (XM_029082501); M. musculus (NM_008238); Homo sapiens (NM_001369369); Foxn4: B. lanceolatum (AJ252025); Branchiostoma belcheri (XP_019621093); Phallusia mammillata (LR785254); R. typus (XM_020536376); A. radiata (XM_033043787); C. milii (NM_001292643); L. oculatus (XM_006640210); E. calabaricus (XM_028824927); D. rerio (NM_131099); Latimeria chalumnae (XM_014493627); Rhinatrema bivittatum (XM_029571604); X. laevis (BC142562); X. tropicalis (NM_001102862); G. seraphini (XM_033957269); P. muralis (XM_028705153); G. gallus (NM_001083359); C. cornix cornix (XM_020584559); Phascolarctos cinereus (XM_021008024); M. musculus (NM_148935); H. sapiens (NM_213596). The sequence of Foxn4 of Ciona intestinalis has been retrieved from the ENSEMBL database (ENSCING00000017653). Sequences were aligned using multiple sequence comparison by log-expectation (MUSCLE) (50), and the phylogenetic tree was reconstructed using the neighbor-joining method implemented in the BioNJ software (51) with 1000 bootstrap replicates and the Jones-Taylor-Thornton (JTT) substitution model (52). Both programs are available on the platform (53). The resulting phylogeny was visualized using the Interactive Tree of Life platform (54).

Data analyses

Transcriptomes were analyzed on the Galaxy platform using featureCounts (55) followed by DESeq2 analysis (56).

Statistical analysis

t tests (two-tailed) were used to determine the significance levels of the differences between the means of two independent samples, considering equal or unequal variances as determined by the F test. For multiple tests, the conservative Bonferroni correction was applied or as indicated, using one-way analysis of variance (ANOVA) with Tukey’s multiple comparison test.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank M. Petersen for help with the phylogenetic analyses and critical reading of the manuscript; B. Krauth, S. Birmelin, A. Haas-Assenbaum, and C. Happe for expert technical help; and the MPIIE deep sequencing and bioinformatics units for help with RNA-seq. We are grateful to S. Santini (CNRS/AMU IGS UMR7256) for maintaining the site Funding: This work was supported by the Max Planck Society and the Deutsche Forschungsgemeinschaft (Project 256073931-SFB 1160). R.M. was supported by a JSPS Overseas Research Fellowship. Author contributions: T.B. conceived the study and wrote the paper with input from all authors. J.B.S., A.N., R.M., and D.N. performed experiments. All authors designed the experiments, analyzed data, discussed the results, and reviewed the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The data reported in this paper have been deposited in the European Nucleotide Archive (ENA) database under the study accession number PRJEB38873.

Stay Connected to Science Advances

Navigate This Article