Genetic compositions and comparative diversity of zebra (Dreissena polymorpha) and quagga (D. rostriformis) mussel populations are compared across their three decade-long histories in the Hudson River and Lake Erie. We analyzed 15 nuclear DNA microsatellite loci for the zebra mussel and 10 for the quagga mussel. Results indicate that Hudson River and Lake Erie zebra mussel populations slightly diverge in genetic compositions, and have similar overall genetic diversity levels. The allelic composition of the Hudson River zebra mussel population significantly changed during the middle time period (2003), suggesting genetic replacement. Its temporal genetic diversity levels stayed similar. In contrast, the Hudson River’s quagga mussel population remained genetically consistent over time in both composition and diversity. The Lake Erie zebra mussel populations displayed slight changes in allelic composition and increased in genetic diversity from the earliest timepoint, suggesting supplementation with new propagules. Populations of the quagga mussel remained genetically consistent over time in Lake Erie. The zebra mussel veliger larvae sample in 2016 differed from the adult samples, attributed to gene flow from other areas resulting in genetic admixture. Findings indicate that invasive populations may undergo significant genetic divergence or remain consistent over time, whose patterns may differ across their ranges and between related species. The population dynamics underlying their invasional success thus may be complex.

Sampling

Dreissenid mussels were collected with collaborators (see Acknowledgements) from the mid-river section of the Hudson River (“HR”, Stuyvesant, NY; latitude: 42.39^oN, longitude: -73.78^oW) and from South Bass Island in western Lake Erie (“LE”; 41.63 ^oN, -82.84 ^oW; Fig. 1) across three time periods. These constituted early (E), middle (M), and later (L) stages in their invasions, from which a total of 319 zebra and 202 quagga mussel individuals were analyzed (Table 1). The zebra mussel sample timepoints spanned from 1994 (E), 2003 (M), and 2016 (L) in both Lake Erie and Hudson River. The quagga mussel sample timepoints in Lake Erie spanned from 1998 (E), 2011 (M), and 2016 (L), while it’s invasion within Hudson River is more recent (~2008), and thus quagga mussel sample timepoints were limited from 2010 (E) to 2016 (L) in Hudson River. Whole dreissenid mussels were preserved in 95% ethanol (EtOH) in the field, labeled, transported to the laboratory (where EtOH was changed), and archived at room temperature. For all the historic samples, tissue was stored in 95% EtOH and replaced annually until DNA was extracted and analyzed with microsatellites in this study. A plankton sample containing veliger larvae was collected from LE in July 2016 with a 63µm Wisconsin Cole-Palmer 40-AD50 plankton net (https://www.coleparmer.com) attached to a 15.2m line. The net was thrown from the seawall/shore and upon sinking, was maintained at the maximum depth the line would allow for 60s, then manually retrieved using a hand-over-hand technique. Upon retrieval, the sample was placed into a15mL falcon tube containing 95% EtOH. Cross-polarized light microscopy was used to identify dreissenid mussel veliger larvae from the plankton following Johnson (1995). Individual larvae were isolated in 0.5mL tubes with 95% EtOH for later DNA analysis.

DNA extractions and gene amplifications

Genomic DNA was extracted and purified from adductor muscle and mantle tissues of adults, and from whole individual veliger larvae, with DNeasy^® Blood and Tissue kits (Qiagen Inc., Valencia, CA, USA) following the manufacturer’s protocol. Half of the reagents were used for the larvae. DNA quality and quantity were evaluated with a Nanodrop^™ 2000 spectrophotometer (Thermo Fisher Scientific^™ Inc., Waltham, MA, USA), visualized on 1% agarose mini-gels with ethidium bromide, and stored at 4°C prior to amplification with the polymerase chain reaction (PCR). A one-step multiplex dreissenid-specific PCR assay targeting the mtDNA 16S gene for zebra mussel and the mtDNA cytochrome oxidase subunit I (COI) gene for quagga mussel (Ram et al. 2011) was used to identify each individual to their species (both adults and veligers). Each individual then was genotyped at 15 nDNA μsat loci for zebra mussels (Table S1A) and 10 nDNA μsat loci for quagga mussels (Table S1B).

Microsatellite PCRs were run in 10μL reactions, containing 0.035 units AmpliTaq^® DNA polymerase (ABI; Applied Biosystems^™, Foster City, CA, USA), 1X GeneAmp^® PCR Buffer I (ABI), 100μM dNTPs, 2.8mM MgCl₂, 0.5μM each of the forward and reverse primers, and ~50ng extracted DNA on a C1000^™ Thermal Cycler (Bio-Rad Laboratories, Hercules, CA). Positive (known genotype) and negative (no template) controls were included in each run. Amplification profiles were 3 min initial denaturation at 95°C, followed by 34 cycles of 30 sec at 95°C for denaturation, 30 sec for annealing (temperature in Table S1), and 1 min extension at 72°C, capped by 5 min final extension at 72°C. Amplifications were conducted separately for each locus. Products were diluted 1:50 with ddH₂O, of which 2μL was added to 13μL of a solution containing formamide and ABI Genescan^™-500 LIZ^® size standard, and loaded into single wells on 96-well plates. Samples were denatured for 2 min at 95°C and analyzed on an ABI 3130x1 Genetic Analyzer with GeneMapper^® 4.0 software (ABI). Output profiles were manually checked to confirm allelic size variants.

Genetic data analyses

For each locus in each sample, possible null alleles, scoring errors, and/or large allele dropout were evaluated using Micro-Checker v2.2.3 (van Oosterhout et al. 2004). Temporal samples were tested for conformance to Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) expectations using Genepop v4.3 (Rousset 2008). Heterozygote deficits are common in many mollusks (Li et al. 2003; Plutchak et al. 2006), especially dreissenid mussels (Marsden et al. 1995; Wilson et al. 1999; Elderkin et al. 2001; Astanei et al. 2005; Therriault et al. 2005; Gosling et al. 2008; Imo et al. 2010; Feldheim et al. 2011; Marescaux et al. 2016; Mallez and McCartney 2017; Lindsay et al. 2018), which may be attributed to the Wahlund effect (Marescaux et al. 2016). Therefore, we used FreeNA to correct for null alleles and heterozygote deficiency within our population samples. Frequencies of null alleles were computed with FreeNA (Chapuis and Estoup 2007).  

Numbers of alleles (N_A) and observed and expected heterozygosities (H_O/H_S) were calculated in Arlequin v3.5 (Excoffier and Lischer 2010) and allelic richness (A_R) in Fstat v2.9.3.2 (Goudet 1995). Genetic diversity differences measured in allelic count among populations was analyzed with permutation test using the R package FPTest (Yang and Fu 2017), with significance evaluated after sequential Bonferroni correction (Rice 1989). Genepop was used to calculate inbreeding index (F_IS). Numbers of private alleles (N_PA), i.e., those found exclusively at a single collection site or in a single temporal sample at a given site, were determined with Convert v1.31 (Glaubitz 2004). Relative proportions of private alleles (P_PA) were calculated by dividing the number of private alleles by the total number of alleles in that sample at all loci. Mean frequency of private alleles (F_PA±SE) was calculated using Fstat. Colony v2.0.6.1 (Jones and Wang 2010) was used to test for full and half siblings in each sample, and for parentage analysis between the adults and veliger larvae collected from Lake Erie in 2016. NeEstimator v2.01 (Do et al. 2014) was used to estimate effective population size (N_e) using (1) a single-sample method from the linkage disequilibrium (LDN_e) described by Waples and Do (2008), and (2) and a temporal two-sample method for each population location. Jackknife methods were used to calculate 95% confidence intervals (CI) for all N_e estimates.

To discern whether samples significantly diverged over time or space, Weir and Cockerham (1984) estimates of F_ST were calculated and corrected for null alleles with FreeNA. Exact G-tests were performed in Genepop with significance evaluated after sequential Bonferroni correction. Analysis of Molecular Variance (AMOVA) in Arlequin examined hierarchical partitioning of variation between populations and their component temporal samples. Additionally, three-dimensional Factorial Correspondence Analysis (3dFCA) in Genetix v4.05 (Belkhir et al. 2004) and Multivariate Discriminant Analysis of Principal Components (DAPC; Jombart et al. 2010) using the R package adegenet (Jombart et al. 2018) were used to visualize spatial and temporal trends in population structure.

Bayesian assignment tests in Structure v2.3.4 (Pritchard et al. 2000) analyzed the likelihood of individuals belonging to K hypothetical groups (1–8 for zebra mussels and 1–6 for quagga mussels), with burn-ins of 50,000, 100,000 replicates, and 10 iterations each. As recommended by Puechmaille (2016) and Janes et al. (2017), K values were evaluated using multiple methods (i.e., ΔK, LnP(D), Medmedk, Medmeak, Maxmedk, and Maxmeak) in StructureSelector (Li and Liu 2018). Geneclass2 (Piry et al. 2004) evaluated assignments of individuals to populations, using the “enable probability computation”, default parameters, and 100,000 simulated individuals. A hierarchical approach was taken, with tests run 1) with each temporal sample separate, and 2) with temporal samples grouped for spatial populations.

Dpolymorpha-ARLEQUIN.arp – Microsattelite genotype data in Arlequin input format for Driessenia polymorpha (zebra mussel).

Drostriformis-ARLEQUIN.arp – Microsattelite genotype data in Arlequin input format for Driessenia rostriformis (quagga mussel).