Genetic Data Suggest Multiple Introductions of the Lionfish (Pterois miles) into the Mediterranean Sea

Widespread reports over the last six years confirm the establishment of lionfish (Pterois miles) populations in the eastern Mediterranean. Accumulated knowledge on lionfish invasions in the western Atlantic Ocean has shown that it is a successful invader and can have negative impacts on native species, indirect ecological repercussions and economic effects on local human societies. Here we analysed genetic sequences of lionfish from Cyprus as well as data from the whole distribution of the species, targeting the mtDNA markers cytochrome c oxidase subunit 1 (COI) and the control region (CR). Our results reflect a pattern of repeated introductions into the Mediterranean from the northern Red Sea and a secondary spread of this species west to Rhodes and Sicily. Presented results agree with previously published studies highlighting the genetic similarity with individuals from the northern Red Sea. Nevertheless, some individuals from Cyprus, in addition to those coming via the Suez Canal, were genetically similar to fish from the Indian Ocean, indicating genetic homogeneity among populations of P. miles across its current distribution, possibly facilitated by the ornamental fish trade and/or transport through ballast water.


Introduction
Invasive species are one of the main threats to biodiversity and natural resources, as they can cause severe changes in marine ecosystems [1][2][3]. The negative impacts often include alterations in the structure of marine communities and ecosystem services, with significant economic cost and even impacts on human health [4][5][6]. Monitoring, understanding and predicting the impacts of marine invasions has attracted much scientific interest [7,8]. The Mediterranean Sea is a biodiversity hotspot, as it currently hosts circa 17,000 known species, including more than 600 established alien spp., although its volume and surface are less than 1% of the world's oceans [3,9,10].
Climate change is considered to be one of the main factors that facilitated the establishment of thermophilic species that manage to reach the Mediterranean through the Suez Canal [11,12]. The Suez Canal is the main source of origin of alien species in the Mediterranean Sea; 13 out of the first 14 species that invaded at the beginning of the past century and 64% of the total number of species currently found in the Mediterranean are Lessepsian migrants (reached the Mediterranean Sea through the Suez canal from the Red Sea [9,13,14]). The lionfish, Pterois miles (Bennett, 1828) is one of the most successful marine invaders, exhibiting a set of behavioral and life-history traits that enables it to proliferate and establish alien populations [15,16]. The native distribution of P. miles includes the Indian Ocean, the Red Sea and some parts of Indonesia [17]. It has now spread throughout the Caribbean, the western Atlantic Ocean and recently into the Mediterranean Sea. More than 150 published studies in the last five years have focused on the western Atlantic region and documented impacts of the invasion of congeneric Pterois miles and P. volitans on local biodiversity, ecosystem services and the economy [3,16,18]. The presence of lionfish off the southeast coast of the USA was recorded for the first time in 1985 [19]. Less than 20 years later lionfish abundance was comparable to that of very common native species (e.g., Mycteroperca phenax (Jordan and Swain, 1884 [20])). Currently, lionfish presence has expanded, reaching the southeast coast of Brazil [21]. Based on genetic data, the spread of lionfish in the Caribbean and in the Mediterranean has been attributed to human activities such as the ornamental fish trade and the opening/widening of the Suez Canal, respectively [22][23][24][25][26]. Founder effects, bottlenecks, genetic drift and rapid population expansions are the main factors that determined the genetic profile of P. volitans in the Atlantic Ocean [27]. The close evolutionary relation of the congeneric species P. volitans and P. miles that invaded the western Atlantic Ocean and the Mediterranean, respectively, is highlighted by their ability to hybridize [28]. We still need to evaluate the effects of the lionfish presence in the Mediterranean, including interactions with native taxa such as groupers and moray eels that are also present at its native distribution range [29,30].
P. miles was first recorded in the Mediterranean Sea in 1991, off Israel [31]. More than two decades later, in 2012, two individuals of the species were captured off the coast of Lebanon [32], while the records from Cyprus started in 2013 [29]. Invasion through the Suez Canal is indicated by genetic data [24]. Abundance patterns at the onset of the species establishment around Cyprus' coastline exhibited a gradient of decreasing density from the southeastern part of the island, where it is now very common, to the southwestern part, where it is less abundant [33,34].
A recent study by Bariche et al. (2017) [24] that included individuals from Cyprus found low genetic diversity and the authors asserted that a few founder individuals may have been responsible for the invasion. Stern et al. (2018) [25] suggested there may have been more than one invasion of lionfish into the Levantine Sea. It is still unknown under which conditions this invasion took place, including whether it is ongoing with steady gene flow from the Red Sea to the Mediterranean. For this purpose, this study incorporates, for the first time, data from the whole range of the species' distribution.
Here we investigated the genetic diversity patterns of individuals sampled from Cyprus, enriching our dataset with sequences from other regions of the species' entire distribution in order to elucidate the invasion and establishment history of lionfish in the Mediterranean.

Sampling-DNA Extraction
The examined material was collected in the framework of the research program RELIONMED-LIFE (preventing a lionfish invasion in the mediterranean through early response and targeted removal) during scuba diving expeditions organized by the Marine and Environmental Research (MER) Lab Ltd. and the Enalia Physis Environmental Research Centre between September and October 2018 along the Cypriot coast. Immediately after collection, dorsal fins were taken from the samples and were placed in separate tubes containing 96% alcohol. All samples were stored at −20 • C until further laboratory analyses.
The great majority of the samples were collected from the southeastern coasts of the island where the population density was highest, but individuals from sites along the southern coasts of the island were also included in our analyses. In total, 56 individuals from 4 sampling sites were collected from Cyprus (Table S1). Small pieces of fins were used for total genomic DNA extraction with DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions. The final concentration and purity (A260/A280nm absorption rate) of DNA extractions were determined with NanoDrop 2000/200c (Thermo Fisher Scientific Inc., Waltham, MA, USA). In all cases the final DNA concentration was >100 ng/µL and the purity ratio >1.7.

Amplification and Sequencing
Cytochrome c oxidase subunit 1 (COI) gene and the control region (CR) were targeted following common PCR procedures for the amplification of the two regions using the primer pairs FISH-BCL/FISHCOIHBC [35,36] and universal CR-A/CR-E [37]. These mitochondrial loci were selected taking into account: (1) The availability of data already deposited in NCBI Genbank with which we can enrich our dataset and identify genetic similarities with populations/sources of introduction and (2) that the higher mutation rate of mtDNA compared to nDNA is more likely to provide useful information for tracking the invasion history of the species in the Mediterranean. Targeted loci were amplified in a Veriti thermocycler (Applied Biosystems, Foster City, CA, USA) with an initial denaturation step of 10 min at 94 • C, followed by 39 cycles of 1 min at 94 • C, 1 min at 52 • C and 1 min at 72 • C, with a final extension of 72 • C that lasted for 10 min. The final reaction volume in all cases was 20 µL, 0.1 of which was Kapa Taq DNA polymerase (5 U µL −1 ), 0.2 µL Kapa Taq DNA polymerase (5 U µL −1 ), 1.2 µL 25 mM MgCl 2 , 2 µL Kapa PCR buffer A, 0.6 µL 10mM dNTPs (Kapa, Sigma-Aldrich, St. Louis, MO, USA) 0.6 µL of each primer (10 µM) and > 20 ng DNA template. PCR products were purified with Qiaquick Purification Kit (Qiagen) following the manufacturer's protocol. Purified products were sent to Macrogen (Amsterdam, the Netherlands) for sequencing on both strands.

Data Elaboration and Phylogenetic Analyses
Sequencing results were reviewed and low quality ends were trimmed. The final dataset was enriched with available NCBI Genbank mtDNA sequences (Table S1). Available sequences of P. miles individuals from Genbank, representing the whole range of the species' distribution, were also included ( Figure 1). Additional sequences of the closely related species Pterois volitans, P. russelii and P. lunulata were also taken from the same database and included in the dataset as outgroups. Our final dataset included molecular data from 129 individuals (Table S1). which is the most common haplotype found in Cyprus. The common "Mediterranean haplotype" is also found in the Red Sea as well as in very distant places such as Indonesia and South Africa (Figures 3 and 4). Unique haplotypes were found in Cyprus and the diversity of haplotypes from the Red Sea and the Indian Ocean, where the lionfish is native, was high.  Multiple sequence alignments were performed with MAFFT v.7 [38]. Genetic divergence, as p-distance, which is the proportion of nucleotide sites at which two sequences being compared are different, within and between geographical groups of specimens and species were calculated using MEGA v.6.0 [39]. Likelihood scores for the selection of the best-fit nucleotide substitution model were calculated with jModelTest v.2.1.1 [40] using the following settings: Three substitution schemes, base frequencies estimation, gamma shape (four categories) and invariable sites estimation. Models including both gamma distribution and invariable sites were neglected [41]. The combined information from both sequenced loci were fed to MrBayes v.3.2.6 [42] and RAxML v.8.1.21 [43], where Bayesian inference (BI) and maximum likelihood (ML) phylogenetic analyses were implemented. In this way we tried to ascertain any possible phylogenetic signal provided by our data.
BI analysis was conducted with four independent runs and eight chains per run for 10 7 generations with a sampling frequency of 100. Hence, the analysis was based on 10 5 sampled trees from each run. Convergence among runs was monitored in MrBayes through the average standard deviation of split frequencies and by inspection of generation versus log probability of the data plot viewed in TRACER v.1.5.0 [44]. Finally, the 50% majority rule tree was constructed relying on 75,004 trees, since the 25% of sampled trees were discarded in the burn-in phase.

Geographic Population Structure-Haplotype Networks Construction
Analyses of molecular variance (AMOVA) and F ST values were estimated using Arlequin v3.5.2.2 [45]. Prior to analysis, available sequences were assigned to different local groups according to their geographic origin. Furthermore, three regional groups were formed in order to make comparisons between larger geographic regions. In particular, the Mediterranean group included individuals from Cyprus, Lebanon, Rhodes and Sicily; the Indian Ocean group included individuals from Madagascar, South Africa, Sri Lanka and Indonesia; and the Red Sea regional group included individuals from the local Red Sea and Gulf of Aqaba groups. Two different datasets representing each targeted genetic locus were created. Different haplotype networks for each gene were constructed using PopART v.1.7 [46]. The haplotypes detected by the software were distinguished according to their origin using traits block. Median-joining network haplotype networks [47], presenting pies proportional in size to the frequency of each haplotype and with different colors indicating the origin of each haplotype, were constructed.

Results
Sequences from 56 individuals of P. miles sampled along Cypriot coasts were obtained and deposited in NCBI Genbank (Accession codes MN150185 -MN150294). Possible contamination from Numts (mitochondrial sequences of nuclear origin) for both targeted genes was eliminated by the unambiguous alignment of retrieved sequences and the high similarity with already published GenBank sequences. Furthermore, in the case of the COI gene, no gaps or internal stop codons were detected after translation [48].
Individuals from very distant and geographically separate areas revealed extremely low genetic diversity among examined P. miles samples. The final alignment, including data retrieved from NCBI Genbank, of the concatenated dataset of both sequenced genes consisted of 968 bp coming from 129 individuals. More specifically, out of 612 bp examined for COI including/excluding outgroup, 517/593 were conserved, 41/19 were variable and 31/6 parsimony informative. For CR, out of 356 inspected sites, 96/105 were conserved, 261/251 were variable and 244/242 parsimony informative.
Calculated mean genetic distances for COI between geographic groups were extremely low, varying from 0.00% to 0.58% The range of variation for CR is 0.00% to 2.88%. For both genes within geographic group genetic divergence was in many cases higher than between groups. In the case of individuals from the Mediterranean region in particular, the highest genetic divergence was identified within the Cyprus populations, reaching up to 0.26 and 1.06 for COI and CR, respectively. The highest values of genetic distance for CR were found between USA and the Indian Ocean and for COI between South Africa and the Indian Ocean groups (Tables 1 and 2). The best substitution model selected for the concatenated dataset under Akaike Information Criterion (AIC) [49] was GTR + I with -In = 1131.2262 for COI and HKY + G with -In = 925.0418 for CR. Separate phylogenetic analyses for each gene were also conducted. These analyses resulted in gene trees exhibiting polytomies, unresolved phylogenetic relationships and low statistical support values. Hence, they are not presented herein. Nevertheless, concatenated, partitioned dataset analyses separated available P. miles sequences into two distinct clades ( Figure 2). Clade A is separated in two different sub-clades from which one is statistically supported. In subclade A(b), some individuals from Cyprus are grouped with samples from South Africa and the Indian Ocean ( Figure 2 and Figure S1). Nevertheless, the majority of samples from these regions formed a second unsupported sub-clade. Even though it is grouped with the rest of P. miles specimens included in our analyses, a single individual from the Indian Ocean forms a separate clade, Clade B, of which it is the only representative. Details about the representatives of clades and subclades are given in Table S1.
Statistically significant AMOVA and F ST values were calculated only in the case of CR gene. The great majority of variation (81.53%) was attributed to within population differentiation ( Table 3). The same results were also indicated by the relatively high genetic distances within geographic groups in both genes (Tables 1 and 2). Slightly different F ST values between the Mediterranean-Red Sea and Mediterranean-Indian Ocean regional groups support the genetic similarity of the invasive population in both regions (Table 4).
In total, six and three different haplotypes were identified for COI and CR genes, respectively, from samples collected along Cypriot coasts. Sicily and Rhodes appear to share the same haplotype, which is the most common haplotype found in Cyprus. The common "Mediterranean haplotype" is also found in the Red Sea as well as in very distant places such as Indonesia and South Africa (Figures 3 and 4). Unique haplotypes were found in Cyprus and the diversity of haplotypes from the Red Sea and the Indian Ocean, where the lionfish is native, was high.    Table S2. Hatch marks along edges represent the number of mutations between nodes. Accession numbers for individuals representing each haplotype are given in Table S2. Figure 3. Median-joining haplotype network for CR sequences. Different colors correspond to geographic origin and circle size is proportional to the number of individuals with the same haplotype. Hatch marks along edges represent the number of mutations between nodes. Accession numbers for individuals representing each haplotype are given in Table S2.  Table S3.  Table S3.

Discussion
The increasing number of reports of lionfish around the coasts of Cyprus show that the species is already established in the Mediterranean and is becoming more and more common [33,34]. While the ecological impacts of the invasion in the Atlantic Ocean have proved to be severe, we still do not know the magnitude of impacts in the Mediterranean [50]. In the case of the west Atlantic Ocean, lionfish are thought to have been released or escaped from marine aquaria [26,51,52]. While lionfish is not a common aquarium species in Cyprus, it is very common in neighbouring countries such as Lebanon (Bariche M, pers. comm.). As a result, the introduction of the species to the Mediterranean Sea through aquarium releases is also a possibility.
Our data for both examined genes are in line with the results of previous studies, as the same haplotypes are found in the Mediterranean and the neighboring Red Sea/Gulf of Aqaba. It is worth noticing that Italy, Rhodes and Lebanon share a common haplotype which in both genes is the most common in Cyprus. This fact strengthens the hypothesis that the first populations established in the eastern part of the Mediterranean were the source for the populations that are currently expanding westwards, reaching Rhodes and Sicily. On the other hand, the presence of six different, some of them unique, COI haplotypes in Cyprus may indicate ongoing gene flow between the northern Red Sea and the Levantine Basin. Beyond the presence of different haplotypes for COI and CR, a pattern of repeated introductions is also indicated by the fact that the species was reported for the first time in 1991 and then, after a long gap of 20 years, in 2012. The occasional presence of the species indicates its ability to reach the Mediterranean Sea and, at the same time, an initial failure to become established. Furthermore, the same pattern was documented in other invasive species of the same origin, such as the bluespotted cornetfish Fistularia commersonii [24,25,53]. In line with previously published genetic data from the Mediterranean Sea, the genetic similarity with individuals from the Red Sea is highlighted and also, at the same time, the low genetic diversity among Mediterranean populations [24,25]. Therefore, published studies support the idea that the lionfish, as many other species, invaded the Mediterranean Sea through the Suez Canal [8,11,54].
Nevertheless, as revealed by our phylogenetic analyses and haplotype networks, some of the individuals captured in Cyprus are genetically closer to individuals from the Indian Ocean and South Africa than to each other. This might be due either to underrepresentation of Red Sea populations or to a contribution of long-distance dispersal via the aquarium trade and/or ballast waters.
The presence of a dominant haplotype in the Cyprus population could be associated with a founder effect and rapid population expansion. However, the same level of genetic similarity with populations from the Red Sea and the broad geographic range of the Indian Ocean, as shown by the estimated F ST values, indicate the possibility of alternative pathways of the species' introduction in the Mediterranean Sea. Moreover, the attribution of total variance by AMOVA analysis to within-population differentiation (Table 3) along with the calculated genetic distance values (Tables 1 and 2), reflect the recent invasion of P. miles in the Mediterranean; hence, a limited time available for evolutionary divergence. Given the short period after the lionfish invasion, within the Mediterranean Sea population divergence is more likely to be the result of multiple introductions than of random genetic drift after a unique introduction event.
Dispersal to new areas is more effective when individuals are in the pelagic larval phase (14-17 days old) when they can disperse up to circa 900 km from the spawning area, taking advantage of the presence of local currents [16,55]. Nevertheless, if currents facilitated the transport of larvae from the Red Sea to the Mediterranean, then genetic divergence should correlate with geographical distance, as predicted by a 'stepping stone' model [56]. Our results do not show any correlation of geographic and genetic distances, as haplotypes from the Mediterranean are found also in remote regions, while the unique Mediterranean haplotypes have not been found in the Red Sea by this or previous studies [25]. Taking into account the possibility of introduction from very distant areas, P. volitans may also arrive and hybridize with P. miles in the Mediterranean [28]. This would alter the dynamics of the invasive population and hence the impact of the invasion on the native ecosystem. The introduction of new species into the Mediterranean through shipping ballast water is well documented, as circa 22% of all alien species currently present were introduced through this pathway [9]. The risk of lionfish larval transfer from the western Atlantic Ocean to the eastern Pacific Ocean through the Panama Canal was highlighted by MacIsaac et al. [57]. The present study verifies the genetic similarity of the Mediterranean lionfish populations with individuals from the Red Sea and at the same time points out the possibility of introduction from very distant areas following various pathways. These findings urge the competent authorities towards the direction of targeted measures. Biosecurity measures focusing on ornamental fish trade, ballast water treatment and the Suez Canal, which is one of the main paths that invasive species follow to reach the Mediterranean basin, are a clear priority to prevent Diversity 2019, 11, 149 9 of 12 continued introduction of lionfish and other invasive species. The latter is becoming an urgent need after the recent widening and deepening of the canal since the effects of these actions have not been evaluated yet.
Overall, no geographical pattern of genetic divergence among the lionfish populations studied was detected. It seems that the species is taking advantage of a variety of human activities and is expanding its distribution in the Mediterranean basin; this may have serious impacts on native species, ecosystems and the economy.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/11/9/149/s1, Figure S1: BI phylogram based on concatenated data set including both targeted genes (CR and COI), Table S1: Species, locality of origin, available sequence data from targeted genes, and codes/ Genbank accession numbers of individuals used in the analyses., Table S2: Frequency and Accession numbers of CR haplotypes, Table S3: Frequency and Accession numbers of COI haplotypes.