To characterise abundance, biomass and demographic structure, samples of zebra and quagga mussels were collected at eight locations on the Moselle River (Fig. 1). The first sampling event was carried out between 4^th and 21^st May 2021 at all eight locations. Thereafter, the sampling events were limited to three stations in the downstream stretch of the Moselle (Thionville, Cattenom-upstream and Sierck-les-Bains; Fig. 1), with four sampling events taking place over a period of 12 months, i.e. 4^th May 2021, 5^th July 2021, 22^nd November 2021 and 5^th May 2022.

On each sampling date at each site, mussels were carefully collected from all submerged substrata (e.g. rocks and stones) close to the shoreline at ca. 1m depth until a minimum of 300 mussels had been collected. All mussel samples were then preserved in 96% ethanol for further analysis in the laboratory, where they were morphologically identified to species using key identification characteristics such as flatness of the ventral face, shape of the carina, shape of the ventral junction and position of the byssus (May and Marsden 1992; bij de Vaate and Jansen 2007; Sablon et al. 2010; Teubner et al. 2016). All individuals were then measured using an electronic calliper (± 0.01 mm). and separated into size classes of 1mm for further analysis.

The number of individual zebra and quagga mussels was evaluated for each sample, whereupon the biomass and biovolume of quagga mussels in each sample was calculated based on regression models between length and biomass and length and biovolume.

To assess the relationship between length and biomass, ca. 100 mussels of each species were sampled in the Moselle River at Sierck-les-Bains on 23^rd January 2022 and kept alive in a plastic container with an air bubbler for the duration of the experiment. To feed the mussels, around 1/3 of the water was changed every two to three days, using water from the Rhine River. In the laboratory, each mussel was placed in hot water for a few minutes to open the valves. Once opened, the mussels were identified to species, their shells cleaned with a scalpel and their length measured with an electronic calliper. The shell and soft bodies were then stamped in an absorbent cloth, following which the wet weight of soft tissue (WW_{soft tissue}) and wet weight of the shell (WW_shell) were immediately measured using a pre-weighed aluminium cup and a microbalance (± 0.0001 g). After 24h at 105 °C in a stream room, the aluminium cups were re-weighed on the same microbalance to obtain the dry weight of soft tissue (DW_{soft tissue}) and dry weight of the shell (DW_shell). Finally, the aluminium cups containing the dried soft tissue were placed in an oven at 450 °C for 2h to obtain the ash-free-dry-weight (AFDW) of the tissue, with AFDW representing the organism’s biomass. To prevent rewetting of the dried material, the aluminium cups were stored in a desiccator before weighing.

To assess the relationship between mussel length and biovolume, 100 zebra and quagga mussels were selected from the sample taken at Sierck-les-Bains on 22^nd November 2021 and the shell length of each mussel measured using an electronic calliper. The biovolume of each individual was then calculated by immersing each mussel into a 25ml graduated cylinder half-filled with water, with the change in volume after immersion representing mussel biovolume.

To determine the number of cohorts present, the size-frequency-distributions for zebra and quagga mussels were analysed separately for each combination of site (Thionville, Cattenom-upstream, Sierck-les-bains) and sampling period (4^th May 2021, 5^th July 2021, 22^nd November 2021, and 5^th May 2022), using the graphical method of Bhattacharya (Bhattacharya 1967), giving a total of 12 size-frequency-distributions for each species. All analyses were conducted using a Microsoft Excel macro developed by one of the co-authors (Jean-Nicolas BEISEL). Despite our efforts to sample the widest range of sizes, some cohorts obtained from a site were based on < 10 individuals. Thus, decided to combine the observations of Thionville, Cattenom-upstream and Sierck-les-Bains into a single dataset for each sampling event and each species to ensure we had a statistically relevant number of individuals and mitigate potential sampling bias. While it is possible there may be ecological differences between stations that could influence growth rate or timing of reproduction, the three stations are spatially close and are all located in the navigable stretch of the Moselle, i.e. they have the same channel width and substrate. Furthermore, the cohorts identified for these three stations separately appeared to be very similar in length (i.e. ± 1.5 mm).

Once all cohorts had been identified for each sampling period, we used a growth-at-length model specific to the Moselle (Beisel et al. 2010) to match zebra mussels cohorts across sampling events. Hypothetical growth rates were then calculated by dividing the difference in length between two sampling periods by the number of days between them.

As there is no known growth-at-length model for the quagga mussel as yet, it was not possible to calculate the age of each cohort. Consequently, for each sampling event, we considered that the growth rate of a cohort was negatively correlated with its age. Using the little data available in European ecosystems (i.e. D’Hont et al. 2018; Haringvliet River, the Netherlands), we calculated an assumed growth rate by dividing the length difference of a cohort between two sampling periods in mm.d^-1. For both species, cohorts with < 20 individuals were not considered significant enough to be included in the analysis.

Genomic DNA was extracted from a subsample of 20 randomly picked quagga mussels from each site to evaluate genetic diversity and structure. DNA extraction was performed using the DNeasy Blood and Tissue kit (QIAGEN, Netherlands), following the manufacturer’s instructions, and stored at -20 °C for further analysis.

Each mussel was then genotyped for the 700bp-long cytochrome c oxidase subunit I mitochondrial region (COI) using the universal primers LCO-1490 and HCO-2198 (Folmer et al. 1994). Amplifications were performed in 25 μL reaction mixtures containing 1X GoTaq reaction buffer (1.5 mM MgCl₂), 0.2 mM of each dNTP, 0.5 μM of each primer, 0.5 U of GoTaq DNA Polymerase (Promega) and 2 μL (ca. 10 ng) of genomic DNA. The samples were then subjected to a polymerase chain reaction (PCR), with conditions comprising an initial denaturation at 94 °C for 4 min, followed by 30 cycles of 45 s denaturation at 94 °C, 45 s annealing at 40 °C, and 50 s elongation at 72 °C, with a final elongation step of 10 min at 72 °C. PCR products were Sanger-sequenced with the same primers used for amplification (Genewiz, Leipzig, Germany). Chromatograms were visualised and edited in Geneious R9.1.8.

In addition, each mussel was genotyped using the 10 microsatellite (µsat) loci developed by Wilson et al. (1999) (Dbug) and Feldheim et al. (2011) (Dbu). PCR amplifications were performed as in Marescaux et al. (2016b), with 1 µl of PCR product being mixed with 10 µl of Hi-Di formamide and 0.1 µl of GeneScan-500 (LIZ) (Applied Biosystems) for genotyping using a ABI3730 48 capillary DNA Analyser (IPG, Gosselie, Belgium). Peak Scanner v2.0 (Applied Biosystems) was then used to score alleles and genotype individuals.

To illustrate mitochondrial diversity, the COI dataset was aligned with previously published sequences (GenBank accession number DQ840132, DQ840133 – EF080861, EF080862, EU484436, JN133734–JN133747, JQ756297, JQ756298, KJ881409–KJ881415, KP057252, MF469063–MF469065, MK358469, MK358470, U47650) in MAFFT (E-INS-i method; Katoh and Standley 2013) and a Median-Joining Network (Bandelt et al. 1999) calculated (HaplowebMaker; Spöri and Flot 2020).

The Moselle microsatellite dataset was checked for the presence of null alleles using MICROCHECKER v.2.2.3 (Van Oosterhout et al. 2004). For each combination of locus and site, linkage disequilibrium (LD) and Hardy–Weinberg equilibrium (HWE) was tested with GENEPOP (Markov chain parameters set to 1000 for dememorisations, batches and iterations per batch; Raymond and Rousset 1995). Levels of significance for all multiple tests were adjusted using Bonferroni correction (Rice 1989). Molecular variance (AMOVA) was also analysed to assess the degree of genetic diversity within and among populations using Arlequin. In addition, Bayesian clustering of microsatellite genotypes was performed using STRUCTURE v.2.1 (Pritchard et al. 2000), with the optimal number of clusters (K) inferred as in Marescaux et al. (2016b).

We found that zebra mussel largely dominated the dreissenid populations of the Moselle River, with quagga mussels representing between 28 and 51% in our samples, with no clear upstream/downstream pattern (Fig. 2). On the other hand, the percentage of quagga mussels varied between sampling events, especially at the Thionville site, where we harvested 48, 10, 49 and 40%, respectively, over four successive sampling events. Despite this, the variation in abundance had no effect on the dominance of zebra mussels for all sites and sampling events, with percentages at Cattenom-upstream being 36, 29, 31 and 45%, respectively, and 32, 28, 41 and 45%, respectively, at Sierck-les-Bains, during successive sampling events.

While we recorded significant positive relationships between length and biomass, and between length and biovolume, for both zebra and quagga mussels, there was no significant differences between species (ANCOVA tests; p > 0.05; Table 1). Based on these relationships, we were able to convert the data for individual size frequencies into biomass and biovolume. The results for biomass differed greatly from those for density, with quagga mussels having around twice the biomass of zebra mussels, ranging between 29 and 81% of total biomass (Fig. 2). Maximum biomass was achieved at the La Maxe site, where quagga mussels represented 81% of total dreissenid biomass during the first sampling event. At Thionville, quagga mussels, represented just 29% of dreissenid biomass and 10% by population abundance during the second sampling event. The quagga biovolume results were similar to those for biomass, with values ranging from 24 to 75% of total biovolume (Fig. 2), indicating that quagga mussels often represent more than twice the biovolume of zebra mussels, despite them being the dominant species at a site.

Analysis of the different cohorts on each sampling event provided us with a comparative framework for population structure. The size structure of the two species differed greatly (Fig. 3). For example, young individuals with a shell length < 15 mm were dominant in the zebra mussel population, representing ca. 70% of the population, with 33% being made up of individuals with a shell length < 10–11 mm (i.e. < 1 yr). In contrast, the quagga mussel population was dominated by larger individuals, with just 22% of the population being < 15 mm in length, with those measuring < 10–11 mm representing just 7%. Unlike zebra mussels, where a large proportion of small individuals survived after attachment, quagga mussels of < 15 mm length appear to have undergone large-scale mortality.

Size-frequency analysis for zebra mussels using Bhattacharya’s method revealed seven cohorts between the first (May 2021) and final (May 2022) sampling events (Table 2, cohorts A, B, C, D E, F, G). The theoretical lengths of these cohorts, calculated using a previous growth-at-length model for the Moselle (Beisel et al. 2010), were very close to the observed cohort lengths (< 1 mm difference, yellow square, Table 2). Thus, our model provided theoretical results that were consistent with the observed growth of a cohort between different sampling events, supporting the hypothesis that the zebra mussel growth rate was indeed maximal in spring, with a maximum of 0.081 mm/day between May and June, and almost nil in winter, at between 0.007 and 0.002 mm/day between November and March.

For the quagga mussel, we identified five cohorts for the first and second sampling events, and three for the third and fourth sampling events (Fig. 4). While very few individuals of < 15mm were found during the third and the fourth sampling events (22^nd November 2021 and 5^th May 2022), the data suggests the presence of two additional cohorts for the third sampling events at ca. 3 and 11 mm); however, these were not included in the analysis due to the low number of individuals. Assuming that (i) the growth rate of a cohort is proportional to its age (Stoeckmann 2003), (ii) the growth rate of the quagga mussel is greater than that of the zebra mussel at the same age (Neumann et al. 1993; Stoeckmann 2003; D’hont et al. 2018), and (iii) the quagga mussel grows throughout the year, even in winter (D’hont et al. 2018), we linked the cohorts observed during the different sampling events. This suggests that five different cohorts occured between the first sampling event in May 2021 and the fourth in May 2022 (i.e. A, B, C, D and E in Fig. 4). During the fourth sampling event, two cohorts (i.e. cohorts F and G, Fig. 4) could not be linked to those identified during the previous sampling events, while a third cohort almost certainly links with cohort B as its growth rate was close to that recorded in D’Hont et al. (2018) for the winter period (i.e. 0.02 mm/day). Quagga mussels achieved their highest growth rate between the first and second sampling campaigns, i.e. between May and July, with a maximum rate of 0.104 mm/day in cohort B (Fig. 4).

To effectively compare the growth of quagga and zebra mussels in 2021, we chose to analyse the small quagga mussels represented by cohorts A, B and C (3.9, 7.3 and 15.5 mm, respectively) sampled during the first sampling event. We then used the growth model specific to zebra mussels from the Moselle to calculate the projected length of the three cohorts during the three other sampling events. This model has proved effective in predicting zebra mussel length in previous analyses and thus provides a useful comparison with the length of quagga mussels in cohorts A, B and C observed during the first three sampling events.

Between the first and second sampling event (May to July 2021), cohort A had the smallest difference in length between the two species, with zebra mussels being 8.9 mm and quagga mussels 9.7mm (Fig. 5). The associated growth rate for this period was 0.081 mm/day for zebra mussels and 0.093 mm/day for quagga mussels. For cohorts B and C, however, quagga mussels were much larger (in average about 2 mm) than the zebra mussels, with cohort B showing the greatest difference in growth rate between the two species at 0.10 mm/day (mean) for quagga mussels and 0.071 mm/day for zebra mussels (Table 3). Between the second and third sampling events (July to November 2021), the mean length zebra mussels in cohorts A, B and C stabilised rapidly, reaching 13.7, 15.9 and 21.2 mm, respectively, during the third sampling event (Fig. 5). In contrast, quagga mussel cohorts A, B and C continued to grow to a length of 18.8, 21.6 and 25.8 mm, respectively. Overall, the zebra mussel cohorts displayed a slower growth rate than quagga mussels, with cohort A having growth rates of 0.034 and 0.065 mm/day, cohort B 0.030 and 0.056 mm/day, and cohort C 0.020 and 0.040 mm/day for zebra mussels and quagga mussels, respectively (Table 3). In autumn, quagga mussel growth rates were almost twice as high as those for zebra mussels, with zebra mussel growth rates in cohorts A, B and C being very low (0.004–0.007 mm/day) between the third and fourth sampling event (November 2021 to May 2022), while the maximum length of the only known quagga mussel cohort (B) during the final survey (May 2022) reached 25.8 mm, with a mean growth rate of 0.025 mm/day.

Four distinct haplotypes were identified among the 60 quagga mussels genotyped for the Moselle River (Fig. 6, Q1, Q2, Q4 and Q5), all of which have been retrieved in previous studies (see Marescaux et al. 2016b), including Q1, the most widespread haplotype historically, which was dominant at all sites. At each of the three sites sampled, the Q1 haplotype represented between 70 and 95% of all individuals sampled, while the other three haplotypes corresponded to variants found at lower frequencies in other European rivers in 2016 (haplotype Q2 in the Ijsselmeer, Markemeer, Meuse, Rhine, Main and Danube rivers; haplotype Q5 in the Danube River and Datteln-Hamm canal; haplotype Q4 has also been recorded in North America) (Fig. 6). Interestingly, haplotypes Q1 and Q4 were previously detected in the Moselle at similar frequencies between 2011 and 2013 (Marescaux et al. 2016b).

STRUCTURE analysis performed on the entire microsatellite dataset, which included all sequences observed in this study alongside those of Marescaux et al. (2016b), failed to distinguish genetic clusters (data not shown). This was corroborated by the AMOVA analysis, which indicated that just 0.59% of genetic variation found in the quagga mussel population was retrieved when comparing worldwide groups (Table 4).

Ten years after its first observation in the Moselle River, our results show that the quagga mussel represents between 28 and 51% of all dreissenids sampled, compared with just 0–1.5% (i.e. 0–3 individuals per site) at the same locations in 2011 (bij de Vaate and Beisel 2011). Though the percentages observed between 2021 and 2022 varied between sites, there was no clear upstream­-downstream distribution pattern along the river, suggesting that the variation determinant may be linked to variations in flow rate through the year.

While our findings indicate a 30-fold increase in the quagga mussel population in ten years, such an increase is not always observed in running waters. In the USA for example, quagga mussels made up only 1% of dreissenids after 12 years of coexistence on the Mississippi and Ohio rivers (Grigorovich et al. 2008). In the Albert Canal (Meuse River, Belgium), quagga mussels accounted for 80% of all dreissenids in 2014, just three years after its first observation (Marescaux et al. 2015). Heiler et al. (2013) also showed a dominance of quagga mussels in rivers and canals in Germany, including the Main-Danube Canal (76–100% at some sites), the Main River (> 80%) and the Upper Rhine (80–90%). In Belgium (Meuse River; Marescaux et al. 2015) and Germany (Heiler et al. 2013), quagga mussel dominance was achieved in less than three years following its first observation/introduction (date of introduction modelled in Heiler et al. 2013). In the Moselle, however, quagga and zebra mussels have been coexisting for at least 10 years (bij de Vaate and Beisel 2011) with quagga mussels still not dominant. While Strayer et al. (2019) noted a wide range of temporal dynamics in their analysis of long-term (>10 years) dreissenid population datasets from 47 lakes and three rivers in Europe and North America, a general pattern that applied across many populations was that quagga mussels arrived later than zebra mussels and usually caused large declines in the zebra mussel population over a relatively short period. After 10 years, they felt it unlikely that zebra mussels would continue to decline any further due to the presence of quagga mussels.

Our results showed that, locally, quagga mussels could represent more than twice the biomass or biovolume of zebra mussels in the Moselle, with biomass ranging between 29 and 81% and biovolume between 24 and 75% of the dreissenid population.

As filtration rate is dependent on the size/biomass of individuals, the ecological impacts of invasive species such as dreissenid mussels tends to be dependent on biomass rather than density (Marescaux et al. 2016a; Strayer et al. 2019). Even where quagga mussels are less numerous than zebra mussels, therefore, they can still cause a greater ecological impact. At present, assuming equal sampling effort throughout the year, we are unable to explain the 10% of quagga mussels found at Thionville during the second sampling event. The sites of Thionville, Cattenom-upstream, and Sierck-les-Bains are all located in the same navigated sector and very similar.

The relationship between biomass and length was more variable for quagga mussels (Table 1; R² = 0.72) than zebra mussels (R²= 0.85). Burlakova et al. (2006) has shown that the biomass of a bivalve can vary greatly over the course of a year, especially during the reproduction period when the gonads make up a greater proportion of overall weight. While some studies have shown that quagga mussels are able to reproduce continuously throughout the year (Marescaux et al. 2015; Hesselschwerdt and Teiber-Siesseger 2021), which may account for this difference in biomass between the species, we did not specifically investigate whether quagga mussels could reproduce during the winter months in this study.

In this study, at a given date, the zebra mussel population was composed of five different cohorts, resulting in demographic dynamics of seven cohorts over one year (Table 2). Overall growth rate was highest in spring, reaching 0.081 mm/day between May and July, and very low in winter, at just 0.007–0.002 mm/day between November and March. These results are in accordance with those of D’Hont et al. (2018) for the Haringvliet River (The Netherlands), where highest growth rates occurred during the summer (0.08 mm/day) and lowest during winter (0.01 mm/day). A previous study in the Moselle also recorded a maximum growth rate of between 0.087 and 0.135 mm/day between June and August, though only for mussels < 3 months old (Beisel et al. 2010). In the Lower Rhine (Germany), a similar study found a maximum growth rate of 0.095 mm/day in June, though also for individuals < 3 months of age (Neumann et al. 1993; Jantz and Neumann 1998). Finally, in Lake Maarsseven (The Netherlands), Dorgelo and Smeenk (1988) reported that zebra mussels of 5.2 mm length grew at an average rate of 0.046 mm/day between June and November. In our study, we recorded zebra mussels of 6 mm length at the beginning of July having a mean growth rate of 0.03 mm/day between July and November (Table 2). Compared to previous studies, therefore, the growth rate of zebra mussels in the Moselle is similar to that reported in other European ecosystems, confirming that optimal growth clearly takes place during spring.

As with Zebra mussels, we recorded five separate quagga mussel cohorts at any given time, with seven cohorts observed over the course of the year (Fig. 3). The quagga mussel growth rate was also highest in spring, at 0.104 mm/day between May and July; however, unlike zebra mussels, the growth rate remained relatively high in winter at 0.025 mm/day. While there have been few exhaustive studies on the growth of quagga mussels in running waters in Europe, there have been several in the Laurentian Great Lakes of the USA (Malkin et al. 2012; Casper et al. 2014; Karatayev et al. 2018; Elgin et al. 2021). For example, a study done by Karatayev et al. (2018) in the Lake Erie have shown that a 12 mm quagga mussel have a mean growth rate of 0.018 mm/day at 15m water depth, between May and December (mean value recalculated by us). Another study by Elgin et al. (2021) in Lake Ontario showed that a 12 mm quagga mussel gained 10.2 mm in one year at 15m depth. However, comparing our values to those obtained in Laurentian Great Lakes is not an easy task as these environments are far too different from the Moselle, especially regarding their thermal regime. D’Hont et al. (2018), however, in their study on growth of quagga mussels in the Haringvliet River, found a maximum growth rate of 0.09 mm/day in summer, and an average growth rate in winter of 0.02 mm/day (recalculated by us), which are close to our own values. It should also be noted that the cohorts identified were often overlapping, suggesting that quagga mussels may be reproducing continuously throughout the year.

Our comparison between zebra and quagga mussels (Fig. 4, Table 3) showed that quagga mussels always displayed a higher growth rate than zebra mussels at any time of the year, with the quagga mussel growth rate being 1.4× greater between May and July, 1.5× greater between July and November and 4× greater in winter (zebra mussel = 0.006 mm/day, quagga mussel = 0.025 mm/day, cohort B, Fig. 4). Note, however, as we were only able to follow one cohort for quagga mussels through winter, and as there is only one reference in the literature to compare with (D’Hont et al. 2018), we recommend that further studies are needed to better define quagga mussel growth over the winter period. As a comparison, Casper et al (2014) observed in the St-Lawrence River that quagga mussel had a growth rate almost twice that zebra mussel, between July and September with water temperature between 19 and 27 °C.

Our analysis also showed a significant difference in the size structure of the two species, with very few quagga mussels < 15 mm found, especially during the third and fourth sampling events. A similar lack of smaller size classes was also observed in the Meuse River by Marescaux et al. (2015). One potential explanation is predation by the round goby, Neogobius melanostomus, which is known to occur at high abundances in both the Moselle and Meuse (Verreycken et al. 2011; Manné et al. 2013). Round gobies in the USA are known to feed on dreissenids (Houghton and Janssen 2013), with up to 90% of its food potentially consisting of mussels (Andraso et al. 2011), particularly dreissenids of between 5 and 13–14 mm (Djuricich and Jansen 2001; Andraso et al. 2011; Naddafi and Rudstam 2014), though the zebra mussel’s position on or within the substrate and their shell shape can make it more difficult to predate (Houghton and Jansen 2014). Quagga mussels are also believed to have greater attachment strength than zebra mussels (Peyer et al. 2009), which could also result in differential predation between the two species, and the very particular quagga mussel population structure observed in this study.

Our genetic results appear to confirm that most individuals sampled worldwide correspond to a single haplotype (Q1) originating from the species’ native area in the Pontic region (Marescaux et al. 2016b). This well-established haplotype suggests a first, or a single massive, dispersion event, followed by successive additional invasions by other haplotypes thereafter. Interestingly, our data highlighted two recent dispersions of haplotypes Q2 and Q5, these being retrieved during one sampling event but not before (see Marescaux et al. 2016b). Though haplotype Q2 has previously been detected on the Meuse River, and could have been transported downstream and then upstream again by boat traffic, the Q5 haplotype is mainly restricted to the Danube River and the Datteln-Hamm canal, both of which are directly connected to the Moselle via the Main-Danube canal and thereafter the Main and Rhine Rivers. These results could potentially indicate either new invasion corridors, considerable changes in haplotype distribution since the last comprehensive study or an undervaluation by Marescaux et al. (2016b) due to under-sampling. Given the results obtained from the microsatellite data, one might expect that genetic exchanges are frequent between distinct populations in rivers/canals, and even between groups from different continents. This observation appears to support our hypothesis that new invasion events have occurred repeatedly since the last comprehensive study and that more haplotypes might also be observed in populations from other rivers.