Evolutionary stasis characterizes lineages that seldom speciate and show little phenotypic change over long stretches of geological time. Although lineages that appear to exhibit evolutionary stasis are often called living fossils, no single mechanism is thought responsible for their slow rates of morphological evolution and low species diversity. Some analyses of molecular evolutionary rates in a handful of living fossil lineages have indicated they exhibit slow rates of genomic change. Here, we investigate mechanisms of evolutionary stasis using a dataset of 1,105 exons for 481 vertebrate species. We demonstrate that two ancient clades of ray-finned fishes classically called living fossils, gars, and sturgeons, exhibit the lowest rates of molecular substitution in protein-coding genes among all jawed vertebrates. Comparably low rates of evolution are observed at four-fold degenerate sites in gars and sturgeons, implying a mechanism of stasis decoupled from selection that we speculate is linked to a highly effective DNA repair apparatus. We show that two gar species last sharing common ancestry over 100 million years ago naturally produce morphologically intermediate and fertile hybrids. This makes Gars the oldest naturally hybridizing divergence among eukaryotes and supports a theoretical prediction that slow rates of nucleotide substitution across the genome slow the accumulation of genetic incompatibilities, enabling hybridization across deeply divergent lineages and perhaps slowing the rate of speciation. Our results help establish molecular stasis as a barrier to speciation and phenotypic innovation and provide a mechanism to explain the low species diversity in living fossil lineages.

Molecular evolutionary rates among major lineages of jawed vertebrates.

To test whether slow rates of morphological evolution are paired with low rates of molecular evolution in living fossils like gars, sturgeons, and paddlefishes, we estimated molecular rate variation across 1,105 exons from a sample of 471 jawed vertebrate species. We identified orthologous exon sequences from the genomes of selected species in the NCBI database for the following major jawed vertebrate lineages: Acipenseriformes, Aves, Crocodylia, Chondrichthyes, Lepidosauria, Lissamphibia, Marsupialia, Placentalia, Polypteridae, and Teleostei (Figure S1a). The HMM protocol available in HMMER 3.1 (Wheeler and Eddy 2013) was used to search each of the downloaded genomes for orthologous exons. These exon sequences were extracted using Python scripts from a phylogenomic analysis of ray-finned fishes using these loci (Hughes et al. 2018). For each group of the exon sequences, we aligned them using MAFFT v.7.3 (Katoh and Standley 2013) with default parameters. Exon alignments for the seven living species of gars (Lepisosteidae) used in a phylogenomic study (Brownstein et al., 2023) were included in the comparative analysis. Each exon was separately aligned among the species in a given vertebrate lineage, resulting in a maximum of 1,105 alignments sampled for each lineage. Fourfold degenerate (4D) sites were extracted from all exons and concatenated in every vertebrate clade except Polypteridae and Acipenseriformes. This was because we could not find orthologous 4D site sequences for all the available genome assemblies for these two clades; sampling 4D sites for all three species would be needed to sample the common ancestors of Polypteridae and Acipenseriformes, as we were only able to include three species of each in our exon rate estimate analyses.

We estimated and compared posterior molecular substitution rates at each exon across all major vertebrate clades using fixed input trees in Bayesian molecular clock analyses. We used previously published time-calibrated phylogenies for Teleostei (Hughes et al. 2018), Acipenseriformes (Kumar et al. 2017), Polypteridae (Near et al. 2014), Lepisosteidae ( Brownstein et al. 2023), Chondrichthyes (Kumar et al. 2017), Testudines (Shaffer et al. 2017), Amphibia (Kumar et al. 2017), Lepidosauria (Pyron and Burbrink 2014), Aves (Prum et al. 2015), Crocodylia (Green et al. 2014), Marsupialia (Upham et al. 2019), and Placentalia (Upham et al. 2019). The time tree of vertebrates used in the branch rate analysis of coelacanths and lungfish was taken from the literature (Wang et al. 2021) and timetree.org (Kumar et al. 2017). We used these time-calibrated molecular phylogenies in BEAST 2.5.2 (Bouckaert et al. 2019) by inserting them in Newick format into the ‘Starting Tree’ tab in the BEAUTi terminal. The following operators were turned off to ensure the input tree remained fixed: tree scaler, tree root scaler, uniform operator, subtree slide, narrow and wide exchange, and Wilson-Balding. In turn, we neither estimated tree topology nor divergence times. Custom scripts for inserting the trees, along with XML files containing the tree topologies used, are in the Supplementary Data.

BEAST 2.5.2 (Bouckaert et al. 2019) was used to estimate the Bayesian posterior nucleotide substitution rate for each of the 1,105 exons and fourfold degenerate sites separately from each vertebrate clade. The computer program BEAUTi (Bouckaert et al. 2019) was used to construct individual XML files from each exon alignment and the pooled fourfold degenerate sites for each clade with the clade-specific time-calibrated phylogeny. The time-calibrated phylogeny was fixed such that BEAST did not estimate topology or divergence times. Because of the large number of BEAST analyses, we took advantage of the Yale High-Performance Computing cluster and built custom scripts to produce and run XML files along this pipeline. First, we produced a template XML for each clade that specified all input parameters regarding a single XML. We then used a custom batch script (Supplementary Information) to fit the specifications given by the template XML file with each gene to produce individual XML files for every single exon sampled in a given vertebrate clade. We used a Yule (pure-birth) branching model as the tree prior, a relaxed Lognormal clock model as the clock before allowing independent rates for each branch, and an HKY model of nucleotide evolution to allow for unequal frequencies and transition rates, and ran each analysis ran for 10 million generations. Upon completion of the analyses, we confirmed sufficient MCMC mixing (ESS>200) for each BEAST run using the program Tracer v. 1.7 (Rambaut et al. 2018) and the R package “coda” v.0.19-4 (Makowski et al. 2019).

To test whether our rate estimates were unaffected by tree model choice, we reran exons in lepidosaurs and teleosts with the five fastest, five slowest, and five middlemost rate estimates under a Birth-Death tree model. Theoretically, tree model choice should not affect our results because we fixed the time tree in each analysis, but we chose to test this outright. Parameter choices were otherwise the same as the original runs. We then compared absolute rate estimates between runs using the Yule and Birth-Death models. Next, we tested to see whether the number of species sampled for different clades biased estimated rates. Among the six clades with more than seven species sampled that had average estimated substitution rates higher than the average in gars (Figure 2), we subsampled seven species of lepidosaurs and teleosts (the number of species sampled for gars) that captured the common ancestry of major subclades (i.e., Squamata, Toxicofera, Euteleostei, Acanthomorpha) out of our exon dataset and reran the analyses under the original parameter specification (i.e., under a Yule model). We then compared absolute rate estimates between runs using the full and reduced sampling of lepidosaurs and teleosts.

Estimation of branch-specific molecular evolutionary rates for living fossils.

We investigated the rates of molecular substitution in candidate living fossil lineages including paleognathous birds, the Hoatzin Opisthocomus hoazin, the Tuatara Sphenodon punctatus, the Salamanderfish Lepidogalaxias salamandroides, the African Coelacanth Latimeria chalumnae, the Australian Lungfish Neoceratodus forsteri, and the West African Lungfish Protopterus annectens. We extracted the estimated posterior molecular substitution rates for the corresponding terminal branches in our fixed time trees. For the two lungfishes and African Coelacanth, we constructed a time tree that included these taxa and the species of Lissamphibia for which exon data was available. We included Lissamphibians to ensure that the time tree accommodated the paraphyly of sarcopterygian “fishes.” This was accomplished by taking the fixed tree of Lissamphibia (see ‘Input Tree Selection’) and manually adding Latimeria chalumnae, the Australian Lungfish Neoceratodus forsteri, and the West African Lungfish Protopterus annectens as progressive outgroups to Amphibia (following (Amemiya et al. 2013; Meyer et al. 2021; Wang et al. 2021)) and fixing divergence times following TimeTree.org searches. BEAST analyses were run using the same parameters as in the whole-clade rate estimation runs and extracted branch-specific rates for L. chalumnae and Dipnoi.

Hybridization among deeply divergent gar lineages Lepisosteus and Atractosteus.

Next, we investigated signatures of hybridization across deeply divergent gar lineages to test whether slow rates of molecular evolution might be associated with incomplete reproductive isolation across deep time in living fossil lineages. To check for regions where both extant gar genera are currently sympatric, we downloaded occurrence data for species of Lepisosteus and Atractosteus spatula from FishNet2 (http://fishnet2.net/). We pruned erroneous and duplicate records. The clean occurrence data files are included in the Supplementary Data.

Natural hybrids have been reported among wild populations of Atractosteus spatula and Lepisosteus osseus in Texas and Oklahoma (Bohn et al. 2017; Taylor et al. 2020), and so we examined genomic evidence of hybridization in these populations. We obtained tissue samples from 206 specimens of Atractosteus spatula, Lepisosteus osseus, and hypothesized A. spatula X L. osseus hybrids from across the Gulf of Mexico coastal river systems (Table S1) to test for both the commonality of hybrids and the presence of both F1s and F2s in the Brazos river system. The Brazos and Trinity systems and Choke Canyon Reservoir were targeted because previous studies have convincingly demonstrated the presence of hybrid Atractosteus spatula X Lepisosteus osseus individuals in this region (Bohn et al. 2017; Taylor et al. 2019). We chose a subset of five A. spatula and seven L. osseus isolations with high DNA concentrations based on Qubit fluorometer (Life Technologies, Carlsbad, CA, USA) readings to test for naturally occurring hybrids between species in these genera. These were compared with tissue samples from individuals that showed morphological intermediacy consistent with being a hybrid, as well as their identification in a previous study (Bohn et al. 2017).

We used double digest restriction site associated DNA sequencing (ddRADseq) to obtain a large dataset of single nucleotide polymorphisms (SNPs) to investigate whether sympatric populations of Lepisosteus osseus and Atractosteus spatula show evidence of hybridization. ddRADseq was performed following a modified version of the frequently used protocol (Peterson et al. 2012). Additional details are described in the Supplementary Information. Pooled libraries were size-selected for 300-500 bps using a BluePippin sequenced by the University of Oregon Genomics & Cell Characterization Core Facility on an Illumina HiSeq 4000. We assembled the gar ddRAD dataset using iPyrad v.0.9.68 (Eaton and Overcast 2020) with the Lepisosteus oculatus genome (NCBI accession number: GCF_000242695.1) (Braasch et al. 2016) as a reference, resulting in a total of 256,750 loci shared by at least four samples. For analysis of hybridization, ddRAD SNPs were filtered using VCFTools v.0.1.15 (Danecek et al. 2011) so that only biallelic SNPs with minor allele counts >1 and <5% missing data were retained. To minimize the effects of linkage among markers, one random SNP per 10,000 bp window was retained, resulting in a dataset containing 2,097 SNPs (2.2% missing data). In addition, we used a custom R script to identify 1,223 SNPs (2.1% missing data) that were fixed between the parental species (Atractosteus spatula and Lepisosteus osseus). Demultiplexed ddRAD Illumina reads are deposited at the NCBI SRA (PRJNA1077910).

We assessed genomic signals of hybridization between Atractosteus spatula and Lepisosteus osseus. First, missing genotypes were imputed with the “impute” function (method = “random”) in the R package LEA v.3.4.0 (Frichot and François 2015). To examine patterns of genetic variation, we then performed principal component analysis (PCA) using the 2,097 filtered SNPs with the “dudi.pca” function in the R package ade4 v.2.1.4 (Dray and Dufour 2007; Thioulouse et al. 2018). In addition, we estimated genomic ancestry coefficients for each individual with sparse nonnegative matrix factorization implemented in the R package LEA. Ten replicate analyses were performed with between two and ten ancestral populations (K). We examined cross-entropy scores to determine the optimal value of K (Figure S3a).

Hybrid classification was performed using the dataset of 1,223 fixed SNPs. First, the “find.clusters” function from the Adegenet v.2.1.4 R package (Jombart and Ahmed 2011) was used to assign all individuals to two genetic clusters by performing 10 million search iterations of the K-means algorithm with 1,000 random starting centroids. Starting with the K-means group memberships, we used the Adegenet “snapclust” function (Beugin et al. 2018) to estimate the probability of membership to parental, F1, or backcross classes. A maximum of ten million generations for 1,000 replicate runs of the expectation–maximization were performed. Additionally, we estimated the ancestry index and interspecific heterozygosity for each individual using the R package HIest v.2.0 (Fitzpatrick 2012).

To better understand patterns of hybridization, including skewed sex ratios among parental Atractosteus spatula and Lepisosteus osseus, we interrogated mitochondrial DNA (mtDNA) data for a sample of 201 gars including all living species in Lepisosteidae. Because mitochondrial DNA is maternally inherited, inferring a phylogenetic tree using mitochondrial sequence data can illuminate skewed sex ratios in hybrid crosses. We included two specimens of Amia calva to serve as outgroups in the phylogenetic analysis. Sampling locations of specimens used for genetic analyses are listed in Supplemental File Table_S1.csv. DNA was extracted from 95% ethanol-preserved tissues using a standard DNeasy Qiagen Blood and Tissue Kit (QIAGEN, Valencia, CA, USA). To minimize downstream enzymatic inhibition, we purified DNA extractions with ethanol precipitation: 3M sodium acetate (pH = 5.2) was added equal to 10% of the total volume of the DNA extraction followed by 100% ethanol equal to 2.5 times the total volume of DNA. After mixing, extractions were incubated for 10 minutes at -80C. Samples were centrifuged for 30 minutes at 8,000 RCF, the supernatant was carefully poured off, and the DNA pellet was washed with 250 uL of cold 70% ethanol. Samples were centrifuged again for 5 minutes at 8,000 RCF, supernatant was poured off, the pellet was allowed to air dry for ~15 minutes, and the DNA pellet was resuspended with the desired amount of DNAse-free water.

The mtDNA gene tree of gars (Figure 5f) was inferred with a phylogenetic analysis of the mitochondrial encoded cytochrome b (cytb) gene. The molecular phylogenetic analysis included 33 specimens of Atractosteus spatula, two specimens of A. tristoechus, eight specimens of A. tropicus, 65 specimens of Lepisosteus osseus, five specimens of L. platostomus, 44 specimens of L. oculatus, 23 specimens of L. platyrhincus, 21 specimens of Atractosteus spatula X Lepisosteus osseus hybrids, and a single specimen of both Amia calva and A. ocellicauda to serve as outgroups. The cytb gene was amplified using previously published PCR primers and cycling conditions (Wright et al. 2012). Amplification products were prepared for DNA sequencing using a polyethylene glycol precipitation. Contiguous sequences were assembled from individual DNA sequencing reactions using the computer program Geneious v.7.2 (Kearse et al. 2012). New cytb sequences were aligned by eye to those previously generated in early studies of gar phylogeny (Wright et al. 2012). The optimal data partitioning scheme, among the three codon positions of the cytb gene, and molecular evolutionary models was determined using the Bayesian information criterion in the computer program Partitionfinder v. 2.1 (Lanfear et al. 2017). The mitochondrial gene tree was inferred from the aligned cytb sequences using the optimal molecular evolutionary models and partitioning scheme using the computer program MrBayes v. 3.2 (Ronquist et al. 2012), where posterior probabilities for the phylogeny and parameter values were estimated using Metropolis-couple Markov chain Monte Carlo (Larget and Simon 1999; Huelsenbeck et al. 2001). The MrBayes analysis was run for 10⁷ generations with two simultaneous runs each with four chains. Convergence of the MC3 algorithm and stationarity of the chains was assessed by monitoring the average standard deviation of the split frequencies between the two runs, which was less than 0.005 after 3 x 10⁶ generations. In addition, the likelihood score and all model parameter estimates were plotted against the generation number to determine when there was no increase relative to the generation number in the computer program Tracer v. 1.7 (Rambaut et al. 2018). The first 50% of the sampled generations were discarded as burn-in, and the posterior phylogeny was summarized as a 50% majority-rule consensus tree. All cytb gene sequences generated for this study are available at GenBank PP331004 - PP331204.

Geometric morphometric analyses of gar hybrids.

To quantify how the phenotypes of hybrid individuals of the two extant gar genera compared to those of their parental lineages, we used a dataset sampled from 25 specimens of Atractosteus spatula, Lepisosteus osseus, and A. spatula x L. osseus from the Brazos River system in Texas. The skull and mandible were selected as regions of study because there is morphological variation in these traits among lineages of extant and extinct gars (Wiley 1976; Kammerer et al. 2006; Grande 2010; Brito et al. 2017) and both the skull and mandible contain key apomorphies of both Atractosteus and Lepisosteus (Wiley 1976; Grande 2010). These include features like the orientation and contact of the dentary symphyses and the width of the skull, which differ among Lepisosteus and Atractosteus and might appear distinct in hybrid individuals. A total of 7 meristic counts and proportions were measured.

We included 24 morphometric landmarks: 12 on the skull in dorsal view and 12 on the mandible in ventral view (Figure S6). Landmarks were placed based on previously defined borders between major craniomandibular elements (Grande 2010). We digitized landmark coordinates and defined scales for each skull using tpsUtil64 v. 1.7 and tpsDIG2 v. 2.26. We ran analyses in both the R package geomorph (Adams and Otárola-Castillo 2013) and the program MorphoJ (Klingenberg 2011). In both programs, we applied a generalized Procrustes superimposition to exclude size, positional, and orientation effects before conducting principal components analyses on the data. We checked the resulting Procrustes fit graph for outliers and then ran principal components analyses.

Assessing the morphological disparity of gars in deep time.

To quantify variation in gar morphology over deep time, we collected data for individual specimens of extant and extinct gar species spanning over 75 million years to compare the shape of the skull across the gar crown clade. We opted for a quantitative comparison of these features in crown gars rather than an analysis of the rates of trait evolution through time given the small sample size of this clade (n=12 species; Fig. 1) and focused on two-dimensional measurements and meristic counts to mitigate the effects of post-mortem compression and deformation of gar fossils.

Our measurements and meristic dataset include all currently described extinct species confidently placed in the gar crown group (Fig. 1), as well as a new sample of Atractosteus spatula X Lepisosteus osseus hybrid crosses. Hybrids examined included one alcohol-preserved head (KUI uncatalogued), three alcohol-preserved full fish (KUI 18407, 18408, 18560), and three skeletonized skulls also used in our geometric morphometric analyses (KUI 18558, 18409, 18559). These specimens were collected from the same river systems from which our sequence data establishes the presence of hybrids: the Trinity River ~1.6 km north of the US Route 90 Bridge in Liberty, Texas (KUI 18407), near the mouth of the Trinity River (KUI 18560, 18558, 18559), and at the mouth of the Trinity River near Trinity Bay (KUI 18408, 18409).

Measurements taken included dimensional ratios of all gar skull roof bones, as well as standard length, head length, lateral line scale, and fin ray counts. All measurements were taken using digital calipers. The sample was combined with the total sample of crown gars in Grande (2010) for a total of n=124 specimens, and the dataset was analyzed and plotted in R using ggplot2. All measurements and meristic data are available in the Supplement.

Detection of introgression among living lineages of gars.

Lastly, we assessed whether living species of Lepisosteidae exhibit signatures of introgression over deeper phylogenetic time scales. First, we used MrBayes 3.2.7 (Ronquist et al. 2012) to generate Bayesian gene trees for the subset of 770 orthologous exons identified by Brownstein et al. (2023) to be variant in gars sampled for all extant gars and two teleost outgroups (Megalops cyprinoides and Osteoglossum bicirrhosum). These represent the variable exons in gars sampled to reconstruct our tip-dated Bayesian phylogenomic hypothesis. We used an HKY+I+G molecular evolutionary model as implemented in the computer program MrBayes v. 3.2 (Ronquist et al. 2012) and estimated posterior probabilities for the phylogeny and parameter values using Metropolis-couple Markov chain Monte Carlo (Larget and Simon 1999; Huelsenbeck et al. 2001). MrBayes was run for 1.5 x 10⁶ generations with two simultaneous runs each with four chains for each exon. Convergence of the posteriors was checked in Tracer v. 1.7 (Rambaut et al. 2018). Finally, we used the program DensiTree (Bouckaert 2010) to assess concordance among the resulting gene trees (Fig. 1, Fig. 4a).

Previous studies (e.g., Maddison, 1997; Mallet et al., 2016; Edelman et al. 2019) have noted that evidence of topological discordance across gene tree topologies may be reflective of both incomplete lineage sorting and introgression. As such, we conducted a secondary test of introgressive episodes across Lepisosteidae using PhyloNet 3.7.3 (Wen et al. 2018). Based on the 770 gene trees generated in Mr. Bayes 3.2.7, we inferred phylogenetic networks from maximum pseudolikelihood (MPL) with a reticulation of zero and then compared the log-probability of the reticulations up to six with 200 replications for each reticulation scenario. 

Comparison of hybridizing species pair MRCA ages and whole-clade rates.

We searched the literature for information on the most deeply divergent lineages that still hybridize for the major vertebrate clades analyzed in this study (Figure 2). The ages of the oldest hybridizing divergences were plotted against the mean exon substitution rate estimates generated in this paper. These data are included in the Supplemental File Table_S2.csv.