Genetic models reveal historical patterns of sea lamprey population fluctuations within Lake Champlain

Sample collection and mitochondrial DNA sequencing

To estimate historical population fluctuations within Lake Champlain, fin tissue samples from spawning-phase sea lamprey were obtained from the US Fish and Wildlife Service (USFWS) in May–June 2009. All tissue collection was conducted as part of routine USFWS sampling; the authors of this study received tissue, but did not handle any vertebrate specimens. First, to estimate fluctuations based on mtDNA sequence data, samples were collected at three Lake Champlain tributaries: Great Chazy River (n = 33), Malletts Creek (n = 33), and Beaver Brook (n = 28), representing the northern, central, and southern regions of Lake Champlain, respectively. A broad geographic sampling regime was used in order to test for population structure. All tissue samples were stored in 95% non-denatured ethanol and genomic DNA was extracted using DNeasy Blood and Tissue Kits (Qiagen, Venlo, Netherlands).

The two non-coding regions of the mtDNA genome, NCI and NCII, were amplified in all Lake Champlain samples using PCR with the lamprey-specific primers CR1 (Waldman et al., 2004) and LampR (5′-AATAGACGGTTGGTGGGACA-3′). PCR reactions were performed in 25 µl volumes with the following reagents: 0.2 µl Qiagen taq polymerase (5 units/µl), 2.5 µl 10X PCR buffer (with 1.5 mM MgCl₂), 10 µM each primer, 10 µM dNTP, and 50–100 ng template DNA. Thermal cycler settings were set at an initial denaturation at 95 °C for 5 min followed by 40 cycles of 95 °C for 45 s, 56.5 °C for 45 s, and 72 °C for 1 min; and a final extension at 72 °C for 10 min. PCR products were cleaned using the Wizard SV Gel and PCR Clean-Up Kit (Promega, Madison, Wisconsin, USA) and eluted in 25 µl H₂O.

To sequence the Lake Champlain samples at both non-coding regions, forward extension reactions were run for each individual using the forward primers CR1 (for NCI) and CASSFN (for NCII) (5′- GACCCCTAAGTTCATTGC-3′). All primers designed specifically for this study, including LampR, CR1, and CASSFN, were designed using Primer3 (Untergrasser et al., 2012). Extension reactions were prepared using 1/4 reactions from the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, California, USA). For each 20 µl reaction, the following reagents were combined: 2 µl Ready Reaction Mix, 3 µl 5× sequencing buffer, 3.2 µl primer (1 µM), 1 µl DNA template (20–30 ng), and 10.8 µl water. Thermal cycler settings were set at an initial denaturation at 96 °C for 1 min followed by 30 cycles of 96 °C for 20 s, 50 °C for 20 s, and 60 °C for 4 min. Reactions were cleaned with 2 µl of 2.2% sodium dodecyl sulfate and returned to the thermal cycler at 98 °C for 5 min followed by 25 °C for 10 min. Reactions were then purified in Sephadex size-exclusion columns (GE Healthcare Life Sciences, Little Chalfont, UK) and 12 µl of purified products were loaded into an Applied Biosystems 3130 Genetic Analyzer for sequencing.

Sequences from individuals that successfully amplified at both non-coding regions (n = 54) were concatenated into a single sequence denoted by NC_total. Concatenated sequences were then trimmed to a 652-bp region to exclude repetitive portions of reads prone to slippage (NCI: excluded 14 bp before position 15382; NCII: excluded 28 bp after position 16173). Regions prone to slippage were identified by superimposed sequences on the chromatograms. All trimmed sequences were aligned with the MUSCLE algorithm (Edgar, 2004). Finally, a haplotype network was generated for all NC_total haplotypes to visualize their relationship using TCS v.1.18 (Clement, Posada & Crandall, 2000).

Second, to estimate population fluctuations based on nuclear microsatellite data, we used previously-published allele frequency data (Bryan et al., 2005). Bryan et al. (2005) genotyped individuals from Great Chazy River (n = 40) and Lewis Creek (n = 40) at 8 microsatellite loci. There was no evidence for linkage disequilibrium or significant deviations from Hardy–Weinberg equilibrium (see Bryan et al., 2005 for more details).

Moment-based analyses of historical population changes

Historical population fluctuations were inferred using three moment-based methods. First, the mismatch distribution of pairwise differences between NC_total haplotypes was plotted. The observed distribution of pairwise differences was compared to the expected number of pairwise differences under a model of population expansion with 1,000 bootstrap replicates in Arlequin v.3.1 (Excoffier, Laval & Schneider, 2005). Second, Harpending’s raggedness index (r) was used to test whether the observed distribution was significantly different from the expected theoretical distribution under a model of expansion (Harpending, 1994). Third, an alternative metric, Fu’s F_s statistic, was used to test the selective neutrality of mutations (Fu, 1997). Fu’s F_s can detect an excess (or deficiency) of haplotypes, given the observed haplotype diversity, thereby indicating a population expansion (or contraction).

Coalescent analyses of historical population changes

Two Bayesian coalescent MCMC models were used to estimate historical demographic fluctuations of P. marinus in Lake Champlain over time. First, the program BEAST v.1.6.2 (Drummond & Rambaut, 2007) was used to make inferences based on NCII sequences. This region was selected because of its high concentration of polymorphic sites. BEAST applies a Bayesian coalescent-based procedure, using MCMC to sample the posterior distribution of genealogical trees, demographic parameters over time, and coalescent events given sequence information and a set of priors. For our demographic model, we applied the Bayesian Skyline Plot (BSP). BSP is a change-point model that, assuming a single panmictic population, estimates fluxes in population size through time and uses a smoothing procedure to visualize these changes (Drummond et al., 2005). To determine which nucleotide substitution model fit the data, Akaike information criteria (AIC) values were calculated in jMODELTEST v.0.1.1 (Guindon & Gascuel, 2003). The substitution model selected by AIC was then used as a prior in BEAST. For our baseline study, we assumed a strict molecular clock of 3.6 × 10⁻⁸ substitutions per base per year (i.e., 3.6% substitutions per million years), based on previous estimates of divergence in mtDNA noncoding regions in fishes (Donaldson & Wilson, 1999). The BSP group number (m) was set to 15. The parameter m allows adjacent coalescent intervals to be grouped so that they can have the same N_e—it serves to smooth the resulting BSP. The maximum effective population size was set at 10,000, a high estimate based on preliminary coalescent model runs. MCMC chains were run for 50⁶ iterations in triplicate, sampling the posterior distribution every 1,000 iterations. The sampling distribution of the model was evaluated in TRACER v.1.5, with the first 10% discarded as burn-in (Rambaut et al., 2014). Quality of the MCMC convergence was assessed by the effective sample sizes (ESS): if the ESS value was less than 100, it was assumed that the MCMC chain had not been run long enough to get an accurate representation of the posterior distribution and the trace was discarded (Drummond et al., 2007).

A sensitivity analysis for the BEAST modeling was carried out on two parameters—clock rate and maximum effective population size—to assess whether priors biased parameter estimates. These parameters were selected because the priors were based on our preliminary estimates. For each parameter change, three additional MCMC chains were run keeping all other baseline input values constant. We used two alternative maximum population sizes (20,000 and 100,000) and one alternate clock rate (2.0 × 10⁻⁸ substitutions per base per year, i.e., 2.0% substitutions per million years). To determine how sensitive model results were to the priors, runs with altered priors were compared.

Historical demographic changes were also inferred from coalescent modeling of nuclear microsatellite data. Allele frequencies at eight microsatellite loci, previously published by Bryan et al. (2005), were used as input for the Bayesian MCMC model Msvar v.1.3 (Beaumont, 1999). Msvar uses probable genealogies of allele frequency data to generate posterior probability distributions of four demographic parameters: current effective population size (N₀), historical effective population size (N₁), mutation rate (μ), and time since the demographic change began (t) (Storz & Beaumont, 2002; Beaumont & Rannala, 2004). Broad prior distributions were defined for each parameter to test whether the model could detect true population fluctuations (Table S1). Five independent chains were run for a panmictic sea lamprey population (Bryan et al., 2005; Waldman, Grunwald & Wirgin, 2006) under a model of exponential growth, with an average generation time of six years (Hardisty & Potter, 1971), and a N_e/N_c ratio of 0.2 (Frankham, 1995). Each chain consisted of 8 × 10⁸ iterations, with sample points taken once every 10,000 iterations. The burn-in portion of the chain was excluded by truncating runs to consider only the second half. To estimate each demographic parameter, the five truncated chains were combined into one posterior distribution and peak density values were recorded for each parameter, along with 95% highest probability density (HPD) credibility intervals. A Bayes factor was also calculated to test whether there was more support for a population expansion or contraction using the method described by Storz & Beaumont (2002). Under the Bayesian statistical framework, a Bayes factor is an appropriate metric for comparing alternative models.

Sensitivity analyses were not conducted for Msvar because the model is computationally demanding; however, simulation modeling has revealed that Msvar has the power to detect true population contractions even with broad priors (Girod et al., 2011). These simulation models have also shown that the precision of estimates increases when demographic parameters are scaled using coalescent theory. Thus the marginal posterior distributions of the following scaled parameters were also plotted: θ₀ = 4N₀μ and θ₁ = 4N₁μ (effective population sizes scaled by mutation rate); t_f = t/(2N₀) (time scaled by current effective population size).

Haplotype diversity & population structure

We observed 14 unique haplotypes among the 54 NC_Total mitochondrial sequences (Table 1). The haplotype network of all 14 NC_total haplotypes visualizes their relationships as well as their relative abundance (Fig. 1). Haplotype 1 is predominant in Lake Champlain, with five other rare haplotypes (2, 3, 6, 7, and 9) having only one base pair difference from haplotype 1. Most of the mutations within these five haplotypes are point mutations within NCI or the 5′ end of NCII (Table 1). In contrast, most of the remaining rare haplotypes branch off of haplotype 7 and are characterized by mutations within the A/T-rich, repetitive 3′ region of NCII. In general, these mutations in the repetitive region were retained in each subsequent repeat of the sequence, leading to a larger number of overall mutations in NCII.

An exact test of population differentiation revealed that haplotypes were randomly distributed across the three sampling locations (n = 54; p = 0.265; # dememorization steps = 10,000; # steps Markov chain = 100,000). Thus, the Lake Champlain population is considered to be panmictic, consistent with results from previous studies (Bryan et al., 2005; Waldman, Grunwald & Wirgin, 2006).

Moment-based analyses of historical population changes

The mismatch distribution of NC_total was smooth and unimodal, suggesting that the sea lamprey population has undergone an expansion (Fig. 2). The distribution of observed pairwise differences closely matches the expected distribution of pairwise differences under a model of population expansion, with a raggedness index that was positive, but not significant (r = 0.04; n = 54; p = 0.86). This indicates that there was no significant deviation from the theoretical model of expansion. Additionally, Fu’s F_s statistic was significantly negative (F_s = − 6.61; n = 54; p = 0.02), indicating an excess of rare haplotypes, which would be predicted under a scenario of a recent population expansion. Overall, these moment-based methods of studying historical demography with mtDNA data strongly support a population expansion, without explicit estimates of the timing and magnitude of the event.

Coalescent analyses of historical population changes

For the BEAST analysis of mtDNA NCII data, the best-fit model of sequence evolution determined by jMODELTEST was Hasegawa, Kishino and Yano (HKY) + I, where I means there is a significant proportion of invariable sites. The BSP derived from the NCII data shows a decrease in effective population size starting around 400 years ago that continued until 50 to 100 years ago, at which point effective population size slowly began to increase (Fig. 3). The sensitivity analysis showed that the maximum population size prior was directly proportional to the current N_e estimate; when doubled, the N_e estimates also doubled. However, changes in the maximum population size and clock rate priors did not affect the overall pattern of the BSP (Table S2).

In contrast to the BEAST analysis, coalescent modeling of microsatellite alleles estimated one historical population contraction. The five Msvar chains converged below a critical potential scale reduction factor of 1.2 (Brooks & Gelman, 1998) for all four demographic parameters, indicating that there was good convergence of parameter estimates. The Bayes factor (BF) for a population contraction indicated very strong support for a population contraction over a population expansion (2^∗ln(BF) = 10.3) (Kass & Raftery, 1995), with the vast majority of MCMC iterations estimating a population contraction (Fig. 4A). The posterior density distributions reveal the estimates for each demographic parameter with 95% HPD intervals (Figs. 4B–4D). There was a decrease in effective population size from approximately 2,660 (HPD_95% = 153, 65750) down to approximately 50 currently (HPD_95% = 0.06, 1440) (Fig. 4B). The density peak was higher, with narrower credibility intervals, for historical effective population size compared to current effective population size. This population contraction was estimated to have begun approximately 820 years ago, although the peak density was only 40% for this time parameter, with wide credibility intervals (HPD_95% = 0.81, 71558) (Fig. 4C). Finally, there was a high density peak for the mutation rate estimate, with over 80% of the estimates approaching a modal value of 3.26 × 10⁻⁴mutations/site/generation (HPD_95% = 3.52 × 10⁻⁵, 2.78 × 10⁻³) (Fig. 4D).

The posterior distributions of the scaled demographic parameters also indicate a population contraction (Fig. 5). As predicted by simulation analyses, scaling the parameters increased the precision of parameter estimates. The magnitude of the scaled population size estimates were very similar to the unscaled estimates. Modal θ estimates corresponded with a decline from approximately 2,335 individuals historically (HPD_95% = 359, 29614) to 60 individuals currently (HPD_95% = 0.08, 772). Despite increased precision, there was still some overlap in the 95% credibility intervals for θ (Fig. 5A). In contrast, there was a substantial increase in precision of the time estimate, suggesting a more ancient demographic event: the modal value of the scaled time distribution corresponded with an event beginning approximately 1,230 years ago (HPD_95% = 128, 7522) (Fig. 5B).

Moment-based analyses of historical population changes

Mismatch distribution analysis is based on the assumption that demographic changes leave corresponding genetic signatures in neutral sequence data and gene trees (Rogers & Harpending, 1992; Harpending et al., 1998). The mismatch distribution of this study fits a model of population expansion well, being both smooth and unimodal. However, a major weakness of this method is its poor ability to make specific inferences about the timing of demographic changes. In general, when using moment-based methods, only a general inference as to whether the change was relatively recent or ancient can be made. The presence of many low-frequency mutations is one indication of a ‘recent’ expansion (Schneider & Excoffier, 1999). This pattern was evident in the haplotype network, in which the majority of haplotypes in Lake Champlain were present in only one to three individuals (Fig. 1). Additionally, the mean number of pairwise differences in the mismatch distribution can be used for a crude estimate of timing: a ‘low’ mean can indicate a ‘recent’ expansion while a ‘high’ mean can indicate a more ‘ancient’ expansion (Okello et al., 2005). The haplotypes in the Lake Champlain lamprey population had a mean of 1.74 mismatches, which supports a ‘recent’ expansion. Yet, without explicit guidelines for differentiating between ‘low’ and ‘high,’ or inferring what ‘recent’ and ‘ancient’ indicate on an evolutionary timescale, caution should be used in interpreting these results. Instead, these results are more appropriately used as supplemental evidence to be considered alongside inferences made by coalescent-based analyses.

Coalescent analyses of historical population changes

BEAST and Msvar analyses both detected an initial decline in effective population size. These congruent results between two separate genomes strengthen the evidence for a decline; however, the magnitude and timing of the event(s) remain uncertain. The BSP generated using BEAST provides credibility intervals at every time point that reflect both coalescent and phylogenetic uncertainty. This proved to be important, as the BSP shows that the credibility intervals for effective population size ranged from 51 to 9,500. These wide credibility intervals around N_e could be due to the limited number of segregating sites within the NCII sequence data. The sensitivity analysis also showed that the maximum population size prior affected the BSP results, with the estimated N_e being proportional to the prior. However, the overall pattern of demographic change was consistent regardless of the priors (see Table S2), indicating that the estimated mean N_e always showed a contraction followed by an expansion.

In contrast to BEAST, Msvar estimates four parameters simultaneously and can only detect the single most likely demographic trend over time based on the strongest genetic signal. Bayesian models, such as Msvar, require prior distributions for each parameter (Storz & Beaumont, 2002; Beaumont & Rannala, 2004). Generating priors can be difficult when no accurate estimates of demographic parameters are available, as is the case for sea lamprey in Lake Champlain. Estimates of current sea lamprey population size are poor because census data are inferred from wounding rates on prey species. Moreover, no historical documentation of sea lamprey population size in Lake Champlain exists. In the absence of reliable data, broad and equivalent priors were used for current and historical effective population size (Goossens et al., 2006). Furthermore, through simulation analyses, Girod et al. (2011) showed that Msvar is powerful at detecting population fluctuations even without informed priors if the demographic event was sufficiently large and ancient. The contraction detected in this analysis indicates that effective population size has declined from roughly 2,660 to 50 (a low but plausible estimate given the aggressive population control program and the fish’s high fecundity). It is important to note that these values represent best estimates based on the mode of the posterior density distributions, but that the 95% credibility intervals of historical and current population size overlap due to the uncertainty associated with the current effective population size parameter. Scaling the population size parameters by the mutation rate did increase precision, but there was still overlap in the credibility intervals. Thus, for a comprehensive population history, the Msvar results must be interpreted alongside results from the BEAST and moment-based analyses.

Scenario 1: evidence for fluctuations within a native population

Capitalizing on multiple analytical approaches and two unique genetic data sets used in this study, it is possible to infer multiple demographic fluctuations. While the mismatch distribution detected a population expansion and coalescent modeling using Msvar detected a decline, the results may not be in conflict. Instead, they may be detecting the two different demographic changes that are both evident in the BSP generated by BEAST. By considering known events in the ecological history of Lake Champlain and the approximate time scale estimated by each method, we present the most parsimonious demographic history of a native sea lamprey population in Lake Champlain.

The decline in effective population size detected by both coalescent models may correlate with land-use changes and fishing pressures that began after the arrival of European settlers in the mid-18th century. The BSP from BEAST indicated that the decline began approximately 400 years ago. Given that there is uncertainty around this time estimate, one hypothesis is that a gradual decline in population size began around 1753 when a large number of European settlers arrived in the area. At this time, human impact on the landscape increased with the onset of large-scale agricultural and clear-cutting practices and, later, the construction of mills and dams that would have limited upstream lamprey spawning migrations (Klyza & Trombulak, 1999). Moreover, the acceleration of the decline shown in the BSP correlates with the extirpation of salmonids in Lake Champlain—the primary food source of parasitic-phase sea lamprey. The last documentation of native Atlantic salmon in the basin was 1830, with lake trout extirpation following shortly thereafter in the 1890s (Fisheries Technical Committee, 2009). Thus, it is plausible that the population contraction detected by both coalescent models can be explained by known ecological changes in the region, though some caution is warranted given that the scaled Msvar distributions suggest that the decline was more ancient (ca. 1,230 years ago).

The mismatch distribution analysis may be detecting a more recent population expansion associated with salmonid stocking in Lake Champlain. This recent expansion is supported by the upward trend in effective population detected in the BSP in the last 50 years. A stocking program began for both Atlantic salmon and lake trout in 1972 and now focuses on stocking the main part of the lake with yearlings (Marsden et al., 2003; Fisheries Technical Committee, 2009). The reintroduction of their primary food source could have allowed the sea lamprey population to expand.

While this demographic history is plausible given the ecological context of Lake Champlain over the last 300 years, it is important to consider whether the data represent two real population changes, or whether these changes are artifacts of the models. To address this issue, an important consideration is the difference in mutation rates between the nuclear and mitochondrial markers used in the study, as faster-evolving markers generally have a greater power to detect more recent events. While the mitochondrial genome tends to have a higher mutation rate than the nuclear genome due to the mutagenic properties of respiration by-products and the limited DNA repair mechanisms of the mitochondrial genome (White et al., 2008), nuclear microsatellites have an even higher average mutation rate due to strand slippage. Indeed, Msvar estimated a fast mutation rate for the microsatellites used in the study, with the narrow posterior distribution peaking at 3.24 × 10⁻⁴ mutations/site/generation. Nevertheless, Msvar, a model that uses microsatellite markers exclusively, has been shown to have a bias towards detecting ancient declines in effective population size (Beaumont, 1999). Most recently, this tendency has been rigorously confirmed through simulation studies (Girod et al., 2011). Therefore, even if there were signatures of two demographic fluctuations in the population’s history, Msvar, which can only identify a single event, is predicted to detect the more ancient decline. As to whether or not the mtDNA sequence data show a true signature of expansion within the last 40 years, some insight can be gained from a comparison between BEAST model inputs and outputs. Even though an extremely strict (i.e., slow) molecular clock was used as an input in this analysis, the posterior distribution of effective population size still shows an upward trend beginning approximately 50 years ago. Posterior estimates that deviate from prior predictions generally indicate a strong genetic signal. Taken together, these lines of evidence suggest that the population fluctuations detected are likely to be real.

Scenario 2: evidence for a founder event in the early 20th century

Differentiating between founder events resulting from an invasion versus bottlenecks within native populations can be challenging, as both events are predicted to lead to a decline in effective population size and a loss of genetic diversity (Nei, Maruyama & Chakraborty, 1975). Thus, an alternative interpretation of the data is that the population contraction detected in Msvar represents a bottleneck immediately following an invasion from the anadromous Atlantic Coast population sometime after the completion of the Champlain Barge Canal in 1916, while the expansion detected in the moment-based and BEAST analyses represents the population boom in the 1970s. While the timing of such an invasion falls within the bounds of the BEAST and unscaled Msvar credibility intervals, both coalescent models assume a single panmictic population. Therefore, if the present-day Lake Champlain population is significantly differentiated from the anadromous population, the coalescent models could not accurately estimate historical changes beginning in the differentiated population of origin. Indeed, Bryan et al. (2005) found evidence for significant genetic structure between Lake Champlain and anadromous sea lamprey populations, and long-term vicariance for the Lake Champlain population.

The mitochondrial haplotype network can provide additional insight into founder events: the random sampling associated with a recent founder event should lead to the loss of rare haplotypes. Therefore, a recent founder event is unlikely given the observed network, which features a number of rare haplotypes branching off of a predominant ancestral haplotype (Fig. 1). In sum, the majority of the evidence from the genetic models suggests that the population is native to Lake Champlain, but uncertainty remains and a founder event associated with an invasion in the early 1900s cannot be completely excluded given the genetic data that are presently available.