Discordant introgression in a rapidly expanding hybrid swarm.

The erosion of species boundaries can involve rapid evolutionary change. Consequently, many aspects of the process remain poorly understood, including the formation, expansion, and evolution of hybrid swarms. Biological invasions involving hybridization present exceptional opportunities to study the erosion of species boundaries because timelines of interactions and outcomes are frequently well known. Here, we examined clinal variation across codominant and maternally inherited genetic markers as well as phenotypic traits to characterize the expansion and evolution of a hybrid swarm between native Cyprinella venusta and invasive Cyprinella lutrensis minnows. Discordant introgression of phenotype, microsatellite multilocus genotype, and mtDNA haplotype indicates that the observable expansion of the C. venusta × C. lutrensis hybrid swarm is a false invasion front. Both parental and hybrid individuals closely resembling C. lutrensis are numerically dominant in the expansion wake, indicating that the non-native parental phenotype may be selectively favored. These findings show that cryptic introgression can extend beyond the phenotypic boundaries of hybrid swarms and that hybrid swarms likely expand more rapidly than can be documented from phenotypic variation alone. Similarly, dominance of a single parental phenotype following an introduction event may lead to instances of species erosion being mistaken for species displacement without hybridization.

Introduction

Interspecific hybridization can either promote or constrain biological diversity. By producing novel epistatic interactions and genotypes, hybridization can lead to the formation of new species (Rieseberg 1997; Nolte et al. 2005; Nolte and Tautz 2010). Species boundaries can also erode following hybridization, leading to the rapid loss of genetic and species diversity (Seehausen et al. 1997; Taylor et al. 2006). Hybrid swarms, for example, can form quickly and overwhelm parental species through genetic homogenization or competitive exclusion in as few as five generations (Rhymer and Simberloff 1996; Mooney and Cleland 2001; Wolf et al. 2001; Perry et al. 2002; Hall et al. 2006), making it likely that the progressive stages of interactions and outcomes go undetected. Consequently, the process of species erosion remains poorly documented, including the formation, expansion, and evolution of hybrid swarms (Taylor et al. 2006; Seehausen et al. 2008; Gilman and Behm 2011).

The erosion of species boundaries is a dynamic process involving the conversion of initially discrete parental trait distributions to a single distribution of intermediate, hybrid forms over time, and space (Jiggins and Mallet 2000; Taylor et al. 2006). At the onset of secondary contact, trait distributions typically take the form of steep coincident spatial transitions from one parental species to another (Barton and Hewitt 1985, 1989; Arnold 1997). If reproductive isolation is incomplete and barriers to gene flow are weak, the frequency of later-generation hybrid and backcrossed individuals increases, leading to the amalgamation of distinct phenotypic and genetic entities (Forbes and Allendorf 1991; Taylor et al. 2006; Gilman and Behm 2011). Over time, the initially steep clines established by secondary contact erode (Endler 1977; Barton and Gale 1993), flattening individual genotype and trait distributions and expanding the spatial extent of introgression.

Hybrid zones have come to serve as natural laboratories for examining processes contributing to the maintenance or the erosion of species boundaries because the factors influencing the evolutionary outcome of hybridization can be inferred from distributions of observable traits (Barton 1983; Barton and Hewitt 1985, 1989; Barton and Gale 1993; Harrison 1993; Arnold 1997). Patterns of spatial variation reflect biotic and abiotic influences on the integrity of species boundaries, such as the extent of gene flow, the strength of reproductive isolation, and the relative fitness of hybrid and parental phenotypes in particular environments. Temporal information can provide additional insight into ecological or genetic factors promoting species erosion and how the influence of such factors changes across progressive stages of introgression (Barton and Hewitt 1981; Rieseberg 1991; Wolf et al. 2001; Hall et al. 2006; Lepais et al. 2009).

Most natural hybrid zones are not ideal laboratories for examining the erosion of species boundaries because they tend to be highly stable (but see Blum 2002; Dasmahaptra et al. 2002) and time since the formation is rarely known. Biological invasions involving hybridization, on the other hand, can provide exceptional opportunities to study the erosion of species boundaries in real time. Hybridization between non-native and native congeners is a well-documented outcome of species introductions (Ellstrand and Schierenbeck 2000; Sakai et al. 2001), particularly in aquatic ecosystems where rates of introduction continue to rise, and the time of initial introductions and onset of hybridization are often known (Cohen and Carlton 1998; Ellstrand and Schierenbeck 2000; Walters et al. 2008). Hybrid swarms have arisen rapidly following secondary contact, even when the initial densities of non-native species are low (Ellstrand and Schierenbeck 2000). For example, hybridization between the native crayfish Orconectes propinquus and the invasive congener Orconectes rusticus in lakes across Wisconsin (USA) led to the displacement of O. propinquus by hybrids phenotypically similar to O. rusticus within a decade of introduction (Perry et al. 2001). Similarly, heterosis and preferential choice of male sheepshead minnow (Cyprinodon variegatus) by female Pecos pupfish (Cyprinodon pecosensis) fueled the spread of hybrids across more than half of the range of the Pecos pupfish in <5 years (Echelle and Connor 1989; Wilde and Echelle 1992; Childs et al. 1996; Rosenfield and Kodric-Brown 2003; Rosenfield et al. 2004).

Here, we characterize a recently formed and rapidly expanding hybrid swarm (Walters et al. 2008) involving native blacktail shiner (Cyprinella venusta stigmatura) and introduced non-native red shiner (Cyprinella lutrensis) in the Upper Coosa river basin (Alabama, Georgia, Tennessee, USA). Cyprinella lutrensis thrive under a wide range of ecological conditions, exhibit high fecundity and a short generation time and frequently dominate species assemblages where introduced (Hubbs 1954; Marsh-Matthews and Matthews 2000; Marsh-Matthews et al. 2002; Herrington and DeVries 2008; Walters et al. 2008; Blum et al. 2010). Declines in the abundance or displacement of native stream fishes have been attributed to competition for resources and predation by non-native C. lutrensis on native larvae (Gregor and Deacon 1988; Ruppert et al. 1993; Douglas et al. 1994; Carpenter and Mueller 2008). Hybridization with native congeners is thought to exacerbate the impact of competition and predation on assemblage structure (Page and Smith 1970; Burr and Page 1986; Mettee et al. 1996; Walters et al. 2008). In this study, we examined clinal variation in codominant and maternal genetic markers, as well as quantitative phenotypic traits, to better understand the nature and outcomes of hybridization between C. lutrensis and C. venusta, including whether the expansion of the Coosa River C. lutrensis × C. venusta hybrid swarm is limited to the leading edge of the phenotypic invasion front. By examining the patterns of genetic and phenotypic introgression in the wake of upstream expansion, we also assessed whether the decline of native C. venusta following the introduction of C. lutrensis is best characterized as a process of displacement or loss through introgressive hybridization.

Materials and methods

Collections

Survey records indicate that C. lutrensis were introduced to the upper Coosa River basin in Lake Weiss (Alabama) in 1974 and have since dispersed upstream at rates approaching 31 km/year (Walters et al. 2008). Annual surveys also first documented a large C. lutrensis × C. venusta hybrid swarm in the mainstem river in 1998 (Burkhead and Huge 2002; Walters et al. 2008). We collected 1060 Cyprinella between June 2005 and June 2006 from 29 sites spanning a 320-km transect of the mainstem Upper Coosa River system (Fig. 1). Sampling locations included areas to the south of Lake Weiss (the area of introduction), as well as sites representative of the species’ subsequent northward expansion from Lake Weiss and the Coosa River to the Conasauga River. We collected fish from sites spanning approximately 50 km north of the phenotypic invasion front (Table 1; Fig. 1). Distances between sites over the length of the transect were measured from satellite imagery using Google Earth software. At each site, fish were collected by seine with a sampling reach length scaled to 25× river width. All habitats (e.g., pools, runs, and riffles) were sampled to reduce any bias in catch that could arise owing to habitat preference. Once caught, all fish were euthanized and immediately placed in 95% ethanol for morphological and genetic analysis.

Figure 1

Sampling locations along a 320-km transect of the Coosa River, Oostanaula River, and Conasauga River including the introduction point of invasive Cyprinella lutrensis (Lake Weiss) and upstream spread.

Table 1

Summary data for 29 populations of Cyprinella in the Upper Coosa River Basin. Population numbers correspond with the numbers on Fig. 1

Population	Latitude (N)	Longitude (W)	Phenotype n	mtDNA n	Microsatellite n
1	33.997	−85.881	15	35	35
2	34.113	−85.853	36	176	177
3	34.091	−85.744	17	26	26
4	34.165	−85.396	10	11	11
5	34.251	−85.381	5	5	5
6	34.200	−85.256	19	20	20
7	34.255	−85.178	7	106	106
8	34.315	−85.118	123	129	128
9	34.371	−85.125	23	26	26
10	34.411	−85.107	51	64	64
11	34.450	−85.027	12	19	19
12	34.468	−85.033	22	26	26
13	34.494	−85.011	26	26	26
14	34.510	−84.958	34	39	39
15	34.529	−84.966	25	28	28
16	34.573	−84.945	36	45	44
17	34.538	−84.898	14	14	14
18	34.595	−84.928	21	24	24
19	34.667	−84.931	37	43	43
20	34.667	−84.933	8	12	12
21	34.709	−84.868	0	22	22
22	34.736	−84.857	25	25	25
23	34.783	−84.872	23	39	39
24	34.807	−84.859	26	26	26
25	34.818	−84.856	6	17	17
26	34.853	−84.838	18	22	22
27	34.895	−84.829	8	11	11
28	34.904	−84.828	14	14	14
29	34.920	−84.842	9	10	10
Total			670	1060	1059

Because the spatial structure and extent of hybridization in the mainstem Coosa River was unknown at the time of sampling, we also collected specimens from two reference sites where no instances of co-occurrence or hybridization have been recorded (Walters et al. 2008). Information from the two sites served as phenotypic and genetic references against which patterns of differentiation and introgression in the mainstem Coosa River could be assessed. Specimens of C. lutrensis (n = 45) were collected from one site in Proctor Creek (Fulton County, Georgia; 33.795N, −84.475W) and C. venusta (n = 66) were collected from three sites in Terrapin Creek (Cherokee County, Alabama; 34.107N, −85.660W; 34.119N, −85.672W; 34.123N, −85.677W). All focal and reference specimens were characterized in an identical manner.

Quantitative trait measurements

Individuals sampled from the mainstem river transect and reference sites were measured for four phenotypic traits that reliably diagnose parental C. lutrensis and C. venusta (Boschung and Mayden 2004) as follows: standard length (SL), maximum body depth (BD), number of lateral line scales and the intensity of the caudal spot. All measurements were taken on the left side of each specimen, and only specimens of >30 mm SL were measured because of the difficulty of obtaining accurate lateral line scale counts for small individuals. One site (site 21) was excluded from analyses based on quantitative traits because all of the specimens were <30 cm long. These fish were included in genetic analyses.

Specimens were measured for SL from the tip of the snout to the rounded edge of the caudal peduncle, using digital calipers calibrated to 0.01-mm precision, as well as for maximum BD from the anterior junction of the dorsal fin at the dorsal midline to the anterior junction of the anal fin at the ventral midline. We used the ratio of SL to BD as an index of size in subsequent analyses, obtained as residuals for each specimen from the regression of BD on SL. The intensity of the melanic caudal spot was scored on an increasing linear scale of 0–3, with 0 representing the complete absence of the caudal spot and 3 representing maximum intensity. Because color intensity varies among hybrids (Walters et al. 2008), scores of 1 and 2 were included to characterize weak and intermediate spot expression, respectively. The three phenotypic traits (frequency of lateral line scales, size, and caudal spot intensity) demonstrated a high degree of linear correlation across all focal sampling sites (size/caudal spot intensity: Pearson r = 0.63, P < 0.001; size/lateral line scales: r = 0.60, P < 0.001; lateral line scales/caudal spot intensity: r = 0.79, P < 0.001); we therefore used principal component analysis (PCA) to derive an overall phenotypic score for each individual. Trait decomposition yielded a single principal component that explained a high percentage of quantitative phenotypic variation (78%).

We tested for evidence of morphological differences between individuals collected from Proctor Creek and Terrapin Creek, as well as for individuals collected from the terminal transect sites using manova, with site specified as a fixed factor and individual phenotypic traits (frequency of lateral line scales, size, and caudal spot intensity) specified as dependent variables.

Microsatellite genotyping and mtDNA-RFLP assays

We extracted and amplified DNA following Walters et al. (2008). For each individual, genomic DNA was extracted from approximately 0.05 g of preserved fin tissue using DNeasy kits (Qiagen, Valencia, CA, USA). With a few exceptions as noted below, polymerase chain reaction (PCR) mixtures for the amplification of both the complete cytochrome b gene (cyt b) and seven microsatellite loci included 2.5 mm MgCl2, 2.5 mm of each dNTP, 0.5 units Taq DNA polymerase (Invitrogen, Carlsbad, CA, USA), 0.5 μm of the oligonucleotide primers HA and LA (Schmidt et al. 1998), and 0.5 μm PCR buffer (Invitrogen). For some individuals, PCR mixtures for the amplification of cyt b were prepared using 0.1 units Paq (Stratagene, Santa Clara, CA, USA) instead of Taq, and no MgCl2 was included in these mixtures.

Restriction of cyt b with HinfI (New England Biolabs, Ipswich, MA, USA) generates species-specific fragment size profiles in C. venusta and C. lutrensis that can be reliably scored from fluoresced agarose electrophoresis gels. HinfI restriction of cyt b amplified from C. venusta generates 130-, 480-, and 530-bp fragments versus 95-, 130-, 350-, and 570-bp fragments from C. lutrensis. Following Walters et al. (2008), all individuals were assigned species-level mtDNA ancestry from restriction of cyt b amplicons.

Individuals were also genotyped at seven polymorphic microsatellite markers developed for other species that have been shown to distinguish parental C. lutrensis and C. venusta (Walters et al. 2008). Four of the loci (Nme 25C8.208, Nme 18C2.178, Nme 24B6.191, and Nme 24B6.211) were developed for Notropis mekistocholas (Burridge and Gold 2003). Two loci (Rhca20 and Rhca24) were developed for Rhinichthys cataractae (Girard and Angers 2006). The remaining locus (Can6EPA) was developed for Campostoma anomalum (Dimsoski et al. 2000). PCR annealing temperatures were modified according to Walters et al. (2008) with reactions run using fluorescently labeled forward primers. Microsatellite PCR products were subsequently characterized on an MJ Research Basestation Genetic Analyzer and Cartographer software.

Analysis of admixture and genetic differentiation from microsatellite variation

Pairwise estimates of population differentiation between all mainstem transect sites were calculated using the program microsatellite analyser (MSA) (Dieringer and Schlotterer 2003). The significance of pairwise linearized FST values was assessed via permutation tests based on 10 000 runs (Weir and Cockerham 1984). Significance levels were Bonferroni corrected for multiple comparisons.

To provide a baseline for comparisons among samples taken from mainstem transect locations, we estimated pairwise levels of differentiation between the samples taken from Proctor Creek and Terrapin Creek (the reference sites), using the program GENETIX (Belkhir et al. 1996). FST values were estimated from Wier and Cockerham θ values, and 95% confidence intervals (CI) were estimated by bootstrapping 1000 times.

Multilocus admixture profiles for all mainstem and reference site individuals were generated using the program Structure v2.3.1 (Pritchard et al. 2000). In preliminary analyses, we evaluated the relative contribution of individual loci to admixture profiles by comparing the results of runs generated using all seven loci with results of runs involving sequential removal of individual loci in the analysis. No loci were found to bias the results, and all loci were found to be informative. Exploratory runs were also performed with K iteratively set at values ranging from 1 to 8, to determine (i) whether the hybrids represent a distinct evolutionary lineage and (ii) whether parental species exhibit intraspecific genetic structure. No support was found for either condition (data not shown). Therefore, five independent runs with K = 2 (representing the parental species) were executed. For all runs, data were collected over 100 000 iterations, following a 50 000 iteration burn-in, under an admixture model of co-ancestry and correlated allele frequencies (Falush et al. 2003). Mainstem individuals were considered to be pure when the posterior probability of assignment to a parental class averaged over all runs was >90%. For all hybrid individuals, admixture was assessed based on the average assignment values to the first parental cluster following the criteria outlined in the study by Walters et al. (2008).

To assess whether the patterns of introgression were robust to overestimates of hybridization, analyses were also run with parental classifications set at posterior assignment probabilities of q > 0.8 and hybrid classes appropriately adjusted. Under less restrictive probability thresholds designed to minimize the detection of hybrids, the total proportion of hybrids in the transect decreased by only 5%. In addition, we did not detect any significant differences in the patterns of introgression across the transect or major differences in the relative proportions of parental and hybrid individuals across sites. Therefore, we only present classifications and statistical results obtained with parental assignment estimates of q > 0.9.

Clinal analysis of phenotype, mtDNA ancestry, and microsatellite variation

We used the program Cfit (Gay et al. 2008; available at http://www.cefe.cnrs.fr/en/genetique-et-ecologie-evolutive/cfit) to fit clines to multilocus microsatellite genotype, the relative frequencies of mtDNA haplotype assignment, and dominant phenotype according to PCA scores. Cfit fits clines using a simulated annealing algorithm and has the advantage of permitting the simultaneous comparison of multiple quantitative and genetic traits.

In initial analyses, we estimated the parameters c (cline center) and w (cline width, defined as the inverse of the maximum slope) from clines fit to each trait or locus independently. In addition, the distribution of phenotypic variation among mainstem transect locations was explored according to five models (unimodal, bimodal, bimodal with introgression, trimodal and trimodal with introgression; Gay et al. 2008), and the best-fitting distribution was selected using Akaike’s Information Criterion (AIC) (Akaike 1973). We then performed a maximum likelihood search for a common slope (concordance) and center (coincidence) for the traits that demonstrated evidence of clinal structure across the hybrid zone (Bensch et al. 2009). Following Gay et al. (2008), we compared constrained and unconstrained individual clines using likelihood ratio tests with significance assessed at the 5% level.

Results

Reference populations

Multilocus analysis of genetic differentiation between the reference sites recovered a global FST value of 0.24 (95% CI, 0.11–0.39). Single-locus FST values ranged from 0.06 to 0.64, with four loci demonstrating FST values > 0.15. Bayesian analysis of genetic identity based upon a minimum threshold of q > 0.90 classified 100% of reference individuals from Proctor Creek and Terrapin Creek as pure C. lutrensis and C. venusta, respectively (Fig. 2).

Figure 2

Proportion of parental Cyprinella venusta (BS), Cyprinella lutrensis (BS), and hybrid (F1, F2, BCR, and BCB) individuals at 29 sites along a 320-km region of the Coosa River and in two allopatric reference populations (C. venusta and C. lutrensis). For all individuals, admixture categories were determined based on Structure assignment values averaged over multiple runs.

Phenotypic comparison of parental C. venusta (n = 32) with C. lutrensis (n = 20) from the reference populations confirmed that the two species can be readily distinguished on the basis of morphological traits (manova: F3,48 = 504.68, P < 0.001; see also Fig. 4). In particular, C. venusta are characterized by the presence of a prominent black, subtriangular caudal spot that is absent in C. lutrensis (coded character state means ± SE: 2.9 ± 0.06 vs. 0.0 ± 0.0; F1 = 1453.30, P < 0.001), more lateral line scales (mean ± SE: 41.59 ± 0.29 vs. 33.95 ± 0.32; F1 = 293.11, P < 0.001) and a more fusiform body shape (F1 = 71.06, P < 0.001).

Figure 4

Integrative characterization of genetic and phenotypic data for Cyprinella individuals sampled from a 320-km transect of the Coosa River. Phenotype values represent individual principal components analysis (PCA) scores, and distance is shown as the cumulative distance away from the southern terminus of the transect. (A) Hybrids are phenotypically variable and are distributed across the ranges of parental phenotypic variation. (B) Dominant phenotype is strongly associated with mtDNA haplotype, irrespective of hybrid status.

Microsatellite and mtDNA variation across the Upper Coosa River basin

The seven microsatellite loci exhibited a mean of 21 alleles (Table S1) among the 1059 fish genotyped from 29 mainstem transect locations. The multilocus FST value calculated across all individuals and loci provided evidence of weak but significant genetic population differentiation among the sampled locations (FST = 0.05, P < 0.001). Pairwise comparisons indicated that populations sampled at the terminal ends of the transect were significantly differentiated from one another (FST = 0.17, P < 0.001). Expected heterozygosity and allelic richness did not vary consistently across the transect (Table S1).

We detected significant spatial structure in the relative frequencies of pure parental and hybrid multilocus admixture profiles based on Bayesian assignment values. Parental C. lutrensis were numerically dominant at the southernmost end of the transect and decreased in frequency with increasing distance northward with a cline center located 83 km from the southern terminus. No pure parental C. lutrensis were detected at distances >252 km from the southern terminus (Fig. 2). However, hybrid genotypes were recovered at sites reaching more than 50 km beyond the northernmost extent of parental C. lutrensis (Fig. 2). The proportion of individuals containing both C. venusta and C. lutrensis alleles averaged 15% across the transect, with values ranging from 0% to 36% across sites. The majority of hybrid individuals were later-generation hybrids and backcrossed individuals (84%), with backcrosses observed more frequently in the direction of C. venusta (BCR = 15 individuals versus BCB = 76 individuals).

The frequency of the C. venusta mtDNA haplotype increased with increasing northward distance, with the center of the mtDNA cline estimated at 103 km from the transect’s southern terminus (Fig. 3). Approximately 20% of specimens at the southern terminus exhibited a C.venusta mtDNA haplotype, whereas all individuals exhibited a C. venusta mtDNA haplotype at distances >300 km to the north.

Figure 3

Clinal changes in frequencies of Cyprinella venusta phenotype, C. venusta mtDNA haplotypes, and C. venusta microsatellite multilocus genotype along a 320-km transect of the Coosa River. The vertical line represents the consensus cline center. Clinal parameters were estimated using Cfit and visualized using Analyse (Barton and Baird 1998). For graphical purposes, mean population phenotypic values have been scaled between 0 and 1.

Phenotypic variation and cline comparisons across the Upper Coosa River basin

Phenotypic variation exhibited by parental species and hybrids in the Coosa river conformed to a trimodal distribution (Table 2), with the cline center estimated at approximately 105 km upstream from the southern terminus of the transect (Fig. 3). Phenotypic variance was highest in the geographic center of the transect, where fish resembling both parental C. lutrensis or C. venusta co-occurred over a wide region spanning approximately 200 km (sites 4–22; Fig. 4). At the terminal sites of the transect, individuals differed significantly in all four measured phenotypic traits (manova: F3,21 = 24.28, P < 0.001), with fish at the northern terminus exhibiting a more slender body (F1,23 = 54.09, P < 0.001), a greater intensity in the degree of caudal spot melanism (2.80 ± 0.20 vs. 0.14 ± 0.29; F1,23 = 64.92, P < 0.001) and more lateral line scales (41.90 ± 0.99 vs. 34.71 ± 0.62; F1,23 = 54.09, P < 0.001) than those at the southern terminus. A comparison of trait means at the terminal sites with those of reference populations indicated that individuals at the southern and northern ends of the transect were phenotypically indistinguishable from the pure reference C. lutrensis (manova: F3,31 = 1.20; P = 0.33) and C. venusta (F3,38 = 0.28; P = 0.84), respectively.

Table 2

Comparison of five models of the distribution of phenotypic traits in the zone of hybridization between Cyprinella lutrensis and Cyprinella venusta in the Coosa River

	n	Deviance	AIC
Unimodal	7	1699.54	1713.54
Bimodal	6	1378.84	1390.84
Bimodal (with introgression)	8	1411.90	1427.90
Trimodal	12	1327.98	1351.98
Trimodal (with introgression	14	1684.12	1712.12

AIC, Akaike’s Information Criterion.

The number of parameters estimated by each model is given as n. The best model (corresponding with the lowest AIC value) is highlighted in bold.

The clines describing genetic and phenotypic distributions along the mainstem transect exhibited both concordance and discordance (Fig. 3). The phenotypic mtDNA, and microsatellite clines statistically converged on a mean cline center located approximately 93 km from the southern terminus, near the northeastern shore of Lake Weiss (close to site 4). However, the center of the cline describing the distributions of co-dominant multilocus genotypes was situated 20 km closer to the southern terminus of the transect than the centers of the unconstrained phenotypic and mtDNA haplotype clines, and 10 km closer to the southern terminus than the mean cline center based on all three traits (P = 0.06; Table 3). The width of the mtDNA haplotype and microsatellite genotype clines was not statistically distinguishable from a mean cline width of 327 km. However, at a width of 250 km, the phenotypic cline was significantly narrower than the mean cline width based on all three traits (Table 3). Overall best-fit comparisons of the non-constrained model with those constrained for both center and width using likelihood ratio tests indicated that the unconstrained model fit the data better overall than the model constraining center, width, or both (differed by >2 AIC points; Table 3), which further suggests that the clines are discordant and incongruent.

Table 3

Comparison of genetic and phenotypic clines across the Cyprinella hybrid zone between unconstrained models (parameters are allowed to vary freely), and models wherein the cline center (c) or the cline width (w) have been constrained. All measures are expressed in kilometers. Cline centers are expressed as the fluvial distance from the southern terminus of the transect, corresponding to the area of Cyprinella lutrensis introduction

	Unconstrained			Center constrained c = 93		Width constrained w = 327		Both constrained
	Center	Width	Deviance	Deviance	P	Deviance	P	Deviance	P
Phenotype	105	250	1328.00	1329.68	0.19	1335.78	0.01	1336.14	<0.01
mtDNA	103	323	1104.26	1106.28	0.16	1104.28	0.89	1106.84	0.11
Multilocus genotype	83	355	1945.42	1948.74	0.06	1947.94	0.11	1949.28	0.05
n			16	14		14		12
AIC			4409.68	4412.70		4416.00		4416.26

AIC, Akaike’s Information Criterion.

The number of parameters is given by the number n. The best model is given in bold.

Discussion

The propensity of species boundaries to erode can vary according to the nature and extent of interactions between parental species and hybrids following secondary contact (Harrison 1993; Arnold 1997). Hybridization and introgression may disrupt co-adapted gene complexes and can lead to the formation of a hybrid swarm characterized by increased genetic and phenotypic variation (Endler 1977; Jiggins and Mallet 2000). Variation in fitness among hybrid genotypes and the rise of advantageous traits in admixed populations can subsequently hasten the erosion of species boundaries or genetic assimilation of one parental species (Arnold 1997; Barton 2001; Coyne and Orr 2004; Hall et al. 2006).

The introduction of C. lutrensis into the Coosa River has resulted in a biological invasion involving introgressive hybridization (Walters et al. 2008). Most of the hybrids that we collected were later-generation and backcrossed individuals, suggesting that continued upstream expansion of the hybrid swarm (Walters et al. 2008) is unlikely to be restricted by reduced hybrid viability and fertility (Blum et al. 2010). Clinal discordance in the distribution of multilocus microsatellite genotypes, mtDNA haplotypes, and phenotypic traits suggests that the rate of spread within the mainstem Coosa River is more geographically expansive than can been documented from observable patterns of phenotypic variation. In addition, most pure parental and hybrid individuals collected in the southern region of the transect exhibited a C. lutrensis phenotype (Fig. 4). These findings indicate that the rise and retention of a single parental species’ phenotype following an introduction event can mask a history of hybridization and can lead to the appearance of one species being displaced by another without hybridization.

Spatial patterns of genetic differentiation: coincidence and concordance of clines

Loci under different selection regimes are expected to introgress across species boundaries at different rates (Harrison 1990; Mallet 2005; Yuri et al. 2009), with genes under divergent selection introgressing to a lesser extent than neutral or positively selected loci (Gay et al. 2009; Maroja et al. 2009). Consistent with this idea, we found wide clines in multilocus microsatellite admixture profiles and mtDNA haplotypes relative to the cline describing variation at phenotypic traits. This suggests that microsatellite alleles and mtDNA haplotypes are introgressing more readily than phenotypic traits (Endler 1977). Comparison with theoretical values of introgression of neutral traits also suggests that diffusion of alleles, haplotypes, and traits is constrained (Endler 1977). Assuming a 1-year generation time (Lee et al. 1980) and a maximum dispersal distance of 31 km/year (Walters et al. 2008), the width of neutral clines would be approximately 490 km (Endler 1977). As C. lutrensis can exhibit two generations per year, it is possible that neutral clines could be as broad as 700 km (Farringer et al. 1979). As all three clines are narrower than expected assuming free diffusion, there must be factors structuring introgression across the hybrid swarm (Endler 1977; Gay et al. 2007).

Phenotypic model comparisons and examination of individual phenotypic scores indicate that variation across the transect can be partitioned into two modes that broadly resemble one or the other parental species, with variable hybrids intermediate to, and distributed across, both phenotypic parental classes (Fig. 4A). Most lutrensis-like and venusta-like hybrid individuals are associated with the corresponding parental mitochondrial haplotype (Fig. 4B), suggesting that a hybrid’s phenotype is influenced by the maternal contribution (see Garrett 1988; Skulason et al. 1989; Holtmeier 2001 for similar examples in other fishes). This could be giving rise to the coincident centers of the clines describing phenotypic variation and the frequency of mitochondrial haplotypes. Similarly, it could explain why the cline centers are located at the transition from sites composed almost exclusively of parental C. lutrensis and hybrids that are phenotypically indistinguishable from C. lutrensis, to sites encompassing elevated phenotypic variability across the transect.

Maternal effects can influence the direction and rate of introgression if maternal phenotype confers a fitness advantage to hybrid progeny through competitive superiority, increased mate success, or both (Mousseau and Fox 1998; Burgess and Husband 2004).The patterns of clinal concordance and discordance that we have observed suggest that selection may be maintaining some degree of phenotypic species integrity despite extensive admixture of other genomic regions (Bensch et al. 1999; Balloux et al. 2000; Rees et al. 2003; Vallender et al. 2007; Brelsford and Irwin 2009). The relative abundance of invasive parental C. lutrensis and C. lutrensis-like hybrids, and the corresponding decline or absence of parental C. venusta or C. venusta-like hybrids in the southern region of the transect, which is the initial area of introduction and the region with the longest history of hybridization, suggests that C. lutrensis traits may be selectively favored. Parental C. lutrensis can overcome native congeners following introduction because of an aggressive disposition, short generation times, and high fecundity (DeVivo 1995; Fuller et al. 1999; Burkhead and Huge 2002). Hybrids that inherit these characteristics would also likely become comparatively more abundant.

The retention of the non-native C. lutrensis parental phenotype in admixed populations could also be promoted by hybrids preferentially mating with individuals most closely resembling maternal parental taxa. Shape and color pattern often function as cues for species recognition (Ptacek 2000) and likely also influence patterns of assortative mate choice in shiners (Blum et al. 2010). The rise and retention of the C. lutrensis phenotype following an introduction event might therefore reflect both prezygotic and postzygotic reproductive advantages, where the benefits of increased fecundity and short generation time are buoyed by preferences for parental traits (Lande 1981; West-Eberhard 1983; Panhuis et al. 2001; Turelli et al. 2001).

Reproductive isolation and hybrid fitness

Patterns of association between mtDNA haplotypes and nuclear genotypes suggest that the direction of hybridization is roughly symmetrical between C. lutrensis and C. venusta. Of the 29 F1 hybrids recovered from mainstem reaches, 55% had a C. lutrensis haplotype and 45% had a C. venusta haplotype. This result contrasts with laboratory-based experiments that found initial hybridization events between female C. venusta and male C. lutrensis occur more frequently than the corresponding reciprocal crosses (Blum et al. 2010). However, in partially isolated groups such as Cyprinella minnows, the rarity of conspecific mates following an introduction event can lead to an increased frequency of heterospecific mating between introduced and more abundant native species (Hubbs 1955; Currat et al. 2008; Lepais et al. 2009; Blum et al. 2010). Most of the hybrids recovered in the mainstem Coosa River are later generation and backcrosses predominantly in the direction of native C. venusta. This finding also suggests that postzygotic barriers to gene flow between the species are weak, which is consistent with controlled laboratory crosses, demonstrating that C. venusta × C. lutrensis F1 and backcross hybrids are capable of reproducing (Hubbs and Strawn 1956) and that by some measures (e.g., larval mortality), hybrid fitness exceeds that of parental species (Blum et al. 2010).

Hybrid success has been identified as a key factor facilitating species erosion via introgressive hybridization (Perry et al. 2001; Behm et al. 2010). For example, the collapse of ecologically divergent limnetic and benthic sticklebacks in Enos Lake (British Columbia, Canada) into a single intermediate hybrid swarm is associated with a shift from diet specialization to generalist feeding tendencies. Intermediate-shaped fish show no decrease in growth or survivorship compared with parental species, implying that the loss of ecologically based postmating isolation has contributed to the collapse of the species pair (Taylor et al. 2006; Behm et al. 2010). Similarly, it is likely that the upstream expansion of the C. lutrensis × C. venusta hybrid swarm has been facilitated by yearly increases in the frequencies of successfully reproducing later-generation hybrids (Endler 1977; Huxel 1999; Hall et al. 2006).

Introgressive hybridization and the success of invasive species

Biological invasion can lead to native species decline via displacement by competitively superior introduced species or through introgressive hybridization leading to genetic homogenization or assimilation (Levin et al. 1996; Wolf et al. 2001). Native and introduced species remain distinct during species displacement, where spatial distributions of phenotypic and molecular traits clearly reflect the geographic transition from one parental species to another. In contrast, biological invasion facilitated by hybridization leads to high frequencies of hybrid genotypes and introgression of phenotypic and molecular traits. In the Coosa River, the frequency of parental and hybrid individuals with phenotypes closely resembling invasive C. lutrensis has increased in the wake of the expanding hybrid swarm, resulting in the appearance of one species being displaced by another without hybridization.

In the absence of selection for or against hybrids, the decline of pure native species in response to hybridization is largely influenced by demographic factors including the relative of abundances hybridizing groups and the dispersal potential of the colonizing species (Endler 1977; Lepais et al. 2009). The establishment of a hybrid swarm can facilitate invasion success by reducing the frequency of pure native parental individuals and increasing the frequency of non-native alleles in admixed progeny (Levin et al. 1996; Wolf et al. 2001). Introgression can accelerate the rate that native alleles or traits are lost if non-native genotypes have a fitness advantage (Huxel 1999; Hall et al. 2006). Physiological factors such as higher rates of growth and reproduction in invasive species (Sousa et al. 2008; Alonso and Castro-Diez 2008) and heritable behavioral attributes such as generalist feeding tendencies, high aggressiveness, and successful predator avoidance (Pennuto and Keppler 2008; Balshine et al. 2005) can all increase the likelihood that non-native alleles or traits will be maintained or favored, leading to the decline of native species despite often initially large asymmetries in relative abundance (Mooney and Cleland 2001).

The expanding C. lutrensis × C. venusta hybrid swarm in the Coosa River serves as an example of how the geographic rate and extent of invasions can be underestimated if surveys are limited to observable patterns of phenotypic variation. Distributional records and other observations solely based on phenotype may not necessarily disclose the nature of species interactions following secondary contact. Our findings indicate that the rate of spread within the mainstem Coosa River has likely been more rapid and geographically more expansive than has been documented from observable patterns of phenotypic variation. In addition, the rise and retention of the C. lutrensis phenotype in the southern region of the transect superficially creates the appearance of native species displacement without any hybridization. Because the loss of native species has often been characterized as a process of displacement following the introduction and spread of a non-native species, it is possible that other instances of native species decline involving introgressive hybridization are being mistaken for species displacement. Re-evaluating other well-studied invasions (resulting from introductions or range expansions) might therefore be necessary to better understand how ecological and evolutionary factors promote the loss of native populations following secondary contact.