Myriad anthropogenic factors have led to substantial declines in North America's freshwater mussel populations over the last century. A greater understanding of mussel dispersal abilities, genetic structure, and effective population sizes is imperative to improve conservation strategies. This study used microsatellites to investigate genetic structure among mussel beds and estimate effective population sizes of a common North American mussel species, Amblema plicata, in the Little River, Oklahoma. We used five microsatellite loci to genotype 270 individuals from nine mussel beds distributed throughout the river and one of its tributaries, the Glover River. Our results indicate that subpopulations of A. plicata in the Little River are genetically similar. Upstream subpopulations had less genetic diversity than sites located downstream of the confluence of the Glover and Little rivers. Downstream subpopulations were primarily assigned to the same genetic group as upstream subpopulations, but they were admixed with a second genetic group. Low flows during droughts likely influenced the observed genetic structuring in A. plicata populations in the Little River. Additionally, downstream subpopulations may be admixed with a genetically distinct population of A. plicata, which may account for the increased genetic diversity. Estimates of effective population sizes (Ne) of large mussel beds were low compared to the total abundance (N) of A. plicata. While our data have limitations, they provide important information on the spatial scale at which conservation plans should focus and the population sizes that should be sustained through relocation and restocking programs.
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.
Study Area, Species, and Sampling
Amblema plicata is a common, wide-ranging mussel species found throughout central and eastern North America, and it is one of the most abundant species in the Little River (Vaughn and Taylor 1999). Amblema plicata have been estimated, based on growth rate, to reach sexual maturity at six years of age (Haag 2012). Male A. plicata, like all unionids, broadcast their sperm into the water column to fertilize females downstream (Haag 2012). After fertilization and larval development, female A. plicata release larval threads, which are mucus threads with attached glochidia, into the water column to infect host fish (Haag 2012). Amblema plicata is a host generalist but typically uses fish in the sunfish (Centrarchidae) and perch (Percidae) families (Freshwater Mussel Host Database 2017).
During the summers of 2015 and 2016, we collected A. plicata tissue samples from eight mussel beds in the Little River and one mussel bed in the Glover River, a tributary of the Little River, in southeast Oklahoma (Fig. 1). Three high-gradient sites were above the confluence of the Little and Glover rivers, while six low-gradient sites on the Little River were below the confluence of the rivers. The Little River is influenced by two large impoundments and two small low-head dams. Its main stem is impounded by Pine Creek Dam (constructed in 1969), while the Mountain Fork River, a major tributary of the Little River, is impounded by Broken Bow Dam (constructed in 1968), a hypolimnetic release dam (Vaughn and Taylor 1999; Matthews et al. 2005). Cold water from Broken Bow Dam has eliminated most mussels in the lower Mountain Fork River and the lower Little River below the confluence of these two rivers (Vaughn and Taylor 1999). The Glover River is unimpounded and enters the Little River approximately 30 km downstream of Pine Creek Dam. The two low-head dams are located on the main stem of the Little River, one between the outflow of Pine Creek Reservoir and the confluence with the Glover River, and the other between the Glover and Mountain Fork rivers' confluences with the Little River (Fig. 1).
Tissue samples from 30 individual A. plicata were collected from each site for a total of 270 samples. Tissue samples were collected from the first 30 individuals found within quadrats at large sites; if fewer than 30 individuals were collected within quadrats, then the remaining individuals were located through semiquantitative timed searches. At small mussel beds, tissue samples were collected randomly from 30 A. plicata individuals that were located during hour-long semiquantitative time searches. Five of the nine sites (sites: 1, 5, 6, 7, and 8) were large mussel beds (>50 m long). These sites were quantitatively sampled with quadrats (Vaughn et al. 1997). Twenty 0.25 m2 quadrats were placed randomly along transects throughout the mussel bed and excavated to a depth of 15 cm. The density (mussels/m2) of A. plicata was calculated using the quadrat data, and the total abundance of A. plicata in each large mussel bed was estimated by multiplying the density by the area of each bed. Semiquantitative timed searches (Vaughn et al. 1997) were conducted for an hour at four small mussel beds (<50 m long; sites: 2, 3, 4, and 9). Mussels were located tactilely or visually while snorkeling or scuba diving over the mussel bed. We collected approximately 20 mg of mantle tissue from each mussel and stored it in 95% ethanol. We also measured the shell length of every mussel sampled.
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).
We used GenAlEx 6.502 to calculate expected (He) and observed (Ho) 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 FST 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 (FST 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 FST 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.
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. Ho = observed heterozygosity. He = 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 (Ne) 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 (Ne/N) in each large mussel bed by dividing the effective population size by the total abundance.
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.
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 FST values ranged from –0.0035 to 0.0596, with significant (FST ≠ 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 (FST) 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 (FST ≠ 0) genotypic differentiation between upstream (1–3) and downstream (4–7) subpopulations (FST = 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.
Pairwise geographic distances in river kilometers and FST values above and below the diagonal, respectively. Bold font indicates significant genetic differentiation between subpopulations (FST ≠ 0, df = 10, P < 0.05).
Large mussel beds had A. plicata densities ranging from 3.5 to 9.4 individuals/m2 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 (Ne/N) among three large mussel beds was 0.071, but values were highly variable.
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.
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. Ne = effective population size. Ne/N = proportion of individuals breeding. Negative Ne values can be explained by sampling error and interpreted as an infinite Ne (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 (Ne/N). Mean estimates of Ne/N ranged from 0.10 to 0.11 from 192 published estimates across 102 species (Frankham 1995). We found that estimates of Ne 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 Ne was not infinite, were highly variable, ranging from 0.002 to 0.209 (mean Ne/N = 0.071). Other broadcast-spawning species have widely variable Ne/N ratios. A freshwater mussel species (Popenaias popeii) endemic to the Rio Grande drainage in the United States and Mexico had an Ne/N ratio of 0.12 in the Black River in New Mexico (Inoue et al. 2015). The estimated Ne/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 Ne/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 Ne. For Lampsilis cariosa from three river drainages in Maine, Ne ranged from 150 to 1,850 individuals across nine sites (Kelly and Rhymer 2005), while Lampsilis cardium exhibited Ne 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 Ne 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.
We thank Katie Murphy, Kyle Broadfoot, and Traci DuBose for assisting with field work and Ann Harris, James Cureton, and Billy Culver for help with laboratory work and data analysis. Lawrence Weider and Richard Broughton provided helpful feedback on the manuscript. This research was funded by the Oklahoma Department of Wildlife Conservation (SWG grant no. F1401225), the U.S. Fish and Wildlife Service (Cooperative Agreement no. F14AC01058, Project no. A15-0044), and the University of Oklahoma.
Appendix A: PCR Reaction Mixes and Conditions
PCR reaction mixes and conditions for all nine loci. aJames Cureton, University of Oklahoma, personal communication; bGalbraith et al. (2015).