The Fox Gene Repertoire in the Annelid Owenia fusiformis Reveals Multiple Expansions of the foxQ2 Class in Spiralia

Abstract Fox genes are a large and conserved family of transcription factors involved in many key biological processes, including embryogenesis and body patterning. Although the role of Fox genes has been studied in an array of model systems, comprehensive comparative studies in Spiralia—a large clade of invertebrate animals including molluscs and annelids—are scarce but much needed to better understand the evolutionary history of this gene family. Here, we reconstruct and functionally characterize the Fox gene complement in the annelid Owenia fusiformis, a slow evolving species and member of the sister group to all remaining annelids. The genome of O. fusiformis contains at least a single ortholog for 20 of the 22 Fox gene classes that are ancestral to Bilateria, including an ortholog of the recently discovered foxT class. Temporal and spatial expression dynamics reveal a conserved role of Fox genes in gut formation, mesoderm patterning, and apical organ and cilia formation in Annelida and Spiralia. Moreover, we uncover an ancestral expansion of foxQ2 genes in Spiralia, represented by 11 paralogs in O. fusiformis. Notably, although all foxQ2 copies have apical expression in O. fusiformis, they show variable spatial domains and staggered temporal activation, which suggest cooperation and sub-functionalization among foxQ2 genes for the development of apical fates in this annelid. Altogether, our study informs the evolution and developmental roles of Fox genes in Annelida and Spiralia generally, providing the basis to explore how regulatory changes in Fox gene expression might have contributed to developmental and morphological diversification in Spiralia.


Introduction
Forkhead box-containing proteins (i.e. Fox genes) form one of the largest families of transcription factors in animals, displaying a remarkable functional diversity in many morphogenetic processes (Kaufmann and Knöchel 1996;Carlsson and Mahlapuu 2002;Hannenhalli and Kaestner 2009). Fox genes are characterized by a conserved DNA-binding domain of approximately 100 amino acids -the Forkhead or winged helix domain-that folds into two stereotypical large loops or "wings" (Clark et al. 1993;Li and Tucker 1993). Since the discovery of the Forkhead domain in the fruit fly Drosophila melanogaster (Weigel et al. 1989), Fox genes have been studied in a wide range of traditional developmental systems, mostly vertebrates (Mazet et al. 2003;Lee and Frasch 2004;Hannenhalli and Kaestner 2009;Jackson et al. 2010;Golson and Kaestner 2016). The initial description of 15 Fox gene classes in chordates, each identified by a letter (Kaestner et al. 2000), and the establishment of a unified nomenclature facilitated phylogenetic analyses and comparisons with other major invertebrate clades, such as hemichordates (Fritzenwanker et al. 2014), molluscs (Yang et al., 2014;Wu et al., 2020), platyhelminthes (Pascual-Carreras et al. 2021), panarthropods , cnidarians (Magie et al. 2005), and sponges (Larroux et al. 2006), as well as animal outgroups (Larroux et al. 2008). This uncovered a complex evolutionary history for this large family of transcription factors. Today, Fox genes are classified into as many as 27 classes belonging to two major clades (Kaestner et al. 2000;Mazet et al. 2003;Larroux et al. 2008;Hannenhalli and Kaestner 2009;Benayoun et al. 2011;Pascual-Carreras et al. 2021;Schomburg et al. 2022), where ancestral duplication events (e.g. the former class foxQ split into foxQ1 and foxQ2, foxN into foxN1/4 and foxN2/3, foxL into foxL1 and foxL2/3, and foxJ into foxJ1 and foxJ2/3), gene innovations (e.g. foxR and foxS are unique of vertebrates, foxT is potentially a novelty of panarthropods), expansions and losses (e.g. foxAB in vertebrates, foxQ2 in tetrapods, and foxAB, foxE, and foxI in panarthropods) are common (Mazet et al. 2003;Tu et al. 2006;Wotton and Shimeld 2006;Paps et al. 2012;Schomburg et al. 2022). Moreover, genomic comparisons also uncovered signs of conserved syntenic linkage for some of the classes, such as foxL1-foxC-foxF-foxQ1, in phylogenetically distant lineages of insects, chordates and spiralians (Mazet et al. 2006;Shimeld 2006, 2011;Wotton et al. 2008;Shimeld et al. 2010a;Schomburg et al. 2022). Yet a comprehensive characterization of Fox genes is still lacking in most major animal groups, most notably in members of Spiralia, one of the two main clades of protostomian animals that comprises nearly half of the extant major metazoan groups, including molluscs and annelids (Marlétaz et al. 2019). Therefore, the evolutionary history and developmental roles of this conserved family of transcription factors are still unclear at key nodes of the animal tree of life.
Fox genes typically show tissue-specific expression patterns and play an important role in cell-type determination and differentiation (Hannenhalli and Kaestner 2009;Jackson et al. 2010). Functional studies in human, mouse, zebrafish, and fly have revealed an array of functions of Fox genes in early development, such as axial patterning, germ layer specification, and organogenesis (reviewed in Carlsson and Mahlapuu, 2002). In Spiralia, however, studies on the function of Fox genes are scarce and mostly focused on certain classes, with just a handful of studies encompassing more than one major spiralian clade (Supplementary table 1 and references therein, Supplementary Material online). For example, foxA is consistently expressed in the developing foregut in many spiralians, including annelids, brachiopods, phoronids, and bryozoans (Arenas-Mena 2006; Seaver 2008, 2010;Adler et al. 2014;Martín-Durán et al. 2016;Vellutini et al. 2017;Kwak et al. 2018;Andrikou et al. 2019;Kostyuchenko et al. 2019) and foxJ1, foxQ2 and foxG are expressed in larval specific tissues in the annelid Platynereis dumerilii, the brachiopod Terebratalia transversa and the phoronid Phoronopsis harmerii (Santagata et al. 2012;Marlow et al. 2014;Gą siorowski and Hejnol 2020). Similarly, the clustered classes foxC, foxL1 and foxF show mesodermal expression in all spiralian species studied to date, suggesting that coordinated activation of these Fox genes in a common germ layer might have contributed to the maintenance of their genetic linkage (Shimeld et al. 2010a;Passamaneck et al. 2015;Martín-Durán et al. 2016). Other Fox gene classes, however, have only been studied in individual species, which prevents inferring an ancestral role for these genes in Spiralia. For example, foxL2/3 is a regulator of ovarian differentiation and development in molluscs (Liu et al. 2012;Yang et al. 2014;Teaniniuraitemoana et al. 2015;Li et al. 2016), foxB is expressed during late mesoderm development in the leech Helobdella austinensis (Kwak et al. 2018), foxO controls tissue regeneration and cell death in the planarian Schmidtea mediterranea (Pascual-Carreras et al. 2021) and foxK is involved in ectodermal regeneration in that same planarian species (Coronel-córdoba et al. 2022). Consequently, the repertoire and developmental functions of most Fox genes remain largely unexplored in Spiralia, and thus its study is not only important to discern the evolution of this gene family in animals, but also the contribution of these developmental regulators to the diversification of body plans and embryonic modes in this major animal group.
Here, we mined the genome of the annelid Owenia fusiformis (Martín-Zamora et al. 2022), a member of Oweniidae and sister group to all remaining annelids (Rouse et al. 2022), and eight other annelids with highquality genomic and transcriptomic datasets to infer the ancestral Fox gene complement to Annelida, one of the most species-rich and morphologically diverse groups within Spiralia. Temporal and spatial gene expression analyses offer insights into the potential role of some of the Fox gene classes in O. fusiformis, uncovering conserved and putative new roles for some of the Fox classes. Moreover, our study reveals that the foxQ2 class is largely expanded in Spiralia, with the paralogs being consistently expressed in apical territories and exhibiting signs of possible sub-functionalization in O. fusiformis. Altogether, our study informs the evolution of the Fox gene family in Annelida, providing valuable data to reconstruct the evolution and developmental roles of these genes in Spiralia and Metazoa.

Results
The Fox Gene Complement in O. fusiformis To characterize the Fox gene complement in the annelid O. fusiformis, we searched for annotated gene models containing the conserved Forkhead DNA-binding domain in its reference genome (Martín-Zamora et al. 2022) and in eight other species with high quality, publicly available genomic and transcriptomic datasets (Simakov et al. 2013;Chou et al. 2018;Li et al. 2019;Shao et al. 2020;Martín-Durán et al. 2021;Sun et al. 2021;Zakas et al. 2022). We obtained a total of 35 putative Fox genes in O. fusiformis, and between 28 and 48 in other annelids, which is a higher number than those previously reported in other Spiralian species, in which the number of Fox genes ranges from 21 to 26 genes (Yang et al. 2014;Wu et al. 2020;Pascual-Carreras et al. 2021). To assign the orthology of each of the O. fusiformis Fox genes, we applied maximum likelihood and Bayesian phylogenetic tree inference, obtaining strongly supported orthologs for 21 of the 24 classes of Fox genes, both in clade I and clade II (figs. 1 and 2; Supplementary figs. 1-3, Supplementary Material online). Only the two ancestral bilaterian Fox gene classes foxE an foxI are missing in O. fusiformis, which have been however reported in other spiralian lineages (Yang et al. 2014;Wu et al. 2020). Nonetheless, these two Fox gene classes show a dynamic pattern of loss/expansion in Annelida, with foxE being lost in earthworms and leeches and expanded in the spionid Streblospio benedicti, and foxI being generally lost in Annelida and only retained in the earthworm Eisenia andrei and S. benedicti ( fig. 2). Indeed, earthworms, leeches, and S. benedicti have divergent Fox gene complements, with many losses and expansions, in contrast to most other annelid species that have relatively stable and more complete Material online), which indicates that this recently described Fox gene class in Panarthropoda is more ancient that initially thought ). In addition, O. fusiformis has 14 fast-evolving orthologs that are either related to foxQ2 in clade I or foxN2/3 in clade II ( fig. 1; Supplementary figs. 2 and 3, Supplementary Material online). While the 11 foxQ2-like genes could be confidently assigned to the foxQ2 class (see below), the orthology of the divergent Fox genes belonging to clade II could not be resolved and thus we deem them "orphan" genes (figs. 1 and 2). Similarly, most other annelids show a varying number of orphan Fox genes ( fig. 2), most of them related to foxQ2 and foxN2/3 Fox genes belonging to clade I and clade II additionally differ in whether they lack or contain introns at conserved sites of the Forkhead domain, respectively (Larroux et al. 2008). To further characterize the Fox gene complement and assess this rule in O. fusiformis, we reconstructed the domain architecture and exon-intron composition of all 35 Fox genes in this annelid. Fox genes in O. fusiformis show diverse gene architectures, with lengths ranging from 675 bp ( foxQ2-8) to 3,372 bp ( fox orphan-3), and all but foxQ2-10 contain an intact Forkhead domain ( fig. 3A). One Fox gene-fox orphan-3-shows two additional protein domains, namely a RING finger and RAWUL domains ( fig. 3A), which we confirmed by RNA-seq sequencing, and that might indicate that this divergent Fox gene have acquired new protein functions. Exon numbers in the Fox genes of O. fusiformis range from one to 11, and most genes follow the clade I and clade II distinction based on the number of introns, apart from foxAB-2, foxF, foxL2/3, and foxQ2-6, which contain introns albeit they belong to clade I and might thus represent independent intron gains ( fig. 3A; Supplementary table 2, Supplementary Material online). Altogether, these findings reinforce the previous observation that O. fusiformis contains a relatively well conserved Fox gene repertoire, while they highlight as well that a handful of Fox gene classes and paralogs might have experience faster rates of molecular evolution.
Certain Fox gene classes (e.g. foxL1, foxC, foxF, and foxQ1) show conserved chromosomal linkage across phylogenetically distant bilaterian taxa, such as panarthropods, vertebrates and amphioxus, and there are evidences of this linkage in some spiralian species (Mazet et al. 2006;Shimeld et al. 2010a;Shimeld et al. 2010b;Schomburg et al. 2022). To assess whether this feature was also retained in O. fusiformis, we studied the chromosomal location and microsyntenic relationship of the 35 Fox genes in this annelid. Fox genes are spread across nine of the 12 chromosomes of O. fusiformis, namely chromosomes 1, 2, 3, 5, 6, 7, 8, 11, and 12 ( fig. 3B). While foxF, foxC, foxL1, and foxQ1 are located on the same chromosome-chromosome 11-only foxC and foxL1 show evidence of a tight linkage, with a high number of genes (>1,000) lying between foxQ1 and foxF as well as between foxF and foxC ( fig. 3B). In addition, we observed tandem duplications of foxQ2 genes in

The Expression Dynamics of the Fox Gene Complement in O. Fusiformis
To investigate the expression dynamics of the Fox genes in O. fusiformis and relate each of these genes to major morphogenetic events during the life cycle of this annelid, we used available stage-specific RNA-seq data covering 14 developmental time points, from the unfertilized oocyte to the juvenile stage (Martín-Zamora et al. 2022). In O. fusiformis, the temporal expression dynamics of the Fox genes seem to correlate with their assignment to clade I and clade II, because half (6) of the clade II Fox genes are expressed maternally and during the early cleavage stages, while clade I Fox genes tend to show short peaks of expression at single developmental stages, from the 32-cell stage Genome Biol. Evol. 14(10) https://doi.org/10.1093/gbe/evac139 Advance Access publication 13 September 2022 onwards ( fig. 4A). Clade II genes that escape this trend are foxP, expressed at the juvenile stage (but also briefly at the blastula stage), foxN2/3 and the orphan-3 gene, which peak during gastrulation, and foxJ1, foxK, and foxO, that are expressed during larval development ( fig. 4A). Notably, a set of Fox genes belonging to clade I (foxAB-1, foxG, and many foxQ2 paralogs, see below) are finely expressed at the time of the specification of the embryonic organizer and the establishment of the axial identities in O. fusiformis, as well as during gastrulation (i.e. foxAB-2, foxH, foxA) (

O. fusiformis has an Expanded foxQ2 Class
Previous studies indicated that an ancestral duplication in the foxQ2 class occurred at least in the last common bilaterian ancestor, and maybe even predated the cnidarianbilaterian split (Chevalier et al. 2006;Fritzenwanker et al. 2014;Pascual-Carreras et al. 2021), which was followed by further duplications of this Fox gene class in some spiralian and deuterostomian lineages (Yang et al. 2014;Wu et al. 2020). This observation agrees with the large expansion of 11 foxQ2 paralogs that we identified in the annelid O. fusiformis ( fig. 1), which is mirrored in many other annelid lineages. Moreover, the presence of an Engrailed Homology 1 (EH)-i-like Groucho binding motif in some foxQ2 orthologs and its variable C-and N-terminal position with respect to the Forkhead and FoxQ2 domain has been used to subdivide the foxQ2 class in foxQ2-C and foxQ2-N, or even to define a new sub-class named foxQD (Fritzenwanker et al. 2014;Pascual-Carreras et al. 2021). To clarify the evolution of this Fox gene class, and how the expansions of foxQ2 genes occurred in O. fusiformis and spiralians generally, we mined available databases and the genomes of seven annelid species in search for genes with complete FoxQ2 domains, which we then used for phylogenetic reconstruction and the identification of EH-i-like motifs. This is a stringent approach that supported the ascription of most of the divergent clade I Fox genes to the foxQ2 class (e.g. the 11 foxQ2-like genes of O. fusiformis; Supplementary fig. 2, Supplementary Material online), with those not meeting the criteria being considered as orphan Fox genes ( fig. 2). While the general orthology of all identified foxQ2 genes was robustly supported ( fig. 5A; Supplementary figs. 6-8, Supplementary Material online), we did not recover two separate monophyletic clades with EH-i-like motifs in either the C-or the N-terminus. Instead, foxQ2 orthologs with an EH-i-like motif at the C-terminal end ( fig. 5B) appear to form a relatively well supported monophyletic clade, comprising sequences of both cnidarians and most bilaterian groups, often as single copy genes (yet the annelids O. fusiformis, Dimorphilus gyrociliatus, and Helobdella robusta have two paralogs). The rest of the foxQ2 sequences lack an EH-i-like motif and probably represent more or less divergent copies. Among these fast-evolving foxQ2 copies we found lineage-restricted expansions, such as those of O. fusiformis-for which the phylogenetic relationship between paralogs correlate well with their genomic linkage-and the vestimentiferan annelids Lamellibrachia luymesi and Paraescarpia echinospica, which have a group of foxQ2 genes that independently aquired an EH-i-like motif on the N-terminal end ( fig. 5A). Together, our findings corroborate previous analyses revealing a complex evolutionary history for the foxQ2 gene class, probably dominated by fast rates of molecular evolution and/or frequent independent events of gene duplication in both cnidarian and bilaterian lineages, specially among spiralians, as well as its complete loss in tetrapods.

The Developmental Expression of foxQ2 Paralogs in O. fusiformis
To assess the potential functional implications of the expansions of foxQ2 genes in spiralians ( fig. 6A), we first compared the developmental expression profiles of foxQ2 genes in O. fusiformis and three other spiralian species, namely the molluscs Crassostrea gigas and Mizuhopecten yessoenssis and the annelid C. teleta ( fig. 6B-E). In the pacific oyster C. gigas, the three foxQ2 paralogs (Supplementary table 5  foxQ2-1 having a maximum at the blastula stage, foxQ2-2 being strongly expressed during gastrulation, and both foxQ2-3 and foxQ2-4 peaking at larval stages ( fig. 6C). The six foxQ2 paralogs of C. teleta display four different temporal patterns of expression: foxQ2-1 is expressed in the oocyte, foxQ2-2 and foxQ2-3 peak during early cleavage, foxQ2-4 is expressed at the blastula stage, and foxQ2-5 and foxQ2-6 show higher expression values during gastrulation ( fig. 6D). Consistently, foxQ2 paralogs are also expressed at distinct developmental stages in O. fusiformis, with foxQ2-11 expressed first maternally and during early cleavage, foxQ2-1, foxQ2-4, foxQ2-5, foxQ2-6, foxQ2-9, foxQ2-10 showing a peak of expression at the early blastula stage, foxQ2-2 and foxQ2-3 exclusively expressed at the time of the specification of the embryonic organizer at 5 hours post-fertilization and their expression being under control of the ERK1/2 signaling pathway (Seudre et al. 2022), and finally both foxQ2-7 and foxQ2-8 coming up at larval stages ( fig. 6E). Notably, spiralian foxQ2 orthologs peaking at the blastula stage (i.e. foxQ2-2 and foxQ2-3 in O. fusiformis, foxQ2-4 in C. teleta, foxQ2-1 in M. yessoensis, and foxQ2-2 in C. gigas) belong to the strongly supported clade of foxQ2 genes with a conserved C-terminal EH-i-like motif. Similarly, foxQ2-1 in C. teleta and fusiformis, at the blastula (lateral view), gastrula (lateral to the left and ventral to the right) and larval stages (lateral view). Consistent with the temporal expression data, foxJ1 starts to be expressed at the putative prototroch precursor cells at the gastrula stage and it is later detected in the ciliated cells of the larva (apical organ and prototroch). Asterisks mark the apical/anterior pole. ac, archenteron; ao, apical organ; bp, blastopore; gp, gastral plate; pt, prototroch. (C) Table summarizing the current knowledge of the developmental roles of Fox genes in Spiralia (from Supplementary table 1). For each gene, reported developmental roles are in the light coloured box to the left and the associated species are shown with a dot on the corresponding horizontal line on the right.
Genome Biol. Evol. 14(10) https://doi.org/10.1093/gbe/evac139 Advance Access publication 13 September 2022 GBE foxQ2-11 in O. fusiformis are phylogenetically related and both show maternal expression. Therefore, these findings support that some of the expansions of foxQ2 genes are shared across spiralian lineages and that the expansion of this class of Fox genes may also resulted in the evolution of novel expression dynamics.
We next analyzed whether foxQ2 paralogs also showed varying spatial expression domains in O. fusiformis using whole mount in situ hybridization across targeted developmental stages. Except for foxQ2-8 and foxQ2-7, for which we did not observe any clear expression pattern at the larval stage, all other nine foxQ2 paralogs showed distinct  Genome Biol. Evol. 14(10) https://doi.org/10.1093/gbe/evac139 Advance Access publication 13 September 2022 GBE domains of expression. The paralog foxQ2-11, which has a high expression in the oocyte, is detected in the zygote and appears equally distributed in all blastomeres at the two-cell and four-cell stages ( fig. 6F). The genes foxQ2-1, foxQ2-4, foxQ2-5, foxQ2-6, foxQ2-9, foxQ2-10, which peak at 4 hours post-fertilization ( fig. 6G) and are probably independent expansions in O. fusiformis ( fig. 5A), are expressed at the apical pole, but encompassing varying areas of expression, from ones restricted to just the animal micromeres (e.g. foxQ2-6 and foxQ2-9) to nearly the entire animal hemisphere (e.g. foxQ2-4, foxQ2-5 and foxQ2-6) ( fig. 6H). Finally, the expression of the two foxQ2 paralogs with a C-terminal EH-i-like motif-foxQ2-2 and foxQ2-3-is restricted to the apical-most micromeres from 5 hours postfertilization onwards, consistent with the observed temporal dynamics and its regulation by ERK1/2 activity at that timepoint ( fig. 6G and I). Altogether, our gene expression analyses support that the multiple copies of foxQ2 have retained the evolutionarily conserved expression of this Fox gene class in anterior/apical development (Yaguchi et al. 2008;Santagata et al. 2012;Hunnekuhl and Akam 2014;Marlow et al. 2014;Fritzenwanker et al. 2014;Martín-Durán et al. 2015;Janssen 2017;He et al. 2019;Schacht et al. 2020), yet they have probably undergone temporal, spatial, and regulative sub-functionalization.

Discussion
The Evolution of the Fox Gene Complement in Annelida and Spiralia Taking advantage of the reference genome assembly and comprehensive functional data available for the annelid O. fusiformis (Martín-Zamora et al. 2022), our study identified and characterized the full Fox gene repertoire in a member of Oweniidae, the sister lineage to all remaining annelids (Rouse et al. 2022). The 35 Fox genes of O. fusiformis, belonging to 20 of the 22 classes predicted for the bilaterian ancestor (Kaestner et al. 2000;Benayoun et al. 2011), are amongst the most complete Fox gene repertoires for a spiralian reported to date, and thus helps to clarify the ancestral Fox gene complement in Spiralia and Lophotrochozoa (Marlétaz et al. 2019), as well as the evolution of certain Fox gene classes in this animal group. Indeed, our study demonstrates that O. fusiformis and the brachiopod Lingula unguis have bona fide foxT genes, revealing that this recently described Fox gene class is more ancient that initially thought and not limited or specific to Panarthropoda . Therefore, all 21 Fox gene classes present in O. fusiformis, together with the missing foxE and foxI classes, were present in the last common spiralian and annelid ancestor, thus setting a reference complement for subsequent studies on the diversification of Fox genes in these animal clades.
Our study supports that lineage-specific losses and expansions are common in the evolution of the Fox genes in Spiralia. While a small number of Fox gene classes have been independently lost in some annelid lineages (e.g. foxE and foxI in O. fusiformis, foxI and foxT in C. teleta), our study supports that a full repertoire of Fox genes was present in the last common annelid ancestor ( fig. 2). This contrasts with Mollusca (or at least Conchifera), which has experienced ancestral losses of two classes, namely foxI and foxQ1, as well as independent lineage-restricted losses (e.g. foxG and foxM in C. gigas, foxO in Lottia gigantea) (Shimeld et al. 2010a;Yang et al. 2014). Consistent with its faster rate of molecular evolution, the platyhelminth S. mediterranea has a more degraded Fox gene repertoire, with up to nine Fox classes missing (foxAB, foxB, foxE, foxH, foxL2/3, foxQ1, foxJ2/3, foxM and foxN1/4) (Pascual-Carreras et al. 2021). Notably, and as discussed below, expansions of the foxQ2 class are common in Annelida and Mollusca, with some lineages, including O. fusiformis, exhibiting large lineage-specific expansions. Future work on the expression and function of Fox genes across spiralian lineages will help to clarify how these expansions and losses impact developmental programs and the diversification of body plans in Spiralia.
Gene architecture differences between clade I and II and the ancestral chromosomal linkage of the foxC, foxF, foxL1 and foxQ1 classes characterize the evolution of Fox genes (Larroux et al. 2008;Shimeld et al. 2010a;Schomburg et al. 2022). While our gene architecture data generally support the previous observations that the clade I of Fox genes are intronless, O. fusiformis-as also observed in the annelid C. teleta-does not exhibit a Fox gene cluster involving foxC, foxF, foxL1 and foxQ1 (Shimeld et al. 2010a). This contrasts with the overall conservation of gene macrosynteny in O. fusiformis (Martín-Zamora et al. 2022) and suggests that some intra-chromosomal rearrangements might have happened in this species. Despite the lack of a cluster organization, however, these genes retain a mesodermal expression in O. fusiformis (at least for foxC, foxF, foxL1; see below) (Martín-Durán et al. 2016), as observed in C. teleta and other spiralians (Shimeld et al. 2010a). Therefore, our study challenges the hypotheses that the concerted mesodermal co-expression of foxC, foxF, foxL1 and foxQ1 genes was the selective pressure for maintaining their cluster integrity (Shimeld et al. 2010a;Shimeld et al. 2010b), suggesting instead that more local, gene-specific regulation is responsible for their expression dynamics.

The Expression Dynamics of the Fox Gene Complement in O. fusiformis
By combining new and previously reported temporal and spatial expression data in O. fusiformis, our study provides a comprehensive view of the developmental dynamics of Fox genes in this annelid species, suggesting conserved and potentially new developmental roles for certain Fox classes. In most animals studied to date, including annelids, molluscs, brachiopods, phoronid, bryozoans, planarians, and nemertean species, foxA is a key effector of foregut formation and a marker of endodermal tissues (Lartillot et al. 2002;Arenas-Mena 2006;Seaver 2008, 2010;Jr Kai Yu et al. 2008;Fuchs et al. 2011;Adler et al. 2014;Fritzenwanker et al. 2014;Martín-Durán et al. 2015, 2016Perry et al. 2015;Vellutini et al. 2017;Kwak et al. 2018;Andrikou et al. 2019;Kostyuchenko et al. 2019;Pascual-Carreras et al. 2021). In O. fusiformis, foxA is first detected in the vegetal macromeres at the blastula stage and later observed in the gastrula endoderm (peak of expression) and in the mouth and midgut of the developing larvae (Martín-Durán et al. 2016), thus supporting the role of the foxA class in endoderm and gut formation in animals. Similarly, the gene foxJ1 is involved in ciliogenesis in a number of animals (Xianwen Yu et al. 2008;Choksi et al. 2014), including during the formation of the prototroch in the annelid P. dumerilii (Marlow et al. 2014). In O. fusiformis, foxJ1 is expressed soon after gastrulation in the presumptive prototroch precursors and later on in the heavily ciliated cell types of the larva (the apical organ and prototroch) (Carrillo-Baltodano et al. 2021), reinforcing a conserved role of this Fox gene in the development of ciliated organs in Annelida and metazoans in general.
Existing gene expression data also support a conserved role of several Fox gene classes in mesoderm development in O. fusiformis. The foxH class regulates mesoderm development and embryonic organizing activity in vertebrates (Pogoda et al. 2000;Hoodless et al. 2001;Kofron et al. 2004), and it is downstream of the embryonic organizer and likely involved in mesoderm and posterodorsal development in O. fusiformis (Seudre et al. 2022), with a similar temporal expression dynamic reported in the oyster C. gigas (Yang et al. 2014). Similarly, the temporal and spatial expression domains of foxL1, foxC and foxF support their role in mesoderm formation in O. fusiformis (Martín-Durán et al. 2016), as described in other spiralians and metazoans (Wotton et al. 2008;Shimeld et al. 2010a;Wotton and Shimeld 2011;Fritzenwanker et al. 2014), albeit they do not retain their linked chromosomal position in this annelid species. Notably, foxQ1 which is absent from all molluscan species studied to date, is expressed in the stomodeum and pharynx in C. teleta (Shimeld et al. 2010a;Yang et al. 2014;Wu et al. 2020). However, foxQ1 is only expressed at the juvenile stage in O. fusiformis, suggesting that the role of this Fox gene class might differ between annelid and spiralian species.
The expression of Fox genes in other annelids and spiralians combined with the temporal dynamics of orthologs in O. fusiformis provide evidence of the potential roles of certain Fox gene classes in this annelid species. For instance, foxD is involved in myogenesis and ventral patterning in the annelid P. dumerilii, the platyhelminth S. mediterranea and the brachiopod T. transversa (Lauri et al. 2014;Passamaneck et al. 2015;Pascual-Carreras et al. 2021) and it becomes expressed after gastrulation and specially during organogenesis and initiation of myogenesis in O. fusiformis (Carrillo-Baltodano et al. 2021). The foxO class is a regulator of cell death in the planarian S. mediterranea (Pascual-Carreras et al. 2021) and an effector of cell division during early cleavage in the annelid H. austinensis (Kwak et al. 2018). In O. fusiformis, foxO is also expressed maternally and during early cleavage, as well as during embryonic periods with active cell turnover (Carrillo-Baltodano et al. 2021). The expression of the foxAB class has been only studied in the annelid C. teleta, which has a single ortholog that becomes expressed in a unique D-quadrant cell during early cleavage and is later involved in ectoderm differentiation and foregut formation (Boyle et al. 2014). Owenia fusiformis has instead two foxAB paralogs, one (foxAB-1) exclusively expressed at the time of the specification of the organizer cell at 5 hours post-fertilization and a second one ( foxAB-2) that peaks at the gastrula stage before gradually fading away during larval development. Although further expression analyses are needed, we speculate that the two paralogs in O. fusiformis could play similar roles than those described in C. teleta, with foxAB-1 acting during the axial body patterning and foxAB-2 being involved in ectoderm and foregut formation later in embryogenesis.
Our comprehensive developmental time course of the Fox genes in O. fusiformis also uncovered dynamics of expression for Fox gene classes for which there is little understanding of their roles during annelid and spiralian embryogenesis. For instance, foxT is mostly expressed in the oocyte, early development and pre-gastrulation, which contrast with the generally late developmental expression reported in arthropods (Lin et al. 2021;Janssen et al. 2022). The genes foxJ2/3, foxM, and foxN1/4 also show maternal expression, and for foxM and foxN1/4 their expression lasts until the blastula stage (Supplementary fig.  5, Supplementary Material online), where they are either expressed ubiquitously ( foxM) or at the apical pole (foxN1/4), suggesting that they might be potential regulators of early cleavage and/or cell fate specification in this species. The genes foxG, foxN2/3 and foxP are expressed at the time of the specification of the embryonic organizer in O. fusiformis (Seudre et al. 2022) (fig. 4A), and could therefore be involved in the establishment of the embryonic polarity and body plan in this annelid, as suggested by the expression of foxN2/3 and foxP in the endoderm and apical/ anterior ectoderm, respectively (Supplementary fig. 5, Supplementary Material online). Finally, some Fox genes are restricted to either the larva (foxK) or the juvenile (foxB, foxL2/3 and foxQ1), with expression data suggesting that at least foxK might be involved in the development of antero-ventral and oral ectoderm (Supplementary fig. 5, Supplementary Material online). Altogether, our study sets the stage for further expression and functional studies of Fox genes in O. fusiformis and spiralian embryogenesis, which ultimately will help to better understand the plasticity and development roles of this major family of transcription factors in animal development and evolution.

The Complex Evolutionary History of foxQ2 Genes
The foxQ2 class often comprises at least two paralogs in many cnidarian and bilaterian lineages studied to date, except in Tetrapoda, which lost this Fox gene class (Mazet et al. 2003;Yu et al. 2003;Chevalier et al. 2006;Santagata et al. 2012;Sinigaglia et al. 2013;Fritzenwanker et al. 2014;Marlow et al. 2014). Notably, foxQ2 genes show a remarkable conservation of their expression patterns across phylogenetically distant animal lineages. In deuterostomes like the cephalochordate Branchiostoma floridae, the echinoderm Strongylocentrotus purpuratus and the hemichordate Saccoglossus kowalevskii, foxQ2 genes are expressed in the apical pole during embryogenesis (Fritzenwanker et al. 2014;Range and Wei 2016), and foxQ2 genes also play central roles in anterior development in the insects D. melanogaster and Tribolium castaneum, and the spider Parasteatoda tepidariorum (Lee and Frasch 2004;Kitzmann et al. 2017;Schacht et al. 2020). In the cnidarians Nematostella vectensis and Clytia hemisphaerica, foxQ2 genes are expressed in and required for the proper development of the aboral pole, in support for the-still debated-homology between the cnidarian aboral pole and the bilaterian anterior pole (Chevalier et al. 2006;Sinigaglia et al. 2013). In Spiralia, foxQ2 genes share an apical expression in the annelid P. dumerilii and the brachiopod T. transversa (Santagata et al. 2012;Marlow et al. 2014), which, as our study shows, is consistent with the majority of expression domains for foxQ2 genes observed in the annelid O. fusiformis. Therefore, foxQ2 genes appear to participate in ancient and broadly conserved genetic programs for apical and axial patterning in metazoans.
The consistent expression of foxQ2 genes in apical territories contrasts, however, with the complex phylogenetic pattern of evolution of this class of Fox genes. As our study reveals, expansions of the foxQ2 class are common in Spiralia, and specially in Annelida, with 11 paralogs in O. fusiformis, and 10 and 13 in the vestimentiferans L. luymesi and P. echinospica, respectively. While many of these paralogs probably emerged from species-specific expansions (e.g. foxQ2-1, foxQ2-4, foxQ2-5, foxQ2-6, and foxQ2-9 in O. fusiformis; and the cluster of vestimentiferan foxQ2 paralogs with an N-terminal EH-i-like motif) other paralogs might have a more ancient origin, tracing back to Annelida (e.g. foxQ2-1 in C. teleta and foxQ2-11 in O. fusiformis, both expressed maternally) and even Spiralia (e.g. foxQ2-2 in C. teleta and foxQ2-4 in M. yessoensis). Notably, however, nearly all bilaterian and cnidarian lineages retain at least a copy (two in the annelids O. fusiformis and D. gyrociliatus) of a subclass of foxQ2 genes with a C-terminal EH-i-like motif. In O. fusiformis, these two genes (foxQ2-2 and foxQ2-3) are controlled by the ERK1/2 signaling that establishes the axial polarity of the embryo and consequently show a narrow peak of expression at the animal pole at the time of the specification of the organizer cell at five hours post fertilization (Seudre et al. 2022). Based on these observations, we propose a model in which an ancestral foxQ2 gene containing a C-terminal EH-i-like motif originated in the last common ancestor to Cnidaria and Bilateria, followed by independent expansions in certain bilaterian and cnidarian lineages and fast divergence of the new copies, which tended to lose the EH-i-like motif. Despite these duplications, however, the new paralogs did not generally acquire radically different functions (i.e. neofunctionalization) but rather retain a role in aboral/apical/anterior development, and evolved temporal, spatial and regulative specialization of their expression (i.e. subfunctionalization), as observed for example in O. fusiformis and the hemichordate S. kowalevskii (Fritzenwanker et al. 2014). Further analyses of the gene regulatory network associated with foxQ2 genes in a broader range of metazoans, especially in spiralians with multiple copies, and role of the EH-i-like motif will contribute to uncover the evolutionary history and developmental consequences of foxQ2 expansions during animal diversification.
Altogether, our study informs the evolution, temporal, and spatial expression of the largely conserved Fox gene repertoire in the oweniid annelid O. fusiformis. Our findings provide valuable information to reconstruct the ancestral complement of these core developmental regulators in Spiralia and to continue unravelling the embryological role and contribution of this major family of transcription factors to the evolution of animal body plans.

Animal Husbandry and Embryo Collection
Adult specimens of Owenia fusiformis Delle Chiaje, 1844 were collected and shipped to London from the coast near the Station Biologique de Roscoff (France) during their reproductive season (May to July). In the lab, animals were kept in aquaria with mud and artificial seawater (ASW) at 15°C. In vitro fertilizations were conducted as previously described (Martín-Durán et al. 2016;Carrillo-Baltodano et al. 2021) and embryos were kept in glass bowls at 19°C until they reach the desired developmental stage. Larval stages were relaxed in 8% MgCl2 and all embryonic samples fixed GBE in 4% formaldehyde in sea water (or MgCl2, for larvae) for 1 h at room temperature. After washing the fixative with phosphate buffer saline (PBS) supplemented with 0.1% Tween-20, embryos and larvae were dehydrated to 100% methanol and stored at −20°C.

Identification of Forkhead Genes and Orthology Assignment
Candidate Fox genes for O. fusiformis were initially retrieved from the functional annotation of its genome assembly (European Nucleotide Archive, accession number: GSE184126) (Martín-Zamora et al. 2022), and the foxC sequence was obtained from a previous study (Martín-Durán et al. 2016 (Katoh and Standley 2013) with the L-INS-i strategy, trimming the forkhead domain as reported in the HMM profile. Maximum likelihood trees were constructed with IQTree v.2.2.0-beta (Minh et al. 2020) with automatic identification of the model of protein evolution and 1000 rapid bootstraps. Bayesian reconstructions in MrBayes v.3.2.7a (Ronquist and Huelsenbeck 2003) were also performed using the CAT model of protein evolution and two runs with four chains (one cold, three hot) run for 50,000,000 generations. Resulting trees (Supplementary figs. 2 and 3, Supplementary Material online) were visualized and edited with FigTree (https://github.com/ rambaut/figtree/).

Manual Curation and Genomic Structure of Fox Genes
Short Illumina RNA-seq reads from a developmental time course aligned to the genome and a de novo assembled transcriptome (Martín-Zamora et al. 2022) were used to manually curate the automatic annotation and produce full length transcripts for foxN1/4 (genes OFUSG25429.1 and OFUSG25430.1 were merged into a single gene), foxK (genes OFUSG10014.1 and OFUSG10013.1 were merged into a single gene), foxQ2-10 (genes OFUSG24642.1 and OFUSG24641.1 were merged into a single gene) and the Fox orphan-2 gene (genes OFUSG25458.1 and OFUSG25457.1 were merged into a single gene). The exon and intron positions, as well as the chromosomal location for each Fox gene were determined based on the gene annotation of the genome assembly of O. fusiformis (Martín-Zamora et al. 2022) using the Integrative Genomics Viewer (Robinson et al. 2011). The position of the Forkhead domain within each Fox gene was determined using the Conserved Domain Database (CDD/SPARCLE) (Lu et al. 2020). The genomic architecture of the Fox genes was visualized using the online software IBS (Liu et al. 2015) and transferred to Illustrator Creative Cloud 2022 (Adobe).

Phylogenetic Analysis of the foxQ2 Family
To study the evolution of the FoxQ2 family in Spiralia, we retrieved sequences from spiralians (C. teleta, C. gigas, D. gyrociliatus, E. andrei, H. robusta, L. gigantea, L. luymesi, M. yessoensis, P. dumerilii, P. echinospica, S. benedicti), ecdysozoans (D. melanogaster, N. vetripennis, P. humanus, T. castaneum), deuterostomes (Oryzias latipes, S. kowalevskii, S. purpuratus, B. floridae) and cnidarians (C. hemisphaerica, Hydra vulgaris, N. vectensis) from published genomes, transcriptomes and databases. For all genes, membership to the foxQ2 family was manually confirmed using the conserved domain database CDD/SPARCLE (Lu et al. 2020) and identifying the presence of a full FoxQ2 specific domain (Accession number CD20035) within the sequences. Some genes previously misannotated and assigned to other Fox gene classes were renamed as foxQ2. To identify the EH-I like domain in the foxQ2 genes, we generated a EH-1-like motif position-specific scoring matrix using STREME v 5.4.1 (Bailey 2021), with a motif width of 7 to 10 amino acid and retaining motifs with P-value < 0.05, and confirmed the identified motif by comparison with results obtained by (Yaklichkin et al. 2007) (fig. 4B). The position and the presence of the EH-like motifs within the FoxQ2 sequences were determined by inputting the EH-i-like domain matrix into FIMO v 5.4.1 (Grant et al. 2011), matching motifs with a q-value < 0.05. For phylogenetic analyses, multiple protein alignments were performed with MAFFT v.7 as explained above and trees were reconstructed from a set of sequences selected with a Q.insect amino acid replacement matrix (Minh et al. 2021) to account for transition rates, the gamma distribution with four categories (G4) (Yang 1994) to describe sites evolution rates, and an optimization of amino acid frequencies using maximum likelihood in IQ-TREE v.2.1.2 (Minh et al. 2020). A thousand ultrafast bootstraps were used to extract branch support values and posterior probabilities were