Natural or Naturalized? Phylogeography Suggests That the Abundant Sea Urchin Arbacia lixula Is a Recent Colonizer of the Mediterranean

We present the global phylogeography of the black sea urchin Arbacia lixula, an amphi-Atlantic echinoid with potential to strongly impact shallow rocky ecosystems. Sequences of the mitochondrial cytochrome c oxidase gene of 604 specimens from 24 localities were obtained, covering most of the distribution area of the species, including the Mediterranean and both shores of the Atlantic. Genetic diversity measures, phylogeographic patterns, demographic parameters and population differentiation were analysed. We found high haplotype diversity but relatively low nucleotide diversity, with 176 haplotypes grouped within three haplogroups: one is shared between Eastern Atlantic (including Mediterranean) and Brazilian populations, the second is found in Eastern Atlantic and the Mediterranean and the third is exclusively from Brazil. Significant genetic differentiation was found between Brazilian, Eastern Atlantic and Mediterranean regions, but no differentiation was found among Mediterranean sub-basins or among Eastern Atlantic sub-regions. The star-shaped topology of the haplotype network and the unimodal mismatch distributions of Mediterranean and Eastern Atlantic samples suggest that these populations have suffered very recent demographic expansions. These expansions could be dated 94–205 kya in the Mediterranean, and 31–67 kya in the Eastern Atlantic. In contrast, Brazilian populations did not show any signature of population expansion. Our results indicate that all populations of A. lixula constitute a single species. The Brazilian populations probably diverged from an Eastern Atlantic stock. The present-day genetic structure of the species in Eastern Atlantic and the Mediterranean is shaped by very recent demographic processes. Our results support the view (backed by the lack of fossil record) that A. lixula is a recent thermophilous colonizer which spread throughout the Mediterranean during a warm period of the Pleistocene, probably during the last interglacial. Implications for the possible future impact of A. lixula on shallow Mediterranean ecosystems in the context of global warming trends must be considered.


Introduction
The European black sea urchin Arbacia lixula (Linnaeus, 1758) is currently one of the most abundant echinoids in shallow rocky habitats of the Mediterranean [1], where it has the potential to greatly influence benthic communities with their grazing activity [2][3][4]. A. lixula has a considerable trophic plasticity, ranging from omnivory to strict carnivory [5] and its scraping predatory behaviour can bulldoze the substrate bare of erect and encrusting algae and sessile animals. A. lixula broadly overlaps its habitat with the common edible sea urchin Paracentrotus lividus (Lamarck, 1816). Both species are traditionally thought to have the ability to trigger the development of subtidal barren zones of reduced benthic productivity and diversity [6][7][8][9]. However, new and increasing evidence suggests that A. lixula could actually be playing the principal role in producing and maintaining these barrens [10] and that this trend could be worsening in the near future due to foreseeable climatic changes [11].
Arbacia lixula is commonly regarded as a typical native species in the Mediterranean fauna [12], since it is currently found in shallow rocky shores all along the Mediterranean, often at high densities, and has been so since historical times. However, its tropical affinities have been suggested for a long time. Based on the lack of Mediterranean fossil record, Stefanini [13] and Mortensen [14] stated that A. lixula (reported as A. pustulosa), probably originated at the Tropical Atlantic region, from where it spread into the Mediterranean. Kempf [15], Tortonese [16] and Fenaux [17] also considered that A. lixula was a thermophilous species.
In NW Mediterranean, increasing abundances over time have been reported for this species. In 1950, Petit et al. reported that Arbacia lixula had become abundant in Marseilles during the previous 30 years [18], despite Marion had described it as rare in the same area in 1883 [19]. More recently, Francour et al. reported a 12-fold increase in the abundance of A. lixula in Corsica over a period of nine years (1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992) and speculated that a long term rise in the water temperature could have been the cause for this proliferation [20]. In the same period (1982 to 1995), a 5-fold increase in A. lixula densities was reported at the Port-Cros Marine Reserve (France) [21]. On the other hand, in a recent 5-year follow-up (2003-2008) at Ustica Island (Southern Thyrrenian Basin), a positive correlation was found between the gonadosomatic index of adult A. lixula and summer surface water temperature, suggesting increased reproductive potential with temperature [10].
Arbacia is an ancient genus with a fossil record that dates back to the Paleocene [22] whose distribution is mainly Neotropical. Unlike other sea urchin genera, Arbacia has a history of latitudinal shifts [23], and the five extant species inhabit mainly temperate and tropical shallow waters [24], being mostly allopatric. Only one species, A. dufresnii, is able to live in cold Subantarctic waters. A. lixula is the only species in the genus that lives in the Old World. Its present distribution includes Brazil, the African Atlantic coast from Morocco to Angola, the East Atlantic archipelagos of Cape Verde, Canaries, Madeira and Azores, and the whole Mediterranean basin, excluding the Black Sea. It has never been reported from the Atlantic European coast north of Gibraltar (J. Cristobo, X. Troncoso, N. V. Rodrigues; pers. comms.), probably due to the low sea surface temperature originated by the southward Portugal Current [25].
Recently, Lessios et al. [26] presented an exhaustive phylogenetic study of genus Arbacia, using sequences of the mitochondrial COI (cytochrome c oxidase I) and the nuclear gamete recognition protein bindin, which has clarified many interesting questions on inter-specific relationships within this remarkable genus. Notably, the sequence of speciation events was consistently reconstructed and their divergence times were reliably estimated. Thus, the splitting between A. lixula and its sister species, the NW Atlantic A. punctulata, was estimated to have taken place some 2.2-3.0 Mya (millions years ago) based on COI sequences, or 1.9-3.3 Mya based on bindin sequences. The phylogeny of bindin sequences also allowed these authors to infer that Brazil populations separated from the rest of A. lixula some 1.8-3.4 Mya; i.e. very early in the evolution of this species (however, only 5 individuals from Brazil were used in the analysis, and no estimation could be inferred for the same event from mitochondrial sequences, due to the unresolved position of the Brazilian clade within other A. lixula haplotypes).
Yet, many questions remain open about the intra-specific relationships of Arbacia lixula. Considering its unusually wide present distribution area, which ranges from equatorial waters to temperate Mediterranean, the great colonizing potential shown by this species, including the ability to cross trans-oceanic barriers to gene flow [26], and the massive potential impact of its behaviour on coastal ecosystems, further research on its phylogeography and population genetics is necessary in order to elucidate the history and ongoing processes that shape the distribution of the species. In this work, we present a phylogeographic study using the mitochondrial marker COI, based on a representative sample of individuals covering most of the distribution area of Arbacia lixula. Our goals were to answer relevant questions concerning the history and present-day distribution of the species: What are the relationships between the main geographic areas where the species is found? Do the main geographic barriers to gene flow, that are known to regulate the genetic structure of many other marine organisms, affect the present-day genetic structure of this species? Can recent geographic and/or population expansion events be traced and reconstructed by analysing the signature left in sequence data of this species?

Ethics Statement
Field sampling required for this work involved only invertebrate species which are neither endangered nor protected. All necessary permits for sampling at localities placed inside protected areas (Cabrera National Park, Columbretes Islands Marine Reserve & Scandola Nature Reserve) were previously obtained from the competent authorities. Non-destructive sampling techniques (external soft tissue biopsy) were used in these localities in order to minimize impact on the ecosystems.

Sampling
Between April 2009 and July 2011, we obtained samples from 24 localities belonging to three predefined regions: West Atlantic, East Atlantic and Mediterranean (see Fig. 1

DNA Amplification and Sequencing
Total DNA was extracted using REDExtract-N-Amp Tissue kit (Sigma-Aldrich, www.sigma.com) from either one tube foot or a tiny portion (5-10 mg) of gonad. A fragment of the COI gene was amplified and sequenced using specific primers designed using the complete genome sequence of A. lixula mitochondrion [27] with PRIMER 3.0 [28], as follows: COIARB-F: 59-TTC TCT GCT TCA  AGA TGA C-39, COIARB-R: 59-CTA TAA TCA TAG TCG  CTG CT-39, COIAL-R: 59-GCT CGG GTA TCT AGG TCC  AT-39. Most individuals were amplified using the COIARB-F/ COIARB-R pair, but some individuals belonging to Atlantic populations had to be amplified using COIARB-F/COIAL-R instead. PCR amplification reactions were performed in a 20 ml total-reaction volume with 10 ml of REDExtract-N-Amp PCR reaction mix (Sigma-Aldrich), 0.8 ml of each primer (10 mM), 6.4 ml of ultrapure water (Sigma-Aldrich) and 2 ml of template DNA. A single denaturing step at 94uC for 5 min was followed by 40 cycles (denaturation at 94uC for 40 s, annealing at 43uC for 45 s and extension at 72uC for 45 s) and a final extension at 72uC for 5 min in a S1000 dual thermal cycler (BioRad, www.bio-rad. com). The PCR products were purified and both strands sequenced in Macrogen (www.macrogen.com) using the same primers for the sequencing reaction.

Genetic Diversity Analyses
All the sequences were edited in BIOEDIT [29] and aligned using CLUSTALW as implemented in MEGA 5 [30]. The single nucleotide mutations found were double-checked by contrasting the agreement and quality of forward and reverse sequencing chromatograms. The Nei & Gojobori procedure with the Jukes & Cantor correction [31][32] implemented in MEGA 5 was used for detecting positive natural selection. Sequences of the haplotypes found have been deposited in GenBank (accession numbers from JQ745096 to JQ745256).
Number of haplotypes (N h ), haplotype diversity (H d ) and nucleotide diversity (p) were computed with DNASP v. 5.10 [33]. Haplotype richness was calculated with CONTRIB v. 1.02 [34] using a rarefaction size equal to the smallest sample size (n = 15) and Student's t-test was used for comparing its values between regions having more than two sampled locations (i.e., Eastern Atlantic and Mediterranean).
We used BAPS v. 5.2 (Bayesian Analysis of Population Structure) [35][36] for clustering the sampled haplotypes into monophyletic clusters of haplotypes (haplogroups). We ran five replicates for every value of the maximum number of clusters (k) up to k = 10.
Haplotypes were assigned to one of the clusters by admixture analysis, performing 50 simulations from posterior haplotype frequencies. The assigned haplotype names reflect the haplogroup they belong.

Phylogeography and Phylogeny
Relationships and geographical distribution of the haplotypes were analysed in a haplotype network constructed with NETWORK v. 4.6.0.0 (http://www.fluxus-engineering.com/sharenet.htm), which implements the median-joining method, in the absence of recombination [37]. The network was optimized using maximum parsimony criterion and the obtained loops were solved using criteria derived from coalescent theory [38][39]. In order to determine the putative ancestral haplotypes, the outgroup weights based on haplotype frequency and connectivity [40] were calculated for each haplotype using the TCS v. 1.21 program [41].
For phylogenetic analysis of the haplotypes obtained, we included a sequence of Strongylocentrotus purpuratus from GenBank (Acc. number NC_001453 [42]). Though the use of an outgroup sequence for rooting intraspecific genealogies has been shown to have little resolution [43], we nevertheless used it since the resulting tree is coherent with the outgroup weights calculated using TCS. We used JMODELTEST v. 0.1.1 [44], based on a hierarchical series of likelihood ratio tests [45] and the Bayesian Information Criterion (BIC), to assess the most appropriate nucleotide substitution model for our data. This condition was satisfied by the Tamura & Nei model [46] with a gamma correction (a = 0.240) (TrN + G). This evolution model was fed into MRBAYES software v. 3.1.2 [47] and the haplotype tree was estimated under the BIC after 1 million generations of 8 MCMC chains with a sample frequency of 100 (10,000 final trees). After verifying that stationarity had been reached, the first 2,000 trees were discarded, an independent majority-rule consensus tree was generated from the remaining (8,000 trees), and it was drawn using MESQUITE v. 2.75 [48].

Population Structure Analyses
Pairwise genetic distances between populations (F st ) were calculated with ARLEQUIN v. 3.1 [49] considering the genetic distance between haplotypes, and their significances were tested by performing 40,000 permutations. The level of significance for these multiple tests was corrected by applying the B-Y false discovery rate (FRD) procedure [50][51]. Kruskal's non-metric multidimensional scaling (MDS [52]) of F st values was performed with RSTUDIO [53] to graphically visualise these results. In order to have a different differentiation measure based only on haplotype frequencies, Jost's D [54] was calculated using SPADE [55]. Negative values for D were corrected to zero. We calculated a confidence interval around the obtained values by 1,000 bootstrap replicates. We set this confidence interval, using the normal approximation, at the appropriate P-value following the B-Y correction as explained above. Significant differentiation was inferred when this confidence interval excluded zero.
Analyses of molecular variance (AMOVA) were performed to assess population structure, using conventional F-statistics (i.e. only with haplotype frequencies), and their significances were tested running 90,000 permutations in ARLEQUIN [56]. AMOVAs were performed using different population sets in order to test the significance of population structure among regions, or among subbasins within regions. These AMOVAs were repeated also considering genetic distances between haplotypes, in order to check the robustness of the results.
The effect of isolation by geographical distance was assessed, for the whole dataset or separately for different populations sets, by the correlation of linearized genetic distances (F st /1-F st ) [57] with geographical distances between localities. Though ideally the oceanic current patterns should be included in the geographical distances calculation, currently we do not know of any reliable method for accurately quantifying this, so we used the shortest distance by sea on GOOGLE EARTH 6 (http://www.google.com/ earth). The significance of the correlation was tested by the Mantel test procedure [58], implemented in ARLEQUIN, with 20,000 permutations for each analysis.

Demographic History Inference
Demographic history was inferred for the three studied regions and for each sub-basin by analysing the mismatch distributions. Populations that have recently experienced a sudden demographic growth show unimodal distributions, whereas those at demographic equilibrium show multimodal distributions [59]. The expected mismatch distributions under a sudden expansion model were computed in ARLEQUIN using Monte Carlo simulations with 10,000 random samples. The sum of squared deviations (SSD) between observed and expected distributions was used as a measure of fit, and the probability of obtaining a simulated SSD greater than or equal to the expected was computed by randomisation. If this probability was .0.05, the expansion model was accepted, and its parameters h 0 , h 1 and t were calculated. For those populations showing large values for the final effective population size h 1 , this method does not usually converge and flawed results could be obtained. In this case, we kept the value of t calculated by this method, which is consistently robust [60], and used DNASP to calculate the value of h 0 which minimized the SSD, letting h 1 have an arbitrary large value of 1000 [61]. In the case that the mismatch distribution was not unimodal, the data were fitted to a constant population size model [62][63] for graphical representation.
To estimate the approximate time of a demographic expansion (t) from coalescence methods, the relationship t = 2 mkt was used [59] where t is the mode of the mismatch distribution, m is the mutation rate per nucleotide and k is the number of nucleotides of the analysed fragment. A range of mutation rates from 1.6% to 3.5% per million years was used for the COI gene, as calculated previously for echinoids [64][65].
In order to add more statistical support for population expansions, Tajima's D test of neutrality [66], Fu's F s [67], and Ramos-Onsins & Rozas' R 2 [68] indices of population expansion were calculated using DNASP. The confidence limits of Tajima's D were obtained assuming that it follows the beta distribution [66], while statistical tests and confidence intervals for F s and R 2 were based on a coalescent simulation algorithm implemented in DNASP, with 20,000 simulations. Harpending's raggedness index r [69] was calculated using ARLEQUIN and its significance was tested using parametric bootstrapping (10,000 replicates). These indices were calculated for the three regions and the six predefined subregions.

Genetic Diversity
We sequenced 635 bp of the mitochondrial gene COI from 604 Arbacia lixula individuals from 24 localities ( Fig. 1 and Table 1). We

Haplotype Network and Phylogenetic Inference
The haplotype network (Fig. 2) showed a strikingly star-shaped topology with a high ratio of singletons (81.4% of all haplotypes), which is typical of populations that have suffered a recent demographic expansion. The three most abundant haplotypes (A2, A17, B6) occupy central positions. All initial loops obtained by the MP criterion could be resolved using coalescent theory, except one, comprising 2 of the most frequent haplotypes (A2, A17), plus haplotypes, A4 & A20, which is therefore left unresolved in the figure. The outgroup weights calculated by the TCS program identified A2 as the ancestral haplotype (Table S1). This is the second most frequent haplotype and the only which is  (Fig. 3) is coherent with the topology of the haplotype network. Haplotypes belonging to haplogroup A were collapsed at the base of the phylogram, indicating that this group is paraphyletic and ancestral, in accordance with the results of the outgroup weights analysis. Haplotypes of group B form a homogenous clade from which haplogroup C derives. The collapsed comb-like shape of haplogroups A and B suggests a recent demographic expansion. Interestingly, Brazilian haplotypes B44, B45 & B46 formed a monophyletic clade with haplogroup C, supported by a PP value of 0.81. This is consistent with previous results by Lessios et al. [26] which found that the samples from Brazil included in their analysis formed a clade nested within Eastern Atlantic (and Mediterranean) sequences.

Population Structure
The analyses of population pairwise genetic differentiation (F st and Jost's D, Table 3) reflected a lack of population structure within both Eastern Atlantic and Mediterranean regions, but a clear differentiation between them and a complete differentiation (no alleles shared) of both regions from the Brazilian samples. Results from F st and D were largely consistent. No significant differences could be found between any pair of localities from Cape Verde and Macaronesia, suggesting a high level of genetic flow among these Eastern Atlantic sub-regions. Likewise, no significant differences were found between any pair of Mediterranean localities (out of 136 possible pairs), with the exception of Torremuelle (the westernmost Mediterranean locality) where F st analysis showed significant differences with two other Mediterranean localities, though these differences were not significant when D measures were analysed. Between Eastern Atlantic and Mediterranean, however, 38 (D) and 31 (F st ) comparisons (out of 85) were significant. Remarkably, the localities of Carboneras (Western Mediterranean), Crete and Kos (Eastern Mediterranean) did not show any significant difference to any other Eastern Atlantic or Mediterranean population, despite the large geographical distances involved in the case of the two latter localities. The MDS analysis (Fig. 4) graphically expresses the relationships among populations obtained from F st measures. Brazilian localities are widely separated in the first dimension from Eastern Atlantic and Mediterranean populations, whereas the Mediterranean and Eastern Atlantic populations were separated along the second axis. The lack of structure between sub-regions within the Eastern Atlantic and the Mediterranean is also apparent in the graphical arrangement. The same analysis using D measures (not shown) reflected the same overall structure.
Consistent with the pairwise differentiation analysis, the AMOVA found significant differences between the three regions (Table 4), which remained significant when only Eastern Atlantic vs. Mediterranean regions were compared (Table 5). Conversely, and again in agreement with the pairwise differentiation analyses, no significant differences within regions between Eastern Atlantic sub-regions (Table 6) or among the three Mediterranean subbasins (Table 7) were detected by AMOVA. The same results were obtained when these AMOVAs were repeated considering genetic distances between haplotypes (data not shown).
The Mantel test showed significant isolation by distance when the whole dataset was analyzed (Fig. 5A). This result remained significant when populations from Brazil were excluded (Fig. 5B). Contrarily, no significant correlation between genetic differentiation and geographical distance was found when populations within just one region, either East Atlantic or Mediterranean, were analyzed ( Fig. 5C & 5D).

Historical Demography
The mismatch distribution of Arbacia lixula populations from the Brazilian region (Fig. 6A) did not fit the sudden expansion model (Table 8). Conversely, the mismatch distribution for the Eastern Atlantic region (Fig. 6B) was remarkably unimodal. This indicates that a recent demographic expansion has occurred in this population. Similar results were obtained when only the Macaronesian sub-region was analyzed (Table 8). However, the distribution for the Cape Verde sub-basin did not fit the sudden expansion model, as reflected by a high SSD (Table 8). Nevertheless, this result may be an artefact due to small sample size (n = 27). The demographic expansion in the Eastern Atlantic populations could be dated, from the value of t and the known mutation rate for the COI of Echinoidea, between 30.6-66.9 kya (thousand years ago), which is a surprisingly recent time.
The mismatch distribution obtained for the Mediterranean region (Fig. 6C) was also typically unimodal. The parameters of the theoretical curves calculated individually for each Mediterranean sub-basin had all similar values, comparable to those of the whole Mediterranean region (Table 8), reinforcing the idea that all the Mediterranean Arbacia lixula populations belong to the same genetic pool. The demographic expansion in the Mediterranean could be dated between 93.8-205.2 kya. This estimation is a little older than that obtained for the Eastern Atlantic expansion, but is still a recent time.
The neutrality and population expansion tests calculated for the different regions and sub-basins (Table 9) were largely coherent with the results inferred from the mismatch distributions. Tajima's D detected significant differences from neutrality in all cases, except for Brazil and the Eastern Mediterranean sub-basin. Fu's F s test for demographic expansion was significant in all cases (though just marginally so in the case of Brazil). Ramos-Onsins & Rozas' R 2 was significant for all cases except the Eastern Mediterranean sub-basin, and the raggedness value r was consistent with unimodal distributions, except for Brazilian and Cape Verdean populations.

Discussion
COI and other mitochondrial markers have proven to be the most useful tool for tracing both intraspecific and intrageneric genealogies of many echinoid species [26,[64][65][70][71][72][73] and usually yield easily interpretable results which are consistent with  Table 3. Genetic differentiation between Arbacia lixula populations, F st (below the diagonal) and Jost's D (above the diagonal). Consistently significant differences obtained by both methods after false discovery rate correction are represented in bold.
Significant P values for F st obtained from randomization. *: significant after false discovery rate correction (P,0.0085).
Significant P values for D indicate that confidence interval obtained by bootstrapping excludes 0.
*: significant after false discovery rate correction (P,0.0085). doi:10.1371/journal.pone.0045067.t003 those of other nuclear markers. Nevertheless, our analyses are based on a single mitochondrial marker (COI). Thus, these results must be taken with caution, and further analyses using nuclear markers would be desirable. On the other hand, previous works in Echinoidea have shown that other nuclear markers were mainly used only to confirm the evolutionary history depicted by mtDNA [26,72] or else displayed too much diversity to produce interpretable results [73]. The Arbacia lixula populations sampled showed high values of haplotype diversity and haplotype richness, but relatively low values of nucleotide diversity. The lowest diversity was found in Brazilian populations and, specifically, in the westernmost locality (Itaipu), which is close to the distribution limit of the species and separated from the other Brazilian locality by the Cabo Frio upwelling. In contrast, the highest diversity was found in the East Atlantic, as expected if this region is the geographical origin of the species [16,26]. We detected three haplogroups in A. lixula. One of them (Group A) seems to be ancestral and is found only in Eastern Atlantic and Mediterranean populations, while another (Group B) is present at both sides of the Atlantic. The third one (Group C) is derived from Group B and found only in Brazil.
In a recent work, Lessios et al. [26] concluded that Arbacia lixula split from a common ancestor with A. punctulata ca. 2.6 Mya, and attributed this split to the mid-Atlantic barrier, separating the western A. punctulata from the eastern A. lixula, which would later have crossed back this barrier to establish itself, as an isolated  clade, in the coast of Brazil. A problem with this view is that the mid-Atlantic barrier was fully in place long before the estimated date of the split, so the separation of the two species could not be a vicariance event but a range expansion event (on the part of the lineage that would become A. lixula), and two crossings of the barrier are required to fully explain the present-day distribution of the species (though the second crossing could be facilitated by the South Equatorial Current system [74]). An alternative scenario would be that the two Atlantic species diverged in Western Atlantic, after the rise of the Panama isthmus isolated their ancestor from the eastern Pacific region (the possible origin of the genus Arbacia [26]), and that A. lixula crossed the Atlantic ridge only once to colonize the Eastern Atlantic. Our results favour the first (Lessios') view, as the haplotypes from Brazil formed a derived monophyletic group nested within the amphi-Atlantic Group B, rather than the opposite. This indicates a derived lineage in Western Atlantic, old enough to have had time to evolve forming the haplotype Group C. A more thorough sampling of the whole range of the Western Atlantic distribution and the inclusion of more data from Western Africa, are necessary before firm evidence can be obtained about the historical whereabouts of the main lineages of A. lixula.
Overall, the pattern of distribution of genetic variability (as shown in F st , Jost's D, MDS and AMOVA analyses) showed three groups of populations that differed significantly from each other (Brazilian, Eastern Atlantic and Mediterranean), while little structure could be found within these groups. It is remarkable that the F st measures based on sequence distance metrics and the differentiation measure D based on haplotype frequencies yielded essentially the same results. This is attributable to the prevalence of close haplotypes separated by small number of mutations (hence the low nucleotide diversity in general) that are widespread among populations. Thus, haplotype genetic differences had relatively little weight and most population structure derives from haplotype frequency differences.
Another striking pattern resulting from our molecular analyses is that recent demographic phenomena have shaped the present-day genetic structure of Arbacia lixula populations in the Eastern Atlantic and the Mediterranean. This does not seem to be the case of the Brazilian population but, given the small sample size, it is unclear if the resulting mismatch distribution (Fig. 6A) is either multimodal or L-shaped in this population. Multimodal curves are typical of populations at demographic equilibrium, but L-shaped distributions may result from very recent demographic bottlenecks [75]. More extensive sampling would be required to get the full picture of the demographic processes that have shaped the Brazilian populations of A. lixula.
The lack of an exclusively Mediterranean mitochondrial lineage of Arbacia lixula is remarkable. Other Atlanto-Mediterranean echinoderms such as Marthasterias glacialis [76], Holothuria mammata [77] or Paracentrotus lividus [73,78] do have lineages exclusive of the Mediterranean. These species have been probably present in the Mediterranean for several million years and their populations may have suffered several episodes of impaired gene flow during the Pleistocene glaciations. The genetic structure shown by A. lixula probably reflects a different demographic history from these other species.
Even if there is no phylogenetic break in the Mediterranean (as also found by Lessios et al. [26]) and alleles are widely shared at both sides of the Gibraltar boundary, this barrier seems nevertheless to restrict gene flow in Arbacia lixula, so as to establish significant differences in terms of haplotype frequencies between Mediterranean and Eastern Atlantic populations. The AMOVA (and pairwise comparisons) detected significant genetic differentiation between these groups of populations (Table 5), suggesting a reduced gene flow through the Strait of Gibraltar. Differently to what can be found in other marine organisms [79], the Strait itself, and not the Almeria-Oran Front (some 350 Km east of Gibraltar), is the place of the phylogeographic break, as the populations from the Alboran Sea are undistinguishable from other Mediterranean populations, but are significantly differentiated from most Atlantic populations (Fig. 4, Tables 3 and 7). Thus, A. lixula does not show any genetic differentiation among populations throughout the whole Mediterranean Sea. This could be due to recurrent gene flow, but oceanographic barriers such as the Almeria-Oran Front or the Siculo-Tunisian Strait [79] are strong enough to maintain genetic differentiation among different sub-basins in the case of other echinoderms of similar larval dispersive capacity [73,[77][78].  We favour the alternative explanation (for the lack of genetic structure) that the colonization of the Mediterranean by A. lixula is so recent (see below) that populations in the different Mediterranean sub-basins have not had yet enough time to diverge from each other.
In the case of Macaronesian and Cape Verdean populations (Table 6), it seems likely that the present-day genetic similarity could be the result of a recent demographic expansion (see below), which could have swamped any trace of previous differentiated lineages potentially formed during periods of restricted gene flow among archipelagos.
Brazilian populations of Arbacia lixula are completely differentiated from Eastern Atlantic and Mediterranean populations (Tables 3 & 4). In addition, they showed the lowest genetic diversity and did not show any signature of demographic expansion. Nevertheless, our sample size is small, and Northern and Central Brazilian populations of A. lixula have never been sampled for phylogeographic studies. More extensive sampling along the Brazilian coast would be required for a full understanding of factors shaping the genetic structure of the West Atlantic populations of A. lixula.
The almost complete lack of fossil record for Arbacia lixula in the Mediterranean is most revealing. At present, the species is highly abundant and occurs in areas that have been thoroughly sampled by palaeontologists. Other Mediterranean echinoids currently cooccurring in the same habitats are commonly found in assemblages of the Pleistocene and have been abundantly reported in the paleontological literature [80][81][82][83]. In contrast, only one fossil individual of A. lixula from the Mediterranean has ever been reported in the literature [13]. It was found in very young deposits from Livorno (Italy) whose recency led Stefanini to speculate that A. lixula had an exotic origin and had entered the Mediterranean in recent times [13]. A. lixula is consistently absent from fossil assemblages of the so-called ''Senegalese fauna'' that characterize the warmer periods from the Tyrrhenian stage (ca. 260-11.4 kya), which have been extensively sampled and thoroughly described [84][85][86][87][88][89].
As for the Atlantic archipelagos, recent work on the fossil echinoid fauna of Azores Islands [90] has revealed the presence of A. lixula, providing several tens of pieces of individuals, including the oldest known record of this species. These deposits are currently dated to 130-120 kya [91], which corresponds to the last interglacial or Riss-Würm (also called MIS 5e, ca. 130-114 kya). These specimens add up to the only other Atlantic A. lixula fossil specimen known from the Pleistocene of Madeira [13] whose dating is more uncertain.
Thus, there is scarce paleontological evidence of the occurrence of Arbacia lixula in the Mediterranean, and somewhat more, but still scarce, evidence of the colonization of the Atlantic archipelagos of Azores and Madeira, which probably occurred during the last interglacial period of the Pleistocene (MIS 5e). These observations are in agreement with the genetic signatures we observed in the mismatch distributions, which clearly show that recent sudden expansions have occurred in the Mediterranean and Macaronesian populations (Fig. 6). This is also supported by the strikingly star-shaped topologies of the haplotype network (Fig. 2) and by the comb-like clades in the BI phylogenetic tree (Fig. 3). Our temporal estimation for the demographic expansion in the Mediterranean (93.8-205.2 kya) is coherent with the only available fossil record [13]. This is considerably younger than the times for expansion events found in other Mediterranean echinoderms using the same estimation method, which vary from 300 to 600 kya [73,77] and fits with the possibility that the colonization of the Mediterranean by A. lixula took place as recently as during the last interglacial period (MIS 5e). This period was also the longest of all interglacial warm periods of the Pleistocene. The minimum winter surface temperature of the Mediterranean Sea stayed warmer than 19uC for several thousands of years [89]. This probably enabled tropical Atlantic populations of A. lixula to cross the Strait of Gibraltar and colonize the Mediterranean.  In the case of Eastern Atlantic populations, the exponential demographic expansion is even more apparent, since the mismatch distribution follows a sharp unimodal curve which fits to a sudden expansion model with a very high value for h 1 . This expansion probably occurred more recently than in the Mediterranean (31.3-68.4 kya). This estimation falls within the Late Pleistocene, an epoch generally dominated by the last glaciation (Würm), during which the mean sea level dropped down to 80 m below the present level [92][93]. Changes in ocean circulation related to this sea level drop can be related to the population expansion of A. lixula in the Eastern Atlantic. Contrary to what happens in the Mediterranean, the fossils available show that the species was present in Macaronesia before this expansion [90], so the demographic history of the Atlantic populations of A. lixula seems to be more complex than that of the Mediterranean populations. To complete the picture of the colonization of Atlantic archipelagos, data from continental African shores would be highly valuable.
An invasive species can be defined as a ''species that threatens the diversity or abundance of native species, the ecological stability of infested ecosystems, economic activities (e.g., agricultural, aquacultural, commercial, or recreational) dependent on these ecosystems and/or human health'' [94]. Although the term is generally applied to species introduced as a result of human activities, it should not be necessarily so. Moreover, ecosystem engineer species such as Arbacia lixula, that have shaped contemporary communities as the result of a colonization event that took place many years ago, can be falsely viewed as native [95]. According to our molecular data, A. lixula has indeed colonized the Mediterranean recently and complies with the terms of the former definition, even if it is usually viewed as native because its colonization took place following natural climatic changes, without human intervention.
Whether considered as an ''old natural invader'' or as native, the present trend of global warming can potentially boost the negative impact of A. lixula in Mediterranean ecosystems, thus possibly turning a ''natural'' colonization into an ecological problem related (at least partially) to human intervention. The ongoing warming [96] may facilitate population blooms of A. lixula in Northern Mediterranean, by releasing the constraint to larval development due to low water temperature. Warnings have been issued about its potential population increase and the generation of barren grounds in sublittoral habitats [10][11].
Thus, genetic data are in agreement with the consideration of Arbacia lixula as a thermophilous species that has recently colonised the Mediterranean and whose densities may increase in the foreseeable future. Monitoring of populations seems highly recommendable as a management tool in the near future for protecting the threatened Mediterranean shallow water ecosystems.

Supporting Information
Table S1 Haplotype frequencies of Arbacia lixula COI for all sampled localities. Haplotypes shared by two or more localities are represented in bold, while numbers not in bold correspond to private haplotypes. Background colours correspond to the three different haplogroups. Outgroups weights calculated by TCS are also displayed for each haplotype, and that with the highest outgroup weight (A2) is highlighted in green background. (XLS)