An atlas of fish genome evolution reveals delayed rediploidization following the teleost whole-genome duplication

Construction of TGD-aware teleost gene trees

We collected a data set of 101 vertebrate genomes, including 74 teleost fish, two nonteleost fish (bowfin and spotted gar, which did not undergo the teleost genome duplication or TGD), 20 other vertebrates of which six are mammals, and five nonvertebrate genomes (Supplemental Fig. S1; Supplemental Table S1). We used TreeBeST to reconstruct the phylogenetic relationships of 26,692 gene families across those 101 genomes (Methods; Vilella et al. 2009; Herrero et al. 2016). We then applied SCORPiOs to correctly place the TGD in these gene trees, using the bowfin and the spotted gar as reference outgroups (Parey et al. 2020). Briefly, SCORPiOs leverages synteny information to distinguish WGD-descended orthologs from paralogs when sequence information is inconclusive. After a WGD, orthologous genes are expected to be embedded in orthologous neighborhoods. For each individual gene tree, SCORPiOs “crowd-sources” additional information from local syntenic genes to identify errors in orthology relationships. SCORPiOs then reorganizes those gene trees to accurately position the WGD duplication node, if the synteny-consistent solution is equally supported by the sequence alignment. These 26,692 WGD-aware teleost gene trees predict that the ancestral genome of teleost fish contained 46,206 genes after the duplication event, in line with the latest estimates from the Ensembl database (49,255 ancestral teleost genes in release v100; Methods).

A high-resolution atlas of the TGD duplication across 74 teleost genomes

Teleost fish genomes are mosaics of duplicated regions, formed through rearrangements of duplicated ancestral chromosomes (Supplemental Fig. S2A). Long-standing efforts have been made to reconstruct the ancestral teleost karyotype before the whole-genome duplication (Jaillon et al. 2004; Kasahara et al. 2007; Nakatani and McLysaght 2017). According to the state-of-the-art reference, this ancestral teleost karyotype comprised 13 chromosomes (Nakatani and McLysaght 2017). Nakatani and McLysaght delineated between 353 and 690 megabase-scale genomic regions that descend from these 13 ancestral chromosomes in zebrafish, tetraodon, stickleback, and medaka.

We combined this pre-TGD ancestral karyotype with our gene trees to identify regions and genes that descend from sister duplicated chromosomes across all 74 teleosts in our data set (Methods; Supplemental Fig. S2). First, we transformed the reference segmentations of the zebrafish, tetraodon, stickleback, and medaka genomes (reference species) from 13 to 26 ancestral chromosomes after the duplication (1a, 1b, …, 13a, 13b). In each genome, we iteratively grouped regions from an ancestral chromosome into two postduplication copies by minimizing intragroup paralogs (Methods; Supplemental Fig. S2B). To assess robustness, we performed 100 groupings with random restart and found that genes were assigned to the same chromosome copy in 80% of iterations on average. We then identified orthologous ancestral chromosomes across all four species using gene orthologs, and arbitrarily named one of each pair “a” and “b” (Supplemental Fig. S2C).

Next, we propagated these ancestral chromosome annotations to the other 70 teleost genomes using gene orthology relationships (Supplemental Fig. S2D). For 1303 gene trees (7%), orthologs from the four reference species were not assigned to the same ancestral chromosome, and we resolved these inconsistencies by assigning all orthologous genes, including those of the reference species, to the most represented ancestral chromosome. This process results in a 74-species comparative genomic atlas with genes annotated to postduplication chromosomes, along with fully resolved orthology and paralogy links between all included species (Supplemental Fig. S2E).

This comparative atlas integrates 70%–90% of each genome into 24,938 postduplication gene families (Fig. 1A,B), and greatly improves upon the state-of-the-art both in terms of species and genomic coverage. The atlas reveals that teleost fishes vary substantially in their retention of duplicated gene copies (ohnologs) since the TGD, which make up 33% of the arowana genome but only 19% of the cod genome (Supplemental Table S2). In general, Osteoglossiformes, Otomorpha, and Salmoniformes species have retained more ohnologs from the TGD than other Euteleosteomorpha clades that diverged later. The atlas is available on the Genomicus database web server (see “Data availability and implementation”).

Quality checking the teleost genome comparative atlas

As a quality check, we assessed discrepancies between the pre-TGD karyotype reconstruction by Nakatani and McLysaght (2017) and our ancestral chromosome assignations. Both should be globally congruent because the ancestral karyotype provides the groundwork for the comparative atlas, but discrepancies can arise when orthologous genes are assigned to different preduplication chromosomes between the four reference species. In this case, we resolve the inconsistency by reassigning all orthologs to the most represented preduplication chromosome, which will be different from the original assignation in Nakatani and McLysaght (2017) for at least one species. We found that 9% of zebrafish genes differ in ancestral chromosome assignation between our comparative map and the pre-TGD karyotype, versus 2%–3% in medaka, tetraodon, or stickleback. This likely reflects small-scale rearrangements in zebrafish captured by our individual gene trees but missed by the macrosyntenic approach of Nakatani and McLysaght (2017). Alternatively, zebrafish genes may be more frequently placed incorrectly in our gene trees and assigned to the wrong orthology group, which may lead to erroneous ancestral chromosome reassignations. To explore this possibility, we identified gene trees that remain synteny-inconsistent after correction by SCORPiOs (i.e., whose orthology relationships are inconsistent with those of their surrounding genes), and therefore potentially contain orthology errors. Zebrafish genes reassigned to a different ancestral chromosome in our comparative map are not overrepresented in synteny-inconsistent trees (7% vs. 14% for all zebrafish genes), suggesting that their orthology relationships and chromosomal reassignations are overall well-supported by their sequences and syntenic gene neighborhoods.

Additionally, we used random noise simulations to show that the comparative atlas is robust to potential uncertainty or errors in the original ancestral genome reconstruction. The delineation of genomic regions descended from each preduplication chromosome from Nakatani and McLysaght (2017) relies on the identification of conserved synteny blocks and their breakpoints between teleost genomes and outgroups. Because breakpoint locations are challenging to determine—as also attested by lower posterior probabilities close to breakpoints in Nakatani and McLysaght (2017)—the region boundaries can vary in precision. We shifted the positions of these boundaries with increasingly large errors, mimicking situations in which up to 25% of genes change preduplication chromosome assignations in each of the reference genomes (Methods). We found that even large errors in region boundaries did not majorly affect the final atlas, with only 11% of genes changing ancestral chromosome assignations at the highest noise settings (Supplemental Fig. S3).

Finally, we report that correcting gene trees with SCORPiOs had a decisive impact on the establishment of the comparative atlas, enabling the inclusion of a significantly larger fraction of teleost genes (83% vs. 61%; Supplemental Fig. S4). The teleost comparative genomic atlas represents therefore a reliable, comprehensive, and robust resource for fish genomics.

The teleost duplication was followed by delayed rediploidization

The comparative genomic atlas is the first resource that allows an in-depth, genome-wide analysis of the TGD, and in particular of the early genome evolution processes that followed after the TGD. Previous work has revealed that auto- and allotetraploids significantly differ in their early evolution (Stebbins 1947; Garsmeur et al. 2014; Mason and Wendel 2020): specifically, autotetraploids initially harbor four near-identical chromosome sets, allowing multivalent meiotic recombination and tetrasomic inheritance to occur until preferential bivalent pairings evolve (i.e., cytological rediploidization). Previous studies in salmonids and Acipenseriformes have revealed that this process is highly dynamic, occurring in distinct temporal waves for different genomic regions, separated by as much as tens of millions of years (Robertson et al. 2017; Gundappa et al. 2022; Redmond et al. 2022). Sequence exchanges caused by prolonged multivalent recombination are detectable because they delay the divergence of duplicated regions until after cytological rediploidization occurs, resulting in different phylogenetic expectations regarding ohnolog sequence evolution. Depending on the respective timings of rediploidization and speciation, ohnologs can either follow the ancestral ohnolog resolution (AORe) or lineage-specific ohnolog resolution (LORe) models (Fig. 2A), introduced in Robertson et al. (2017). In the AORe model, cytological rediploidization occurs before lineages split, thus initiating ohnologous sequence divergence before speciation. In the LORe model, because cytological rediploidization has not been resolved before speciation occurs, ohnologs share more sequence similarity within clades than across clades, and can therefore be misidentified as clade-specific duplications.

To investigate the polyploidization mode and rediploidization processes of teleost fishes, we established whether meiosis was fully resolved with bivalent chromosomal pairing in the teleost ancestor before extant teleost lineages diverged, or whether some duplicated genomic regions still recombined during meiosis at that time. We focus our analysis at the Osteoglossiformes/Clupeocephala lineage split, which is the earliest possible speciation point in the teleost phylogeny (Parey et al. 2022), dated approximately at 267 Mya (Kumar et al. 2017). We selected a subset of six teleost genomes (Paramormyrops, arapaima, Asian arowana, zebrafish, medaka, and stickleback) to ensure an equal representation of Osteoglossiformes and Clupeocephala genomes in the data set, along with outgroup genomes (Methods; Supplemental Fig. S5). We developed an extension to SCORPiOs named LORelEi for “Lineage-specific Ohnolog Resolution Extension”. SCORPiOs LORelEi is built around two modules (“diagnostic” and “likelihood tests”) to identify delayed rediploidization in gene trees. The LORelEi diagnostic module leverages the gene tree correction applied by SCORPiOs to identify sequence-synteny conflicts. Sequence-synteny conflicts arise when the orthology relationships of a gene family are inconsistent with conserved synteny information, but correcting the WGD duplication node to rectify this inconsistency induces a significant drop in likelihood estimated from sequence divergence. We grouped gene trees with sequence-synteny conflicts according to their ancestral chromosome of origin in the comparative atlas. Three anciently duplicated chromosome pairs (3, 10, and 11) are significantly enriched in sequence-synteny conflicts, thus hinting toward prolonged recombination between these homeologs after the TGD (Fig. 2B; P < 0.05, hypergeometric tests corrected for multiple testing).

We next sought to confirm that the identified sequence-synteny conflicts were consistent with the phylogenetic expectations of delayed rediploidization. We used the LORelEi likelihood test module to explicitly compare the likelihoods for the tree topologies expected under the AORe and LORe rediploidization models, in which ohnolog resolution occurs before or after the Osteoglossiformes/Clupeocephala lineage split, respectively (Fig. 2A; Methods). We performed these likelihood tests on 5557 gene trees, which retain both ohnologs in at least one of the descending lineages. For 638 gene trees, the early resolution topology (AORe) was significantly more likely, whereas the late resolution topology (LORe) was favored for 1361 trees (likelihood AU-tests, alternative topology rejected at α = 0.05; no significant differences for the remaining 3558 trees). For example, the col12a1a - col12a1b ohnologs had stopped recombining and accumulated independent substitutions by the time Osteoglossiformes and Clupeocephala diverged (Fig. 2C). On the other hand, the map1aa - map1ab TGD ohnologs presumably still recombined when the taxa diverged, and both ohnologs diverged independently later on in each lineage (Fig. 2D). We mapped the location of these genes on the medaka karyotype, which remains close to the ancestral teleost karyotype, to identify chromosomal regions of ancestral and lineage-specific rediploidization (Fig. 2E; Methods). This visualization revealed that ancestral chromosome pairs 3, 10, and 11 had not yet rediploidized when teleosts started diversifying, with large runs of LORe ohnologous families spanning the entire length of their descendant chromosomes in medaka (Chromosomes 2, 3, 6, 21, and 23). The other homeologs appear as a mix of localized AORe and LORe regions, suggesting that rediploidization was ongoing.

This snapshot of the rediploidization status at the time of the Osteoglossiformes/Clupeocephala lineage split shows that rediploidization was not uniform across the ancestral teleost genome, consistent with results in salmonids (Robertson et al. 2017; Gundappa et al. 2022). In salmonids, different temporal waves of rediploidization have been linked with variations in gene functions (Gundappa et al. 2022). In teleosts, we found no evidence of functional enrichments in either the AORe genes, which rediploidized early, or the LORe genes carried by late-rediploidized ancestral Chromosomes 3, 10, and 11. Interspersed LORe genes on other homeologs were functionally enriched for different pathways, including “TGF-beta signaling pathway,” “regulation of transferase activity,” and “regulation of intracellular signal transduction” (corrected P-values < 0.05; Supplemental Tables S3, S4; Methods). One possible explanation to the functional differences between the two sets of LORe genes might be that homeologous Chromosomes 3, 10, and 11 remained tetraploid owing to specificities in their nuclear organization and topological features, which prevented rediploidization, whereas sequence exchanges at LORe loci on other homeologs may have been locally maintained as a result of selection, for example by removing deleterious mutations or transferring favorable alleles at all gene loci. In this speculative model, meiosis resolution would largely occur in waves because some chromosomes or chromosomal regions are mechanistically “easier” to rediploidize than others.

Altogether, we provide the first evidence that entire chromosomes experienced delayed rediploidization in teleosts and continued to exchange genetic material between homeologs for at least 60 million years after the teleost whole genome duplication. This prolonged exchange of genetic material between duplicated chromosomes after the TGD strongly suggests that the teleost ancestor was an autotetraploid.

The teleost genome experienced biased gene retention after duplication

Previous observations on a limited subset of genes have suggested that the ancestral teleost genome has experienced biased gene retention after the TGD (Conant 2020). Biased gene retention is generally considered a hallmark of allopolyploidization (Garsmeur et al. 2014), but has mostly been observed in plants, and the mechanisms of gene retention in polyploid vertebrate genomes may differ. We tested whether the ancestral tetraploid teleost underwent biased gene retention, with one duplicated chromosome copy systematically retaining more genes than the other. We computed gene retention in nonoverlapping 10-gene windows along duplicated chromosome copies in several teleost genomes, using an outgroup fish as proxy for the ancestral gene order (Fig. 3A,B; Supplemental Figs. S6, S7; Methods) and total gene retention on homeologs in the ancestor (Osteoglossocephalai) of all 74 teleosts (Supplemental Table S5). We uncover a consistent, significant bias in gene retention on ancestral chromosome pairs 3, 4, 7, 11, and 13, in which genes were preferentially retained on one homeolog over the other, regardless of study species. We find additional but weaker evidence for biased gene retention on ancestral chromosome pairs 2, 5, 6, 8, and 9 in some combinations of genome comparisons. We however do not find evidence of retention bias between ancestrally duplicated copies for chromosome pairs 1, 10, and 12. This preferential gene retention is not an artifact of the atlas construction because of genes unassigned to either duplicated chromosome: indeed, conservatively assigning such genes to the homeolog with the lowest retention did not offset the retention imbalance (Supplemental Table S6). In summary, we find that genes were preferentially retained on one duplicated chromosome in a subset of ancestral chromosomes, a character more frequently observed in allopolyploids (Garsmeur et al. 2014; Cheng et al. 2018). Of note, the chromosome copy with highest gene retention is typically annotated as “a” in the atlas as a consequence of the construction process. As the TGD is likely an autopolyploidization event, “a” and “b” correspond to different chromosome copies and are interchangeable for each chromosome pair—“a” and “b” chromosomes should not be interpreted as distinct parental subgenomes of allopolyploidization.

Differences in selective pressures and gene expression do not explain the observed bias in gene retention

Previous observations have suggested that biased gene retention in allopolyploid genomes reflects partial epigenetic silencing of one parental subgenome (Freeling et al. 2012; Bird et al. 2021). Genes on silenced chromosome copies are expressed at lower levels, and therefore subjected to lower selective pressures and ultimately to pseudogenization. We therefore investigated whether anciently duplicated chromosome copies that retained fewer genes have decayed because of unequal selective pressure. We estimated d_N/d_S ratios for 1263 pairs of TGD-derived paralogs conserved in at least 40 species, and tested whether ohnologs on one of the two anciently duplicated chromosomes systematically underwent more relaxed sequence evolution (Methods). On ancestral chromosome pairs 8 and 13, genes of one chromosome copy experienced a significant relaxation of selection compared with their ohnologous counterparts (Fig. 3C). Genes under relaxed selection were located on the homeolog with lower retention, as predicted if low-selective pressure promotes gene loss. However, homeologs with significant differences in selective pressure correspond only to two out of six chromosomes showing strongly biased gene retention (n = 13; P = 0.1923; Fisher's exact test), with nonsignificant differences in selective pressure in the same direction as the gene retention bias for two chromosome pairs (4 and 11) and in the opposite direction for the two others (3 and 7). Thus, differences in selective pressures do not explain the observed bias in gene retention.

In addition, we investigated whether ancestral chromosome copies show differences in gene expression, which would have been established after the TGD and have been maintained since. We find no bias in ohnolog gene expression between ancestral chromosomes from each pair. This result is consistent whether looking at average gene expression across 11 tissues in medaka (Fig. 3D), zebrafish (Supplemental Fig. S8), or tissue by tissue in either medaka or zebrafish (Supplemental Figs. S9, S10).

In conclusion, our findings are consistent with previous reports that gene retention on ancestral chromosomes was biased following the teleost duplication (Conant 2020). However, we find that biased gene retention in teleosts is not correlated with other features classically observed in allopolyploids, and our results suggest that biased gene retention can occur following autopolyploidization, possibly driven by distinct and underappreciated factors.

The comparative atlas enhances teleost gene and genome annotations

Finally, we investigated how the comparative atlas may improve crucial resources for fish evolutionary, ecological, and functional genomics. The Zebrafish Information Network (ZFIN) provides manually curated, high-quality reference annotations for zebrafish genes and implements rigorous conventions for gene naming (Ruzicka et al. 2019). These gene names are then propagated to orthologous genes in other teleost genomes, providing the basis of the entire teleost gene nomenclature and functional annotation transfers. In an effort to convey evolutionary meaning within the gene names, zebrafish paralogs descended from the TGD are identified with an “a” or “b” suffix. ZFIN guidelines recommend that adjacent genes should carry the same suffix when they belong to a continuous syntenic block inherited since the TGD, sometimes called syntelogs (Zhao et al. 2017). This aspiration has however been difficult to implement in the absence of a high-resolution map of zebrafish duplicated regions and their ancestral chromosomes of origin. To assess the consistency in consecutive zebrafish gene names, we extracted and compared zebrafish “a” and “b” gene suffixes along duplicated regions from the comparative atlas (Fig. 4). We find that despite previously mentioned efforts, zebrafish “a” and “b” gene suffixes are not consistent with the polyploid history of the zebrafish genome (Fig. 4A): 43% of gene suffixes would have to be reassigned in order to reflect the shared history of syntenic genes, which would be impractical to implement (Methods). Gene annotations are therefore unhelpful to study large-scale processes such as chromosome evolution or genome-wide rediploidization. Additionally, the ZFIN nomenclature does not impart a suffix to singleton genes, which correspond to TGD-duplicated genes in which one of the copies was eventually lost. As a result, only 26% of zebrafish genes are annotated with an “a” or “b” suffix in ZFIN.

Using the comparative atlas, we annotated which paralogs were retained for 84% of all zebrafish genes, providing a complementary resource that significantly extends the evolutionary annotation of gene histories in this species (Fig. 4B). The comparative atlas also includes similar annotations for all 74 studied teleosts, as shown for medaka, stickleback, and tetraodon (Supplemental Fig. S11). How ancient tetraploids return to a mostly diploid state is an active area of research, where distinguishing which paralog gene has been retained can be essential (Inoue et al. 2015; Robertson et al. 2017; Conant 2020; Simakov et al. 2020; Gundappa et al. 2022). The comparative atlas opens the possibility to formally identify and investigate which ancestral copies have been retained and lost during teleost diversification, and transfer functional gene annotations between model and nonmodel fish genomes. As such, the comparative atlas constitutes a biologically meaningful, historically accurate insight into reference gene annotations to support further investigations of teleost genome evolution.

Libraries and packages

Scripts to build the comparative atlas were written in Python 3.6.8 and assembled together in a pipeline with Snakemake version 5.13.0 (Köster and Rahmann 2012). The ete3 package (version 3.1.1) was used for phylogenetic gene tree manipulation and drawing. Other Python package dependencies used for plots and analyses include Matplotlib (version 3.1.1), seaborn (version 0.9.0), numpy (version 1.18.4), pandas (version 0.24.2), pingouin (version 0.3.4) for Wilcoxon paired tests and multiple testing correction, and SciPy (version 1.4.1) for Fisher's exact tests. Chromosome painting and synteny comparisons were drawn with the RIdeogram R package (Hao et al. 2020).

Genomic resources

Genome assemblies and annotations for the 101 vertebrate genomes were downloaded from various sources, including the NCBI, Ensembl version 95 and GigaDB. The source and assembly version used for each genome are listed in Supplemental Table S1. The gene coordinates files for the 74 teleost genomes in the comparative atlas have been deposited to Zenodo (https://doi.org/10.5281/zenodo.5776772).

Phylogenetic gene trees

Initial gene trees were built using the Ensembl Compara pipeline (Vilella et al. 2009). Briefly, starting from the sets of proteins derived from the longest transcripts in each genome, we performed an all-against-all BLASTP+ (Altschul et al. 1990), with the following parameters “-seg no -max_hsps 1 -use_sw_tback -evalue 1e-5”. We then performed clustering with hcluster_sg to define gene families, using parameters “-m 750 -w 0 -s 0.34 –O”. We built multiple alignments using M-Coffee (Wallace et al. 2006), with the command “t-coffee -type = PROTEIN -method mafftgins_msa, muscle_msa,kalign_msa,t_coffee_msa -mode = mcoffee”. Next, we conducted phylogenetic tree construction and reconciliation with TreeBeST, using default parameters (Vilella et al. 2009). Although TreeBeST remains the most efficient method to build gene trees for large data sets (Noutahi et al. 2016), it systematically infers a number of gene duplication nodes that are overly old and poorly supported (Hahn 2007). We therefore edited the TreeBeST gene trees: we used ProfileNJ (Noutahi et al. 2016) to correct nodes with a very low duplication confidence score (duplication score < 0.1, computed as the fraction of species that retained the genes in two copies after the duplication). Specifically, ProfileNJ rearranges subtrees below these poorly supported nodes to make them more parsimonious in terms of inferred duplications and losses. Finally, we ran SCORPiOs (version v1.1.0) to account for several whole-genome duplication events in the species phylogeny and correct gene trees accordingly: the teleost duplication (TGD), using bowfin and gar as outgroups. SCORPiOs reached convergence after five iterative rounds of correction. We identified 8144 teleost gene subtrees out of 17,493 (47%) that were inconsistent with synteny information, of which 5611 could be corrected (32% of all subtrees). We note that the corrected-to-inconsistent tree ratio (69%) is comparable to our previous application of SCORPiOs to fish data. Similarly, the proportion of errors in initial sequence-based gene trees is in line with our previous application of SCORPiOs to a data set of 47 teleost species (Parey et al. 2020). We also applied SCORPiOs to correct the nodes corresponding to the carp 4R WGD, using zebrafish as outgroup; and the salmonid 4R WGD, using Northern pike as outgroup. In the presence of LORe in salmonids, as is the case for teleosts, SCORPiOs will attempt to correct gene trees on the basis of synteny information but the correction will be rejected for being inconsistent with sequence evolution. As a result, the salmonid 4R WGD will be placed as lineage-specific duplications in LORe trees, consistent with sequence evolution.

Ancestral gene statistics

We calculated the predicted number of genes in the postduplication ancestral teleost genome using our set of 26,692 gene trees, and compared this to 60,447 state-of-the art gene trees stored in Ensembl Compara v100. Specifically, to calculate the number of genes inferred in the teleost ancestor (Osteoglossocephalai), we counted ancestral gene copies assigned to Osteoglossocephalai in the 26,692 and 60,447 TreeBeST phylogeny-reconciled gene trees. This ancestral gene number is an indirect but accurate approximation of the quality of inferred gene trees, because the major challenge is to accurately position duplications at this ancestral node.

Comparative genomic atlas

The FishComparativeAtlas pipeline annotates teleost genes to post-TGD duplicated chromosomes. It takes as input genomic regions annotated to preduplication chromosomes for four reference species (zebrafish, medaka, stickleback, and tetraodon), and the set of gene trees described above. Segmentation of the four teleost species with respect to the 13 ancestral chromosomes were extracted from Nakatani and McLysaght (2017) and genomic interval coordinates converted to lists of genes. All genomes were then reduced to ordered lists of genes. We first converted input genomic intervals from zebrafish genome assembly Zv9 to GRCz11 and from medaka genome assembly MEDAKA1 to ASM223467v1, by transferring boundary genes between assemblies using Ensembl gene ID histories. Next, we identified putative TGD sister regions within each of the four reference species, as regions sharing a high fraction of TGD-descended paralogs (Supplemental Fig. S2B). We grouped regions descended from the same preduplication ancestral chromosome, and iteratively annotated pairs of regions into internally consistent sets of “a” and “b” postduplication sister regions as follows: for each ancestral chromosome, we started from the largest descendant region and arbitrarily defined it as “a”. All regions sharing 50% ohnologs or more with this “a” region are identified as the “b” paralog region(s). Additional search rounds were then conducted to extend the “a” and “b” annotations in a stepwise manner to all regions descended from this ancestral chromosome. The required ohnolog fraction was decreased at each round, and the search was stopped when remaining blocks shared <5% ohnologs with previously annotated regions. Because this identification of duplicated regions was performed independently in each of the four species, “a” and “b” region annotations were homogenized to ensure consistency across species (Supplemental Fig. S2C). The homogenization step uses orthology relationships from the gene trees and stickleback as an arbitrary guide species: annotations were switched in other species when “a” segments shared more orthologous genes with the “b” region of stickleback. Where conflicts remained for individual gene annotations, we used a majority vote procedure across all four species to define the postduplication chromosome. Finally, annotations were propagated to all teleost genomes in the gene trees using orthologies.

Simulation of ancestral chromosome boundary shifts

To simulate uncertainty in interval boundaries in the original ancestral genome reconstruction from Nakatani and McLysaght (2017), we randomly drew new boundaries in the vicinity of their original location according to a Gaussian distribution with standard deviation σ varying in [5, 10, 25, 50, 75, 100] genes. These boundary shifts were independently generated for each of the four reference species, with n = 100 random noise simulations for each σ value and each reference species. In total, simulations generated 600 noisy input data sets that were processed with the FishComparativeAtlas pipeline to assess its robustness to noise. Each of the 600 produced outputs were then compared to the comparative atlas, by counting the proportion of gene families with a reassigned ancestral chromosome (Supplemental Fig. S3).

Early and late rediploidization gene tree topologies

Phylogenetic gene trees were built for the reduced set of 33 genomes as follows. We first filtered the CDS multiple alignments previously computed for the 74 teleosts and 27 nonteleost outgroup genomes set as described above, to retain only genes of all 27 outgroups and 6 teleost genomes, including three Osteoglossiformes: paramormyrops (Paramormyrops kingsleyae), arapaima (Arapaima gigas), and arowana (Scleropages formosus) and three Clupeocephala: zebrafish (Danio rerio), medaka (Oryzias latipes), and stickleback (Gasterosteus aculeatus) (Supplemental Fig. S5). We used these reduced multiple alignments to build phylogenetic gene trees with TreeBeST (Vilella et al. 2009), using default parameters and the option “-X 10”. This resulted in a set of 14,391 gene trees containing genes of the 33 retained genomes. We then ran SCORPiOs to correct trees for the teleost-specific duplication using the gar and bowfin genomes as outgroups.

To investigate the occurrence of delayed rediploidization, we implemented an extension to the SCORPiOs pipeline (SCORPiOs “LORelEi” for “Lineage-specific Ohnolog Resolution Extension”) to (1) identify gene trees characterized by sequence-synteny prediction conflicts and (2) perform likelihood AU-tests (Shimodaira and Hasegawa 2001) to evaluate how AORe and LORe rediploidization tree topologies are supported by gene sequence evolution. For (1), we identify gene trees that SCORPiOs attempts to correct based on syntenic information, but whose correction is rejected because of low sequence-based likelihood. For (2), we selected the 5557 gene families containing a gene in at least one reference outgroup (bowfin or gar) and resulting in distinct tree topologies under AORe and LORe. In practice, these topologies can be distinguished when at least one teleost group (Osteoglossiformes or Clupeocephala) retained both ohnologs, although not necessarily in the same species. For each of these 5557 families, we built three gene trees using RAxML 8.2.12 (Stamatakis 2014), with the GTRGAMMA model: the unconstrained maximum likelihood (ML) tree, the constrained AORe topology, and the constrained LORe topology (Fig. 2B). We then used CONSEL (Shimodaira and Hasegawa 2001) to test for significant differences in log-likelihood reported by RAxML (Stamatakis 2014). A tree topology was rejected when significantly less likely than the ML tree at α = 0.05.

We used Circos version 0.69–9 (Krzywinski et al. 2009) to visualize AORe and LORe ohnologs on the medaka genome. Before the Circos construction, we used the “bundlelinks” tool available in the circos-tools suite version 0.23 to bundle together consecutive genes with the same rediploidization mode, using 50 genes as the maximum distance parameter (-max_gap 50). Using this setting, isolated AORe and LORe are less visible (i.e., have thinner links) than high-confidence regions of consecutive genes displaying the same rediploidization mode. On the Circos, we annotate ancestral chromosomes corresponding to each medaka chromosome with color labels. For clarity purposes, we only add a label if more than 17.5% of genes of a given medaka chromosome are annotated to the ancestral chromosome.

Functional enrichment tests

For each of the AORe, LORe whole-chromosomes, and LORe interspersed sets, we extracted high-confidence ohnologs list in medaka. Specifically, we retained only ohnologs falling in high-confidence AORe and LORe regions defined by “bundled” gene families in the Circos representation (i.e., regions formed from neighboring genes with the same rediploidization mode, using 50-gene sized windows). Finally, we used the zebrafish orthologs of these genes (n = 248 zebrafish genes for AORe, n = 193 for LORe whole-chromosomes and n = 215 for LORe interspersed) to perform Gene Ontology and KEGG pathway enrichment analyses through the WebGestalt web server (Liao et al. 2019).

Gene retention on homeologous chromosomes

Because ancestral gene order is particularly difficult to reconstruct in teleosts because of an elevated rate of microsyntenic rearrangements and gene copy losses (Inoue et al. 2015; Nakatani and McLysaght 2017), we take advantage of a nonduplicated outgroup fish genome as a proxy for the ancestral gene order. Here, we make the assumption that consecutive genes on the outgroup genome, all assigned to the same preduplication chromosome, represent the ancestral gene order. Using the gene trees, we identify 10,629 1-to-1 orthologies between spotted gar genes and teleost preduplication gene families, and 11,599 with bowfin genes. We then used the gene order in the outgroup genomes as an approximation of the ancestral teleost gene order. We reduced outgroup genomes to these genes, and extracted all blocks of consecutive genes annotated to the same preduplication chromosome. We computed the percentage of gene copies retained on “a” and “b” homeologs in each extant duplicated genome, using nonoverlapping windows of 10 genes along these blocks.

Tests for relaxed selective pressures

We considered 1263 teleost gene families annotated in the comparative atlas for the d_N/d_S analysis, selecting all families that contained exactly two ohnologous copies in at least 40 teleost genomes (excluding salmonids and carps, which underwent additional WGDs), and exactly one orthologous copy in the bowfin and gar outgroups. We pruned the trees from any species with only one ohnolog copy or where additional gene duplications were present, to obtain “a” and “b” clades with the exact same species, for informative d_N/d_S comparisons. For each gene family, we used translatorX vLocal (Abascal et al. 2010) to (1) translate coding sequences, (2) align resulting amino acid sequences with MAFFT v7.310 (Katoh and Standley 2013) using option “‐‐auto”, (3) trim poorly aligned regions with Gblocks version 0.91b (Castresana 2000) using parameters “-b4 = 2 -b5 = h”, and (4) back-translate the sequences into codon alignments. We used the RELAX model in HyPhy (Wertheim et al. 2015) to estimate d_N/d_S ratios and test for significant relaxation or intensification of selection on branches of the “a” subtree compared with branches of the “b” subtree. Briefly, RELAX estimates d_N/d_S distributions across sites on the sets of “a” and “b” branches and fits a selection intensity parameter “k”, which captures the extent of selection intensification (k > 1) or relaxation (k < 1). Likelihood ratio tests are conducted to compare the alternative model with the k parameter to a null model without. We identified relaxation or intensification of selection on “a” versus “b” branches when the null model was rejected (P-values < 0.05, corrected for multiple testing using the Benjamini–Hochberg procedure). Finally, we performed Fisher's exact tests corrected for multiple testing using the Benjamini–Hochberg procedure, to identify chromosome pairs with significantly higher numbers of genes under relaxed evolution on one chromosomal copy.

Expression level of ohnologous genes

We used RNA-seq data sets across 11 tissues (bones, brain, embryo, gills, heart, intestine, kidney, liver, muscle, ovary, and testis) in zebrafish and medaka from the PhyloFish database (Pasquier et al. 2016), normalized into TPM (transcripts per million transcripts) and quantile normalized across tissues as previously described (Parey et al. 2020). Gene IDs were then converted from Ensembl version 89 to version 95 using conversion tables downloaded from BioMart (Smedley et al. 2009). Average expression across tissue (Fig. 3C; Supplemental Fig. S8) and by-tissue expression (Supplemental Figs. S9, S10) were calculated for ohnologs grouped by their ancestral chromosome of origin.

ZFIN gene names

Zebrafish ZFIN gene names were extracted using BioMart (Smedley et al. 2009) from the Ensembl database (version 95). We extracted the last letter of gene names, which represents “a” and “b” ohnology annotations in ZFIN. We then computed the minimal number of “a” and “b” ZFIN gene name reassignments that would be necessary to be consistent with the comparative atlas. In the comparative atlas, “a” and “b” labels are arbitrarily assigned to duplicated chromosomes; that is, genes descended from Chromosomes 1a and 1b could be swapped to 1b and 1a. To not artificially overestimate discordances, we first swapped such arbitrary “a” and “b” annotations to minimize differences with ZFIN. Finally, we counted the remaining number of “a” and “b” disagreements for zebrafish genes in the comparative atlas that were also annotated in ZFIN.

Software availability

All code for the FishComparativeAtlas pipeline is publicly available at GitHub (https://github.com/DyogenIBENS/FishComparativeAtlas). An archive containing a stable version of the code along with all input data (including gene trees) and the final atlas has been deposited in Zenodo (https://doi.org/10.5281/zenodo.5776772), to reproduce the generation of the comparative atlas or directly inspect it. The SCORPiOs LORelEi extension is available at GitHub (https://github.com/DyogenIBENS/SCORPIOS) and has also been deposited in Zenodo (https://doi.org/10.5281/zenodo.6913688), along with all input data, environments, and outputs. In addition, both the FishComparativeAtlas pipeline and SCORPiOs LORelEi extension are available in the Supplemental Material.

Data access

Gene homology relationships and local synteny conservation between teleosts can be interactively browsed through the Genomicus web server (Nguyen et al. 2022), accessible at https://www.genomicus.bio.ens.psl.eu/genomicus-fish-03.01/cgi-bin/search.pl. The homology relationships from the comparative atlas can be downloaded as flat or HTML files via an ftp server (ftp://ftp.biologie.ens.fr/pub/dyogen/genomicus-fish/03.01/ParalogyMap), along with the gene trees in New Hampshire eXtended (NHX) format (ftp.biologie.ens.fr/pub/dyogen/genomicus-fish/03.01/protein_trees.nhx).