Geopolitical species revisited: genomic and morphological data indicate that the roundtail chub Gila robusta species complex (Teleostei, Cyprinidae) is a single species

G. robusta Baird & Girard 1853a

Gila robusta Baird & Girard 1853a is the type species for the genus. The type series of the species is cataloged as USNM 246, but a note included with the type specimens states, “These, the types of Gila robusta B.+G., are cat. No. 246. They were reentered by error as 47,983 and attributed to nos. 276 + 273, which are cods! Nos. 276 + 277 were attributed as type nos. of this species, by error, by Jordan & Evermann (1896: 227). R.R. Miller III: 1945” (Fig. S1). This is one of a number of clerical errors uncovered in the taxonomic history of this fish and is indicative of the historic pattern of confusion surrounding the systematics of this complex.

The original description reported that the collection locality for the syntypes was the Zuni River, New Mexico, but Smith, Miller & Sable (1979) suggest that this locality is a clerical error, based on the argument that the Zuni River was unsuitable habitat for G. robusta during the time at which the type specimens were collected (1851). Smith, Miller & Sable (1979) suggested that the specimens were actually from another collection site of the Sitgreaves expedition, the Little Colorado River, below Grand Falls, Coconino County, Arizona. Sublette, Hatch & Sublette (1990) dispute the assertions of Smith, Miller & Sable (1979), and contend that the syntypes were collected from the Zuni River in 1851, and note that additional specimens were subsequently collected twice in 1873 and once in 1879 on the Zuni River by different collectors. It is highly unlikely that multiple clerical errors on different expeditions would have occurred at this locality with this species. Sublette, Hatch & Sublette (1990) suggest that the Zuni River represented marginal habitat that may have received recruits from the Rio Pescado and that they have since been extirpated from the Zuni River. The journal of SW Woodhouse, the naturalist who originally collected these fishes, clearly states that on Saturday the 6th of September 1851, he received these fishes while camped at the Zuni Pueblo, which were “collected from the creek.” The argument that the syntypes were collected on the Little Colorado River rather than the Zuni is therefore dubious and we maintain the type locality as the Zuni River for G. robusta.

G. intermedia (Girard 1856)

Girard (1856) described Tigoma intermedia from specimens collected on the Rio San Pedro in Arizona and noted that it was morphologically intermediate between T. pulchella (Baird & Girard 1854) and T. purpurea Girard 1856 (both of which are now regarded as valid species within the genus Gila). Evermann & Rutter (1895) placed T. intermedia within the genus Leuciscus Cuvier 1816. Jordan & Evermann (1896) synonymized Tigoma and Richardsonius Girard 1856 with Leuciscus, and suggested that L. intermedius and L. niger (Cope in Cope & Yarrow 1875) (now Gila nigra Cope in Cope & Yarrow 1875) may be conspecific. Fowler (1924) retained this genus placement for L. intermedius. Snyder (1915) regarded T. intermedia and G. nigra as synonyms of Gila gibbosa Baird & Girard 1854, which he placed within the genus Richardsonius Girard 1856. Jordan & Gilbert (1883) asserted that gibbosus Baird & Girard 1854 was unavailable due to homonymy resulting from placement of Leuciscus gibbosus Ayres 1854 and Gila gibbosa Baird & Girard 1854 within the genus Squalius Bonaparte 1837 (Gila gibbosa Baird & Girard 1854 being the junior homonym). Jordan & Evermann (1896) placed both of these species within the genus Leuciscus Cuvier 1816, creating another case of homonymy, this time within Leuciscus. In both cases, gibbosus Ayres 1854 was identified as the senior homonym; however, that name was itself preoccupied by Leuciscus between gibbosus Storer 1845 and gibbosus Ayres 1854. In any case, G. gibbosa Baird & Girard 1854 is not an available name. Miller (1945) treated intermedia as a subspecies of G. robusta Baird & Girard 1853a, which was followed by subsequent authors (Miller, 1946, 1961; Uyeno, 1961; La Rivers, 1962; Sigler & Miller, 1963; Miller & Lowe, 1964, 1967; Uyeno & Miller, 1965; Barber & Minckley, 1966; Cole, 1968; Minckley & Alger, 1968; Minckley, 1969; Lee et al., 1980; Robins, Bailey & Bond, 1980).

Rinne (1969) recognized two distinct species within the Lower Colorado River Basin: G. robusta for the more broadly distributed form, and G. intermedia for populations principally distributed in central and southern Arizona. Within the synonymy of G. intermedia, he also included G. gibbosus, G. nigra, and G. lemmoni, noting that the former was unavailable due to homonymy, thereby asserting G. intermedia as the next available name. Rinne (1969) does not appear to have examined any of the type series of G. intermedia, but states “in all respects, they (populations he labels as G. intermedia) correspond to the original description of Tigoma (=Gila) intermedia Girard 1856.” The original description of G. intermedia (Girard 1856) consisted of the following text: “Intermediate between T. pulchella and T. purpurea, more closely related however to the former than to the latter. The fins are much less developed, the inferior fins especially are quite small.” The original description did not include enough diagnostic characters to confidently conclude that the populations Rinne (1969) defines as G. intermedia are conspecific with the type specimens of G. intermedia, and there is no evidence that Rinne examined either pulchella or purpurea before asserting that G. intermedia was in fact, the correct name for these populations. Nevertheless, most subsequent authors followed Rinne (1969) in treating G. intermedia as a valid species (Stout, Bloom & Glass, 1970; Minckley, 1971, 1973; Rinne, 1976; Hocutt & Wiley, 1986; Minckley, Hendrickson & Bond, 1986; Sublette, Hatch & Sublette, 1990; Robins, 1991; Winfield & Nelson, 1991; Espinosa-Pérez, Gaspar-Dillanes & Fuentes-Mata, 1993; Gilbert, 1998; Minckley & DeMarais, 2000; Norris, Fischer & Minckley, 2003; Nelson et al., 2004; Scharpf, 2005; Miller, 2005; Minckley & Marsh, 2009; Page & Burr, 2011; Page et al., 2013), despite the known complex taxonomic history and a lack of any discrete identifying characteristics.

G. grahamii Baird & Girard 1853b

Although Gila grahamii is currently recognized as a synonym of G. robusta, due to the complicated systematic history of this group, particularly with respect to G. nigra, we include this nomenclatural account to promote clarity. Baird & Girard (1853b) described G. grahamii (often misspelled in the literature as grahami) from specimens collected in the Rio San Pedro, Gila River basin. Günther (1868) placed it within the genus Leuciscus Cuvier 1816. Cope & Yarrow (1875) placed the species back in the genus Gila and recognized it as distinct from both G. robusta Baird & Girard 1853a and G. nigra Cope in Cope & Yarrow 1875. Evermann & Rutter (1895) treated G. grahamii as a synonym of G. robusta, and subsequent early authors followed this assignment. This synonymy remained stable until Rinne (1969), who regarded grahamii as a subspecies of G. robusta, recognized populations collected from the tributaries of the Verde River and the upper Gila River system as distinct from the subspecies G. r. robusta in the main stem Verde and Gila Rivers. Rinne (1969) recognized current San Pedro populations (type locality of G. grahamii) as belonging to the species he referred to as G. intermedia, even though he accepted G. r. grahamii for his “headwater” form (again, apparently without examining type specimens). In the years between 1969 and 2000, there was not consistent recognition of G. r. grahamii as a valid subspecies, but no authors treated it as a valid species, or as a synonym or subspecies of any species other than G. robusta (Rinne, 1976; Lee et al., 1980; Robins, Bailey & Bond, 1980; Holden & Minckley, 1980; Mayden, 1992; La Rivers, 1994; Gilbert, 1998). Minckley & DeMarais (2000) regarded the populations referred to by Rinne (1969) as G. r. grahamii as representing a distinct species, even though they noted that it is likely of hybrid origin and paraphyletic. However, they also noted that the syntypes of G. grahamii belong to what Rinne (1969) regarded as the subspecies G. r. robusta (citing personal communication between RR Miller & WL Minckley). Therefore, they recognized G. nigra Cope in Cope & Yarrow 1875, as the earliest available name for the species previously referred to by Rinne (1969) as G. r. grahamii.

G. nigra Cope 1875

Gila nigra Cope in Cope & Yarrow (1875) was described from specimens collected in Ash Creek and at San Carlos, Arizona. Jordan & Gilbert (1883) placed it in the genus Squalius Bonaparte 1837, and later Jordan & Evermann (1896) placed it in the genus Leuciscus. They also regarded Gila gibbosa Baird & Girard 1854 as conspecific, but unavailable due to homonymy (see above). Gilbert & Scofield (1898) synonymized G. nigra with T. intermedia, which they placed in the genus Leuciscus. Snyder (1915) regarded G. niger as a synonym of G. gibbosa (within the genus Richardsonius), failing to recognize that the latter species name was not available due to homonymy. Fowler (1924) followed Gilbert & Scofield (1898) in treating G. nigra as a synonym of intermedius, within the genus Leuciscus. Subsequent treatments placed nigra in synonymy with intermedia (see above). This synonymy was broadly followed until Minckley & DeMarais (2000) recognized G. nigra as the earliest available name to refer to the species treated by Rinne (1969) as G. r. grahamii.

Synonymies now considered valid

Along with the considerable synonymy of Gila robusta are species that were at one time considered synonyms of G. robusta and are now considered valid. Gila elegans Baird & Girard 1853a was treated as a synonym of G. robusta by Ellis (1914) and as a subspecies by Miller (1945) and La Rivers (1994). However, G. elegans was treated as valid by Vanicek (1967) and subsequent authors. Gila jordani Tanner 1950 was treated as a subspecies of G. robusta by Rinne (1976), Lee et al. (1980), La Rivers (1994), and Gilbert (1998) but was recognized as valid by Minckley & Marsh (2009) and subsequent authors. Clinostomus pandora Cope 1872 was treated as a synonym of G. robusta by Ellis (1914) but valid as Gila pandora (Cope 1872) by subsequent authors. Tigoma egregia Girard 1858 was treated as a synonym of G. robusta by Ellis (1914) but regarded as valid as Richardsonius Girard 1856 egregius (Girard 1858) by subsequent authors. Finally, Gila seminuda Cope & Yarrow 1875 was treated as a subspecies of G. robusta by Ellis (1914), Snyder (1915), Rinne (1976), and Lee et al. (1980) but regarded as valid by Gilbert (1998) and subsequent authors.

Mitochondrial genome

Raw Illumina reads were paired, trimmed and mapped to the mitochondrial genome of G. robusta (GenBank DQ536424.1) using Geneious v.6.1.8 (Kearse et al., 2012). Five mapping iterations were completed at high sensitivity. For each sample, consensus sequences of all contigs that successfully mapped to the reference genome were extracted, and all consensus sequences were aligned, manually inspected and low coverage regions removed in Geneious.

Reduced representation nuclear genome

Reads were trimmed, assembled, and genotyped using the dDocent pipeline (Puritz, Hollenbeck & Gold, 2014). Any loci appearing in less than 85% of individuals were excluded from these analyses. The paired-end overlapping read algorithm was used for de novo assembly. Clustering similarity of 0.9 and mapping parameters of A (match score) = 2, B (mismatch score) = 3, and O (gap penalty) = 4 were used. dDocent processing recovered 89,896 loci with an average read depth of 227 in all 19 individuals. Complex variants were decomposed using vcflib (Garrison & Marth, 2012) to deconstruct haplotypes and insertions and deletions were removed using VCFtools (Danecek et al., 2011). Any contigs that mapped to the mtDNA were removed from the dataset using VCFtools, and the remaining contigs were considered to represent only the nuclear genome. The resulting data were then collapsed into haplotypes with the rad_haplotyper pipeline (https://github.com/chollenbeck/rad_haplotyper; Willis et al., 2017), which uses read alignments to record combinations of SNPs present across paired-end reads. For each individual, rad_haplotyper removes complex loci, missing data, paralogs, and sequencing errors. Any locus that is not present in at least 14 of 19 individuals with a depth of coverage of at least 20 is not included in the final dataset. The rad_haplotyper method was employed to overcome many of the problems that can arise with SNP data in the absence of a reference genome, such as inflated homozygosity, artifacts, or inflated heterozygosity (Willis et al., 2017). Contigs were then collapsed into genotypes for final analyses. PGDspider v.2.1.1.3 (Lischer & Excoffier, 2011) was used to convert the dataset to the required file types for further analysis.

Phylogenetic analysis

The optimal model of sequence evolution was selected using the Akaike Information Criterion in JModelTest v.2.1.4 (Posada, 2008). GTR + G was found to be the best-fit substitution model for both the mtDNA and total evidence datasets. The HKY model was found to be the best-fit model for the nuclear dataset. To calculate posterior probabilities of clades, MrBayes v.3.2.6 (Huelsenbeck & Ronquist, 2001; Ronquist & Huelsenbeck, 2003) was used to run a 1,000,000-generation Markov chain implementing the best-fit model for the dataset. We used flat Dirichlet prior probability densities with an initial burn in of 250,000 generations. Trees were saved every 500 generations for a total sample size of 1,500 trees. A majority rule consensus tree calculated from the 2,000 sampled trees was used to determine the posterior probabilities of clades. Under these parameters standard deviations between independent runs stabilized and were all less than 0.01. Maximum likelihood (ML) analyses were conducted using the randomized accelerated maximum likelihood (RAxML) software v.8 (Stamatakis, 2014). Best-fit models and 30,000 bootstrap replicates were implemented for all datasets. Uncorrected pairwise divergence times for mtDNA were estimated using Mega v.7.0 (Kumar, Stecher & Tamura, 2016). Phylogenetic trees were constructed and visualized using FigTree v.1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/). Structure analysis and Discriminate analysis of principal components (DAPC) were also implemented (see Supplementary Materials).

SNAPP coalescent analysis

Species trees were estimated from the nuclear dataset with the SNAPP package in BEAST v.2.3.2 (Bouckaert et al., 2014). Polymorphic loci were extracted from the rad_haplotyper output and outgroups were removed. BEAUti (Drummond et al., 2012) was executed with the following parameters: all taxa were treated as distinct species/populations, the mutation rates u and v were calculated from the data, and default values were used for the exponential priors. The data were then analyzed in BEAST (Bouckaert et al., 2014). The MCMC chain was run for 10,000,000 generations, sampling every 1,000 generations with 300,000 preBurnin. The results of this analysis were visualized using DensiTree (Bouckaert & Heled, 2014), and presented as a cloudogram.

Molecular clock estimation

To estimate the time of coalescence, a Bayesian MCMC approach was implemented in BEAST on the mtDNA dataset, under a coalescent constant population model with a strict clock of 2% per million years (Brown, George & Wilson, 1979; Bowen et al., 2001; Reece et al., 2010). Simulations were run with default priors under the GTR + G model of mutation. Simulations were run for 10 million generations, sampling every 1,000 generations following a 30,000 preBurnin. A total of 10 independent runs were computed to ensure convergence and log files were combined using TRACER v.1.7 (Rambaut et al., 2018).

Tests of introgression

Tests of introgression were performed on the total dataset using ABBA–BABA statistics implemented in HybridCheck (Ward & Van Oosterhout, 2016) to assess whether introgression with G. cypha could be detected within the G. robusta clades. To ensure that our data met the assumptions of ancestry ((P1, P2)P3), required for this test, we only compared in-groups (P1, P2) separated by the distinct clades of our phylogenetic trees. Standard ABBA–BABA calculations use Patterson’s D to infer introgression, but this statistic may not be sufficient to separate introgression from ancestral population structure. Therefore, we also include an alternative statistic F_d (Martin, Davey & Jiggins, 2014), which estimates the fraction of the genome shared through complete introgression between P2 and P3, and P1 and P3 and is not subject to the same biases as the D statistic. We also report the statistical significance (p-value) that expressed the deviation from an equal number of ABBA–BABA sites. Significant admixture was determined by Z scores of 3 or higher.

Type material

Examination of the type series of G. robusta (N = 2), G. intermedia (N = 4), G. nigra (N = 5), and G. grahamii (N = 1) reveal differences between the types within this complex. However, similar or greater morphological dissimilarity was observed for the G. grahamii type (now considered a synonym of G. robusta; see above) when compared to each of the other type series (Table 1). Furthermore, these differences apply to only these restricted type series and differences were not supported when the fresh material was added (see below). Unfortunately, the taxonomic key (Minckley & DeMarais, 2000) used to assign the names G. intermedia and G. nigra to the populations that Rinne (1969, 1976) and Minckley & DeMarais (2000) recognize as distinct species fails to correctly assign the type material to the correct species.

Fresh material

Analysis of fresh material of specimens assigned to the three species reveals extensive overlap in characters, prohibiting any definable difference between groups (Table 2). There is no single diagnostic character that can be used for species identification of fresh material, with considerable overlap among species in every morphological character. Likewise no suite of characters can distinguish the fresh material by species unambiguously.

Comparisons of fresh specimens to type material

Morphological comparisons of type material (Table 1) to fresh specimens (Table 2) also failed to resolve the species as currently recognized (Table 3). Type specimens as well as fresh material exhibit as much or more variation within species as between species. As such, it is impossible to assign any of the fresh specimens back to the type material and thus to species, without location information. No character in putative G. robusta specimens could be uniformly assigned back to the type of G. robusta but instead each character was assigned to multiple type series. Only 28% of the time did a morphological character (Tables 1 and 2) align correctly to the types of G. robusta, while 63% of the time the characters aligned to G. nigra types and 51% to G. intermedia, with many of the characters aligning with multiple type series (Table 3). Similar patterns are observed with both G. nigra and G. intermedia.

Molecular analysis

Mapping of the mitochondrial genome resulted in recovery of 14,892 bp of the 16,595 bp reference mtDNA genome. Filtering and haplotyping of the nuclear DNA resulted in a final dataset of 1,292 RAD contigs containing 4,821 haplotypes across all individuals, which consisted of 6,658 polymorphisms across 52,483 total bp of nuclear DNA. Analysis of the mitochondrial genome and nuclear datasets revealed high concordance between tree topologies constructed for each of the datasets using both ML and Bayesian tree building methods. Only a single discrepancy was observed between the nuclear dataset and mtDNA dataset: the O’Donnell Canyon specimen is assigned to clade 2 in the nuclear tree and clade 4 in the mtDNA tree (Figs. 2 and 3), which may be the result of accidental mixing of captive stocks. The results of these molecular analysis are consistent with the morphological finding; G. intermedia and G. nigra are not distinct evolutionary units and not distinguishable from G. robusta. Likewise, the SNAPP coalescent analysis results plotted as a cloudogram reveal high concordance between tree topologies (Fig. 4). Similarly, Structure (Pritchard, Stephens & Donnelly, 2000) analysis and discriminate analysis of principle components (DAPC; Jombart, 2008) fail to discriminate the nominal taxa (Fig. S2).

Morphological analysis of phylogeny

Comparison of morphological characters to the phylogenetic trees revealed nearly 100% overlap of morphological characters between each clade (i.e., there was no differentiation in any morphological character between phylogenetic clades; Table 4) and no diagnostic morphological character was identified that could align with the phylogenetic lineages resolved in this study.

Test of introgression

The ABBA–BABA test between the G. robusta complex and G. cypha revealed significant introgression between many of the locations (Table S2) using Patterson’s D statistic. However, as expected, the alternative, F_d suggests far fewer locations exhibiting significant introgression. There is evidence for introgression in populations from the Verde River (133), Aravaipa Creek (141), Hot Springs Canyon (125), and Spring Creek-Tonto (132) for both D and F_d.