Population genetics

To study the population genetics of Haminoea japonica in the native and non-native ranges, DNA was extracted, amplified, and sequenced for the mitochondrial gene cytochrome c oxidase subunit I (COI) as described in Hanson et al. [8]. Thirty-four individuals from Japan and eleven individuals from Europe were added to the samples studied by Hanson et al. [8] and Crocetta et al. [9]. This resulted in a total sample size of 142 individuals, spanning the entire known range of the species. Sequences were assembled and edited using the software Geneious Pro 4.8.3 [14]. Geneious Pro 4.8.3 was also used to extract the consensus sequences and to align them using the default parameters. A haplotype network was constructed using TCS 1.21 [15] with a 95% connection limit. To assess the potential for missed haplotypes due to sample size in both the native (Japan and Korea) and non-native (North America and Europe) range, rarefaction analyses were run in EstimateS 8+ [16] with each location considered a single sample, even though collections may have been made at each location in more than one occasion. North American and European samples were run together under a single analysis since the low number of samples from Europe would make the identification of the rarefaction curve asymptote impossible. Chao2 means and 95% lower and upper confidence intervals were estimated.

Based on geographic distribution, five groups (Northeastern Japan, Southern Japan, Sea of Japan, North America, and Europe) and 27 populations (see Table 1) were tested for genetic structure using the analysis of molecular variance (AMOVA) implemented in Arlequin 3.5 [17]. Three AMOVA analyses were run for all localities, native range localities and non-native range localities. Significance of the AMOVAs was tested using 16,000 permutations of individuals between groups. Arlequin 3.5 was also used to calculate F_ST values as a measure of pairwise differences between all population groups (Table 1) as well as the nucleotide diversity (π) and haplotype diversity (h) of each population. The significance of the pairwise F_ST value was estimated by performing 16,000 permutations. The F_ST analysis was conducted with population groups rather than populations because of two reasons, the large number of populations would result in lower confidence values after Bonferroni’s correction and the presence of many populations with one or a small number of haplotypes would render the results uninformative or difficult to interpret.

Phylogenetic analysis

To gain insight into the geographic origins of H. japonica, the species was placed within a larger phylogeny of the genus Haminoea. For this part of the study DNA was extracted, amplified, and sequenced for the mitochondrial gene cytochrome c oxidase subunit I (COI) and nuclear gene 28S rRNA, as described by Malaquias et al. [18]. Sequences were verified by both forward and reverse comparisons and assembled and edited using Sequencher 5.01 (Gene Codes Corp.). All sequences were deposited in GenBank (Table 2). Sequences were aligned using Clustal_X [19] with a slow-accurate method and default settings. The alignments were further optimized by eye using MacClade 4.06 [20]. The best-fit models of evolution (TVM + I + G for COI and GTR+ I + G for 28S rRNA) were chosen using the Akaike information criterion [21] implemented in ModelTest 3.6 [22]. Aligned sequences were trimmed to exclude ambiguous reads at each end and gap-rich regions where homology could not be confidently assessed. A total of 522 bp of COI and 720 bp of 28S rRNA remained for use in phylogenetic analyses comprising 44 sequences. A combined analysis of the two genes (total 1242 bp) was conducted using Bayesian inference analysis in MrBayes 3.1.2b with default priors [23,24] for 1.5 x 10⁶ generations, with a sampling frequency of 100 and two separate runs to ensure that independent analyses were converging on the same tree. Convergence of runs was diagnosed using the program Tracer 1.4 [25]. For each analysis the first 3750 trees were discarded (‘burn-in’ period). Robustness of each node was assessed using Bayesian posterior probabilities calculated by MrBayes 3.1.2b. Two members of the family Bullidae (Bulla occidentalis A. Adams, 1850 and Bulla striata Bruguière, 1792) were chosen as outgroup taxa.

Assessment for parasite infection

In order to test for a possible Japanese origin of the parasites associated with H. japonica in San Francisco Bay, 108 individuals collected in Uchiura Bay, Japan and 30 individuals collected in Matsushima Bay, Japan were assessed for schistosome parasite infections. No permits are necessary to study and/or collect Haminoea japonica on the two sites. A permit from the local fishermen’s association is required to collect commercially valuable species, but H. japonica is not commercially valuable or an endangered species.

After collection, specimens were kept in containers filled with seawater and exposed to bright, direct sunlight for 1-2 hours. The specimens were then removed and the water in the container was examined for cercaria under a stereomicroscope or with the naked eye. Additionally, a subset of 72 specimens from the 108 specimens collected in Uchiura Bay, were dissected and examined for schistosome infection of the digestive gland.

Phylogeny of Haminoea

The phylogeny of Haminoea shows H. japonica as monophyletic and sister to Haminoea natalensis (Krauss, 1848), a species described from the Indian Ocean coast of South Africa. These two species are placed within a highly supported (posterior probability = 0.96) clade of Southeast Asian species (Figure 2). Sister to this is a clade containing seven species, six of which are Atlantic or Mediterranean in origin, and one Eastern Pacific (Haminoea vesicula).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Bayesian molecular phylogeny produced from the combined mitochondrial COI and nuclear 28S rRNA genes.

Numbers above branches are posterior probabilities. Outgroup taxa removed for clarity.

https://doi.org/10.1371/journal.pone.0077457.g002

Population genetics

The AMOVA results (Table 3) revealed that when all populations are analyzed together the majority of genetic variation (50.9%) occurred within populations, whereas the variation among groups was 36.74%. The AMOVA run based on native range populations only gave similar results, with most of the variation (52.87%) within populations and only 32.7% among groups. However, the AMOVA run for non-native populations only, calculated much higher genetic variation among populations within groups (46.23%) and very little among groups (15.43%). The F_ST analysis results (Table 4) confirmed that all population groups examined are genetically distinct from one another with the exception of Northeastern Japan, which is not significantly different from North America and Europe.

Thirty-three haplotypes were identified in the native range, and four in the non-native range (Figure 3). Two of the four non-native haplotypes were found in Japan, whereas two (H24 – from a single individual collected from Tomales Bay and H25 – widespread in Spain and Italy) were not detected in the native range. Rarefaction analysis of the non-native range haplotypes predicted that four haplotypes should exist in the non-native range (Chao2 mean and 95% CI lower bound), with a maximum of 4.46 (95% CI upper bound). This suggests that all the haplotypes present in the non-native range have been recovered, although further sampling may reveal one additional haplotype. Because both the sample data and the Chao2 estimate curves reached an asymptote (Figure 4C), these results can be considered reliable. The same analysis run for the northeastern Japan populations suggested that a minimum of 12.59 (lower 95% CI) and a maximum of 176.22 (upper 95% CI) with a Chao2 mean of 33.5 haplotypes should be present. Although this maximum is not as reliable as that for the non-native range analysis, since the rarefaction curves did not reach an asymptote (Figure 4B), the minimum 12.59 estimate reveals that a substantial sampling effort is necessary to completely recover the genetic diversity existing in this region. Finally, a rarefaction analysis for the rest of the native range revealed that between 51.45-484.5 (135.25 Chao2 mean) haplotypes should be present in the rest of Japan and Korea, meaning at least 31 more haplotypes could be detected with further sampling effort (Figure 4A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Haplotype network and distribution map of records of H. japonica in Japan, Korea, Europe and North America.

Numbers in the map indicate the localities listed in the key to localities. Numbers in the haplotype network indicate haplotype number as in Table 1. Haplotypes marked with an asterisk have been newly found in this study.

https://doi.org/10.1371/journal.pone.0077457.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Rarefaction analyses of COI haplotype diversity with sampling effort in H. japonica.

Black circles = Chao2 mean, Black squares = Sobs. A. Native range excluding Northeastern Japan. B. Northeastern Japan. C. Non-native range.

https://doi.org/10.1371/journal.pone.0077457.g004

Parasite infection rates in Japan

None of the individuals examined for cercariae, either through shedding or dissection, were found to be infected with schistosome parasites.

Origin and population structure of the non-native populations of H. japonica

The phylogenetic analysis of Haminoea suggests an Indo-Pacific origin of H. japonica [7] as this species is nested in a clade containing only Indo-West Pacific species (Figure 2). However, the placement of H. japonica as sister of H. natalensis is probably an artifact of the limited taxon sampling within the Indo-pacific species examined. Additionally, the F_ST analysis results (Table 4) confirmed the population differentiation and structure revealed with a smaller sample size [8]. This is consistent with the hypothesis of a northeastern Japan origin for the non-native populations.

The results of the AMOVA analyses are indicative of substantial gene flow between populations in the native range resulting in most of the genetic variation occurring within populations (52.87%) and very little among populations (14.43%). However, there is substantial variation among groups (32.70%), suggesting a certain degree of isolation between the Sea of Japan, Southern Japan and Northeastern Japan populations of H. japonica.

As suggested below, the movement of bivalves between regions in Japan may have resulted in the introduction of southern haplotypes of H. japonica to Mangokuura in Northeastern Japan. Introductions of H. japonica into the non-native range have been associated with commercial imports of bivalves, particularly oysters and clams [7,8,11], but this new evidence suggests the same vector have resulted in movements of H. japonica within the native range, but between distinct populations. Such movements may have also brought haplotypes of diverse origins to Tokyo Bay, resulting in lower percentages of variation among groups than within populations. However, considering the high intensity of historic and recent transport of bivalves between regions of Japan and Korea [26], the impact of these within-Japan introductions in the present genetic structure of H. japonica is surprisingly low, as revealed by the relatively high percentage of variation among groups. This suggests that introductions of H. japonica are probably rare.

The genetic structure in the non-native range is very different, with very little variation among groups (15.43%) and high percentages of differentiation among populations and within populations (46.23% and 38.34% respectively). The genetic uniformity between Europe and North America suggest that these populations are the result of either a single or a small number of introduction events, whereas the high percentage of variation among populations is likely the result of founder effects.

Sample size and rarefaction curves

The rarefaction analysis confirms that this study has recovered all or most of the COI haplotypes found in the non-native range of H. japonica (only one potentially remains unrecovered). However, a number of haplotypes from the native range remain undetected by our sampling effort. This would explain the fact that haplotypes H24 and H25 are missing from our sample in the native range. The present analysis confirms that the reduced genetic diversity in the non-native range is not an artifact of a small sample size, but it is indicative of severe population bottleneck, which is consistent with a relatively recent introduction from a single or a small number of sources. Further evidence of population bottleneck in the non-native populations is provided by the low levels of both nucleotide and haplotype diversity in North America and Europe (Table 5).

Invasion pathway

With the available evidence it appears that the introductions of H. japonica in North America and Europe are independent events. The presence of one unique haplotype in each North American and European populations (H24 and H25 respectively) suggest independent introductions. If the European populations were the result of a secondary introduction from North America it would be difficult to explain the presence of the common European haplotype H25, which has not been detected in North America. Although these two haplotypes (H24 and H25) have not been detected in Japanese samples, suggesting they are rare in the native range, very similar haplotypes were recovered in Northeastern Japan (H29 and H32). The European haplotype H25 is common (frequency 36%) and widespread (found in Spain and Italy), this is likely the result of founder effect. On the contrary the North American haplotype (H24) is rare (frequency 4%).

Secondary invasions within Japan

The presence in Tokyo Bay of several haplotypes commonly found in other regions of Japan was suggested to be the result of introductions of slugs from other populations [8]. In this paper for the first time we report the presence of haplotypes commonly found in Southern Japan in Northern Japan samples, but only in the coastal lagoon of Mangokuura in Miyagi Prefecture. Mangokuura is a major hub for bivalve culture and oyster spat exports and during the 1980-1990s received imports of oyster spat from a few locations in the Seto Inland Sea (Southern Japan) including Hiroshima and Okayama (Kenji Okoshi, pers. comm.). It is a distinct possibility that southern Japanese haplotypes of H. japonica could have been introduced into Mangokuura as a consequence of these activities. These secondary introductions are difficult to detect in the absence of molecular data and could be problematic, as alleles from non-native populations can spread through hybridization potentially reducing fitness of native populations.

Potential for further spread

Previous work suggested that cold winter ocean temperatures could limit the potential range of non-native populations of H. japonica [8]. This hypothesis was based on two observations: 1) The non-native populations originated in a region influenced by the Oyashio cold current; 2) All known non-native records by the time of the publication of the paper were from regions with relatively cold water temperatures such as northern California and Washington, USA, the Atlantic coasts of Spain and France, and the northern Adriatic Sea. This was corroborated by the assumption that the Mediterranean populations have become extinct [11] due to absence of recent records even in well-sampled areas such as the northern Adriatic Sea lagoons. However, evidence of established populations of H. japonica in coastal lagoons along the western coast of Italy [9] suggested that the potential spread of this species could be much greater. In the present paper we report H. japonica for the first time in the Mediterranean coast of France as well as in three new localities in Italy. It appears that H. japonica is now well established in the western Mediterranean and central Adriatic and that milder winter temperatures do not appear to have prevented further spread. At this point all coastal lagoons in the Mediterranean as well as in Southern California should be considered vulnerable to invasion by H. japonica. It is possible that the range of H. japonica is already much larger than currently known, but it remains undetected due to difficulties in identification. The main external diagnostic characteristic of H. japonica is the presence of a deep notch in the cephalic shield [7], which is also present in the Atlantic and Mediterranean species H. hydatis [27], making these two species difficult to distinguish in the absence of molecular or anatomical data. These morphological similarities hamper proper identification and early detection of H. japonica.

Potential ecological consequences of the invasion of H. japonica

Although based on anecdotal evidence [8], H. japonica seems to have displaced the Pacific Northwest species H. vesicula in British Columbia, suggesting it could be an aggressive invader. This is paralleled by a similar observation held in Laguna di Sabaudia [10] where H. japonica has become a dominant organism. The spread of H. japonica could have severe consequences for the biodiversity of Haminoea in the Mediterranean, where a number of species with restricted ranges occur in coastal lagoons. For example, Haminoea fusari is only known from Lago Fusaro, Italy, a coastal lagoon currently occupied by an established population of H. japonica [9]. Other species such as Haminoea hydatis, H. navicula, H. orbignyana, H. orteai and H. templadoi have been reported in southern European coastal lagoons, often in association with green algae (Ulva) and sea grasses [27,28,29,30,31] and are potentially vulnerable to the spread of H. japonica, which occupies a similar niche. Neither the low genetic diversity of non-native populations of H. japonica, nor their cold-water origin, seem to have hindered the rapid spread of this species in warmer Mediterranean waters and costal lagoons. Although H. japonica is currently considered a non-native species, it is possible that it has already become an aggressive invasive species and efforts to control it could become increasingly difficult.

Potential human health consequences of the invasion of H. japonica

Examination of specimens in Japan did not reveal the presence of schistosome parasites. These results do not reveal anything new about the origin of the cercarial parasite associated with H. japonica in San Francisco Bay, but fail to provide evidence of a Japanese origin. The apparent absence of parasites associated with H. japonica in Washington State (USA), Canada and Europe seems to confirm that the San Francisco case is unique. It is thus likely that the parasite has traveled independently from H. japonica into San Francisco Bay [12] resulting in a serendipitous association. Regardless of the origin of the cercarial parasite, this association makes northern California populations of H. japonica susceptible to spread the parasite to nearby locations. Since water temperatures did not appear to prevent the spread of H. japonica in the Mediterranean, it is unlikely to do so in California. Thus, southern Californian estuaries, some of which are heavily dependent on water-related tourist activities, may be vulnerable to invasion by both the slug and the parasite, potentially resulting in economic loss.