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.

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₂-free buffer (Invitrogen, UK), 3 mM (for COI gene)/2.5 mM (for 16S, 28S and ITS genes) MgCl₂, 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₂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).

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⁷ 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₂) (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 demographic 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₂ 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.

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₂) 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.

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.

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). 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).

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.

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).

Figure 2.

Geographical distribution of the 39 mtDNA haplotypes (Hp) of Paracaprella pusilla in the populations sampled. Each site is represented by a pie chart showing population composition and relative haplotype frequency. Number of analysed individuals per population appears in brackets. White-shaded areas are the cumulative proportion of private haplotypes per location. Sites are coded as in Tables 1 and 2. The legend gives information about the existing haplotypes across all locations and the number of individuals carrying each haplotype. For the Cadiz marina (ESCAD) population, the change in haplotype frequency is shown in the four years when the species was recorded.

Overall, haplotype (Hd = 0.879) and nucleotide (π = 0.00418) diversities were high (Table 3). Both presumed native and non-native populations showed similar values of haplotype diversity (Table 3). Highest haplotype diversities were found for the presumed native populations of Paranaguá marina (BRPAR) and Niteroi (BRRIO), and the supposedly non-native population of Australia (AUAUS) (Hd = 0.810, 0.800 and 0.790, respectively). The lowest haplotype diversity was found for the presumed native population of Mexico (MXSIS) (Hd = 0.600; Table 3). Among Eastern Atlantic and Mediterranean populations, Palma marina (ESBAL) displayed the highest value of haplotype diversity (Hd = 0.879) and TNGGB the lowest (Hd = 0.667; Table 3). Nucleotide diversity ranged from 0.00066 in TNGGB to 0.00598 in ESCAD.

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 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).

Figure 3.

Multidimensional scaling plot (MDS) based on F_ST values for Paracaprella pusilla. For Cadiz marina (ESCAD) population, four points are represented, each one corresponding to one of the four years when the species was recorded. Populations are coloured according to the region they belong: Northeast Atlantic + Mediterranean (red); South Pacific (green); and presumed native region (blue).

Neutrality tests, Tajima’s D, Fu’s FS and Ramos-Onsis and Rozas’ R₂, 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.

Figure 4.

Mismatch distribution of Paracaprella pusilla for each region. a) Europe (Northeast Atlantic + Mediterranean), b) Australia (South Pacific), and c) presumed native region. Blue bars show the observed frequency distributions and the orange lines represent the expected ones under the sudden expansion model.

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 (Supp. 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 diversities of haplotype (Hd = 0.764) and nucleotide (π = 0.00598) were high (Suppl. material 4, Table S3). Haplotype diversity showed a linear pattern of decrease over time (Suppl. material 4, Table S3), with maximal values occurring after the species’ initial discovery at Cadiz marina in 2010 (Hd = 0.885), and lowest values occurred in the last monitoring year (Hd = 0.343). Nucleotide diversity ranged from 0.00169 in 2017 to 0.01515 in 2016 (Suppl. material 4, Table S3).

The F_ST statistics showed intermediate levels of divergence between years, with significant values ranging from 0.131 (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₂ were negative and not significant for all years (Table 7), which indicated that the Cadiz marina population was not under an expansion phase.

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).

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.

Figure 5.

Global distribution of Paracaprella pusilla and the genus Paracaprella a Current worldwide distribution of P. pusilla including its introduced range and the proposed native range. Information based on Ros and Guerra-García (2012) and references therein; Ros et al. (2013c); Ros et al. (2014); Ros et al. (2016a); Alarcón-Ortega et al. (2015); Alfaro-Montoya and Ramírez-Alvarado (2018); and Fersi et al. (2018) b Number of Paracaprella species recorded per marine ecoregions. Information based on Winfield and Ortiz (2013); Mauro and Serejo (2015); Sánchez-Moyano et al. (2015); Lacerda and Masunari (2014). Marine ecoregions follow Spalding et al. (2007).

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.

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 (Supp. 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 non-indigenous 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.