Stable genome structures in living fossil fishes

Assembly, annotation, and description of the chromosome-level gar genomes

We assembled high-quality, chromosome-level genome assemblies of the alligator gar (Atractosteus spatula) and longnose gar (Lepisosteus osseus) by combining long- and short-read sequencing data with Hi-C. A k-mer survey (k = 17) revealed estimated genome sizes of 1.18 Gb for A. spatula and 1.19 Gb for L. osseus (Supplemental Fig. S1A). We generated long-read sequences to assemble the genomes into contigs, resulting in 1316 sequences with a contig N50 of 5.16 Mb for A. spatula, and 747 sequences with a contig N50 of 28.18 Mb for L. osseus. After polishing and quality improvement, we used Hi-C data to construct the chromosome-level genome assemblies, oriented and joined contigs into super-scaffolds, and manually curated them (Supplemental Table S1).

Our final genome assembly of Atractosteus spatula was 1.02 Gb, with a scaffold N50 of 56 Mb. Approximately 98.12% of the assembled genome was anchored into 28 chromosomes (2n = 56), consistent with karyotyping (Supplemental Fig. S1B,C). These chromosomes consisted of 16 macrochromosomes (≥20 Mb) and 12 microchromosomes (<20 Mb) (Supplemental Fig. S1D). The discrepancy between the estimated and final assembled genome sizes was likely explained by an overestimation of repetitive content (Wang et al. 2024), as indicated by the prominent repeat peak in the k-mer analysis (Supplemental Fig. S1A). The final genome assembly of Lepisosteus osseus was 1.2 Gb with a scaffold N50 of 55 Mb. Approximately 90% of the assembled genome was anchored into 28 chromosomes (2n = 56), consistent with previous karyotyping (Ráb et al. 1999). The new assemblies consist of 15 macrochromosomes (≥20 Mb) and 13 microchromosomes (<20 Mb) (Supplemental Fig. S1E). These chromosome numbers differ from a previously reported number for the L. oculatus (n = 29) (Braasch et al. 2016). Whole-genome and gene synteny maps show that Chromosomes 14 and 17 in the L. oculatus (Braasch et al. 2016) are fused in A. spatula and L. osseus (Fig. 1A; Supplemental Fig. S1F). Through sequencing depth analysis, we confirmed the accuracy of the identified chromosome fusion events. By aligning the sequencing data of L. oculatus and A. spatula to Chromosome 3 of A. spatula, we found that, at the fusion region, the sequencing depth of a 2-kb fragment in L. oculatus was zero, whereas A. spatula showed a relatively high sequencing depth in the fusion region. Further genomic feature analysis of the fusion region revealed a high proportion of repetitive sequences. We speculate that this region may represent the centromeric region, which likely contributed to the chromosomal changes (Supplemental Fig. S2). By incorporating homology-based, transcriptome and de novo annotation, we identified 21,055 genes in the A. spatula genome and 20,178 genes in the L. osseus genome and successfully annotated 96.94% of the A. spatula genes and 97.8% of the L. osseus genes to gene function databases (Supplemental Tables S2, S3).

Figure 1.

Chromosome-level genome synteny and characterization of the macro- and microchromosomes in gar species. (A) Comparison of whole-genome alignment synteny segments among A. spatula, L. osseus, and L. oculatus. (B) Comparison of interaction decay exponents between macrochromosomes and microchromosomes in A. spatula and L. osseus. (C) Comparison of substitution rates between macrochromosomes and microchromosomes in gars. (D) Comparison of macrosynteny conservation across diverse phylogenetic lineages. The conservation index (CI), Oxford dot plots, and gene synteny plots all indicate that A. spatula microchromosomes exhibit higher conservation than macrochromosomes throughout evolution. Colored highlights indicate microchromosomes.

Open in new tabDownload PowerPointLink to figure

Our new genome assemblies are more complete than previously published gar genomes (Braasch et al. 2016; Mallik et al. 2023). The scaffold length distributions of our gar genome assemblies skew toward substantially longer sequences, and the genome evaluation metrics markedly surpass those of earlier published genomes (Supplemental Fig. S1G–I; Supplemental Table S4). Using Benchmarking Universal Single-Copy Orthologs (BUSCO) analysis, we identified higher numbers of single-copy orthologs in our Atractosteus spatula and Lepisosteus osseus assemblies compared to previously published gar and bowfin genomes (Supplemental Fig. S1J). These analyses establish that the new A. spatula and L. osseus genomes are the most complete holostean assemblies currently available.

Gars have a stable chromosomal architecture

Despite diverging more than 100 million years ago (Grande 2010; Brownstein and Lyson 2022; Brownstein et al. 2023), Atractosteus spatula and Lepisosteus osseus exhibit an extraordinary degree of genomic similarity. The three gar genomes exhibit high synteny (Fig. 1A) with bowfins as well (e.g., Amia calva), the living sister lineage of gars. Gar genomes maintain higher synteny with tetrapod genomes (e.g., Xenopus laevis) than with teleost fish (e.g., Danio rerio), despite sharing a more recent common ancestry with teleosts (Supplemental Fig. S1K). Nonreference whole-genome alignment using Progressive Cactus (Armstrong et al. 2020) revealed up to 83.35% sequence identity between A. spatula and the most recent common ancestor of L. osseus and L. oculatus (Fig. 1A; Supplemental Table S5). Comparative genomic analysis revealed distinctive characteristics of macrochromosomes and microchromosomes in A. spatula, L. osseus, and L. oculatus. Macrochromosomes contain a higher total number of genes but demonstrate lower gene density per megabase, whereas microchromosomes exhibit higher GC content (Supplemental Figs. S1L–P, S3A–C). The proportion of repeat sequences is 31.24%, 37.39%, and 20.23% for A. spatula, L. osseus and L. oculatus, respectively (Supplemental Table S6).

Chromosomes 14 and 23 of A. spatula have exceptionally high repeat content, and large segments of these chromosomes maintain no synteny to the genomes of L. osseus and L. oculatus. Analysis of repeat sequence divergence among the chromosomes of A. spatula showed that Chromosomes 14 and 23 have lower average repeat divergence compared to other chromosomes (Supplemental Figs. S3A, S4). This suggests that the high number of GC repeats along Chr 14 and Chr 23 of A. spatula appeared recently in the evolutionary history of Atractosteus.

Some vertebrate genomes, including those of birds, turtles, squamates, sharks, and certain ray-finned fishes, consist of both macrochromosomes and microchromosomes (Waters et al. 2021). To investigate the features and functions of gar microchromosomes, we compared the characteristics of macrochromosomes and microchromosomes in A. spatula and L. osseus. The new assemblies revealed significant differences between these two chromosome types in gars: microchromosomes tend to have higher GC content, lower repeat sequence density, and higher gene density than macrochromosomes (Supplemental Fig. S5A–F). Differences in genomic features between macro- and microchromosomes are not unique to gars. They are also exhibited by other species with both chromosome types, such as birds (e.g., chicken) and cartilaginous fishes (Waters et al. 2021).

Hi-C data analysis showed that gar microchromosomes exhibit a slower decay rate and a higher degree of interchromosomal contact than macrochromosomes (Fig. 1B; Supplemental Fig. S5G,H), suggesting that microchromosomes are more structurally stable. Additionally, we infer that gar microchromosomes have a significantly lower substitution rate than macrochromosomes (Fig. 1C).

We also compared macrosynteny conservation across vertebrates. Despite millions of years of divergence (Grande 2010; Brownstein et al. 2022a), gar microchromosomes exhibit more conserved synteny than macrochromosomes when compared with other species through whole-genome synteny analysis. Macrochromosomes showed frequent chromosomal rearrangements in vertebrates, particularly during the teleost genome duplication (TGD). By calculating the conservation index (CI), which measures the ratio of one-to-one ortholog pairs on homologous chromosomes to one-to-one ortholog pairs on nonhomologous chromosomes (Wang et al. 2017; Simakov et al. 2020, 2022; Li et al. 2022), we found that microchromosomes exhibit higher conservation across the whole genome (Fig. 1D; Supplemental Table S7; Supplemental Figs. S6–S8).

To investigate gene function on gar microchromosomes, we calculated Ka/Ks ratios (the ratio of the number of nonsynonymous substitutions per nonsynonymous site Ka to the number of synonymous substitutions per synonymous site Ks) (Hurst 2002; Zhang 2022) of paralogous genes between macro- and microchromosomes. Ka/Ks values for genes on gar microchromosomes are normally distributed (Ka/Ks < 1) around 0.5, suggesting that these genes are under purifying selection (Supplemental Fig. S5I; Hurst 2002; Roth and Liberles 2006; Fuselli et al. 2023). We also compared gene expression levels and found that genes on microchromosomes in gars have higher expression levels than those on macrochromosomes (P < 0.05) (Supplemental Fig. S5J). We also identified 91 positive selection genes for DNA repair and genome stability in gar microchromosomes (see Supplemental Results; Supplemental Fig. S12; Supplemental Table S10). Finally, we performed GO functional enrichment analysis on gar microchromosome genes, which showed significant enrichment in “death receptor activity” and “tumor necrosis factor receptor activity” (Supplemental Fig. S5K). The enriched genes are members of the tumor necrosis factor receptor superfamily and play a crucial role in apoptosis (programmed cell death) (Ashkenazi and Dixit 1998). For example, the CD40 gene is essential for mediating a wide range of immune and inflammatory responses and can induce cell death (Hess and Engelmann 1996). Similarly, the NGFR gene can mediate both cell survival and cell death in neural cells (Bruno et al. 2023). This suggests that microchromosomes may be crucial for maintaining homeostasis in gars by regulating cell death and inflammatory processes.

Evolutionary stasis of gar genome sequences and structures

Gars and bowfins are emerging as a model system to understand the molecular evolution of living fossil lineages (Braasch et al. 2016; Thompson et al. 2021; Brownstein et al. 2022a, 2024; Simakov et al. 2022). To better understand rates of genomic evolution in gars, we estimated nucleotide substitution and genome structural change rates using the new Atractosteus spatula and Lepisosteus osseus assemblies. We constructed a phylogenetic tree from 43,420,318 loci after conducting reference-free whole-genome alignment using Progressive Cactus (Armstrong et al. 2020). Whereas the coelacanth lineage displayed the shortest branch length, indicative of an exceptionally slow substitution rate (Nikaido et al. 2013), our primary focus, the gars, also clearly exhibited relatively short branch lengths and lower substitution rate compared to other sampled chordates (Fig. 2A).

Figure 2.

Evolutionary rate of gars compared with diverse phylogenetic lineages. (A) Comparison of whole-genome nucleotide substitution rates using reference-free whole-genome alignments. (B) Heat map showing microsynteny conservation within major vertebrate clades, with gene copy numbers of orthogroups labeled in different colors. (C) Box plot compares the conservation index of major vertebrate clades to that of gars. (D) Chromosomal rearrangements across representative species of major vertebrate clades. Phylogenetic branches are colored according to the average rate of intrachromosomal rearrangements. The heat map displays the number of intrachromosomal rearrangements per chromosome for each species compared to Narcine bancroftii, with varying colors indicating different frequencies. The lollipop chart indicates that gars, Amia calva, and Acipenser ruthenus exhibit a higher syntenic length ratio with Narcine bancroftii compared to other species. This suggests that gars have experienced more significant sequence loss.

Open in new tabDownload PowerPointLink to figure

We also observed that gar genomes possess fewer nucleotide insertions and deletions than other species of chordates, with differences in single base pair deletions among the three gar species possibly accounting completely for their varying genome sizes (Supplemental Fig. S9A,B). We also examined whether evolutionary rates differ meaningfully among gars and other canonical living fossil lineages. To do this, we calculated the relative evolutionary rate in 11 species of putative living fossil animals using single-copy orthologs. Our results revealed that gars exhibit a similar evolutionary rate to most other “living fossils” and that relative evolutionary rate is positively correlated with phylogenetic branch length (Supplemental Fig. S9C; Supplemental Tables S8, S9).

Because gars diverged prior to the whole-genome duplication event (Braasch et al. 2016; Thompson et al. 2021) that occurred along the teleost stem lineage (Bi et al. 2021; Davesne et al. 2021), we suspected that gars have a relatively conserved genome architecture compared with teleosts. Using the whitespotted bamboo shark (Chiloscyllium plagiosum) as a reference species (Zhang et al. 2020), we performed syntenic gene retention analysis, which revealed a strong negative correlation between synteny retention rates of gars and teleost fishes and their corresponding branch lengths along our inferred phylogeny. Gars, which are subtended by the shortest branch length, had the highest number of collinear genes (Supplemental Fig. S9D). Additionally, analysis of syntenic block size (measured by the number of conserved anchor genes within blocks) indicated that gars exhibit the slowest decay rate among representative teleost fishes, suggesting that gars have a stable genome architecture accompanied by a slow sequence substitution rate (Supplemental Fig. S9E).

We also profiled microsynteny conservation among major vertebrate groups to further understand the nature of structure change in gar genomes. The results showed that the overall syntenic conservation of homologous genes in gars is higher than other vertebrates and far higher than in teleost fishes (Fig. 2B). We also calculated values of the macrosynteny conservation index (Wang et al. 2017; Simakov et al. 2020, 2022; Li et al. 2022) for gars and other jawed vertebrates. The resulting CI values showed that teleost fishes show the lowest degree of microsynteny with gars, whereas birds and nonbird reptiles show the highest (Fig. 2C; Supplemental Fig. S9F). Our analysis of chromosome rearrangement among gars and other vertebrate species showed that, compared to other species, gars have the lowest number of intrachromosomal rearrangements, with fewer rearrangements in gar microchromosomes than in macrochromosomes. Gars have experienced an average of 0.5 intrachromosomal rearrangement events per million years (Fig. 2D; Supplemental Fig. S10). Collectively, these results demonstrate that gars have maintained a higher degree of genomic structural stability and sequence conservation than other vertebrates.

Abundance and activity of transposable elements in gar genomes

Transposable elements (TEs) are DNA sequences in the genome that can replicate themselves or move between chromosomes. TE activity can lead to genomic instability (Bourque et al. 2018; Bhat et al. 2022). We compared the abundance and activity of TEs in gars and other vertebrate species to examine what role TE evolution might play in the stagnant evolution of gar genomic sequences and chromosomal structures. Based on the Kimura distance (K2P) of TE divergence from consensus sequences, the TE landscape in Atractosteus spatula shows two expansion peaks. The oldest peak, located on the right side of the graph (Fig. 3A) represents highly divergent copies resulting from mutation accumulation, corresponding to inactive TEs. We compared these older peaks in A. spatula with those in other species and found that gars have the largest proportion of inactive TEs, including DNA, LINE, SINE, and LTR families. To quantify the activity and expansion rate of younger TEs, we focused on families with <10% divergence from their consensus sequences. Using the species phylogeny, divergence time, and genome size, we reconstructed the ancestral state of young TE insertions. The results showed that A. spatula has the lowest insertion rate among examined vertebrates (Fig. 3A; Supplemental Fig. S11). Additionally, we analyzed the copy number distribution of young TE families, revealing that gars exhibit the lowest copy numbers among vertebrates (Fig. 3B–D). Together, these results indicate that TE activity in gars is lower than in other species and likely contributes to their genomic stability.

Figure 3.

Repeat sequences evolutionary rate of gars compared with diverse chordate lineages. (A) Phylogenetic branches are colored according to the evolutionary rate of recently accumulated transposable elements (TEs) with <10% Kimura divergence from consensus TEs across diverse species. The plot on the right shows the distribution of transposable elements based on Kimura distance analysis in these species. Kimura divergence of 25 marks the major ancient peak in A. spatula and is used for comparison. (B–D) Comparison of DNA/LINE/LTR transposable element copy numbers with <10% Kimura divergence from consensus sequences across diverse species.

Open in new tabDownload PowerPointLink to figure

Consequences of genomic stasis on gar population structure

Hybridization of Atractosteus spatula, Lepisosteus osseus, and L. oculatus in the wild has been confirmed using microsatellite and ddRAD single nucleotide polymorphism (SNP) data (Herrington et al. 2008; Bohn et al. 2017; Taylor et al. 2020; Brownstein et al. 2024). The existence of F1 and F2 A. spatula × Lepisosteus spp. crosses is surprising because these genera diverged from one another during the Cretaceous period, between 97 and 120 million years ago (Grande 2010; Brownstein et al. 2022a). The slow rate of genomic evolution in gars may enable gars to hybridize across deep divergences that would otherwise correspond to increased genetic divergence and resulting hybrid inviability from the emergence of incompatibilities (Bolnick 2005; Matute et al. 2010; Coughlan and Matute 2020). The maintenance of species barriers in gars has resultantly been tied to differences in reproductive behavior and fecundity (Brownstein et al. 2024). Because gar species in the two genera have existed in sympatry in North America for at least 50 million years (Brownstein et al. 2024) and there is currently wide overlap between most living species (Grande 2010), it is also surprising that little evidence exists so far for historical episodes of introgression between Atractosteus and Lepisosteus (Brownstein et al. 2024).

To investigate how population structure and species divergence among gars might be affected by their slow rates of genomic evolution, we conducted principal component (PCA) and phylogenetic analyses, which assigned Lepisosteus oculatus, L. osseus, and Atractosteus spatula populations to three distinct clusters and inferred that these lineages are reciprocally monophyletic, respectively (Fig. 4A; Supplemental Fig. S13A). Admixture analysis of 15,305,348 SNPs grouped our sample of 69 gars into three clusters (K = 3) corresponding to the three included species. However, we also find evidence for admixture among species (Fig. 4B). Identity-by-descent (IBD) analysis revealed that some individuals of L. osseus and L. oculatus also share haplotype blocks (Supplemental Fig. S13B), and F3 (Reich et al. 2009; Patterson et al. 2012) statistics indicate that L. osseus has experienced admixture with both A. spatula and L. osseus (F3 < 0, Z-score < −3) (Fig. 4D; Supplemental Table S11). Treemix identified one migration event between A. spatula and L. osseus across values of m (1 to 5), in which m is the weight of the migration matrix and specifies the number of potential gene flow events (Fig. 4E). We also inferred gene flow between these gar species using simulations in Momi2 (Supplemental Fig. S13C–E). Nonetheless, calculated pairwise fixation index (F_ST) values among A. spatula, L. osseus, and L. oculatus, which measure the number of alleles shared by different populations as proportions of alleles shared in a single population, range from 0.384 between L. osseus and L. oculatus to 0.785 between A. spatula and L. oculatus (Fig. 4C); values >0.25 are classically considered to be indicative of comparisons among species (Wright 1965). These results indicate that, despite strong differentiation and long divergence times among the three gars, ancient genetic introgression events have still occurred over their long history of sympatry.

Figure 4.

Population structure and ancient introgression of living gar genera. (A) PCA analysis of population data from three gar species. (B) ADMIXTURE results with K values ranging from 2–4 with three (K = 3) found as the optimal number of populations. (C) Average F_ST and π values among the three gar species. (D) F3 statistics for L. oculatus, L. osseus, and A. spatula. (E) Treemix trees with one migration events (m = 1).

Open in new tabDownload PowerPointLink to figure

Ethics statement, sampling, and sequencing

All sampling and experimental procedures were approved by the Ethics Committee of the Institute of Hydrobiology, Chinese Academy of Sciences. We collected A. spatula and L. osseus specimens from an aquarium market in Wuhan, Hubei Province, China. We subjected the obtained samples to high-throughput Illumina sequencing, Pacific Biosciences (PacBio) HiFi, ultralong Oxford Nanopore Technologies (ONT), Hi-C, and RNA sequencing. For details, see the Supplemental Methods.

Genome assembly and evaluation

High-throughput Illumina data were employed for k-mer-based genome size estimation. PacBio HiFi and ONT ultralong reads were used to generate draft assemblies, whereas Hi-C data anchored contigs into chromosome-level genomes. For details, see the Supplemental Methods.

We evaluated the completeness of the A. spatula and L. osseus genome assemblies using BUSCO (v4.1.2) (Simão et al. 2015) and visualized the synteny of BUSCO genes using ChrOrthLink (https://github.com/chulbioinfo/chrorthlink) (Yoo et al. 2022). We assessed genomes by mapping Illumina short reads to the assembled genomes with BWA (Li and Durbin 2009) using default parameters. We also mapped all transcripts assembled from RNA-seq reads to the genome using BLAST (v35.1) (Kent 2002).

Repeat and protein-coding gene annotation

For repeat sequence annotation, we used a combination of a custom de novo repeat sequence library and a public repeat sequence database to predict repeat elements in assembled genome. For gene prediction, we used a combination of ab initio gene prediction, homology-based gene prediction, and transcript-based prediction to identify protein-coding genes in gars. For details, see the Supplemental Methods.

Whole-genome alignment and genome synteny

The genomes of three gar species (A. spatula, L. osseus, L. oculatus) were aligned using Progressive Cactus (v2.8.0) (Armstrong et al. 2020). To prepare this alignment, we soft-masked genomes of the three gars with RepeatMasker (Tarailo-Graovac and Chen 2009). The guide tree and divergence times used as input for Progressive Cactus were retrieved from the TimeTree database (Kumar et al. 2017). Based on the alignment results, we calculated nucleotide similarity levels among the three gar species using A. spatula as the reference. For chromosomal-level synteny analysis, we used NGenomeSyn (v1.41) (He et al. 2023).

Gar macro- and microchromosome analysis

To compare differences between gar macrochromosomes and microchromosomes, we first calculated the GC content, repeat sequence content, and gene density using a 1-Mb sliding window in BEDTools (v2.31) (Quinlan and Hall 2010). We then analyzed chromosome interaction frequency and interaction decay as a function of increasing distance. Hi-C paired-end reads were processed using HiC-Pro (v3.1.0) (Servant et al. 2015), and insulation scores, Pearson's correlation matrices, and eigenvectors were calculated and visualized using FAN-C (v0.9.26) (Kruse et al. 2020). The observed-to-expected interchromosomal matrix was calculated by counting normalized interchromosomal interactions in the 1-Mb KR-normalized matrix. Interaction frequency and distance decay were analyzed with a 500-kb sliding window, and the decay curve was plotted using ggplot2 (https://ggplot2.tidyverse.org) in R (v4.2.2) (R Core Team 2014).

To compare substitution rates between macrochromosomes and microchromosomes, we aligned the genomes of the three gar species using Progressive Cactus (v2.8.0) (Armstrong et al. 2020) with A. spatula as the reference genome. Based on alignment results, we removed repeat regions and coding regions, then split the chromosomes into 100-kb segments and extracted region sequences from the whole-genome alignment results. We calculated substitution rates for each region using r8s (v1.81) (Sanderson 2003).

To examine the conservation of microchromosome organization in gars, we conducted macrosynteny analysis between gars and other chordate species (Branchiostoma floridae, Petromyzon marinus, Amblyraja radiata, Acipenser ruthenus, Oryzias latipes, Naja naja, Chelonia mydas, Gallus gallus, Latimeria chalumnae, Polypterus senegalus, and Amia calva) using LAST (v1542) (Frith 2011) and PanSyn (Yu et al. 2024), and visualized macrosynteny with Circos (v0.69) (Krzywinski et al. 2009). To quantify the degree of chromosome conservation, we calculated the conservation index for each comparison between A. spatula and other species by dividing the number of one-to-one orthologous gene pairs on homologous chromosomes by the number of one-to-one orthologous genes on nonhomologous chromosomes. We used Fisher's exact test (P < 0.05) to assess the significance of chromosome conservation. The conservation index ranges from 0 (no conservation) to 1 (highly conserved) (Simakov et al. 2013; Li et al. 2022).

We used KaKs_Calculator (v2.0) (Zhang 2022) to calculate the Ka/Ks ratios of paralogous gene pairs between macrochromosomes and microchromosomes, and gene expression values for L. osseus RNA-seq data were calculated using featureCounts (v2.0.2) (Liao et al. 2014). The positive selection analysis was designed to identify genes putatively enriched on microchromosomes and was accordingly restricted to taxa possessing them. Thus, the bowfin, despite being a sister group to gars, was excluded due to the absence of microchromosomes (Thompson et al. 2021).

Genomic nucleotide substitution and chromosomal rearrangement rates

To estimate the genome-wide nucleotide substitution rates for gars, we first conducted whole-genome alignment to compare the evolutionary rates of gars with those of other animals (including Amphioxous, Agnathans, Chondrichthyans, Nonteleost actinopterygians, teleostei, Aves, Suqmata, Testudines, and Mammals) using Progressive Cactus (v2.8.0) (Armstrong et al. 2020). We then processed the alignment results using the commands halStats and halSummarizeMutations (Hickey et al. 2013) and plotted the phylogenetic tree, branch lengths, and substitution rates per site using ggtree (v3.14) (Yu et al. 2017). We calculated the relative evolutionary rate of several “living fossils” (including Ornithorhynchus anatinus, Ailuropoda melanoleuca, Crocodylus porosus, Dromaius novaehollandiae, Acipenser ruthenus, Polypterus senegalus, Carcharodon carcharias, Branchiostoma floridae, Tachypleus tridentatus, Eptatretus atami, Nautilus pompilius, and Latimeria chalumnae) using two-cluster analysis and Tajima's relative rate test (Lü et al. 2021). In the two-cluster analysis, we assessed faster or slower evolutionary rates in specific taxa using Z-statistics and the tpcv command in LINTRE (Takezaki et al. 1995). Using Tajima's relative rate test, we measured lineage-specific substitutions and validated faster evolutionary rates using χ² tests.

To investigate synteny loss rates in gars, we aligned 35 ray-finned teleost fishes to Chiloscyllium plagiosum using MCScan in the JCVI toolkit (Tang et al. 2024). We identified synteny blocks by the number of anchor gene pairs in each block, with the relative synteny retention rate calculated as the total number of genes in the reference species that have syntenic relationships with the compared species (Shi et al. 2022). We also calculated Pearson's correlation between collinear genes and branch lengths based on the phylogenetic tree. We compared the microsynteny conservation of gars with other vertebrates (including 14 teleostei, 11 Aves, five Lissamphibia, five Testudines, seven Squamata, one Alligator, and 14 Mammalia using Syntenet (Almeida-Silva et al. 2023).

For genome rearrangement analysis, we retained only chromosome sequences from each species and performed whole-genome alignment using LAST (v1542) (Frith 2011), with the genome of Chelonia mydas as the reference. We merged small syntenic blocks into larger blocks and discarded blocks shorter than 50 kb. We identified genome rearrangements (inversions, translocations, and inverted translocations) based on the orientation and position of retained blocks (Liu et al. 2021). To infer the rearrangement rate at each phylogenetic tree node, we estimated ancestral chromosome rearrangements using a maximum-likelihood approach along our time-calibrated phylogeny in the R package phytools (Revell 2012; R Core Team 2014).

Transposable element evolution

To identify transposable elements in the genomes of the three studied gar species (L. osseus, L. oculatus, A. spatula), Danio rerio, Acipenser ruthenus, Amia calva, Polypterus senegalus, Thamnophis elegans, Xenopus laevis, Amblyraja radiata, Petromyzon marinus, and Branchiostoma floridae, we created a de novo TE database for each genome using RepeatModeler (Flynn et al. 2020) with the -LTRStruct module. We combined this custom database with a public repeat database and ran RepeatMasker (Tarailo-Graovac and Chen 2009) with parameters -a -s -c -gccalc to calculate Kimura values for all identified TEs using the script calcDivergenceFromAlign.pl. We then generated TE landscapes using the script createRepeatLandscape.pl. In addition, because our repetitive sequence annotation method only relies on RepeatModeler and RepeatMasker, it may lead to incomplete annotation of repetitive sequences and omission of certain lineage-specific TEs.

TE landscape peaks represent the most abundant and active bursts of TEs in each species (Gozashti et al. 2023; Decena-Segarra and Rovito 2024). TEs located on the left side of the landscape indicate insertions with low divergence from their respective consensus sequences, suggesting recent or active TEs. In contrast, TEs on the right side show higher divergence, indicating an accumulation of mutations, corresponding to older or inactive TEs. To quantify TE activity, we selected TE insertions with a K2P genetic distance of <10% compared to their respective consensus sequences. We retrieved a phylogeny and divergence times from the TimeTree database (Kumar et al. 2017) and inferred active TE insertion states at ancestral nodes using maximum-likelihood analysis with the phytools package (Revell 2012), applying the fastAnc function to reconstruct ancestral states according to divergence times and genome size. Additionally, we counted the active TE copies for each species and generated a heat map for comparison.

Population structure and admixture among gars

To examine population structure and infer episodes of historical introgression among living gar species, we collected single nucleotide polymorphism data for individuals of A. spatula, L. osseus, and L. oculatus. For details on population analysis methods, see the Supplemental Methods.