INTRODUCTION

Freshwater mussels (Bivalvia: Unionoida, hereafter “mussels”) are a highly diverse and imperiled group of animals. With approximately 300 species, North America has the highest diversity of mussels. Roughly 70% of endemic species are considered threatened, vulnerable, endangered, or extinct (Lydeard et al. 2004), and even common species are decreasing in abundance (Anthony and Downing 2001; Vaughn et al. 2015). Declines in mussel populations can be attributed to a variety of factors, such as habitat destruction and fragmentation, introduction of invasive species, pollution, and commercial exploitation of shells for the pearl and pearl button industries (Lydeard et al. 2004; Strayer et al. 2004). Conserving the remaining mussel fauna is a priority, but without understanding the basic population biology of mussels, developing successful conservation plans may be impossible. One emerging conservation tool is the propagation and restocking of mussels (FMCS 2016). To use this tool successfully, we need to understand the spatial scale of genetic structuring and effective population sizes of mussel beds.

Mussels often occur in dense aggregations or beds that are separated by stretches of river with no or very few mussels (Strayer et al. 2004). Mussel larvae (glochidia) are obligate ectoparasites on fish, but adults are sedentary (Barnhart et al. 2008). Thus, gene flow among and within beds must occur either via movement of glochidia attached to host fish, or through the downstream drift of sperm, juvenile mussels, or glochidia before they are attached to fish (Schwalb et al. 2011; Ferguson et al. 2013; Irmscher and Vaughn 2015, 2018). Mussel beds that are connected through gene flow can be considered one large metapopulation with individual beds serving as subpopulations (Vaughn 2012). Thus, the overall genetic structure in mussels should be the sum of connectivity among subpopulations due to gene flow and the isolation of subpopulations due to dispersal barriers such as impoundments or stretches of unsuitable habitat (Galbraith et al. 2015).

Genetic structure can indicate the spatial scale at which gene flow is occurring. Habitat disturbances and fragmentation can influence the genetic structure of mussel beds by blocking gene flow among beds (Watters 1996; Strayer et al. 2004; Newton et al. 2008; Schwalb et al. 2011; Galbraith et al. 2015). Multiple studies have evaluated genetic structure in unionid mussels (Berg et al. 1998; Kelly and Rhymer 2005; Elderkin et al. 2006; Berg et al. 2007; Zanatta and Wilson 2011; Galbraith et al. 2015; Jones et al. 2015; Inoue and Berg 2017; Hoffman et al. 2018), but few studies have estimated effective populations sizes of mussel beds and compared these estimates to total bed populations (Inoue et al. 2015).

Our goal was to better understand connectivity between mussel beds and the effective population sizes of these beds. We used microsatellites to evaluate population genetic structure and effective population size in a common, widespread mussel species, the threeridge mussel (Amblema plicata), in a medium-sized south-central U.S. river known for its diverse and relatively healthy mussel and fish populations (Vaughn and Taylor 1999; Vaughn 2003; Matthews et al. 2005). The Little River is fragmented by both large and low-head dams (Vaughn and Taylor 1999; Allen et al. 2013), which might restrict gene flow and result in distinct genetic clusters of individuals upstream and downstream of these dams. Our objectives were to assess the population genetic structure of A. plicata in the Little River and estimate the effective population size of each sampled mussel bed.

METHODS

DNA Extraction and Genotyping

DNA was extracted using the methods of the Qiagen DNeasy blood and tissue kit (Qiagen, Hilden, Germany). We successfully amplified nine microsatellite loci using primers developed for Amblema neislerii: Anec101, Anec114, Anec122, Anec126, Aned103, Aned104, Aned108, Aned126, and Aned140 (Díaz-Ferguson et al. 2011) and a variety of different polymerase-chain-reaction (PCR) conditions described in Table A1. Additionally, we tested two other loci that did not successfully amplify with PCR (Anec103 and Aned132). We used the ILS600 red size standard (Promega, Madison, Wisconsin, USA) and analyzed the PCR products on an ABI 3130xl genetic analyzer (Applied Biosystems, Foster City, California, USA). Alleles were binned and scored in GeneMapper V3.7 (Applied Biosystems).

Figure 1.

Sampling sites in the Little River drainage, Oklahoma, USA.

Data Analysis

We used GenAlEx 6.502 to calculate expected (H_e) and observed (H_o) heterozygosities and number of alleles per locus and site (genetic diversity) and to check for deviations of genotype frequencies from Hardy-Weinberg expectation (HWE; Peakall and Smouse 2006, 2012). We checked for linkage disequilibrium within and among mussel beds with GENEPOP V4.6 (Raymond and Rousset 1995; Rousset 2008). We estimated null-allele frequencies with MICRO-CHECKER (van Oosterhout et al. 2004). Subpopulation pairwise F_ST was calculated with GENEPOP V4.6 (Raymond and Rousset 1995; Rousset 2008). We ran exact G-tests to check for significant allelic differentiation and genotypic differentiation (F_ST values) in GENEPOP V4.6. Due to the genetic similarity and geographic proximity (Table 3) of sites 6, 7, and 8, we combined them into a single subpopulation with 90 individuals before we tested for isolation-by-distance (IBD) and evaluated genetic structure. We ran paired Mantel tests with 9,999 permutations in GenAlEx 6.502 (Peakall and Smouse 2006, 2012) comparing pairwise F_ST and river distance between sites (measured with the path function in Google Earth Pro) to analyze IBD across all sites. Pairwise geographic distances from the combined subpopulation of sites 6–8 were taken from the midpoint of the distance between sites 6 and 8.

We used STRUCTURE (Version 2.3.4; Pritchard et al. 2000), which uses a Bayesian clustering method to assign individuals to populations and infer genetic structure, to evaluate population genetic structure. Across all runs, we assumed independent allele frequencies and allowed for individuals to be admixed among subpopulations. We used the sampling location of each individual as prior information to assist clustering (LOCPRIOR model). Each run had an initial burn-in period of 50,000 and was followed by an additional 100,000 Markov Chain Monte Carlo (MCMC) replicates. We ran 10 iterations for each value of K (genetic clusters). Values of K ranged from 1 to 7 and were based on the number of subpopulations. STRUCTURE HARVESTER (V0.6.94; Earl and vonHoldt 2012) was used to determine the number of genetic clusters (K) that best fit the data. The value of K that corresponds to the greatest P(X|K) value was identified as the number of genetic clusters in the study area, which according to Evanno et al. (2005) is a good predictor of the real number of genetic clusters. We used the FullSearch algorithm in CLUMPP (Version 1.1.2; Jakobsson and Rosenberg 2007) to find the optimal alignment of 10 replicate cluster analyses from STRUCTURE with K = 2, and we used DISTRUCT (Version 1.1; Rosenberg 2004) to graphically represent the individual assignment scores of all 270 individuals across the seven subpopulations. Upstream beds (sites 1–3) and downstream beds (4–9) were grouped together to form two populations and checked for evidence of recent population declines using a Wilcoxon test in BOTTLENECK (Version 1.2.02; Cornuet and Luikart 1996) under the infinite-allele and two-phase models with 1,000 iterations.

Table 1.

Genetic diversity metrics per locus and site for A. plicata in the Little River. Figure 1 shows the locations of each site. n = number of individuals genotyped per locus. H_o = observed heterozygosity. H_e = expected heterozygosity. Bold font indicates departures from Hardy-Weinberg expectation. Values in the rightmost column are for means across all nine sites.

We used NeEstimator (Version 2.01; Waples and Do 2008; Do et al. 2014) to calculate the effective population size (N_e) of each mussel bed using the linkage disequilibrium (LD) method with a critical value of 0.05. We calculated the proportion of reproductively active A. plicata individuals (N_e/N) in each large mussel bed by dividing the effective population size by the total abundance.

Table 2.

Null allele frequencies per locus across the nine sites. Negative null allele frequencies indicate a heterozygote excess at a given locus and site. Bold font indicates the presence of null alleles at a given locus and site due to a significant excess of homozygotes (P < 0.05), which is calculated using Fisher's combined probability test.

RESULTS

We genotyped 30 A. plicata individuals from each bed for a total of 270 individuals. The number of alleles ranged from four to 31 across loci and beds (Table 1). Genetic diversity was higher across the four downstream subpopulations (mean number of alleles per locus [±SE] = 18.33 ± 2.78) than at the three upstream subpopulations (14.11 ± 2.53). Mean observed heterozygosities ranged from 0.50 to 0.63 and mean expected heterozygosities ranged from 0.66 to 0.79 among sites (Table 1). Because deviations from HWE due to heterozygote deficiencies occurred at six or more sites for four loci (Anec101, Aned103, Aned104, and Aned108), and these loci also had null alleles at high frequencies (Table 2), they were not included in subsequent analyses. The remaining five loci deviated from HWE at three or fewer sites with null alleles present at low frequencies at two or fewer sites. We found no evidence of large allele dropout or scoring errors due to stuttering. There was no evidence of linkage disequilibrium between loci across all subpopulations. Linkage disequilibrium between two or fewer loci was detected within six sites. Microsatellite genotypes of all 270 individuals can be obtained by contacting the authors.

Pairwise F_ST values ranged from –0.0035 to 0.0596, with significant (F_ST ≠ 0, df = 10, P < 0.05) allelic and genotypic differentiation at 13 of the 21 subpopulation pairs, while pairwise geographic distances between sites ranged from 12.40 km to 155.80 km (Table 3). The paired Mantel test did not find a significant relationship between genetic (F_ST) and geographic (river km) distance within the Little River (R = 0.51, P = 0.09), suggesting a lack of isolation-by-distance. Although there was significant (F_ST ≠ 0) genotypic differentiation between upstream (1–3) and downstream (4–7) subpopulations (F_ST = 0.0102, df = 10, P < 0.001), analysis of population genetic structure revealed a single genetic cluster (K = 1) among the seven A. plicata subpopulations. Downstream subpopulations exhibited almost no genetic structure among sites, and they are genetically similar to upstream subpopulations (Fig. 2). Individual-based population assignment scores indicated that downstream subpopulations had a higher degree of admixture between two genetic groups (blue and orange); however, both the upstream and downstream subpopulations were overwhelmingly assigned to the same genetic group (blue). There was no evidence of recent population bottlenecks in upstream and downstream populations.

Table 3.

Pairwise geographic distances in river kilometers and F_ST values above and below the diagonal, respectively. Bold font indicates significant genetic differentiation between subpopulations (F_ST ≠ 0, df = 10, P < 0.05).

Large mussel beds had A. plicata densities ranging from 3.5 to 9.4 individuals/m² and estimated total abundance ranging from 1,572 to 61,776. Small beds had catch per unit effort ranging from 34 to 82 individuals/hr (Table 4). The effective population sizes of the five large beds ranged from 81.4 (95% CI: 28.7–Infinite) at site 7 to Infinite at sites 5 and 6; and the effective population sizes of small beds ranged from 100.8 (95% CI: 21.9–Infinite) to Infinite (Table 4). The mean proportion of individuals breeding (N_e/N) among three large mussel beds was 0.071, but values were highly variable.

Figure 2.

Bayesian clustering analysis of A. plicata genetic structure among seven subpopulations in the Little River, with K = 2.

DISCUSSION

Our results indicate that upstream and downstream subpopulations of A. plicata in the Little River are genetically similar. The three subpopulations upstream of the confluence of the Glover and Little rivers were overwhelmingly assigned to one genetic group, while the four downstream subpopulations were admixed between two genetic groups with 70–80% of each individual-based population assignment score being assigned to the same genetic group as upstream subpopulations. While most studies have found little or no within-river genetic structuring of mussel populations where stream flows are consistent and unfragmented (Szumowski et al. 2012; Ferguson et al. 2013; Galbraith et al. 2015; Jones et al. 2015; Inoue and Berg 2017), our study and one other have found genetic structuring at microsatellite loci among mussel populations within a stream. Inoue et al. (2015) found genetic differences in upstream and downstream populations of Popenaias popeii in the Black River of New Mexico. Although we found no evidence of a recent population bottleneck at the upstream sites, the low mean number of alleles across loci at upstream sites suggests that these sites have lower genetic diversity than downstream sites. Two possible mechanisms underlying these differences in genetic diversity are (1) restricted gene flow between upstream and downstream subpopulations during periods of drought and (2) loss of rare alleles by genetic drift associated with decreases in upstream population sizes during droughts. Droughts are common and cyclical in the south-central USA and have been shown to lead to decreases in mussel population sizes in rivers in this region (Galbraith et al. 2010; Atkinson et al. 2014; Vaughn et al. 2015). Upstream reaches of the Little River are smaller and higher gradient than downstream reaches, and during droughts riffle areas can become dry or very shallow (Vaughn 2003; Matthews et al. 2005). Gene flow among mussel beds requires sufficient flow for the movement of fish hosts, juveniles, sperm, and/or unattached glochidia. Irmscher and Vaughn (2015) found that the movement of fish hosts in the Little River was restricted during droughts. Thus, low-water conditions during droughts may restrict gene flow between upstream and downstream populations or decrease upstream population sizes and exacerbate genetic drift. We did not observe genetic structuring among sites in the lower river (below the confluence with the Glover River), and this is likely because there is sufficient gene flow among these sites. Downstream sites were more genetically diverse than upstream sites, which may be due to admixture from another genetically distinct population of A. plicata from further downstream.

Table 4.

Demographic metrics for A. plicata at each site. Area, density, total number of A. plicata individuals, and proportion of individuals breeding were estimated only for large beds where quantitative sampling using quadrats was completed. Small beds are indicated with an asterisk by the site number; relative abundance (as CPUE = catch per unit effort) was measured for these. N = total number of individuals. N_e = effective population size. N_e/N = proportion of individuals breeding. Negative N_e values can be explained by sampling error and interpreted as an infinite N_e (Do et al. 2014).

Although the Little River is fragmented by both small low-head dams and large dams, as well as the reservoirs formed by them, we did not see evidence of interrupted gene flow, but this could be due to the long generation times of mussels. Pine Creek Dam (constructed in 1969) impounds the river itself and thus impedes host-fish dispersal between beds upstream and downstream of the dam. Broken Bow Dam (constructed in 1968) is a hypolimnetic release dam on a major tributary of the Little River, the Mountain Fork River. Cold water constantly flowing into the Little River via the Mountain Fork has caused significant declines in mussel abundance downstream from the release (Vaughn and Taylor 1999), along with preventing host-fish dispersal. Finally, small low-head dams on the Little River main stem also may restrict fish movement. However, we did not find distinct genetic clusters upstream and downstream of any of the dams; rather, we found that sites upstream and downstream of Pine Creek Dam (sites 1–3) assigned to the same genetic group. These populations are now isolated, but we likely did not see the genetic signal yet because the time of isolation is relatively short given the long generation times of mussels. Many mussel species, including A. plicata, are long-lived organisms with long generation times (Haag and Rypel 2010). In a study of growth and longevity of mussels in southeast Oklahoma using dendrochronological techniques, maximum ages of adult A. plicata from three rivers ranged from 53 to 79 years old (Sansom et al. 2016). Thus, there have likely not been enough generations for differentiation to occur upstream and downstream of large dams through the loss of alleles due to genetic drift. Low-head dams may not completely block gene flow because the Little River experiences frequent high flows (Matthews et al. 2005) and fish hosts may be able to freely migrate over them during floods. Other studies also have failed to show the isolating effects of dams on genetic structure in mussels, again, likely a consequence of the long generation times of mussels (Kelly and Rhymer 2005; Szumowksi et al. 2012).

Few studies have compared the effective population sizes of mussel beds to the estimated total population size (N_e/N). Mean estimates of N_e/N ranged from 0.10 to 0.11 from 192 published estimates across 102 species (Frankham 1995). We found that estimates of N_e for A. plicata were small relative to N estimated by quantitatively sampling mussels. Proportions of breeding mussels in the three large beds where N was estimated and the estimated N_e was not infinite, were highly variable, ranging from 0.002 to 0.209 (mean N_e/N = 0.071). Other broadcast-spawning species have widely variable N_e/N ratios. A freshwater mussel species (Popenaias popeii) endemic to the Rio Grande drainage in the United States and Mexico had an N_e/N ratio of 0.12 in the Black River in New Mexico (Inoue et al. 2015). The estimated N_e/N ratio for Pacific oysters (Crassostrea gigas) was less than 10^–6 (Hedgecock et al. 1992). Another broadcast-spawning species, sea bass (Atractoscion nobilis), had N_e/N ratios ranging from 0.27 to 0.40 (Bartley et al. 1992). Our results are corroborated by other studies of mussels that have found relatively low values of N_e. For Lampsilis cariosa from three river drainages in Maine, N_e ranged from 150 to 1,850 individuals across nine sites (Kelly and Rhymer 2005), while Lampsilis cardium exhibited N_e from 1.5 to 205.8 individuals across eight sites in Ohio (Ferguson et al. 2013). The effective population size of Quadrula fragosa in the St. Croix River was estimated to be roughly 150 individuals (Roe 2010). Estimates of effective population sizes are more informative when compared to total population sizes. Restocking programs should be designed to ensure that N is sufficiently large to lead to values of N_e that are high enough to offset the effects of genetic drift on target mussel populations.

This study contributes to our understanding of the population genetics of a common mussel species, but there are limitations to the results. While only nine of 11 loci were successfully amplified with PCR, an additional four loci consistently deviated from HWE due to heterozygote deficiencies, likely due to high frequencies of null alleles. Additional loci would provide more resolution when evaluating genetic structure and estimating effective population sizes. Furthermore, to date, all microsatellite studies of A. plicata genetic structure have used primers developed for A. neislerii (Díaz-Ferguson et al. 2011). Primers developed specifically for A. plicata may amplify more successfully.

This study provides important information on the genetic structure and effective population size of a common mussel species, which can be used to help manage and conserve not only common species but rare ones as well. Galbraith et al. (2015) found that mussel genetic structure did not vary as a function of rarity. Because sampling for common species is less time-intensive and thus less expensive than sampling for rare ones, common species could be used as surrogates for rare species when attempting to understand the population genetic structure of mussels.

We found that A. plicata subpopulations within a large extent (156 km) of the Little River were genetically similar, but genetic structuring was present within this reach and is likely influenced by flow conditions and possibly admixture of downstream subpopulations with a genetically distinct subpopulation. Although our data have limitations, our results provide useful information on the spatial scale at which conservation plans should focus and the population sizes that should be sustained through relocation and restocking programs. In stretches of river with genetically similar beds, individuals could be translocated from healthy beds to beds that are declining (Galbraith et al. 2015). Additionally, managers could use individuals from stable beds to propagate mussels for stocking into beds that are suffering from declines. Restocking programs should be designed to ensure that total population sizes are high enough to lead to effective population sizes high enough to offset the effects of genetic drift. For any mussel conservation program to be successful, there must be high-quality, unfragmented habitat into which to translocate or restock mussels. An understanding of the mussel population genetics of a region is important before implementing conservation strategies, including propagation and relocation.