Global invasion genetics of two parasitic copepods infecting marine bivalves

Invasive species, and especially invasive parasites, represent excellent models to study ecological and evolutionary mechanisms in the wild. To understand these processes, it is crucial to obtain more knowledge on the native range, invasion routes and invasion history of invasive parasites. We investigated the consecutive invasions of two parasitic copepods (Mytilicola intestinalis and Mytilicola orientalis) by combining an extensive literature survey covering the reported putative native regions and the present-day invaded regions with a global phylogeography of both species. The population genetic analyses based on partial COI sequences revealed significant population differentiation for M. orientalis within the native region in Japan, while introduced populations in North America and Europe could not be distinguished from the native ones. Thus, M. orientalis’ invasion history resembles the genetic structure and recent spread of its principal host, the Pacific oyster, Crassostrea gigas, while M. intestinalis lacks population genetic structure and has an overall low genetic diversity. Therefore, the native origin of M. intestinalis remains unclear. With this study, we demonstrate that even highly related and biologically similar invasive species can differ in their invasion genetics. From this, we conclude that extrapolating invasion genetics dynamics from related invasive taxa may not always be possible.

Cross-infections of hosts and parasites from separate invasion fronts where the parasite had invaded both locations at a similar time (i.e., similar lengths of coevolutionary time) showed that parasite invasion led to different evolutionary trajectories, with differences in the parasite's infectivity and the host's resistance and tolerance 13 . These trajectories were characterized by specific molecular interactions between host and parasite 14 . Mytilicola orientalis was co-introduced with the invasive Pacific oyster, Crassostrea gigas, through aquaculture transports 15,16 , and now also infects a suite of native hosts, including blue mussels, M. edulis [16][17][18] . Thus, the native host, M. edulis, is now shared between M. intestinalis and M. orientalis.
The reported invasion histories already suggest substantial differences in the mode of invasion between these two Mytilicola species that could, in turn, lead to dissimilarities in ecological effects and evolutionary responses in M. edulis. However, while the major invasion route of the principal host of M. orientalis, the Pacific oyster, is known 19,20 , the invasion route of its hitchhiking parasite has yet to be validated. Similarly, the dispersal routes of Mytilus galloprovincialis are better known than those of its parasitic copepod M. intestinalis 21,22 . To fully understand the parasites' invasion histories and their consequences for ecological and evolutionary processes, the genetic structure in the native ranges and along invasion routes of both parasite species needs to be investigated.
Mytilicola is a genus of cyclopoid copepods that have a direct life cycle involving one pelecypod host. Their free-living larval stages are short and have mainly been studied in M. intestinalis (e.g. 10,23,24 ). The duration of the free-living and parasitic larval stages of the other Mytilicola species are thought to be similar 25 . The two free-living larval stages, nauplius and metanauplius, last 36 to 48 hours at a laboratory temperature of 18 °C 23 . The next larval stage, the first copepodite, is infective and will search for a suitable host to infect. It may live for 11 to 15 days under laboratory conditions, but it is doubtful whether it lives that long in the field 23 . Once Mytilicola has infected a host, its dispersive capacity depends on its bivalve host, which is sessile and does not move across long distances unaided. Therefore, the dispersal capacity of Mytilicola was considered to be limited 12,26 .
In this study, we aimed to (1) use the available literature to reconstruct the invasion history of both parasites, and (2) test whether the reported invasion routes could be confirmed by a global phylogeography based on partial cytochrome-c-oxidase I (COI) sequences. We specifically wanted to verify their putative areas of origin (i.e. Mediterranean Sea Basins for M. intestinalis and Japanese seas for M. orientalis), and determine the genetic diversity and population differentiation between putative native and invaded regions for both invasive Mytilicola species.  Supplementary Information). The literature included peer-reviewed publications, fishery leaflets, reports, book sections and dissertations. A total of 41 publications with potentially relevant information could not be read because the full texts were not available online, neither accessible via Dutch and German national libraries, nor by attempts to contact the original authors, and thus, could not be included.

Literature-based invasion history.
The distribution of Mytilicola intestinalis before the 1930s was limited to the Mediterranean Sea Basin, where the parasite was recorded from a few locations near large ports and aquaculture areas (Figs 1A, S1 and Online Resource 1). The hypothesized native region of M. intestinalis is, therefore, the Mediterranean Sea Basin. In the 1930s, its distribution extended to a few specific locations at the German and English coasts, from which, in the following years, the parasite spread to coasts of neighboring countries. In the late 1940s, the first sightings in the Netherlands and Ireland were reported, followed by the Atlantic coasts of France in the early 1950s. In 1951, researchers collaborated to record the invasion and potential effects of M. intestinalis along the Atlantic coasts of Europe, which resulted in the publication of a special issue in Revue des Travaux de l'Office Scientifique et Technique des Pêches Maritimes on M. intestinalis 27 . Ship hull fouling was seen as the main vector of M. intestinalis introductions, after which further spread was aided by the relaying of mussels (i.e., from the German Wadden Sea to aquaculture areas in Zeeland, the Netherlands). It is, thus, unclear from which Mediterranean basin (or basins) M. intestinalis was introduced, and whether this happened repeatedly, which might lead to different expectations regarding signatures of bottlenecks and genetic diversity in the introduced region.
The distribution of M. orientalis in its putative native range has not been the subject of many studies, which is reflected in our literature survey (Figs 1B, S2 and Resource 2). The co-introduction of M. orientalis with the Pacific oyster in North America was better recorded. Where young Pacific oysters or infected stock were imported, M. orientalis was also found, but the parasite did not spread quickly to adjacent natural areas 26 . In the 1970s, Pacific oysters from North America were imported into France, and, once more, M. orientalis was co-introduced across half the globe. In 1993, the species was also found in Ireland in imported oysters from France 28 . There are many gaps in the distribution of M. orientalis in Europe, which presumably reflect the lack of studies in these areas ( Supplementary Fig. S2, Online Resource 2). new records and infection levels of Mytilicola spp. A total of 5604 Pacific oysters (Crassostrea gigas), 2658 Mediterranean mussels (Mytilus galloprovincialis), more than 1595 blue mussels (Mytilus edulis), 104 bay mussels (Mytilus trossulus) and 90 mussels (Brachidontes pharaonis) were dissected in this study. Dissections of Pacific oysters and mussels from 20 locations along Mediterranean coasts, 13  Mytilicola intestinalis was detected in M. galloprovincialis and M. edulis. Mytilicola orientalis was also detected in these bivalves, but also in C. gigas and M. trossulus. Prevalences of M. intestinalis ranged from 0 to 82.0% in the Mediterranean Sea Basins, the putative native region of M. intestinalis, and from 0 to 93.8% in the invaded regions ( Table 1). Prevalences of M. orientalis ranged from 0 to 10.6% in Japan, the putative native region of M. orientalis, and from 0 to 73.8% in the invaded regions (Table 1). invasion genetics of Mytilicola spp. A total of 424 M. intestinalis and 444 M. orientalis were successfully sequenced, which resulted in the detection of 18 and 26 haplotypes, respectively (Tables S1 and S2) (Genbank accession No. MN334483-MN334526). The minimum spanning network of M. intestinalis had a star-like structure with one central, highly abundant haplotype from which several haplotypes were separated by a few mutational steps ( Fig. 2A). The star-like structure of the network suggests recent demographic expansion, as is also suggested by the negative Tajima's D and Fu's F S values which were significant for most locations (Table 2), and by the mismatch distributions (Fig. 3). Nine of the ten non-singleton haplotypes were shared between the putative native and invaded range (Figs 2A, 4 and Supplementary Table S1); only haplotype Mi16 was solely found in the invaded range (La Rochelle, Supplementary Table S1).
In contrast, the minimum spanning network of M. orientalis was relatively complex with two common frequent haplotypes that both showed associated derived haplotypes (Fig. 2B). With the exception of one location in Japan, all locations had one of the main haplotypes as its most frequent haplotype. A small group of related singleton haplotypes was found only in North America (light purple colored haplotypes Mo11-Mo13, Figs 2B, 3 and Supplementary Table S2). This type of network is more typically associated with stable demographics, also suggested by the non-significant Tajima's D and Fu's F S values (Table 2), and by the mismatch distributions (  (Table 3).
For M. orientalis, we adopted the default expectation of high genetic diversity in the whole native region and genetic differentiation among populations in the absence of detailed knowledge of the distribution in the native range. Because of the many transport occurrences of Pacific oysters to North America and Europe, we did not expect to find a reduction in M. orientalis genetic diversity in the invaded regions. The overall level of population differentiation for M. orientalis was higher and estimated at Φ ST = 0.05242 (p < 0.00001) and F ST = 0.03002 (p < 0.00001), with an overall population differentiation in Japan of Φ ST = 0.10287 (p < 0.00001) and F ST = 0.07019 (p < 0.00001), and in the introduced region of Φ ST = 0.02243 (p = 0.0303) and F ST = 0.01005 (n.s.) ( Table 3).
The small, overall genetic differentiation in M. intestinalis could not be attributed to any single location or group of locations based on the pairwise population comparisons because none of the comparisons were significant after Bonferroni correction (Supplementary Tables S3, S4). The hierarchical AMOVAs did not clarify this further, as they did not detect significant groupings on any of the hierarchical levels. Including Cassis into the Western Mediterranean group led to negative Φ CT values (Table 4). When grouping Cassis into the Ligurian Sea, the percentage of variation among groups was 0.38%, among populations within groups 1.24%, and within populations 98.38%, resulting in Φ CT = 0.0038 (n.s.), Φ SC = 0.01247 (n.s.) and Φ ST = 0.0162 (p = 0.0480) ( Table 4). Using classical F-statistics, the percentage of variation among groups was greater at 2.07%, smaller among populations within groups at 0.84%, and similar within populations at 97.12%, resulting in F CT = 0.02069 (n.s.), F SC = 0.00825 (n.s.) and F ST = 0.02878 (n.s.) ( Table 4). In both hierarchical AMOVAs using mitochondrial statistic Φ and conventional F-statistic, there was no significant difference between the Gulf of Lion, the Ligurian Sea and the Adriatic Sea. Therefore, the small amount of population structure that was detected in the overall test cannot clearly be attributed to one of the Mediterranean Sea Basins, nor to one or a few individual locations.
In contrast to M. intestinalis, there was strong and significant genetic differentiation between locations in the native region of M. orientalis. Iwaya, located in the Sea of Japan, was strongly differentiated from nearly all other sampling locations along the Seto Inland Sea (Φ ST = 0.215-0.290, Supplementary We did not detect significant differences in haplotype diversity (h) between the putative native regions and the introduced regions (ANOVA, F 1, 26 = 1.1081, p = 0.302) (Fig. 5A). Haplotype diversity for M. intestinalis ranged from 0.131 to 0.569 in the putative native region and from 0.260 to 0.608 in the invaded region (Supplementary Table S1), while for M. orientalis, h ranged from 0.627 to 0.820 in its putative native region,   www.nature.com/scientificreports www.nature.com/scientificreports/ Japan, (all locations with >15 individuals) and from 0.699 to 0.839 in its introduced region (Supplementary  Table S2). Similarly, no significant differences were detected for nucleotide diversity (π) between the putative native region and the invaded regions (ANOVA, F 1,26 = 0.0284, p = 0.867). Nucleotide diversity for M. intestinalis ranged from 0.000276 to 0.00179 in the putative native region, and from 0.000572 to 0.00267 in the introduced region (Fig. 5B, Supplementary Table S1), and M. orientalis' nucleotide diversity ranged from 0.00423 to 0.00550 in Japan and from 0.00366 to 0.00541 in the introduced regions (Fig. 5B, Supplementary Table S2). However, both haplotype diversity and nucleotide diversity were significantly larger for M. orientalis than for M. intestinalis (ANOVA h: F 1, 26 = 99.2662, p < 0.00001; ANOVA π: F 1, 26 = 251.7502, p < 0.00001). The interaction term in the linear models (species × region) was not significant for haplotype diversity h or nucleotide diversity π (ANOVA h: F 1, 26 = 0.5634, p = 0.460; ANOVA π: F 1, 26 = 1.3178, p = 0.261), indicating that genetic diversity differed consistently between species.

Discussion
Our study utilized the consecutive invasions of two congeneric invasive parasites to align the invasion history derived from their reported spread in the literature with their phylogeography based on molecular markers. Based on mitochondrial haplotype diversity, we found that the native range and the invasion route of Mytilicola orientalis match the classical pattern of a recent invasion described in the literature, while the phylogeography of Mytilicola intestinalis did not match our literature-based expectations due to low mtDNA genetic diversity. While several scenarios could explain the low genetic diversity, we cannot identify the origin of the species, which potentially also makes it a cryptogenic species in the Mediterranean Sea Basins. www.nature.com/scientificreports www.nature.com/scientificreports/ M. orientalis phylogeography supports invasion history. Our study confirmed the coastal seas of Japan as the native region of M. orientalis, because of the significant population genetic structure between the Seto Inland Sea and Iwaya in the Sea of Japan, and the stronger overall population differentiation in Japan than in the invaded regions. This was expected based on the invasion history derived from our literature survey. The invasion of M. orientalis was extremely long-distance, crossing both the Pacific and Atlantic Ocean via aquaculture transports, followed by a relatively rapid spread through the invaded regions (Fig. 1B, Resource 2) that is still ongoing 34 (Feis, unpubl. data).
Pacific oysters have been imported from Japan to locations across the globe, with Miyagi Prefecture as the main overseas export area and Kumamota and Hiroshima Prefectures including the Seto Inland Sea 35 contributing smaller amounts 35 . Oyster seed is regularly transported between the regions largely homogenizing the genetic diversity and restricted gene flow has only been reported for the Nagasaki and Kumano regions 36 . Mytilicola orientalis hitchhiked with oyster transports to North America 16 and to Europe 15,37,38 , and our population genetic analyses support these patterns in the native range and invaded regions: the North American and European samples could not be distinguished from the Seto Inland Sea samples from Japan, while the sample from Iwaya in the Sea of Japan was significantly different.  Fig. 4. Figure 2 was drawn in Adobe Illustrator based on the Arlequin output list of OTU differences.
Within the invaded region in Europe, M. orientalis showed more resemblance in the MDS plots between North America and the northern Wadden Sea invasion: Belfair (representing North America) and Puan Klent (representing the northern Wadden Sea oyster invasion). This pattern represents a small mismatch with the phylogeography of the host 19,20 and might suggest that other vectors than oyster transports contributed to M. orientalis' secondary spread in the Wadden Sea. For example, the wider host range of M. orientalis and spillover to native bivalves 18 could have further sped up and modified its spread. In the Mediterranean, the distribution of M. orientalis is disjunctive, with presence at the Thau lagoon in France since 1979 39 and in the northern Adriatic 40 , but an absence in between (Table 1). This patchy distribution pattern probably reflects the long-distance transports of oysters and M. orientalis' ongoing invasion, and furthermore demonstrates that Mytilicola spp. do not colonize rapidly over large distances unaided by anthropogenic factors.
In conclusion, M. orientalis was likely co-introduced in large numbers or via multiple introductions, based on the equally high genetic diversity in native and invaded regions. Thus, the invasion history from the native region along the reported invasion route of M. orientalis from our literature survey matches well with the resulting phylogeography based on COI haplotypes.
No differentiation in the putative native range of M. intestinalis. Based on our literature survey, the invasion of Atlantic coasts of Europe by M. intestinalis can be characterized by local introductions at seaports and harbor areas, and a slow natural spread afterwards. While M. intestinalis was thought to originate from the Mediterranean Sea Basins, our genetic data based on mtDNA could not confirm this as a source of the invasion. Even though the overall level of population genetic structure in the Mediterranean was significant, it was very low. Zooming in on this small difference, however, we did not detect a significant difference among Mediterranean Sea Basins based on hierarchical AMOVAs nor between any pairwise comparison between Mediterranean samples. The absence of clear population genetic differentiation in the Mediterranean is striking. If the Mediterranean Sea Basins were the native distribution range of M. intestinalis, the long geological history of the basins should have resulted in significant population structure. Although no uniform phylogeographic pattern in other marine taxa in the Mediterranean Sea Basins has been found 41 , for many species, the topography and oceanographic conditions of the Strait of Otranto and the Siculo-Tunisian Strait represent a break that typically causes population genetic differentiation between the Adriatic, Ionian and Ligurian Seas (e.g. [42][43][44]. Likewise, genetic breaks were found near Marseille, separating the Ligurian Sea and Gulf of Lion, especially for species with lower dispersal ability [44][45][46] . As the short-lived pelagic larvae of Mytilicola limit the potential for natural dispersal 12  www.nature.com/scientificreports www.nature.com/scientificreports/ A bottleneck in M. intestinalis may have occurred prior to its invasion, which would have reduced standing genetic variation [47][48][49][50] . Such a scenario is supported by the star-like haplotype network, relatively low genetic diversity and significant negative Tajima's D statistic that indicates recent population expansion or recovery from an ancient bottleneck after which diversity slowly recovered, creating a star-like haplotype network at the start. However, such a bottleneck must have happened fairly recent to prevent the build up of genetic differentiation over extended periods of time. Alternatively, anthropogenic shuffling of native genetic diversity might explain the lack of population genetic structure in M. intestinalis within the Mediterranean. Like M. orientalis' invasion, aquaculture transports might also be a cause here. However, the typically small size of translocated mussels ('mussel spat') contests this. The number of parasites is usually proportional to the shell length of mussels 51,52 and very young mussels are rarely infected 12 . Moreover, Mediterranean commercial mussels are usually grown in hanging cultures, which are considerably less vulnerable to parasite infections than mussels grown on the benthos. These mussels are also replaced yearly by a new parent stock, which makes it difficult for M. intestinalis to build up a dense population (reviewed in 12 ). Nevertheless, transportation of spat between Mediterranean regions does occur and could be contributing to or have contributed in the past to the homogenization of its parasite fauna. Other human activities related to mussel dispersal are ballast water and ship hull fouling 53 . Ship hull fouling was seen as the primary route by which M. intestinalis invaded the German and British coasts (e.g. 23 ) and explained its disjunctive distribution there (reviewed in 12 ). Higher resolution of more independent genetic markers might allow for detecting differentiation. Nevertheless, while anthropogenic shuffling could explain the lack of population genetic structure in the Mediterranean and Adriatic Seas, it cannot explain the overall low genetic diversity compared to M. orientalis.
Finally, M. intestinalis may be an introduced species in the Mediterranean Sea Basins, which would explain the lack of population genetic structure and the indications for population expansion. A few hints supporting this scenario can be found in literature. Mytilicola intestinalis was first described in 1902 and 1903 from two important ports for shipping, Trieste and Gravosa (both Adriatic Sea) 54,55 -around 40 years after the opening of the Suez Canal (November 1869), an event that led to the invasion of many marine species into the Mediterranean (Lessepsian migration, e.g. 56,57 ). In 1902, nearly all dissected mussels at those Adriatic sites were infected by M. intestinalis, with up to 50 parasites in one host 54,55 . Prevalences and intensities that high are typically seen in mussel populations that have recently invaded -e.g., maximum intensity of up to 14 parasites in Germany 10 , up to 59 in the UK 58 and up to 70 in the Netherlands 59 . In 1953, M. intestinalis could not be found in the same locations where A. Steuer initially described the species 60 , and we also did not find M. intestinalis in > 200 mussels from Trieste (Table 1). Similarly, in Mare Piccolo, Taranto (Italy), M. intestinalis had been a sudden pest before 1932, but authorities successfully eradicated it 61 . At the same site, a prevalence of 6.7% and mean intensity 1.8 (n = 106) was found in 1953 62 but in 1996, M. intestinalis was absent (n = 31) 63 . The cases of Trieste and Taranto both may be interpreted as outbreaks upon initial invasion followed by low infestation over longer time scales, resembling the typical boom-bust dynamics of biological invasions before populations enter the comparatively stable adjustment phase 64 . Moreover, Adriatic harbors that were rarely visited by large vessels had no or low M. intestinalis prevalences, while mussels at large ports were all heavily infested 65 . www.nature.com/scientificreports www.nature.com/scientificreports/ Because our data may not be consistent with the idea that M. intestinalis is native in the Mediterranean Sea Basins, we suggest that M. intestinalis should be regarded as a cryptogenic species (i.e., a species whose origin cannot readily be determined with the available data, sensu Carlton 66 ). Along European Atlantic and North Sea coasts, however, M. intestinalis is clearly an invasive species, as there are records of times when M. intestinalis was not present and later invaded. The absence of such data in the Mediterranean neither proves nor disproves presence or absence of M. intestinalis in the Mediterranean Sea Basins before its taxonomic description. The higher prevalence around shipping ports and aquaculture areas could indicate that M. intestinalis is invasive, or that its distribution has been influenced by humans also in the Mediterranean. Such a scenario would also imply that its native host species may not be Mytilus galloprovincialis. Solving the puzzle whether M. intestinalis is native or is an invader in the Mediterranean would need explorative, basic parasitological work in potential native regions such as the southern and eastern Mediterranean Sea Basins, the Red Sea and potentially further, in combination with a population genetic study.
Our study shows that despite the large body of literature existing on these two common marine parasites of commercially important bivalves, hypotheses derived from literature-based invasion histories reports do not necessarily match conclusions based on population genetic demography. This indicates that some basic understanding of the invasion processes is lacking, and highlights the usefulness of population genetic studies to test anecdotal evidence of invasion histories. The native region and invasion route of M. orientalis were validated with the sampling and choice of marker in the current study: the native region of M. orientalis is Japan, and the parasite hitchhiked with its host to North America and Europe. Whether M. orientalis was introduced to Europe via Canada, or directly from Japan, or both, and how M. orientalis spread throughout Europe remains to be answered with the use of higher resolution nuclear genetic markers. In contrast, the native region of M. intestinalis could not be validated, because of low genetic diversity and lacking population structure within and among the www.nature.com/scientificreports www.nature.com/scientificreports/ Mediterranean Sea basins, which may point to M. intestinalis being an invader in the Mediterranean Sea Basins (cryptogenic species). The invasion histories and demographics thus differ between the two Mytilicola species, and this may have profound impacts on their further invasion and evolutionary trajectories, as well as their competition in the invaded regions where they now co-occur.

Materials and Methods
Literature review. To describe the distribution of both invasive Mytilicola species through time in putative native and invaded regions, we searched the Web of Science/Web of Knowledge and Google Scholar with the keywords "Mytilicola intestinalis" and "Mytilicola orientalis". Further publications, including grey literature, were found by checking the references of each relevant publication (i.e., "snowballing"). All publications were screened for recorded absence or presence of Mytilicola species per host species, and for coordinates of sampling sites. This resulted in data sets for Mytilicola intestinalis and for Mytilicola orientalis (available upon request) which were used for plotting occurrence maps with the R packages rworldmap 67 and ggplot 68 . If latitude and longitude were not reported in a publication, we inferred approximate coordinates from maps so that the data could be plotted.   A possible caveat in the literature survey is that M. orientalis' invasion in Europe has been cryptic because of difficulties with reliable species identification 69 . In the literature survey, we used the identification presented by the authors, but that may be incorrect when studies are based on morphology, as the distinguishing characteristics for M. intestinalis and M. orientalis are not fully reliable 69 . For instance, M. intestinalis has been reported from Crassostrea gigas 70,71 , although it was shown not to be a suitable host for M. intestinalis 72,73 ; these observations are more likely to be of M. orientalis instead.
Sampling of host and parasite species. Different    www.nature.com/scientificreports www.nature.com/scientificreports/ intertidal zone during low tide at 61 locations along European coasts (Mediterranean Sea Basins, Atlantic, North Sea, and Baltic Sea), along American and Canadian Pacific coasts (Puget Sound, Vancouver Island) and along Japanese coasts (Seto Inland Sea, Sea of Japan, and Oyashio Current; Table 1). Mussel species were identified based on morphology only. There are no reports of the mussel species, B. pharaonis, being dissected for Mytilicola spp. Nevertheless, we included a sample from one location because if it is present in this species, this could extend the occurrence range of the parasite. Bivalves were dissected in the laboratory under a dissecting microscope (10-30×) or with the help of a magnifying glass (10×) in order to screen for the presence of Mytilicola individuals. Parasites were isolated from their host and stored separately by individual host in 100% ethanol. DNA extraction, amplification and sequencing. DNA was extracted from individual Mytilicola using either the QIAgen 96 DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany) or the Zymo Research genomic DNA kit (the latter by Baseclear B.V., Leiden, the Netherlands). One to five M. intestinalis and one to eight M. orientalis per host were sequenced. All extractions were performed according to the protocols provided by the manufacturers. DNA concentrations were measured on a Nanodrop to confirm DNA quality and quantity. A 534 base pair (bp) fragment of the mitochondrial cytochrome-c-oxidase I region (COI) was amplified using the primer pair MOICOIf and MOICOIr for M. intestinalis 69  Data analyses. We report parasite prevalence as the number of individual mussels or oysters infected with M. intestinalis or M. orientalis divided by the total number of mussels or oysters examined per location, and mean intensity as the mean number of M. intestinalis or M. orientalis within an infected mussel or oyster per location (sensu Margolis et al. 75 ). Part of the M. orientalis sample from Tavira (Portugal) was acquired through an anti-Mytilicola treatment of mussels using Dichlorvos Pestanal (DDVP, 30 mg per 1 liter of seawater) 13,76 . Fifteen live Mytilicola specimens escaped the mussels during the recovery period shortly after the 4 h treatment. Prevalence and mean intensity (Table 1) for Tavira were calculated based on dissected mussels only.
Quality control of the sequences was performed in CLC Genomics Workbench 8.5.1 (for M. orientalis) or BioEdit 7.2.5 77 (for M. intestinalis), after which forward and reverse reads were assembled. We screened all chromatograms for double peaks and at least one of the forward and reverse sequences had no double peaks at any base in the sequence. Some of the North American samples failed to amplify or sequence entirely and were not considered further. COI sequences were cropped to 483 base pairs (bp) for M. intestinalis and to 476 bp for M. orientalis. Unique haplotypes were determined using DnaSP 5 (M. orientalis) 78 or using a custom Python 3.6 script (M. intestinalis) (GitHub: https://github.com/pluttik/collarl). Arlequin 3.5.1.2 79 was used to calculate pairwise population comparisons (10,000 permutations), population differentiation based on the mitochondrial F-statistic Φ ST and conventional F-statistics F ST (10,000 permutations), Analyses of MOlecular VAriance (AMOVA), molecular diversity indices, Tajima's D statistic 80 and to perform Fu's F s tests 81,82 for selective neutrality (10,000 simulated samples), mismatch distributions (10,000 permutations) and minimum spanning networks for both species. For neutral markers, Tajima's D and Fu's F s can be used to detect changes in population size. Significantly negative D and F s values can be interpreted as signatures of population expansion or, unless another independent marker has the same pattern, a past selective sweep 81,82 . We calculated both Φ ST and F ST to detect any signs of possible population differentiation, as well as for easier comparison to other studies that use either or both fixation indices. In our analyses we only included locations with more than 15 Mytilicola individuals, with the exception of calculations for minimum spanning networks, in which all individuals were included. In addition to overall AMOVAs, we also calculated AMOVAs within putative native and introduced regions for both species, and hierarchical AMOVAs for Mediterranean M. intestinalis with two designs to test whether the overall population differentiation could be attributed to different regions within the Mediterranean Sea Basins. The first design tested the Gulf of Lion (Cassis, Istres, Rouet and Bouzigues) vs. the Ligurian Sea (Livorno) vs. the Adriatic Sea (Malinska and Marjan). In the second design, the location Cassis was added to the Ligurian Sea group instead of the Gulf of Lion group. This was done to account for the possibility that the location Cassis, which is east of Marseille and a known location for genetic breaks [44][45][46] , belongs to the Italian group rather than to the group with more proximate locations within France. To test for differences in genetic diversity between species and regions, we fitted a linear model testing haplotype diversity (h) and nucleotide diversity (π) as a function of species (levels: M. intestinalis and M. orientalis), region (levels: putative native region, introduced region) and the interaction term species × region in the R statistical environment version 3.4.3 83 .

Data Availability
Sequencing data is available on Genbank under the accession numbers MN334483-MN334508 for M. orientalis and MN334509-MN334526 for M. intestinalis.