A Complicated History of Mitogenome Introgression Among Luxilus Species (Teleostei, Family Leuciscidae) in the Ozark Highlands

1. Introduction

Hybridization and introgression occur frequently in nature and can be an important source of genetic variation in species (Arnold, 1997; Dowling & Secor, 1997; Edelman & Mallet, 2025; Hedrick, 2013; Mallet, 2005). Discordant patterns of mitonuclear introgression, which in the most extreme cases can result in mitogenome replacement even in the absence of measurable nuclear introgression, are well documented (Toews & Brelsford, 2012). Furthermore, the evidence of historical hybrid zones that have moved or disappeared is often apparent from genetic traces of introgression in the extant populations of impacted species (Wielstra, 2019).

Introgression may be driven by deterministic factors. Natural selection may favor introgression of locally adapted mitogenomes when species encounter novel or extreme environments or drive genetic rescue when mitogenomes have accumulated mutation load (Bonnet et al., 2017; Hill, 2019; Sloan et al., 2017). Mitochondrial replacement also may be nonadaptive and driven by sex-biased fitness of hybrid offspring (Dowling & Hoeh, 1991; Gerber et al., 2001; MacPherson et al., 2023). Alternatively, introgression may be driven as an outcome of range shifts and species invasion. When range shifts lead to invading species encountering congeneric species and engaging in hybridization, the outcome can lead to massive introgression in the genome of the invading species as a purely neutral demographic process (Currat et al., 2008; Wielstra & Arntzen, 2012).

Hybridization is common in fishes (Hubbs, 1955), and particularly prevalent in minnows and shiners of the family Leuciscidae (Scribner et al., 2001; Zbinden et al., 2023). Multi-species unguarded nesting behaviors are common among minnows and shiners and associated with high rates of hybridization (Corush et al., 2021). Discordant mitonuclear introgression patterns have been extensively reported among fish species, with examples of complete mitogenome replacement at the population level reported in charrs (Doiron et al., 2002; Wilson & Bernatchez, 1998), cichlids (Nevado et al., 2009), darters (Bossu & Near, 2009), pupfish (Carson & Dowling, 2006), topminnows (Duvernell & Schaefer, 2014), and shiners in the genus Luxilus (Duvernell & Aspinwall, 1995).

The genus Luxilus is comprised of nine North American shiner species (Dowling & Naylor, 1997). Studies have demonstrated the frequent occurrence of hybridization among sympatric Luxilus species (Dowling et al., 1989; Dowling & Moore, 1984; Meagher & Dowling, 1991), and members of the genus commonly spawn in heterospecific groups that can include multiple congeners (Pflieger, 1997). Within the zone of extant sympatry between L. chrysocephalus and L. cornutus, mechanisms of reproductive isolation and the dynamics of hybridization vary among populations isolated in separate drainages (Dowling et al., 1997; Dowling & Hoeh, 1991; Gleason & Berra, 1993). In the Ozark Interior Highlands, outside the zone of extant species sympatry, historical hybrid introgression has resulted in L. chrysocephalus populations in drainages throughout the region appearing monomorphic for L. cornutus mitogenomes (Duvernell & Aspinwall, 1995).

Recently, it was discovered that L. zonatus populations in some Ozark drainages exhibit L. chrysocephalus mitogenomes (Halas, 2011). What makes this particularly noteworthy is that in the river drainages where L. zonatus exhibit L. chrysocephalus mitogenomes, the L. chrysocephalus populations themselves are monomorphic for L. cornutus mitogenomes, suggesting a complex history of hybridization, introgression, and serial mitogenome replacement among Luxilus species in the Ozarks (Halas, 2011). The objective of our study was to determine the geographic extent of mitogenome introgression among L. zonatus populations, determine which drainages were impacted, and interpret these introgression events in a historical phylogeographic context.

2. Materials and Methods

The genus Luxilus is a relatively common minnow found throughout the eastern United States, including the cornutus and zonatus species groups. The cornutus group is comprised of four species, L. chrysocephalus, L. cornutus, L. cerasinus, and L. albeolus. Luxilus chrysocephalus is one of the most widely distributed members of this group with a range that spans the southern Great Lakes and the eastern two-thirds of the Mississippi River basin, as well as south into Gulf Coast drainages. Luxilus cornutus, has a similarly large range but is more northern in its distribution, including the northwestern portion of the Mississippi River basin, through the Great Lakes basin, and extending east to Atlantic drainages from Nova Scotia to the mid-Atlantic coast. A zone of extant sympatry between the species spans drainages of the Great Lakes, upper Illinois River, and portions of the upper Ohio River basin. The other two species of the group are more narrowly distributed in mid-Atlantic drainages with L. albeolus ranging from the Roanoke River of Virginia, southward to the Cape Fear River of North Carolina while L. cerasinus is found primarily in the Roanoke and New rivers.

There are three allopatric species in the L. zonatus species group, restricted mostly to the Ozark Highlands. They include L. cardinalis in the Arkansas River drainage, L. pilsbryi in the White River drainage, and L. zonatus in Missouri and Mississippi River tributaries that drain the Ozarks, including the Osage, Gasconade, and Meramec River drainages in the northern Ozark Highlands, and the Black and St. Francis River drainages in the southeastern Ozark Highlands (Mayden, 1988). All three species in the L. zonatus complex are sympatric with L. chrysocephalus.

We sampled L. zonatus from eleven sites in each of five major Ozark river drainages where they are found, along with L. pilsbryi from three sites in the White River drainage, and L. cardinalis from two sites in the Arkansas River drainage (Table 1). Specimens were collected by seining moderate to strong flowing current in the main channels of creeks and rivers. We also collected L. chrysocephalus from ten sites within the Ozark Highland study area (Table 1) and combined current with previously reported results (Duvernell & Aspinwall, 1995) for the species in this region. Additional reference sequences of Luxilus species outside the Ozark Highlands study area were included from previous published and unpublished sources (Dowling et al., 1992; Dowling & Naylor, 1997).

Fin clips were preserved in 95% ethanol, and voucher specimens were preserved in 10% formalin. All specimens were identified to species based on diagnostic morphological features (Pflieger, 1997). Ozark Highlands samples collected for this study were extracted using Blood and Tissue extraction columns (Qiagen). Subunit 2 of the NADH dehydrogenase gene (ND2) of mitochondrial DNA (mtDNA) was isolated by polymerase chain reaction (PCR) using primers designed from an L. chrysocephalus reference sequence. The primers were LuxND2_4010F 5’-GCCCATACCCCGAACATGAC-3’ and LuxND2_5145R 5’-TCTGCTTAGAGCTTTGAAGGCT-3’. PCR was conducted using Thermo Scientific PCR Master Mix (2X) (Thermo Fischer), and cycling conditions of an initial denaturing step at 95^oC for 4 min, 35 cycles of 95^oC for 30 sec, 58^oC for 30 sec, and 72^oC for 45 sec, followed by 72^oC for 4 min. PCR products were prepared for Sanger Sequencing using HighPrep PCR magnetic beads (Magbio Genomics). Sanger Sequencing of forward and reverse strands was conducted at the MU Genomic Technology Core (mugenomicscore.missouri.edu/) using the PCR primers. Chromatograms were assembled and aligned using Geneious Prime (Dotmatics). Additional reference samples of Luxilus species were prepared using similar PCR and Sanger sequencing methods as previously described (Dowling & Naylor, 1997).

A phylogeny was generated using the Maximum Likelihood (ML) method employed in MEGA11 (Tamura et al., 2021) with the GTR+I nucleotide substitution model and 1000 bootstrap replicates. Ericymba dorsalis, Notropis atherinoides, and N. buccata were included as outgroups. Nucleotide diversity (π), and nucleotide divergence (D_xy) parameters were estimated from ND2 sequence alignments using DnaSP 6.12.03 (Rozas et al., 2017).

Molecular divergence times of mitochondrial lineages were estimated using the suite of programs from BEAST 2.7.3 (Bouckaert et al., 2019) and reported as mean divergence time +/- 95% HPDs. Fossil-based calibration is confounded in Leuciscids by a poor fossil record. Therefore, we conservatively employed a uniform prior distribution of divergence rates as recommended by Smith and Larson (pers. comm.) to generate estimates of divergence time. Divergence rates for minnows and suckers correlate with size. Luxilus minnows (3g – 100g) lie at the low end of the size range for the family, with a corresponding expected rate of 2.5 – 3.5 times the pairwise divergence rate per million years. Input files were generated using BEAUti 2 and *BEAST was run with a species tree relaxed clock with a log normal rate distribution spanning the rate range reported by Smith and Larson. Analyses were run for 20 million generations with parameters logged every 1,000 generations. Multiple runs were conducted, and outputs were examined in Tracer 1.7.2 (Drummond & Rambaut, 2007) to verify that key parameters had effective sample sizes of 200 or greater. Log and tree files from four runs were combined using LogCombiner 2.7.3 (Drummond & Rambaut, 2007). with a 10% burnin. Mean age estimates and topology were based on the maximum clade credibility using common ancestor heights from our combined trees using TreeAnnotator 2.7.3 and the tree was visualized in FigTree 1.4.4.

3. Results

A mtDNA ML phylogeny (Figure 1) was consistent with previously reported species relationships in the genus Luxilus (Dowling & Naylor, 1997). Samples of the L. zonatus group (L. pilsbryi, L. cardinalis, and L. zonatus) formed a monophyletic clade (98% of bootstrap replicates). Within the L. zonatus species group, L. cardinalis and L. pilsbryi were most closely related (86% bootstrap replicates), followed by L. zonatus as expected (Dowling & Naylor, 1997; Mayden, 1988). However, among morphologically identified L. zonatus populations, the L. zonatus mtDNA haplogroup was only found in specimens from a portion of its range, in the Black/Current River drainage (Table 2, Figure 2a).

The cornutus group excluding L. cerasinus (ie. L. cornutus, L albeolus, and L. chrysocephalus) formed a second monophyletic clade (95% of bootstrap replicates). Luxilus cerasinus was sister to the rest of the cornutus group, though without bootstrap support. There were two mtDNA lineages within the cornutus group, one contained the L. cornutus mtDNA lineage including L. albeolus (100% of bootstrap replicates) while the other contained the L. chrysocephalus mtDNA lineage (100% of bootstrap replicates).

Comparison of mtDNA results with morphological species assignments identified considerable introgression among these taxa. A well-defined mtDNA lineage of morphologically identified L. zonatus (100% of bootstrap replicates) was nested within the range-wide L. chrysocephalus mtDNA lineage, indicating that this distinct lineage (hereafter ‘introgressed L. zonatus’) resulted from mitochondrial introgression from L. chrysocephalus into L. zonatus. These morphologically identified introgressed L. zonatus were from the St. Francis, Meramec, Gasconade, and Osage river drainages. Most individuals sampled from these populations exhibited the introgressed L. zonatus mtDNA haplogroup (Table 2). The paraphyletic L. chrysocephalus range-wide mtDNA clade, exclusive of introgressed L. zonatus, exhibited nucleotide diversity (π) of 0.033. Nucleotide diversity for the monophyletic introgressed L. zonatus lineage was nearly an order of magnitude lower at 0.004.

The L. cornutus mtDNA haplogroup also exhibited evidence of introgression as previously reported (Duvernell & Aspinwall, 1995). Sampled populations of morphologically identified L. chrysocephalus from drainages throughout the northern Ozark Highlands, as well as the White River drainage in the southern Ozark Highlands all exhibited a monophyletic clade of introgressed L. cornutus mtDNA (98% support, known hereafter as ‘introgressed L. chrysocephalus’) (Table 2, Figure 2b) and were the sister lineage to a range-wide L. cornutus clade. The range-wide L. cornutus and introgressed L. chrysocephalus clades appeared reciprocally monophyletic, though with lower bootstrap support (86%) and small sample size for the L. cornutus clade exclusive of the introgressed L. chrysocephalus clade. The L. cornutus clade exhibited nucleotide diversity of 0.010 compared to 0.006 for the introgressed L. chrysocephalus clade. Samples of morphologically identified L. chrysocephalus taken from the St. Francis River and Whitewater River, exhibited range-wide L. chrysocephalus haplotypes while the Little St. Francis River sample included both range-wide L. chrysocephalus and introgressed L. zonatus haplotypes (Table 2, Figure 2b).

Among all sampled L. zonatus individuals, the only ones that did not exhibit the L. zonatus haplotype or the here-defined introgressed L. zonatus haplotype were one individual in the Little Piney, a tributary of the Gasconade River, which exhibited an introgressed L. chrysocephalus haplotype, and one L. zonatus individual from the Little St. Francis River exhibiting an L. chrysocephalus haplotype (Table 2, Figure 2a).

A calibrated phylogeny (Figure 3) exhibited very large HPD confidence intervals in part due to calibration uncertainty built into the analysis. Nevertheless, confidence intervals placed the coalescent ages of most clades of interest within the Pleistocene (2.6 to 0.01 Ma). Our analysis placed the most likely age of the most recent shared ancestor between L. cornutus and L. chrysocephalus as pre-Pleistocene, at 5.05 Ma (0.61 - 11.77 Ma), and the L. zonatus group (L. zonatus, L. pilsbryi, L. cardinalis) at 4.75 Ma (0.59-11.07) (Figure 3). The coalescent age of all samples in the L. chrysocephalus mtDNA lineage was 1.98 Ma (0.24 – 4.62 Ma), and the introgressed L. zonatus clade within the L. chrysocephalus mtDNA lineage was 0.30 Ma (0.03 – 0.72 Ma). Nucleotide divergence (D_xy) between the introgressed L. zonatus clade and the next most closely related L. chrysocephalus haplotype was 0.049 with an estimated age of 1.52 Ma (0.16-3.57 Ma). The age of the ancestor of the L. cornutus and introgressed L. chrysocephalus mtDNA lineage (i.e. excluding L. albeolus) was 0.73 Ma (0.07 – 1.72 Ma), with an age for the nested range-wide L. cornutus clade estimated at 0.42 Ma (0.04 –1.02 Ma), and the introgressed L. chrysocephalus clade estimated at 0.31 Ma (0.03 – 0.74 Ma), respectively.

4. Discussion

Discordant mitonuclear introgression has been widely documented in many organisms (Bonnet et al., 2017; Hill, 2019; Sloan et al., 2017; Toews & Brelsford, 2012). The features that make our study exceptional are first that mitogenome replacement appears to have occurred twice in succession involving three congeneric shiner species. Second, the ranges of impacted populations in the two affected species are largely congruent and cover multiple isolated drainages that span the northern and eastern Ozark Highlands. The phylogeographic patterns of mitogenome introgression of both L. chrysocephalus and L. zonatus require at least two separate incidences of mitonuclear introgression that are Pleistocene in age, leading to mitogenome replacement in ancestral populations of both species. The broad, congruent geographic overlap of impacted populations, involving numerous Ozark Highland drainages, may reflect the impact of shifting distributions of both species associated with Pleistocene glacial cycles, and coincident with Ozark Highlands geomorphological processes.

Glaciation in North America west of the present day Mississippi River valley was marked by the Laurentide Ice Sheet advancing as far south as the vicinity of the Missouri River Valley at the northern edge of the unglaciated Ozark Highlands. Five glacial advances occurred in this region prior to the Illinoian stage at approximately 2.4, 1.3, 0.8, 0.4, and 0.2 Ma (Rovey & Balco, 2011). These events modified the course of rivers in the region, pushing the range of affected species southward (Robison, 1986). Periglacial lakes along the glacial front would have increased in size as glaciers melted. This would have also increased connection of periglacial lakes along the glacial front, and overflow of meltwater may have allowed for east-west dispersal of fishes along the northern border region of the Ozark Highlands (Robison, 1986).

The major drainages of the Ozark Highlands radiate out of the Salem Plateau. Biogeographically, drainages north of the Salem Plateau (e.g. Osage, Gasconade, Meramec) have similar ichthyofaunal species assemblages that differ in notable respects from drainages south of the Salem Plateau (e.g. St. Francis, Black, White). These biogeographic differences are evident in some Leuciscids (Campostoma spp., Cyprinella spp., Notropis spp.), Centrarchids (Ambloplites spp., Lepomis spp.), Ictalurids (Noturus spp.) (Pflieger, 1997), and intraspecific phylogeography of the Etheostoma caeruleum complex (Ray et al., 2006). There is also evidence that connections occurred occasionally between drainages that were otherwise isolated by the Salem Plateau. A population genetic study of the topminnow, Fundulus olivaceus, identified gene flow between the Big and St. Francis rivers, the Meramec and Black rivers, the Big Piney and North Fork rivers, and the Meramec and Gasconade rivers (Duvernell et al., 2019). Study of fluvial geomorphology in the Ozark Highlands offers evidence that headwaters on opposite sides of the Salem Plateau have commingled historically through numerous stream capture events (Beeson et al., 2017) creating opportunities for faunal exchanges among drainages.

The L. zonatus species group (L. cardinalis, L. pilsbryi, and L. zonatus) diversified in isolated drainages, likely beginning in the Pliocene (common ancestor estimated here at 4.7 Ma), and the phylogeographic distribution of these three species is consistent with their origin in the southern Ozark Highlands followed by a range expansion of L. zonatus into the more northern drainages within its range. Those northern L. zonatus populations appear monophyletic for the introgressed mitogenome acquired from L. chrysocephalus. Among L. zonatus populations, mitogenome introgression has impacted at least four of the five major Ozark drainages encompassed by the extant species distribution such that, range-wide, far more L. zonatus individuals exhibit introgressed than non-introgressed mitogenomes. The distribution of mitogenome replacement in northern L. zonatus populations is consistent with massive introgression predicted by an invasion-with-hybridization model driven by demographic inequality at a northward invasion front (Currat et al., 2008). According to this prediction, as L. zonatus expanded its range into eastern and northern portions of the Ozark Highlands where L. chrysocephalus was already present, the demographic inequality alone between resident and invading populations could have been sufficient to drive asymmetric mtDNA introgression (Currat et al., 2008; Wielstra & Arntzen, 2012). The reduced nucleotide diversity, and monophyly of the introgressed L. zonatus clade relative to other range-wide L. chrysocephalus haplotypes may indicate that this mitogenome replacement event coincided with a founder event initially within a small geographic region (e.g. within one drainage) and then expanded into other drainages through continued L. zonatus range expansion. The broad extant syntopy of L. zonatus and L. chrysocephalus throughout the Ozark Highlands suggests a northward expansion of L. zonatus without competition or species replacement.

The age of the split between L. cornutus and L. chrysocephalus is similarly dated to the Pliocene (5.05 Ma). With a northern L. cornutus range that was almost completely formerly glaciated, the presence of introgressed L. cornutus mitogenomes in Ozark L. chrysocephalus populations provides evidence for range overlap within or in the vicinity of the Ozark Highlands at least sometime during the Pleistocene (Duvernell & Aspinwall, 1995). The divergence between introgressed haplotypes found in Ozark L. chrysocephalus and the range-wide L. cornutus species haplogroup is dated to approximately mid-late Pleistocene (0.73 Ma). Evidence that hybrid introgression resulted in L. chrysocephalus mitogenome replacement from an L. cornutus donor spans all drainages originating in the Ozark Highlands, with most drainages appearing monophyletic for the introgressed haplogroup, and only drainages in the southeastern Ozark Highlands (Black, St. Francis, and Whitewater Rivers) exhibiting predominantly non-introgressed L. chrysocephalus mitogenomes (Duvernell & Aspinwall, 1995).

Our phylogenetic analysis shows that, after introgression, range-wide L. cornutus and introgressed Ozark L. chrysocephalus populations achieved reciprocal monophyly. The similar levels of nucleotide diversity between range-wide L. cornutus and introgressed Ozark L. chrysocephalus suggests that an L. cornutus range-wide expansion was similarly Pleistocene in age and occurred on a similar time frame to when the Ozark L. chrysocephalus introgression event occurred. The relatively limited range-wide nucleotide diversity in the northerly distributed L. cornutus compared to the more southerly distributed L. chrysocephalus is consistent with contrasting and disproportionate impacts of Pleistocene glaciation on these species in their respective extant latitudinal distributions.

The notable feature of these overlapping patterns of introgression is that in multiple drainages where L. zonatus appears monomorphic for an introgressed haplogroup from an L. chrysocephalus donor, sympatric L. chrysocephalus populations appear monomorphic for a different introgressed haplogroup from an L. cornutus donor. This dual mitogenome replacement pattern was observed in all drainages north of the Salem Plateau, including all samples in the Osage, Gasconade and Meramec River drainages. The coalescent ages of the introgressed L. chrysocephalus and L. zonatus haplogroup clades are both pre-Illinoian late Pleistocene, and with essentially the same age estimates (310 and 297 kya, respectively), suggesting that introgressed populations of both species may have established nearly concurrently via either range expansions or gene flow. However, L. zonatus introgression would have had to preceded L. chrysocephalus introgression or else introgressed L. zonatus populations would exhibit L. cornutus-lineage haplotypes instead of L. chrysocephalus-lineage haplotypes.

The phylogeographic distribution of mitogenome replacement events across multiple isolated drainages may suggest a deterministic driver operating similarly in separate drainages throughout the Ozark Highlands. However, the concurrent prevalence of introgressed mitogenomes in populations of both L. chrysocephalus and L. zonatus presents challenges to the invocation of natural selection as a driving factor. If mitogenome replacement was prevalent across multiple Ozark drainages because natural selection favored introgression of locally adapted mitogenomes in both species, then this does not account for why natural selection would favor different non-native mitogenomes in sympatric populations of each species. Similarly, if the driver of selection favoring introgressed mitogenomes was to alleviate some form of accumulated mutation load, as might have been caused by founder effects resulting from shifting ranges driven by glacial cycles, then that would not address the multi-drainage pattern of introgression or explain why introgressed mitogenomes were more fit than native mitogenomes in both species. The alternative explanation for the broad distribution of introgressed mitogenomes in both Luxilus species across very similar geographic areas that encompass most of the same Ozark drainages is range expansions and drainage invasions that could have occurred repeatedly for both species. The shared history of mitogenome haplogroups points to a complicated history of hybridization and introgression between these species. This calls for a better understanding of the phylogeographic histories of both species that will only be possible with further sampling and nuclear genome data.

Distributions of mitochondrial haplogroups between L. zonatus and L. chrysocephalus are further complicated in the St. Francis River drainage in the southeastern Ozark Highlands. The L. chrysocephalus sampled in the St. Francis River drainage in this study exhibited exceptional mitochondrial sequence diversity that reflected the entire range-wide mitogenome diversity of L. chrysocephalus (Figure 1). Results reported here are consistent with those reported by Duvernell and Aspinwall (Duvernell & Aspinwall, 1995) except that the earlier study relied on a PCR-RFLP assay that only resolved L. cornutus and L. chrysocephalus mtDNA lineages, and did not delineate differences between the introgressed L. zonatus haplogroup and all other range-wide L. chrysocephalus haplotypes.

In the Little St. Francis River, both L. zonatus and L. chrysocephalus were polymorphic for L. chrysocephalus haplotypes, including a predominance of the introgressed L. zonatus haplogroup in both species. The prevalence of the introgressed L. zonatus haplogroup could be indicative of an ancestral polymorphism in the Little St. Francis L. chrysocephalus population that has been shared with L. zonatus through introgressive hybridization. The St. Francis River L. chrysocephalus could be the source of the introgressed L. zonatus haplogroup that subsequently expanded into other Ozark Highland drainages. Alternatively, the shared mitochondrial polymorphisms between L. zonatus and L. chrysocephalus in the Little St. Francis River could indicate more recent, ongoing hybridization between the species. Further sampling within this drainage, with the inclusion of nuclear genetic data, is needed before conclusions can be drawn. Either way, this level of interspecies haplogroup sharing is distinctive within this system and was not encountered to the same extent in any other drainages. In the Osage, Gasconade and Meramec drainages where L. chrysocephalus and L. zonatus are syntopic and both essentially monomorphic for introgressed haplogroups, there is little evidence of recent or ongoing hybridization between the species except for one L. zonatus individual from Little Piney Creek in the Gasconade River drainage that exhibited the introgressed L. chrysocephalus mitogenome.

5. Conclusions

The concurrent patterns of mitogenome introgression among Luxilus species within the Ozark Highlands offer an exceptional example of mitonuclear discordance impacting multiple species and provide unique insights into Pleistocene-era impacts of gene flow and range shifts driven or facilitated by glacial advances and geomorphological processes occurring in the Ozark Highlands during this period.