Invasion Pathway of the Ctenophore Mnemiopsis leidyi in the Mediterranean Sea

Gelatinous zooplankton outbreaks have increased globally owing to a number of human-mediated factors, including food web alterations and species introductions. The invasive ctenophore Mnemiopsis leidyi entered the Black Sea in the early 1980s. The invasion was followed by the Azov, Caspian, Baltic and North Seas, and, most recently, the Mediterranean Sea. Previous studies identified two distinct invasion pathways of M. leidyi from its native range in the western Atlantic Ocean to Eurasia. However, the source of newly established populations in the Mediterranean Sea remains unclear. Here we build upon our previous study and investigate sequence variation in both mitochondrial (Cytochrome c Oxidase subunit I) and nuclear (Internal Transcribed Spacer) markers in M. leidyi, encompassing five native and 11 introduced populations, including four from the Mediterranean Sea. Extant genetic diversity in Mediterranean populations (n = 8, N a = 10) preclude the occurrence of a severe genetic bottleneck or founder effects in the initial colonizing population. Our mitochondrial and nuclear marker surveys revealed two possible pathways of introduction into Mediterranean Sea. In total, 17 haplotypes and 18 alleles were recovered from all surveyed populations. Haplotype and allelic diversity of Mediterranean populations were comparable to populations from which they were likely drawn. The distribution of genetic diversity and pattern of genetic differentiation suggest initial colonization of the Mediterranean from the Black-Azov Seas (pairwise F ST = 0.001–0.028). However, some haplotypes and alleles from the Mediterranean Sea were not detected from the well-sampled Black Sea, although they were found in Gulf of Mexico populations that were also genetically similar to those in the Mediterranean Sea (pairwise F ST = 0.010–0.032), raising the possibility of multiple invasion sources. Multiple introductions from a combination of Black Sea and native region sources could be facilitated by intense local and transcontinental shipping activity, respectively.


Introduction
Introduction of non-indigenous species (NIS) beyond their native range is considered a principal threat to marine ecosystems worldwide [1]. The rate of such introductions accelerated in the past few decades in conjunction with increased maritime shipping and global trade [2][3]. Maritime traffic often involves use of ballast water loaded in source ports and later discharged in destination ports, resulting in mass transfer of organisms between distant regions [4][5][6]. Species with planktonic life stages have a high chance of interfacing with a shipping vector when ballast water is loaded, and thus of being moved around the world to new locations [7].
In recent years, gelatinous zooplankton outbreaks have raised concerns regarding the health of aquatic ecosystems [8]. A number of biological traits of gelatinous zooplankton may contribute to global outbreaks by this group. For example, many gelatinous zooplankton have a broad diet, high growth rate, high fecundity, high regeneration, encystment, and even reverse development potential [9][10][11], which enable them to overcome harsh conditions associated with the transport vector (i.e. ballast tanks) and successfully reach and establish in new environments [12][13][14].
The Mediterranean Sea has an enormously rich native biodiversity, though it is also the world's most invaded marine ecosystem [15][16] and is considered at very high risk of future invasions from ballast water discharges [17] and, especially, canal connections [18][19]. A total of 986 NIS have been recorded in the Mediterranean Sea, including 48 new species since 2011 [19]. The eastern section of the Sea has accumulated a disproportionate number of these NIS, principally due to Lessepsian invaders [18][19][20] that colonized following opening of the Suez Canal with its link to the Indian Ocean.
Knowledge of the source and pathways of NIS introductions is essential for developing management strategies to prevent invasions. A focus on areas at high risk of biological invasions is crucial and should be considered a management priority [17,[21][22]. In this paper, we explore the spread of the ctenophore Mnemiopsis leidyi A. Agassiz 1865 to the Mediterranean Sea. Mnemiopsis leidyi is native to the western Atlantic Ocean from Massachusetts, USA to Argentina. The species is a simultaneous hermaphrodite capable of self-fertilization, may reach maturity at two weeks of age, and can release up to 10,000 eggs per day [23]. Over the past 30 years, the species spread across Europe in a remarkable series of invasions, first entering the Black Sea (and Azov Sea) in early 1980s [24], the eastern Mediterranean in early 1990s (mainly Aegean Sea where an established population was not reported, [25][26]), followed by the Caspian Sea in 1999 [27].
Blooms of M. leidyi were reported throughout the Mediterranean Sea in 2009, from eastern to western coastal areas [28][29][30][31]. Previous studies have addressed invasion pathways of M. leidyi from its native region to Eurasia excepting the Mediterranean Sea [32][33]. These studies suggested that M. leidyi was introduced to Eurasia via at least two pathways. The first invasion occurred from the Gulf of Mexico to the Black Sea, followed by secondary spread to the Caspian Sea [32][33]. The second invasion was from the northern distribution of this species in the western Atlantic (possibly Narragansett Bay) to the Baltic and North Seas in northern Europe [32][33]. However, the source of the M. leidyi population in the Mediterranean Sea remains unclear. Several possibilities can be envisaged. It is possible the species has spread exclusively from the Black Sea [34] or other south Eurasian Seas in currents or in discharged ballast water. Alternatively, the species may have spread in discharged ballast water that originated in the North or Baltic seas, from the western Atlantic Ocean, or via a combination of the above pathways. To clarify the invasion pathway(s) of this species into the Mediterranean Sea, here we explore the population genetic structure of native and introduced populations using both mitochondrial (Cytochrome c Oxidase subunit I; COI) and nuclear ribosomal (Internal Transcribed Spacer; ITS) genes.

Ethics Statement
No specific permits were required for the described field studies in Eurasia, North America or South America. The species collected is an invasive pest in Eurasia and is not protected throughout its range. Sampling points did not include any protected or private lands.

Sample Collection and DNA Extraction
A total of 286 M. leidyi individuals were sampled from five native (Narragansett Bay, Massachusetts; York River, Virginia; Morehead, North Carolina; Tampa Bay, Florida; Peninsula Valdes coast, Argentina) and 11 introduced populations (two from the eastern Black Sea; Sea of Azov; north and south Caspian Sea; Baltic Sea; Limfjorden Fjord, Denmark; and Spain, France, Italy and Israel in the Mediterranean Sea). Individuals were preserved separately in 95% ethanol prior to genetic analysis.
Genomic DNA was isolated from gelatinous lobe tissue of the ctenophores using the automatic extraction protocol described by Elphinstone et al. [35], and DNeasy Blood and Tissue Kit (Qiagen Inc., ON, Canada). A fragment of the COI gene was amplified using the species-specific primers (Ml-COIF: 59-TGTCGCCCAAATTACTGTTTC-39 and Ml-COIR: 59-TGACGGGGTAAACCTCATAAA-39). Primers were designed in this study according to the available sequenced M. leidyi mitochondrial genome (GenBank accession no: NC016117). The universal primer pair, (ITS5F and ITS4R) [36] was used to amplify the ITS-1, 5.8 S gene, and ITS-2. We conducted PCR amplifications in a 40-mL reaction volume, with about 50 ng of genomic DNA, 1 unit of Taq DNA Polymerase (QIAGEN), 1 x PCR buffer, 2.5 mM of MgCl 2 , 0.2 mM of dNTPs, and 0.4 mM of each primer. PCR was performed with an initial denaturing step at 95uC for 1 min, followed by 35 amplification cycles (95uC for 30 s, 50uC for 30 s, 72uC for 50 s), and a final elongation step at 72uC for 7 min.

Sequencing and Cloning Protocol
We purified PCR products, which were then sequenced for both COI and ITS markers with forward (Ml-COIF) and reverse primers (ITS4R), respectively, using Big Dye terminator sequencing chemistry with an ABI 3130XL genetic analyzer (Applied Biosystems). Sequences were inspected, manually edited, and aligned using Codon Code Aligner 2.0 (Codon Code Corporation, Dedham, MA). Sequence of alleles containing double nucleotide calls (overlapping peaks) were cloned using Cloning and Amplification Kit (pSMART GC HK, Lucigen) according to Ghabooli et al. [33].

MtDNA Analysis
We assessed diversity indices within populations, such as the number of haplotypes (n), haplotype diversity (h) and nucleotide diversity (p) [37] using DnaSP v5 [38]. We constructed phylogenetic relationships among haplotypes using the neighborjoining algorithm in MEGA version 4 [39]. We used a fragment of COI from a cydippid ctenophore, Pleurobrachia pileus (GenBank accession no JF760211) as an outgroup. We generated a parsimony network of haplotypes using TCS 1.0 [40].

Nuclear marker (ITS) Analysis
Using the protocol described above, we processed four new populations from Mediterranean Sea (Spain, France, Italy and Israel) as well as one more from the native range (MH from North Carolina) in addition to our previously published dataset which consisted of 190 individuals analyzed for ITS marker [33]. We measured genetic diversity within populations with number of alleles (N a ), observed (H o ), and calculated expected heterozygosity (H e ) using GENEPOP (online version http://genepop.curtin.edu. au) and Arlequin version 3.1 [41]. We used the Markov chain method to estimate the probability of significant deviation from Hardy-Weinberg equilibrium using GENEPOP. We determined genetic differentiation among populations from pairwise F ST using Arlequin.
To estimate the sufficiency of our sampling, we generated rarefaction curves using ECOSIM and 5000 random iterations [42] for both haplotypes and alleles found in native region, and the Black-Azov and Mediterranean Seas. We estimated Chao-1 diversity [43] using SPADE software version 3.1 [44], based on the number of rare haplotype/allele present in sampled populations.

Results
Analysis of a 656-bp fragment of COI obtained from 241 individuals resulted in 17 different haplotypes in surveyed populations (GenBank accession nos KF435105-KF435121). In total, we detected 29 variable sites (4.42%), 16 of which were specific to the divergent haplotype Ml01 from Peninsula Valdes, Argentina (2.43%). Ml03 and Ml09 were the most common haplotypes. We found haplotype Ml03 in all populations except in Peninsula Valdes, while Ml09 was not recovered from Peninsula Valdes, Limfjorden, or the Baltic Sea.
We found twelve different haplotypes in native populations, all of which were present in introduced populations except for Ml01 from Peninsula Valdes, the single private haplotype at this site. We detected a total of 16 haplotypes among the introduced populations. Black-Azov Sea populations contained 11 haplotypes, which was higher than in all other introduced regions. Mediterranean Sea populations contained eight haplotypes, while those from the Caspian and Baltic seas had four haplotypes each. Out of eight haplotypes observed in the Mediterranean Sea, only Ml11 was not recovered from native populations in North America. Six haplotypes including Ml11 were detected in Black-Azov Seas. Two haplotypes from Mediterranean Sea populations were not found in either the Black or Azov Sea, though they were present in the native region, mainly in Florida and Morehead ( Figure 1, Table 1).
The Black-Azov Seas shared six haplotypes with native populations, while the other five haplotypes from this region were either private for one population (Ml04, Ml05, Ml10, and Ml13) or shared with France in Mediterranean Sea (Ml11). All four haplotypes found in Caspian Sea populations were present in both the Black Sea and North America. The Baltic Sea and Limfjorden (Denmark) shared all of their haplotypes with the native region, mainly Narragansett Bay, and only one haplotype with other introduced populations ( Figure 1, Table 1).
The introduced population (BL) from Black Sea contained the highest number of haplotypes (n = 7) (Table 1). Among introduced populations, those from Limfjorden and the south Caspian Sea had the lowest number of haplotypes (n = 2 and 3, respectively). Native populations from Morehead and Peninsula Valdes exhibited the highest (n = 7) and lowest (n = 1) number of haplotypes, respectively (Table 1).
Mean COI haplotype diversity (h) and nucleotide diversity (p) in all introduced populations were 0.70460.059 and 0.002060.0003, respectively. Comparable values in Mediterranean Sea populations were nearly identical, 0.70260.008 and 0.001960.0001, respectively. Native populations exhibited higher values for each of these indices (h = 0.85060.043 and p = 0.003860.0007, respectively). We excluded the non-diverse individuals of Peninsula Valdes of South America from this calculation.
The reconstructed phylogenetic relationship for the mtDNA haplotypes supported three main groups. The first group consists of the unique and highly divergent Ml01 haplotype restricted to South America, whereas the second one includes haplotypes Ml02, Ml07 and Ml08, which were common in northern areas of the distribution range in North America and Europe (Narragansett  Bay, Baltic Sea and Limfjorden). The rest of the haplotypes formed the third group ( Figure 2A). The complex parsimony haplotype network was star-shaped for the third group, with Ml03 in the middle. There were one or a few mutation steps between haplotypes, except for Ml01, which was separated from Ml03 by 19 mutation steps ( Figure 2B).
Chao-1 COI haplotype richness estimates were moderately higher than obtained values in Black-Azov Sea populations (15.2 vs. 11, respectively), indicating undersampling of these regions, although the lower 95% confidence interval limit (11.7) was marginally higher than observed diversity in the Black-Azov Seas ( Figure 3A). Chao-1 estimates for native region were also higher than the observed diversity (16 vs. 12), with the lower 95% confidence interval limit of 12.6 suggesting moderate undersampling of native region ( Figure 3A). However, Chao-1 estimates for the Mediterranean Sea were similar to the observed diversity (8.1 vs. 8), with the lower 95% confidence interval limit of 8 suggesting sampling was sufficient ( Figure 3A). The percentage of singletons for the Chao analyses of the native region, Black-Azov and Mediterranean Seas was 33, 45, and 13, respectively.
Analysis of the 619 bp DNA fragment comprising the complete ITS1, 5.8 S rRNA and ITS2 regions obtained from 286 individuals of M. leidyi -including the 190 individuals analyzed in our previous study [33] -resulted in 18 different alleles. We found five new alleles (GenBank accession nos KF435100-KF435104) in the Mediterranean Sea and Morehead (Figure 1) which were not previously identified. Alleles N and O were the most and least common, respectively. Alleles A and B were the most common in all populations (Figure 1), consistent with the previous survey of Ghabooli et al. [33].
We detected thirteen different alleles in native populations, all of which were recovered from introduced populations, except for the private allele G from Peninsula Valdes (Figure 1). Mediterranean Sea populations had 10 alleles, eight of which were present in native region. There was one private allele (O) in Haifa, Israel (Figure 1). Only five of 10 alleles found in Mediterranean populations were shared with Black-Azov Sea populations. In total, we recovered seven alleles in the Black-Azov Seas, six of which were also obtained from North America. Alleles C and J in Baltic Sea, Limfjorden, and Narragansett Bay were not present in Mediterranean populations, consistent with Black, Azov and Caspian Seas ( Figure 1, Table 1, see Ghabooli et al. [33]). The Chao-1 allele richness estimate for the Black-Azov Seas (Chao-1 estimator = 8; lower 95% confidence interval = 7.1; allele richness = 7) indicates reasonably comprehensive sampling of this region ( Figure 3B). For the native region, the estimated Chao-1 allele richness was 19.3, while the observed richness was 13, indicating undersampling of this area ( Figure 3B). For Mediterranean Sea populations, the Chao-1 estimates were similar to the observed diversity (10.2 vs. 10) with the lower 95% confidence interval of 10 indicating sufficient sampling in this region ( Figure 3B). The percentage of singletons for the native region, and Black-Azov and Mediterranean Seas was 38, 14, and 10, respectively.

Discussion
In this study, we build upon our previous study to explore genetic diversity, and determine the source(s) of, Mnemiopsis leidyi populations in the Mediterranean Sea using both mitochondrial (COI) and nuclear (ITS) markers. Our results support a multiple source model, composed by at least two different introduction pathways. One source of M. leidyi in the Mediterranean appears to have originated from Black Sea, consistent with the view of Bolte et al. [34] and with natural flows between the basins. However, we propose a second possible invasion pathway, originating from North America (Gulf of Mexico).

Genetic diversity and population differentiation
Introduced populations in the Mediterranean Sea exhibited lower values of haplotype diversity (0.70260.008) and observed heterozygosity (0.5860.05) relative to native ones (0.85060.043 and 0.6260.19, respectively). However, none of the Mediterranean populations exhibited erosion of genetic diversity for either of the analyzed markers relative to their putative source populations. This pattern could be driven by repeated introductions from the native range as well as from the adjacent Black Sea area, given intense vector activity between these regions and the high diversity of source populations [17], [45].
Two Mediterranean Sea populations (Spain and France) exhibited deviation from Hardy-Weinberg equilibrium (Table 1). Both populations exhibited lower than expected heterozygosity, which can be explained by possible inbreeding and/or population admixture (i.e. Wahlund effect) [33], [46]. We did not detect a heterozygosity deficit in other newly analyzed populations in Morehead and Limfjorden (Table 1).
Mediterranean populations had the lowest F ST with populations from the Black-Azov Seas (F ST = 0.001-0.028). However, Mediterranean populations also exhibited low genetic differentiation with those from Florida and Morehead (F ST = 0.010-0.032) in the native range. The highest genetic differentiation occurred among introduced populations in Mediterranean or Black-Azov-Caspian seas and those in the Baltic Sea and Limfjorden (Denmark), ranging from 0.177 to 0.417 (Table 2). High genetic divergence between introduced populations implies very low or lack of genetic connectivity and gene flow among these locations, implying that northern populations were not responsible for invasion of the Mediterranean Sea. As well, initial reports of invasion of the Mediterranean Sea occurred prior to those from the Baltic or North Seas [25][26].
Our results suggest that the Black-Azov Seas are a likely source of M. leidyi in the Mediterranean Sea, in accordance with Bolte et al. [34]. It is important to note, however, that the presence of similar alleles and haplotypes in the Mediterranean Sea and native populations -specifically those in the Gulf of Mexico and North Carolina -suggest a possible invasion pathway from North America. Namely, two COI haplotypes (Ml16 and Ml17) found in Mediterranean populations were not recovered from Black or Caspian Seas, but were present in native populations ( Figure 4A) in North America (Florida and Morehead). Similarly, our ITS survey revealed five new alleles for this species which were not recovered from populations in Sea of Azov, Black or Caspian Seas ( Figure 4B). Although the absence of the above alleles/haplotypes in Black and Caspian Seas populations may be explained by insufficient sampling from these regions or by seasonal, variation in frequency of genotypes/haplotypes, or other ecological and Phylogeography of Mnemiopsis leidyi PLOS ONE | www.plosone.org evolutionary processes, the possibility of introduction of M. leidyi from the native source region cannot be excluded. This conclusion is supported by our Chao-1 diversity estimates and rarefaction curves for native and Black-Azov Seas populations. These analyses indicate that our sampling recovered most of the diversity present in native and especially in the Black-Azov Seas and, therefore, the Black Sea as a sole source seems less likely.
The Mediterranean Sea receives an enormous flow of global shipping [16][17]. The tropical Western Atlantic Ocean is a source of trade to the Mediterranean Sea, and places it at risk of future invasions from discharged ballast water [17]. Moreover, high shipping activity within the Mediterranean Sea itself poses additional risk of translocation of M. leidyi and other NIS throughout the basin [17].
Despite of M. leidyi's dynamic invasion history, we observed geographic structure with some haplotypes/alleles being restricted to particular latitudes. The geographic distribution of genetic diversity is clearly not random and appears to reflect adaptation to specific biogeographic conditions. It is likely that this association is not only due to vector directionality but also to ecological and evolutionary processes [47]. The three haplotypes forming the second group in the NJ tree are very common in the northern region and less prevalent elsewhere (Figure 2A). The rest of the haplotypes that form the star in the parsimony network are distributed mainly in warmer waters and some were not found at all in northern regions ( Figure 2B). Shifts in haplotype/allele frequencies are expected due to selection to local conditions. Some haplotypes/alleles could become dominant in several generations if they are strongly favored by selection or linked to regions favored by selection [48][49]. Genetic differentiation among native populations was relatively high (F ST = 0.036-0.437), suggesting some structuring and limited gene flow in the native region. The private haplotype Ml01 was separated from other haplotypes by at least 15 mutation steps. All other haplotypes had only one or a few mutation steps between them. Pairwise genetic differentiation, parsimony network analysis, and phylogenetic reconstruction of haplotypes demonstrate high genetic divergence between South America and all other locations, notwithstanding the paucity of samples available from the former region. Long-term isolation of populations could explain the observed divergence. Pleistocene glacial periods in the northern hemisphere could drive high genetic divergence between South America and North America, resulting in population fragmentation [50]. However, further studies and more comprehensive sampling of the region, especially South and Central America, could shed light on the degree of isolation between populations along the western Atlantic coast. Our present, albeit very limited analysis does not support an introduction pathway for M. leidyi between South America and Eurasia.

Introduction pathways
Genetic analyses have revealed pathways of M. leidyi introduction into major Eurasia Seas [32][33][34]. M. leidyi entered the Black Sea via ships' ballast water from the Gulf of Mexico region. Spread of M. leidyi into the Sea of Azov occurred via the natural connection between these basins [26]. Secondary introduction into Caspian Sea likely occurred through ballast water discharged by a vessel after transiting the Volga-Don canal [51]. A second pathway from a port in New England, possibly Narragansett Bay, was likely responsible for the translocation of M. leidyi into the Baltic Sea, with subsequent spread into the North Sea [32]. The Mediterranean Sea was the most recent European basin invaded, with the eastern portion of the basin colonized first. Water flow between the Black and Mediterranean seas could account for this invasion, with subsequent transfer within the latter accommodated by a combination of ballast transfer and natural spread. Bolte et al. [34] used six microsatellite data to suggest a Black Sea source of M. leidyi in the Mediterranean Sea. However, in this study, genetic differentiation of North American and Mediterranean Sea populations was only slightly greater than that with populations from the Black Sea (Table 2). In addition, there were more haplotypes/alleles present in the Mediterranean Sea that were not shared by Black-Azov populations than with those from the Gulf of Mexico region (Table 1). Finally, there exists substantial ballast water movement from the Gulf region to the Mediterranean Sea [17]. Each of these lines of evidence supports the view that North America could have been an additional source of the introduced population in the Mediterranean Sea. The analysis of ITS and COI data in this study are consistent with the hypothesis of multiple introductions, with both native and Black Sea populations serving as sources of M. leidyi in the Mediterranean Sea.