Genetic trajectories of zebra and quagga mussel invasions across three decades: Lake Erie versus Hudson River populations

Genetic compositions and comparative diversity of zebra (Dreissena polymorpha) and quagga (D. rostriformis) mussel populations are compared across their three decade-long histories as invasive species in the Hudson River and Lake Erie of North America. We analyze 15 nuclear DNA microsatellite loci for the zebra mussel and 10 for the quagga mussel. Results indicate that the Hudson River and Lake Erie zebra mussel populations slightly diverge in genetic compositions, and possess similar overall genetic diversity levels. The allelic composition of the Hudson River zebra mussel population significantly changed during the middle time period (2003) analyzed, suggesting genetic replacement. Yet, its overall levels of genetic diversity levels have stayed similar. In contrast, the Hudson River’s quagga mussel population has remained genetically consistent over time in both composition and diversity. Lake Erie’s zebra mussel population underwent slight change in allelic composition and increased in genetic diversity from the earliest timepoint, suggesting allelic supplementation from newly arriving propagules. In contrast, Lake Erie’s quagga mussel population has remained genetically consistent over time. The genetic composition of Lake Erie zebra mussel veliger larvae sampled in 2016 differed from its adult samples, attributable to gene flow from other areas and genetic admixture. Overall 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 successes thus may be complex.


Introduction
Some of the most invaded North American aquatic habitats have included the Hudson River system (> 122 invasive species; Mills et al. 1996;Strayer 2012) and the Laurentian Great Lakes (188 invasive species; NOAA 2019; Sturtevant et al. 2019). In the Great Lakes, the majority of these invasive species originated from Eurasia (> 100) and most were associated with ballast water releases from commercial shipping vessels (~ 70) (Sturtevant OPEN ACCESS. et al. 2019), resulting in significant ecological and economic impacts (Kelly et al. 2009). Two of the most "notorious" introductions have been the Eurasian zebra mussel Dreissena polymorpha (Pallas, 1771) and the quagga mussel D. rostriformis (Andrusov, 1897), due to their tremendous fouling and filtering capabilities, transformation of soft to hard benthos, and ecological community changes (Berg et al. 1996;Vanderploeg et al. 2002;Stewart et al. 2009;Marshall and Stepien 2019).
Lake Erie, the shallowest and most eutrophic of the Great Lakes, has the longest history of Dreissena invasion in North America dating to the zebra mussel's first appearance in 1986 (Hebert et al. 1989;Carlton 2008). The zebra mussel was followed by its close relative, the quagga mussel, in 1989 (Mills et al. 1993). The zebra mussel populations initially grew very quickly, while the quagga mussel remained rare for a decade. Zebra mussels reached maximum density by 1992, and then cycled through density-dependent biomass fluctuations (Karatayev et al. 1997(Karatayev et al. , 2002. Quagga mussels experienced a longer lag time, not reaching maximum population densities until 1998 (Karatayev et al. 2014). Quagga mussels have outcompeted zebra mussels throughout the central and eastern Lake Erie basins, although a mixed species community remains in the western basin, where shallower waters provide refugia for zebra mussels (Karatayev et al. 2014).
The zebra mussel was discovered in the Hudson River in 1991 , followed by the quagga mussel in 2008 (Strayer and Malcom 2013). By 1992, zebra mussels accounted for > 50% of the Hudson River's entire benthic biomass . Survival of adult zebra mussels initially was high (~ 50%/year) but steadily fell over time (< 1%/year since 2005), accompanied by decreasing mean body size (Strayer et al. 2011). Quagga mussels quickly spread throughout the Hudson's entire freshwater reach, but still remain < 10% of the entire Dreissena community (Strayer and Malcom 2013;Strayer et al. 2020). Quagga mussels tend to outcompete zebra mussels in areas where they co-occur (Mills et al. 1996;Karatayev et al. 2011), but not in all systems (Strayer et al. 2019). Such replacement can take a decade or longer (Karatayev et al. 2011), and might be in the early stages within the Hudson River.

Invasion genetics and hypotheses
Monitoring invasive populations across time is necessary for understanding their successes Pérez-Portela et al. 2012;Rius et al. 2015), although sample archives and long-term temporal datasets are rare and often difficult to obtain and maintain. Population genetic patterns relating to the relative success of invaders may follow a variety of pathways and alternatives for different species, possibly even differing among populations and across spatio-temporal scales. In brief, Stepien et al. (2018) outlined potential invasion genetic hypotheses, with the null hypothesis being that the genetic composition of an invasive population would be expected to remain relatively consistent over time (hypothesis I). The alternative hypothesis is that the invasive population may genetically change over time, which might be due to selection, drift, gene flow, and/or the arrival of new genotypes (hypothesis II). Under these hypotheses, three scenarios are possible: the overall level of genetic diversity (i.e., measures of heterozygosity, allelic richness, etc.) might remain statistically similar (A), decrease (B), or increase (C) over the time course of an invasion (Stepien et al. 2018).
Under scenario IA, a newly established invasion might experience high reproductive success and retain all or most of its genetic variation over time (Stepien et al. 2018). For example, invasive round goby Neogobius melanostomus (Pallas, 1814) populations in the Great Lakes that were founded from various Eurasian sources (Dillon and Stepien 2001;Stepien and Tumeo 2006;Brown and Stepien 2009) have remained relatively genetically consistent over 25 years (Snyder and Stepien 2017). Following scenario IIA, an invasion might change in allelic composition over time, yet possess approximately the same overall genetic diversity level. For example, the Eurasian ruffe Gymnocephalus cernua (Linnaeus, 1758) invasion in the Great Lakes displayed significant temporal changes in genetic composition over 30 years, yet had relatively consistent genetic diversity levels (Stepien et al. 2018). In scenario IIB, genotypic frequencies may change and diversity levels decline temporally due to small population sizes and drift (i.e., founder effects and bottlenecks). In scenario IIC, genotypic variation and diversity might increase if new propagules with additional alleles arrive, become established, and augment the population. In some cases, variation across the time course of an invasion might reflect a combination of these scenarios. For example, the spiny waterflea Bythotrephes longimanus (Leydig, 1860) initially exhibited a founder effect after invading the Great Lakes (IIB), which then was supplemented over time with new introductions from multiple sources (IIC; Berg et al. 2002).

Research objectives
The present study compares and contrasts the relative genetic trajectories and diversities of two concurrent invasions by closely-related congeneric species in two North American populations across their invasional time courses. We test the hypotheses of temporal (I) genetic stasis versus (II) changes in allelic composition, with the scenarios of overall genetic diversity levels (A) remaining similar, (B) declining (with drift), or (C) increasing (as might stem from arrival and establishment of new allelic variants) over time. We compare their differentiation in genetic composition (using pairwise and hierarchical comparisons, and Bayesian assignment tests), and genetic diversities (heterozygosity, allelic richness, and private alleles) at 15 nuclear DNA microsatellite (μsat) loci for zebra mussel and 10 for quagga mussel. Population fluctuations presumably have impacted the genetic compositions, diversity levels, and population structure of these dreissenid mussel species. Moreover, little is known about whether new introductions of their propagules have occurred during the last 30 years, possibly adding or replacing alleles. Interpretation of the genetic similarities, differences, and interactions between zebra and quagga mussels aims to advance understanding the population dynamics underlying invasional success.

Sampling
Dreissenid mussels were collected with collaborators (see Acknowledgements) from the mid-river section of the Hudson River ("HR", Stuyvesant, NY; latitude: 42.39°N; longitude: −73.78°W) and from South Bass Island in western Lake Erie ("LE"; 41.63°N; −82.84°W; Figure 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 S1). The zebra mussel samples spanned from 1994 (E), 2003 (M), and 2016 (L) in both Lake Erie (LE) and the Hudson River (HR). The quagga mussel samples in LE extended from 1998 (E), 2011 (M), and 2016 (L). Due to the more recent invasion of the quagga mussel within the Hudson River (invaded ~ 2008), the Hudson River samples included two timepoints, 2010 (E) and 2016 (L). Whole dreissenid mussels were preserved in 95% ethanol (EtOH) in the field, labeled, transported to the laboratory (where the EtOH was changed), and archived at room temperature. For all historic samples, tissues were stored in 95% EtOH and replaced annually, until their 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 that the line would allow for 60s, then manually retrieved using a hand-over-hand technique. Upon retrieval, the sample was placed in a 15 mL falcon tube containing 95% EtOH. Cross-polarized light microscopy was used to identify dreissenid mussel veliger larvae from the plankton sample, following Johnson (1995). Individual larvae then were isolated in 0.5 mL tubes containing 95% EtOH, until later DNA analysis.

DNA extractions and gene amplifications
Genomic DNA was extracted and purified from the 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. Reagent concentrations were halved 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). Onestep multiplex dreissenid species-specific PCR assays that respectively targeted the mtDNA 16S gene for the zebra mussel and the mtDNA cytochrome oxidase subunit I (COI) gene for the quagga mussel (Ram et al. 2011) were used to identify each individual to its species (both adults and veligers). Each individual then was genotyped at 15 nDNA μsat loci for zebra mussels (Supplementary material Table S2A) and 10 nDNA μsat loci for quagga mussels (Table S2B).
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.8 mM MgCl 2 , 0.5 μM each of the forward and reverse primers, and ~ 50 ng 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 S2), 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 2 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.
Numbers of alleles (N A ) and observed and expected heterozygosities (H O /H E ) 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 in allelic counts among populations were analyzed using permutation tests using the R package FPTest (Yang and Fu 2017), with significance evaluated after sequential Bonferroni correction (SBC; Rice 1989). Genepop was used to calculate the 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 frequencies of private alleles (F PA ± SE) were calculated using Fstat. Colony v2.0.6.1 (Jones and Wang 2010) was used to determine whether full or half siblings were present in samples, and for parentage analysis between the adults and veliger larvae from LE in 2016. NeEstimator v2.01 (Do et al. 2014) estimated effective population size (N e ) using (1) a single-sample method from linkage disequilibrium (LDN e ) by Waples and Do (2008), and (2) with a temporal two-sample method per population location. Jackknife methods were employed 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 SBC. 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) in the R package adegenet (Jombart et al. 2018) were used to visualize spatial and temporal trends of population structure.

Genetic diversity and composition
A total of 250 alleles occurred among 15 μsat loci for the zebra mussel samples (Table S2A), averaging 16.7 alleles per locus (ranging from 6-36). For quagga mussel samples, we identified 231 alleles among 10 μsat loci (Table S2B), averaging 23.1 per locus (ranging from 12-36). Although heterozygote deficits were indicated in population samples of both species (seven loci per species did not conform to HWE), those did not occur at any given locus among most/or all samples or across all loci in any sample(s).
Colony analyses identified few full siblings in zebra or quagga mussel samples, and some half siblings (Table S3). There were no noticeable patterns in these occurrences between populations or among sampling time periods. For zebra mussels, the HR's early temporal sample and LE's later temporal sample each contained a single sibling pair. The sole quagga mussel sample with a full sibling pair was in the early temporal sample from the HR. Analyses thus were run with and without the full sibling pairs (removing a single individual per sibling pair), which did not affect the results (< 0.001) for F ST or the assignment tests. Thus, results from the complete data sets are presented here.

Population comparisons between the Hudson River and Lake Erie
Genetic diversity values were found to be similar between the HR and LE populations for zebra mussels (HR: A R = 10.74 and H E = 0.70, LE: A R = 11.20 and H E = 0.73; Table S1A) and quagga mussels (HR: A R = 14.71 and H E = 0.85, LE: A R = 14.50 and H E = 0.85; Table S1B). The amount of inbreeding did not differ between the HR and LE for either species (Table S1A, B). LE populations of both species contained greater numbers of private alleles Table 1. Pairwise genetic divergences (F ST ; below the diagonal), and comparisons of allelic distributions (above diagonal), for (A) zebra mussel and (B) quagga mussel, based on 15 and 10 microsatellite loci, respectively. Significance: *0.01 < p < 0.05, **0.001 < p < 0.01, ***p < 0.001 after SBC, NS = p > 0.05 (or without * below diagonal). Dashed boxes denote temporal comparisons within one location.  Table S1A, B). Both zebra and quagga mussels showed slight yet significant genetic divergences between the HR and LE populations when the temporal samples were grouped (ZM: F ST = 0.005*, QM: 0.002*; Table 1A, B). Geneclass2 assignment tests depicted high respective self-assignments for populations of each species (i.e., LE to LE, and HR to HR; Table S4).

Temporal comparisons within populations
For the zebra mussel, no significant changes in overall genetic diversity levels occurred across time in the HR population, whereas the LE population underwent a significant increase in numbers of alleles (E-M: p < 0.001***; Table 1A). Both zebra mussel populations changed in their respective genetic compositions over time. This difference was greatest for the middle time period in the HR population (E-M: F ST = 0.037***; M-L: F ST = 0.037***; Table 1A; results from pairwise exact G-tests were equivalent). A small, yet significant genetic divergence occurred for the LE samples (E-M: F ST = 0.005*; M-L: F ST = 0.008***; Table 1A). Furthermore, hierarchical analyses in Structure and DAPC revealed their slight genetic divergence from the later temporal sample in LE (Figure 2A, Figure S1). Geneclass2 results showed appreciable self-assignments among the temporal samples of zebra mussels, with the greatest self-assignment being from the HR for the middle time period (0.78; Table S4A2). For the quagga mussel, no significant temporal genetic diversity changes were observed in the HR or LE populations (Table 1B). However, significant temporal genetic divergence occurred in the LE quagga mussel population (E-L: F ST = 0.008**; Table 1B). Hierarchical analyses in Structure and in DAPC displayed overall similar genetic structure throughout the time course of the LE invasion ( Figure 3A, Figure S2). Geneclass2 results showed appreciable self-assignments among the temporal samples of quagga mussels, with the greatest self-assignment being for the early time period in the HR (0.83; Table S4B2).
AMOVA analyses of the two population regions (HR and LE) revealed small yet significant partitioning of genetic variation across the temporal courses for both species (ZM: 1.75%, p < 0.001; QM: 1.02%, p < 0.001**; Table 2). Single sample estimates of effective population size based on linkage disequilibrium yielded large values for both species, with CI ranging to infinity, except in the HR early zebra mussel population (LDN e = 280.9, CI = 163.4-872.8; Table 3). The three temporal two-sample NeEstimator Structure results illustrating hierarchical genetic structure of quagga mussel populations, with membership to K clusters (each cluster is represented by a different color, and each vertical bar represents an individual). Structure initially was run with all samples in the dataset (K = 2 shown in grey and white), followed by subsequent runs for the Hudson River (K = 2 shown, green and rose) and Lake Erie (K = 2 shown, yellow and red). B. Three-dimensional factorial correspondence analysis (3dFCA) illustrating the mean relationships among quagga mussel populations from the Hudson River and Lake Erie. C. Multivariate Discriminant Analysis of Principal Components (DAPC) for quagga mussel, indicating all individuals. Circles = Hudson River and Triangles = Lake Erie. E = Early, M = Middle, L = Later (time periods).  Nei and Tajima (1981), F k of Pollak (1983), and F s of Jorde and Ryman (2007)) produced similar results, and thus just F k is reported here. The temporal two-sample estimates displayed the trend of decreasing F k over the time courses of the zebra and quagga mussel invasions alike in LE and HR (Table 3). Furthermore, the zebra mussel temporal F k in LE was 1.27-2.94X larger than in HR over the invasion courses. However, since the CI values tended to overlap, comparisons between these F k values are problematic.

Temporal comparisons between population regions
No temporal samples differed in allelic counts between the HR and LE population regions (Table 1A). However, all temporal samples revealed significant genetic divergences between the HR and LE populations, with the greatest divergence occurring during the middle timepoint (F ST = 0.036***; Table 1A). Quagga mussel populations retained similar genetic diversity values across all temporal samples for both LE and the HR (Tables 1B, S1). However, the HR and LE quagga mussel populations had significantly different genetic compositions at their later timepoints (F ST = 0.011***; Table 1B).
For both species, Bayesian Structure analyses varied in the numbers of K clusters that were best-supported, according to the analysis method used (i.e., ΔK, LnP(D), Medmedk, Medmeak, Maxmedk, and Maxmeak; see Methods). However, each displayed similar patterns of population sample clustering, and results for the ΔK method thus are reported here. The complete zebra mussel dataset supported K = 2, with the HR middle timepoint (colored white; Figure 2A) grouping separately from the other population samples (colored grey; Figure 2A). The 3dFCA and DAPC similarly depicted the zebra mussel HR middle temporal sample as being genetically distinct ( Figure 2B, C). The hierarchical Structure analysis showed this same temporal separation for the HR (K = 2), along with discerning further structuring in LE (K = 4), with the veliger larvae sample differing from the nearby sample of adults (Figure 2A).
The quagga mussel dataset supported K = 2 clusters in Structure, without clear separation among population samples ( Figure 3A). Likewise, the hierarchical Structure analysis revealed little differentiation within the HR (K = 2) or LE (K = 2; Figure 3A). The quagga mussel 3dFCA and DAPC similarly demonstrated little difference between samples, except for separation of the later temporal samples ( Figure 3B, C).

Divergences between life stages
The LE zebra mussel veliger larval sample significantly diverged from the adults in genetic composition (L-Vel: F ST = 0.015***; Table 1A). Furthermore, A R values were slightly, but not significantly, lower and possessed four private alleles (Tables S1, S2). The larvae sample displayed self-assignment (0.57), with its second greatest assignment to the LE 2016 adult population (0.20; Table S4A3). The Structure analysis also indicated its unique genetic composition (Figure 2A), along with the 3dFCA (Figure 2B, C).

Genetic diversity levels over the temporal courses of the invasions
Lower genetic diversity in invasive populations relative to their native sources is hypothesized to stem from founder effects and bottlenecks (Roman and Darling 2007;Handley et al. 2011;Blackburn et al. 2015). Our findings indicate that neither dreissenid mussel species experienced significant founder effects in either population (LE and HR) throughout the > 25 years of their invasions, which concurs with previous genetic analyses of their North American and European introductions (Brown and Stepien 2010;Stepien et al. 2013;Marescaux et al. 2016), as well as findings from some other invasive bivalve species (e.g., Limnoperna fortunei in South America (Duarte et al. 2018) and Asia (Ghabooli et al. 2013), Mytella charruana along the southeastern North American coast (Calazans et al. 2017), Rangia cuneata in the Baltic Sea (Voroshilova et al. 2018), and Semimytilus algosus in South Africa (Zeeman et al. 2020)). The high genetic diversity of North American dreissenid mussel populations is believed to stem from large and genetically variable introductions, which were enhanced by admixture from several Eurasian source populations (Stepien et al. 2005(Stepien et al. , 2013Brown and Stepien 2010). Their ballast water introductions likely involved huge numbers of individual propagules introduced from multiple ships and sources (Stepien et al. 2013). Large numbers of propagules from multiple sources are assumed to facilitate the introduction success, establishment, and spread of invasive species (Lockwood et al. 2005;Roman and Darling 2007;Simberloff 2009).
Here, populations of both dreissenid mussel species possessed high genetic diversity levels (H E and A R ), consistent with their general global invasional trends (Astanei et al. 2005;Therriault et al. 2005;Brown and Stepien 2010;Marescaux et al. 2016), including populations stemming from secondary spread (i.e., inland lakes; Mallez and McCartney 2017). Furthermore, the LE populations of both zebra and quagga mussels contained greater number of private alleles. This result is not surprising, since the Lake Erie/Lake St. Clair system was the initial site of introduction discerned for both species in North America (Mills et al. 1993;Carlton 2008); Lake Erie thus is likely to possess and have retained more rare alleles. Additionally, estimates of effective population size in the present study depict greater population sizes of both species in LE than in the HR.
Over the time course of invasions, it has been hypothesized that continued propagule pressure may increase genetic diversity (e.g., allelic richness and heterozygosity), by introducing new alleles to augment the population (Lockwood et al. 2005;Simberloff 2009). Our study found a significant initial temporal increase in allelic count values in the LE zebra mussel population (partially supporting hypothesis scenario C), but not for the quagga mussel population (supporting scenario A). The HR genetic diversity levels for both species have remained consistent over the time courses of their respective invasions (supporting scenario A). Although their genetic diversity has been stable, temporal estimates of their effective population sizes reveal decreasing populations over the time course at each location for both species. These results are congruent with known population dynamics for dreissenid mussels, which tend to reach abundances beyond their carrying capacity, and then subsequently decline (Karatayev et al. 2011).

Temporal divergences in allelic compositions
The HR's zebra mussel population has undergone significant changes in allelic composition over the past three decades. Divergences were most pronounced during the middle time period (2003), which significantly differed from both the earlier (1994) and later (2016) time periods for the HR population (supporting hypothesis II; illustrated in the Structure, DAPC, and AMOVA analyses). Furthermore, overall population genetic diversity levels in the HR have remained temporally consistent (supporting hypothesis scenario IIA). Employing a mostly different suite of microsatellite loci, Brown and Stepien (2010) likewise found discernable change in genetic composition from 1991 to 2003 at the HR site of Catskills, NY, located ~ 23 river km downstream from the present collection site. By incorporating further analysis of a later temporal sample, we found a second genetic transition shift towards more similarity with the earlier period, which both significantly diverged from the middle invasion stage. Population ecology studies of dreissenid mussels in the HR showed that larval recruitment typically has occurred in large year classes in ~ 2-4 year cycles, with annual population sizes fluctuating by an order of magnitude due to carryingcapacity related adult mortality (Strayer and Malcolm 2006;Strayer et al. 2020). Such cyclical patterns associated with population expansions and declines may provide opportunity for allelic contributions from more distant sources, resulting in allelic turnovers in the population. The population patterns of genetic replacement that were observed here match the ecological population trends of zebra mussels in the HR (Strayer and Malcolm 2006;Strayer et al. 2020).
In contrast, the genetic composition of the HR's quagga mussel population has varied over time only slightly, remaining consistent in overall genetic diversity (supporting hypothesis IA of genetic stasis). The initial HR quagga mussel invasion (2010) appeared genetically most similar to the middle temporal sample from LE (2011), suggesting that the Erie Canal connection might have provided a pathway for genetic contribution to the HR. However, additional suspected source populations would need to be investigated for verification. By 2009, the quagga mussel had replaced most of the zebra mussel population throughout the Erie Canal (Stewart 2013), and likely continued into the HR.
The quagga mussel commonly outcompetes the zebra mussel in many habitats, as has occurred across most of the Great Lakes, reaching proportions of > 80% abundance (Karatayev et al. 2014). Species replacement also has occurred in many European rivers, with the quagga mussel increasing by an estimated 26% per year in Germany and the Netherlands (Heiler et al. 2013). Species replacement of the zebra mussel by the quagga mussel sometimes takes many years (Mills et al. 1999;Stoeckmann 2003;Nalepa et al. 2009;Karatayev et al. 2011). The eventual predominance of the quagga mussel often further alters benthic and pelagic ecology (Nalepa et al. 2009;Karatayev et al. 2014), as well as influences the population genetic structure of both species (see Stepien et al. 2013).
The LE zebra mussel population that we examined underwent apparent small, yet significant temporal changes in genetic composition, whose number of alleles significantly increased between the early and middle temporal stages (supporting hypothesis IIC). Separation of the temporal samples in the hierarchical Structure and DAPC results depict genetic divergence occurring from the early to later time periods. This genetic divergence and increased genetic diversity suggest supplementation with new propagules (and their alleles) from other invaded areas or overseas, primarily between the early (1994) and middle (2003) time periods. In contrast, relative genetic stasis has characterized the population of LE quagga mussels over time (supporting hypothesis IA). The invasive marine charru mussel (Mytella charruana) similarly displayed no temporal changes in genetic composition from 2006 to 2010, despite high genetic diversity and continued propagule release (Calazans et al. 2017); that study, however, was only four years long, much shorter than ours. However, another Great Lakes species introduced via ballast water, the round goby, has retained relative temporal stasis in genetic composition and diversity alike over 25 years (Snyder and Stepien 2017).
The zebra mussel population in western LE fluctuated from > 80% of the dreissenid biomass in 1998, to < 20% in 2001 (Stoeckmann 2003), and has experienced long-term population declines over the past decade, primarily due to competition with the quagga mussel (Strayer et al. 2019). Our temporal estimates depict this decreasing effective population size over time, which is congruent with its known population dynamics (Karatayev et al. 2011;Strayer et al. 2019). The present findings indicate that the zebra mussel's genetic composition has accrued small changes throughout the course of its invasion, while maintaining relatively consistent levels of overall genetic diversity since 2003 (supporting hypothesis IIA from the middle to later timepoints). The steady levels of genetic diversity since 2003 might have stemmed from reduction in arrivals of overseas propagules due to increased ballast water regulations implemented by Canada (Canada Shipping Act 2006) and the U.S. Coast Guard (2004Guard ( , 2018, to stem invasions (Bailey et al. 2011). New invasions in the Great Lakes have declined significantly following the "No Ballast on Board" regulations (Sturtevant et al. 2019).

Genetic relationships between life stages
The LE veliger larvae sample from 2016 had its second greatest assignment (Geneclass2 results) to its 2016 LE adults; however, the two samples significantly diverged in the hierarchical Structure and DAPC plots, as well as with F ST . Similarly, Haag and Garton (1995) found significant genetic divergence between zebra mussel larvae and adult samples in western LE (from the same location as our LE samples), analyzing variation at a single allozyme locus. Those researchers attributed the divergence to possible selection pressures on the different life-stages, although they might have mistakenly sampled both zebra and quagga mussel veligers, since the quagga was cryptically present at the time of their study. The possibility of species misidentification was eliminated in our investigation, since we used a molecular assay (Ram et al. 2011) to accurately determine the species of each individual adult and larva. It is important to note that each female mussel may spawn > 1,000,000 eggs (Borcherding 1991), and thus there is a magnitude greater abundance of larvae compared to adults. Furthermore, dreissenid mussels may undergo multiple spawning events over a reproductive season (Haag and Garton 1992;Ram et al. 1996), and thus, there may be appreciable intraspecific genotypic variation related to spatiotemporal patterns of their spawning across the lake. Thus, results from a single-snapshot with a limited sample size (e.g., 40 individual larvae) may be biased by the underlying substructure and stochasticity. These findings need to be further examined by incorporating sampling throughout the reproductive season, with genotypic comparisons made between the species and across their life stages, while incorporating metapopulation sampling.

Conclusions
The present study discerned different population genetic patterns for two sympatric conspecific invasions across a three-decade span, which both were characterized by high diversity levels. The HR's zebra mussel population underwent significant change in genetic composition over time, while maintaining its overall net level of variation, following hypothesis IIA, whereas its quagga mussel population has remained consistent in both composition and diversity (hypothesis IA). The LE zebra mussel population significantly changed over time, and initially increased in genetic diversity (hypothesis IIC), which then declined through the later timepoint (hypothesis IIA). The quagga mussel LE population sampled has remained in genetic stasis throughout its time course (hypothesis IA). The genetic composition of the veliger larvae population sample in LE significantly diverged from its adult counterpart, suggesting dispersal and gene flow from other areas. To fully understand the invasion dynamics and evolutionary patterns of zebra and quagga mussels, it is important to jointly analyze their demographic and genetic compositions in concert over time. This is especially critical as both species continue to expand their invasive ranges (e.g., Naranjo-García and Castillo-Rodríguez 2017; Prié and Fruget 2017). Their population genetic histories reveal different trajectories between the species and across their temporal and spatial courses, which may significantly influence their relative successes and interactions.
Calazans CSH, Walters LJ, Fernandes FC, Ferreira CE, Hoffman EA (2017) Genetic structure provides insights into the geographic origins and temporal change in the invasive charru mussel The following supplementary material is available for this article: Table S1. Comparative genetic variation of (A) zebra mussel and (B) quagga mussel populations based on 15 and 10 nDNA microsatellite loci, respectively. Table S2. Loci analyzed for (A) Zebra mussels and (B) Quagga mussels. Table S3. Colony analysis results for (A) zebra mussel and (B) quagga mussel populations. Table S4. Assignment test results, using Geneclass2, for (A) zebra mussels and (B) quagga mussels. Appendix 1. References to Table S2. Figure S1. Multivariate Discriminant Analysis of Principal Components (DAPC) for zebra mussel populations from Lake Erie. Figure S2. Multivariate Discriminant Analysis of Principal Components (DAPC) for quagga mussel populations from Lake Erie.