{'fraction_non_alphanumeric': 0.036688083796771236, 'fraction_numerical': 0.037012040386588195, 'mean_word_length': 5.142099154369093, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Taxonomic Uncertainty and Phylogenomics: Rescuing a Contentious Species Complex from the Anomaly Zone (6438)

Species are an indisputable unit for biodiversity conservation, yet their delimitation is fraught with both conceptual and methodological difficulties. A classic example is the taxonomic controversy surrounding the Gila robusta complex in the lower Colorado River of southwestern North America. Three species were originally defined according to subtle morphological differences with weak diagnostic power, with traditional genetic analyses failing to support the nominal species. Consequently, the complex was re-defined as a single polytypic unit, with the proposed 9threatened9 status of two putative species withdrawn at the federal level. Here, we utilized dense spatial and genomic sampling (N=387 and >22k loci) to re-evaluate the status of the complex, based on SNP-based coalescent and polymorphism-aware phylogenetic models. In doing so, all three species were supported as evolutionarily independent lineages, despite widespread phylogenetic discordance. To understand this discrepancy with past studies, we categorized evolutionary mechanisms driving discordance. We tested (and subsequently rejected) prior hypotheses suggesting that phylogenetic discord in the complex was hybridization-driven. Instead, we found the G. robusta complex to have diverged within the 9anomaly zone9 of tree space and, as such, have accumulated inconsistent patterns of diversity which have confounded prior studies. After extending these analyses with phylogeographic modeling, we propose that this is reflective of a rapid radiation promoted by Plio-Pleistocene tectonism. Our results not only support resurrection of the three species as distinct entities, but also offer an empirical example of how phylogenetic discordance can be categorized in other recalcitrant taxa.

Importantly, phylogenetic uncertainty also translates to taxonomic uncertainty. This is 30 because modern systematic taxonomy fundamentally describes homology [i.e. Darwin's (1859) 'propinquity of descent' (Simpson 1961)], which, by definition, requires a phylogenetic context. 32 Phylogenetic uncertainty in this sense can manifest itself as a soft polytomy (= 'honest' uncertainty), the erroneous promotion of non-monophyletic clades, or controversial 'splitting' 34 versus 'lumping.' Incomplete or biased sampling is often a driver of this disparity (Ahrens et al. 2016;Reddy et al. 2017). Here, narrow taxon sampling may introduce substantial ascertainment 36 bias (=systematic deviations due to sampling). On the other hand, a broader yet sparse sampling regime often fails to sample cryptic lineages (Heath et al. 2008) -with subsequent impacts on 38 both the delimitation of species (Pante et al. 2015;Linck et al. 2019) and study of their traits (Beaulieu and O'Meara 2018). 40 These sources of uncertainty culminate in topologies that often fluctuate with regard to sampling designs or methodologies, and this results in taxonomic uncertainty [e.g. Ctenophora 42 versus Porifera as sister to all other animals (Pisani et al. 2015;Whelan et al. 2015;Simion et al. 2017)]. Access to genome-scale data has alleviated some of these issues by offering a level of 44 precision not possible with single-gene phylogenies (Philippe et al. 2005). However, their inherent complexity and heterogeneity introduces new problems, and consequently, additional 46 sources of phylogenetic uncertainty.
Gene tree heterogeneity is a ubiquitous source of discordance in genomic data, and 48 "noise" as a source of this variance must consequently be partitioned from "signal" (where "noise" is broadly categorized as systematic or stochastic error). Large genomic datasets can 50 reduce stochastic error (Kumar et al. 2012), yet it still remains a prevalent issue when individual genes are examined (Springer and Gatesy 2016). On the other hand, systematic error in 52 phylogenomics may represent a probabilistic bias towards incongruence that is inherent to the evolutionary process itself (Maddison 1997). This, in turn, exemplifies the complications 54 introduced by genomic data: As genomic resolution increases, so also does the probability of sampling unmodeled processes (Rannala and Yang 2008;Lemmon and Lemmon 2013). This 56 potential (i.e., simultaneously decreasing stochastic error as systematic error increases) produces the very real possibility of building a highly supported tree that is ultimately incorrect. 58 Certain demographic histories are more predisposed to systematic error than others. For instance, when effective population sizes are large and speciation events exceptionally rapid, 60 time between divergence events may be insufficient to sort ancestral variation, such that the most probable gene topology will conflict with the underlying species branching pattern. This results 62 in what has been coined an "anomaly zone" of tree space (i.e., dominated by anomalous gene trees (AGTs); Degnan and Rosenberg 2006). Inferring species trees is demonstrably difficult in 64 this region (Liu and Edwards 2009), and exceedingly so if additional sources of phylogenetic discordance, such as hybridization, are also apparent (Bangs et al. 2018). 66 In clades with such complex histories, it is often unclear where the source of poor support and/or topological conflict resides. Yet, to analytically account for gene tree conflict, it is 68 necessary to categorize these sources and select approaches accordingly. Failure to do so promotes a false confidence in an erroneous topology, as driven by model misspecification 70 (Philippe et al. 2011). The overwhelmingly parametric nature of modern phylogenetics insures that imperative issues will revolve around the processes being modeled, and what they actually 72 allow us to ask from our data (Sullivan and Joyce 2005). However, the selection of methods that model processes of interest requires an a priori hypothesis that delimits which processes are 74 involved. Yet, diagnosing prominent processes is difficult in that a phylogenetic context is required from which to build hypotheses. Fortunately, a wealth of information can be parsed 76 from otherwise "non-phylogenetic" signal (sensu Philippe et al. 2005). For example, many statistical tests diagnose hybridization via its characteristic signature on the distribution of 78 discordant topologies (e.g. Pease and Hahn 2015). Theoretical predictions regarding AGTs and the parameters under which they are generated are also well characterized (Degnan and Salter 80 2005; Degnan and Rosenberg 2009). Thus, by applying appropriate analytical approaches that sample many independently segregating regions of the genome, empiricists can still derive 82 biologically meaningful phylogenies, despite the presence of complicated species-histories (McCormack et al. 2009;Kumar et al. 2012). 84 Here, we demonstrate an empirical approach that infers species-histories and sources of subtree discordance when conflict originates not only from anomaly zone divergences but also 86 hybridization. To do so, we used SNP-based coalescent and polymorphism-aware phylogenetic methods (Chifman and Kubatko 2014; Leache et al. 2014;De Maio et al. 2015) that bypass the 88 necessity of fully-resolved gene trees. We combine coalescent predictions, phylogenetic network inference (Solís-Lemus and Ané 2016), and novel coalescent phylogeographic methods (Oaks 90 2018) to diagnose the sources of phylogenetic discordance and, by so doing, resolve a rather conflicted complex of study-species. We then contextualize our results to demonstrate the 92 downstream implications of 'problematic' tree-space for threatened and endangered taxa, as represented by our study complex. 94 Gila 96 Few freshwater taxa have proven as problematic in recent years as the Gila robusta complex (Cyprinoidea: Leuciscidae) endemic to the Gila River basin of southwestern North 98 America (Fig. 1). The taxonomic debate surrounding this complex exemplifies an inherent conflict between the traditional rigidity of systematic taxonomy versus the urgency of decision-100 making for conservation and management (Forest et al. 2015). Our study system is the Gila River, a primary tributary of the lower basin Colorado River that drains the majority of Arizona 102 and ~11% of New Mexico. The critical shortage of water in this region (Sabo et al. 2010) is a major geopolitical driver for the taxonomic conflict surrounding our study species. As an 104 example, the lower Colorado basin is responsible for approximately half of the total municipal and agricultural water requirements of the state of Arizona, and nearly two-thirds of its total 106 gross state product (GSP) (Bureau of Reclamation 2012; James et al. 2014). This disproportionate regional reliance creates tension between the governance of a resource and its 108 usage (e.g. Huckleberry and Potts 2019) which in turn magnifies the stakes involved in conservation policy (Minckley 1979;Carlson and Muth 1989;Minckley et al. 2006). 110 We focused on three species (Roundtail chub, Gila robusta; Gila chub, G. intermedia; and Headwater chub, G. nigra) that comprise a substantial proportion of the endemic 112 ichthyfaunal diversity of the Gila Basin [=20% of 15 extant native species (excluding extirpated Gila elegans and Xyrauchen texanus); Minckley and Marsh 2009]. Historically, the focal taxa 114 have been subjected to numerous taxonomic rearrangements (Fig. 1). Until recently, the consensus was defined by Minckley and DeMarais (2000) on the basis of morphometric and 116 meristic characters. These have since proven of limited diagnostic capacity in the field, thus provoking numerous attempts to re-define morphological delimitations (Brandenburg et al. 2015;118 Moran et al. 2017;Carter et al. 2018). Genetic approaches have to date been unproductive (Schwemm 2006;Copus et al. 2018), leading to a recent taxonomic recommendation that 120 subsequently collapsed the complex into a single polytypic species (Page et al. 2016(Page et al. , 2017. For the sake of clarity, we employ herein the nomenclature of Minckley and DeMarais (2000). 122 Methods (1514)  124 Taxonomic sampling A representative panel of N=386 individuals (Table S1; Fig. 2) was chosen from existing 126 collections (Douglas et al. 2001;Douglas and Douglas 2007), to include broad geographic sampling of the complex as well as congeners. For the sake of clarity, we retained species-level 128 nomenclature for all members of the Gila robusta complex, and additionally discriminate between G. robusta from the upper and lower basins of the Colorado River ecosystem (Chafin et 130 al. 2019).
One species, Gila elegans, is extirpated from the wild, thus samples were provided by the 132 Southwestern Native Aquatic Resources and Recovery Center (Dexter, NM). The genus Ptychocheilus served to root the Gila clade within the broader context of western leuciscids 134 (Schönhuth et al. 2012(Schönhuth et al. , 2014(Schönhuth et al. , 2018.

Reduced-representation sequencing
Genomic DNA was extracted using either PureGene® or DNeasy® kits (Qiagen Inc.) and 138 quantified via fluorometer (Qubit™; Thermo-Fisher Scientific). Library preparations followed the published ddRAD protocol (Peterson et al. 2012). Restriction enzyme and size-selection 140 ranges were first screened using an in silico procedure (

Data processing and assembly
Raw Illumina reads were demultiplexed and filtered using the PYRAD pipeline (Eaton 156 2014). We removed reads containing >1 mismatch in the barcode sequence, or >5 low-quality base-calls (Phred Q<20). Homologs assembly was then performed using de novo clustering in 158 VSEARCH (Rognes et al. 2016)

Phylogenetic inference 168
We formulated two simple hypotheses with regards to independent evolutionary subunits. If populations represented a single polytypic species, then phylogenetic clustering should 170 reflect intraspecific processes (e.g. structured according to stream heirarchy; Meffe and Vrijenhoek 1988). However, if a priori taxon assignments are evolutionarily independent, then 172 they should be recapitulated in the phylogeny. Given well-known issues associated with application of supermatrix/ concatenation approaches (Degnan and Rosenberg 2006;Edwards et 174 al. 2016) and pervasive gene-tree uncertainty associated with short loci (Leaché and Oaks 2017), we also employed SNP-based methods that bypassed the derivation of gene trees (Leaché and 176 Oaks 2017).
We also used a polymorphism-aware method (POMO; Schrempf et al. 2016) in IQ-TREE 188 (Nguyen et al. 2014). POMO considers allele frequencies rather than single nucleotides, thus allowing evaluation of change due to both substitution and drift. To provide POMO with 190 empirical estimates of polymorphism, we used the entire alignment, to include non-variable sequences. We filtered liberally using individual occupancy thresholds of 10% per-locus so as to 192 maximize individual retention and per-population sample sizes. We then deleted populations that contained <2 individuals, and loci with >=90% missing data per-population. This yielded a 194 dataset of 281,613 nucleotides and 40 tips. Non-focal outgroups were excluded due to their disproportionate effect on missing data. 196 We also calculated concordance factors (CFs) using a Bayesian concordance analysis in BUCKY (Larget et al. 2010), parallelized across all quartets via an adaptation of the TICR 198 pipeline (Stenz et al. 2015). To prepare these data, we sampled all non-monomorphic full gene alignments for which at least 1 diploid genotype could be sampled per population. We excluded 200 outgroups and non-focal Gila so as to maximize number of loci retained. This yielded 3,449 genes across 31 sampled tips. Gene-tree priors were generated using MRBAYES v.3.2.6 (Ronquist 202 et al. 2012) with 4 independent chains, each of which was sampled every 10,000 iterations, with a total chain length of 100,000,000 iterations and 50% discarded as burn-in. BUCKY was then 204 run in parallel to generate quartet CFs across 31,465 quartets, using a chain length of 10,000,000, again with 50% burn-in. Quartet topologies were used to generate a population tree using 206 QUARTETMAXCUT (Snir and Rao 2012), using the get-pop-tree.pl script from TICR (Stenz et al. 2015; https://github.com/nstenz/TICR). 208

Comparing phylogenies and estimating site-wise conflict 210
To evaluate the performance of SVDQUARTETS, TICR, and PoMo, we first computed site-wise log-likelihood scores (SLS) for each topology by performing a constrained ML search 212 in IQ-TREE. For comparison, we also generated an unconstrained concatenated tree. All ML analyses employed a GTR model with empirical base frequencies and gamma-distributed rates, 214 and were assessed across 1,000 bootstrap pseudoreplicates. Analyses were also reduced to a subset of tips common across all variably filtered datasets. We quantified the phylogenetic signal 216 supporting each resolution as the difference in site-wise log-likelihood scores (ΔSLS) between each population tree and the concatenation tree (Shen et al. 2017). We then calculated site-wise 218 concordance factors (sCF) as an additional support metric (Minh et al. 2018).

Tests of co-divergence
The contemporary course of the Colorado River resulted from the Pliocene erosion of the 240 Grand Canyon (McKee et al. 1967;Minckley 1986). Gila in the lower Colorado River basin then differentiated following one or more colonization events (e.g. Rinne 1976). Subsequent work 242 (Douglas et al. 1999) supported this conclusion by examining contemporary phenotypic variation among all three species as a function of historical drainage connectivity, with the conclusion that 244 body shape was most readily explained by Pliocene hydrography.
We tested if divergences were best explained by a model of in situ diversification 246 following a single colonization event, or instead by multiple, successive colonizations. We compared divergence models using a Bayesian approach (program ECOEVOLITY, Oaks 2018) that 248 used a coalescent model (Bryant et al. 2012) to update a prior expectation for the number of evolutionary events across independent comparisons. Four independent MCMC chains were run 250 with recommended settings and a burn-in that maximized effective sample sizes. Event models followed a Dirichlet process, with the concentration parameter exploring four alternative gamma 252 distributed priors (i.e. α=2.0, β=5.70; α=0.5, β=8.7; α=1.0, β=0.45; and α=2.0, β=2.18).
We randomly sampled 2000 full-locus alignments, then examined potential co-254 divergences in the lower-basin complex by selecting a series of pairwise comparisons: Gila elegans x G. robusta (lower); G. seminuda x G. robusta (lower); G. jordani x G. robusta 256 (lower); G. intermedia x G. robusta (lower); and G. intermedia x G. nigra (lower). These targeted nodes represent H, G, F, I, and N in the SVDQUARTETS topology (Fig. 3A). 258 Results (1147)  260 Phylogenetic conflict in Gila We found tree reconstructions across all three population methods relatively congruent 262 (SVDQUARTETS = Fig. 3; TICR, and POMO = Fig. 4). The concatenated supermatrix tree ( Fig.   S1) was also largely congruent with the population trees, but with two major disparities 264 (discussed below). Bootstrap support was variable and declined with decreasing nodel depth in the SVDQUARTETS analysis (Fig. 3), whereas the vast majority of nodes in POMO were supported 266 at 100% (Fig. 4A).
All analyses consistently supported the monophyly of a clade consisting of G. intermedia, 268 G. nigra, and lower basin G. robusta (hereafter the 'lower basin complex'). This clade had high bootstrap support in both SVDQUARTETS and POMO, and was universally placed as sister to G. 270 jordani. Gila robusta was unequivocally polyphyletic in all analyses, forming two distinct clades separated by the Grand Canyon. The lower basin clade was monophyletic in all cases, save the 272 concatenated tree, where it was paraphyletic (Fig. S1). It was also consistently recovered as sister to a monophyletic G. nigra + G. intermedia, with the exclusion of a single sample site (Aravaipa 274 Creek) that nested within G. intermedia in the POMO tree. Of note, this population had been previously diagnosed as trending towards G. intermedia in terms of morphology (Rinne 1976;276 Demarais 1986), although hybridization was not supported by D-statistics (Table 1).
Topology within the G. nigra + G. intermedia clade was less consistent. Both were 278 reciprocally monophyletic in the SVDQUARTETS tree (albeit with low support; Fig. 3), whereas POMO yielded a monophyletic G. intermedia, with but one population (Spring Creek) contained 280 within G. nigra (Fig. 4A). The POMO tree also conflicted with the other methods in its paraphyletic placement of upper basin G. robusta. We suspect this represents an artefact of well-282 known hybridization with sympatric G. cypha (Dowling and DeMarais 1993;Gerber et al. 2001;Douglas and Douglas 2007;Chafin et al. 2019). 284

Discriminating among sources of phylogenetic conflict 286
Phylogenetic conflict was variably attributable to either hybridization or rapid divergence. We found support for a single reticulation event connecting G. seminuda and G. 288 elegans, an hypothesis consistent with prior interpretations (DeMarais et al. 1992). This  Table 1).
Multiple internode pairs were observed in the anomaly zone (Fig. 5). In all cases, 298 internode branches separating G. nigra and G. intermedia, and those separating their constituent lineages, reflected coalescent lengths that would yield anomalous gene trees. Not surprisingly, 300 the internode separating G. jordani from the lower basin complex, and that of G. robusta from G. intermedia/ G. nigra (Fig. 5C; tan branches) also fell within the anomaly zone, per TICR and 302 concatenated topology results.

Relative performance of species-tree methods
Change in site-likelihoods among constrained and unconstrained IQ-TREE searches in all 306 cases suggested that our recovered species-trees were supported by a minority of sites (Fig. S3), an observation consistent with tree regions being in the anomaly zone. Several discrepancies also 308 reflected idiosyncrasies among the different approaches. For example, the POMO topology has a paraphyletic upper basin G. robusta within which G. elegans, G. cypha, G. seminuda, G. jordani, 310 and the lower basin complex were subsumed (Fig. 4A). However, only ~10% of SNPs supported this resolution (Fig. S4), a value far below the theoretical minimum sCF derived from 312 completely random data (Minh et al. 2018). Of note, hybridization is a well-known artefact when a bifurcating tree is inferred from reticulated species (Sosef 1997;Schmidt-Lebuhn 2012), with 314 concatenation or binning approaches using genomic data being demonstrably vulnerable (Bangs et al. 2018). Thus, we tentatively attribute the observed paraphyly as an artefact of documented 316 hybridization between G. cypha and G. robusta (Chafin et al. 2019), and the inability of POMO to model hybridization. Hybridization also potentially drives the lack of monophyly in G. 318 seminuda, per TICR and the concatenation tree (Fig. S1).
We also explored the impact of matrix occupancy filters on SVDQUARTETS, and 320 bootstrap support and overall topological consistency declining with increasingly stringent filters (Fig. 3b). This corroborates prior evaluations with regard to the impacts of over-filtering 322 RADseq data (Eaton et al. 2017). In all cases, site-wide concordance was significantly predicted by subtending branch lengths, but not by node depths (Fig. S5). This suggests that site-wise 324 concordance was unbiased in our analyses at either shallower or deeper timescales, but was affected instead by the extent of time separating divergences. Some bioinformatic biases such as 326 ortholog misidentification or lineage-specific locus dropout will disproportionally affect deeper nodes (Eaton 2017). However, we interpret the lack of correlation between node depth and site-328 wise concordance as an indication that these processes lack substantial bias.

330
Biogeographic hypotheses and co-divergence ECOEVOLITY used four independent MCMC chains that sampled posterior parameter 332 space over 300,000 iterations. Burn-in was determined as the number of samples removed to optimize effective sample sizes (ESS>500 in all cases following 10-30% burn-in). Model 334 selection was not found to be vulnerable to alternative event priors (Fig. S6). The best-fitting model across all priors consistently demonstrated co-divergence of G. jordani with the lower 336 basin complex (G. robusta x G. intermedia and G. intermedia x G. nigra; Fig. 6). The divergence of G. elegans and G. seminuda from a theoretical lower basin ancestor pre-dates this putatively 338 rapid radiation, although it is unclear if these estimates were impacted by the aforementioned introgression between G. seminuda and G. elegans. 340 Posterior effective population size (N e ) estimates were large (e.g. >20,000) and consistent with previous estimates (Garrigan et al. 2002). Gila jordani was an exception, with a mean 342 posterior N e =6,062. We find this discrepancy rather unsurprising, given the extremely narrow endemism of this species (Tuttle and Scoppettone 1990), and its recent bottleneck (Hardy 1982), 344 although this is still a rather large estimate given the latter. Posterior divergence time estimates suggested a late-Miocene/ early-Pliocene origin of G. elegans. Results for G. seminuda and the 346 lower basin radiation indicated Pliocene and early Pleistocene divergences, respectively. These results are supported in the fossil record (Uyeno 1960;Uyeno and Miller 1963), although we 348 note paleontological evaluations of Gila have been sparse. Thus, we hesitate to interpret these as absolute dates, given our fixed mutation rate for these analyses and an uncertainty regarding the 350 capacity of RADseq methods to yield an unbiased sampling of genome-wide mutation rate variation (e.g. Cariou et al. 2016). 352

Discussion (1962 + 313 = 2275) 354
The goal of our study was to determine if extensive geographic and genomic sampling could resolve the taxonomically recalcitrant G. robusta complex. We applied diverse 356 phylogenetic models and tests of hybridization and predictions of parameter space within the anomaly zone to diagnose sources of subtree discordance. In so doing, we also tested multiple 358 hypothesized hybrid speciation events. We detected a single reticulation (G. seminuda), although other events with a lower component of genomic introgression may have also occurred. We 360 documented rapid co-divergence of lower basin taxa within the anomaly zone and were able to resolve these despite the prevalence of incomplete lineage sorting. This scenario (as outlined 362 below) is consistent with the geomorphology of the region and seemingly represents an adaptive radiation by our study complex, as facilitated by drainage evolution. 364

Methodological artefacts and conflicting phylogenetic hypotheses for Gila 366
Increased geographic and genomic sampling revealed the presence of diagnosable lineages within the G. robusta complex, with both rapid and reticulate divergences influencing 368 inter-locus conflict. Phylogenetic hypotheses for our focal group had previously been generated using allozymes (Dowling and DeMarais 1993), Sanger sequencing (Schwemm 2006;Schönhuth 370 et al. 2014), microsatellites (Dowling et al. 2015), and more recently RADseq (Copus et al.

2018). None could resolve relationships within the lower basin complex. To explain these 372
contrasts, we argue that prior studies suffered from systematic artefacts and ascertainment biases that were overcome, at least in part, by our approach. sampled extensively, including nearly all of the sites included in this study, but was only able to examine a handful of genes. Because anomalous gene trees are most probable under a scenario of 386 rapid radiation (as documented herein), the reduced number of loci used by Schwemm (2006) could not recover a consistent species tree. Copus et al. (2016Copus et al. ( , 2018) examined a dataset 388 containing 6,658 genomic SNP loci (across 1,292 RAD contigs), but only did so across a sparse sample of 19 individuals. A bioinformatic acquisition bias also likely impacted this study, in the 390 form of strict filtering that disproportionately excluded loci with higher mutation rates (Huang and Knowles 2016). 392 A necessary consideration when validating phylogenetic hypotheses across methods (and datasets) is to gauge compatibility between the underlying evolutionary processes and those 394 actually being modeled. In this sense, the consideration of statistical support metrics alone can be not only misleading, but also promote false conclusions. For example, bootstrapping is by far the 396 most prevalent method of evaluating support in phylogenetic datasets (Felsenstein 1985). While bootstrap concordances may be appropriate for moderately-sized sequence alignments (e.g. 398 Efron et al. 1996), they can be meaningless when applied to sufficiently large datasets (Gadagkar et al. 2005;Kumar et al. 2012). This is apparent in the high bootstrap support displayed for 400 anomalous relationships in our own analysis (Fig. S1). Phylogenetic signal also varies among loci, such that in many instances, relatively few loci drive contentious relationships (Shen et al. 402 2017). Likewise, not all methods are equal with respect to their simplifying assumptions. Given this, we feel it imperative to consider the biases and imperfections in both our data, and the 404 models we apply.

Complex evolution and biogeography of the Colorado River
The taxonomic instability in Gila is not uncommon for fishes of western North America, 408 where confusing patterns of diversity were generated by tectonism and vulcanism (Minckley et al. 1986;Spencer et al. 2008). This issue is particularly dramatized when viewed through the lens 410 of modern drainage connections (Douglas et al. 1999). Historic patterns of drainage isolation and intermittent fluvial connectivity not only support our genomic conclusions but also summarize 412 the paleohistory of the Colorado River over temporal and spatial scales.
The earliest record of fossil Gila from the ancestral Colorado River is mid-Miocene 414 (Uyeno and Miller 1963), with subsequent Pliocene fossils representing typical 'big river' morphologies now associated with G. elegans, G. cypha, and G. robusta (Uyeno and Miller 416 1965). The modern Grand Canyon region lacked any fluvial connection at the Miocene-Pliocene transition, due largely to regional tectonic uplifts that subsequently diverted the Colorado River 418 (Spencer et al. 2001;House et al. 2005). Flows initiated in early Pliocene (c.a. 4.9 mya; Sarna-Wojcicki et al. 2011), and subsequently formed a chain of downstream lakes associated with the 420 Bouse Formation (Lucchitta 1972;Spencer and Patchett 2002). Evidence suggests 'spillover' by a successive string of Bouse Basin paleolakes was episodic, and culminated in mid-Pliocene 422 (House et al. 2008), with an eventual marine connection via the Salton Trough to the Gulf of California (Dorsey et al. 2007). Prior to this, the Gila River also drained into the Gulf (Eberly 424 and Stanley 1978), and sedimentary evidence indicated that it was isolated from the Colorado until at least mid-Pliocene by a northward extension of the Gulf (Helenes and Carreno 2014). 426 This geomorphology is reflected in a broader phylogeographic pattern that underscores marked differences between resident fish communities in the upper and lower basins (Hubbs and Miller 428 1948).

Intra-basin diversification also occurred as an addendum to hydrologic evolution. 430
Although the course of the pluvial White River is now generally dry, it may have been a Pliocene-early Pleistocene tributary of a paleolake system when the proto-Colorado first 432 extended into the modern-day lower basin (Dickinson 2013). This may represent an initial colonization opportunity for upper basin fishes, an hypothesis that coincidentally aligns well 434 with our rudimentary age estimate for Virgin River chub, G. seminuda (Fig. 6). This early isolation, as well as the continued contrast between the spring-fed habitats therein, and the high 436 flows of the ancestral Colorado River, provide an explanation for the unique assemblage of Gila and other fishes therein (Hubbs and Miller 1948). with a relatively mesic Pliocene as precursor to a protracted monsoonal period extending through early Pleistocene (Thompson 1991;Smith et al. 1993). The latter, in turn, may have resulted in 444 relatively unstable drainage connections (Huckleberry 1996). The potential for climate-driven instability, and the complex history of intra-drainage integration of Gila tributaries during the 446 Plio-Pleistocene (Dickinson 2015), lends support to the 'cyclical-vicariance' model proposed by Douglas et al. (1999). These periods of isolation may have promoted an accumulation of 448 ecological divergences that persisted post-contact, and were sufficient to maintain species boundaries despite contemporary sympatric distributions and weak morphological 450 differentiation. This hypothesis is also supported by the non-random mating found among G. robusta and G. nigra, despite anthropogenically-induced contact (Marsh et al. 2017). 452 This study provides a much needed resolution to this debate by defining several aspects: First, our study reinforced the recognition of G. robusta as demonstrably polyphyletic, with two 464 discrete, allopatric clades corresponding to the upper and lower basins of the Colorado River (Dowling and DeMarais 1993;Schönhuth et al. 2014). These data, together with the geomorphic 466 history of the region that promoted endemic fish diversification (as above), clearly reject 'G.
robusta' as a descriptor of contemporary diversity. This underscores a major discrepancy in the 468 taxonomic recommendations for the lower basin complex (Page et al. 2016). Given that the type locality of G. robusta is in the upper basin (i.e., the Little Colorado River), we note a pressing 470 need either to determine taxonomic precedence for the lower basin 'G. robusta,' or to provide a novel designation. The potential resurrection of a synonym is a possibility, necessitating a 472 detailed examinations of the type specimens prior to a formal recommendation. This may be appropriately adjudicated by the AFS-ASIH Names of Fishes panel, as a follow-up to their 474 earlier involvement.
The situation with G. intermedia and G. nigra is slightly more ambiguous. The short 476 internodes and anomaly zone divergences identified herein explain previous patterns found in population-level studies, with elevated among-population divergence but scant signal uniting 478 species (Dowling et al. 2015). We also unequivocally rejected the previous hypothesis of hybrid speciation for G. nigra (Minckley and DeMarais 2000;Dowling et al. 2015). 480 Rather, intermediacy in the body shape of G. nigra reflects differences accumulated during historic isolation (Douglas et al. 1999) and/ or the retention of an adaptive 482 ecomorphology (Douglas and Matthews 1992). These hypotheses warrant further exploration, with provisional results employed in future management decisions (Forest et al. 2015). With 484 regards to taxonomy, we confidently recommend that G. intermedia be resurrected, and that additional studies be implemented to dissect the potential distinctiveness of G. nigra. For 486 management purposes, we echo a conservative, population-centric approach (previously argued for by Dowling et al. 2015;Marsh et al. 2017). Conclusion (313) The intractable phylogenetic relationships in Gila were resolved herein through improved 502 spatial and genomic sampling. Our data, coupled with polymorphism-aware methods and contemporary approaches that infer trees, yielded a revised taxonomy for Gila in the lower 504 Colorado basin. The geomorphic history of the Colorado River served to explain many anomalous patterns seen in this and previous studies, wherein opportunities for contact and 506 colonization were driven by the tectonism that is characteristic of the region. The signal of rapid diversification is quite clear in our data, as interpreted from patterns inherent to phylogenetic 508 discord. We emphasize that discordance in this sense does not necessarily represent measurement error or uncertainty, but rather an intrinsic component of phylogenetic variance 510 that is not only expected within genomes (Maddison 1997), but also a necessary component from which to build hypotheses regarding the underlying evolutionary process (Hahn and Nakhleh 512 2016). Ignoring this variance in pursuit of a 'resolved phylogeny' can lead to incorrect inferences driven by systematic error. Similarly, insufficient spatial or genomic sampling may also promote 514 a false confidence in anomalous relationships, particularly when character sampling is particularly dense, whereas taxon sampling is sparse. 516 We reiterate that phylogenetic hypotheses, by their very nature, cannot exhaustively capture the underlying evolutionary process. One approach is to categorize phylogenetic (and 518 "non-phylogenetic") signals in those regions of the tree that are refractive to certain models (as done herein). We acknowledge that attempting to reconstruct the past using contemporary 520 observations is a battle against uncertainty and bias, with the evolution of phylogenetic/ taxonomic hypotheses expected as additiona data are accrued. As such, we urge empiricists that 522 engage in taxonomic controversies (such as this one) to interrogate their results for transparency.
The task of sorting through conflicting recommendations invariably falls to managers, with 524 unreported biases (be they methodological or geopolitical) only confounding those efforts.   for divergence times are plotted on corresponding nodes. Units are in years, using a static 908 mutation rate of 1.2 e -08 substitutions per year. Posterior probabilities for divergence models (insert) suggest the co-divergence of Gila jordani, G. robusta, G. nigra, and G. intermedia. 910  and TICR phylogenies. Panels show (left to right): Linear regression of subtending branch lengths (log-transformed) with sCF; node height (cumulative branch lengths from root to focal 940 node); and densities of sCF across nodes as compared to the discordance factors for the two conflicting quartet resolutions (sDF1 and sDF2). 942 944 Figure S6: Prior and posterior probabilities for number of independent divergence events in ECOEVOLITY co-divergence models for Gila. Parameters across all runs were identical, except 946 for the shape (α) and scale (β) of the gamma-distributed prior on the Dirichlet process concentration. 948