Taxonomic Uncertainty and the Anomaly Zone: Phylogenomics Disentangle a Rapid Radiation to Resolve Contentious Species (Gila robusta Complex) in the Colorado River

Abstract

Species are indisputable units 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. Nominal species designations were originally defined according to weakly diagnostic morphological differences, but these conflicted with subsequent genetic analyses. Given this ambiguity, the complex was re-defined as a single polytypic unit, with the proposed "threatened" status under the U.S. Endangered Species Act of two elements being withdrawn. Here we re-evaluated the status of the complex by utilizing dense spatial and genomic sampling (n = 387 and >22 k loci), coupled with SNP-based coalescent and polymorphism-aware phylogenetic models. In doing so, we found that all three species were indeed supported as evolutionarily independent lineages, despite widespread phylogenetic discordance. To juxtapose this discrepancy with previous studies, we first categorized those evolutionary mechanisms driving discordance, then tested (and subsequently rejected) prior hypotheses which argued phylogenetic discord in the complex was driven by the hybrid origin of Gila nigra. The inconsistent patterns of diversity we found within G. robusta were instead associated with rapid Plio-Pleistocene drainage evolution, with subsequent divergence within the "anomaly zone" of tree space producing ambiguities that served to confound prior studies. Our results not only support the resurrection of the three species as distinct entities but also offer an empirical example of how phylogenetic discordance can be categorized within other recalcitrant taxa, particularly when variation is primarily partitioned at the species level.

Keywords: RADseq; anomaly zone; conservation genomics; hybridization; phylogenomics; taxonomic uncertainty.

Fig. 1.

Timeline of the conservation status of Gila species endemic to the lower Colorado River basin [*See Copus et al. (2018) for a detailed overview of taxonomic synonymies; †“The Center” refers to the Center for Biological Diversity (501c3), Tuscon, AZ; ‡“DPS” = Distinct Population Segment as referenced in the United States Endangered Species Act (ESA 1973; 16 U.S.C. § 1531 et seq), here referring specifically to a lower basin sub-unit of Gila robusta]. Note that the timeline is not to scale.

Fig. 2.

Sampling localities for Gila (n = 380 individuals) within the Colorado River Basin, southwestern North America. Locality codes are defined in Supplementary table S1, Supplementary Material online. Sympatric locations (R14 and C2) are slightly offset for visibility purposes. Map insert increases the viewing scale for sampling sites within the lower basin G. robusta “complex” (Bill Williams and Gila rivers).

Fig. 3.

(A) Majority-rule consensus cladogram of SVDQuartets across 12 variably filtered SNP data sets varying from 7,357–21,007 SNPs and 256–347 individuals representing 16 Gila OTUs from the Colorado River Basin; Ptychocheilus spp. used as outgroups. (B) Binned bootstrap concordance values are reported for major nodes in the majority-rule consensus tree (A) labeled as (A–P). Supports are partitioned by data set, coded by the matrix occupancy threshold per individual (“i”) and per column (“c”; e.g., i50_c50 = 50% occupancy required per individual and per column). Dashed terminal branches indicate positions for taxa missing from >50% of data sets. For detailed locality information, refer to supplementary table S1, Supplementary Material online.

Fig. 4.

(A) PoMo phylogram across 12 Gila OTUs from the Colorado River basin. Branch lengths reflect the number of substitutions and inferred number of drift events per site. Branch supports (only shown for those <100%) represent concordance among 1,000 bootstrap replicates, inferred using a data set consisting of 281,613 nucleotides and 40 tips (i.e., populations); (B) Corresponding TICR phylogram reporting branch lengths in coalescent units, calculated from 31,465 quartets evaluated across 3,449 full alignments of ddRAD loci. For detailed locality information, refer to supplementary table S1, Supplementary Material online.

Fig. 5.

Diagram comparing internode pairs within the anomaly zone, as determined using coalescent-unit transformed branch lengths mapped onto the (A) SVDQuartets, (B) PoMo, (C) TICR, and (D) concatenated trees (displayed here as cladograms). Paired internodes are color-coded: two successive internode branches of the same color = a pair of coalescent branch lengths falling within the anomaly zone; branches bicolored = branches involved in two separate significant anomalous divergences.

Fig. 6.

Admixture results for lower basin Gila robusta, G. intermedia, and G. nigra at K = 10–12. Individuals were obtained from 21 sites in the Gila River (fig. 2B). Results were generated for n = 140 individuals having <50% missing data for unlinked SNPs with a minor allele frequency ≥0.05 (=5,118 SNPs). Individuals are arranged according to the SVDQuartets individual-level phylogeny and are represented by stacked bar plots where colors are proportional to assignment probabilities as aggregated by Clumpak for 20 replicate runs of Admixture.

Fig. 7.

Posterior probability distributions for node ages and effective populations sizes (N_e) for six Gila OTUs from the Colorado River Basin. Estimates were derived from EcoEvolity and 2,000 randomly sampled full-length ddRAD locus alignments. Branches are annotated with a mean (std. dev.) N_e and posterior probabilities for divergence times are plotted on corresponding nodes. Units are in years, using a static mutation rate of 1.2 e⁻⁰⁸ substitutions per year. The inset figure shows posterior probabilities for the total number of divergence events, as contrasted with a prior distribution weighted against codivergence (i.e., with all 5 nodes having different ages).