Introgression between ecologically distinct species following increased salinity in the Colorado Delta- Worldwide implications for impacted estuary diversity

Sample collection

Specimens of the genus Colpichthys were collected by seine from estuarine habitats along the coast of the Northern Gulf of California and within the Colorado River Delta in the years 2005, 2006, 2007, and 2011 (Table 1). Collections were carried under Mexican Federal Collecting Permit (Permiso de Pesca de Fomento) DGOPA 14253. 101005.6950, and its extension DGOPA 06435.210606.2640, issued to Findley and Jacobs by the Comisión Nacional de Acuacultura y Pesca of the Secretaría de Agricultura, Ganadería, Desarrollo Rural, Pesca y Alimentación (SAGARPA). Collected specimens were preserved in 95% ethanol in the field and stored at −20 °C upon returning to the laboratory for use in both molecular and morphological analyses. For this study, we divided the coastline of the Northern Gulf of California into 4 regions: (1) the Delta, native to populations of C. hubbsi; (2) the Baja Coast and (3) the Sonora Coast, native to populations of C. regis; and (4) the Delta Edge, the putative hybrid zone. A full list of collection sites along with the abbreviated designation for each site are shown in Table 1. In addition to the field-collected specimens, formalin-fixed specimens from the Natural History Museum of Los Angeles County (LACM) and Scripps Institution of Oceanography were examined for their morphology to confirm species identification. A list of these “historical” specimens is included in the supplementary materials (Table 2).

Phenotypic analyses

The dorsal-lateral scale counts, vertebral counts, hemal spine morphology, and the position of the first dorsal fin were used as species-distinguishing features for this study; these features were chosen based on the species description of Colpichthys hubbsi (Crabtree, 1989). The first dorsal fin position is assessed based on the linear length from the snout to the first dorsal fin origin, normalized by the standard length. Length measurements were taken with dial calipers. Morphometric and meristic counts were performed according to Hubbs, Lagler & Smith (2004). Radiographs of specimens were taken for examination of their vertebral counts and hemal spine morphology. Six specimens from location NWC and three specimens from location YAV were too small to accurately assess their hemal spine morphology; they were thus excluded from the hemal spine data set.

Molecular procedures

Genomic DNA was isolated from caudal peduncle muscle tissues using the DNeasy Blood and Tissue Kit (Qiagen, Inc., Valencia, CA, USA). The mitochondrial cytochrome b gene was amplified by PCR using primers AJG15 and H5 (Akihito et al., 2000) under the following cycling conditions: 95 °C for 2 min; 40 cycles of 95 °C for 30 s, 51 °C for 30 s, followed by 72 °C for 90 s; then finally 72 °C for 10 min. PCR products were size-checked on 1.5% agarose gel and cleaned with ExoSap before Sanger sequencing with BigDye Terminator v3.1. Samples were then submitted to DNA Analysis Facility on Science Hill at Yale University for capillary electrophoresis.

Amplification of the nuclear RAG1 gene was carried out using primers RAG1-2533F and RAG1-4090R (López, Chen & Ortí, 2004) under the following cycling conditions: 95 °C for 2 min; 40 cycles of 95 °C for 1 min, 51 °C for 1 min with an incremental increase of 0.5 °C per cycle until 55 °C , followed by 72 °C for 90 sec; then finally 72 °C for 10 min. Post-PCR processing and sequencing of the RAG1 gene were performed as described above. From the sequences, we were able to identify four polymorphic sites; several individuals display ambiguity in base calling at such sites, and these were tagged as putative heterozygotes. A total of four RAG1 variants were inferred from the sequence data. Homozygous individuals were used to confirm the existence of three RAG1 variants; the remaining RAG1 variant was inferred from a heterozygous individual. To confirm the gametic phase of these RAG1 sequences, we employed allele-specific polymerase chain reaction (ASPCR) as detailed in Wu et al. (1989). Three allele-specific primers (T717F: 5′-CTACAAAATCTTCCAGGAT-3′, T770R: 5′-TTTATCTAAGGCTGCCCTCCAGA-3′, and A909F: 5′-AACTGGTGCCCTCAGAAGAA-3′) were designed to uniquely amplify one of the two diploid copies of RAG1 in the putative heterozygotes. The 3′ terminating nucleotide of these primers are designed to specifically anneal at the predetermined polymorphic site, thus enabling unique amplification of alleles and determination of the RAG1 genotypes of the heterozygotes.

Six microsatellite loci (Odont08, Odont09, Odont11, B18, B19, B39) were amplified for genotyping. Primer sequences for loci Odont08, Odont09, and Odont11 were obtained from Beheregaray & Sunnucks (2000). Primer sequences for loci B18, B19, B39 were obtained from Byrne & Avise (2009). Loci which can be easily differentiated by lengths were multiplexed two at a time using QIAGEN Multiplex PCR kits with 6-FAM labeled M13 primer (Boutin-Ganache et al., 2001). Only the forward primers of each locus are 5′ tagged with the M13 sequence, and reactions were performed under the following cycling conditions: 95 °C for 15 min; 25 cycles of 94 °C for 30 s, 55 °C for 90 s, 72 °C for 1 min; then another 25 cycles of 94 °C for 30 s, 50 °C for 90 s, 72 °C for 1 min; and finally 60 °C for 30 min. Diluted PCR products were submitted to the UCLA GenoSeq core facility for genotyping on an ABI3730 (Applied Biosystems, Foster City, CA, USA).

Molecular analyses

All sequences were aligned, trimmed, and analyzed using Geneious Pro ver. 5.5.6 (Drummond et al., 2011). A Maximum Likelihood tree of Cyt b sequence was made with MEGA7 (Kumar, Stecher & Tamura, 2016) and rooted using Atherinops affinis as the outgroup. A time-calibrated tree was inferred using the Reltime method (Tamura et al., 2012) and the Tamura-Nei model (Tamura & Nei, 1993). The timetree was computed using one calibration constraint at the node of the split between C. regis and C. hubbsi; the calibration point of 4.83Ma was obtained from the age of the ash bed within the Bouse Formation, the rock unit from which two fossilized C. regis were documented (Todd, 1976; Spencer et al., 2013).

Microsatellite data were obtained from 101 out of 103 specimens and imported into Geneious for allele scoring and analyzed using STRUCTURE ver. 2.3.4 (Pritchard, Stephens & Donnelly, 2000). STRUCTURE is a program that uses a Bayesian clustering approach to infer population structure from genotypic data. Under an admixture model of K populations or genetic groupings, STRUCTURE computes the q-value, a quantity between 0 and 1, that reflects the proportions of an individual’s genome originating from a certain population. Hybridization can then be inferred based on the computed q-value; for example, first-generation (F1) hybrids of two populations (K = 2) are expected to have a q-value close to 0.5 (Vähä & Primmer, 2006). We ran STRUCTURE analyses using an admixture model with 1,000,000 MCMC generations after a burn-in period of 100,000. Sampling locations were incorporated as prior parameters in the LOCPRIOR model to improve clustering for our data set; the LOCPRIOR model is preferred when there are weak signals of structure due to low number of markers analyzed (Hubisz et al., 2009). Both correlated and independent allele frequency models were used with number of groups set from K = 1 to K = 4; for each K value, analysis was repeated ten times. Of the two allele frequency models, the correlated model provides greater power in differentiating closely related populations (Falush, Stephens & Pritchard, 2003). The optimal K value was informed by ΔK (Evanno, Regnaut & Goudet, 2005) as computed using Structure Harvester (Earl & VonHolt, 2012). Output of the STRUCTURE analyses were summarized in CLUMPP ver. 1.1.2 (Jakobsson & Rosenberg, 2007) and visualized in DISTRUCT ver. 1.1 (Rosenberg, 2004).

Multivariate analysis

We further analyzed our phenotypic and microsatellite data by discriminant analysis of principal components (DAPC) using the R-package adegenet (Jombart & Ahmed, 2011). This analysis is used to help describe the phenotypic and genetic clustering of our specimens. Cluster priors for DAPC were identified by k-means clustering; based on BIC scores we selected the 2-cluster model to use for the DAPC. To perform DAPC, the program transforms the data into principal components, which is then analyzed through discriminant analysis. Repeated cross-validation was used to evaluate the optimal number of principal components to retain for the discriminant analyses to avoid overfitting of data. For the phenotypic analysis, we used measurements of the four species-distinguishing features as our data and retained one PC (91.28% of the cumulative variance) for generating the morphology discriminant function. For the microsatellite analysis, we used the allele frequencies of the six loci as our data and retained 20 PCs (66.35% of the cumulative variance) for generating the microsatellite discriminant function.

Demographic histories

To explore the demographic histories of Colpichthys, Cyt b data sets were separated by species according to the Cyt b gene tree (Fig. 1B). Forty-four sequences from C. hubbsi and 52 sequences from C. regis were analyzed for their haplotype diversity. Median-joining haplotype networks were drawn in POPART ver. 1.7 (Leigh & Bryant, 2015). Cyt b haplotypes were examined by various neutrality tests—Tajima’s D (Tajima, 1989), Fu and Li’s D* and F* (Fu & Li, 1993), Fu’s Fs (Fu, 1997), and R ₂ (Ramos-onsins & Rozas, 2006)—in DnaSP v. 5.10.01 (Librado & Rozas, 2009). P-values for these test statistics were obtained by coalescent simulations of 1,000 replicates, conditional on the number of segregating sites. Pairwise differences were analyzed by mismatch distributions in Arlequin ver. 3.1 (Excoffier, Laval & Schneider, 2005); coalescent simulations were run with 10,000 bootstrap replicates. Past demographic expansions can be inferred from haplotypes based on the mismatch distribution of the pairwise differences; an uneven, multi-modal distribution is expected for populations under demographic equilibrium, whereas a “smooth” uni-modal distribution is expected for populations that had undergone demographic expansion (Rogers & Harpending, 1992). The observed mismatch distributions of the two species were compared to the expected mismatch distribution under a sudden expansion model using the sum of square deviation (SSD) and the raggedness index.

Pure-bred Specimens of C. regis and C. hubbsi are morphologically distinct

Consistent with the findings of Crabtree (1989), our morphological analysis based on both modern and historical specimens, summarized in Fig. 2, shows that C. regis (green and gold) differs from C. hubbsi (pink) in having higher dorsal-lateral scale counts, higher vertebral counts, and shorter snout to first dorsal fin origin distance. In addition, the symphases of the hemal spines of the anterior-most caudal vertebrae of C. regis are modified into an expanded process (Fig. S1), forming a funnel-like structure antero-posteriorly. The number of modified hemal spines within C. regis outside of the Delta Edge varies between 5 to 11; the hemal spines of C. hubbsi in the Delta appears to show no signs of such modification (Fig. 2D). Outside the putative hybrid zone at the Delta Edge, individuals are relatively easily classified to species based on morphological differences. However, more continuous variation is evident within the Delta Edge region.

Introgression of Cytochrome b and RAG1 haplotypes

Maximum-Likelihood phylogeny reconstruction of Cyt b sequence reveals ∼6% divergent clades corresponding to the two species of Colpichthys (Fig. 1B). C. hubbsi haplotypes were recovered from the Delta region, while Baja and Sonoran Coast samples contained only C. regis haplotypes. Both species were expected to cohabit the localities at the southwestern border of the Delta (PRI, SGU, and TRC), but none of these locations contain a mixture of both haplotypes. PRI and SGU fish are all C. regis haplotypes, while the more Delta internal TRC sample is exclusively C. hubbsi haplotypes, albeit with limited sample sizes. Within the Delta, five specimens collected from predominantly C. hubbsi populations—ELZ, CHY, and SFG—were found to possess the C. regis mitochondrial haplotype; these specimens are morphologically indistinguishable from other C. hubbsi specimens.

We identified four RAG1 variants (labeled as H₁, H₂, H₃, and R; Table 3) based on four polymorphic sites in the RAG1 sequence data (975 bp). We find that H₁, H₂, and H₃ are common in populations of C. hubbsi in the Delta, while R is fixed in all populations of C. regis in the Baja Coast and Sonora Coast regions. Thus, the two species are well differentiated by their RAG1 genotypes, with the diagnostic SNPs at positions 717 and 909. Six specimens, heterozygous at both of these positions, are flagged as putative hybrids. Four of the six hybrids occur at the Delta Edge, while the remaining two are found within the Delta proper (Table 4). No RAG1 heterozygotes are detected in any of the Baja and Sonoran collection sites.

As noted above, morphology of the anterior-most caudal hemal spines provides a phenotypic difference between the two species of Colpichthys. The non-modified hemal spine phenotype (Fig. S1A) is observed in all specimens possessing any of the three C. hubbsi RAG1 alleles (H ₁, H ₂, or H ₃).

Bayesian analyses of microsatellites show evidence of admixture

Based on the evaluation of ΔK, the STRUCTURE analyses on the six microsatellite loci recovered two genetically distinct clusters corresponding to the two species of Colpichthys (Fig. 3; for K = 3 and K = 4, see Fig. S2). Under the admixture model, the q-value represents admixture proportion of each individual. We used a q-value threshold to distinguish putative purebred and hybrid individuals; for optimal efficiency, this threshold is set to 0.1 (Vähä & Primmer, 2006). Populations in the Delta proper—ELZ, SFG, and CHY—were identified as purebred C. hubbsi. Several populations away from the Delta such as ANI, NWC, and YAV were identified as purebred C. regis. All individuals within or near the western Delta Edge—PRI, SGU, and TRC—show substantial degrees of genetic admixture, having q-values between 0.1 and 0.9. Moreover, among these populations there appears to be a gradient in q-values as a function of geographic distance from the mouth of the Colorado River. TRC has an average q-value of 0.17 while PRI, which is farther from the river mouth, has an average q-value of 0.72 (Fig. 3A). A relatively weak signal for admixture in the ADR population beyond the eastern edge of the delta was detected exclusively in the correlated allele frequency model; this signal is absent in the analysis with the independent allele frequency model (Fig. 3B).

Consilience of phenotypic and nuclear data

Results from separate DAPC analyses were pooled to compare phenotypic and microsatellite variation. Two clusters are supported, consequently only one discriminant function was generated for each of the DAPC analyses. These results show two well separated species with congruent cluster assignment from both data sets (Fig. 3C). Delta Edge specimens preferentially situate between the species clusters. Of the six RAG1 heterozygotes (triangles on Fig. 3C): two recovered in the Delta Edge area fall near the midpoint in morphology between the two species; two recovered in the Delta proper are near the centroid of C. hubbsi and the last two from the Delta Edge are phenotypically more similar to C. hubbsi than C. regis.

Genetic diversity and demographic history

The Cyt b median-joining haplotype networks of the two species are displayed in Figs. 4A and 4B. Twenty-Eight C. hubbsi haplotypes are recovered from the four Delta populations; in contrast, only ten haplotypes are found in C. regis over a larger, more disjunct geographic range. C. regis haplotype diversity differs substantially between the Baja and the Sonoran Coasts; with the exception of one individual from PRI, all non-major haplotypes are found exclusively in Sonora Coast populations. The major C. regis haplotype is present in all locations outside the Delta except the two southernmost Sonora Coast sites, SLD and YAV. All five individuals within the Delta proper with the C. regis Cyt b sequence possess the major haplotype. Neutrality tests were conducted with the exclusion of the mitochondrial-introgressed individuals; all neutrality tests report significant departure from neutrality. Significantly negative values for Tajima’s D and Fu’s Fs in both species are consistent with either population growth or selective sweep (Table 5). The sum of square deviations and raggedness indices for the mismatch distributions are close to zero and not significant (Table 6), thus we are unable to reject the sudden expansion model for both species. Figures 4C–4D shows the nucleotide mismatch distributions for both species; C. hubbsi has a bimodal distribution, while C. regis has a unimodal distribution that matches closely with the expected curve under sudden expansion.

The hybrid zone

Crabtree (1989) recognized that the two Colpichthys species co-occur in the region just north of San Felipe, designated the Delta Edge in this study. While external morphology of specimens collected in 1968 confirms the co-occurrence of the two species at this region (Table S1), there is no indication from these samples of historical hybridization. Our study presents the first evidence for the modern hybridization of the sister species. STRUCTURE and DAPC analyses strongly support our hypothesis that the Delta Edge region currently represents a hybrid zone. STRUCTURE results limited to the correlated allele frequency model raise the possibility that hybridization may also be occurring on the eastern border of the Delta at Bahía Adair (ADR). The fact that weak admixture signals at ADR are not detected using our other approaches suggests that hybridization at this region, if present, is likely infrequent.

Directional introgression

In addition to hybridization at the Delta Edge, we detected C. regis mitochondrial and RAG1 haplotypes near the mouth and upstream along the Colorado Estuary within the Delta. This could be explained by either the dispersal of first generation hybrids from the hybrid zone to the Delta or by multigenerational backcrossing of hybrids with parental C. hubbsi in this case far from the Delta Edge (i.e., introgressive hybridization). Our data lend support to the introgressive hybridization interpretation. Unlike the RAG1 hybrids at the Delta Edge, specimens that possess C. regis mitochondrial and/or RAG1 haplotypes within the Delta do not exhibit intermediate morphology. Morphological and microsatellite data for these specimens are typical for C. hubbsi. Thus, they likely represent the product of multiple backcrosses. Thus C. regis haplotypes have been introduced into the center of the C. hubbsi range via introgressive hybridization. Conversely, we failed to detect any C. hubbsi haplotypes outside of the Delta and Delta Edge regions. Thus, introgression does not appear to be reciprocal, although more sampling south of the Delta Edge would allow a stronger conclusion that C. hubbsi haplotypes are absent in predominantly C. regis populations.

Parental sex asymmetry in introgression

Phenotypic and microsatellite data (Fig. 3C, Fig. S3) reveal a narrow geographic cline across the 10 km wide hybrid zone which spans three estuaries at the southwestern Delta border—Estero Primero (PRI), Estero Segundo (SGU), and Estero Tercero (TRC). Yet despite morphological and nuclear evidence of hybridization all thirteen specimens from PRI and SGU possess the C. regis mitochondrial haplotype. Given matrilineal mitochondrial inheritance, this pattern may indicate hybridization between female parental C. regis and male parental C. hubbsi. Little is known about the life history and spawning behavior of Colpichthys. However, Crabtree (1989) noted that the ovipositor of C. hubbsi is twice the size of that in C. regis, suggesting that differences in reproductive organs could contribute to the directionality of introgression. Behavioral or ecological differences may also play a role, but these hypotheses have yet to be assessed.

Timing of divergence

Due to the lack of available genetic data from historical samples, the onset of hybridization cannot be easily constrained. However, the ∼6% mitochondrial sequence divergence between C. regis and C. hubbsi suggests that the two species have long been separated with hybridization a relatively recent phenomenon. This mitochondrial divergence is comparable to that of another fish sister species pair, Gillichthys mirabilis and Gillichthys detrusus (Swift et al., 2011), which share distributional patterns with C. regis and C. hubbsi. Like C. hubbsi, G. detrusus is narrowly restricted to the silty tidal channels of the Colorado River Delta; its congener, G. mirabilis, can be found in the estuaries and lagoons on both coasts of the Gulf. Swift et al. (2011) estimated a ∼5 Ma divergence time between G. detrusus and G. mirabilis. Given that both C. hubbsi and G. detrusus are Delta endemics, their divergences from their sister species may have been linked to the formation of the Delta estuaries, which occurred after the opening of the Gulf of California. Fossils of C. regis have been discovered in the Bouse Formation at the location of Cibola Lake, Arizona (Todd, 1976). The Bouse Formation records the changing depositional environments associated with Late Miocene/Early Pliocene progression of the Colorado River to the Gulf (Poulson & John, 2003; McDougall & Martínez, 2014; Spencer & Patchett, 1997; Roskowski et al., 2010). Stratigraphic and tectonic reconstructions suggest that the Colorado River reached the Gulf of California by ∼4–5 Ma (Winker & Kidwell, 1986; Dorsey et al., 2007; Spencer et al., 2013; Crossey et al., 2015; Howard et al., 2015). Based on this chronology, speciation leading to the Delta endemics presumably followed the establishment of the Delta habitats around this time (Swift et al., 2011; Ellingson, 2012).

Selection and ecological speciation in response to salinity conditions

The parapatric distribution of Colpichthys in the Northern Gulf suggests that the two species may have diverged due to ecological differentiation reinforced by divergent selection on local adaptations. The divergence is likely linked to the gradient in salinity in the Northern Gulf. While virtually no physiological data exists for C. hubbsi, C. regis is known to tolerate hypersaline conditions as it is frequently observed in lagoons with salinity values reaching 50 psu (Castro-Aguirre & Espinosa Pérez, 2006). Prior to the completion of Hoover Dam in 1935, the water output at the Delta was estimated to be 16 to 18 billion cubic meters per year (Stockton & Jacoby Jr, 1976), supporting a substantial historic fresh to salt water gradient through much of the Delta. Development of an inverse hypersaline estuary followed in the 20th Century damming and water removal from the system (Lavín & Sánchez, 1999). Although we are not aware of any studies of salinity tolerance in captivity for Colpichthys, as noted above, silversides have developed ecological species in response to salinity gradients in Brazil and Australia, the Mediterranean and the Eastern US (Bamber & Henderson, 1988; Beheregaray & Sunnucks, 2001; Fluker, Pezold & Minton, 2011; Johnson, Watts & Black, 1994; Klossa-Kilia et al., 2007; Francisco et al., 2006; Olsen, Anderson & McDonald, 2016; Trabelsi et al., 2002). Thus, it is likely that historic salinity gradients in the Delta played a role in reproductive isolation between the Colpichthys sister species. Moreover, selection maintaining these species is prone to breakdown following the loss of freshwater.

Due to the difference in density of fresh and saltwater, salinity is known to influence the morphology of the gas bladder, which maintains neutral buoyancy in teleost fishes. Among its closest relatives in Atherinopsinae, C. hubbsi is the only species where the gas bladder does not extend posterior to the visceral cavity (Crabtree, 1989), suggesting that salinity related selection impacts the species. Other potential selective factors may involve the unusual silty conditions and tides in the Delta and the limited depth in the Delta. Relative to other examples of species boundaries particularly maintained by selection and differences in salinity in fishes such as those associated with the transition into the Baltic, which developed during the Holocene (Barrio et al., 2016); or those in Brazilian coastal silversides where divergence of freshwater forms was in the Pleistocene (Beheregaray, Sunnucks & Briscoe, 2002), our inferred early Pliocene timing of speciation based on sequence differences between C. regis and C. hubbsi is substantially more ancient.

Evolutionary implications of introgression for colpichthys

Introgressive hybridization potentially has a number of distinct causes, including anthropogenic introduction and habitat modification. It also has a number of different evolutionary outcomes ranging from sharing of advantageous alleles between species (Martin, Bouck & Arnold, 2006; Whitney, Randell & Rieseberg, 2010; Pardo-Diaz et al., 2012; Hedrick, 2013) to extinction (Rhymer & Simberloff, 1996; Todesco et al., 2016). Of particular interest here are examples of habitat alteration that lead to ”speciation reversal”, or the introgressive loss of ecological speciation. Eutrophication, for example, has led to numerous examples of species loss in fishes (Vonlanthen et al., 2012; Seehausen, Van Alphen & Witte, 1997). Given that salinity is associated with ecological species in a number of silversides, and the loss of freshwater in the Colorado Delta, introgression into the more freshwater adapted C. hubbsi from the marine adapted C. regis would be expected. Ultimately, introgression of C. regis alleles adapted for higher salinity could leave C. hubbsi at risk of genetic swamping. In the long run, the risk of extinction by speciation reversal for C. hubbsi will depend on the extent to which the loss of salinity gradient has eliminated the habitat distinction in the Northern Gulf. While this study only identified backcrossed hybrids based on the limited genetic loci sampled, it is possible that many of the morphological “pure bred” C. hubbsi are in fact “cryptic” hybrids. Depending on the frequency and relative abundance of such “cryptic” hybrids, interspecific gene flow may be more pervasive than suggested by our data. To assess whether “speciation reversal” could pose a threat to C. hubbsi, future whole-genome analyses will be required to quantify the amount of genomic introgression as well as identify loci that may remain under divergent selection in the face of gene flow.

Demographic expansions

Given the geographic restriction of C. hubbsi, one immediate concern is whether the reduced Colorado River outflow has directly contributed to a population decline. The rise in salinity at the Delta has been implicated in the local decline of the clams in the genus Mulinia (Rodriguez, Flessa & Dettman, 2001). Based on our results from neutrality tests and mismatch analyses, we find significant departure from neutrality, but no evidence to support a recent bottleneck in C. hubbsi. On the other hand, contrary to initial expectation, C. regis exhibits surprisingly low genetic diversity compared to its geographically restricted congener. This pattern is seen in all three types of molecular markers (Cyt b: Fig. 4; RAG1: Table 3; microsatellites: Table S2). The results from this study are consistent with demographic expansions in the recent past. Other factors, such as selective sweep or fine-scaled population structure (Ptak & Przeworski, 2002), could produce similar mtDNA patterns, but these mechanisms seem unlikely to have generated the similar patterns of diversity in all three data sets. We offer the following scenario as a possible explanation for the observed data. In the mismatch analysis, the model parameter τ estimates the age of demographic expansion. Based on this, our results (Table 6) suggest that C. hubbsi may have attained larger populations in the Delta estuarine habitats than C. regis and that C. hubbsi populations were either more continuous or experienced a more ancient expansion relative to C. regis. One possible explanation for this could come from Pleistocene climate and sea-level fluctuations. Recent work (Dolby et al., 2016) demonstrated that interactions between coastal geomorphology and Pleistocene sea-level fluctuations limited estuary habitats to refugia during glacial sea-level minima, impacting modern population structures of estuarine fish. In addition, during the Pleistocene, wetter climates increased Colorado River outflow into the Gulf, while the gentle sloping bathymetry of the Northern Gulf together with the lowering of sea-level likely expanded habitats for C. hubbsi much further south, allowing for population growth. Meanwhile, the steeper outer shelves on both sides Gulf to the south could have restricted suitable refugial habitats for C. regis during sea-level lowstand (Dolby, 2015). The lack of genetic diversity in C. regis populations may represent rapid range expansion from these limited refugia following sea-level rise at the end of the Pleistocene. Crucial to this hypothesis is the precise timing of the demographic and range expansions. More complete sampling of C. regis within each of the modern localities or genomic approaches are needed in future studies to further explore this issue.

Potential implications for Delta endemics

Of the numerous Gulf endemics (Palacios-Salgado et al., 2012), some are confined to the northernmost Gulf, and a smaller subset are very locally endemic to the Delta itself, which is an area protected as a Biosphere Reserve (Reservade la Biosfera del Alto Golfo de California y Delta del Rio Colorado). Of these northernmost endemics, some are relicts of broader distributions, while others appear to be ecological species that evolved in association with the Colorado Delta (Jacobs, Haney & Louie, 2004). Apparent relictual species include the critically endangered Vaquita, the world’s smallest cetacean, whose closest relative is found on the coasts of South America south of Peru (Munguia-Vega et al., 2007). The totoaba (Totoaba macdonaldi) was historically an important fishery resource that appears to have depended on freshwater input from the Colorado River for rapid growth early in development (Rowell et al., 2008). This genus, broadly endemic to the Gulf, is known from fossils in the Miocene of California Central Valley (Huddleston & Takeuchi, 2007); its closest relatives are in the Atlantic and IndoPacific (Lo et al., 2015). In addition to these relicts which may have been directly impacted by salinity changes or other anthropogenic activities such as overfishing, three taxa in the region appear to have evolved ecological species pairs associated with the Delta. Apart from the fish genera Colpichthys and Gillichthys as discussed above, the fiddler Crab Uca monilifera and U. princeps conform to the same pattern (Brusca, 1980). Thus, out of a suite of five northernmost Gulf endemic estuarine species, three are potentially at risk of introgression due to hydrologic modification. While a number of studies have argued that the northern Gulf fauna has been directly impacted by salinity (Rodriguez, Flessa & Dettman, 2001; Rowell et al., 2008), others have emphasized the minimal change in Northern Gulf biological productivity associated with the anthropogenic loss of freshwater input from the Colorado River (Brusca et al., 2017). However, there may nevertheless be strong impacts on the endemic diversity of this biologically unique region, especially with regards to introgression in endemics specific to the Colorado Delta proper.

Global change context

Early 20th century water removal tended to impact larger river systems primarily in semi-arid settings, where freshwater is at a premium and where economic development is well advanced; this is well exemplified by the Southwestern US. With increasing population growth and economic development in a wider range of geographic settings, damming and diversion efforts have continued to expand from rivers that traverse arid regions, impacting the Deltas estuaries of rivers such as the Nile (Stanley & Warne, 1998) and the Euphrates (Isaev & Mikhailova, 2009), to many of the largest river systems of Asia, such as the Mekong (Kondolf, Rubin & Minear, 2014). Despite this, the impact on estuarine diversity, especially as regards to the potential loss of diversity due to introgression between ecological species, has not been widely discussed or assessed. In the case of the Colorado River Delta, we have been able to establish that introgression is ongoing in the context of the ecologically altered Delta environment. This is the first step in establishing the risk of introgressive species loss in river estuaries that have suffered from water extraction and loss of freshwater input. This is yet another factor in a long list of ecological and biodiversity impacts consequent to altering riverine hydrologic processes, and should be systematically considered in a global context so that the full ecological consequences and diversity costs of such development can be properly understood.