Clarifying the taxonomic status of the alien species Branchiomma bairdi and Branchiomma boholense (Annelida: Sabellidae) using molecular and morphological evidence

This study was performed to analyse the genetic and morphological diversity of the sabellid annelid genus Branchiomma, with special emphasis on a taxon so far identified as Branchiomma bairdi. This species, originally described from Bermuda, has frequently been reported as an invader in the Mediterranean, the Atlantic and the Eastern Pacific, but recent observations have raised some taxonomic questions. Samples of this taxon were collected from five sites in the Mediterranean Sea, two sites in the original distribution area of B. bairdi in the Gulf of Mexico and four localities in the east Pacific and Atlantic Oceans where B. bairdi has been reported as invasive. The molecular results revealed a conspicuous genetic divergence (18.5% K2P) between the sampled Mediterranean populations and all the other ones that led to a re-evaluation of their morphological characters. The latter showed that the Mediterranean and extra-Mediterranean populations also differ in some discrete morphological and reproductive features. Consequently, the Mediterranean samples were re-designated as B. boholense, another non-indigenous species originally described from Philippines. Branchiomma bairdi and B. boholense differ in body size, development and shape of micro and macrostylodes, size of radiolar eyes and body pigmentation. Genetic diversity was high in B. boholense from the Mediterranean as well as in B. bairdi from the Gulf of Mexico, but low in B. bairdi populations outside their native range. The phylogenetic analysis revealed the presence of connections between the Mediterranean localities as well as between native and introduced B. bairdi populations that focus the attention on the Panama Canal as important passage for the introduction of the species from the Gulf of Mexico to the north-east Pacific Ocean.


Introduction
In this regard, in the last decades, molecular tools have been proven particularly useful to improve our understanding of the causes and consequences of biological invasion [50][51][52][53][54][55]. Genetic methods offer the possibility to gain insights into many aspects of the invasion process, including a better understanding of sources and routes of invasions, the detection of cryptic diversity and the connectivity among native and introduced populations [56][57][58][59][60][61][62]. However, despite all the benefits provided by the molecular approach, one of the best ways to fully characterize aquatic invasions is to combine genetic and morphological approaches, as the availability of taxonomic expertise is considered the first crucial issue for dealing with the assessment and management of NIS [11,54,[63][64][65][66][67].
This synergistic approach might be particularly useful when the species show high intraspecific phenotypic plasticity or when distinct evolutionary lineages do not show clear morphological diversity as it has been observed for the species of the genus Branchiomma [18].
The present study was first designed to investigate pathways of introduction of B. bairdi in the Mediterranean Sea as well as in other extra-Mediterranean invaded regions. To this end, specimens of the taxon so far identified as B. bairdi were sampled in different Mediterranean and extra-Mediterranean sites, including some localities in the original distribution range of the species and examined via amplification of the mitochondrial cytochrome c oxidase subunit I (COI). However, the molecular results suggested the presence of two different taxa between the sampled Mediterranean and extra-Mediterranean specimens. Therefore, the purpose of this study was reoriented towards a clarification of the taxonomic status of the collected Branchiomma species and, consequently, results from molecular analysis were integrated with a thorough morphological examination that was performed on several individuals of each studied population.  Table 1). For Mexican material (Veracruz and Mazatlán) the sampling permission was granted by the Comisión Nacional de Acuacultura y Pesca (Authority: Lic. Aldo Gerardo Padilla Pestaño). Tampa Bay material was collected during the fouling project specifically designed to search for non-indigenous species [17,20,68]. The permission was issued by the Smithsonian Environmental Research Center (SERC) in Edgewater, Maryland. This material was part of the SERC reference collection (Authority: Dr. Gregory Ruiz). The permission for samples collection in Galapagos was allowed by the Ministry of Environment of Ecuador (MAE) and by the Galapagos National Park Directorate, while the permit for sampling specimens from Hawaii was granted by Hawaii Department of Land and Natural Resources. Sampling in Madeira island was allowed by the administrator of the Funchal Marina without any specific permission. Specimens collected from Madeira Island in our study are not considered an endangered or protected species. Taranto samples were collected in private land (aquaculture facilities Cantiere Navale Greco) and the permission was granted by the owner Dr. Francesco Grego. Samples  Table 1.

Sample collection and preparation
https://doi.org/10.1371/journal.pone.0197104.g001 Collected worms used for the molecular analysis were preserved in absolute ethanol and stored in the laboratory at -20˚C until processing. Some specimens were fixed in 4% buffered formaldehyde seawater solution and then preserved in 70% ethanol for morphological analysis.
In addition, seven specimens of B. boholense [22] from Indonesia (Komodo Island, Sumba and Banda Sea) were morphologically examined and measured for comparative purposes. This material was collected on sandy calcareous bottom from 10 to 90 m depth in 1984 during the Indonesian-Dutch Snellius II Expedition (Zoological Museum of Amsterdam).

Molecular analysis
DNA extraction was carried out following the procedure used by [18] using tissue removed from the posterior end of the worms. Genomic DNA was extracted using standard protocols for the DNeasy Blood and Tissues Kit (QIAGEN Pty Ltd, Dusseldorf, Germany) and stored at -20˚C until further processing. The primer cocktail designed by [69] was used for the amplification of the COI gene for all the individuals of the Mediterranean populations (MM, TA, LA, SA, CA) and for some specimens collected from the Mazatlán populations (MZ). The amplifications were carried out in a 25 μl reaction containing 1 μl of the genomic DNA, 0.5 μl of Hot-StartIt Taq DNA polymerase, 1.5 μl of MgCl 2 (25 mM), 1 μl of dNTPs (50 mM each), 5.5 μl of ddH 2 O, 12.5 μl of 10% trehalose, 2.5 μl of 10X buffer HotStartIt and 0.25 μl of each primer (10 mM). The PCR thermal cycling profile was 94˚C for 1 min, 5 cycles of 94˚C for 30s, 45˚C for 40s, 72˚C for 1min followed by 35 amplifications cycles of 94˚C for 30s, 50˚C for 40s, 72˚C for 1min with a final extension step at 72˚C for 10 min.
As the PCR using the COI primer cocktail [69] was unsuccessful for the amplification of the populations of Madeira, Veracruz, Galapagos, Tampa Bay and Hawaii, a new set of primers PolyBr_COIF/PolyBr_COIR was designed: the forward primer 5' TCWATWAGWGT WATTATTCGKGCTG 3' and the reverse 5' CMGCAGGATYAAARAACCTAGTA 3'. PCR mixture (20 μl in total volume) contained 0.5 of the genomic DNA, 0.08 μl f Platinum Taq DNA polymerase, 0.6 μl of MgCl 2 (50 mM), 0.4 μl dNTPs (50 mM), 15.42 μl of ddH2O, 2 μl of 10X buffer and 0.5 μl of each primer (10 mM PCR products were verified by means of electrophoresis on 1% agarose gel with ethidium bromide. The products of successful amplifications were purified using the ExoSAP-IT PCR purification system (USB Corporation, Cleveland, Ohio, USA) and then bidirectionally sequenced using the forward and reverse primer M13 for the Mediterranean and Mazatlán populations. For the sequencing of the samples collected from all the other populations, the new set of primer was utilized. Cycle sequencing was performed using BigDye Terminator version 3.1 in a 10 μl reaction containing 5 μl of 10% trehalose, 2 μl of 5X BigDye buffer, 1 μl of BigDye, 1 μl of each primer (or primer cocktail) and 1 μl of the PCR products. After cleanup of the sequence reactions using Zymo Sequencing Clean-Up Kit, the sequences were analysed with an ABI PRISM 1 3130 Genetic Analyzer [http://dx.doi.org/10.17504/protocols.io. nx6dfre]. A total of 123 sequences were generated using either of the two primer combinations (S4 Table).
Forward and reverse sequences were combined and manually edited with Sequencher v. 4.8 (Gene Codes Corporation, Ann Arbor, MI, USA). The sequences were aligned using the algorithm CLUSTAL W implemented in MEGA 5 [70] and trimmed to a uniform length of 493 bp.
The levels of genetic polymorphism for each studied population were estimated by the number of haplotypes (Nh), haplotype diversity (h) and nucleotide diversity (π) using DNAsp 5.10.01 [71].

Morphological analysis
For this study, detailed observations of the external morphology of specimens collected from different localities were undertaken using light microscopy. For each locality, 15 adult individuals were examined and for all the specimens, the total body length (including radiolar crown, thorax and abdominal regions), the trunk length (only thorax and abdominal areas) and the radiolar crown lengths were measured. Moreover, some morphological features such as length of radiolar tips, shape of stylodes, dorsal lips and number of rows of teeth above the main fang of thoracic and abdominal uncini were analysed. The number of thoracic and abdominal chaetigers and the ratio of crown-trunk lengths were also measured.
All the analysed morphological features were photographed using a digital camera attached to the stereo-and compound microscopes (Nikon, Coolpix 990 and Carl Zeiss AxioCam Erc 5s, Zen 2012, Carl Zeiss Microscopy GmbH, Jena, Germany). Thoracic uncini were drawn with the aid of a camera lucida. Due to the variation in size and morphology of uncini along the body segments and within the tori, their morphology was analysed considering the dorsal most uncini of second thoracic chaetiger.
A one-way analysis of variance (ANOVA) was performed by using STATSOFT STATIS-TICA v. 6.0 to test for significant differences in the mean value of the crown/trunk length ratio between Mediterranean and extra-Mediterranean samples. Pearson correlation analysis was carried out on the number of chaetigers and the values of body length in order to evaluate the relationship between the two measurements taken.

Genetic analysis
The sequencing of the 493 base pair fragment of the COI gene in 123 individuals identified 11 haplotypes. These haplotypes derived from 44 polymorphic sites of which 43 were parsimony informative. Most polymorphisms were the results of synonymous changes (43 of the 44 sites) while five mutations produced amino acid replacements (non-synonymous changes). Seven of the 11 identified haplotypes were location private, six of which were represented by a single individual. Four haplotypes were shared by individuals collected from distinct populations and sampling sites, two of which resulted the most common haplotypes: haplotype 1 was present just in the Mediterranean Sea and was shared by 35 individuals belonging to all the sampled Mediterranean populations, while haplotype 9 was shared by 67 individuals collected from different populations in both the Pacific Ocean, the Gulf of Mexico and the Atlantic Ocean (S1 Table).
Haplotype (h) and nucleotide (π) diversity were quite different within populations ( Table 1). Haplotype diversity was high in most of the Mediterranean populations (mean h ± standard deviation = 0.806 ± 0.218), while samples from Madeira, Mazatlán, Galapagos and Hawaii populations showed low values of haplotype diversity ranging from 0.000± 0000 to 0.100±0.088. Similarly, the mean nucleotide diversity of the Mediterranean samples was high (0.0031 ± 0.001), while Mazatlán, Madeira and Hawaii populations showed minimum values of nucleotide diversity, ranging from 0.000 ± 0.000 to 0.00155 ± 0.00067, highlighting the higher genetic diversity observed in the Mediterranean samples. By contrast, the maximum values of nucleotide diversity were detected in Veracruz populations with h = 0.01328 ± 0.00159 ( Table 1).
Values of pairwise F ST and average K2P genetic distance analysis evidenced a great molecular divergence between the Mediterranean populations and those sampled both in the Gulf of Mexico, Atlantic Ocean and Pacific Ocean (S2 and S3 Tables). Indeed, all the comparison between the sampled Mediterranean populations and the Pacific, Mexican and Atlantic ones showed high and significant pairwise F ST values, ranging from 0.93957 to 0.99684. By contrast, low and non-significant genetic divergence was observed in all the comparisons among the Mediterranean populations, as well as between the Pacific populations and the Atlantic ones, including populations from the Gulf of Mexico. Moreover, no molecular differentiation was observed in all the comparisons between Madeira, Mazatlán, Hawaii and Tampa Bay populations, while slightly higher values were obtained for the comparisons between Veracruz population and Madeira, Mazatlán, Hawaii, Galapagos and Tampa Bay (S2 Table). Comparably, values of the average Kimura two-parameter distance calculated between the different populations, evidenced a strong genetic variation (18.5% ± 0.000614) between the Mediterranean populations compared to all the other ones, while comparisons between the Mediterranean populations as well those independently calculated within the Pacific, Mexican and Atlantic group of populations showed no significant average genetic distance (S3 Table).
The phylogenetic tree constructed using the Bayesian approach evidenced two distinct monophyletic groups of samples (posterior probability, pp. 1.00) (Fig 2). Group 1 was well supported (pp. 1.00) and enclosed samples belonging to populations collected in the Gulf of Mexico, Pacific Ocean and Atlantic Ocean. On the other hand, Group 2 comprised exclusively COI sequences obtained from individuals collected in the Mediterranean basin (pp. 1.00). The COI sequence identified as Branchiomma bairdi [78] and used as outgroup did not cluster inside the studied clades of samples, probably because it refers to another Branchiomma species.
The MJ network of haplotype showed results highly concordant with the phylogenetic analysis. The analysis distinguished two groups of haplotypes, separated by 33 mutations (Fig 3). Haplogroup 1 included three haplotypes (Hap_9, Hap_10, Hap_11), of which one (Hap_9) was shared by the 89% of the individuals sampled from different populations in both the Gulf of Mexico, Pacific and Atlantic areas. Haplotype 10 was present in seven individuals sampled from Veracruz, Galapagos and Tampa Bay populations, while haplotype 11 was represented by OTUs names consist of the abbreviated location name according to Table 1 and a unique identifier. For Mediterranean populations, dark green indicates samples collected from carbon dioxide vents, light green samples from non-acidified sites. For extra-Mediterranean populations, dark blue indicates samples from sites where the species was introduced, light blue indicates samples from the native range. The tree was rooted using sequences of Bispira manicata, Bispira serrata, Branchiomma sp. and Branchiomma bairdi as outgroups (taxon names followed by GenBank 1 accession a single individual and was location private. Haplogroup 2 contained eight haplotypes exclusively present in the Mediterranean basin. Haplotype 1 resulted the most frequent haplotype inside the basin, as it was shared by 35 individuals belonging to all the sampled Mediterranean populations. Haplotype 6 differentiated from haplotype 1 by one mutation and was present in four individuals, two sampled in the locality of Lacco Ameno and two belonging to Mar Menor population. Haplotype 2 was present just in the Taranto population while all the other haplotypes were location private and represented by a single individual (Fig 3, S1 Table).

Morphological analysis
The two groups of populations that molecular analysis showed to be genetically divergent, also showed slight morphological differences. As shown in Table 2, one of the main differences numbers). Only nodes supported by posterior probabilities ! 0.95 are reported. ÃÃ , posterior probability = 1; Ã , posterior probability ! 0.95 and < 1.
https://doi.org/10.1371/journal.pone.0197104.g002   Table 3). The values of crown/trunk length ratio resulted higher in Mediterranean populations than in extra-Mediterranean samples (Fig 4). Results from Pearson correlation analysis evidenced a linear relationship between body size of the worms and the number of chaetigers both in Mediterranean ( Fig  5A) and extra-Mediterranean samples (Fig 5B).
The morphological features of Mediterranean and extra-Mediterranean specimens are listed in Table 4. Individuals from Mediterranean material showed a brownish colouration with small dark spots both in preserved and live specimens (Fig 6A). Their radiolar crown showed to bear from 20 to 25 radioles with short radiolar tips (as long as equivalent space of 6-8 pinnules), in some specimens abruptly tapering with pinnules extending up to the radiolar tips ( Fig 7B). Dorsal lips were as long as 1/3 of the crown (Fig 7E). Microstylodes appeared  digitiform in all the examined individuals, small, as long as rachis width or slightly longer (Figs 7A, 7C and 8A, 8B) and not covering the radiolar eyes that appeared small (covering less than 1/2 of side of radiolar rachis at mid-radiolar length) along the radiole and brownish coloured (Figs 7A, 7C and 8A, 8B). Flattened, tongue-like macrostylodes, 2-3 times longer than microstylodes, were present in number of 1-4 pairs per radiole (Figs 7A, 7C and 8A, 8B). However, high variability in the morphology of macrostylodes was observed according to the position of the radiole within the radiolar crown, indeed, some macrostylodes appeared strap-like (cylindrical). Variation was also viewed in tips of macrostylodes, sometimes being distally rounded others being acuminate (Fig 8A and 8B). The basal microstylodes appeared small (as long as rachis width or slightly longer) and unpaired (Fig 6B). Thoracic and abdominal uncini with 1-2 rows of teeth over the main fang (Fig 7D and 7F). Specimens from all the extra-Mediterranean material showed a green (in alive specimens) and brownish colouration (in fixed specimens) with interramal small brownish spots along the body and an orange-red pigmented crown in some specimens (Fig 6A). Their radiolar crown showed from 15 to 18 pairs of radioles terminating with long radiolar tips, as long as about 10-15 pairs of pinnules (Fig 9B). Several specimens showed the presence of the eggs within the crown, suggesting brooding ( Fig 6C). Dorsal lips were as long as 1/3 of the crown (Fig 9E). The microstylodes appeared long (notoriously longer than rachis width), thin and digitiform in all the examined individuals, but not covering the radiolar eyes that appeared highly developed (covering 3/4 of side of radiolar rachis at mid-radiolar length) and orange coloured (Figs 8D and 9A, 9C). Macrostylodes were strap-like, 2-3 times longer than microstylodes and 1-4 pairs per radiole (Figs 8C, 8D and 9C). The basal microstylodes appeared unpaired (Fig 6C). Thoracic and abdominal uncini with 2-3 rows of teeth over the main fang (Fig 9D and 9F).
In Fig 10, the morphological features of B. boholense specimens from its original area of distribution (Indonesia) are shown. Microstylodes appeared digitiform and long as rachis width or slightly longer, while macrostylodes were 2-3 times longer than microstylodes and tonguelike (flattened) in shape (Fig 10B, 10D, 10E and 10F). The basal microstylodes appeared small (as long as rachis width or slightly longer) and unpaired (Fig 10A). Radioles showed short radiolar tips (as long as equivalent space of 6-8 pinnules) (10C) and radiolar eyes appeared small (covering less than 1/2 of side of radiolar rachis at mid-radiolar length) (Fig 10D, 10E and 10F). Thoracic and abdominal uncini showed two rows of teeth over the main fang (Fig 10G).

Taxonomic status of the Mediterranean taxon
This study revealed significant genetic differentiation in the mitochondrial gene COI between Mediterranean and extra-Mediterranean Branchiomma populations, strongly suggesting that they represent two different species. Observed values of pairwise F ST and average K2P genetic distance denoted that the two groups of populations belong to two different congeneric species, indeed the levels of genetic distance (K2P) observed between the two groups of populations were similar to those reported for other congeneric polychaete species [18, [79][80][81][82][83][84][85]. The results obtained from the phylogenetic analysis emphasized the strong separation between the two forms as well: sequences obtained from extra-Mediterranean samples and Mediterranean ones clustered in two distinct and well supported clades.

Morphological features Mediterranean taxon (presumably B. boholense) Extra-Mediterranean taxon (B. bairdi)
numbers of radioles 20-25 pairs 15-18 pairs basal stylodes small (as long as rachis width or slightly longer), unpaired long (notoriously longer than rachis width), unpaired shape of macrostylodes flattened, tongue-like mainly but strap-like may be present cylindrical, strap-like shape and size of microstylodes short (as long as rachis width or slightly longer), digitiform long (notoriously longer than rachis width), digitiform compound eyes small (covering less than 1/2 of side of radiolar rachis at midradiolar length) /brownish big (covering 3/4 of side of radiolar rachis at midradiolar length)/orange radiolar tips short (as long as equivalent space These results led to a careful re-examination of the morphology of the two groups of samples that were molecularly divergent. The morphological features observed for extra-Mediterranean specimens correspond to those reported for the species B. bairdi [86]. By contrast, the morphology of the Mediterranean taxon led to presumably reconsider the species B. boholense for our Mediterranean material originally identified as B. bairdi. Indeed, Mediterranean specimens showed morphological features comparable to those reported for specimens from Indonesia, reviewed in this paper (Fig 10) and identified as B. boholense [22] as well as to the syntype examined by Cepeda & Rodríguez-Flores [41]. Besides, the taxonomic characters observed for our Mediterranean material did not coincide with those considered diagnostic for the identification of any other Branchiomma taxa, including species native of the Mediterranean area. However, lacking molecular data for the material from Indonesia, it was not possible to corroborate our morphological results with the molecular data during this phase. This approach is highly suggested in the future to validate our findings. It is worth mentioning that biological data on reproduction and population dynamics concerning the Mediterranean form previously referred as B. bairdi [29,30] could probably be attributed to B. boholense. Significant biological differences found between the Mediterranean populations and Mazatlán (Gulf of California) population were previously explained as adaptation to a new environment [29,30]. Results from the present study led to reconsider these biological differences due to the fact that the two examined taxa are actually different species. Indeed, the two species differ in the reproductive features, with a different period and mode of reproduction. Brooding was reported for Malta and Mazatlán populations [26,45], while external fertilization occurs in Mediterranean specimens [30]. Moreover, no sign of asexual reproduction was found in the Mediterranean taxon contrary to what has been reported for B. bairdi [45]. During the morphological analysis, the presence of the eggs among radioles in the crown of B. bairdi were often detected, and never observed in the Mediterranean taxon presumably identified in this study as B. boholense.  Morphological distinctiveness. As already reported by different authors [17,22,25,26,37,41,42], Branchiomma bairdi and B. boholense are difficult to distinguish and it is likely that in the past the two species have been frequently misidentified. According to Tovar-Hernández et al. [43], the main distinctive characters separating these two species are the thoracic uncini dentition and the different morphology of macrostylodes: strap-like macrostylodes and 2-3 rows of teeth on thoracic uncini are present in B. bairdi, whereas B. boholense has flattened tongue-like macrostylodes and thoracic uncini with one large tooth. However, recently, Keppel et al.
[17] and Cepeda & Rodríguez-Flores [41] showed that also in B. boholense two rows of teeth can be present, suggesting that this character cannot be considered for the separation of the two taxa.
Our morphological observations on Mediterranean and extra-Mediterranean specimens as well as on individuals from Indonesia integrated the results reported by Cepeda & Rodríguez- Flores [41] that compared some Mediterranean material identified as B. bairdi with the syntype of B. boholense giving some useful suggestion about the main differences separating the two taxa.
Results from this study suggested that the main discrete morphological differences between the two taxa are the shape of macrostylodes, the development of microstylodes and radiolar eyes, coupled with differences in body size and in the number of teeth over the main fang of thoracic and abdominal uncini. Specimens from Mediterranean populations, presumably to be considered as B. boholense, showed mainly flattened, tongue-like macrostylodes, short microstylodes (as long as rachis width or slightly longer) and small brownish radiolar eyes (covering less than 1/2 of side of radiolar rachis). By contrast, specimens of B. bairdi had cylindrical, strap-like macrostylodes, long microstylodes (notoriously longer than rachis width) and big red-pigmented radiolar eyes (covering more than 3/4 of side of radiolar rachis). Moreover, both species showed more than one row of teeth over the main fang. Nevertheless, in B. bairdi the main fang appeared followed by 2-3 rows of teeth, while in the Mediterranean taxon as well as in specimens from Indonesia (Fig 10G) reviewed in this paper, only two rows of teeth over the main fang seemed to be present.
The slight difference observed between the morphology of the Mediterranean material and the syntype of B. boholense examined by Cepeda & Rodríguez-Flores [41], could be due to the fact that the syntype was a juvenile (small body size and short, blunt radiolar tips typical of juveniles or regenerating worms), and it was in very bad preserved conditions. In particular, the difference between the two taxa concerning the first, basal stylodes that were reported paired for B. boholense by Cepeda & Rodríguez-Flores [41], was probably a mistake due to bad preservation, because the distance between the base of the crown and their appearance suggests that they probably referred to the second pair of stylodes that obviously are paired. In addition, well preserved material of B. boholense from Indonesia also examined in this study (4 adult and 3 juvenile specimens), only present unpaired basal stylodes (Fig 10A). However, Capa et al. [18] already emphasized the presence of basal stylodes being paired in some few taxa (B. lucullanum, B. bombyx and B. galei) and showing plasticity in the majority of the other studied Branchiomma species. But the relation of presence-absence of this feature according to juvenile or adult stages, or regenerating worms have not been determined.
As this study shows, two invasive species of Branchiomma in the Mediterranean present different reproductive modes and discrete morphological differences, including body size. Capa et al. [18] already highlighted that the traditional diagnostic morphological features of species in the genus Branchiomma are greatly homoplastic and polymorphic, and Keppel et al. [17] emphasized that other characters including body size, radiolar crown, dorsal lips and thorax length are not informative, as they are impacted by anaesthetics and fixation methods, or by the regeneration process. Moreover, during the process of morphological identification other difficulties can arise as it is not always possible to work using live specimens or well fixed or preserved samples, and sometimes may not be possible to corroborate identifications using type materials. Six species of Branchiomma have been reported as alien worldwide, but it appears that many additional NIS remain undetected for this family, in multiple geographic regions around the world. Consequently, for early detection programs of invasive species we recommend the combined use of molecular and morphological approaches for Branchiomma species. It would work at least for B. bairdi, B. boholense, B. bombyx, B. luctuosum and other seven unnamed, but alien species included in Capa et al. [18]. However, a major revision of type materials is still needed in order to provide the available names for those alien species.
Species distribution within the Mediterranean basin. On the basis of both molecular and morphological analysis the present knowledge led to reconsider the distribution of Branchiomma bairdi and B. boholense within the Mediterranean basin.
Both species are present in the Mediterranean area, although, following the results of this study, the previous records of B. bairdi reported for the Italian coast should probably represent B. boholense [25,26,27,29,30,87]. Even though B. boholense seems to have been present in the Eastern basin since the early 1900s [22], recently the species has shown to spread and establish successfully in the central and western Mediterranean basin as shown by its recent reports along the Italian coast [25][26][27] and in the south-eastern coast of Spain [37]. By contrast, B. bairdi was actually collected in the islands of Malta and Formentera [26,41], and along the Turkish and Tunisian coasts [38,40].
The presence of B. bairdi in the western side of the basin and its records at Canary and Madeira Islands [26,42] is consistent with the probable introduction of B. bairdi from the Atlantic Ocean throughout the Gibraltar Strait. By contrast, considering the Indo-Pacific origin of B. boholense and that the first findings of the species are from the eastern side of the basin [22,34], the species might have entered the Mediterranean trough the Suez Canal.

Genetic structure
The mitochondrial genetic diversity observed within populations of the examined species are comparable to those reported for other polychaetes [88,89]. However, different genetic scenarios are displayed for the two introduced Branchiomma species, with the Mediterranean taxon (presumably B. boholense) showing high levels of genetic diversity, and, otherwise, the introduced Pacific and Atlantic B. bairdi populations exhibiting low levels of genetic diversity.
Different ranges of genetic diversity have been documented for invasive species according to their life history traits and to the different type of introductions [57,[90][91][92][93]. As proposed by Voisin et al. [94], invasive species typically present two distinct patterns of genetic diversity: a reduced genetic variability associated to stochastic introductions of a small number of individuals that are more likely subjected to successive bottleneck, or increased levels of genetic diversity due to multiple introductions.
The pattern of genetic diversity of (presumably) B. boholense observed in the Mediterranean basin is in accordance with the second hypothesis. Indeed, the expansion of this species across the Mediterranean is very recent with possible recurrent introduction events recorded through different years and in localities that are wide apart across the basin [25,26,34,37]. A similar pattern of high genetic variability associated to continuous introductions from different sources was observed in some invasive ascidian species such as Styela clava Herdman, 1881 [57] and Styela plicata (Lesueur, 1823) along the Italian coasts [92].
This pattern of increased within-population genetic diversity may be also associated to the life history traits of B. boholense, which showed a short life span, associated with a fast growth, early maturity and a rapid generation turnover ( [29] as B. bairdi). Moreover, this species displays a great ecological plasticity and adaptability in colonizing different habitats and substrates. In particular, among the examined sites around Ischia Island the species was collected also at the Castello Aragonese vent's system, a venting area where intense CO 2 emissions lower the seawater pH up to 6.6 values [95] and where the few polychaete species collected in this site respond to stress either by acclimating via phenotypic plasticity (e.g the sabellid Amphiglena mediterranea) or through genetic adaptation (e.g. the nereidid Platynereis spp) [85,96,97]. Branchiomma boholense seems to acclimate to the CO 2 vent conditions like A. mediterranea, as the vent-inhabiting populations are not genetically distinct from nearby non-acidified populations. The ecological plasticity of the species coupled with the high genetic diversity observed might in turn enhance its potential of invasion. A connection between high intrinsic genetic diversity and high ecological plasticity was observed also for introduced populations of the Mediterranean species Sabella spallanzanii (Gmelin, 1791) in Australian waters [98], where it is considered a pest species [66].
A different scenario of within-population genetic diversity was observed for the studied native and introduced populations of B. bairdi. As expected, populations sampled from the original distribution range of the species (Tampa Bay and Veracruz) showed very high values of diversity mainly linked to their higher size and stability compared to the introduced ones. Indeed, these populations belong to a taxon with a huge Caribbean distribution, described for Bermuda in 1885, but successively found in Florida [99] and in Veracruz [100][101][102].
By contrast, the low values of diversity observed in the studied introduced B. bairdi populations (Madeira, Mazatlán, Galapagos, Hawaii) and their monomorphic condition denoted the occurrence of recent bottlenecks from a single introduction event that accounted for the total loss of the native genetic variability observed. This is also consistent with the recent age of introduction of B. bairdi in the southern Gulf of California [43] as well as in Galapagos, Hawaii [20,68] and Madeira Island [42]. The recent introduction of this species within Madeira fauna can be inferred considering that Madeira is a site whose fauna is well-known and highly investigated over time [103].
Moreover, Madeira, Galapagos and Hawaii are island populations and, as already observed in other studies [104,105], such populations show lower levels of genetic variation than their mainland counterparts as stochastic introduction events counteract the arrival of new alleles in the gene pool. However, the pattern of diversity observed in the introduced populations of B. bairdi can be due to the presence of asexual reproduction [45] as well, with the propagation of clones that could have genetically homogenized the population [106,107]. Nevertheless, this low genetic diversity does not seem to affect the successful invasion of the species as already observed for other invasive taxa [57,108].
Pair-wise F ST and K2P analysis evidenced a genetic similarity among the Mediterranean populations as well as among the extra-Mediterranean ones, suggesting the presence of a certain relatedness among the Mediterranean populations and between native and introduced B. bairdi populations as well.
The presence of connections between the different populations of the same species was demonstrated also in the median joining network of haplotypes. The network showed two distinct haplogroups representing Mediterranean and extra-Mediterranean detected mitochondrial haplotypes. In both haplogroups the occurrence of just one haplotype shared by the majority of the individuals confirmed the degree of connections between the Mediterranean populations of presumably B. boholense as well as between the native and introduced populations of B. bairdi. As regards the Mediterranean basin, the occurrence of two haplotypes shared by different populations suggested the presence of a low but existing gene flow between these populations likely mediated by local vessels rather than by natural dispersion of the species as supported by the in vitro short pelagic larval phase of the species [30]. The Mediterranean basin is highly urbanized and a dense network of harbors and marinas are present along its coasts [109,110], likely acting as sites of secondary dispersal of non-indigenous species [93,111]. Likewise, the relatedness among native populations of B. bairdi and individuals introduced in the Pacific Ocean and in Madeira Island was demonstrated by the fact that these populations share the same haplotype detected in native populations. This genetic structure corroborates the hypothesis that the introduction of B. bairdi both in the East Pacific and in the East Atlantic occurred from the Caribbean as already suggested by Tovar-Hernández et al. [45,46] and Ramalhosa et al. [42] for the port of Mazatlán and Madeira Island, respectively. Moreover, this scenario brings attention to the Panama Canal. The Panama Canal is an important point of passage for the maritime traffic with thousands of ships transit the Canal annually [112] and, as already reported by several studies [113][114][115][116][117] it is a potential corridor for the movement of marine species between the Pacific and Atlantic Ocean. It is worth mentioning that haplotype 9 is clearly distinctive from the remaining B. bairdi samples (average genetic distance, K2P = 0.015) (Figs 2 and 3). No morphological differences were observed and we therefore include this haplotype under B. bairdi. It comprises samples from Veracruz and Tampa Bay (native range) but also a single individual from the Galapagos Islands (introduced range). We have no information about the underlying causes for the genetic split, but the wide geographic distribution of this haplotype underlines its invasive potential.
In this study, the combined use of both molecular and morphological data allowed to clarify the taxonomic status of the two alien Branchiomma species and the genetic structure of the studied populations, providing reliable information on the distribution range of B. bairdi outside its native range and on the invasion scenario of (presumably) B. boholense in the Mediterranean Sea.
Further studies with the use of other molecular markers and a bigger pool of samples could corroborate our findings and provide a better understanding of two species invasion process for a complete and exhaustive management program.
Supporting information S1