Unravelling the origin and introduction pattern of the tropical species Paracaprella pusilla Mayer, 1890 (Crustacea, Amphipoda, Caprellidae) in temperate European waters: first molecular insights from a spatial and temporal perspective

Paracaprella pusilla Mayer, 1890 is a tropical caprellid species recently introduced to the Eastern Atlantic coast of the Iberian Peninsula and the Mediterranean Sea. In this study, we used direct sequencing of mitochondrial (COI and 16S) and nuclear (28S and ITS) genes to compare genetic differences in presumed native and introduced populations in order to infer its introduction pattern and to shed light on the native range of this species. The temporal pattern of genetic diversity at the westernmost limit of the Copyright M. Pilar Cabezas et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. NeoBiota 47: 43–80 (2019) doi: 10.3897/neobiota.47.32408 http://neobiota.pensoft.net RESEARCH ARTICLE Advancing research on alien species and biological invasions A peer-reviewed open-access journal


Introduction
Non-indigenous species (NIS) are a fundamental component of global change and are currently considered one of the most important drivers of biodiversity alteration in marine ecosystems worldwide (Bax et al. 2003;Molnar et al. 2008). Some NIS successfully establish themselves, form self-sustaining populations, and spread into new locations, becoming invasive and causing both significant ecological and economic impacts (Molnar et al. 2008). They may out-compete native species and alter community structure and ecosystem processes. They may also threaten ecosystem services, which may result in significant economic losses in fisheries, aquaculture, and tourism sectors Ojaveer et al. 2015;Katsanevakis et al. 2016).
Marine organisms have been spread by human-mediated transport long before the first comprehensive biological studies were carried out (Carlton 1999(Carlton , 2003. However, the ever-increasing magnitude and efficiency of global maritime trade and associated transported vectors, as well as rising seawater temperatures associated with global climate change, are drastically increasing the spread of NIS (Ruiz et al. 1997;Carlton and Cohen 2003;Katsanevakis et al. 2013;Booth et al. 2017).
Europe, where approximately 1500 NIS have been introduced, is the major recipient of marine NIS worldwide (Katsanevakis et al. 2014;AquaNIS 2015;Tsiamis et al. 2018). Consequently, legislation mostly rooted in the Marine Strategy Framework Directive (MSFD) (EC 2008) and the Biodiversity Strategy (EC 2014) has been adopted to deal with NIS and, thus, protect, conserve, or enhance marine ecosystems. These strategies aim to mitigate or reverse the impacts of existing NIS and prevent future introduction and the establishment of new ones by identifying and managing introduction pathways, among other things. In this regard, genetic data have been recognized as a powerful and useful tool (Holland 2000;Geller et al. 2010;Rius et al. 2015). The study of the genetic structure and degree of gene flow within and between native and non-native populations specifically provides crucial insights into the pattern of introduction, colonization, and spread of introduced taxa (Geller et al. 2010;Rius et al. 2015). Such studies help to determine source populations (Rius et al. 2015) and provide insights into the invasive potential of species (Roman and Darling 2007), and, thus, can lead to a better understanding of the mechanisms and dynamics underlying introduction and invasion. Such information can be used to prioritize management strategies, to prevent further introduction events, and to assess the abundance and status of NIS.
Nonetheless, most genetic studies on NIS have focused on terrestrial and freshwater organisms (Lowry et al. 2013;Sherman et al. 2016;Viard et al. 2016), likely due to the logistics of sampling in the marine environment, which make it difficult to obtain reliable neutral population estimates (Sherman et al. 2016). Therefore, more research is needed to improve our knowledge of marine NIS and better understand the patterns of their introduction and invasion. On this subject, hotspots and stepping-stone areas for these species (e.g. marinas, aquaculture installations) and zones of special interest such as marine reserves or NATURA 2000 sites, should be prioritized (Olenin et al. 2010). In European waters, the Mediterranean Sea and the region of the Strait of Gibraltar deserve special attention, as they are considered hotspots for both biodiversity and biological invasions (Drake and Lodge 2004;Molnar et al. 2008;Boudouresque et al. 2017). Indeed, the Mediterranean Sea hosts the highest documented number of marine NIS globally, with around 900 NIS recorded so far (Ulman et al. 2017;Zenetos et al. 2017;Galil et al. 2018). Both the Mediterranean Sea and the Strait of Gibraltar have important shipping links to other areas worldwide, intense aquaculture activity, and are among the most important destinations for tourism worldwide, with many vessels arriving from America and northern Europe as well as from the Indo-Pacific via the Suez Canal (Streftaris et al. 2005;Galil et al. 2017;Tsiamis et al. 2018).
Crustaceans are among the most introduced taxa worldwide (Carlton 2011). In the Mediterranean Sea and nearby areas, they account nearly 20% of all reported NIS (Zenetos et al. 2012;Ulman et al. 2017). Within crustaceans, caprellid amphipods associated with fouling communities of artificial habitats are considered as prime candidates for introduction and establishment in regions where they are not native (Ashton et al. 2010;Ros et al. 2016a). The great abundances that some caprellids attain in these communities and their ability to survive on floating objects and vessel hulls (Thiel et al. 2003;Ashton et al. 2010) make them good models for understanding marine introductions and invasions. Thus, the number of introduced species belonging to this group have been documented worldwide with increasing frequency (e.g. Ros et al. 2016a;Gillon et al. 2017;Marchini and Cardeccia 2017).
Three non-indigenous species of caprellids have been recorded in temperate European seas: Caprella mutica Schurin, 1935(Ashton 2006Almón et al. 2014), Caprella scaura Templeton, 1836 (Sconfietti andDanesi 1996;Martinez and Adarraga 2008), and Paracaprella pusilla Mayer, 1890 (Ros andGuerra-García 2012;Ros et al. 2016a). While the dynamics of invasion by C. mutica and C. scaura have been explored through molecular tools (Ashton et al. 2008;Cabezas et al. 2014), there are no similar studies for P. pusilla. Paracaprella pusilla is a tropical caprellid species first described from Brazil (type locality: Rio de Janeiro) (Mayer 1890). Nonetheless, the native range of this species is not entirely clear (Farrapeira et al. 2011;Rocha et al. 2013). Some authors have considered the Atlantic coast of Central and South America as the most likely native range of this species (Mayer 1903;McCain 1968;Carlton and Eldredge 2009;Rocha et al. 2013). Paracaprella pusilla is frequently found in the Caribbean (Carlton and Eldredge 2009), with most records coming from the Gulf of Mexico and the coasts of Venezuela and Colombia (Díaz et al. 2005;Guerra-García 2006;Winfield et al. 2006). However, other authors have considered P. pusilla to be cryptogenic, sensu Carlton (1996), in this region (Serejo 1998;Farrapeira et al. 2011). Since its original description, P. pusilla has been reported from numerous other tropical and subtropical areas around the world, including both the East and West African coasts, the Indian peninsula, Australia, and Hawaii (see Ros and Guerra-García 2012 and references therein), mainly on fouling communities associated with artificial structures. Today, the known distribution of P. pusilla also includes the Pacific coasts of Panama (Ros et al. 2014), Mexico (Alarcón-Ortega et al. 2015, and Costa Rica (Alfaro-Montoya and Ramírez-Alvarado 2018), in addition to temperate European waters. Its first recorded occurrence (2010) in Europe was at Cádiz on the Atlantic coast of southwestern Spain (Ros and Guerra-García 2012). Soon afterwards, this species was first found in the western Mediterranean at Mallorca (2011) and Ibiza in August 2012 (Ros et al. 2013c). In 2014, several individuals of P. pusilla were found at Zikim on the southern coast of Israel, the first record of this species in the eastern Mediterranean Sea (Ros et al. 2016a). Most recently, P. pusilla has been reported in the central Mediterranean, in the Gulf of Gabès in tidal channels of the Kneiss archipelago (Tunisia) (Fersi et al. 2018). Thus, by its presence in western, central, and eastern regions of the Mediterranean, we suggest that P. pusilla might be present more generally throughout the Mediterranean and that it might have been overlooked due to its small size or temporal instability.
Two main pathways have been suggested for the introduction of P. pusilla to European waters. Ship fouling is the most probable vector for the introduction and dispersion of this species (Ros and Guerra-García 2012;Ros et al. 2013b, c), either through the Strait of Gibraltar, from source populations in Central and South America, where it is supposedly native, or via the Suez Canal, from the Indo-Pacific. This second alternative is less likely because P. pusilla has not yet been found in the Red Sea (Zeina and Guerra-García 2016) and has only been reported once in the Suez Canal (Schellenberg 1928) despite recent studies (see El-Komi 1998;Emara and Belal 2004;Zeina and Guerra-García 2016). These pathways are only assumptions based on historical records of P. pusilla. No molecular studies have been conducted to elucidate the introduction and dispersion patterns of this species.
Another question that remains genetically unexplored is whether P. pusilla is indeed a cosmopolitan species or if populations across its presumed large range belong to different cryptic species. In the order Amphipoda, molecular evidence supports the existence of cryptic species among widely distributed marine NIS, such as Ampithoe valida Smith, 1873 and Jassa marmorata Holmes, 1905, two biofouling species introduced to the Northeast Pacific (Pilgrim and Darling 2010). Nevertheless, other NIS are actually widely distributed (no cryptic species); examples include the caprellids Caprella mutica and Caprella scaura sensu stricto (Ashton et al. 2008;Cabezas et al. 2014). Morphological evidence supports the conspecificity of populations of P. pusilla (Ros et al. 2014), but molecular evidence is still needed to confirm P. pusilla as a neocosmopolitan species (introduced species that have achieved a widespread distribution through anthropogenic dispersal; sensu Darling and Carlton 2018).
In this study, we analysed the genetic diversity, population structure, and levels of differentiation of populations of P. pusilla from its presumed native and introduced distribution ranges. We sequenced mitochondrial and nuclear genes of P. pusilla in order to (i) provide the first molecular evidence to verify the conspecificity of populations; (ii) shed light on this species' native range, and (iii) to infer its introduction pattern in temperate European waters, particularly on the Iberian Peninsula. In addition, we analysed the temporal pattern of genetic diversity at Cadiz marina, which is the westernmost limit of the range of P. pusilla in Europe, beginning soon after its first detection and for a period of eight years. We use the Cadiz marina as a model for understanding how genetic diversity influences the introduction process of this tropical NIS into new areas where it previously could not survive. This information is crucial to better understanding the initial phases of marine introductions and identifying the factors associated with it. Additionally, this information allows for the better understanding of possible future invasions to other localities on the Atlantic coasts of Europe in the scenario of global warming, and, thus, it provides valuable information for the effective management of introduced species. As far as we know, this is the first study of temporal change in genetic diversity of an introduced marine amphipod.

Sample collection
Spatial sampling. A total of 230 specimens of P. pusilla were collected from 12 localities across its presumed native and introduced geographic ranges, including from the type locality at Rio de Janeiro and the whole of its introduced range in Europe (Table  1). Unfortunately, the Indo-Pacific region was represented only by a single sample from Australia. The greater number of individuals collected at Spanish localities is the result of our continuous monitoring over eight years, when at least two samplings per year were conducted. Samples were collected mostly from fouling communities predominantly comprised of hydroids and macroalgae, attached to floating pontoons, ropes, buoys, and ship hulls. At each locality, individuals were removed by hand and immediately preserved in 96-100% ethanol. In the laboratory, using a stereomicroscope, male individuals (see Guerra-García 2006) were identified as P. pusilla (see Mayer 1903: pl. 2, figs 36, 37;Ros et al. 2013c: fig. 2).
To compare the levels of intra-and interspecific genetic diversity, four individuals of the congeneric Paracaprella tenuis Mayer, 1903 from Celestún, Mexico (Table 1) were also included. Paracaprella tenuis is very similar to P. pusilla (see morphological characters by McCain 1968; Winfield and Ortiz 2013: table 1) and these species occur in sympatry in the northern Gulf of Mexico (Foster et al. 2004).
Temporal sampling. Paracaprella pusilla was first recorded in Europe in September 2010 on a floating pontoon at Cadiz marina, southern Spain, during a survey of peracarid crustaceans from harbours along the Strait of Gibraltar (Ros and Guerra-García 2012). At this locality, P. pusilla is associated with Eudendrium sp., a hydroid and a common component of fouling communities; these species seem to have a mutualistic relationship (Ros and Guerra-García 2012;Ros et al. 2013a). Samples of the Cadiz population were collected annually from 2010 to 2017 (Table 1), mostly during summer or early autumn, when caprellid abundance was generally greatest. Individuals of P. pusilla were removed by hand from samples of Eudendrium collected from the sides of floating pontoons, near the water surface. Seventy-five specimens were collected from this site and immediately preserved in 96-100% ethanol.

DNA extraction, amplification, and sequencing
Genomic DNA was extracted from gnathopods, pereopods, antennae and gills along one side of the body of each specimen sampled. We used the commercial kit PureLink Genomic DNA Mini Kit (Invitrogen, UK) according to the manufacturer's protocol.
The DNA was eluted in 120 µl of elution buffer and stored at −20 °C. Fragments of two mitochondrial (COI and 16S rRNA) and two nuclear (28SrRNA and ITS) genes were amplified by polymerase chain reaction (PCR), the latter two genes only for a subset of representative individuals of each population. PCR amplifications consisted of 25 µl reaction volumes containing 3 µl of template DNA, 10× MgCl 2 -free buffer (Invitrogen, UK), 3 mM (for COI gene)/2.5 mM (for 16S, 28S and ITS genes) MgCl 2 , 0.2 mM dNTPs, 1 µM of each primer, 0.1 µg µl -1 Bovine Serum Albumin (BSA, Promega, Madison, WI), 0.3 U Platinum Taq DNA polymerase (Invitrogen, UK), and double-distilled H 2 O to volume. Primers for amplification and PCR conditions are listed in Table 2. PCR product purification and unidirectional or bidirectional Sanger sequencing were provided by a commercial company (GENEWIZ, London, UK).

Sequence analysis
The resulting sequences were checked and edited using SEQUENCHER version 5.4.6 (Gene Codes Corporation, Ann Arbor, MI, USA). Mitochondrial COI sequences were translated into amino acids to search for stop codons that are indicative of the presence of pseudogenes. All sequences were thereafter deposited in GenBank (Suppl. material 1, Table S1). For mitochondrial (COI and 16S) and ITS genes, all sequences were aligned using MUSCLE (Edgar 2004) as implemented in MEGA version 7 (Kumar et al. 2016). Sequences of both mitochondrial genes were subsequently concatenated using the APE package (Paradis et al. 2004) in RStudio (RStudio Team 2016). For 28S gene, sequences were aligned using the MAFFT algorithm (Katoh and Standley 2013) and highly variable regions were eliminated from the analyses using GBLOCKS (Castresana 2000) with default parameters and allowing all gap positions. Uncorrected pairwise distances among haplotypes were calculated using MEGA version 7 (Kumar et al. 2016).

Spatial analysis
Phylogenetic reconstruction. Phylogenetic relationships were estimated using two model-based methods of phylogenetic inference to verify whether alternative topologies were supported by different tree-building approaches: Bayesian inference (BI) in MrBayes version 3.2.6 (Ronquist et al. 2012) and maximum likelihood (ML) in RAxML version 7.2.8 (Stamatakis 2008). These analyses were carried out for three sequence datasets: one using the concatenated sequences of the mitochondrial genes (COI + 16S), other using sequences of the nuclear 28S ribosomal gene, and the last using sequences of the nuclear ITS gene. Analyses were conducted using data partitions by codon (1+2+3) for the mitochondrial COI gene, to minimize saturation effects of codon positions on phylogenetic reconstructions (Salemi 2009) and to account for different rates of evolution of each one (Pond et al. 2009). Only one individual (or sequence) per haplotype was included in the phylogenetic analyses. Caprella liparotensis Haller, 1879 and Caprella danilevskii Czerniavski, 1868 were used as outgroups ( Table  1). The best-fit model of sequence evolution for the three datasets was estimated using PartitionFinder version 2.1.1 (Lanfear et al. 2016). According to the corrected Akaike Information Criterion (AICc) (Akaike 1974), the best models for the mitochondrial concatenated dataset were GTR+I (1 st partition), GTR (2 nd partition), GTR+G (3 rd partition), and GTR+G (4 th partition = 16S), and for the 28S and the ITS genes, the models GTR+I and GTR+G were selected, respectively. For BI analyses, two independent runs, of four chains each, were conducted for 2 × 10 7 generations (runs converged with average standard deviation of the split frequencies below 0.01). Trees and parameters were sampled every 1000 generations, with the heating parameter set to 0.25. The consensus (majority-rule) tree was estimated combining results from duplicated analyses, after discarding 25% of total samples as burn-in, determined from plotting log-likelihood values against generation time in Traver version 1.7.1 (Rambaut et al. 2018). For ML analyses, phylogenetic tree was calculated using the GTRGAMMAI model for the mitochondrial dataset, GTRCATI for the 28S gene, and GTRGAMMA for the ITS gene, and bootstrap supports were calculated using 1000 replicates. Consensus tree inferred for each molecular dataset was visualized and rooted using FigTree version 1.4.3 (Rambaut 2017).
Furthermore, relationships among mitochondrial haplotypes (using the concatenated dataset) were examined via a haplotype network using statistical parsimony method (Templeton et al. 1992) in TCS version 1.21 (Clement et al. 2000) with a 95% connection limit. The network was plotted with tcsBU (Santos et al. 2016).
Estimates of genetic diversity and population structure. Two measures of mtDNA diversity, haplotype (Hd, Nei 1987) and nucleotide diversity (π, Nei 1987), were estimated for each P. pusilla locality and region, using DnaSP version 6 (Rozas et al. 2017). Three distinct regions were considered: Northeast Atlantic + Mediterranean, South Pacific, and Western Atlantic (presumed native region). The single individual from Israel (ILZIK) was excluded from the analysis.
The genetic differentiation among populations was determined by means of the statistics F ST (Weir and Cockerham 1984) with Arlequin version 3.5.1.2 (Excoffier and Lischer 2010), using the pairwise differences distance method. Pairwise F ST values were calculated for the mitochondrial dataset, excluding the population with less than three individuals (i.e. ILZIK). Statistical significance was assessed through 10000 permutations, and a multidimensional scaling (MDS) analysis was performed on the matrix of F ST values for a graphical depiction of the structure using TIBCO STATISTICA version 13 (TIBCO Software Inc., CA, USA). Additionally, a hierarchical analysis of molecular variance (AMOVA, Excoffier et al. 1992) was conducted in Arlequin version 3.5.1.2 (Excoffier and Lischer 2010) to study the distribution of genetic variability between presumed native and introduced P. pusilla populations and to explore differentiation across geographic locations. In this regard, two groups were used for the AMOVA tests: (i) presumed native (Brazil + Mexico) vs non-native (Spain + Tunisia + Australia), and (ii) regions (Northeast Atlantic + Mediterranean; South Pacific; Western Atlantic). Statistical significance of variance components was tested with 16000 permutations.
Finally, to test if the selection of demographic events (population expansion or contraction) affected the genetic structure of non-native and potentially native populations, neutrality tests (Tajima's D, Fu's FS and Ramos-Onsis and Rozas' R 2 ) (Fu 1997;Tajima 1989;Rozas and Ramos-Onsins 2002) and mismatch distribution were performed for the mitochondrial dataset. Neutrality tests provide trends with respect to equilibrium and non-equilibrium conditions and indicate recent population expansion when the null hypothesis of neutrality is rejected due to significant negative values. They were assessed for each region with the statistical significance obtained by 10000 coalescent simulations. The distribution of frequencies of observed numbers of differences between pairs of haplotypes for each region is shown in the mismatch distribution. It uses tree shape to provide a rough estimate of population expansion or contraction because of a bottleneck. Populations that have experienced rapid demo-graphic growth in the recent past exhibit unimodal distributions, whereas populations that have been constant over time (demographic equilibrium) have bimodal or multimodal distributions (Rogers and Harpending 1992;Haydar et al. 2011). To test the goodness-of-fit between the observed and the expected distributions under the sudden expansion model, the sum of squared deviations (SSD) (Schneider and Excoffier 1999) and Harpending's raggedness index (Rg) (Harpending 1994) were also computed using 10000 bootstrap replicates. DnaSP version 6 (Rozas et al. 2017) was used to calculate R 2 statistic, and the remaining estimates and respective significance tests were obtained with Arlequin version 3.5.1.2 (Excoffier and Lischer 2010). In all analyses, localities with fewer than three individuals (i.e. ILZIK) were excluded.

Temporal monitoring
In the Cadiz marina (ESCAD) population, genetic diversity over time was assessed by estimating the haplotype (Hd) and nucleotide (π) diversity (Nei 1987) for each year sampled using DnaSP version 6 (Rozas et al. 2017). Frequencies of haplotypes per year were also calculated with this program. In addition, to test whether variation in genetic diversity (haplotype diversity) was linearly related to time (in years), a linear regression analysis was performed using the ggplot2 package (Wickham 2016) in RStudio (RStudio Team 2016).
Estimates of population differentiation over time were obtained from pairwise F ST calculations for the mitochondrial dataset, and neutrality tests (Tajima's D, Fu's FS and Ramos-Onsis and Rozas' R 2 ) were also estimated. All these analyses were conducted as described for the spatial analysis. The MDS analysis based on the matrix of F ST values was performed together with the data from the spatial analysis.

Sequence variation
The mitochondrial markers COI and 16S rRNA were successfully amplified for 236 caprellid individuals: 230 Paracaprella pusilla, four P. tenuis, and the outgroups Caprella liparotensis and C. danilevskii (Suppl. material 1, Table S1). Overall, 44 haplotypes were observed: 39 for P. pusilla, three for P. tenuis, and one for each of the outgroup species (Table 3). The complete alignment of the COI dataset had a total length of 612 bp. No insertions or deletions were detected in any of the sequences, and when they were translated into proteins, no stop codons were found. However, for the 16SrRNA (alignment of 408 bp), some indels were identified. Most of these correspond to insertions or deletions in sequences for P. tenuis or Caprella spp. Interestingly, among P. pusilla, a one bp insertion of a thymine (T) was observed at position 45 in six individuals from Australia and Brazil (all corresponding to the haplotype H25). The alignment of the concatenated dataset of these two genes (COI + 16S) had a total of 1022 bp.
The nuclear marker 28S was amplified for 60 P. pusilla individuals and the two outgroups species (Suppl. material 1, Table S1). Unfortunately, we were not able to amplify this gene for any individuals of P. tenuis. The total alignment length was 1265 bp, but only 1135 bp were selected using the software GBLOCKS. Some insertions and deletions were found, most of them distinguishing between P. pusilla and the outgroup species. Among P. pusilla sequences, a lack of sequence variation was observed: only two haplotypes, differing by the presence of an indel at position 443-444 of the alignment, were retrieved for the 60 individuals sequenced.
Finally, the alignment of the nuclear ITS marker had a total of 518 bp and included 73 P. pusilla and four P. tenuis individuals, plus the two outgroup species (Suppl. material 1, Table S1). Some indels were found between P. pusilla and the other species. However, all P. pusilla sequences were identical.

Spatial analysis
Phylogenetic reconstruction. Phylogenetic analyses of the mitochondrial dataset using the two different approaches (ML and BI) rendered trees with similar overall topologies, with main clades receiving high bootstrap or posterior probabilities support (Suppl. material 2, Fig. S1). All analyses revealed that P. pusilla and P. tenuis are monophyletic and formed highly supported clades. Within P. pusilla, no clear genetic structure was apparent and all haplotypes from the presumed native and non-native ranges appeared mixed, matching the results from the haplotype network (see below).  H6, H12, H13, H14, H15, H16, H17, H18 Nevertheless, haplotypes 6 and 7 appeared a little more differentiated from the remaining haplotypes. Pairwise divergence (uncorrected p distances) between most P. pusilla haplotypes was small, not exceeding 1.1%, with the exception of H6 and H7, which differed from other haplotypes by 2.0-2.4% (Suppl. material 3, Table S2). Interspecific divergence within the genus Paracaprella were much larger that intraspecific variation, ranging from 16.6 to 17.9%, which is similar to the values found between the two Caprella species included in this study (21.5%) (Suppl. material 3, Table S2). For the 28S gene, the ML and BI analysis produced a tree with identical topologies. Paracaprella pusilla was found to be monophyletic in both analyses (Suppl. material 5, Fig. S2A). Sequence divergence between them was 0%, whereas divergence between Caprella species was 4.5%. Divergence between the two genera exceeded 25%. Finally, for the ITS, the ML and BI analyses also rendered trees with identical topologies (Suppl. material 5, Fig. S2B). For this gene, sequence divergence between P. pusilla and P. tenuis was 16.9%, a higher value than that found between Caprella species (10.7%).
The haplotype network reconstruction for all sequenced mtDNA data retrieved two separate networks that could not be connected using the 95% parsimony connection limit (Fig. 1). For the first network, a star-like phylogeny was observed, with one very common haplotype surrounded by several low-frequency and some medium-frequency haplotypes with a maximum distance of eight mutation steps (corresponding to the haplotype H25). Interestingly, this most different haplotype corresponds to those individuals from Australia and Brazil that presented a one bp insertion at position 45 of the alignment. The central haplotype (H2) accounts for ~30% of P. pusilla individuals sequenced and was detected at all locations except the Gulf of Gabès, Tunisia (TNGGB; only three individuals sequenced) and Paranaguá Bay (BRPAB) ( Table 3). The remaining haplotypes found in Europe were within one to four point mutation steps from the central haplotype, those from Australia differed by one to eight steps (Fig. 1). The second network included only two haplotypes (H6 and H7) separated by seven mutation steps; these were detected in 10 individuals of Cadiz (ESCAD and ESSFN) and Baleares (ESBAL) populations (Fig. 1).
Genetic diversity and population structure. The spatial distribution of the 39 mitochondrial haplotypes of P. pusilla did not show any clear pattern (Table 3; Fig. 2). Private haplotypes were present in almost all populations. Eighteen haplotypes were identified in the presumed native range at the Atlantic coast of America, twelve of them were private (Table 3; Fig. 2). Among populations from other regions that are considered non-native (Europe and Australia), 25 haplotypes (14 private) were identified and only four of them (H2, H20, H22, and H25) were shared with the presumed native locations. Only H2 was shared between European non-native populations and the presumed native ones. Eight haplotypes (H1-H4, H6, H7, H9, and H19) were shared among non-native populations, most of them between ESCAD and ESSFN (i.e. Cadiz populations). Among these haplotypes, only one (H2) was present in almost all non-native populations (except TNGGB), but also in all presumed native locations except BRPAB (Table 3; Fig. 2).
The estimates of pairwise F ST values showed mostly low and intermediate levels of divergence between populations, with significant values ranging from 0.067 (ESCAD-AUAUS) to 0.538 (TNGGB-MXSIS) ( Table 4). Despite the great geographic distances, F ST values between presumed native and non-native populations were not high. However, they revealed that the Gulf of Gabès (TNGGB), Paranaguá Bay (BRPAB), and Sisal (MXSIS) were genetically differentiated from most other populations ( Table   Figure 1. Mitochondrial (COI+16S) haplotype network of Paracaprella pusilla from its presumed native and non-native range. Haplotypes 6 and 7, corresponding with 10 individuals of Cadiz (ESCAD and ESSFN) and Baleares (ESBAL) populations, were grouped in an independent network. This network could not be connected using the 95% parsimony connection limit to the main haplotype network which includes most of the haplotypes found in P. pusilla. Haplotype circles are proportional to haplotype frequency and numbers represent haplotype identities (Table 3). Non-observed haplotypes (extinct or unsampled haplotypes) are represented by small white circles. Each line connecting haplotypes represents a single mutational step. 4), although in the case of the Gulf of Gabès population this might be an effect of a low sample size (N = 3). These patterns are reflected in the MDS plot, which did not show any clear separation between non-native and presumed native populations, but TNGGB, BRPAB and MXSIS were slightly separated from the others (Fig. 3). Hierarchical AMOVA tests revealed significant genetic differences within populations, and among populations within groups at all geographical levels (native vs non-native, and  regions) ( Table 5). Intrapopulation variance explained most (over 80%) of the genetic variation found in P. pusilla (Table 5). Neutrality tests, Tajima's D, Fu's FS and Ramos-Onsis and Rozas' R 2 , were negative for all regions but not statistically significant ( Table 6); note that according to Fu (1997), FS statistic should be considered as significant if its p-value is less than 0.02. Additionally, the observed mismatch distribution was nearly bimodal for all regions (Fig. 4), which, thus, disproves the sudden expansion model and suggestes possible diminishing or structured population sizes. Regarding the sum of the square deviations (SSD), statistically significant differences were observed (p < 0.05) in presumed native populations and in the introduced Australian population (Table 6), which further support no recent population expansion. Nevertheless, these results contrast with the non-significant values of the Harpending's raggedness index (Rg) (Table 6), which indicated that a recent population expansion may have occurred in these populations. In addition, both SSD and non-significant values of Rg suggested goodness of fit between the observed and the expected distributions in East-Atlantic and Mediterranean introduced populations (Table 6; Fig. 4a), and, thus, the null hypothesis of recent population expansion should not be totally rejected.

Temporal monitoring
Paracaprella pusilla was monitored in Cadiz marina (ESCAD) soon after its first detection and for a period of eight years (2010-2017). However, the species was not found during the surveys carried out from 2012 to 2015. Therefore, we considered only four years (2010, 2011, 2016 and 2017) in our study.
Nine haplotypes (same as in the spatial study; Table 3), were obtained from the 75 individuals sequenced (Suppl. material 4, Table S3; Fig. 2). Interestingly, haplotype H4 was the only one present in all monitoring years, and it was the most frequent haplotype found (Suppl. material 4, Table S3). It was only detected in three individuals in 2010, but its frequency increased over time (Suppl. material 4, Table S3; Fig. 2). The remaining haplotypes were only detected in one or two specific years. For instance, haplotypes H1, H2, and H6 (the second the most common haplotype in the spatial study; Table 3; Fig. 2), were found in 2010 and their frequency increased in 2011, but they disappeared afterwards. Haplotype H7 was also present in 2010, not detected the next year, but detected again in 2016 and disappearing again in 2017. Finally, haplotypes H3, H5, H8 and H9 were only detected in one of the years. Overall, the     2010-2011) to 0.261 (2010-2017). Significant differentiation was found between years 2010, 2011 and 2017. Only the year 2016 did not show genetic differences from the other years during the monitoring period, but this could be an artefact due to the low sample size (N = 3). In the MDS plot, the year 2017 appeared more separated from the remaining monitoring years carried out in the Cadiz marina population (Fig. 3).
Finally, Tajima's D, Fu's FS and Ramos-Onsis and Rozas' R 2 were negative and not significant for all years (Table 7), which indicated that the Cadiz marina population was not under an expansion phase.

Conspecificity of Paracaprella pusilla populations
Unlike other caprellid taxa with a wide distribution, such as Caprella penantis (Cabezas et al. 2013a) or C. andreae (Cabezas et al. 2013b), the absence of population genetic structure ( Fig. 1; Suppl. material 2, Fig. S1), the small variation of studied mitochondrial markers, and no differentiation (except one indel in the 28S gene) in nuclear markers reveal that the populations of Paracaprella pusilla that we studied did not harbour any cryptic species. Moreover, our analyses confirm that P. pusilla and the morphologically close P. tenuis are monophyletic and formed highly supported clades (Suppl. material 2, Fig. S1; Suppl. material 5, Fig. S2B). Therefore, our findings support the assumption that anthropogenic dispersal is responsible for the broad geographic distribution of P. pusilla and confirms that this is a neocosmopolitan species (see Darling and Carlton 2018).

Native range of Paracaprella pusilla
The Atlantic coast of Central and South America has been postulated as the most likely native range for P. pusilla (Mayer 1903;McCain 1968;Carlton and Eldredge 2009;Rocha et al. 2013) (Fig. 5a). In our study, the six populations sequenced for this region and including the type locality (Rio de Janeiro) accounted for a higher percentage of private haplotypes (66.7%) than that found for all non-native populations sequenced (56.0%). This could be considered an indicator of long-term residency far exceeding the time-frame of human introductions (Wares 2002). However, the non-significant values of the Harpending's raggedness index (Rg) found for the native region, that may also imply a recent expansion of P. pusilla to this region, or the high genetic diversity found in some introduced areas, show the complexity to determine with accuracy the native region of the species through isolated approaches. For example, in an increasingly interconnected world, where maritime traffic continuously connects very distant areas, it is difficult to keep the native region isolated. Secondary introductions from populations introduced in remote areas, and, even more importantly, among sites within the native range, occur, which increase the connectivity and possibly also diversity within particular populations. This is particularly true in fouling species, such as P. pusilla, which can be found in both natural and artificial habitats in the native region. There are, however, several aspects that point to the Atlantic coast of Central and South America as the most likely native area for P. pusilla. First, most records of P. pusilla, both recent and old, come from this area (Ros and Guerra-García 2012). Second, while most records of P. pusilla from putative introduced areas are located in artificial habitats (such as those from India, Europe, Australia, Hawaii, and Pacific Mexico and Panama), in the putative native region P. pusilla is also common in natural habitats (Ros et al. 2016b). Third, the biogeographic distribution of species of Paracaprella (Fig. 5b) reveals that the Atlantic coast of Central and South America has a high diversity of recorded species, which infers that the centre of diversity for this genus may lie in this area. The other area with a high diversity of recorded Paracaprella species is a small region of the Pacific coast of Central America. However, all records of P. pusilla from this region are recent and, unlike other records of Paracaprella species, are located in artificial habitats (Alarcón-Ortega et al. 2015). For all these reasons, we believe that the Atlantic coast of Central and South America is the most likely native range for P. pusilla.

Introduction pattern in temperate European waters
Genetic studies have shown that introduced populations are generally much less diverse than the native ones because of the founder effects and post-introduction demographic bottlenecks (Holland 2000;Rius et al. 2015;Viard et al. 2016). However, many introduced populations do not exhibit reduction in genetic diversity and may even exceed native diversity as a result of admixture or high propagule pressure from multiple introductions events (Holland 2000;Roman and Darling 2007;Rius et al. 2015;Viard et al. 2016). The existence of multiple introductions has been widely reported in the marine environment (see Rius et al. 2015), including the two invasive caprellid species which have been genetically studied in Europe (Ashton et al. 2008;Cabezas et al. 2014). In our study, the high genetic diversity found within introduced European populations, similar to that seen in the native ones, coupled with the presence of numerous private haplotypes (Table 3; Fig. 2) suggest that the introduction of P. pusilla in temperate European waters likely occurred from multiple introduction pathways and source populations.
Our results support the existence of one of the two main introduction pathways previously suggested by Ros and Guerra-García (2012) and Ros et al. (2013b, c), that is, through the Strait of Gibraltar, from native populations of the Atlantic coast of America. Shipping routes have existed across the Atlantic for more than 500 years (Carlton 1989). Moreover, Europe and the Mediterranean, in particular, are characterized by large clusters of ports with intermediate to high levels of trade (Drake and Lodge 2004;Seebens et al. 2013), and it is estimated that approximately 8% of vessels that travel through the Strait of Gibraltar come from the Caribbean, Gulf of Mexico and Central America (Dobler 2002;Kaluza et al. 2010;Tsiamis et al. 2018). Our mitochondrial dataset suggests that populations from Brazil could be the source of European introduced populations, because their haplotypes grouped closely (Fig. 1) and also because of the lower levels of divergence among them (Table 4; Fig. 3). Moreover, our results indicate that there could exist at least two different introduction pathways through the Strait of Gibraltar: one responsible for the introduction of P. pusilla in the Iberian Peninsula (ESCAD, ESSFN, ESBAL) and Israel (ILZIK), and another responsible for the introduction of this species in Tunisia (TNGGB). The presence of haplotype H2, the most common and possibly the ancestral one given its central position in the network (Fig. 1), on both the native Western Atlantic coast and on the non-native European sites (Table 3; Fig. 2), as well as the close relation between haplotypes (Fig. 1) and the low level of genetic divergence among these populations (Table 4; Fig. 3), indicate a clear link between these two regions. According to our results, any of the Brazilian populations, except for Paranaguá Bay (BRPAB), could be the source of P. pusilla in European waters, as all of them shared some haplotypes with these introduced populations (Table 3; Fig. 2). In addition, although only three individuals of the Gulf of Gabès (TNGGB) population were sequenced, the exclusive presence of the haplotype H20 in this population but not in the others (where a high number of individuals were sequenced), the absence of haplotype H2 and the lack of any shared haplotypes (Table 3; Fig. 2), as well as the genetic divergence found among this and the European introduced populations (Table  4; Fig. 3) indicate that an independent introduction through the Strait of Gibraltar could have happened. Haplotype H20 only occurred in the native population of Paranaguá marina (BRPAR). This particular Brazilian area might, thus, have been the source for P. pusilla in the Gulf of Gabès.
On the other hand, although P. pusilla has not been reported in the Suez Canal since Schellenberg (1928) nor in the Red Sea (Zeina and Guerra-García 2016), its recent record in the Israeli coast (Ros et al. 2016a) and the fact that is one of the most abundant caprellid species along the coast of India (Guerra-García et al. 2010), lead us to think of the possibility of these locations as potential sources. Therefore, some individuals of P. pusilla could have been introduced to the Mediterranean region from the Indo-Pacific through the Suez Canal (Tsiamis et al. 2018). Future studies including samples from the Indo-Pacific region are necessary.
Paracaprella pusilla was reported for the first time in European waters in the fouling community of a marina on the Atlantic coast of southwest Spain (Ros and Guerra-García 2012) (Cadiz marina, ESCAD, in the present study), and only one year later (2011) the species was found for the first time in the western Mediterranean (ESBAL in the present study) (Ros et al. 2013c). So, according to historical records, the Cadiz marina population could represent the first step in the introduction pathways of this species in this region. However, our molecular results are not in general agreement with this hypothesis. The Palma marina (ESBAL) population had the greatest genetic diversity together with the greatest number of private haplotypes in the introduced range (Table 3; Fig. 2), which indicates that it, and not the Cadiz marina, could be the initial entry point of P. pusilla in European waters, and, thus, the source population for subsequent range expansion of this species in this region. Palma marina, in Balearic islands, is the largest port and an important point for commercial cargos, recreational boating, and commercial fishing; this port is one of the most important cruise destinations in the entire Mediterranean (Minchin et al. 2006) and a potential hot-spot of marine bioinvasions (Drake and Lodge 2004;Ros et al. 2013b). As far as we know, the studies of Ros et al. (2013a, c) were the only ones focused on caprellids associated with fouling communities in marinas and ports of Mallorca. So, it is possible that P. pusilla was present in Palma marina before its first record in the Cadiz marina. Our data show that European introduced populations are closely related. They shared haplotype H2 (Fig. 2), and the level of divergence between populations was relatively low (Table 4; Fig. 3), indicating that these populations are most likely stepping stones along the same introduction pathway, which is consistent with the scenario of transport by small vessels (Wasson et al. 2001). The stepping-stone invasion pattern is characteristic of many marine invasions and has been reported for other caprellids, such as Caprella mutica (Ashton 2006) and C. scaura (Cabezas et al. 2014). Thus, P. pusilla could have spread from Palma marina (ESBAL) to Cadiz marina (ESCAD), and from there to San Fernando (ESSFN), where the genetic diversity was less (Table 3). Many small vessels of Mallorca overwinter in marinas in southern Spain (Minchin et al. 2006). This, together with the high use of recreational boats on this island (Balaguer et al. 2011), represent a suitable vector for the secondary spread of P. pusilla from one location to another (Ros et al. 2013b, c). Unfortunately, only one individual of the eastern Mediterranean population (ILZIK) could be sequenced, which is insufficient to draw any conclusions about the source population and introduction pathway of P. pusilla at this locality.
Interestingly, the presence in Palma marina (ESBAL) population of one haplotype (H19) also found in Australia (AUAUS) ( Table 3; Fig. 2) suggests that the same pathway or source population may have been responsible for the introduction of P. pusilla at these localities. Australian population did not show significant genetic differences from the native region (except with MXSIS) (Table 4; Fig. 3), and half of the haplotypes detected among sequenced Australian individuals were shared with some of the Brazilian populations ( Fig. 2; Table 3). One of these haplotypes was H25, a highly distinct haplotype (Fig. 1) that has a specific insertion, which was otherwise observed only in Brazilian populations. All this indicates that the Atlantic coast of South America, namely Brazilian populations, could be the most likely origin for P. pusilla in Australian waters.

Temporal monitoring: loss of genetic diversity over time
Our monitoring of the Cadiz marina (ESCAD) population showed a progressive loss of genetic diversity over time (Suppl. material 4, Table S3; Fig. 2). This is consistent with a temporal instability of P. pusilla at this location at the westernmost limit of the geographic range of this species in Europe. This species was found in high densities, including ovigerous females, in September 2010, which somewhat refutes the presence of an initial bottleneck due to founder effects, as reported for other marine invertebrates (Pérez-Portela et al. 2012;Bayha et al. 2015). In fact, the high genetic diversity found in this population, comparable with the intrapopulation diversity observed in the presumed native range, together with the presence of private haplotypes (Table 3), indicates that a high number of colonizers arrived, probably from multiple source populations. After its first detection in September 2010, P. pusilla was recorded within the following two months. Then, it was not recorded until it was recorded again in the summer of 2011, associated with the absence and presence, respectively, of its main host, Eudendrium racemosum (Ros and Guerra-García 2012). In this year, some haplotypes (H1, H2, H4, and H6) found on the initial discovery of the species, were present with higher frequency, but other haplotypes disappeared (H5 and H7) (Fig. 2). These results could indicate that, even during periods when P. pusilla was not observed, some individuals could have persisted but remained undetected, due to low abundances (Carlton 2009), and re-established the population when favourable conditions (higher temperatures in summer and presence of E. racemosum) returned. Moreover, the presence of new haplotypes (H3 and H8) (Fig. 2) and the weak but significant differences observed between 2010 and 2011 (Fig. 3) could both indicate the existence of a previous bottleneck that unmasked the presence of these haplotypes or that new introductions from nearby populations also occurred. After 2011, P. pusilla was not observed during a period of several years until a few individuals appeared in December 2016. Interestingly, these individuals were detected after the finding of this species in a nearby marina (San Fernando, ESSFN in the present study) three months earlier. The presence of the haplotype H4 and the reappearance of the haplotype H7 in ESCAD population, both present in the ESSFN population (Table 3; Fig. 2), clearly suggest a link between these populations. In 2017, more individuals were found at the Cadiz marina and the presence of a new haplotype, only present in ESSFN population (H9), was observed (Fig. 2). This indicates that the reappearance of P. pusilla in the Cadiz marina is likely due to the arrival of new propagules from the ESSFN population that resulted in a successful establishment. The decrease in genetic diversity observed (Suppl. material 4, Table S3) consistent with the increasing dominance of one haplotype (H4) (which was not the commonest in the introduced range), indicates that the Cadiz marina population was re-established by a small number of founding individuals ("founder effect"; see Novak 2007;Pérez-Portela et al. 2012;Bayha et al. 2015;Rius et al. 2015). The founder effect is expected to influence the likelihood of long-time survival of NIS, either by reducing the evolutionary potential for adaptation to novel habitat conditions (Sakai et al. 2001;Willi et al. 2006) which inhibit adaptative potential, or by exposing populations to the negative effects of inbreeding (Sakai et al. 2001;Charlesworth and Willis 2009), and the success of the introduction and invasion can be significantly compromised as a result (Sakai et al. 2001;Novak 2007;Wellband et al. 2017). However, some studies have shown that this is not always true, and that low levels of genetic diversity do not prevent the success and spread of nonindigenous species (Roman and Darling 2007;Pérez-Portela et al. 2012;Bariche et al. 2017). Results from the neutrality tests (Table 7) indicate that the Cadiz marina population was not under an expansion phase. This, together with the absence of P. pusilla for five consecutive years, suggest that the long-term establishment and success of this species may be compromised. Rather than founder effect, the instability of P. pusilla at its westernmost limit in Europe, could be the result of ecological and environmental factors, one of them being the water temperature, as P. pusilla is a tropical species and has been mostly found in summer months in its introduced range (Ros and Guerra-García 2012; Ros et al. 2013a). However, the high recreational boating pressure that occurs in this area could increase genetic diversity over time, increasing the likelihood of local adaptation and therefore allowing the expansion of its invaded range. It would be interesting if further temporal genetic analysis could be addressed, preferentially monthly, in all European introduced populations to determine the current status of genetic diversity in this species, and thus, to fully understand how genetic diversity is influencing its introduction process.
Together with the increase in maritime traffic, climate change directly or indirectly increases the spread of NIS into new areas (Carlton 2000;Molinos et al. 2016;Hulme 2017), some of them establishing populations where they previously could not survive (Carlton 2000;Hellmann et al. 2008;Mellin et al. 2016). The increase in the average surface seawater temperatures of the Mediterranean during the last two decades has affected the distribution and abundance of native and non-native species, leading to an enlarged pool of non-native species that have become established and expanded their distributions (Coll et al. 2010;Ulman et al. 2017). This explains why tropical species, such as P. pusilla, are penetrating into temperate ecosystems. In fact, the high number of private haplotypes and the star-shaped haplotype network observed in the present study, seems to be a result of the propagule pressure from the species' range in the tropical regions. The occurrence of P. pusilla inside marinas, its association with the fouling communities of ships, its ability to spread locally by rafting on detached fragments of these fouling communities (Ros et al. 2016a), and its fecundity (greater than another caprellid, Caprella scaura, introduced in this area; Ros et al. 2013c), suggest possible future introductions to other Mediterranean and adjacent localities.

Conclusions
Our study constitutes the first molecular approach to verify P. pusilla as a neocosmopolitan species, which has been introduced in European waters from multiple introduction pathways likely including at least populations from Brazil. Molecular, ecological and biogeographic evidences point to the Atlantic coast of Central and South America as the likely native range of P. pusilla. While the species appears to be expanding in the Mediterranean, populations from the westernmost distribution limit in Europe (the Atlantic coast of southern Spain) showed a temporal instability. This may indicate that P. pusilla is not fully adapted to the environmental conditions in this area, with a water temperature cooler than in the Mediterranean. Further intensive sampling including both native (especially Caribbean populations) and non-native populations of this species, as well as temporal genetic studies, are still necessary to improve knowledge on the diversity of this species in its native and introduced range, confirm the introduction pattern suggested here, and understand the ecological and evolutionary process involved in the introduction success or failure of this species in European waters.