Hagfish genome elucidates vertebrate whole-genome duplication events and their evolutionary consequences

Polyploidy or whole-genome duplication (WGD) is a dramatic genomic event commonly invoked causally in organismal evolution¹. The generally accepted ‘2R hypothesis’^2,3 suggests that two rounds of WGD occurred during early vertebrate evolution (referred to as 1R and 2R); however, their timing and macroevolutionary consequences remain unclear^4–6. Most studies agree that 1R occurred before the divergence of living vertebrates, but debate centres on whether 2R predated^7,8 or postdated^9–12 the divergence between cyclostomes and gnathostomes (Fig. 1c). Reconstruction of the ancestral vertebrate karyotype is fundamental to unravel the timing of 2R^8,12–15, but this goal has been stymied by a dearth of cyclostome genomes. The recently described genome of the sea lamprey (Petromyzon marinus) has been interpreted to support 2R occurring before⁸ or after¹² the gnathostome–cyclostome split, or not at all (with the karyotype diversity explained as the result of large-scale segmental duplications^16,17). Analysis of the Arctic lamprey (Lethenteron camtschaticum) genome has suggested that 2R occurred in the gnathostome lineage while independent WGD event(s) might have occurred in the lamprey lineage^11,18, perhaps shared with the hagfish^11,19 (Fig. 1c). However, the lack of a hagfish genome assembly, the only major vertebrate group without a reference genome, has challenged attempts to constrain the number and phylogenetic timing of ploidy events in early vertebrate evolution. Here we describe the outcome of sequencing and comparative analysis of the genome of the inshore hagfish, Eptatretus burgeri (Fig. 1a,b).

Similar to the lamprey²⁰, the hagfish genome undergoes somatic programmed DNA rearrangement in the way of chromosome elimination²¹, making it crucial to obtain a reference assembly from a germline source. We sequenced DNA extracted from the testis of a single sexually mature male of E. burgeri and generated a preliminary draft assembly using ~240X of short-read Illumina data assisted by a Chicago in vitro proximity ligation assay at Dovetail Genomics²² (Supplementary Table 5). We estimated the genome of E. burgeri at 3.12 Gb on the basis of k-mer frequency distribution (Extended Data Fig. 1a and Supplementary Table 4) in line with other hagfish species (~2.2–4.5 Gb)²³. Chromosome conformation capture (Hi-C) data obtained from the testis DNA of a second individual were used to further scaffold the genome into a final assembly (v.4.0) containing 19 contact clusters, which we consider as chromosomes for subsequent analyses (Fig. 1d), and 9,295 unplaced scaffolds and contigs (Methods and Supplementary Tables 6–10). The genome was annotated following the Ensembl annotation pipeline²⁴, assisted by RNA-seq from 9 different adult tissues, and previous embryonic and juvenile transcriptomics data¹⁹ (Supplementary Table 11 and Methods). We generated a final gene dataset of 16,513 protein-coding genes (with 27,960 transcripts), 446 long intergenic non-coding (linc)RNAs and a minority of other classes of non-coding RNA genes (Extended Data Fig. 1b). A total of 180 microRNA (miRNA) genes were found in the E. burgeri genome conserved with the hagfish Myxine glutinosa²⁵ belonging to 77 miRNA families and catalogued at MirGeneDB.org²⁶. The germline haploid number of E. burgeri is 26. However, chromosome elimination occurs in somatic tissues of the hagfish, by which, in the case of E. burgeri, 8 pairs of microchromosomes are eliminated during development²¹ (somatic n = 18). Cluster 19 and unplaced contigs/scaffolds probably correspond to these difficult-to-assemble microchromosomes, which presumably consist mainly of highly repetitive sequences and contain almost no protein-coding genes²¹. Consistently, 98.3% (16,240/16,513) of annotated genes are located in clusters 1–18.

BUSCO analyses show high levels of completeness of the hagfish genome (96.0 and 94.2% of single orthologues are present in the assembly and annotation, respectively; Fig. 1e and Extended Data Fig. 1c). GC-content distribution pattern analysis of the hagfish and other chordate genomes shows that the E. burgeri genome represents an intermediate condition between the lamprey and other chordates (Extended Data Fig. 1d), although having an overall GC content similar to that of the lamprey (46.7% and 48.1% for the hagfish and lamprey, respectively). While lamprey protein-coding gene sequences have been demonstrated to pose difficult challenges for comparative analyses due to their high GC content²⁷ (64.0%), the lower content in hagfish coding sequences (50.4%) is within the typical range of most gnathostomes and non-vertebrate chordates (42.5%–53.4%; Extended Data Fig. 1e,f and Supplementary Table 12). Lamprey represents an outlier in terms of both codon usage bias and amino acid composition, while the hagfish is more similar to other vertebrates (Fig. 1f,g). The hagfish genome contains, on average, significantly longer introns and intergenic regions than other vertebrates (P < 2.2 × 10⁻¹⁶, two-sided Wilcoxon rank-sum test), while the average length of coding sequences is similar to that of other chordates (Extended Data Fig. 1g, Supplementary Fig. 4 and Supplementary Tables 18–22). This might explain why hagfish genomes are larger than lamprey genomes²³. Altogether, the hagfish genome provides essential, complementary information to lamprey genomes, especially in analyses such as gene tree reconstruction and comparative genomics.

Whether hagfish form a clade with lampreys (Cyclostomata) or represent the sister group to all other vertebrates (including lampreys) has depended on whether molecular or morphological evidence are considered (Fig. 1c). Morphological studies historically supported cyclostome paraphyly but more recent analyses have recovered cyclostome monophyly (reviewed in ref. ²⁸). Phylogenies inferred from molecular evidence have almost exclusively recovered cyclostome monophyly (reviewed in ref. ²⁵). We used Bayesian inference to reconstruct the phylogeny of vertebrates from an alignment of 190 single-copy genes in all taxa analysed (84,017 sites), strongly supporting a monophyletic Cyclostomata (Fig. 2 and Extended Data Fig. 2a). We calculated the likelihood of gene duplication and loss patterns under the competing phylogenetic hypotheses²⁹ (see Methods), finding that patterns of gene gains and losses better fit cyclostome monophyly. To further compare the two alternative hypotheses of hagfish relationships, an approximately unbiased (AU) test³⁰ was performed, which strongly rejected cyclostome paraphyly (log likelihood difference = 7,947.7, AU = 0.004, multiscale bootstrap probability < 0.001). These results corroborate previous molecular analyses and recent morphological studies²⁸, supporting the view that cyclostomes are monophyletic.

To better understand the genomic changes accompanying major transitions in chordate evolution, we used a phylogeny-aware comparative genomic approach^31–33 to infer ancestral gene complements and gene family gains and losses across the vertebrate tree (Fig. 2, Extended Data Fig. 3 and Supplementary Tables 23–26). We observed two peaks of gene novelty in both the vertebrate and gnathostome stem-lineages (novel genes: +560 and +771, respectively), also characterized by a very low amount of gene losses (−341 and −382, respectively) when compared with other deuterostome and chordate nodes (Fig. 2 and Extended Data Fig. 3a). Furthermore, the fraction of highly retained novel gene families (also known as novel core genes, that is, genes that are not lost in descendant lineages and by convention indicated by ++) is the highest in the last common ancestors of vertebrates, gnathostomes and cyclostomes (novel core genes: ++81, ++86 and ++98, respectively; Fig. 2 and Extended Data Fig. 3a). These are notably larger than those observed in other major evolutionary episodes in metazoan evolution^31,33, but generally similar to a recent study using more chondrichthyan but only two invertebrate genomes³⁴, suggesting that the emergence of new gene families played important roles in the origin and diversification of early vertebrates. Gene Ontology (GO) enrichment analyses demonstrate that the origin of vertebrates was characterized by the appearance of genes involved in signalling pathways, cell communication and transcriptional regulation (Supplementary Tables 24 and 25), while novel core genes involved in immunity played an important role in the origin of gnathostomes (Supplementary Tables 24 and 26). Consistently, gnathostomes and cyclostomes convergently evolved independent adaptive immune systems, based on immunoglobulins in the former, and in variable lymphocyte receptors in cyclostomes³⁵ (Supplementary Fig. 6 and Supplementary Table 27). The largest fraction of gene losses occurred in the ancestral cyclostome lineage (Fig. 2), suggesting that a strong asymmetric reduction of gene complements accompanied the early evolution of the group. For instance, the hagfish genome lacks several vision and circadian rhythm-related genes, probably associated with its vestigial eyes (Supplementary Fig. 7 and Supplementary Table 28). Inferred rates of gene duplication (irrespective of the duplication mechanism) across Metazoa identify widespread duplications associated with the vertebrate and teleost stem-lineages (Extended Data Fig. 2b,c), probably reflecting the 1R, 2R and teleost 3R WGD events³⁶. We also inferred high duplication rates in each of the lineages leading towards crown-gnathostomes and lampreys (Extended Data Fig. 2b,c) which might suggest large-scale duplications associated with these groups, consistent with the WGD events proposed recently^12,18. This type of analysis, however, cannot discriminate between WGD and other large-scale gene duplication mechanisms.

The number of Hox clusters and ancestral WGD events are usually correlated; hence, the former has been used as a genomic marker of the latter. The presence of 6 Hox clusters in lamprey genomes^17,18 has been interpreted to indicate the possibility that more than two WGD events occurred in this lineage¹⁸. We have extended previous observations in E. burgeri¹⁹ to confirm the presence of 40 Hox genes arranged in 6 complete Hox clusters (Fig. 3). Two of the hagfish clusters are located in the same chromosome (cluster 3), separated by >80 Mb, probably the result of chromosomal shuffling due to the intense reorganization of the hagfish genome from ancestral chromosomes (ACs; see below). Phylogenetic analyses of Hox coding sequences have long proven inconclusive to determine the orthology relationship between lamprey and hagfish Hox counterparts¹⁹. We thus applied a microsynteny conservation approach using extended Hox loci which, together with phylogenetic analyses of selected non-Hox syntenic genes, allowed us to establish clear one-to-one orthologous correspondences between hagfish and lamprey Hox clusters, named α to ζ after the lamprey clusters¹⁸ (Fig. 3, Extended Data Fig. 4a–d, Supplementary Fig. 8 and Supplementary Tables 29 and 30). This suggests that the crown-cyclostome already possessed 6 Hox clusters, distinct from the ancestral crown-gnathostome, which possessed 4 clusters (Supplementary Fig. 9). This observation provides further evidence of cyclostome monophyly, by suggesting that lampreys and hagfish share a genome history exclusive of gnathostomes. This implies that the events suggested from the different analyses of the Arctic lamprey genome, two extra WGDs¹⁸ or a triplication¹¹, might have occurred in early cyclostome evolution, probably before the lamprey and hagfish divergence^11,19.

The reconstruction of the pre-WGD vertebrate proto-karyotype by means of macrosynteny analysis stands as the most robust approach to test the 2R event and its phylogenetic position³⁷. Earlier attempts at reconstructing the ancestral vertebrate karyotype have yielded widely disparate outcomes, indicating 10–13 (refs. ^14,16) or 17–18 (refs. ^8,11–13,38) ancestral pre-duplicative chromosomes. These reconstructions have also unveiled a perplexing scenario where lampreys’ divergence from gnathostomes occurred either before or following the 2R event^8,11,12,18. To shed new light on early vertebrate genome evolution, we performed a macrosynteny conservation analysis between gnathostomes, cyclostomes and selected invertebrate deuterostomes. First, to minimize noise from lineage-specific fusion and fission events, we reconstructed ancestral chicken³⁹ and spotted gar⁴⁰ genomes using elephant shark⁴¹ as an outgroup, obtaining an almost perfect one-to-one chromosome orthology (Supplementary Table 31). Next, to infer the ancestral vertebrate karyotype, we elaborated a map of homology relationships between the genes of these slow-evolving gnathostome genomes with the chromosome-level genome assembly of the sea cucumber Apostichopus japonicus (echinoderm) as a pre-duplicative outgroup species⁴². With this, we inferred a proto-vertebrate karyotype of 17 ACs (Supplementary Fig. 15 and Supplementary Table 32).

We mapped genes from the Belcher’s lancelet (Branchiostoma belcheri) genome⁴³ to each AC (Methods and Supplementary Fig. 21) using very stringent criteria, requiring homology relationships of an amphioxus gene with both a sea cucumber gene and several chicken and/or spotted gar genes, and all anchored to the same AC. In total, we mapped 5,065 Belcher’s lancelet genes to AC1–17 (ranging from 115 to 534 genes in AC1–AC16; AC17 consisted of only 20 genes and was thus excluded from several subsequent analyses; Supplementary Table 34). With these in hand, we corroborated our ancestral vertebrate karyotype reconstruction through comparisons with the chromosome-scale genome of the amphioxus Branchiostoma floridae¹² (Supplementary Tables 35 and 36).

Our inference of an ancestral karyotype with 17 ACs matches a previous study⁸, has minor differences with other 17-chromosome inferences^12,13 and depicts one less chromosome than more recent studies^11,38 (Supplementary Tables 37 and 38). In our model, 4 of the 17 ACs (AC1, AC2, AC3 and AC6) each correspond to 2 or 3 linkage groups (putative chromosomes) of the sea cucumber genome (Supplementary Fig. 15 and Supplementary Table 32) as well as to distinct homologous chromosomes in B. floridae (Supplementary Table 36), suggesting that these 4 ACs probably originated via fusions of ancestral chordate chromosomes in the vertebrate lineage before 1R³⁸. The difference from previous 18-chromosome models is that while we consider that the vertebrate AC3 is a single chromosome resulting from a pre-1R fusion event of two ancestral chordate chromosomes (Supplementary Fig. 23a), others^11,12,38 consider that these two chromosomes remained separate through 1R (Nakatani’s Pvc8 and 9, or Simakov’s CLGQ and CLGI, respectively). While Pvc8/CLGQ and Pvc9/CLGI are consistently co-located in gnathostome chromosomes, they remain separate in invertebrate karyotypes^11,12,38. We did not find any signals of linkage between Pvc8/CLGQ and Pvc9/CLGI in the lamprey and the hagfish genomes (Supplementary Figs. 25 and 26). Therefore, there exist two alternative scenarios: (1) the 18-chromosome model implies that two independent pairwise fusions occurred after 1R in a stem gnathostome, mimicking a single pre-1R fusion event (Supplementary Fig. 23a); and (2) our 17-chromosome model requires symmetric fissions of two AC3-derived post-1R chromosomes occurring in an ancestral cyclostome (Supplementary Fig. 23b). Although in silico simulations show that a scenario of pairwise post-1R fusions would not be extremely rare (30% of cases expected by chance; Supplementary Table 40 and Methods), we believe the pairwise fissions to be more plausible given the higher level of reorganization found in cyclostome karyotypes (see next section). Altogether, while we propose a scenario involving 17 ancestral vertebrate chromosomes, a scenario with 18 chromosomes^11,38 is also possible.

Importantly, all ACs correspond to sets of four (11/17) or three (6/17) paralogous chromosomes in the gnathostomes chicken and gar, a strong genome-wide pattern consistent with 2R^2,11–13. We stringently selected 701 sets of orthologous genes (Methods and Supplementary Table 39) between the sea cucumber, chicken, spotted gar and an AC gene (from B. belcheri), and built robust chromosome-level phylogenies with a median of 38 concatenated gene sets across each of the ACs (Extended Data Fig. 5a and Supplementary Table 39). The highly supported, clear-cut topologies further support the existence of 2R in gnathostomes and depict the exact evolutionary trajectory from each AC to their modern chicken and spotted gar descendants (Extended Data Fig. 5a,b). Our reconstruction of gnathostome karyotype evolution involved 8 fusion events that took place after 1R but before 2R (Extended Data Fig. 5b), similarly to what has been previously found by others^8,11,12 (Supplementary Table 41). Furthermore, we found a significant gene retention asymmetry after 2R, with a median of 1:2.28 genes per ohnologous (duplicates that originate through WGD, after ref. ³) chromosome pair, but not after 1R (median 1.16; P = 3.4 × 10⁻⁷, Wilcoxon rank-sum test; Extended Data Fig. 5b). This pattern is consistent with previous studies suggesting that 1R was an autotetraploidization event and 2R an allotetraploidization event^11,12 (but see ref. ⁴⁴ on asymmetric gene retention after teleost 3R autotetraploidy).

We next tested hypotheses of WGD timing relative to cyclostome divergence. We assessed the phylogenetic signal of hagfish and lamprey genes anchored to 661 orthologous gene sets (Supplementary Table 42), including elephant shark orthologues as a control for the 2R signal and amphioxus genes as outgroups. Approximately 73.2%, 79.1% and 75.7% of trees including hagfish, lamprey or both hagfish and lamprey orthologous genes, respectively, are compatible with shared 1R (Fig. 4a). However, while 99.5% of elephant shark gene tree topologies are 2R-compatible, only 19.1% of hagfish, 10.6% of lamprey and 8.2% of cyclostome (including both lamprey and hagfish) gene trees are compatible with a 2R history (Fig. 4a and Supplementary Files 5–8). Thus, we find strong support only for 1R as shared among cyclostomes and gnathostomes.

To further confirm the timing of the 1R event, we investigated whether signals of the inferred four pre-1R and eight post-1R fusion events are present in cyclostomes. When assessing how hagfish and lamprey chromosomes descended from the 17 ACs, we found that the hagfish genome displays a large amount of rearrangement (at least 52 fusions detected), making any signal of hypothetically shared events unreliable (Extended Data Fig. 6b). However, most lamprey chromosomes are descendants of single ACs¹² (Extended Data Fig. 6a), making the lamprey a better model to investigate these rare genomic changes. We found that the sea lamprey genome^17,45 bears signals of three (AC1, AC2 and AC6) and the hagfish genome of two (AC1 and AC2) of the four pre-1R fusions (Supplementary Figs. 25 and 26). On the other hand, similar to previous studies^11,12 we did not find any reliable signal of the eight post-1R fusions detected in gnathostomes, suggesting that the lamprey and hagfish diverged after the 1R but before all eight post-1R/pre-2R fusions^11,12. Taken together, our comprehensive phylogenetic analysis and the constraints given by pre- and post-1R chromosomal fusions provide strong evidence in favour of a pan-vertebrate 1R event, but constrains 2R to the gnathostome lineage as recently suggested in similar analyses^11,12.

It has been suggested that the lamprey genome has been shaped by either three duplicative events¹⁸ or a hexaploidization¹¹. The presence of six orthologous Hox clusters in both the lamprey and the hagfish (Fig. 3) implies that this is the ancestral condition for cyclostomes and supports the triplication event^11,19. Although we find that multiple chromosomes and large chromosomal sections are descendant copies of each AC in both cyclostome groups, the extensive rearrangements observed in the hagfish and the large haploid number in the lamprey impede chromosome-level macrosynteny conservation analysis to distinguish intraspecific ohnologous and interspecific orthologous relationships. To confidently infer karyotype evolution in cyclostomes, we developed a new metric, the ‘overlapping ratio’ (OR), to measure the similarity of gene retention profiles of any two chromosomes hypothetically descending from a common AC (Fig. 4b and Supplementary Fig. 29; gene-poor AC17 was excluded from this analysis, which required at least 20 genes retained in each descendant chromosome). A retention profile is defined by a vector listing the presence or absence of genes on a modern vertebrate chromosome from their corresponding AC. Therefore, we expect the OR of chromosomes deriving from a duplication event to be significantly higher than that of chromosomes deriving from an ancestral fission followed by gene translocations. As proof of concept, we applied this metric to gnathostomes: knowing that their genomes have been shaped by the 2R event, we found that the median OR of ohnologous chromosome pairs in chicken or spotted gar was 0.49 (interquartile range, IQR: 0.44–0.56) and 0.54 (IQR: 0.47–0.65), respectively (Fig. 4c), while OR value for simulated fission-derived chromosome pairs was never larger than 0.15, indicating that ohnologous chromosomes indeed share more retained genes (Fig. 4c and Methods).

We then applied the OR metric to the sea lamprey (after correcting misassemblies, assisted by a meiotic map of the Pacific lamprey Entosphenus tridentatus¹⁷ and confirmed by the recent chromosome-level genome assembly⁴⁵) (Supplementary Table 50 and Supplementary Fig. 24) and the hagfish, defining ohnologous chromosome pairs as those with OR > 0.15. We found that the median OR between putative ohnologous chromosomes was 0.30 (IQR: 0.23–0.36) and 0.29 (IQR: 0.23–0.37) for the lamprey and hagfish, respectively (Fig. 4d, Extended Data Fig. 7a,b and Supplementary File 9). Using this approach, we found that most ACs analysed (12/16 or 75%) have descended into three or more mutually ohnologous chromosomes in both the lamprey and hagfish (Fig. 4f), suggesting that at least a second WGD might have occurred in cyclostomes. In both genomes, at least five chromosomal regions are direct descendants of each of the same five ACs (1, 2, 6, 10 and 14), with 3 ACs contributing to 6 chromosomes each in the lamprey (Fig. 4f). We do not find more than 6 descendant copies from any AC, supporting a whole-genome triplication in the cyclostome lineage as previously proposed in the analysis of the lamprey genome¹¹. The distribution of multiplicity across ACs is highly correlated across the two species (Spearman ρ = 0.91), suggesting that this triplication event is conserved between the lamprey and hagfish and thus occurred in an ancestral cyclostome. It is expected that OR will decrease with each WGD (Supplementary Fig. 29), hence the lower value in cyclostomes is consistent with the occurrence of this larger polyploidy event. While our data do not definitively rule out the possibility of two cyclostome-specific WGD events followed by extensive chromosome losses, this scenario is less plausible than a single triplication event, particularly given the absence of instances with eight copies of any chromosomal region in the lamprey or the hagfish.

To further confirm that this proposed triplication event is shared by hagfish and lamprey, we extended the use of the OR to detect putative orthologous chromosomes. During the process of diploidization after a polyploidy, two descendant chromosomes diverge and fix their mutations independently, hence it is expected that interspecific orthologous chromosomes will have more similar gene retention profiles than intraspecific ohnologous chromosomes (as long as rediploidization precedes speciation⁴⁶). Accordingly, orthologous chromosomes of chicken and spotted gar have a median OR = 0.96 (IQR 0.95–0.98; Fig. 4e and Supplementary Table 44) and clustering-based analysis based on gene retention profiles places chicken and gar orthologous chromosomes closer to each other, completely reflecting the phylogenetic signal (Extended Data Fig. 7c and Supplementary File 10). When we applied this approach to cyclostome genomes, we found the median OR = 0.84 (IQR: 0.74–0.91; Fig. 4e and Supplementary Table 45) for 52 (~87%) chromosome pairings between lamprey and hagfish that putatively represent 1:1 orthologues (higher than that of ohnologous chromosomes; Supplementary Table 45) and only 8 (~13%) one-to-two or two-to-one ambiguous relationships, probably due to secondary independent chromosome losses in either group. Clustering analysis of retention profiles recovers orthologous relationships between lamprey and hagfish (Fig. 4g and Supplementary File 11). Overall, intra- and interspecific gene retention profile analyses indicate that a triplication event took place in the cyclostome stem-lineage; we refer to this as CR, to avoid confusion with the gnathostome-specific 2R event.

To investigate the immediate consequences of the independent CR event on cyclostome genome evolution, we first asked whether retained duplicates (ohnologues) are especially associated with developmental functions in the hagfish as in their gnathostome counterparts^13,47. GO enrichment analysis shows that hagfish gene ohnologues are also significantly enriched for functions associated with developmental processes (Extended Data Fig. 8a,b). Gnathostomes have increased their regulatory complexity (higher number of regulatory regions per gene), particularly of developmental ohnologues⁴⁷. We identified accessible chromatin regions (ACRs) as putative non-coding regulatory elements in the hagfish genome with an assay for transposase-accessible chromatin coupled to sequencing (ATAC-seq), using a total of two embryos of E. burgeri, each at a different stage⁴⁸ (45 and 53; Supplementary Fig. 32). We found a significantly higher number of ACRs per gene than in the cephalochordate amphioxus, similar to what has been observed in gnathostomes⁴⁷ (Fig. 5a), particularly in distal regions from transcriptional start sites (Fig. 5b,c and Extended Data Fig. 8c). This pattern is especially evident in developmental genes (Fig. 5d and Extended Data Fig. 8d–f), implying that their higher retention after a cyclostome CR event is underlain by a more complex regulatory landscape of developmental genes, as in gnathostomes⁴⁷.

In gnathostomes, retained duplicates can evolve via expressional specialization (reduction of expression domains of one of the ohnologues)⁴⁷, probably coupled to neofunctionalization rather than subfunctionalization (differential erosion of enhancers)⁴⁹. Taking advantage of adult transcriptome data across nine organs (see Methods), we next analysed the putative fates of hagfish ohnologues after CR (Supplementary Fig. 31). Hagfish duplicates also tend to reduce their expressional domains: over 68% and 71% of gene families subfunctionalized or specialized in the hagfish and chicken, respectively (Fig. 5e). Hagfish ohnologues that have potentially restricted their expression domains (subfunctionalization or specialization) are associated with a larger amount of regulatory elements and a higher sequence evolutionary rate than those that have maintained the ancestral patterns (Extended Data Fig. 8g,h), similar to gnathostomes (Extended Data Fig. 8i)⁴⁷. Furthermore, the largest portion of ohnologues with strong specialization (one or two ancestral expression domains) are expressed in the brain (Fig. 5f), mirroring the pattern observed in gnathostomes⁴⁷ (Extended Data Fig. 8j). In summary, our results indicate that cyclostomes and gnathostomes followed parallel evolutionary pathways after their independent WGD events. Genes gained a larger regulatory complexity, mostly on distal regions and especially in duplicates with developmental functions, which tend to be retained more often. Furthermore, specialization is a common fate of ohnologues associated with faster sequence evolution and the acquisition of novel regulatory elements that drive their tissue-specific expression. Alternatively, the possibility that a decrease in the number of regulatory elements took place in the amphioxus cannot be confidently ruled out.

Hypotheses on the role of WGD events in the origin and elaboration of the vertebrate bodyplan range from deterministic to permissive⁴. There can be no doubt that many vertebrate and gnathostome novelties are contingent on gene paralogues that are the product of the 1R and 2R events, although whether WGD played a causal role remains unclear. We employed two tests of a causal relationship: (1) absolute timing of the WGD events and the clades with which they are causally associated and (2) contrast in morphological phenotypic diversity before and after the WGD events. Using a dataset of 177 genes and 33 fossil calibrations, we provide estimated times at which duplicated subgenomes diverged. We adopted a sequential Bayesian approach in which the posterior clade age estimates from our species timescale were used as prior on the speciation nodes in our concatenated gene tree; to achieve statistical consistency, this requires that the same molecular loci are not used in the two analyses⁵⁰. A parallel analysis in which the concatenated gene tree was calibrated using the original fossil calibrations yielded results that are less precise but otherwise not materially different. Interpretation of the results depends on the nature of the ploidy event; in the case of autopolyploidy, we estimate the minimum timing of rediploidization, when the two subgenomes derived from WGD stopped homologously recombining, a process that can be asynchronous and span dozens of millions of years^{44,46,51–54}. In the case of allopolyploidy, the age estimate represents the speciation event that isolated each of the 2 subgenomes that later came together to form an allopolyploid genome⁵⁵. Our relaxed molecular clock analyses estimate the 1R event to have occurred 535.3–524.8 Ma (early Cambrian), 14.3–29.2 Myr before the divergence of crown-vertebrates (510.5–506.1 Ma; middle Cambrian) (Fig. 2); the CR event is dated to 500–492 Ma (late Cambrian), 23.5–36.5 Myr before the divergence of crown-cyclostomes (468.5–463.5 Ma; Middle Ordovician) (Fig. 2); and the 2R event is dated to 498.4–485.2 Ma (late Cambrian–earliest Ordovician), 35.1–53.3 Myr before the divergence of crown-gnathostomes (450.1–445.1 Ma; Late Ordovician) (Fig. 2).

To characterize morphological disparity across WGD events, we compiled a phenotype matrix composed of 577 traits for 278 living and fossil chordates encompassing all aspects of morphology, which we subjected to pairwise distance analysis followed by ordination using non-metric multidimensional scaling (NMDS) (Fig. 6a). This multivariate approach groups organisms with similar suites of characteristics while separating organisms with dissimilar traits, providing us with a relative measure of bodyplan diversity. The results of this analysis show that each genome duplication is followed by an increase in morphological disparity through occupation of novel regions of morphospace (Fig. 6b), but the majority of chordate disparity (88–97% of the morphospace encompassed by a vertebrate convex hull) emerged subsequent to the 2R event (Fig. 6b). Thus, while 2R and CR are of comparable antiquity, there is a stark contrast in terms of bodyplan evolution and species diversity between the descendants of 2R and the other WGD events.

The sequencing of a hagfish genome has enabled us to better understand the evolution of early vertebrates. First, our analysis of vertebrate genomes, including the hagfish, now establishes a robust and accurate history of WGD events in early vertebrates, corroborating the idea that cyclostomes diverged from gnathostomes after the 1R but before the 2R event. This is consistent with early^9,10,13 and recent studies on the matter that included the lamprey^11,12. We think that debate over the timing of 2R⁵ can now be concluded. Our hagfish genome also confirms an additional genome-wide duplicative event in stem-cyclostomes, CR, which most probably was a triplication, as previously suggested¹¹. Thus, key vertebrate innovations (for example, elaborate tripartite brain, neural crest cell-derived tissues among other novelties⁵⁶) originated in a stem-vertebrate. However, at this point we cannot reliably establish whether these innovations pre- or postdate the 1R event. This basic vertebrate bodyplan was further elaborated independently in cyclostomes and gnathostomes as a result of their lineage-specific genome duplications, for instance, facilitating the evolution of different adaptive immune systems (immunoglobulin-based in jawed vertebrates, variable lymphocyte receptor-based in cyclostomes³⁵), or the appearance of key morphological innovations, such as the jaw and paired appendages in gnathostomes. Interestingly, these independent WGD events shaped their ancestral genomes in similar ways by permitting an increase in regulatory complexity, especially of genes with roles in development. Duplicates of developmental genes are indeed more likely to be retained in both lineages, highlighting the crucial role of development in evolution of novel complex traits.

The contrasting morphological phenotypic consequences of 2R versus the other WGD events might suggest that there is no direct causal relationship or that there should be no general expectation of macroevolutionary consequences from WGD events despite their clear impact in increasing the regulatory potential of the genome. Another possibility is that the 2R event is different in nature from the 1R and CR events. Indeed, a number of recent studies together with our own analyses (Extended Data Fig. 5b) have suggested that while 1R was probably an autopolyploidy event, 2R was an allopolyploidy^11,12,57. This is significant since it impacts our interpretation of the absolute timing of 2R, with the age estimate reflecting the divergence of the two lineages that later hybridized, not the allopolyploidy event itself⁵⁵. Thus, the event occurred later than 498.4–485.2 Ma (late Cambrian), potentially coinciding with or even postdating the acquisition of gnathostome novelties that accrued among ostracoderms^4,58 before the divergence of crown-gnathostomes (which definitively postdated 2R) at 450.1–445.1 Ma (Late Ordovician). By the same token, the macroevolutionary consequences of allopolyploidy are expected to be more immediate than those of autopolyploidy, resulting in chromosomal rearrangements, changes in chromatin structure, DNA methylation, gene expression and the activation of transposable elements^59–61, extensive and immediate changes that promote species, and ecological diversification^62,63 as well as evolutionary novelty^64–67. This may go some way to explain why the evolutionary consequences of the 2R WGD are so much greater, leading to the profound diversification of gnathostome bodyplans that have dominated vertebrate communities since the early Palaeozoic.

No statistical methods were used to predetermine sample size. The experiments were not randomized and investigators were not blinded to allocation during experiments and outcome assessment.

This research was supported by grants from the Ministry of Science and Innovation of Spain (PID2021-125078NA-I00) and a KAKENHI Grant-in-aid for Scientific Research C from the Japan Society for the Promotion of Science (19K06798) to J.P.-A.; the National Key R&D Program of China (2019YFA0802600), the Chinese Academy of Sciences (ZDBS-LY-SM005), the National Natural Science Foundation of China (31970565) and the Open Research Program of the Chinese Institute for Brain Research (Beijing) to Y.E.Z.; the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB13000000) to Y.E.Z and W. Wang; the Wellcome Trust (to F.J.M., Grant number 108749/Z/15/Z); the Leverhulme Trust (RF-2022-167 to P.C.J.D); the Biotechnology and Biological Sciences Research Council (BB/T012773/1 to P.C.J.D); the Natural Environment Research Council (NERC) grant (NE/P013678/1 to P.C.J.D. and D.P.), part of the Biosphere Evolution, Transitions and Resilience (BETR) programme, which is co-funded by the Natural Science Foundation of China (NSFC); the John Templeton Foundation (Grant 62220 to P.C.J.D. and D.P.; the opinions expressed in this publication are those of the author(s) and do not necessarily reflect the views of the John Templeton Foundation); Wellcome Trust Seed Award (210101/Z/18/Z to J.P.); the Tromsø Research Foundation (Tromsø Forskningsstiftelse, TFS) to B.F. (20_SG_BF ‘MIRevolution’); a National Science Foundation (NSF) grant (1755418) to M.H.; and National Institutes of Health (NIH) grants (R01AI072435 and R35GM122591) to M.D.C. Computing was jointly supported by the HPC Platform of BIG and that of the Scientific Information Centre of IOZ. We thank O. Kakitani and the members of the Fishery Association in Gotsu City, Shimane Prefecture, Japan, for assistance in collecting hagfish; the technical support staff of the Laboratory for Phyloinformatics in RIKEN Kobe for sequencing data production, and RIKEN aquarium staff; T. de Dios Fernández for animal illustrations. For the purpose of open access, the authors have applied for a CC BY public copyright licence to any author-accepted manuscript version arising from this submission.

D.Y., Y.R., J.S., Y.L., W. Wan, W. Wang, Y.E.Z. and J.P.-A sequenced, assembled and evaluated the hagfish genome. F.J.M. generated and coordinated the annotation process at Ensembl. M.M., C.C., M.P. and W.A. generated Ensembl Compara data. A.d.M., D.Y., A.J.S.B. and J.W.C. generated orthogroups. R.F., J.F.F., F.S., S.D’A., C.B., C.A., G.S., B.F., K.J.P., S.D., M.H., J.P.R. and M.D.C. contributed to gene family curation. A.J.S.B., J.W.C., D.P. and P.C.J.D. performed phylogenomics and dating analyses. J.P. performed analysis of gene family evolution. D.Y., J.S. and Y.E.Z performed macrosynteny analysis. I.M. and M.I. advised on macrosynteny and ohnologue evolution analyses and discussed the data. J.P.-A, I.S., F.S. and S.K. obtained hagfish embryonic and adult material. J.P.-A performed ATAC-seq and RNA-seq experiments and ensured sequencing project management with BGI and RIKEN. M.U. and J.S. carried out regulatory genome profiling and ohnologue fate evolution analyses. J.N.K., E.M.C., R.P.D., S.G., E.R., R.S.S. and P.C.J.D. obtained fossil and calibration data, and performed morphological disparity analysis. W. Wang, S.K., N.I. and J.P.-A. conceived the project. S.K., N.I., P.C.J.D., W. Wang, Y.E.Z and J.P.-A contributed to project design. J.P.-A. coordinated the project. Y.E.Z., P.C.J.D and J.P.-A. wrote the paper, with inputs from all authors. All authors revised and approved the paper.

The Eptatretus burgeri (inshore hagfish) v.4.0 genome is available in NCBI GenBank under accession number GCA_900186335.3↗. Raw genome sequencing data together with adult RNA-seq data have been deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB21290↗. ATAC-seq data have been deposited in Gene Expression Omnibus (GEO) under accession number GSE247552↗. Supplementary files are available at FigShare (https://figshare.com/projects/Hagfish_Genome_Project/163186↗). Gene annotation used in this study is available at https://www.ensembl.org/Eptatretus_burgeri↗. A mirror of the UCSC Genome Browser containing hagfish assembly and annotations is available at http://ucsc.crg.eu/↗.

M.D.C. is a cofounder and shareholder of NovAb, Inc., which produces lamprey antibodies for biomedical purposes. J.P.R. is a consultant for NovAb. The other authors declare no competing interests.

Fergal J. Martin, Email: fergal@ebi.ac.uk

Wen Wang, Email: wwang@mail.kiz.ac.cn.

Philip C. J. Donoghue, Email: Phil.Donoghue@bristol.ac.uk

Yong E. Zhang, Email: zhangyong@ioz.ac.cn

Juan Pascual-Anaya, Email: jpascualanaya@gmail.com.