Bithynia tentaculata is believed to have been extirpated from North America during the last glacial maximum. It was reintroduced into North America via the Great Lakes basin in the 1800's and has recently been expanding its geographic range. This snail serves as intermediate host for three trematodes that cause extensive recurring morbidity and mortality events in migratory water birds along the Mississippi River. Using twelve microsatellite loci for ∼200 individual snails from 11 populations in North America and Europe, we examined one of the three major geographic regions from which founding populations into the Great Lakes typically originate. Our data supports a single recolonization of North America into the Great Lakes Basin followed by subsequent introduction events from the Great Lakes to other large watersheds in North America. However, additional watersheds in Europe require sampling to confirm this result. No populations with genetic signatures indicative of North American glacial relics were found. The initial invasion of North America was likely not from the Ponto-Caspian basin, the usual source of freshwater invasive species to the Laurentian Great Lakes.
INTRODUCTION
The Laurentian Great Lakes of North America have been a hotspot for invasion by exotic species. Many ecologically damaging aquatic invasive species have been introduced into the United States (U.S.) via this route (Mills et al. 1993). Molecular data have been used to determine the source of invasion of various aquatic invaders. For example, using the mitochondrial cytochrome oxidase I gene, Gelembiuk et al. (2006) concluded that the source of invasion of zebra (Dreissena polymorpha) and quagga mussels (D. bugensis) into the Great Lakes was the Ponto-Caspian Sea basin (the Black, Caspian, and Azov Seas and their surrounding watersheds). This is congruent with other studies that have shown that the Ponto-Caspian Sea basin has been an important source of many aquatic invaders into the Great Lakes (Lee & Bell 1999, Ricciardi & MacIsaac 2000). Up to 70% of recent invaders in the Great Lakes (1985-2000) trace their source population to this region (Ricciardi & Maclsaac 2000, Brown & Stepien 2009, Keller et al. 2010). However, this is not the only possible source of aquatic invasive species. Shipping routes from European waters to the Great Lakes commonly originate in three locations: (1) the Black/Mediterranean Seas, (2) the North Sea, (3) or the Baltic Sea (Ricciardi & Maclsaac 2000, Grigorovich et al. 2003, Brown & Stepien 2009). These usual source regions for invaders into the Great Lakes provide an excellent starting point for comparing the genetic structure of invasive and potential source populations.
Bithynia tentaculata (L., 1758) had a Holarctic distribution prior to the last glacial maximum with shells found in Pleistocene fossil deposits from Lake Michigan, Illinois, U.S.A (Baker 1928). It is believed to have been extirpated from North America by glaciation events with subsequent recolonization through human-mediated introduction, Following the last glacial maximum, the first North American record of B. tentaculata was in Lake Michigan (presumably by passage through the Great Lakes Waterway, via the Hudson River) in 1871 (Baker 1928, Mills et al. 1993). It was speculated at that time that the snail was carried into Lake Michigan through ballast of timber ships arriving from Europe (Baker 1928). The species then spread throughout the Great Lakes region and into other U.S. waterways. It is now widespread in the Great Lakes, Northern Atlantic Coast drainages, isolated lakes in Montana, and most recently, in the Upper Mississippi River and Wolf River drainages, WI (Sauer et al. 2007)
Since the snail's introduction into the Mississippi River, first recorded in 2002 (National Wildlife Health Center, Madison, WI, unpublished data), parasites carried by B. tentaculata have caused recurring morbidity and mortality events in water bird populations during spring and fall migrations. The intestinal trematodes Cyathocotyle bushiensis, Sphaeridiotrema globulus, S. pseudoglobulus and Leyogonimus polyoon (Sauer et al. 2007, Mitchell & Cole 2008) cause intestinal hemorrhage and extensive mucosal damage. One snail can be infected with hundreds of infectious larval trematodes (Cole, unpublished data) and thus, by eating a small number of snails, a bird can receive a lethal infection in a short period of time (Sauer et al. 2007). From 2002-14 over 135,000 water birds consisting of 17 species have died in mortality events in Wisconsin, Minnesota and Illinois. These mortalities have been attributed to the four trematodes transmitted by B. tentaculata. The majority of these events have occurred in navigation pools (long stretches of river between dams) 7-11 of the Mississippi River and were predominately lesser scaup (Aythya affinis) and American coot (Fulica americana) (National Wildlife Health Center, Madison, WI, unpublished data).
Negative interactions between invasive species and native species are a leading cause of animal extinctions (Claver & Garcia-Berthou 2005) and freshwater gastropods are a highly threatened freshwater fauna with 74% of species categorized as imperiled or extinct (Johnson et al. 2013). In eutrophic lakes in upstate New York, B. tentaculata contributed to the decline of populations of the pleurocerids Elimia livescens, E. virginica, and Pleurocera acuta (Harman 1968k, Harman 1968, Harman & Forney 1970, Jokinen 1992). The proposal that the ability of B. tentaculata to both graze and filter-feed contributed to their competitive ability was supported by a finding that B. tentaculata adds biomass approximately 10 times faster than pleurocerids (Harman & Forney 1970) due to the higher efficiency of carbon and nitrogen assimilation associated with filter-feeding (Tashiro & Colman 1982). While these pleurocerids are relatively common – this is indicative that this introduced species is a potential competitor with other native pleurocerids. Furthermore, an initial study indicated native snails could suffer negative consequences from B. tentaculata invasion, largely due to increased exposure to trematode parasite larvae transmitted at high densities of B. tentaculata (Sandland et al. 2013). A further study in an experimental setting found several native snail species and B. tentaculata were equally infected with the larval stage (metacercariae) of an echinostome parasite suggesting a potentially positive effect of the invasive snail on natives may occur by diluting the parasite load of the entire snail community (Gladosky & Sandland 2014); however, this does not consider the ability of B. tentaculata to form very high population densities, which would serve to enhance overall parasite abundance and pose a threat to native snails.
Understanding the history of invasion and parent populations of B. tentaculata may lead to precautionary steps to be implemented to limit the spread of this species. Use of microsatellite data can be helpful in understanding the routes of introduction and pinpointing parent populations (Stepien et al. 2005, Brown & Stepien 2009). In this study, we used microsatellite data to determine the colonization route of the invasive populations of B. tentaculata into and throughout North America. We distinguish between four alternative hypotheses of potential colonization routes: (1) a single population of B. tentaculata was introduced into the Great Lakes from a single source population and has since dispersed; (2) Bithynia tentaculata were introduced multiple times into the Great Lakes from multiple sources and have since dispersed; (3) there were multiple introductions of B. tentaculata from Europe into geographically distinct locations within North America; (4) while some invasive populations may have been introduced from Europe, some populations may be glacial relics that persisted in North America.
METHODS
Bithynia tentaculata samples (Figure 1) were stored at —20°C in 70-100% ethanol after collection and were deposited in the Field Museum of Natural History (F numbers 344681-344697). This study was carried out in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. As an invertebrate, this species is exempt from the approval process of the Institutional Animal Care and Use Committees at UWL and University of Texas Rio Grande Valley. Total cellular DNA was extracted from a snip of foot tissue using the CTAB (Cetyltrimethylammonium bromide) method (Saghai-Maroof et al. 1984). Primers that target twelve loci were used to determine allele frequencies as described by Henningsen et al. (2010). DNA was amplified via PCR with a CAG tagged primer along with the associated primers, with and without a tag. PCR was performed in a 20uL reaction volume with the amplification mixture at concentrations as follows: Taq DNA Polymerase Thermopol Buff-2 0.05 U/µL, 0.15 mM dNTP, 1.5 mM MgCl2, 25 µg/mL BSA, and 1X Taq Thermopol buff-2 2 Tagged primers were included at a concentration of 0.025 µM, untagged primers at a concentration of 0.25 µM, and 2 µL of DNA template was added to each reaction. Cycling conditions consisted of 4 min at 94°C; followed by 32 cycles of: 94°C for 30 s, 55°C for 30 s, 65°C for one min followed by a final extension at 65°C for 3 min. Samples were then diluted with water 1:10 and genotyped at the University of Wisconsin Biotechnology center, on an ABI 3730xl Genetic Analyzer. Output files were analyzed using the auto rum setting in GeneMarker® (Holland et al. 2008) with a GS500 size standard and ABI template, to determine the size of alleles present at each locus.
Figure 1.
Collection sites of Bithynia tentaculata are indicated with filled black squares. Stars represent potential origination shipping routes for colonization from Europe. Locality information presented in Table 1.
A Bayesian analysis in STRUCTURE v 2.2 (Pritchard et al. 2000, version 2. from http://pritch.bsd.uchicago.edu/sofware/structure2_1.html.) was used to infer the number of populations (K) using Markov Chain Monte Carlo (MCMC), with five independent rims of 100,000 steps following a 100,000 step burn-in for each K from K=1 to K=12 (these represent the maximum and minimum possible values for K based on the number of populations sampled plus one). A test run of 200,000 steps and a 200,000 bum in was conducted to see if this level of iterations were required. The model assumed correlated allele frequencies among populations, sampling locations were informative about ancestry (LOCPRIOR), and followed an admixture model with a single value of lambda (λ=1.0) inferred for all populations. K was estimated based on the log likelihood score and posterior probability of K, Ln P(D) also known as L(K) (Falush et al. 2007) as well as the rate of change in the log likelihood score (ΔK) (Evanno et al. 2005). The log likelihood score was calculated following Evanno et al. (2005): ΔK = m|L(K + 1) — 2L(K) + L(K – 1)|/s[L(K)] for each K. Three values were used to estimate K, L(K), ΔK, and α. The best estimate of K is identified as the maximal value of L(K), yet as the true K is reached, L(K) at larger K values will plateau or even increase slightly (Evanno et al. 2005). The rate of change in the log likelihood score, ΔK, will be the highest at the true K. Finally, the lowest value for α indicates that most individuals are essentially from one population or another. The posterior probability of K, L(K), Ln(K), and α were output directly from the program, and ΔK was calculated using the equation above for K=1 to K=12. Once K was estimated, 5 runs of STUCTURE was used to calculate FST, HE, and HO for each of the K populations. Average expected and observed heterozygosity (HO, NA) for each of the three populations clusters determined by STUCTURE were calculated in Microsoft Excel. Finally, the North American populations were run separately to determine if the greater European diversity masked internal structure in North American populations.
Table 1.
Eleven populations examined for 12 microsatellite loci with statistics summarizing genetic variation within populations. Data presented are latitude, longitude, number of individuals (n), number of microsatellite alleles (NA), observed heterozygosity (HO) and expected heterozygosity (HE) averaged across loci. Dunaremete, Baracska, and Botanical Garden are all collection sites in the Danube River drainage. The order of populations matches the population number in Figure 2.
Genepop v 4.0.10 (Raymond & Rousset 1995), was used to estimate number of alleles (NA) for each population. GenAlEx 6.501 (Peakall & P.E. 2006, Peakall & Smouse 2012) was used to perform an analysis of molecular variance to calculate FST and F'ST of the three population clusters (regions) determined by STRUCTURE. We used these combinational, regional groups rather than individual populations to increase the sample size of each group. This conforms with the findings of Hale et al. (2012) that 25-30 individuals are needed per “group” for accuracy of microsatellite data. When grouping by the regions Danube (n=27), Lake Balaton (n=26), and North America (n=109) we have sufficient sampling for comparison among regions, although not comprehensive sampling for any region. GenAlEx v. 6.501 was also used to perform a genetic distance based Principal Coordinates Analysis (PCoA) on alleles from all populations.
As a final examination of patterns of population genetic structure we used Bottleneck 1.2.02 (Cornuet & Luikart 1996) to examine each of the 3 regions for signs of a recent genetic bottleneck using a Wilcoxon sign-rank test, to accommodate our limited loci and sampling under the two-phase-model of evolution as recommended for microsatellite loci (Luikart & Comuet 1997). In a genetic bottleneck, reduced population size results in loss of alleles and a decline in heterozygosity, recent bottlenecks should appear to have higher than expected genetic diversity compared to expectations from Hardy-Weinberg equilibrium.
RESULTS
Eleven populations of Bithynia tentaculata were sampled for 12 microsatellite loci from the native European range and from North America (Table 1). Structure runs with 200,000 burn-in and iterations were not different from 100,000 burn-in and 100,000 iterations all further analyses were carried out using the latter settings.
All three ad hoc estimates for K from our analysis using STRUCTURE, considered together, suggest a best estimate of K=3 (data not figured). All eleven populations grouped into one of these 3 clusters and provided the rationale for combining populations in further analyses. The three clusters include 1) the Danube population (all populations in the Danube River Basin, Hungary), 2) Lake Balaton population (includes Lake Balaton and River Zala—which flows into Lake Balaton, Hungary), and 3) North America. The FST value for the combined Danube populations is much lower than that for North America or Lake Balaton (Table 2). The North American populations analyzed without the European data resulted in K=1.
In Figure 2 which illustrates all eleven populations as well as the three combined populations, Lake Balaton stands out as being the least intermixed (also Table 2 with the highest among population FST value, 0.3379). The Danube has some contribution from Lake Balaton and from the population that is the source of the North American populations. The North American populations in Canada and Montana are the most heterogeneous, although lacking unique alleles. Examination for genetic signatures of a recent population bottleneck found no signature of this event in the Danube (Wilcoxon sign-rank test for two-phase-model, P =0.57), and Lake Balaton regions (P =0.688), however the North American region shows the signature of a recent bottleneck (P =0.012).
Table 2.
Results of AMOVA on populations grouped by regions. Genetic differentiation among populations, FST (+SD) and expected heterozygosity, HE (+SD) for K=3.
Figure 2.
STRUCTURE output of Q (or proportion of each individual attributed to each cluster). Clusters indicated by color, population by number: (Lake Balaton Cluster: 1-River Zala HUN, 2-Lake Balaton HUN; Danube River Cluster: 3-Dunaremete HUN, 4-Baracska HUN, 5-Botanic Garden HUN; North American Cluster: 6-Lake Winnibigoshish MN, 7-Georgetown Lake MT, 8-Rattlesnake Reservoir MT, 9-Upper Mississippi River, 10-Lake Winnebago WI, 11-Ottawa River, Canada).
The number of alleles at 12 microsatellite loci ranged from 1 - 39. In most populations, observed heterozygosity was lower than expected under Hardy-Weinberg equilibrium (Table 1). Allele-frequency divergence among these three populations is shown in Table 3. The two European populations are more similar to each other (0.1160) than to the North American population, which is roughly equally different from the two European populations (0.1484, 0.1547). The Danube River population had the largest number of private alleles (21), followed by North America (11) and the fewest in Lake Balaton (8).
Most loci not only had different frequencies, but usually unique alleles in each population. In general, allelic diversity is highest in the Danube populations, followed by North America, then Lake Balaton (Table 1). For example, a single well-sampled locus, Bt03, displays very different allele frequencies across all populations, each population also has unique alleles at this locus (not figured). This pattern is repeated in most other loci. However, some population structure was observed. Allele frequencies are more similar in populations such as the Montana Lakes and Ottawa River (i.e. those outside the areas adjacent to Lake Michigan) and have fewer alleles, none unique to those populations (not figured).
A principal coordinates analysis (PCoA) of the microsatellite allele data across all populations (excluding Baracska, a Hungarian population represented by a single individual) resulted in three significant PCoA axes, axis one explains 17% of the variation present, axis two 13.19 %, and axis three 9.28 %. A scatterplot showing all individuals grouped by population on axis one and two is shown in Figure 3. On both axis one and two (Figure 3) the Upper Mississippi River population and Lake Winnibigoshish populations have the widest variation in allelic diversity. The other populations are restricted to the lower left quadrant of the graph, encompassed within the diversity of those two populations. Axis three (not figured) distinguishes the European populations from Lake Winnebago, WI and the Montana populations.
Table 3.
Allele frequency divergence among regions (net nucleotide distance) calculated using STRUCTURE.
DISCUSSION
The goal of this study was to distinguish among four alternative hypotheses for introduction of B. tentaculata into North America. Hypothetical scenarios were proposed considering what is known of the possible invasion route and history of the species. The expected genetic consequences of each scenario are proposed based on the patterns observed in a review of genetic consequences of invasion in 80 species of animals, plants, and fungi by Dlugosch and Parker (2008) and glacial refugia by Maggs et al. (2009).
Hypothesis One incorporated a single introduction scenario: B. tentaculata was introduced once into the Great Lakes and dispersed into other North American watersheds. If this hypothesis were supported we would expect a low FST in North American populations compared to European populations, similar allele frequencies to the source, and very few private alleles in North America compared to European source. Hypothesis Two was developed around a multiple introduction scenario into the Great Lakes; after which B. tentaculata dispersed. If this hypothesis were correct, we would expect a higher FST, and greater heterogeneity in North America and many alleles from all across the founding European populations (albeit at lesser frequencies) and many private alleles if comparing North American population against a single source population. Hypothesis Three described a scenario with introductions from more than one European source population into the different North American watersheds in which B. tentaculata has now been found. If this hypothesis were correct, we would expect a high FST, very different allele frequencies both within North American populations and between North American and European populations, as well as private alleles unique to different North American populations. Hypothesis Four incorporated a European source of introduction for some of the North American populations, while other North American populations were assumed glacial relics. If this hypothesis were correct, we would expect a signal similar to Hypothesis Three but with private alleles in the glacial relic populations that are different from other North American and European populations.
Populations in Hungary compared to populations in North America
The B. tentaculata samples collected in Hungary are all part of the Ponto-Caspian Basin which contributes to the Black sea colonization route, and has contributed >70% of Great Lakes invaders (Ricciardi 2001, Gelembiuk et al. 2006). The concordance of the ad hoc estimates given in STRUCTURE (Figure 2) grouped all 11 populations into one of three clusters (K). However, Evanno et al. (2005) warn that these estimates must also coincide with scenarios that are biologically significant. Considering the a priori knowledge of sampling sites, K=3 corresponds with what would be predicted based on sampling localities from North America and two different Hungarian watersheds, the Danube River (Dunaremete, Baracska, and the Botanic Garden in Budapest sites) and Lake Balaton (Lake Balaton and River Zala sites) watersheds.
Figure 3.
Genetic comparisons between populations using a Principal coordinates analysis of the genetic distance matrix from the microsatellite allele data. Similar shapes represents populations from the same region or drainage basin. Axis 1 represents 17% and axis 2 represents 13.19% of the variation present.
Examination of allele frequencies and variability (Figure 2, Figure 3, Table 3) from the Danube River and Lake Balaton supports the idea that the populations could be genetically distinct with the North American samples grouped as a third distinct population. There appears to be some shared alleles with alleles found in the Hungarian populations also present in the North American populations (Figure 2). This may indicate some contribution to ancestry of or a shared common ancestry with the North American populations, however additional European populations must be sampled for a robust comparison to be made. Some preliminary estimates can be made from the allele-frequency divergences of each population (Table 3). Even though the two Hungarian populations are distinct (Figure 3), they are more similar to each other than either is to the North American population (Table 2). This may indicate that the Hungarian populations are not the source of the North American populations. This inference is also supported by the number of private alleles, as there are more private alleles in samples from Hungary than in samples from North America, supporting the idea that Hungary is not the source of the North American populations. However, this is not definitive; the source could be from further downstream on the Danube, which would explain some of the same alleles present seen in the STRUCTURE figure (Figure 2) and as overlap in the PCoA (Figure 3).
These data do not support the Black sea colonization route as the source of the North American invasion, however this route cannot be definitely ruled out considering its large range and our limited sampling. Samples from further downstream on the Danube River, or within the Black Sea, will be necessary before this entire watershed can be definitively excluded as a source for the North American invasion. Sampling of the other two likely colonization routes through the North Sea and the Baltic Sea is also needed to determine if those routes may be the source population of the North American invasion.
Populations within North America
The North American populations were all combined into a single population by the ad hoc estimates in STRUCTURE which points toward a recent shared ancestry of all the populations within North America. Further when analyzed separately to determine if the European diversity masked variation within North American, STRUCTURE found the most support for a single North American population. This evidence supports Hypotheses one or two that B. tentaculata was introduced into the Great Lakes and dispersed from there. This is supported by few private alleles in North America (compared to the European populations) and that the North American FST is relatively low compared to the native European populations. Another source of data offers some support for this scenario. In a relatively recent M.S. thesis Whalen (2011) used 11 microsatellite loci of which 4 overlap with this study and several of the same populations but also including 2 population not included in this study which are from Eastern Wisconsin near the Great Lakes. They found the populations near the Great Lakes (e.g. Lake Winnebago, WI) were probably “parent” populations to the Lake Onalaska and Lake Winnibigoshish populations included in both studies.
The STRUCTURE analysis does not appear to support Hypothesis Three, that there are multiple European source populations for the North American populations, unless Europe has very homogeneous populations, which is unlikely with the amount of divergence we observed in just the two Hungarian watersheds. Hypothesis Four also appears unlikely given our data set, though it is hard to distinguish private alleles from potential glacial relics without more extensive sampling of specimens from Europe for comparison. It is difficult at this point to confirm whether Hypothesis One or Two is more likely, as there were not enough data from European specimens to compare source populations. However, by comparing allele frequencies among populations, we can get a hint of which hypothesis (One or Two) is more likely.
With the exception of two loci, Bt22 and Bt40, the allele frequencies across the three populations differ, once again supporting the hypothesis of a different source population for invasion than those sampled in Hungary (data not figured). There were alleles present in snails from North America that were absent in snails collected in Hungary. While these may be attributed to new alleles arising in the population, such a scenario is unlikely across so many loci. There were also alleles present in the Hungarian populations that are not in snails from North America. Founder effects could account for this; however, if this were the case, the other alleles present would most likely be at similar frequencies, which they were not. It is possible that some snails from Hungary were mixed with other European source populations that then gave rise to the North American populations. If this were the case, there would likely be more heterogeneity in the North American populations (Figure 2). Most of the North American loci are dominated by one allele indicating a significant recent genetic bottleneck occurred, this is also supported by the Bottleneck analysis finding the signature of a recent bottleneck only in the North American population. This signature could be due to a relatively small initial invasive population from a single source population and subsequent bottlenecks with colonization of additional watersheds. This coincides with other invasions into the Great Lakes that have been shown to be from a single source population (Ricciardi & MacIsaac 2000, Brown & Stepien 2009). At this point it is still speculative, but Hypothesis One is more likely, and the higher FST (Table 3) may be an artifact of founder effects which can elevate FST levels (Weir & Cockerham 1984).
Our data suggest B. tentaculata has dispersed across the U.S. from a single initial colonization, not from multiple invasions from different sources. It also appears that all sampled North American populations are recent recolonizations, not glacial relic populations. Given the few populations sampled in the European range of B. tentaculata, the European source of the introduction into North America is still unknown, and will require further study of the European range. However, it does appear, based on the data available, that the most common route for invasion into the Great Lakes, from the Ponto-Caspian Sea Basin through the Black sea, is not the likely source of introduced Bithynia tentaculata in North America.
ACKNOWLEDGEMENTS
The authors would like to thank individuals that contributed their insight, time, and specimens to the project: Russell Minton, Anindo Choudbury, Jenny Sauer, Roger Haro, Mark Sandheinrich, Ben Walker, Greg Sandland, Rachelle Amundson, Adam Prise, Kirk Gallant, and Josh Laurila. This project was funded by the UWL River Studies Program, United States Geological Survey, Upper Mississippi Environmental Research Center, National Great Rivers Research and Education Center, University of Wisconsin Institute for Race and Ethnicity, the UWL faculty research grant program, and the UWL College of Science and Health Dean's Fellowship program. Several anonymous reviewers contributed greatly to the improvement of the manuscript, we appreciate their critical insights. Use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the US Government.
