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

Abstract

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 specificities, we generated a chromosome-scale genome assembly for the sea urchin Paracentrotus lividus and an extensive gene expression and epigenetic profiles of its embryonic development. We found that, unlike vertebrates, sea urchins retained ancestral chromosomal linkages but underwent very fast intrachromosomal gene order mixing. We identified 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 identified 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.

Graphical abstract

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.

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 3

Intrachromosomal gene order rearrangement in sea urchins and vertebrates (A) The relationship between divergence time and gene order collinearity. Hsap, Homo sapiens; Ggal, Gallus gallus (chicken); Locu, Lepisosteus oculatus (spotted gar); Pliv, P. lividus; Spur, S. purpuratus; Lvar, L. variegatus. (B) Oxford plot between human and mouse, showing interchromosomal rearrangement but long colinear segments between the two species. (C) Oxford plot between the two sea urchin species, showing similar chromosomal architecture but reshuffled gene orders within chromosomes.

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

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 μm. (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

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 7

cis-Regulatory landscape conservation during sea urchin development (A) Sequence conservation scores in OCRs showing stage-specific activation (distinct from the TSS and consecutively expressed). Arrowheads indicate the stages’ highest non-coding conservation. Const., constitutive. Dev1 and Dev2 correspond to populations of OCRs that show broad activation domains. (B) Minimal Jensen-Shannon distance between staged transcriptomes of P. lividus, S. purpuratus, and B. lanceolatum. For each P. lividus stage, stages with minimal distance are highlighted by an arrowhead. (C) Endomesoderm GRN recovered by regulatory network analysis at the late blastula stage. (D) Enrichment scores in ATAC-seq footprints across the developmental stages of P. lividus. TFs playing a putative role in zygotic genome activation are highlighted (MZT TFs), with arrowheads pointing to Ets4 and SoxB1. (E) Regulatory activity around the He2 gene with RNA-seq (red) and ATAC-seq signals (blue), where regulatory elements located at the 5′ end of the gene include several TFBS footprints for MZT TFs, including Ets4 and SoxB1.