Hybridization drives genetic erosion in sympatric desert fishes of western North America

Abstract

Many species have evolved or currently coexist in sympatry due to differential adaptation in a heterogeneous environment. However, anthropogenic habitat modifications can either disrupt reproductive barriers or obscure environmental conditions which underlie fitness gradients. In this study, we evaluated the potential for an anthropogenically-mediated shift in reproductive boundaries that separate two historically sympatric fish species (Gila cypha and G. robusta) endemic to the Colorado River Basin using ddRAD sequencing of 368 individuals. We first examined the integrity of reproductive isolation while in sympatry and allopatry, then characterized hybrid ancestries using genealogical assignment tests. We tested for localized erosion of reproductive isolation by comparing site-wise genomic clines against global patterns and identified a breakdown in the drainage-wide pattern of selection against interspecific heterozygotes. This, in turn, allowed for the formation of a hybrid swarm in one tributary, and asymmetric introgression where species co-occur. We also detected a weak but significant relationship between genetic purity and degree of consumptive water removal, suggesting a role for anthropogenic habitat modifications in undermining species boundaries or expanding historically limited introgression. In addition, results from basin-wide genomic clines suggested that hybrids and parental forms are adaptively nonequivalent. If so, then a failure to manage for hybridization will exacerbate the long-term extinction risk in parental populations. These results reinforce the role of anthropogenic habitat modification in promoting interspecific introgression in sympatric species by relaxing divergent selection. This, in turn, underscores a broader role for hybridization in decreasing global biodiversity within rapidly deteriorating environments.

Fig. 1

Sampling localities for Gila cypha (blue) and G. robusta (red) within the Colorado River Basin, western North America. Locality codes are defined in Table 2. Sympatric locations (BKR, DES, WWC, and YAM) are slightly offset for visibility purposes. Inset cartoons the respective morphologies of each species

Fig. 2

Results of a discriminate analysis of principal components (DAPC) analysis depicting Gila robusta (red), G. cypha (blue), and their respective populations (as colored). a Discriminant function axes 1 and 2 (=DF1 × DF2) showing discrimination among both species; b axes 2 and 3 (=DF2 × DF3) reflecting the manner by which populations of each species (grouped within ellipses) are distributed in discriminant space. The relative percent variance captured by each discriminant function is presented in parentheses. Sample localities are defined in Table 1. An asterisk denotes sympatric localities

Fig. 3

Assignment results for Admixture and Structure analyses involving Gila robusta, G. cypha, and G, elegans. K-values range from Structure optimum K = 5 (see Fig. S2) to Admixture optimum K = 8 (see Fig. S3). Locality abbreviations are as defined in Table 1. An asterisk denotes sympatric localities

Fig. 4

Genealogical assignment for individual Gila robusta and G. cypha, as compiled from NewHybrids analysis. Individuals are represented by colored bars, with proportion of color indicating posterior probability of assignment per genealogical class. Prior “parental” allele frequencies for G. cypha were derived from the Little Colorado River (HLCR) and from the Little Snake River (RLSR) for G. robusta (alternative prior assignments had no significant affect; see Fig. S4 and S5). Colors are as follows: Red = pure G. robusta; Blue = pure G. cypha; Purple = F₁ hybrid; Light purple = F₂ hybrid; Light blue = cypha-backcrossed hybrid; Light red = robusta-backcrossed hybrid. An asterisk denotes sympatric localities

Fig. 5

Genomic cline analyses for populations of Gila robusta and G. cypha, presented as: a Per-locus clinal relationships for 50 SNPs with δ > 0.8 (all significantly non-neutral at α = 0.001) compared with the neutral expectation (shaded gray region); b Log-likelihood ratio distribution of site-wise per-locus clines compared to the global pattern, where higher log-likelihood ratio indicates greater discordance; c Per-locus incongruence in genomic clines in the Gila robusta samples from the San Rafael River, partitioned by year (2009 versus 2017). Locality codes for populations of each species are defined in Table 1