Introduction

Biological invasions (whether recognized or cryptic) can alter ecosystem functioning and services1, and can have ecological and evolutionary impacts on native communities and ecosystems (e.g.2,3,4,5). Invasive species can also be parasites, which may be introduced with or without their hosts. Parasite invasions can have wide-ranging effects on host-parasite interactions with strong ecological and evolutionary implications6. To examine the effects of invasive parasites, a detailed understanding of their origins and invasion pathways is important. To this extent, the traditional approach of using literature records can be fruitfully complemented by population genetic analyses. Phylogeographic studies of species in their putative native and invaded regions (“invasion genetics studies”)7 can provide more detail to clarify the origin, invasion routes, population genetic structure and genetic diversity of invasive parasites8. However, invasive parasites are greatly understudied in this respect, as well as with regard to their ecological and evolutionary impacts9.

Two exceptions are the consecutive invasions of the parasitic copepods Mytilicola intestinalis Steuer 1902 and Mytilicola orientalis Mori 1935 along the European Atlantic shores. Both of these natural experiments have been proven to be particularly useful to study their ecological consequences and coevolution with the blue mussel, Mytilus edulis, one of their host species. Mytilicola intestinalis invaded European waters first in the 1930s presumably through ship hull fouling and relaying of infected mussels (e.g.10,11,12). Its invasion into the western European Wadden Sea was recently used as a model system to study host-parasite coevolution in the wild13,14. 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 tolerance13. These trajectories were characterized by specific molecular interactions between host and parasite14. Mytilicola orientalis was co-introduced with the invasive Pacific oyster, Crassostrea gigas, through aquaculture transports15,16, and now also infects a suite of native hosts, including blue mussels, M. edulis16,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 known19,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. intestinalis21,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 similar25. The two free-living larval stages, nauplius and metanauplius, last 36 to 48 hours at a laboratory temperature of 18 °C23. 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 field23. 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 limited12,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.

Results

Literature-based invasion history

A search for “Mytilicola intestinalis” resulted in 72 hits on the Web of Science/Web of Knowledge (all databases) and 934 hits on Google Scholar (on 15 December 2017). A search for “Mytilicola orientalis” resulted in 15 hits on the Web of Science/Web of Knowledge (all databases), and 373 hits on Google Scholar (on 15 December 2017). More literature was found by snowballing (see Methods). Our literature survey led to the inclusion of 173 publications on occurrences of M. intestinalis (2187 data points including those of this study; reference list in Supplementary Information) and of 47 publications on occurrences of M. orientalis (416 data points including those of this study; reference list in 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.

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. intestinalis27. 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.

Figure 1
figure 1

Temporal development of the known distributions and invasions of (A) Mytilicola intestinalis and (B) Mytilicola orientalis from time of first observation and species description (i.e., 1902 for M. intestinalis and 1935 for M. orientalis) until December 2017. Locations where Mytilicola spp. presence has been recorded are indicated with red dots, showing the species distribution, invasion route and spread, as well as research effort through time. For temporal development of Mytilicola spp. invasions decade-by-decade, including species absence records, see Supplementary Figs S1 and S2. For animated temporal development of Mytilicola spp. year-by-year, see the .gifs in Online Resources 1 and 2. This figure is based on data from this study and from data acquired through our literature surveys. See Supplementary Information for reference lists.

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 areas26. 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 France28. 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 locations along Japanese coasts, 7 locations along Pacific American and Canadian coasts, and 21 locations along European Atlantic, North Sea and Baltic coasts resulted in detection of M. intestinalis and M. orientalis at 20 and 30 locations, respectively. In Japan, M. orientalis has only been previously recorded in Hiroshima29 and Takehara30, both in the Hiroshima Prefecture, and furthermore in Hiwasa, the Tokushima Prefecture31 – the other locations, including those outside the Seto Inland Sea are, thus, new records. Also, Tavira (Portugal) is a new location for M. orientalis. The Mediterranean locations Malinska, Marjan (both Croatia) and Cassis (France) are new records for M. intestinalis, although Cassis is close to locations where M. intestinalis was already known from (i.e., Berre Lagoon and Marseille32,33).

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).

Table 1 Parasitic copepods Mytilicola intestinalis (MI) and Mytilicola orientalis (MO) collected from mussels: Mytilus galloprovincialis (Mg), Mytilus edulis (Me), Mytilus trossulus (Mt) and Brachidontes pharaonis (Bp), and from oysters: Pacific oysters Crassostrea gigas.

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 FS 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).

Figure 2
figure 2

Minimum spanning networks for COI haplotypes (circles) in (A) Mytilicola intestinalis and (B) M. orientalis. Branch length is proportional to the number of nucleotide differences (the shortest branches represent one nucleotide difference). The black dots on the branches represent missing haplotypes. Frequency of observation is proportional to circle area. Colors represent large scale differentiation between clusters of haplotypes and correspond with the pie charts in Fig. 4. Figure 2 was drawn in Adobe Illustrator based on the Arlequin output list of OTU differences.

Table 2 Tajima’s D and Fu’s Fs neutrality tests for partial cytochrome-c-oxidase 1 sequences in Mytilicola intestinalis and M. orientalis, based on 10,000 simulations.
Figure 3
figure 3

Mismatch distributions of Mytilicola intestinalis (top row) and Mytilicola orientalis (lower row). From left to right: mismatch distributions based on population differences (M. intestinalis: all locations pooled; M. orientalis: Iwaya separate from the rest of the locations), mismatch distributions in the putative native region and in the invaded region(s). The sharp decline in M. intestinalis is indicative of population expansion after a bottleneck or selective sweep; the bimodal distributions in M. orientalis are indicative of stable demographics.

Figure 4
figure 4

Maps showing pie charts of cytochrome-c-oxidase 1 (COI) haplotype frequencies in the invasive parasite species Mytilicola intestinalis and Mytilicola orientalis. Shown plotted overlaying the maps are pie charts for the different sampling locations, while the large pie charts in the corner of the maps represent the grouped frequencies for all locations in that map. On the side, the large pie charts represent the grouped frequencies for all invaded European locations for M. intestinalis and M. orientalis (“Europe”), and for all invaded locations for M. orientalis that includes the North American samples (“Invaded”). Note that the pie charts from the Mediterranean and Japan represent the overall frequencies of all sampled native sampling locations for M. intestinalis and M. orientalis, respectively. Pie charts are smaller for sampling locations with fewer than 15 individuals, which were not included in our analyses. Colors correspond to large scale differentiation between clusters of haplotypes in the minimum spanning networks presented in Fig. 2. Pie charts were made in R statistical environment version 3.4.383, maps were plotted with R packages rworldmap67 and ggplot68 and Fig. 4 was assembled in Adobe Illustrator.

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 FS values (Table 2), and by the mismatch distributions (Fig. 3). Only one population (Puan Klent, Germany), representing the northern invasion front of M. orientalis at time of collection, stood out with significant negative Tajima’s D and Fu’s FS values (Table 2).

Because of the long geological history of the Mediterranean Sea Basins, its different hydrographic regimes and the short dispersive larval phase of M. intestinalis, we expected a high genetic diversity of M. intestinalis in this region and population genetic structure between the different basins. The overall level of population differentiation for M. intestinalis was estimated at ΦST = 0.01468 (p = 0.0108) and FST = 0.00848 (n.s.), with an overall population differentiation in the Mediterranean Sea Basins of ΦST = 0.01531 (p = 0.0323) and FST = 0.02397 (p = 0.0325), and in the invaded region of ΦST = 0.01081 (n.s.) and FST = −0.00544 (n.s.) (Table 3).

Table 3 Results from analyses of molecular variance (AMOVA) on partial cytochrome-c-oxidase 1 (COI) sequences in Mytilicola intestinalis and M. orientalis with mitochondrial fixation index ΦST and with conventional F-statistics.

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 FST = 0.03002 (p < 0.00001), with an overall population differentiation in Japan of ΦST = 0.10287 (p < 0.00001) and FST = 0.07019 (p < 0.00001), and in the introduced region of ΦST = 0.02243 (p = 0.0303) and FST = 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 FCT = 0.02069 (n.s.), FSC = 0.00825 (n.s.) and FST = 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.

Table 4 Results from hierarchical analyses of molecular variance (AMOVA) based on pairwise differences, on partial cytochrome-c-oxidase 1 (COI) sequences in Mediterranean Mytilicola intestinalis with mitochondrial Φ-statistic and conventional F-statistic.

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 Table S5; FST = 0.165–0.216, Supplementary Fig. S3, Supplementary Table S6). Only the pairwise population comparison between Iwaya (Fukuoka Prefecture) and Tarumi (Kagawa Prefecture) was not significant after Bonferroni correction (ΦST = 0.220, p = 0.00129 > Bonferroni corrected p-level of 0.000476, Supplementary Table S5). However, with conventional F-statistics, also this pairwise comparison was significantly different (FST = 0.191, p = 0.00040 < Bonferroni corrected p-level of 0.000476, Supplementary Table S6). No genetic differentiation was detected between the locations sampled along the Seto Inland Sea, i.e., Takahara, Ondo, Niihama and Tarumi. All pairwise population comparisons between the Sea of Japan (Iwaya) and locations from the invaded regions in North America and Europe were significant after Bonferroni correction and showed strong genetic differentiation (ΦST = 0.228–0.344, Supplementary Table S5; FST = 0.187–0.247, Supplementary Table S6). No significant pairwise population differences were detected between locations in the invaded regions.

We did not detect significant differences in haplotype diversity (h) between the putative native regions and the introduced regions (ANOVA, F1, 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, 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, F1,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: F1, 26 = 99.2662, p < 0.00001; ANOVA π: F1, 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: F1, 26 = 0.5634, p = 0.460; ANOVA π: F1, 26 = 1.3178, p = 0.261), indicating that genetic diversity differed consistently between species.

Figure 5
figure 5

Haplotype diversity h and nucleotide diversity π of Mytilicola intestinalis (MI) and Mytilicola orientalis (MO). Putative native ranges are indicted in grey and invaded ranges in black. Bars represent standard deviations. Only samples with more than 15 individuals were included. See Supplementary Tables S1 and S2 for π and h values of individual locations underlying this figure.

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.

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 ongoing34 (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 Sea35 contributing smaller amounts35. 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 regions36. Mytilicola orientalis hitchhiked with oyster transports to North America16 and to Europe15,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.

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 host19,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 bivalves18 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 197939 and in the northern Adriatic40, 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 found41, 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 ability44,45,46. As the short-lived pelagic larvae of Mytilicola limit the potential for natural dispersal12,23, population differentiation among the long-established Mediterranean populations should have been maintained. Especially, since breaks in the distribution of mussel hosts introduce additional barriers for M. intestinalis larval dispersal. This limited dispersive ability of M. intestinalis can also be observed in the the literature-based invasion data (Resource 1), showing a slow spread upon invasion in the Dutch Wadden Sea and Scottish lochs.

Low mtDNA variability masks the origins of M. intestinalis

There are several possible explanations for the observed absence of population genetic differentiation between the Mediterranean M. intestinalis samples and the overall low genetic diversity compared to M. orientalis.

A bottleneck in M. intestinalis may have occurred prior to its invasion, which would have reduced standing genetic variation47,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 mussels51,52 and very young mussels are rarely infected12. 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 in12). 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 fouling53. 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 in12). 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 host54,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 Germany10, up to 59 in the UK58 and up to 70 in the Netherlands59. In 1953, M. intestinalis could not be found in the same locations where A. Steuer initially described the species60, 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 it61. At the same site, a prevalence of 6.7% and mean intensity 1.8 (n = 106) was found in 195362 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 phase64. 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 infested65.

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 Carlton66). 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 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 rworldmap67 and ggplot68. 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 identification69. 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 reliable69. For instance, M. intestinalis has been reported from Crassostrea gigas70,71, although it was shown not to be a suitable host for M. intestinalis72,73; these observations are more likely to be of M. orientalis instead.

Sampling of host and parasite species

Different mussel species (Mytilus edulis, Mytilus galloprovincialis, Mytilus trossulus and Brachidontes pharaonis), and Pacific oysters (Crassostrea gigas) were collected in the 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. intestinalis69 and a newly developed primer pair for M. orientalis (MoriCOI50f: 5′-TTG ATC GGG CTT AAT TGG-3′ and MoriCOI50r: 5′-GAT CGG TTA ARA GCA TGG T-3′). The M. orientalis primers were developed in CLC Genomics Workbench 8.5.1 (CLCbio, Denmark, https://www.qiagenbioinformatics.com/) and were based on M. orientalis COI sequences from Elsner et al.74.

For M. intestinalis, polymerase chain reaction (PCR) and Sanger sequencing using both forward and reverse primers were performed according to Goedknegt et al.69 at Baseclear B.V. (Leiden, Netherlands). For M. orientalis, PCR was performed with HotStarTaq (Qiagen, Hilden, Germany) in 30 μl reaction volumes, adding 15 μl HotStarTaq Mastermix, 1.5 μl of each 10 μM primer, 11 μl H2O and 1 μl of sample for each reaction. Amplification started with a heat-activation of 15 min at 95 °C, followed by 38 cycles of 95 °C for 30 s, annealing at 50 °C for 30 s and elongation at 72 °C for 60 s. The final extension step was at 72 °C for 10 min. Amplification success of all M. orientalis samples was checked with a QIAxcel electrophoresis system (Qiagen, Hilden, Germany). PCR purification and Sanger sequencing of M. orientalis samples were performed on both the forward and reverse strand at the Institute of Clinical Molecular Biology (Kiel, Germany).

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.577 (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.279 was used to calculate pairwise population comparisons (10,000 permutations), population differentiation based on the mitochondrial F-statistic ΦST and conventional F-statistics FST (10,000 permutations), Analyses of MOlecular VAriance (AMOVA), molecular diversity indices, Tajima’s D statistic80 and to perform Fu’s Fs tests81,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 Fs can be used to detect changes in population size. Significantly negative D and Fs values can be interpreted as signatures of population expansion or, unless another independent marker has the same pattern, a past selective sweep81,82. We calculated both ΦST and FST 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 breaks44,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.383.