Analysis of the P. lividus sea urchin genome highlights contrasting trends of genomic and regulatory evolution in deuterostomes

Sea urchins are emblematic models in developmental biology and display several characteristics that set them apart from other deuterostomes. To uncover the genomic cues that may underlie these speciﬁcities, we generated a chromosome-scale genome assembly for the sea urchin Paracentrotus lividus and an extensive gene expression and epigenetic proﬁles of its embryonic development. We found that, unlike vertebrates, sea urchins retained ancestral chromosomal linkages but underwent very fast intrachromosomal gene order mixing. We identiﬁed a burst of gene duplication in the echinoid lineage and showed that some of these expanded genes have been recruited in novel structures (water vascular system, Aristotle’s lantern, and skeletogenic micromere lineage). Finally, we identiﬁed gene-regulatory modules conserved between sea urchins and chordates. Our results suggest that gene-regulatory networks controlling development can be conserved despite extensive gene order rearrangement.


INTRODUCTION
Sea urchins are benthic marine animals that have attracted the curiosity of scientists since antiquity. 1They are one of the most amenable model systems in developmental biology because they produce abundant gametes and transparent embryos, and they have enabled major discoveries, such as the chromosomal nature of heredity, regulators of the cell cycle, and gene-regulatory networks controlling development. 2,3ea urchins belong to the clade of echinoderms, whose body plan displays an intriguing pentaradial symmetry, acquired secondarily during metamorphosis of their bilaterally symmetrical larvae. 4,5In the tree of animals, echinoderms together with hemichordates form the Ambulacraria clade, the sister group to chordates. 6Beyond their original symmetry, echinoderms possess intriguing novelties, such as their calcium carbonate endoskeleton and their water vascular system running through five canals and associated tube feet. 7Among echinoderms, sea urchins also stand out with some remarkable embryological and adult morphological characteristics.During embryogenesis, they proceed with early specification of a micromere cell lineage that gives rise to the embryonic skeleton, a mechanism absent in other echinoderms. 8,9As adults, they also possess an intricate calcified masticatory organ: Aristotle's lantern.
The first sequenced sea urchin genome, that of the Pacific sea urchin Strongylocentrotus purpuratus, revealed the extensive conservation of its gene repertoire with that of the human genome 10 and served as a support for the elucidation of generegulatory networks at play during development. 11,12However, while comparative genomics has helped us understand functional differences across different vertebrate model species, the genomic differences between distinct echinoderms have not been thoroughly examined. 13A possible clue could come from the Hox genes, a classic locus linking genomes and body plan establishment. 146][17] The subsequent discovery of an intact Hox cluster in the sea star 18 and sea cucumber 19 suggested, however, that this rearrangement is not the primary reason for the novel body plan of echinoderms. 20The sea urchin Lytechinus variegatus has revealed the same rearranged Hox cluster in this species as in S. purpuratus, 21,22 which could indicate that sea urchins are more prone to genomic rearrangements than other echinoderms.
][25] Ancestral linkage groups (ALGs) also appear to have been reshuffled after whole-genome duplication at the origin of vertebrates that underwent extensive chromosomal rearrangements. 26Other deuterostomes, such as tunicates, underwent very fast genomic evolution, 27 while others, such as hemichordates, retained more ancestral traits. 28Particularly, ALGs corresponding to maintained chromosomal units over time appeared at the origin of animals and were conserved at long evolutionary distances in lineages such as cnidarians, molluscs, or other spiralians with very limited fusions or rearrangements. 29,30Little is known, however, about the genomic rearrangements that took place in the sea urchin lineage and might underlie the evolution of sea urchin-specific novelties in embryogenesis and body plan.
Here, we report the generation and analysis of a chromosomescale assembly of Paracentrotus lividus, one of the main sea urchin species in the Mediterranean and North-East Atlantic coast and a reference species in embryology of the sea urchin.This species has been used since the 19th century in marine stations of France and Italy by prominent scientists, such as Ho ¨rstadius, Hertwig, Boveri, and Driesch, whose work has led to the foundation of key concepts in embryology and cell biology. 31,32or instance, Theodor Boveri demonstrated the requirement for a complete chromosome set in each blastomere for development, 33 and Hans Driesch discovered the astonishing capacity of isolated blastomeres from P. lividus to develop into smaller but harmoniously patterned larvae. 34This species is also widely consumed around the Mediterranean as a delicacy since antiquity, as testified by archaeological work conducted in Pompei. 35. lividus is also a keystone species in benthic communities because its grazing activity can control the development of algal populations. 36This genomic resource will thus promote population genomics studies on this species, which will improve the management of this ecologically and economically important species. 37n this study, we investigated how conserved the genomic and regulatory architecture is between P. lividus, other sea urchins, and chordates.We compared the gene order and chromosome organization of P. lividus with that of S. purpuratus and L. variegatus, which diverged $60 mya. 38We profiled genomewide chromatin accessibility throughout embryonic development.We investigated the expression and regulation of genes that originated at the origin of the echinoids clade.We identified a novel family of Pmar1-related transcriptional regulators specific to euechinoids, expressed in the micromere lineage and capable of converting any cell of the embryo into skeletogenic precursors.By integrating genomic and regulatory datasets, we demonstrated how regulatory changes could be associated with the origin of the novel body plan of urchins and other echinoderms.

RESULTS
The genome and developmental regulatory landscape of P. lividus We sequenced the genome of a single male individual of P. lividus using 2303 Illumina sequencing (Table S1).The genome was initially assembled using stringent parameters, and haploid copies were filtered out because the 2.9% polymorphism caused both haplotypes to be present in the initial contigs (Figure S1A).The assembly was then improved using long-read PacBio sequencing for gap closing, followed by scaffolding using proximity ligation (Chicago and Hi-C).The draft assembly of P. lividus amounts to 927.4 Mb with 5% of gaps and only 1.7% missing and 0.7% duplicated BUSCO genes (STAR Methods).Half of the assembly is comprised of 8 scaffolds larger than 41 Mb (N50), and there are 18 main scaffolds (>2 Mb) that represent 88.3% of the assembly (Figure 1A; Table S2).These 18 scaffolds correspond to the chromosome number observed in cytogenetic studies 39 (Figure 1B).Other sea urchin model species possess 21 chromosomes per haploid genome (S. purpuratus) 40 and 19 (L.variegatus). 21While most of the repetitive fraction of vertebrate genomes is usually made of retrotransposons, DNA transposons are dominant in the sea urchin, similar to other deuterostomes, such as amphioxi 24,25 (Figure S1B).In addition, transcriptomes for 17 embryonic stages, from the fertilized egg to the pluteus stage, and 5 adult tissues were generated (Table S3).We annotated 32,041 genes, of which 80.1% have similarity to other species, and 92.1% are expressed (TPM > 1).We also annotated 5,087 high-confidence long non-coding RNAs (lncRNAs).By using gene expression clustering, we determined that 4,965 lncRNAs show stage-specific expression and, hence, play a potential role in the regulation of development (Figures S2C and S2D).
To gain further insight into their developmental regulatory program, we applied assay for transposase-accessible chromatin sequencing (ATAC-seq) to 6 stages from the 16-cell stage to the pluteus larva stage (Figure S3A; Table S4).Genomic regions with higher-than-background chromatin accessibility were labeled as open chromatin regions (OCRs) and classified depending on their location (Figure 1D) and their activity profile during development (Figure 1C).Among 64,701 OCRs distinct from transcription start sites, we found 40% ''dynamic'' regions (peaks) that were specific to one or two successive developmental stages and 57% constitutive regions present in three or more stages (Figures 1C and 1G).To further elucidate the regulatory activity in these regions, we applied Cut&Tag targeting the H3K27ac histone modification that is associated with active enhancers 41,42 (Figures S3D, S3E, and S3G).Many OCRs also exhibit a level of transcription consistent with regulatory activity (Figure S3D). 43Our data efficiently recovered the activity of previously well-characterized regulatory elements (Table S5), such as, for example, the proximal and intronic enhancers of the nodal gene 44 (Figure 1F).We compared the distribution of OCR locations with that of amphioxi and zebrafish, for which developmental ATAC-seq is available 25 (Figures 1D and 1E).It has been proposed that OCRs that play an important regulatory role usually lie either in the vicinity of the promoter (proximal, <5 kb) or within introns (gene body). 45In vertebrates, a large fraction of these elements is located more distally than in amphioxi, which has been associated with more complex gene-regulatory processes (Figure 1D). 25 P. lividus shows a higher proportion of distal OCRs compared with the cephalochordate amphioxi and a lower proportion compared with zebrafish (Figures 1C and 1D).
We noticed that 16-cell-stage embryos (27,501 peaks) as well as pluteus-stage larvae (39,623 peaks) exhibit fewer OCRs than the other sampled stages.The blastula (137,938 peaks) and gastrula (147,188 peaks) stages show many stage-specific active elements, which are likely involved in the extensive gene-regulatory events that control embryo patterning and cell fate specification at these stages (Figure 1C; Table S5).In early 16-cell-stage and pluteus larva, the lower number of peaks could be interpreted either as a generally more relaxed chromatin state, potentially related to the absence of transcription before zygotic genome activation, or a more condensed state (Figure S3C).We hypothesize that the 16-cell stage possesses large open chromatin domains, as seen in other species, 46 while the larval stage (pluteus) likely undergoes a reduction of open chromatin domains, as seen in other lineages, such as annelids 47 (Figure 1C).Similar patterns are also recovered with an H3K27ac enrichment signal (Figure S3G).The sea urchin regula-tory landscape therefore generally resembles the architecture observed in vertebrates or amphioxi in terms of enhancer location and developmental dynamics. 48posite trends of genomic architectural changes in sea urchins and vertebrates Many animal genomes retain chromosomal linkage over time as sets of genes stay localized on homologous chromosomes, even in distantly related species. 49Particularly, comparison of highquality genomes of vertebrates, cephalochordates, and molluscs has revealed the existence of 17 ALGs ancestral to chordates. 26However, the occurrence of chromosomal fusions in some lineages suggests that more genomes and lineages should be taken into account to infer the ancestral bilaterian complement of ALGs.We examined the distribution of these ALGs in P. lividus as well as two other available sea urchin genomes (Figures 2A and 2B). 21,26Such representations pinpoint pairs of chromosomes or scaffolds that exhibit significant mutual ortholog enrichment (Fisher's exact test, p < 0.05).Our comparisons indicate that some of the ALGs actually fused in the vertebrate and amphioxus lineages and suggest the existence of 23 ALGs for bilaterians, in agreement with recent studies 29 (Figures S4A and S4B).Our comparisons of the sea urchin chromosomes further support the validity of these 23 ALGs as independent genomic units and reveal a missing linkage group that was not detected previously (dubbed ''ALG R''), likely dispersed across multiple chromosomes in chordates, which indicates 24 ALGs (Figure 2).This ALG R is merged with ALG Q on P. lividus chromosome 4 (chr4) but remained intact in S. purpuratus and L. variegatus (Figures S4C and S4D).We found that 13 of 18 chromosomes of P. lividus descend from a single ALG, four are derived from the fusion of two ALGs, and one (chr2) receives contributions from three more (Figure 2).
Despite the conservation of ancestral linkages in P. lividus, the distribution of ALGs in the chromosomes of L. variegatus and S. purpuratus indicates some species-specific chromosomal fusions 10,21 (Figure 2B).The fusions observed in P. lividus are absent in the two other species, but, conversely, some independent fusions took place, particularly in L. variegatus, where the two largest chromosomes (chr1 and chr2) derived from a mixture of distinct ALGs.Uniquely, ALG A1 and A2 are fused in P. lividus but not in the other sea urchin species, a fusion event that also took place in amphioxi. 26,50The higher number of chromosomes of S. purpuratus reflects the occurrence of only two fusions, the rest of the ALGs being represented as a single chromosome.In contrast, some chromosomal fusions are more ancient; two fusions are shared by all three urchin species: fusion of the three ALGs C1, B2, and E and fusion of ALGs B3 and J1.
Unlike vertebrates, sea urchins did not undergo a sizable rearrangement of units derived from ancestral linkages in most of their chromosomes, even when secondary fusions are observed (Figures 2C, S4C, and S4D).Incidentally, the distinct sea urchin lineages under scrutiny here have diverged since 48-68 mya, less than for human and mouse. 38In contrast to this echinoid chromosomal stability, we observed an extensive reshuffling of the microsyntenic intrachromosomal gene order, which results in the absence of an observable ''colinear'' gene order visible as linear segments across pairs of homologous chromosomes, as seen when comparing human and mouse genomes (Figures 3B and 3C).To quantify the rate at which gene collinearity is eroded, we compared the retention of microsynteny with the divergence time for selected sea urchin and vertebrate species (Figure 3A) 51 and showed in this way that intrachromosomal gene order appears to evolve at a much slower pace in vertebrates than in sea urchins (Figure 3A).
Consequently, the gene order in the genomes of sea urchin appears to evolve following a trend distinct from that of verte-brates: the rate of interchromosomal rearrangement appears to be very low, while the rate of intrachromosomal gene order change appears to occur much faster.To this end, we performed gene family reconstruction (Table S6) and applied phylogenetic reconciliation to detect duplication events (Figure 4A).Sea urchins, like other echinoderms, do not show a particularly increased occurrence of gene losses in their genomes, such as that observed in the tunicate lineage 53 (Figure 4A).Some losses, however, are shared by all three sea urchin species (Table S7); for instance, some members of the transforming growth factor b (TGF-b) signaling pathway, like BAMBI or BMP9, present in other echinoderms, like the sea cucumber. 54nterestingly, we noticed a burst of gene duplication events in the echinoid ancestor with the second-highest number of duplicated gene families in deuterostomes after the origin of vertebrates and its whole-genome duplications. 55However, in the sea urchin, these duplicates are located in close genomic proximity, present on the same chromosomes (49%) and at close distances (36% closer to 100 kb), indicating a common origin by tandem duplication rather than by large-scale genomic duplications, as reported recently for cephalochordates. 56The gene duplicates show enrichment of Gene Ontology (GO) terms associated with membrane transport and the circulatory system, with, for instance, multiple ABC transporters and solute carriers encoding genes (Figure 4C), which possibly suggests that these genes could play a role in the functioning of the water vascular system (Figure 4C).In the most extreme cases, some gene families underwent large expansions in the sea urchin lineages, as detected by hypergeometric tests (p < 0.01; Table S8): the glycoprotein Kirrel, which plays a role in sea urchin skeletogenesis; 57 the SLC16 transporter family, which is enriched in pigment cells; 58 a class of GPCRs (GPCR135) that acts as a putative neuropeptide receptor; 59 and the muscarinic acetylcholine (ACM) receptors.Some of these duplicates seem to have acquired specific expression profiles or domains during development (Figures S5C and S5D).

The expression of recently evolved sea urchin genes
To further evaluate the expression of genes gained and duplicated in the echinoid lineages, we classified genes according to their expression profiles using network-based clustering 60 (Figure S5A).We used gene family reconstruction to determine whether genes that originated and duplicated at different phylogenetic nodes showed a particular enrichment in some of these expression clusters (Figure 4B).Novel echinoid genes are particularly enriched during early embryonic stages or in adult structures, such as tube feet or the body wall (comprising the water vascular system and skeleton).We also found such increased expression of novel genes in Aristotle's lantern, a calcified buccal apparatus specific to echinoids.We similarly observed preferred expression of genes duplicated in the echinoid lineage during the earliest embryonic stages and among the maternally expressed genes, suggesting that some of them could be involved in cell lineage specification (see below the example of pmar1) (Figure 4B).We then examined the putative regulatory elements (OCRs) that are associated with genes duplicated in echinoids; we identified 3.45 elements per gene (on average) for echinoid duplicates in contrast with 4.13 elements for genes arising in older duplications (Figure 4E).These elements are closer to promoters and less conserved in sequence than the elements associated with older duplicates or even the single-copy genes (Figure 4F).At the gene expression level, very recent duplicates show a higher organ specificity than other genes according to the tau estimate (Figure 4D), indicating that novel, fast-evolving regulatory elements are driving their expression in new organs and stages. 61We assessed whether some transcription factor binding sites (TFBSs) are enriched in the OCRs associated with sea urchin gene duplicates and recovered significant enrichments for 22 TFBSs (hypergeometric test, p < 0.01; Table S9).This list comprises five nuclear receptors, including homologs of HNF4, thyroid receptor, and COUP-TF, as well as other factors known for controlling organogenesis (AP2, ELK, and ATF1), findings compatible with the organ-specific expression pattern recovered for a number of these duplicates (Figure 4D).
By studying the evolution, expression, and regulation of the sea urchin genes, we showed that newly evolved and newly duplicated genes are preferentially expressed in specific embryonic stages and anatomical structures.Some of these structures, such as Aristotle's lantern or tube feet, are novelties of Article ll OPEN ACCESS echinoids and echinoderms, respectively.This observation parallels similar findings regarding the younger transcriptome age of early developmental stages 62 or gene expression in the mollusc shell and radula, 63 but further studies will be necessary to explore their functional significance.
The origin of pmar/hbox12 genes To understand how tandem gene duplications underlie the emergence of novelties, we investigated the origin of the pmar1 genes that play a critical role in specifying micromeres and their skeletogenic primary mesenchyme cell (PMC) derivatives, a specific character of euechinoids. 64We identified 12 copies of genes related to the original pmar1/hbox12 sequence clustered within a 350-kb region in P. lividus (Figure 5A). 65In L. variegatus, a similar number of pmar-related genes has recently been described in two distinct genome locations of chr2. 21Moreover, eight copies of a pmar1-related gene have also been identified in the cidaroid Eucidaris tribuloides, the sister group of euechinoids. 66pmar1 genes do not have clear orthologs in other echinoderms, but the identification of phb as a pmar1-related gene in brittle stars (Ophiuroidea) suggested that pmar1 genes are divergent copies of an ancient class of paired-type homeobox genes in this clade. 8,66To clarify their respective relationships, we cataloged the multiple pmar copies found in sea urchins and copies of other pmar-related echinoderm genes, including the Pplx gene from brittle stars and the Phb gene from sea stars. 8We re-constructed a phylogenetic tree that confirms the relationship of pmar and phb genes but, surprisingly, also shows that the aprd genes found in amphioxi and spiralians are related to these genes. 8,67Our phylogeny suggests an independent expansion of pmar genes in multiple sea urchin lineages because the different pmar genes group into species-specific clades.Intriguingly, we identified a novel family of pmar-related genes that group in a distinct clade, showing a particularly fast rate of evolution, and dubbed them parent of pmar (pop) genes to reflect that they are sister genes of pmar (Figure S6E).The architecture of pmar/pop clusters appears distinct in all four sea urchins: in L. variegatus, pops and pmars are in distant loci, while in S. purpuratus, one locus contains 3 pmars and another 2 pmars and 2 pops.In E. tribuloides, eight pmar1-related genes are present in two different scaffolds, while all of these genes are closely located in P. lividus (Figure 5A).In the cidaroid E. tribuloides, only one set of pmar/pop-related genes is found, which could indicate that the duplication that gave rise to the ancestors of pop and pmar genes took place within the Euechinoidea lineage, followed by an independent expansion of pmar genes in the lineages leading to each sea urchin species.
Because pop genes seem to constitute a distinct clade of pmar1-related genes, we sought to determine whether their expression and function are the same as other pmar1 genes.In P. lividus, the pop1 gene is expressed maternally and ubiquitously, while the two other pop genes show an expression spatially Article ll OPEN ACCESS restricted to the micromere lineage in the late cleavage and early blastula stage, consistent with these genes being, like pmar1, involved in specification of the micromeres and PMC lineage (Figures 5D and 5E).Indeed, overexpression of these pop genes caused massive delamination of PMC-like mesenchymal cells that popped up from the vegetal pole at the onset of gastrulation and that was accompanied by ectopic expression of the PMCspecific marker genes delta and alx1 (Figure 5F).This phenotype is identical to the phenotype caused by overexpression of pmar1, reinforcing the idea that pop genes also act to specify the PMC lineage. 64These observations illustrate how pervasive gene duplication and rearranged gene order between the three sea urchin models are associated with novel cell lineage specification mechanisms.

Conservation and divergence of gene expression modules across deuterostomes
The dynamic gene order across sea urchin species contrasts with their highly conserved embryonic development and corresponding gene expression. 68To understand how the gene-regulatory program is affected by these gene order changes, we assessed the extent of the evolutionary conservation of sea urchin regulatory programs by comparing clusters of genes with similar temporal expression profiles in P. lividus, S. purpuratus, and the cephalochordate Branchiostoma lanceolatum. 69We found significant pairwise enrichment of genes belonging to identical gene families in clusters of co-expressed genes, either between sea urchins or when comparing sea urchins (Figure 6A) and chordates (Figure 6B).This enrichment indicates conservation of genes involved in temporal gene expression modules, as observed previously within chordates 25 and insects. 70Interestingly, the patterns of conservation at short and long divergence times (Figures 6A and 6B) show remarkable similarities, with most highly conserved modules being the stage-specific ones and the ones active during mid-to late development.We noticed that evolutionarily conserved modules can involve gene expression at the same developmental stage (homochronic) or shifted in their expression timing (heterochronic) at short and long evolutionary distances (Figures 6A and 6B).Despite the difference between the pluteus larva and the amphioxus larva, we observed several clusters that share a significant homologous gene set.For instance, the cluster (Pliv18) that shares the highest number of homologous genes appears to be enriched in GO terms associated with sensory perception and neuronal function (ion and amino acid transport) (Figure 6C). 71An examination of TFBSs enriched in OCRs related to genes belonging to these clusters revealed multiple factors associated with circadian regulation of gene expression, such as ARNT (BMAL) and CLOCK, which constitute the ''core clock;'' associated regulators such as USF; downstream circadian effectors such as ATF and CREB; plus genes such as SREB, which is related to nutrient-related circadian adjustments (Figure 6D).These results suggest an evolutionarily conserved mechanism of circadian-based activation of neural activity at the end of embryonic development and the onset of larval life. 72eyond gene expression, we focused on the dynamics of the OCRs during development in sea urchins to determine the key regulatory steps and how evolutionarily conserved they are.We identified conserved non-coding regions in P. lividus by performing an alignment of all three available sea urchin genomes and applying a statistical model to infer non-coding conservation (Phastcons).We found that 39% of OCRs with a putative regulatory role (non-repetitive, non-transcription start site [TSS]) overlapped with the 81,142 evolutionarily conserved regions that otherwise represent 1.94% of the genome (17.9 Mb) (Fisher's exact test, p < 1eÀ9).We further assessed the sequence conservation of putative regulatory elements across development by assessing the evolutionary conservation of OCRs active at specific stages (Figure 1C). 73We identified the OCRs specific to the blastula, gastrula, and prism stages as the most conserved in sequence (Figure 7A).
We also compared the transcriptomic distance (Jensen-Shannon) between embryonic stages in sea urchin and deuterostome species to determine which embryonic stages exhibit the closest transcriptomic proximity 25 (Figures 7B and S6A).The early blastula, gastrula, and prism stages that show the highest sequence conservation are also the ones that have the lowest transcriptomic divergence in the comparison of the two sea urchins (P.lividus and S. purpuratus) and the comparison with the cephalochordate amphioxus (Figure 7A). 25,47,74However, sea urchins appear to show a higher level of conservation at earlier stages (Figures 7B and S7A).This conservation can be explained by an earlier onset of major cell lineage specification at the early blastula stage. 12

Regulatory landmarks during sea urchin development
To determine which TFs control successive phases of sea urchin development, we employed a footprinting approach to detect the most likely occurrences of TF binding events in the ATACseq signal at successive stages 75 (Figure 7).Then, by using reconstructed gene families, we identified 815 TFs in P. lividus and assigned 568 of them to TFBSs derived from the JASPAR database, which were later employed to analyze the footprint of the ATAC-seq signal (Figure 7D; Table S10) that we compared with their gene expression across developmental stages (Figure S7A).To ascertain the validity of the regulatory interactions deduced from our chromatin profiling data, we also Article ll OPEN ACCESS reconstructed, for each stage, the global gene-regulatory network (GRN) by leveraging chromatin accessibility in OCRs, the assignment of OCRs to genes, and gene expression. 76The resulting networks represent genes as vertices and regulatory interactions as edges and can be interrogated to highlight a specific subnetwork.This analysis recovered, for instance, the endomesoderm specification GRN (Figure 7C) 11 and the skeletogenic GRN (Figure S7E) 12 described previously for S. purpuratus, validating the idea that the OCRs we identified are indeed cis-regulatory modules that bind key TFs of these GRNs.While a small number of loci appear to be transcribed during cleavage stages (Figure S7A), we found that the onset of large-scale zygotic genome activation is taking place at the end of cleavage, in the early blastula stage, as pointed out previously 77 (Figures 7D and S8B).Gene expression clustering (WGCNA) identifies small sets of genes specifically expressed in each of the cleavage stages, such as G-protein receptors (Oprx1, 16-cell), ionic channels (Sc5a2, 60-cell), Toll-like receptors (2-cell), or homeodomain genes (Figure S6D).However, zygotic expression of TFs appears to only start at the 60-cell to late cleavage stages and gradually increases to mid-blastula stage (Figure S6E).Transcripts encoding TFs are present as maternal messages and remain detectable in later stages.At the 16-cell stage, most OCRs correspond to constitutively active and not stage-specific elements (Figure 1C), with a limited ATAC signal at TSS, confirming that zygotic expression is not yet taking place (Figure S3C).Accordingly, the TFBSs enriched at the 16-cell stage (Figures 7D and 7E) are also enriched in the early blastula stage, with few of corresponding TFs showing detectable expression at these stages (Figure S7A), which could indicate that some chromatin regions are present in a relatively open configuration before zygotic genome activation.
Conversely, the early blastula stage is characterized by the enrichment in TFBSs for maternally expressed TFs such as Ets4 and SoxB1 (Figures 7D and S7B).9][80] Consistent with this, we detected Ets4 binding motifs in an OCR located in the vicinity of the He2 gene (Figure 7E).PITX2, and MIX-1), nuclear receptor family (HNF4), and bZIP family (CREB1 and Myc) also have their binding sites enriched at the 16-cell and early blastula stages and may likely participate in activation of the zygotic genome (Figure 1D).Later developmental stages showed increased TF motifs for developmental regulators expected to be active by mid-development, such as Sox, Otx, or Gsc factors.For instance, ATACseq footprints of SOXE, expressed in the left coelomic pouch that gives rise to the rudiment, and PITX2, expressed in the right coelomic pouch, are also enriched at the prism stage, consistent with establishment of the left-right asymmetry of the larva at this stage. 81Accessibility to FOX binding sites is restricted to the gastrula stage and subsequent stages.At the prism and then pluteus larva stage, we detected binding motifs for a set of circadian-controlled factors, such as CLOCK and ARNT (see above), HNF factor and bZIP factor, as well as CEBP-related factors.The pluteus larva stage similarly shows fewer OCRs (Figure 1) and limited enrichment of TFBSs.The most characteristic larva-specific factors are Ari5b, which plays a role in organ growth in vertebrates, and Fuk-1, which bears similarity to NFATC2, which a putative role in immune system activation in vertebrates.The role of TFs at the successive developmental stages is confirmed by analysis of their influence in the GRNs reconstructed at each stage, as estimated by their prominence in the network architecture or centrality (Figure 7D).
Our examination of chromatin occupancy and TF binding informs on the mechanism of zygotic genome activation as well as on the subsequent sequence of regulatory events coordinating development, such as establishment of cell lineages, segregation of germ layers, and activation of circadian gene expression, and ultimately should help to reconstruct GRNs controlling developmental processes. 11

DISCUSSION
Sea urchins, as echinoderms, exhibit a unique and derived body plan with pentaradial symmetry, which has originally been related to a reshuffling of GRNs caused by an extensive gene order change.This reshuffling was exemplified by the original description of a rearranged Hox cluster in S. purpuratus. 15,16ere, we performed gene order comparison between the available sea urchin models and between sea urchin and amphioxi and demonstrated that the chromosomal architecture of sea urchins retained the ALGs inherited from a bilaterian ancestor. 26,29hese linkage groups were extensively reshuffled after wholegenome duplication in vertebrates or completely lost in tunicates and nematodes.Interestingly, while chromosomal architecture is very conserved, the local intrachromosomal gene order appears to evolve much faster in sea urchins than in vertebrates (Figure 3).While it could be tempting to link this observation with a faster evolution of gene regulation in sea urchins, most developmental and transcriptomics studies suggest conservation of GRNs and developmental mechanisms between sea urchin models, such as S. purpuratus, L. variegatus, and P. lividus. 68,82The relatively limited number of high-quality echinoderm and invertebrate genomes available at distinct evolutionary distances makes it difficult to determine whether the evolution gene order in sea urchins or in vertebrates constitutes the exception or the norm among bilaterians.This observation also suggests, more practically, that many loci will have a different organization and gene order between models such as P. lividus and S. purpuratus; for instance, the pmar1 gene locus (Figure 5A).
We also identified extensive echinoid-specific gene duplications and expansions and pinpointed that some novel echinoderm and echinoid structures, such as tube feet and Aristotle's lantern, captured the expression of many of these recently evolved genes.A similar enrichment of new genes has, for instance, been witnessed with the mollusc radula 63 or mollusc shell. 83Specialization of gene expression after duplication is similar to that observed for WGD duplicates in vertebrates or other lineage-specific duplicates in chordates. 25,50We also explored the origin of the pmar1 class of homeobox genes involved in cell lineage determination in an echinoid-specific fashion 64 and identified a novel family of pmar1-like homeobox genes, the pop genes, that arose in euechinoids and that share with pmar1 the ability to promote formation of the skeletogenic mesoderm.Remarkably, homeodomain genes acting at early cleavage stages, such as Paired or TALE genes, have shown the propensity to undergo lineage-specific duplication in other animal lineages, such as mammals 84 or spiralians, 85 suggesting a recurrent phenomenon during animal evolution.
Interestingly, by using RNA sequencing (RNA-seq) to monitor gene expression and ATAC-seq to explore the open chromatin domains, we were able to perform evolutionary comparisons of gene expression profiles and gene-regulatory modules across species.Our data suggest possible long-range conservation of regulatory logic between sea urchins and chordates (Figure 6), which is consistent with the idea of deeply conserved regulatory mechanisms controlling animal development, although experimental testing of these regulatory mechanisms will be required to validate this idea. 86,87We showed that early development (cleavage) was the most divergent stage at the gene expression level, even between two sea urchin species that split 75 mya, possibly suggesting that the dynamics of zygotic genome activation might vary across species (Figure 7). 38Finally, we confirmed that tightly regulated and large-scale transcription of the zygotic genome did not occur before the early blastula stage in P. lividus (Figure 7).
In sum, we identified conserved principles of gene regulation across deuterostomes in terms of genomic location and of activity of regulatory elements during development.We showed that these gene-regulatory modules are maintained despite the fast intrachromosomal gene order change in sea urchin genomes.The detected conservation of cis-regulatory modules complements previous observations comparing vertebrates and classic model systems like Drosophila or C. elegans and suggests that, despite their derived body plan, the presence of a larval stage and of mechanisms of early determination of embryonic cell lineages, echinoderms preserved such principles. 88mitations of the study Our analysis of gene order in sea urchins is limited by the extent of available genomes.With more echinoderm and invertebrate genomes, we could generalize observed trends of intrachromosomal and interchromosomal gene order changes and determine which, of sea urchins or vertebrates, is representative of the ancestral state.We also think that functional work targeting some of the candidate genes proposed in this study to be expressed preferentially in urchin novelties would be important to appreciate their role; for instance, in the water vascular system or Aristotle's lantern.Moreover, our characterization of the accessible and active genomic region is limited to the ''bulk'' whole embryo and could be expanded across germ layers and cell lineages using approaches such as single-cell RNA-seq and ATAC-seq.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: Article ll OPEN ACCESS Illumina sequencing libraries were prepared at Genoscope (Evry, France) using TruSeq kit for fragment library and Nextera Matepair Kit for mate-pair libraries.These libraries were sequenced on MiSeq, HiSeq2500 and HiSeq2000 instruments at Genoscope (see Table S1) and reads were quality trimmed using a custom script derived from the fastx toolkit.Long-insert Pacific bioscience libraries were prepared and 20 SMRT cells were sequenced on an RSII instrument at the Genotoul plateforme (INRA, Toulouse, France) with the P6C4 chemistry.Sequencing data was submitted to ENA under the master PRJEB25800 accession.

Genome assembly
Genome size was evaluated using the k-mer spectrum approach as described in. 113All occurences of 31-mer were counted in the fragment data using Jellyfish. 114The haploid genome size was estimated at $845Mb.The presence of two well defined coverage peaks, with the main one at half coverage was indicative of high levels of polymorphism (3%).
The Illumina data was assembled using Meraculous (v2.2.2.2) and using a k-mer of 91 and 'diploid_mode = 2' to enable the split assembly of both haplotypes in stringent conditions. 115Briefly, Meraculous performs a de Brujin graph contig assembly after excluding erroneous low frequency k-mers followed by subsequent greedy scaffolding rounds using mate-pair libraries and a minimum of 2 or 3 links to create a link between contigs.Finally, gaps in the assembly were closed using contig extension and paired-end information.We obtained a 1394Mb diploid assembly with a large fraction of residual gaps (20.7%) (Table S2).
To fill the gaps, and improve the Illumina assembly, we performed gap-filling and local reassembly using the PBJelly tool 91 (v15.8.24).PBJelly aligns the Pacbio reads to the assembly using the Blasr aligner designed to account for underlying read error profile, then collects reads surrounding and spanning gaps and scaffold extremity regions and assemble them with the ALLORA assembler relying on a OLC (overlap-layout-consens) design.Resulting error-prone assembled sequences are used to fill gaps and extend scaffolds.We used the parameters '-minMatch 8 -sdpTupleSize 8 -minPctIdentity 75 -bestn 1 -nCandidates 10 -maxScore À500' for blasr alignment.The resulting patched assembly has a reduced fraction of gaps (4.9%) and incorporates 1,475Mb of assembled DNA for a 1,551Mb total size (Table S2).
To generate a haploid reference genome from our diploid assembly, we employed the Haplomerger (v2, 20151124 build), which relies on a graph of reciprocal LASTZ alignments to extract the best path across haplotype scaffolds. 116Before reciprocal alignments, we masked repetitive regions in the diploid assembly with RepeatMasker using a custom repeat library (see below) and a custom scoring matrix obtained with the script lastz_D_Wrapper.pl.We then applied module A aimed at detecting and splitting misassembled scaffolds, and module B aimed at performing haplotype reconciliation.Finally, residual haplotype sequences smaller to be processed by module A (<5000bp) were screened using module G relying on residual similarity to resolved haplotigs.
To further extend the contiguity of our haploid reference, we used long-range contact information from Chicago and HiC library prepared by Dovetail Genomics (Santa Cruz, USA).The Chicago method relies on the reconstitution of synthetic chromatin in controlled condition followed by chromatin conformation capture and evenly distributed contact information in the 150 kb range.Library preparation is described in detail in. 92Alternatively, HiC chromatin conformation capture provides contact information at a broader range, and particularly helps reconstruct chromosomes which usually do not present much interactions in trans.Tissue was crosslinked in 1% PFA, and chromatin subsequently extracted, immobilised on SPRI beads, washed and digested with DpnII. 117fter end-labelling, proximity ligation was carried out using T4 DNA ligase and cross-linking reversed using Proteinase K, removed from the beads and the DNA fragments were purified again on SPRI beads.Sequencing library was constructed using the NEB Ultra library preparation kit (New England Biolabs, Ipswitch).Chicago and HiC libraries were sequenced for 476M and 210M paired-end reads in 2 3 150bp mode on a HiSeq4000 instrument.Chicago and HiC data were processed through two distinct runs of the HiRise scaffolder.The final assembly shows the following BUSCO statistics (v5.1.3)when using the Metazoa gene set: C:94.7% [S:94.0%,D:0.7%],F:3.6%,M:1.7%,n:954.

Annotation
Transcripts assembled with Trinity (see below) were aligned to the genome using GMAP (version of 2018-03-25). 118These alignments and the merged stringtie assemblies were leveraged using Mikado (v1.2.1) to generate a high-quality reference transcriptome. 97A set of curated splice-junctions generated from RNA-seq alignments using Portcullis (v1.0.2) was also provided to Mikado. 96Putative fusion transcripts were detected by Blast comparison against Swissprot and ORFs were annotated using Trans-decoder (Haas et al. 2008).Transcripts derived from the reference transcriptome were selected to train the Augustus de novo gene prediction tool (Stanke et al. 2006).Exon and intron positions derived from the mikado consensus transcriptome were converted into hints for Augustus gene prediction.
We annotated repetitive regions in the genome by constructing a repeat library using RepeatModeler (v1.0.11) that was subsequently used for masking with RepeatMasker (v4.0.7).Repeat landscape was subsequently inferred by computing Kimura 2-parameters distances with the consensus for each repeat category.Gene models with half or more of their exons overlapping at 50% with repeats were discarded, yielding 41717 filtered gene models.Alternative transcripts and UTRs were subsequently incorporated using the PASA pipeline (Haas et al. 2008).These gene models contain a total number of 4915 distinct PFAM domains.

Long non-coding RNAs
For lncRNA annotation, we used all RNA-seq data from P. lividus available in SRA archive (accessions: PRJEB10269, PRJNA392084, PRJNA376650, PRJNA288758, PRJNA264358, 4787.4M reads in total, 14.7 M reads per sample).After quality control with FastQC v0.11.6 119 120 we mapped the reads to the reference genome using HISAT2 v2.2.1 121 and obtained a transcriptome assembly for each sample using StringTie v2.1.4 122by providing our gene models as a reference annotation file.We obtained 175 GTF files corresponding to the individual assemblies which were merged to obtain a reference GTF using the merge option from StringTie.The merged GTF was used to identify the candidate lncRNAs using FEELnc software 99 in a three steps pipeline: i) filter out transcripts shorter than 200ncl, monoexonic and overlapping protein-coding exons, ii) discard transcripts with coding-potential and iii) classify the candidate lncRNAs according to their relative position in the genome.Given that there is no previous knowledge of lncRNAs in P. lividus, we used two different strategies to compute the coding-potential: i) we took a set of mRNAs and shuffled them while preserving 7-mer frequencies (shuffle approach), and ii) we provided as training set lncRNAs from a closely related species (S. purpuratus, reference approach).The number of candidate lncRNA transcripts obtained was 32,147 and 32,107 genes for the shuffle and reference approaches, respectively, with more than 99% overlap.After filtering the transcripts uniquely annotated by a single approach, we obtained a list of 32,001 candidate lncRNAs genes (56,259 transcripts).Subsequently, we estimated the abundance of the candidate lncRNAs in each SRA sample using the feature Counts function of the Rsubread package from Bioconductor 110 and we classified genes according to their expression range using the filter ByExpr function from EdgeR package 123 for R by setting the min.countparameter to 10.This allowed us to identify a high confidence lncRNAs set that includes genes that are expressed in most of the samples, consisting of 5,087 lncRNA genes.

Synteny and gene family reconstruction
To compute pairwise synteny comparisons, we used mutual-best-hits based on MMSeqs2 (MMSeqs2/12-113e3) comparisons between proteomes after selecting for the longest protein for each locus.After reindexing of gene coordinates, we used Fisher's exact test to determine mutual enrichment of orthologues between chromosomes or scaffolds.To estimate the relationship between gene order and divergence time, we estimated as the fraction of orthologues located in blocks of 2 or more consecutive genes in the same order with no more than one interspersed gene (Figure 3).Corresponding python and R codes upload to https://github.com/paracentrotus/urchinpaper.
We used Broccoli (v1.2) for gene family comparisons using the species specified in Table S5 , 101 and inferred clade specific gains and losses by comparing the content of each family with that of a reference species tree using the ete3 library. 102To assess duplication within gene families, families with less than 500 genes and at least 5 genes in 3 species were subjected to phylogenetic reconstruction: after alignment with Msaprobs (v0.9.7) 103 and alignment trimming using clipkit (v0.1, option -m gappy), 124 a tree was reconstructed using IQTREE (v2.1.1)assuming a LG4X + R model.Then, Generax (v1.2.2) was used to perform genes and species tree reconciliation and to detect duplication events in a maximum likelihood framework. 105

Gene expression analyses
We extracted RNA for successive embryonic stages and several organs (Table S3) using Trizol reagent (Invitrogen).Strand-specific RNA-seq libraries were prepared using the TruSeq RNA Library (Illumina) and sequenced in a 2x150bp layout with an average of 64.6M reads per sample.Reads were aligned to the genome using STAR (v2.5.2b) at an average rate of 65.64% of unique mapping. 107A transcriptome was assembled for each sample using Stringtie (v1.3.3b) 108and sample-specific transcriptome assemblies were merged using Taco. 109The reads from all the samples were also assembled de novo using Trinity. 94ene expression was quantified from reads aligned using STAR (v2.5.2b) and using featureCounts from the Subread package (v1.6.3) 110and counts converted in FPKM.We also evaluate the coverage of OCRs for both strands using featureCounts to evaluation their level of transcription on both strands.
We used WGCNA (v1.7.0) clustering for the full set of embryonic stages and organs.After filtering out genes with limited variance and counts, then, the 'softpower' parameters were estimated and set at 13, and clustering was run with a 'signed' network type. 60For all clustering analyses, FPKM were calculated and replicates if available were merged to obtain a single gene expression value.Then, we used mfuzz (v2.54.0) clustering to compare subsets of 8 embryonic stages in sea urchin, amphioxus and zebrafish using datasets for the two later as described in. 25After filtering genes with low expression or limited variability, expression values were normalised, the fuzzifier parameters (m) was estimated and the optimal number of clusters was determined by computing with minimal distance between cluster centroids (Dmin) for various numbers of clusters.To compare the evolutionary conservation of gene expression modules, a hypergeometric test was performed on the number of genes belonging to distinct gene families shared between each pair of alternate gene expression modules (accounting for possible many-to-many paralogy relationships between individual genes).Gene ontology enrichments were computed using the topGO package using terms transferring from PFAM annotation and Swissprot best hits (evalue 1e-10). 125Comparisons of transcriptomic distances across embryonic stages and were computing as Jensen-Shannon divergence and single-copy orthologues inferred using Broccoli as performed in. 25S. purpuratus data analyzed corresponds to PRJNA81157 and Branchiostoma floridae to PRJNA416866.Table S6.Genomes and proteomes used for gene family reconstruction, related to Figure 4.

Figure 1 .
Figure 1.Genome organization and regulatory landscape of the sea urchin P. lividus (A) HiC contact map of the P. lividus assembly, with the 18 longest scaffolds of higher contact density corresponding to putative chromosomes highlighted.(B) Density of annotated genes (color scale) and repeated elements (ridge plot on the right) with a picture of an adult P. lividus (C.G.).(C) Classification and number of OCRs for the different stages.(D) Number of OCRs located at the transcription start site (TSS), in the proximal region (<5 kb upstream of the TSS), in the gene body, and/or in the distal region (>5 kb of the TSS) in three deuterostome species.(E) Cumulative distance to TSSs of OCRs in the same three species.(F) RNA-seq (red) and ATAC-seq (blue) signals in the region of the nodal gene, showing two well-characterized CREs in the proximal and intronic regions of this gene.

BFigure 2 .
Figure 2. Evolution of sea urchin chromosomal architecture (A and B) Oxford plots visualizing the respective positions of orthologs inferred by reciprocal best blast in the sea urchin P. lividus; the cephalochordate Branchiostoma floridae, where ALG A1 and A2 fused (A); and the mollusc Pecten maximus, where several other ALGs fused (B).(B) Dots located in pairs of chromosomes showing a significant mutual enrichment of orthologs (Fisher's exact test, p < 0.05) are colored by ALG assignment, while others are colored in gray.Axis values represent gene indexes.(C) Synteny between chromosomes of all three available echinoid genomes (P.lividus, L. variegatus, and S. purpuratus), colored by ALG.

Figure 4 .8
Figure 4. Gene and organismal novelties in sea urchins (A) Gene family gains (green), losses (blue), and duplications (orange) on a phylogenetic tree of deuterostomes.(B) Enrichment of genes originated (top) and duplicated (bottom) at different phylogenetic nodes in WGCNA clusters of syn-expressed genes using a hypergeometric test.(C) GO terms enriched in genes duplicated at the echinoid nodes for Biological Process (BP) and Molecular Function (MF) categories.(D-F) For genes duplicated at distinct nodes, we evaluated (D) gene expression tissue/stage specificity (tau), (E) distance to the TSS, and (F) Phastcons conservation score in OCRs associated with the corresponding genes.

Figure 5 .
Figure 5.The evolution of pmar/hbox12 genes in echinoids (A) Genomic organization of the pmar and pop loci in all three echinoid genomes.(B) Regulatory landscape with RNA-seq (red) and ATAC-seq signal (blue) and OCRs in P. lividus.(C) Phylogeny of pmar-related paired homeobox genes using the homeobox residues (IQTREE LG4X+R model).(D) Expression of pmar and pop genes.(E) In situ expression of pmar and pop genes.pop1 is expressed maternally and ubiquitously, while pop2 and pop3 are expressed in the PMC precursors.Scale bar.30 mm. (F) Phenotypes caused by overexpression of pmar or pop genes.Overexpression of pmar1, pop1, pop2, pop3, or pop2 fused to the repressor domain of Engrailed causes massive production of PMC-like mesenchymal cells and ectopic expression of PMC marker genes such as Delta and alx1.Inset: ventral view.vv, vegetal pole view; DIC, differential interference contrast.

Figure 6 .Figure 7 .
Figure 6.Conservation of gene expression modules in deuterostomes(A and B) Gene content conservation between cluster syn-expressed genes (mfuzz) during the development of P. lividus and the sea urchin S. purpuratus (A) or the cephalochordate B. lanceolatum (B), estimated using a hypergeometric test.The arrowheads underneath indicate whether the pair of clusters is homochronic (black filled) or heterochronic (white filled).The side heatmaps indicate average expression for each mfuzz cluster as normalized Z score.(C) GO terms overrepresented in cluster 18, which shows the highest conservation of ortholog content.(D) TF motif enrichment in peaks associated with genes belonging to each cluster in P. lividus, computed using a hypergeometric test (p < 0.01).Only clusters with at least one significant TF motif (p < 0.05) are shown.

Figure S1 .
Figure S1.The genome of P. lividus, related to Figure 1.(A) K-mer spectrum for k=31 showing the two peaks characteristics of highly heterozygous genomes.(B) Repeat landscape showing the distribution of divergence (approximation for age) for distinct class of mobile and repetitive elements.(C) Size of different categories of gene features highlighting the comparatively large size of UTRs and introns in P. lividus.

Figure S3 .
Figure S3.Chromatin profiling of P. lividus embryonic stages, related to Figure 1.(A) PCA based on normalised counts of ATAC-seq data.(B) Validation of OCR classification (Figure 1C) by visualising normalised ATAC-seq read coverage for each stage in the unified OCR for regions showing dynamic activity during development.OCRs are clustered according to their coverage signal across the successive developmental stages to verify the accuracy of the intersection-based classification, which is summarised under the heatmap.The top-row of the heatmap and the legend on the right indicates at which stages unified OCR are overlapping with OCR.A high coverage (red) indicates that the element is active while no coverage (blue) supports limited accessibility.(C) Profile of ATAC-seq signal around predicted OCR (peaks) for distinct location of OCR.(D) Distribution of coverage of ATAC-seq, Cut-and-tag K27Ac and transcription (RNA) for distinct OCR location inferred from ATAC-seq data, as well as a similar population of randomised regions of the same size distribution.(F) Transcription in OCR indicating occurence of bidirectional transcription for some of them.(G) Classification of OCR activity defined by 75% upper coverage quartile of OCR regions across the four stages for which we performed Cut-and-tag.

Figure S4 .
Figure S4.Synteny comparison between representative bilaterians and sea urchins, related to Figure 2. 'Oxford' dotplots representing the genomic coordinates of orthologues in pairs of genomes.Dots located in pairs of chromosomes showing a significant mutual enrichment of orthologues (t-test <0.05) are colored by ancestral linkage group assignment (ALG) while others are colored in grey.(A) Comparison of amphioxus (B.floridae) and purple sea urchin (S. purpuratus) (B) amphioxus and a bivalve mollusc, the sea scallop (Pecten maximus) (C) the sea urchins P. lividus and L. variegatus and (D) the sea urchins P. lividus and S. purpuratus.

Figure S5 .
Figure S5.Expression of recently evolved sea urchin genes, related to Figure 4. (A) Inferred WGCNA modules with their normalised fold-change enrichment in gene expression datasets (with on the rightmost side, the description of the module expression profile), marplot indicate the phylostrata of genes assigned to each module according to their origin (middle) or and last duplication nodes (right).Expression of GPCR135 (B) and Kirrel (C), two expanded gene families in sea urchins.Few duplicates usually retain a strong and ubiquitous expression while other duplicates acquired a specialised expression profile.Stage and gene are arranged according to clustering of euclidean distances and 'complete' method in pheatmap.
Figure S7.Transcription factor activity and gene regulatory networks during P. lividus development, related to Figure 7. (A) Expression of selected transcription factors with an enriched footprint in the same order as in Figure 7D.(B) Pairwise comparisons of footprint enrichment scores displayed as 'volcano' plots.Genes belonging to top centiles of fold-change are displayed in red.(C) Distribution of Network vertices degree (centrality) for non-TFs and TFs showing the high level of connectivity for some TFs.(D) Average centrality (z-score normalised) for various TF classes across developmental stages highlighting for instance the importance of Fox during mid-development (Blastula and Gastrula) and that of Sox during later development (Prism).(E) Skeletogenic GRN captured using Network reconstruction at the early blastula stage.

TABLE
d RESOURCE AVAILABILITY B Lead contact B Materials availability B Data and code availability d METHOD DETAILS B DNA extraction and sequencing B Genome assembly B Annotation B Long non-coding RNAs B Synteny and gene family reconstruction B Gene expression analyses B ATAC-seq B Cut&Tag B In situ hybridization B Cell Genomics 3, 100295, April 12, 2023 e2