Gene flow and species delimitation in fishes of Western North America: Flannelmouth (Catostomus latipinnis) and Bluehead sucker (C. Pantosteus discobolus)

Abstract The delimitation of species boundaries, particularly those obscured by reticulation, is a critical step in contemporary biodiversity assessment. It is especially relevant for conservation and management of indigenous fishes in western North America, represented herein by two species with dissimilar life histories codistributed in the highly modified Colorado River (i.e., flannelmouth sucker, Catostomus latipinnis; bluehead sucker, C. (Pantosteus) discobolus). To quantify phylogenomic patterns and examine proposed taxonomic revisions, we first employed double‐digest restriction site‐associated DNA sequencing (ddRAD), yielding 39,755 unlinked SNPs across 139 samples. These were subsequently evaluated with multiple analytical approaches and by contrasting life history data. Three phylogenetic methods and a Bayesian assignment test highlighted similar phylogenomic patterns in each, but with considerable difference in presumed times of divergence. Three lineages were detected in bluehead sucker, supporting elevation of C. (P.) virescens to species status and recognizing C. (P.) discobolus yarrowi (Zuni bluehead sucker) as a discrete entity. Admixture in the latter necessitated a reevaluation of its contemporary and historic distributions, underscoring how biodiversity identification can be confounded by complex evolutionary histories. In addition, we defined three separate flannelmouth sucker lineages as ESUs (evolutionarily significant units), given limited phenotypic and genetic differentiation, contemporary isolation, and lack of concordance (per the genealogical concordance component of the phylogenetic species concept). Introgression was diagnosed in both species, with the Little Colorado and Virgin rivers in particular. Our diagnostic methods, and the agreement of our SNPs with previous morphological, enzymatic, and mitochondrial work, allowed us to partition complex evolutionary histories into requisite components, such as isolation versus secondary contact.


| INTRODUC TI ON
The delimitation of species (i.e., the process by which boundaries are not only identified but new species discovered; Wiens, 2007) is a fundamental issue in biology, and its mechanics contain aspects both theoretical and applied (Carstens, Pelletier, Reid, & Satler, 2013). It is a requirement not only for effective biodiversity conservation (Frankham, 2010) but also for management, particularly with regard to the Endangered Species Act (ESA) (Waples, 1991). However, distinct boundaries traditionally assumed to characterize species (Baum, 2009;Sloan, 2008) are particularly difficult to identify early in the speciation process (Sullivan et al., 2013) or within groups where extensive reticulation has occurred (Mallet, Besansky, & Hahn, 2015). This has led to an evolving interpretation of the speciation process, now viewed as a continuum from population through various ascending steps, but with genealogical distinctiveness achieved gradually and manifested differentially across the genome (Mallet, 2001).
Another complicating factor is the daunting number of narrow concepts that now encapsulate speciation (Coyne & Orr, 2004).
These species concepts tend to define species based on gaps that can often be in conflict with one another and fail to consider speciation as a "process" and thus contribute little to species delineation (de Queiroz, 2007). Clearly, a more unified approach that takes into account multiple lines of evidence would help in species delineation as well as defining conservation units across the population-species continuum, which in turn promotes effective management strategies to protect genetic diversity (as stipulated in Strategic Goal C of the 2010 United Nations Convention on Biological Diversity <https:// www.cbd.int/sp/targe ts/>).
The most commonly utilized approach for delineating lineages is DNA-based, but with questionable reliance upon a single marker (i.e., DNA barcoding; Ahrens et al., 2016), especially given the semipermeable boundaries now recognized in species. Despite being cofounded by various problems, single-gene methods still predominate in the literature, often with sample sizes that do not capture intraspecific haplotype variability, an issue scarcely parameterized (Phillips, Gillis, & Hanner, 2019). Additionally, single-locus delimitation methods fall under several broad categories, yet each suffers from limitations not easily overcome (Dellicour & Flot, 2018).
However, genomic DNA techniques are increasingly being applied to more formally delineate lineages (Allendorf, Hohenlohe, & Luikart, 2017). SNP (single nucleotide polymorphisms) panels are being applied to not only broaden and extend signals of population and species differentiation, but also to unravel their interrelationships. Two insights have emerged so far from the application of contemporary genomic approaches. First, increased resolution provided by SNPs has identified populations that diverge markedly within taxa and are subsequently characterized as "cryptic" species (Singhal, Hoskins, Couper, Potter, & Moritz, 2018;Spriggs et al., 2019), with population histories not only statistically inferred but also tested against alternative models of divergence. A second insight is that admixture among lineages is not only quite common but also greater than previously thought (Dasmahapatra, 2012;Fontaine et al., 2015;Quattrini et al., 2019), even to the extent of promoting new adaptive radiations (Lamichhaney et al., 2017).
The increased resolution provided by reduced representation genomic approaches has negative connotations as well. For example, elevated lineage resolution (as above) has the propensity to re-ignite earlier debates regarding the oversplitting of species (Isaac, Mallet, & Mace, 2004;Sullivan et al., 2014). Not surprisingly, this process is rife with value judgments, one of which seemingly intuits that species defined by parsing previously identified biodiversity are more problematic than those discovered de novo (Padial & de la Riva, 2006;Sullivan et al., 2014). Similar issues emerge when intraspecific diversity is interpreted for conservation actions (Funk, McKay, Hohenlohe, & Allendorf, 2012). For example, the evolutionarily significant unit (ESU) was conceived as a complement to existing taxonomy (Moritz, 1994;Ryder, 1986), with an intent to quickly identify conservation units worthy of protection without resorting to a laboriously slow and unwieldy taxonomic categorization. Population genomics provides an abundance of neutral loci for delimitation of ESUs. However, despite its intended simplicity, the ESU concept (Frazer & Bernatchez, 2001;Holycross & Douglas, 2007) has seemingly become an either/or categorization (i.e., dichotomized) such that it not only contravenes the continuum F I G U R E 1 Map depicting the Colorado River and Bonneville basins, with adjacent basins or recognized geographic regions; ID, Idaho; MEX, México; NM, New Mexico; NV, Nevada; TX, Texas; UT, Utah; WY, Wyoming through which populations evolve, but also reflects those difficulties that emerge when subspecies are designated arbitrarily from continuous geographic distributions (Crandall, Bininda-Emonds, Mace, & Wayne, 2000;Douglas, Douglas, Schuett, & Porras, 2006).
Similarly, the management unit (MU) is yet another conservation category with traditional roots, generally referred to in fisheries literature as a "stock" (Ryman & Utter, 1986). It now has a more contemporaneous meaning, defined primarily by population genomic data, and represents a conservation unit isolated demographically from other such units (Mussmann, Douglas, Chafin, & Douglas, 2019;Palsbøll, Bérubé, & Allendorf, 2014). While genomic data clearly hold great potential for elucidating the evolutionary process, arguments must still be resolved before they become a de facto diagnostic tool for species delineation (Stanton et al., 2019). For example, genomic techniques were unsuccessful in unraveling a hybrid complex among Darwin's finches in the Galapagos (Zink & Vázguez-Miranda, 2019).
In this study, we applied a large ddRAD dataset (~20,000 loci with ~100,000SNPs) to a contemporary framework for genomic analysis (Leaché & Fujita, 2010), wherein we (a) did Bayesian clustering to detect admixture in our study species so as to detect erroneous species designations derived from interspecific gene flow (Camargo, Morando, Avila, & Sites, 2012;Stewart, Timmer, Lawrence, Pryor, & Peever, 2014), (b) applied phylogenetic methods and hybrid testing to elucidate the clustering results, and then finally (c) test species delineation to examine potential splitting within the species examined.
This approach gains additional power when multiple lines of evidence are integrated, such as life history, geographic distributions, and morphology (Fujita, Leaché, Burbrink, McGuire, & Moritz, 2012;Knowles & Carstens, 2007;Schlick-Steiner et al., 2010); thus, we also tie in the previous published work on the species of concern.
As a result, the complex histories of study species can be more clearly discerned despite difficulties imposed by introgression. This is particularly appealing as herein, when problematic species are a focus of conservation concern (Pyron, Hsieh, Lemmon, Lemmon, & Hendry, 2016).

| The biogeography of our study species
Flannelmouth sucker (Catostomus latipinnis) and bluehead sucker (C. (Pantosteus) discobolus) have complex evolutionary histories that reflect historical introgression (Smith, Stewart, & Carpenter, 2013), as well as contemporary hybridization with various congeners (Bangs, Douglas, Thompson, & Douglas, 2017;Douglas & Douglas, 2010;Mandeville, Parchman, McDonald, & Buerkle, 2015). Both have been relatively understudied, yet their conservation concerns have accelerated due to a prolonged drought in western North America superimposed onto an ever-increasing anthropogenic demand for water (Seager & Vecchi, 2010). Given this, a federal and multistate effort has now coalesced on basin-wide mitigation and recovery of both species (Carman, 2007). Consequently, the accurate delimitation of species, as well as the designation of potential conservation units, is highly relevant, especially given that our study species comprise, in an historic sense, the greatest endemic fish abundance/biomass in the Upper Colorado River Basin (Hubbs and Miller, 1948).
Each species exhibits a different life history (Behn & Baxter, 2019;Sigler & Miller, 1963), with flannelmouth sucker primarily inhabiting the mainstem Marsh, 1998, 2003) and bluehead sucker preferring higher elevation streams that have subsequently become more fragmented over time (Hopken, Douglas, & Douglas, 2013). However, mark-recapture studies (Frazer et al., 2017)  The response of our study species to the geologic history of western North America is an aspect of their life histories (reviewed in Bezzerides & Bestgen, 2002). Vicariant processes (i.e., vulcanism and drainage rearrangements), coupled with episodic drought, have induced long periods of isolation sporadically augmented by more pluvial periods that, in turn, have promoted secondary contact (Smith, Badgley, Eiting, & Larson, 2010). Thus, a comparative study of each species can not only provide insights into the manner by which admixture has influenced their evolution, but also clarify our understanding of the Colorado River Basin itself.
Both species are primarily endemic to the Upper Colorado River Basin, but flannelmouth sucker is also found in the Virgin River of the Lower Colorado River Basin and bluehead sucker in the neighboring Bonneville Basin to the west (Figure 1). The latter may potentially represent a different species (C. (P.) virescens), as judged by morphological (Smith et al., 2013), mitochondrial (Hopken et al., 2013;Unmack et al., 2014), and nuclear phylogenies (Bangs, Douglas, Mussmann, & Douglas, 2018 in one major tributary-the Little Colorado River (Figure 2). One lineage may represent a unique species (i.e., Little Colorado River Sucker), currently grouped with flannelmouth sucker, whereas a second may be a unique subspecies (i.e., Zuni bluehead sucker, C. (P.) d.
yarrowi) now found only in the Zuni River (NM) and Kin Lee Chee Creek (AZ), but with a presumed historic distribution that potentially included the entire Little Colorado River (Minckley, 1973).
The quantification of molecular variability in both catostomids is a key element in delimiting species boundaries, management units, and historic patterns of reticulation. Here, we build upon previous work (Bangs, Douglas, et al., 2018) by applying species delimitation methods, phylogenomic (i.e., concatenated and multispecies coalescence), and population genomic approaches (i.e., Bayesian clustering and hybrid detection) to identify potential species range-wide, but with special focus on the Little Colorado River. In this regard, the impacts of divergent life histories, as well as the role of stream capture and hybridization, are particularly germane with regard to the breadth and depth of differentiation found within each.

| Data collection
DNA was extracted with PureGene ® Purification Kit or DNeasy ® Tissue Kit (Qiagen Inc.) and stored in DNA hydrating solution (same kits). Libraries for double-digest restriction site-associated DNA (ddRAD) were generated following Bangs, Douglas, et al. (2018).

| Filtering and alignment
Illumina reads were filtered and aligned per Bangs, Douglas, et al.

| Clustering algorithm and phylogenetic methods
All analyses were run separately by subgenus and utilized the unlinked SNP file generated from PyRAD, save for the concentrated SNP phylogenetic methods that required all the SNP files. replicates, averaged across replicates to determine final values.
The most likely clusters were resolved by using the estimated log probability of data Pr(x|k) and the Δk statistic (per Evanno, Regnaut, & Goudet, 2005). Bayesian clustering also substantiated that all contemporary hybrids with invasive White Sucker had been eliminated per a priori removal based on morphology and previous genetic work (Bangs, Douglas, et al., 2018;Douglas & Douglas, 2010 However, multispecies coalescent methods perform well in situations where introgression is limited, and thus represent an important consideration when species with admixed ancestries are delimited (Edwards et al., 2016).
Applicability of these methods is limited with regard to SNP data, due to the common requirement of a priori inference of gene trees (see Leaché et al., 2014b). Thus, multispecies coalescent inference was restricted to SVDquartets (Chifman & Kubatko, 2015) as implemented in PAUP* v. 4.0 (Swofford, 2003) that effectively bypasses the gene tree inference step, thereby extending its applicability to SNP datasets. This approach uses a coalescent model to test support for quartets and to calculate frequencies of SNPs for each species. The process does not require concatenation, but does necessitate the a priori partitioning of individuals into species or populations. Because of extensive run times using exhaustive tip sampling, species were instead subdivided into populations based on high support under both concatenated SNP methods. All possible quartets were exhaustively sampled using 1,000 bootstraps.

| Bayesian species delimitation
Species delimitation methods are a popular analytical approach, especially those coalescent-based (Fujita et al., 2012), and applicable to larger datasets. However, these can lead to oversplitting, particularly with respect to integrative taxonomic methods and Bayesian assignment tests (Miralles & Vences, 2013). The response of these methods to the effects of introgression is still tentative and thus should be viewed with caution (Leaché, Fujita, Minin, & Bouckaert, 2014).
Bayes factor delimitation (BFD; Leaché, Fujita, et al., 2014) is another powerful tool for testing proposed taxonomic revisions and to assess whether models are congruent with the patterns of divergence obtained from multilocus genetic data. We applied it to test  (Leaché, Fujita, et al., 2014).

| Hybrid detection
We calculated a hybrid index by mapping against interspecific heterozygosity. This served as a second means of assessing admixture, as well as to assess contemporary hybrid events. The est.h function (R-package

| Phylogeny
Both concatenated SNP methods produced the same topology for each species (Figures 3a and 4a clade in the Pantosteus subgenus (Figure 3b). For the concatenated methods, bluehead sucker was placed outside the remaining species (Figure 3a), whereas for the multispecies coalescent method, Rio Grande Sucker was outside (Figure 3b). The latter reflects the results of previous research, to include morphological phylogenies, fossil evidence (Smith et al., 2013), as well as mitochondrial (Chen & Mayden, 2012;Unmack et al., 2014), and nuclear phylogenies (Bangs, Douglas, et al., 2018). One impact of introgression was to obscure lineage-level topologies in Catostomus (Bangs, Douglas, et al., 2018).
Given this, we elected to emphasize full-lineage concatenated topologies, and note in so doing that topological discrepancies among the employed methods are both minimal, and reflective of processes that have been reviewed elsewhere. MRBAyeS, but less so by SVDquARtetS (<70% bootstrap support). Also, the split between Upper Colorado and Little Colorado rivers was only moderately supported (at 86%). It should be noted that Wenima was not included in the SVDquARtetS phylogenetic analysis, due to hybridization with sonora sucker. However, its removal had no effect on topology or supports.

| Structure
The optimum number of supported clusters for Pantosteus was k = 6 (per log probability and Evanno methods, Figure S1)

| Hybridization
Individuals (  the rest of the Little Colorado and Colorado rivers, were ranked lower than a model that split the Little Colorado River (to include Zuni River, Defiance Plateau, and Upper Little Colorado) from the Colorado River. However, the highest ranked model was one that split all three groups in the Little Colorado River (Table 2).

| D ISCUSS I ON
Contemporary hybridization is problematic for freshwater conservation and management, particularly with regard to invasions (Bangs, Oswald, et al., 2018;Hargrove, Weyl, Zhao, Peatman, & Austin, 2019) and translocations (Bruce & Wright, 2018). Yet, these situations can most often be resolved through proper application of genomic approaches. However, it is much more difficult in a deep history context, in that phylogenetic relationships may be obscured as a result. Similarly, introgression is difficult to detect given genetic recombination (Wallis et al., 2016). Interestingly, freshwater fishes show particularly high levels of hybridization, due in large part to the occurrence of numerous sympatric species with small population numbers that are subsequently fragmented by environmental perturbations (Dowling & Secor, 1997). These issues have clearly impacted western North American freshwater fishes and, in particular, the genera evaluated herein (Dowling et al., 2016;Mandeville, Parchman, Thompson, Compton, Gelwicks, Song, & Buerkle, 2017).
The six states that encompass the Colorado River Basin signed a "Range-wide Conservation Agreement Plan" (2004) to adaptively manage our two study species basin-wide [as well as a third species, Roundtail Chub (Gila robusta)]. This, in turn, was a pre-emptive mechanism for these states to avoid potential listing under the Endangered Species Act (Carman, 2007). All three species exhibit distinct life histories and habitat preferences that may have driven their potential divergences across the basin.
Since speciation is a gradual process with biodiversity elements scattered along its continuum (Sullivan et al., 2014), potential incongruence would be expected when different species delimitation methods are employed. Introgression would further complicate this process, yet its impacts on most species delimitation methods remain unknown; thus, confounding any attempt to decipher results (Camargo et al., 2012). As such, the guideline of Carstens et al. (2013) is important considerations in this process, that is, be conservative and employ multiple lines of evidence, given that a failure to delineate is expected. This includes the use of multiple algorithms for analyses of multilocus data and alternative lines of evidence that include (when possible) the life histories, morphologies, distributions, fossil histories, and behaviors of the biodiversity elements under study.
Here, we explore different species delimitation approaches for two species, Flannelmouth and Bluehead sucker, to include the recent listing of the endangered Zuni bluehead sucker under the ESA.
Our purpose was to evaluate similarities and differences in patterns of divergence in these two largely sympatric species with different life histories and to diagnose (if appropriate) the potential for taxonomic revisions. In doing so, we also examined the impacts of introgression as a mechanism to disentangle their complex evolutionary histories that have evolved in lockstep with the geomorphology of the basin.

| Life history and its effects on differentiation
Comparative phylogenomics of Catostomus and Pantosteus subgenera (per Smith et al., 2013) revealed parallel patterns throughout much of the Colorado River and neighboring basins (Bangs, Douglas, et al., 2018). However, the scale of divergence varied greatly between these groups, as emphasized within the Upper Colorado River Basin (Figure 1).
Although three distinct clades were identified in flannelmouth sucker, they are relatively recent as underscored by the lack of distinct clustering ( Figure 4) and having less than 1% fixed SNP sites as compared to >1.9% for all other comparisons (Table 3). These levels of differentiation fit with recent events, including volcanic barriers that appeared during in the last 20 kya, such as Grand Falls on the Little Colorado River (~20 kya).
Lineages of bluehead sucker, on the other hand, reflect temporally deeper origins as underscored by the distinct clustering, branch lengths (Figure 3), and number of fix SNPs (1.9%-3.3%; Table 3) that are equal to or greater than well-established species pairs represented in our analyses, as well as by previous mitochondrial dating We suggest the contrasting timescales for these clades may stem from life history differences, particularly with regard to subgeneric habitat preferences. Pantosteus is commonly designated as "mountain sucker," due to its predilection for cooler habitats within higher elevation streams, whereas Catostomus is physically larger, omnivorous, and restricted to larger rivers that form lower-elevation components of basins (Sigler & Miller, 1963;Smith, 1966). Although Bluehead and flannelmouth sucker largely co-occur, their habitat preferences are profound and must, in turn, influence diversification rates. For example, Douglas et al. (2003) (Hopken et al., 2013).

| Bonneville Basin
Although both species are sympatric in the Colorado River Basin, the bluehead sucker also occurs in the Bonneville and Upper Snake River basins (Figure 1). Therein, it may represent a unique species (originally described as C. (P.) virescens Cope & Yarrow, 1875;Snyder, 1924) that was subsequently collapsed into C. (P.) discobolus (Smith, 1966). The split between C. (P.) virescens in the Bonneville Basin/Snake River and C. (P.) discobolus in the Colorado River Basin is supported in all of our analyses. This includes population clustering, three different phylogenetic methods, and BFD analyses (Figure 3; Table 2). The convergence of all methods, along with recent morphological (Smith et al., 2013) and mitochondrial phylogenies (Hopken et al., 2013;Unmack et al., 2014), supports the reclassification of the Bonneville bluehead sucker. Furthermore, the chronology for the split between these two species (i.e., ~4.8 mya per mtDNA time-calibrated phylogenies) exceeds that found in other catostomid species (Unmack et al., 2014) and emphasizes the deep divergence.

| Little Colorado River Basin
Our phylogenetic analyses also separate the Little Colorado River

| Zuni bluehead sucker
When Pantosteus was first described (Cope & Yarrow, 1875), the Zuni bluehead sucker was designated as a separate species. Subsequent allozymic and morphological data (Smith, Hall, Koehn, & Innes, 1983) not only recalibrated it to subspecies, but also suggested a hybrid origin that encompassed Bluehead and Rio Grande sucker. However, results from our studies now refute this hypothesis by demonstrating alleles from Rio Grande Sucker are found only within a single population (i.e., Rio Nutria) (Figures 3c and 5d). This result is consistent with more contemporary analyses of allozymes (Crabtree & Buth, 1987) as well as single-gene sequencing data (Hopken et al., 2013;Turner & Wilson, 2009).
Zuni bluehead sucker seemingly originated in the mountains of northeast Arizona and northwest New Mexico, to include the Zuni River and Kin Lee Chee Creek of the Defiance Plateau (Smith et al., 1983). However, phylogenetic analyses render populations in Kin Lee Chee Creek and the Defiance Plateau as paraphyletic with the Zuni River and the remainder of the Little Colorado River (Figure 3a,b). In addition, the entire Little Colorado River Basin clade is a monophyletic group sister to the remainder of the Colorado River (Figure 3a,b). This suggests that Zuni bluehead sucker spread into the Little Colorado River following its integration with mountain streams (per Minckley, 1973;Smith et al., 1983). The current hypothesis (Smith et al., 1983) suggests that it was replaced by bluehead sucker in all Little Colorado River drainages, save Zuni River and Kin Lee Chee Creek.
However, population-clustering analyses (Figure 3c) yielded a clade unique to the Little Colorado River, within which only Zuni River populations were assigned. All other populations were assigned to a composite representing this cluster and the remainder of the Colorado River Basin, with proportions for the latter ranging from 0.5% to 38.6%. This admixture was also detected in hybrid index analyses, suggesting the remainder of the Little Colorado River Basin may represent an admixture of these two lineages (Figure 5c).

| Little Colorado River Sucker
In contrast to the Zuni bluehead sucker, the Little Colorado River Sucker did not cluster separately, despite its representation as a monophyletic group in all phylogenetic analyses (Figure 4). This may reflect its recent origin, concomitant with formation of Grand Falls ~20 kya. This vicariant break effectively separated the Upper Little Colorado River from the rest of the Colorado River and prevented contemporary upstream gene flow (Duffield et al., 2006).
Although similar contemporary phylogeographic patterns are found in Zuni bluehead sucker and Little Colorado River Sucker, different evolutionary histories are apparent, as driven by habitat preference.
This process ultimately resulted in levels of divergence that differ, but within similar contemporary ranges. This underscores the chaotic fluvial history of the Desert Southwest, as well as the need for comparative studies that can disentangle the organismal histories that coexist there.
Hybridization was also detected between Sonora and flannelmouth sucker in Wenima Wildlife Area of the Little Colorado River ( Figure 4c). These admixed individuals are presumably due to a recent hybrid event, as gauged by the variation in q-scores found in sonora sucker (Figure 4c), as well as hybrid index values (Figure 5b), high interspecies heterozygosity (Figure 5b), and the presence of four second-generation hybrids. Regardless, further sampling is needed to confirm this assumption.

| Virgin River
Despite forming a monophyletic group, the Little Colorado River Sucker fell within a paraphyletic flannelmouth sucker. This was due largely to the placement of the Virgin River (Figure 4), also suggested as potentially unique due to an elevated morphological variability stemming from potential hybridization with sonora sucker (C. insignis) and razorback sucker (Xyrauchen texanus) (Minckley, 1980).

Indeed, historic introgression with sonora sucker was detected in all
Virgin River samples, as reflected in the elevated hybrid index and low interspecies heterozygosity (Figure 5b). Although the sonora sucker proportion is reduced, it is nevertheless significant based on previous D-statistic tests (Bangs, Douglas, et al., 2018) and hybrid index values for all samples (Figure 5b).
Although the phylogenetic splitting of the three flannelmouth sucker groups (i.e., Upper Colorado, Little Colorado, and Virgin River) was also supported in BFD (Table 2), they grouped as a single cluster in StRuctuRe (Figure 4c) and the splits could not be replicated in cluster analyses, even at higher k-values. This, in turn, suggests a recent origin for these groups, further supported by their short branch lengths (Figure 4a). There is also a lack of fixed differences between these lineages in a previous mitochondrial analysis (Douglas et al., 2003). These considerations fit well with the previous assumption that the Virgin River population may have separated recently, that is, Late Pleistocene, most likely due to climatic oscillations that alternately connected and separated Grand Canyon and Virgin River as recently as 7.5 kya (Douglas et al., 2003). The support in BFD for the splitting of these groups may be due to an increased sensitively in defining recent splits, or may instead be biased by differential introgression with sonora sucker, particularly given the unknown capacity of this method to discern introgression (Leaché, Fujita, et al., 2014).

| CON CLUS IONS
Flannelmouth and Bluehead sucker are recognized as "species of concern" in the Colorado River Basin (Carman, 2007). Proposed taxonomic revisions will not only impact the management of these species, but also the basin as a biogeographic unit. Both species reflect similar phylogenomic patterns, yet their levels of divergence underscore evolutionary histories that differ significantly, and which impact their species delimitations. Three lineages of bluehead sucker were detected in all phylogenetic and population genetic methods, with C. (P.) virescens in the Bonneville and Upper Snake River elevated as a species separate from C. (P.) discobolus in the Colorado River (per Smith et al., 2013;Unmack et al., 2014). Results also support the Zuni bluehead sucker as a unique form. However, the current designation of Kin Lee Chee Creek as congruent with the Zuni River is erroneous, as they are instead paraphyletic. This situation can be resolved by including the Little Colorado River bluehead sucker, or by removing Kin Lee Chee Creek from the listing of the Zuni bluehead sucker under the ESA. The situation is further complicated by hybridization with Rio Grande Sucker and bluehead sucker from the Colorado River.
The Little Colorado River Sucker falls within a paraphyletic flannelmouth sucker and can only be resolved by designating the Virgin River population as a unique lineage. However, these three clades are of recent origin, based on population genetic analyses (herein) and the lack of resolution found in mitochondrial analyses (Douglas et al., 2003). Thus, all three flannelmouth sucker lineages are more accurately represented as evolutionary significant units (ESUs), as reflected by their reduced phenotypic and genetic differentiation.
They thus lack concordance under the genealogical component of the phylogenetic species concept.

ACK N OWLED G M ENTS
Numerous agencies contributed field expertise, specimens, technical assistance, collecting permits, funding, or comments for the completion of this study. We are in debt to students, postdoctorals, and faculty who contributed to the development of our research: S. ing through the revision process. We would also like to thank the reviewers and editors that provided guidance on the revision of this manuscript.

CO N FLI C T O F I NTE R E S T
None declared.