Variation in habitat across the landscape can limit dispersal between populations and influence the diversity and degree of genetic structure of populations. This effect may be larger in small aquatic species with low vagility, specific habitat requirements, and that inhabit lotic systems where dispersal may be biased by the unidirectional movement of water downstream. Darters offer a unique perspective to understand these influences because they have variable life histories and habitat requirements, and evidence suggests that certain life history characters may affect gene flow. Data were collected from a combination of microsatellite loci from 825 individuals to compare genetic structure and differentiation of three species of darters (Ammocrypta beanii, Etheostoma swaini, and Percina nigrofasciata) with differences in life histories and habitat requirements in the Pearl River basin. Our analyses reveal low, but interesting, patterns of population subdivision. Genetic subdivision among populations was highest for E. swaini, intermediate for P. nigrofasciata, and lowest for A. beanii. This pattern suggests a difference in dispersal, influenced by the variation in available habitat required for each species. Further investigations are warranted for small, benthic, and especially threatened species in the Pearl River basin as the main channel has undergone a significant shift in geomorphology and habitat diversity since the construction of reservoirs and low-head dams along its course.
ASSESSING population genetic structure through space and time can provide insight into the factors that impact contemporary populations (Templeton, 2006). Restrictions to gene flow among local populations can lead to subdivision of a population into genetically distinct groups. Such restrictions may be caused by physical barriers, life history characteristics, or by habitat preferences that limit or bias the dispersal ability of a species (Frankham et al., 2010; Alp et al., 2012; Allendorf et al., 2013; Caputo et al., 2014; Thomaz et al., 2016). Genetic subdivisions may be more pronounced in aquatic systems, because unlike terrestrial environments, the dendritic arrangement of fluvial habitats restrict the movement of individual organisms in a complex and asymmetrical pattern (Meffe and Vrijenhoek, 1988; Gornall et al., 1998; Bilton et al., 2001; Fagan, 2002; Morrissey and Kerckhove, 2009; Hutsemékers et al., 2013; Paz-Vinas et al., 2015).
Patterns of dispersal and gene flow among obligate freshwater species inhabiting lotic systems may be biased by the unidirectional movement of water downstream (Congdon, 1995; Fagan, 2002; Miller et al., 2002; Morrissey and Kerckhove, 2009; Pollux et al., 2009). Aquatic species also may be susceptible to variation in habitat across the landscape, which can also limit dispersal between populations and influence the diversity and degree of genetic structure of populations (Piller et al., 2004; Nikula et al., 2011; Wang and Bradburd, 2014; Waters et al., 2015; Waters and Burridge, 2016). This effect may be more pronounced in small and/or non-migratory aquatic species with specific habitat requirements (Tibbets and Dowling, 1996; Vanormelingen et al., 2015; Waters et al., 2015; Hintz et al., 2016; Stoll et al., 2016; Waters and Burridge, 2016).
Darters (Teleostomi: Percidae) are a taxonomically rich group of morphologically and ecologically diverse species (Page, 1983; Ross, 2001). Most species occupy lotic systems and have strict habitat requirements, including preferences for a particular type of benthic substrate or flow regime, which make them particularly sensitive to fragmentation and habitat degradation (Page, 1983; Karr et al., 1986; Tipton et al., 2004; Fluker et al., 2014). Inherent variation in intrinsic traits (i.e., habitat preference) coupled with extrinsic habitat availability can also influence contemporary gene flow and population genetic structure in darters (Roberts et al., 2013). Therefore, darters are an ideal candidate to understand how environmental variation (i.e., habitat differences) can influence dispersal and genetic structure of closely related sympatric species in a fluvial system.
We investigated the contemporary population genetic structure of three sympatric etheostomatine darters (Etheostoma swaini, Percina nigrofasciata, and Ammocrypta beanii) with three different habitat preferences in the Pearl River basin (LA/MS; Figs. 1, 2). The Gulf Darter, Etheostoma swaini, primarily occurs in small- to medium-sized creeks and rivers and occasionally in the main channel of the Pearl River (Piller et al., 2004; Geheber and Piller, 2012), in habitats comprised of mud, sand, gravel, and vegetation (Ruple et al., 1984; Ross, 2001). Based on several life history accounts across its range (Mathur, 1973; Mayden, 1992; Ross, 2001; Hughey et al., 2012), the Blackbanded Darter, P. nigrofasciata, can be found in small- to medium-sized rivers and streams, and movement distances of adults is considered to be low to moderate (Freeman, 1995). Its abundance is positively correlated with amount of instream cover (Rakocinski, 1988). The Naked Sand Darter, Ammocrypta beanii, is widespread throughout the basin and primarily occurs in the main channel and larger tributaries of the Pearl River basin (Ross, 2001). It is found almost exclusively over sandy substrate (Heins and Rooks, 1984; Ross, 2001) and often in large numbers. Although all three species occur in the main channel of the Pearl River, A. beanii occurs in the greatest abundance relative to the other two species which are much less common in the main channel (Piller et al., 2004; Geheber and Piller, 2012; K. Piller, pers. obs.).
Our objective for this study is to understand the degree of genetic connectivity of three sympatric species of darters with different habitat preferences using contemporary genetic data. The hypothesized relationships found in previous studies between gene flow and certain intrinsic and extrinsic factors in other aquatic species allow us to make predictions in regard to gene flow and genetic structure (Caputo et al., 2014; Waters et al., 2015; Waters and Burridge, 2016). It may be predicted that the darter species in this study will show varying levels of gene flow and genetic structure as a result of their variation in preferential habitats. We hypothesize that A. beanii will display high levels of gene flow and minimal genetic structure due to its preference for medium- to large-sized sandy bottom streams like the Pearl River. Comparatively, we expect that P. nigrofasciata and E. swaini will show lower levels of gene flow and a greater degree of genetic structure relative to A. beanii because of their occupancy of small- to medium-sized streams with woody structure or vegetation. The main channel is not the preferred habitat of these species.
MATERIALS AND METHODS
Study area and specimen collection.—The Pearl River is a Gulf Coastal Plain system that forms the southeastern border of Louisiana and Mississippi, flowing in a southwesterly direction from central Mississippi to the Mississippi Sound (Gulf of Mexico). It has been subjected to several major anthropogenic modifications, including the construction of the large Ross Barnett Reservoir and two low-head dams (the Bogue Chitto and Pools Bluff sills). Unfortunately, prominent changes in fish community and habitat have been amplified by human-induced perturbations caused by channel instability, increased erosion, and severe bank destabilization in the Pearl River (Piller et al., 2004; Tipton et al., 2004; Geheber and Piller, 2012). Species of darters included in this study commonly occur in the Pearl River basin. Specimens of each species were collected in the Pearl River and its tributaries, including the Bogue Chitto River and Hobolochitto Creek, from March to June 2011 (Fig. 1). A total of nine sites were sampled for each species, with the exception of E. swaini in which collections occurred only at eight total sites (see Tables S1–S3 in Supplemental Material for sampling locality details, associated population codes, sample sizes, and museum catalogue numbers; see Data Accessibility). Populations were coded based on the associated river system or tributary (BC = Bogue Chitto, UP = upper Pearl River, or LP = lower Pearl River, taken from the Hobolochitto creeks) and from the relative location sampled (3 = upper, 2 = middle, and 1 = lower). Lower Pearl River populations were designated with a ‘w' or ‘e' if the population was sampled from West or East Hobolochitto Creek, respectively. Total DNA was extracted using the Wizard® SV 96 genomic DNA purification system and Vac-Man® 96 vacuum manifold kit, following the manufacturer's protocol (Promega Inc., Madison, WI; https://www.promega.com/products/dna-purification-quantitation/genomic-dna-purification/wizard-sv-96-genomic-dna-purification-system/). Voucher tissues and specimens were deposited in the Southeastern Louisiana University Vertebrate Museum.
Microsatellite amplification and genotyping.—A total of 26 published microsatellite loci (Table S4; see Data Accessibility) were individually screened for each species using primers designed for Ammocrypta pellucida (Ginson, 2012; Ginson et al., 2015), Etheostoma blennioides (Beneteau et al., 2007), E. caeruleum (Tonnis, 2006), E. chermocki (Khudamrongsawat et al., 2007), E. osburni (Switzer et al., 2008), E. scotti (Gabel et al., 2008), and Percina rex (Dutton et al., 2008). Each marker was reviewed for quality of amplification and polymorphism. Loci that failed to amplify or were monomorphic were thrown out of subsequent analyses. Ultimately, eight loci were used for A. beanii, 11 loci for E. swaini, and seven loci for P. nigrofasciata. Loci were individually amplified and then combined post-PCR into two groups (A. beanii, P. nigrofasciata) or three groups (E. swaini) before fragment analysis (Table S5; see Data Accessibility). Individual PCR amplifications were carried out in a reaction volume of 12.5 μl each containing 1.0 μl DNA extract (10–50 ng); dNTPs (0.2 mM each); 1.2 μl of primer pair mixture (2 μM each); MgCl2 (2 mM); Taq DNA Polymerase (0.06 U/μl); 1.25 μl of 10X Standard Taq Buffer; and 7.65 μl of distilled water. Amplification protocols for all loci were as follows: initial denaturation of DNA at 95°C for 30 s, followed by 30–35 cycles of denaturation at 94°C for 30 s, annealing at 57–60°C for 1.5 min, extension at 72°C for 1 min; and a final extension at 72°C for 10 min. PCR products were submitted to the DNA Analysis Facility on Science Hill at Yale University Sequencing Center for genotyping using a 3730xl 96-Capillary Genetic Analyzer (Thermo Fisher Scientific, Inc., Waltham, MA; https://www.thermofisher.com/order/catalog/product/3730XL). Alleles were scored using the program GeneMarker® v1.3 (SoftGenetics, LLC, State College, PA; http://www.softgenetics.com/GeneMarker.php).
Microsatellite characters and genetic diversity estimates.—Conformance to Hardy-Weinberg equilibrium (HWE) was evaluated for each locus in each sample population using ARLEQUIN v3.5.2.2 (Excoffier and Lischer, 2010). Significance was assessed using exact tests with 1,000,000 steps in Markov chain and 100,000 dememorization steps at a significance level of 0.05. Tests for linkage disequilibrium (LD) between loci in each population were performed in ARLEQUIN v3.5.2.2 (Excoffier and Lischer, 2010). Significance was assessed at the 0.05 level using 99,999 permutations and an initial condition set at 5. A false discovery rate (FDR) of 0.10 was applied to HWE and LD tests in order to correct for multiple testing (Benjamini and Hochberg, 1995). To detect genotyping errors and to screen for null alleles, species genotypic matrices were analyzed by population with the program MICROCHECKER v2.2.3 (Van Oosterhout et al., 2004). Loci that consistently failed tests for HWE, LD, and/or null alleles were removed from subsequent analyses. Unbiased gene diversity (expected heterozygosity, HE), observed heterozygosity (HO), allelic frequencies, richness, and number were calculated for each population across loci (Table 1) and for each locus across populations (Table S6; see Data Accessibility) using ARLEQUIN v3.5.2.2 (Excoffier and Lischer, 2010) and the R package PopGenReport v.2.2.3 (Adamack and Gruber, 2014).
Table 1.
Observed genetic diversity at each subpopulation across all loci of each darter species in the Pearl River basin. Subpopulations include Bogue Chitto River (BC3, BC2, BC1), upper Pearl River (UP3, UP2, UP1), and lower Pearl River (LP3w, LP2w, LP3e, LP2e, or LP1e). N is total number of individuals sampled; NT is total number of alleles; NA is mean number of alleles; AR is mean allelic richness; HE is expected heterozygosity; HO is observed heterozygosity; NL is total number of microsatellite loci
Population genetic structure and differentiation.—To estimate levels of population differentiation, the program FSTAT v2.9.3.2 (Goudet, 2002) was used to calculate pairwise FST estimator θ (Weir and Cockerham, 1984) between populations, assuming Hardy-Weinberg equilibrium based on aforementioned Hardy-Weinberg conformance tests (Tables 2–4). Pairwise significance was assessed using 1000 permutations, and P values were adjusted after standard Bonferroni corrections for multiple testing (Rice, 1989). Pairwise comparisons were performed between each sample site in the Pearl River system. Other pairwise statistics, including a separate estimate of pairwise FST, and associated bias corrected 95% confidence intervals were calculated using the diffCalc function in the R package diveRsity (Keenan et al., 2013). The program STRUCTURE v.2.3.4 (Pritchard et al., 2000) was used to assess population genetic structure without a priori group designations using a Bayesian clustering method (Fig. 3). STRUCTURE determines the number of distinct genetic clusters (K) among samples and probabilistically assigns individuals to a specific cluster or multiple clusters if there is admixture (Pritchard et al., 2000; Falush et al., 2003). Because different runs may produce variation in likelihood values, the number of MCMC repetitions to be used after burn-in for each analysis was determined by multiple runs at different values of K (number of population clusters) to ensure consistent results between runs. Once sufficient parameters were determined, final analyses of 20 independent runs for each value of K (100,000 burnin; 300,000 after burnin) were performed under a model allowing correlated allele frequencies between populations and admixture of genotypes (Falush et al., 2003). Total values of K considered for each analysis was equivalent to total number of sample sites plus 3 in order to encompass all sample sites as potential populations and to include the possibility of intra-population structure (Evanno et al., 2005). Appropriate values of K were estimated using the ad hoc statistic ΔK (Evanno et al., 2005) through the web-based program STRUCTURE HARVESTER v0.6.93 (Earl and vonHoldt, 2012). Estimated cluster membership coefficient matrices of all 20 runs at each value of K was further analyzed in the program CLUMPP v1.1.2 (Jakobsson and Rosenberg, 2007) in order to find optimal alignments of R (20) replicated clusters at a given K by employing the full search method. The mean of the permuted matrices across replicates was used to create graphical displays the R-based package ggplot2 v.2.2.1 (Wickham, 2009). The ΔK method can obscure substructure in the data, particularly if there is stronger subdivision between groups. Substructure was suspected within the Pearl River groups for E. swaini (see Results section). A hierarchical structure analysis was performed for this species, removing Bogue Chitto River data in order to visualize underlying substructure within Pearl River populations.
Table 2.
Pairwise estimates of FST (θ; Weir and Cockerham, 1984) between populations of Ammocrypta beanii. Pairwise significance values were adjusted after Bonferroni corrections. Significant values are denoted with asterisks.
Table 3.
Pairwise estimates of FST (θ; Weir and Cockerham, 1984) between populations of Etheostoma swaini. Pairwise significance values were adjusted after Bonferroni corrections. Significant values are denoted with asterisks.
Table 4.
Pairwise estimates of FST (θ; Weir and Cockerham, 1984) between populations of Percina nigrofasciata. Pairwise significance values were adjusted after Bonferroni corrections. Significant values are denoted with asterisks.
To further test for genetic differentiation, an analysis of molecular variance (AMOVA) was utilized in the program ARLEQUIN v3.5.2.2 (Excoffier and Lischer, 2010). An AMOVA is related to pairwise FST estimates, but instead of comparing allele frequencies between populations, an AMOVA examines the amount of mutational differences among loci. Unlike STRUCTURE, analyses of molecular variance allows genetic variation to be partitioned into groups that are designated a priori. For each species, populations were categorized into three groups: Bogue Chitto River (BC1-3), Upper Pearl River (UP1-3), and Lower Pearl River (LP1-3 east or west, depending on the species). Groupings were chosen to test if differences between systems (groups) were driving genetic differences. A hierarchical AMOVA was run for each species using 10,000 permutations, and variance fixation indices were calculated among groups, among populations within groups, among individuals within populations, and within individuals (Table 5). Results of the analysis are based on averages over all loci within groups for each species.
Table 5.
Global AMOVA results averaged over all loci for each species. Local populations were categorized into three groups: Bogue Chitto River, Upper Pearl River, and Lower Pearl River. Significant P values are denoted by asterisks. (Va, FCT): among groups; (Vb, FSC): among populations within groups; (Vc, FIS): among individuals within populations; (Vd, FIT): within individuals.
Isolation-by-distance.—In order to examine a possible isolation-by-distance (IBD) pattern, the R packages adegenet v2.1.0 (Jombart, 2008; Jombart and Ahmed, 2011) and lmodel2 v1.7-2 (Legendre, 2014) were utilized (Fig. 4; Table 6). Mantel tests (Mantel, 1967) were performed with 100,000 randomizations to estimate correlations between matrices of genetic distance (θ) estimates and geographical distances (river km) between pairwise populations and to calculate associated P values. All geographical distances for each species were log-transformed. A reduced major axis regression analysis was performed to calculate the intercept, slope, and R2 value of the genetic distance versus log-transformed geographic distances.
Table 6.
Results from the isolation by distance analysis for each species. Correlation coefficients and associated P values were assessed using a Mantel test. Linear model parameters (intercept and slope), associated 95% confidence intervals, and coefficients of determination (R2) were estimated using a reduced major axis regression method. Significant correlations are designated by an asterisk.
RESULTS
Ammocrypta beanii.—A total of 287 tissue samples were collected, averaging about 32 individuals per site (Table 1; Table S1). From a total of 72 tests for Hardy-Weinberg equilibrium, 39 (54%) showed significant deviations from HWE after corrections. However, 27 of these 39 tests (69.2%) involved loci EosD107, EosD10, and EosC112, and these loci were significant across all populations. Other loci were significant at only some of the populations, showing no consistent patterns. Out of 252 tests for linkage disequilibrium, 85 (34%) were significant. Of those 85 tests, the majority (64; 74%) were significant at loci EosD107 and EosD10. Tests using MICROCHECKER were consistent with HWE deviation and linkage disequilibrium tests. Based on these results, three loci (EosD107, EosD10, and EosC112) were removed from subsequent analyses. Ammocrypta beanii had a global mean number of alleles of 20.82±0.95 and a global mean allelic richness of 18.27±0.77 (Table 1; Table S6; see Data Accessibility). The global mean observed heterozygosity (HO) was 0.854±0.013 (HE = 0.922±0.007).
Pairwise FST estimates showed little genetic differentiation overall for A. beanii among Pearl River populations (Table 2). The only significant estimates occurred between BC3 and BC1 (FST = 0.012), BC1 and LP2w (FST = 0.006), and BC1 and LP1e (FST = 0.012). For the STRUCTURE analysis, the predicted number of population clusters for A. beanii in the Pearl River was K = 5, but individuals across all populations were assigned to either cluster almost equally suggesting no population structure (Fig. 3A). Additionally, the ΔK estimator cannot estimate populations with K = 1 (Evanno et al., 2005). When the mean log likelihood values were compared, K = 1 was much larger than any other estimates of K.
Results of the AMOVA reveal all variance components were significant, with the exception of within groups among populations (Vb, FSC; Table 5). The percentage of variation was lowest for between groups and within groups among populations estimates (0.201% and 0.191%, respectively). Individual level components had the largest percentage of variation (about 7.4% among individuals within populations and 92.2% within individuals). For measurements of isolation by distance, the correlation between genetic distances and log-transformed geographic distances was low among populations of A. beanii (r = 0.010), and this association was not significant (P = 0.477, α = 0.05), indicating that the data did not fit a model of isolation by distance (Table 6; Fig. 4A). The model only explained 0.009% of the genetic variation within this species (R2 = 9.39e−5).
Etheostoma swaini.—A total of 251 tissue samples were collected, averaging about 31 individuals per site (Table 1; Table S2; see Data Accessibility). Just 20 out of 88 tests (22%) for Hardy-Weinberg equilibrium were significant, but only locus EosD107 was significant across all populations. Likewise, only 10 out of 440 (2%) tests for linkage disequilibrium were significant, and 4 out of the 10 tests (40%) involved locus EosD107. There was no pattern of LD or HWE deviations among other loci. Tests using MICROCHECKER were consistent with HWE deviation and linkage disequilibrium tests. Based on these results, locus Eos107 were removed from subsequent analyses. The global mean number of alleles was 13.68±1.02 and the global mean allelic richness was 11.80±0.85 (Table 1; Table S6; see Data Accessibility). Global mean HO was 0.641±0.037 (HE = 0.688±0.039).
Etheostoma swaini showed significant levels of differentiation among populations, with significant pairwise FST values ranging from 0.016 to 0.069 (Table 3). Only populations within the Bogue Chitto River (BC 1-3) showed non-significant pairwise estimates, suggesting these populations showed little genetic differentiation relative to other populations within the Pearl River system. STRUCTURE predicted the number of population clusters was only K = 2 using the ΔK estimator. Individuals from sites in the Bogue Chitto River (BC 1-3) were primarily assigned to a single cluster, and individuals from the upper and lower Pearl River populations were primarily to assigned to another cluster (Fig. 3B). However, all individuals across all populations showed indications of some admixture. It should be noted that despite an estimation of K = 2 clusters, the mean estimated log likelihoods of K = 2 to 5 were relatively similar. This may be indicative of masked substructure within the data, particularly within the Pearl River groups (see methods section). Once the Bogue Chitto populations were removed, there was still no indication of substructure within the Pearl River groups after multiple tests (data not included). The STRUCTURE results corroborate results from pairwise FST estimates (Table 3), which indicated a much higher level of genetic differentiation between groups (global average FST = 0.047) and the highest differentiation between Bogue Chitto and Pearl River populations (mean FST = 0.054). When comparing pairwise estimates within groups, values were low for Pearl River groups (mean FST = 0.027) and were very low for Bogue Chitto groups (mean FST = 0.005).
All components in the AMOVA were significant for E. swaini (Table 5). The lowest percentages of variation occurred among groups and among populations within groups (about 3.1% and 1.7%, respectively). The highest percentages of variation occurred among individuals within populations and within individuals (about 6.5% and 88.7%, respectively). Isolation by distance measures for populations of E. swaini revealed significant correlation between genetic and geographic distances (r = 0.660, P = 0.004), indicating that these data may fit an isolation by distance model (Table 6; Fig. 4B). However, the isolation by distance model only explained 43.6% of the genetic variation (R2 = 0.436).
Percina nigrofasciata.—A total of 287 tissue samples were collected, averaging about 32 individuals per site (Table 1; Table S3; see Data Accessibility). Only six out of 63 tests (9.5%) deviated from Hardy-Weinberg expectations, with significant deviations occurring only with locus Ebl4. Only 8 out of 189 tests (4%) for linkage disequilibrium were significant. There were no consistent patterns of deviations across populations. Tests using MICROCHECKER were consistent with HWE deviation and linkage disequilibrium tests. Results were consistent when subsequent analyses were run with and without the locus Ebl4. Therefore, results presented here represent analyses without removing loci. For this species, the global mean number of alleles was 18.86±0.44, with a global mean allelic richness of 16.66±0.34. The global mean HO was 0.883±0.011 (HE = 0.933±0.002).
Percina nigrofasciata showed variable levels of genetic differentiation, with significant values ranging from FST = 0.003 to 0.022 (Table 4). Interestingly, there was no significant differentiation among populations within the upper Pearl River (UP 1-3), among populations within the Bogue Chitto River (BC 1-3), and among populations within the lower Pearl River (LP3w, LP2w, LP1e). Likewise, populations in the Bogue Chitto River system seemed to be more genetically similar to lower Pearl River populations than upper Pearl River. Similar to A. beanii, predicted population clusters from STRUCTURE for P. nigrofasciata was K = 4 and all individuals assigned to either cluster (Fig. 3C). Although individuals from each sample site were assigned to all four clusters with varying probabilities, individuals from the upper Pearl River (UP 1-3) showed intermediate structure with a higher percentage of individuals being assigned to a single cluster. These results corroborate pairwise FST estimates.
All variance components of the AMOVA for P. nigrofasciata were significant at the α = 0.05 level (Table 5). Similar to the other species, the lowest percentage of variation occurred among groups (0.6%) and among populations within groups (0.3%), while the highest percentages occurred among individuals within populations (5.5%) and within individuals (93.7%). Like E. swaini, isolation by distance estimates for populations of P. nigrofasciata showed a significant correlation between genetic and geographic distances (r = 0.439, P = 0.003), indicating that these data may fit an isolation by distance model (Table 6; Fig. 4C). However, the isolation by distance model only explained 19.3% of the genetic variation (R2 = 0.193).
Species comparisons.—Ammocrypta beanii had the largest mean number of alleles across loci, ranging from 19.60±2.9 (LP2e) to 21.8±3.3 (UP3), despite the lowest mean total number of alleles between populations (104.1±1.3) and the smallest, viable number of microsatellite loci included in the study (NL = 5; Table 1). Etheostoma swaini had the smallest mean number of alleles between populations ranging from 11.20±2.5 (UP3) to 15.50±3.2 (BC2), with a mean total number of alleles of 136.75±5.2. For P. nigrofasciata, the mean number of alleles ranged from 17.86±1.1 (BC2) to 20.57±1.5 (LP2w), with a mean total number of alleles of 132.0±2.0. Similarly, mean allelic richness ranged from 17.54±2.4 (LP2e) to 19.35±3.0 (UP1) for A. beanii, from 15.90±0.9 (BC2) to 17.97±1.1 (LP2w) for P. nigrofasciata, and from 9.83±2.1 (UP3) to 13.14±2.6 (BC2) for E. swaini. Global mean observed heterozygosity did not significantly differ from expected heterozygosity for any species once suspect loci were removed (α = 0.05 level). Global mean observed heterozygosity were similar among populations of A. beanii (HO = 0.854±0.013; HE = 0.922±0.007) and P. nigrofasciata (HO = 0.883±0.011; HE = 0.933±0.002). Comparatively, E. swaini had a much smaller mean observed heterozygosity (HO = 0.641±0.037; HE = 0.688±0.039).
Pairwise FST estimates showed varying degrees of significant differentiation for each species (Tables 2–4). Pairwise FST estimates ranged from FST = 0.000 to 0.012 for A. beanii, FST = 0.003 to 0.069 for E. swaini, and FST = 0.000 to 0.022 for P. nigrofasciata. Despite relatively low pairwise FST estimates across species, estimates for E. swaini were about an order of magnitude larger compared to the other two species. The STRUCTURE results were similar to pairwise estimates (Fig. 3). Ammocrypta beanii and P. nigrofasciata had low signals of population genetic structure, indicating significant admixture among populations for both species. Populations of Etheostoma swaini showed a stronger signal of genetic structuring, particularly between Bogue Chitto populations and Pearl River populations.
In the AMOVA, all components for each species were significant, with the exception of tests among populations within groups for A. beanii (Table 5). Percentage of genetic variation among groups was highest for E. swaini (about 3.1% variation for E. swaini, 0.6% for P. nigrofasciata, and 0.2% for A. beanii), but this component was fairly low relative to other components. Therefore, genetic divergence between major groups was not that high among species, but significantly higher in E. swaini relative to other species. Etheostoma swaini also had higher genetic variation among populations within groups relative to other species (about 1.7% variation for E. swaini, 0.3% for P. nigrofasciata, and 0.2% for A. beanii). Although minimal, these results reflect similar findings in the pairwise FST estimations and the STRUCTURE analysis. The majority of genetic variation occurred within individuals among populations (about 99.6% for A. beanii, 99.1% for P. nigrofasciata, and 95.2% for E. swaini).
Populations of Ammocrypta beanii did not fit a model of isolation by distance (r = 0.010, P = 0.477; Table 6; Fig. 4). In contrast, populations of E. swaini and P. nigrofasciata showed a significant correlation between genetic and geographic distances (E. swaini, r = 0.660, P = 0.004; P. nigrofasciata, r = 0.439, P = 0.003), indicating that these data may fit an isolation by distance model. However, the isolation by distance model only explained 43.9% and 19.3% of the genetic variation for E. swaini and P. nigrofasciata, respectively. This suggests that something other than distance may contribute to genetic diversity among populations of both species, particularly for P. nigrofasciata.
DISCUSSION
The aquatic environment has long been viewed as homogenous with respect to the dispersal, migration, and concomitant levels of gene flow for many organisms. However, a multitude of studies have noted that population subdivision in riverine habitats can be complex as a result of numerous historical and contemporary factors such as vicariance events, dynamic hydrological patterns, and species extinction and re-colonization, along with ecological and biological variation among species' dispersal capabilities, life history differences, and habitat preferences (Meffe and Vrijenhoek, 1988; Costello et al., 2003; Whiteley et al., 2004; Zickovich and Bohonak, 2007). Our samples were collected above and below low-head impoundments (Bogue Chitto Sill and Pools Bluff Sill). Our data indicate that the sills, particularly Bogue Chitto Sill, may inhibit gene flow, especially between Bogue Chitto and Pearl River populations of E. swaini. However, we believe there is insufficient evidence of direct impediment to gene flow across the dams for any of the species included in our study, but may have contributed indirectly by facilitating geomorphological changes. The microsatellite data clearly showed divergent patterns of contemporary genetic structuring among the three species, with Ammocrypta beanii displaying the highest level of gene flow between populations and no genetic structuring, Etheostoma swaini showing the lowest level of gene flow between populations and the greatest amount of structuring, and Percina nigrofasciata being generally intermediate between these two species in both genetic structure and levels of gene flow. The results from our study suggest that habitat preferences of our focal species may have played a major role in dictating the levels of gene flow and genetic differentiation in the Pearl River basin, thereby supporting our original predictions.
Ammocrypta beanii, like all other species in the genus, is associated with sandy bottomed habitats of medium- to larger-sized rivers in which they bury themselves for predator avoidance, energy conservation, and enhancement of prey-capturing abilities (Daniels, 1989; Ross, 2001). Based on historic surveys of the main channel Pearl River over the last 29 years, A. beanii was one of the most abundant darters collected (Geheber and Piller, 2012; Piller and Geheber, 2015); however, prior to 1988, A. beanii was only the third most abundant darter species (Piller et al., 2004). The increase in abundance (post-dam and reservoir construction) may be attributed to an expansion of preferential sandy habitats (Tipton et al., 2004) and may also be responsible for the recovered genetic patterns and larger effective population sizes. Ginson et al. (2015) also recovered a similar lack of within-river genetic structure for the Eastern Sand Darter, A. pellucida, a northerly distributed congener. As suggested for A. pellucida by Ginson et al. (2015), the lack of genetic structure for A. beanii (our study) is also likely caused by a substantial level of short and long dispersal among populations across continuous, preferred sandy habitats, which essentially homogenize genetic structure across large geographic areas resulting from high levels of movement.
Relative to A. beanii, populations of E. swaini within the Pearl River fell on the opposite end of the spectrum with regard to genetic structure and gene flow. The average genetic differentiation between population groups (BC, UP, and LP) were much higher than within group values. Populations within the Bogue Chitto River showed the greatest genetic structure and the lowest average values of genetic differentiation. While the geomorphology and habitat of the Pearl River has been altered through impoundment and channelization (Piller et al., 2004; Tipton et al., 2004), the Bogue Chitto River has been designated as a Louisiana Scenic River and has remained relatively undisturbed. The Natural and Scenic Rivers Regulations (Title 76) prohibits activities such as channelization, channel realignment, clearing and snagging, impoundments, and commercial clearcutting within 100 feet of the low-water mark (Louisiana Scenic Rivers Act, 1988; http://www.wlf.louisiana.gov/scenic-rivers). Resulting genetic differentiation estimates for E. swaini among three populations from other Gulf Coastal drainages (FST = 0.058–0.083; Fluker et al., 2010) were similar to our significant estimates between population groups (FST = 0.016–0.069). Although only three populations of E. swaini were examined in the Fluker et al. (2010) study and a comparison between our study and theirs requires a caveat for interpretation, the similarity in pairwise FST estimates on a similar spatial scale in a region that has also undergone impoundment and habitat alteration is noteworthy. Pairwise FST estimates were a magnitude larger for E. swaini compared to A. beanii and P. nigrofasciata in our study. Also, genetic data for E. swaini fit a model of isolation by distance, although distance only accounted for 44% of the genetic variation within the model. Furthermore, geographic distance could be a symptom of habitat fragmentation (Wang and Bradburd, 2014; Vanormelingen et al., 2015; Waters et al., 2015; Pilger et al., 2017). Geographic distances between populations were comparable between all species. Therefore, the recovered patterns are likely the result of the Gulf Darter's inherent preference for smaller streams and muddy, sand-gravel, or gravel substrates near vegetation (Ruple et al., 1984; Ross, 2001). The main channel of the Pearl River lacks heavily vegetated areas, thereby minimizing dispersal throughout the basin because of the lack of preferential habitats. Additionally, E. swaini displays reduced larval drift, with juveniles often found near spawning riffles, further supporting reduced movement of this species (Simon and Wallus, 2006). However, there are no studies examining larval drift of A. beanii and P. nigrofasciata, and direct comparisons of larval drift cannot be achieved. The Eastern Sand Darter (Ammocrypta pellucida), a close relative to A. beanii, has a primarily nonbenthic larval stage that may aid in downstream gene flow (Simon and Wallus, 2006; Ginson et al., 2015). Furthermore, previous studies have suggested that darter egg size is positively correlated with larval size and hatchling morphological development, and may predict larval dispersal ability (Paine, 1984, 1990; Turner, 2001). Based on previous life history studies (Mathur, 1973; Heins and Rooks, 1984; Ruple et al., 1984; Ross, 2001; Hughey et al., 2012), E. swaini has a larger egg size (1.1–2.0 mm) relative to A. beanii (0.92–1.0 mm) and P. nigrofasciata (0.81–1.60 mm), suggesting that larval drift may be more limited in E. swaini compared to the other two species. Gene flow, however, may not always exhibit a significant, positive correlation to larval transport, and the relationship can be obscured by ecological or habitat requirements of the species being compared (Turner, 2001). In our study, patterns of gene flow match suggested predictions of larval drift, but a more complete comparative dataset is warranted before acknowledging a definitive relationship between gene flow and larval drift in our species.
Percina nigrofasciata is more of a generalist in terms of its habitat occupancy, often inhabiting shallow riffles and deeper pools with silt, sand, or gravel substrates. The abundance of P. nigrofasciata seems to be positively correlated with the occurrence of submerged woody debris and aquatic plants (Mathur, 1973; Rakocinski, 1988; Ross, 2001). As previously mentioned, aquatic vegetation is relatively uncommon in the main channel, but is often present in tributaries. However, as a result of ongoing geomorphic and bank-destabilization issues in the basin, woody debris is abundant, but not continuously distributed, in both the main channel and tributaries throughout the Pearl River basin. The availability of some woody material in the main channel, however, suggests that habitat is not limited for this species, and movement throughout the basin may be less restrictive. Both pairwise FST estimates and STRUCTURE results indicate minimal levels of genetic differentiation between upper Pearl River populations and all other populations. Interestingly, a miniscule difference is seen between upper Pearl River, Bogue Chitto River, and lower Pearl River populations in the AMOVA results. IBD results suggest an isolation by distance pattern, but the model only explains 18% of the variation within the data. It has been hypothesized that body size can increase species dispersal and, thus, gene flow (McCallister et al., 1986). Of the three species in this study, P. nigrofasciata has a much larger body size (Ross, 2001). Therefore, despite patchiness in preferential habitat, Blackbanded Darters may be more capable of migrating between populations, which may be attributed, in part, to the presence of a reduced swim bladder, a structure not possessed by the other species in this study (Page and Swofford, 1984; Evans and Page, 2003). Most species of Percina have a reduced swim bladder that can vary in volume depending on the habitat and swimming behavior, thereby potentially increasing dispersal ability (Evans and Page, 2003).
Darters are the second most species-rich group in North America and display a wide variety of morphological, ecological, and life history characteristics (Page, 1983; Paine, 1990; Mayden, 1992; Ross, 2001; Boschung and Mayden, 2004). Despite the diversity of the group, there have been a limited number of comparative studies addressing the correlation between life history variation and gene flow on darters (Zimmerman, 1987; Turner et al., 1996; Heithaus and Laushman, 1997; Turner and Trexler, 1998; Faber and White, 2000). Specifically, previous studies have suggested that certain life history character differences (i.e., egg size and fecundity) in darters may correlate to differences in gene flow (Turner et al., 1996; Turner and Trexler, 1998; Faber and White, 2000). These types of life history-gene flow correlations have been recovered for other non-darter taxa as well (Whiteley et al., 2004; Neville et al., 2006; Palstra et al., 2007; Reid et al., 2008; Galarza et al., 2009). Neutral divergence among populations of stream species may be impacted by other life history and ecological characteristics as well and should vary depending on these characteristics (Whiteley et al., 2004). Likewise, genetic divergence may also be dependent on equilibrium between migration and genetic drift and the availability of suitable habitat between populations. Using a dendritic ecological network architecture, Pilger et al. (2017) found that patterns of genetic diversity or connectivity could be predicted by network structure for five of the nine species included in their study. These species were primarily small-bodied or habitat-specialists, similar to the darter species in our study. Their results were consistent with expectations of gene flow constrained by stream structure in fish species with smaller bodies and/or specialized habitats, and less constraint in larger-bodied and/or generalist species. We see similar patterns within our dataset. Although our species are all small-bodied habitat specialists, patterns of genetic differentiation, structure, and isolation by distance suggest that an expansion of sandy habitats and a more channelized structure in the Pearl River (Tipton et al., 2004) has increased suitable habitat and, thus, dispersal corridors for A. beanii, while limiting suitable habitat for E. swaini and, to an extent, P. nigrofasciata.
Bohonak (1999) conducted a meta-analysis of more than 300 animal species and noted that dispersal ability is strongly correlated with population genetic structure. Dispersal ability, in part, is limited by the availability of a species' preferred habitat. Understanding and quantifying habitat preferences for a particular species is relatively straightforward and known for many darter species. However, directly measuring dispersal ability is a difficult task and remains a mystery for many aquatic species. Furthermore, landscapes with homogenized habitat may provide challenges to species' movement, which in turn directly impacts gene flow and genetic structure. Our study highlights the importance of understanding basic life history parameters for inferring genetic structure, echoing conclusions from Turner and Trexler (1998). Clearly, the habitat preferences of the three study species in this study differ from one another, which in turn, have played a role in genetic structure. Possessing an understanding of basic life history information can help predict the temporal and spatial genetic structure of species within and among basins, and this is particularly relevant in regard to the current and future challenges that are impacting riverine systems and their inclusive aquatic faunas. Population connectivity is a key component to the persistence of fishes in riverine environments, and this is particularly true for species with specific habitat preferences, like darters. As habitats change as a result of anthropogenic influences and climate change, population structure of these species may also change, potentially isolating or connecting populations over time. As a whole, darters are an environmentally sensitive group, which serve as indicator species in many aquatic systems in eastern North America (Page, 1983; Mayden, 1992). Understanding the relationship between gene flow and life-history features in darters and other species will allow researchers to make predictions about current and future environmental change, and guide conservation management efforts in disturbed aquatic systems.
ACKNOWLEDGMENTS
This study was supported with funds from the Louisiana Department of Wildlife and Fisheries, State Wildlife Grants program to K. R. Piller. Additional support was provided to K. R. Piller from the Edward Schlieder Foundation through an endowment to Southeastern Louisiana University. We would also like to acknowledge the field assistance provided by Daniel Powell, Elizabeth Marchio, Mollie Cashner, and Malorie Hayes, and for the darter illustrations provided by Elizabeth Marchio.