Genetic diversity and population structure of the rockpool shrimp Palaemon elegans based on microsatellites: evidence for a cryptic species and differentiation across the Atlantic–Mediterranean transition

The rockpool shrimp Palaemon elegans is an ecologically important crustacean species within the European coastline fauna. In the present study, genetic diversity and population structure and connectivity were assessed by examining 21 polymorphic microsatellite loci at 13 sampling sites located along the natural distribution range. All localities showed similar levels of genetic variability. Significant deficits of heterozygosity were recorded, most likely due to the presence of null alleles. Genetic structure analyses revealed two clearly genetically distinct groups within P. elegans but without following any geographical or oceanographic basis. Thus, our results provided nuclear evidence for the existence of a Mediterranean cryptic species within P. elegans, highlighting the need to revise its taxonomic status. Regarding P. elegans sensu stricto, population structuring was reported across the Atlantic–Mediterranean transition area, where the Almería-Orán Front restricts the gene flow between the Atlantic and the Mediterranean population. Moreover, while population connectivity was suggested between all Mediterranean localities, some substructure was found within the Atlantic group. Canary Islands exhibited a weak but significant genetic differentiation from all Atlantic mainland localities, consistent with the isolation-by-distance pattern detected throughout the Atlantic population. Overall, all these findings provided new insights into the population biology of P. elegans complex.

www.nature.com/scientificreports/ Sea nowadays. Similarly, P. elegans was introduced in the Baltic Sea, replacing the native counterpart P. adspersus in that region 5,6 . Palaemon elegans has been also recently reported to occur on the northeast coast of the United States 7 , so is currently distributed beyond European waters. Thus, even though P. elegans is considered an important species within the European coastline fauna due to its broad ecological niche and ongoing geographic expansion 8 , population genetics analyses are scarce. Genetic breaks occur when there is a marked genetic differentiation in a continuously distributed species 9 . For many years, high levels of gene flow associated with little genetic structuring were assumed for marine species with planktonic larval phases 10 as P. elegans. However, it has been proved that historical and oceanographic processes could drive genetic differentiation between marine populations, even in species with high dispersal abilities 11 . This is the case of the shrimp P. elegans since Reuschel et al. 8 carried out a phylogeographic analysis using two mitochondrial DNA (mtDNA) markers, revealing the existence of three haplogroups, one of them in the Atlantic localities (type I) and two from the Mediterranean localities (types II and III). In that study, a genetic break was found in the Atlantic-Mediterranean transition area, restricting the gene flow between these basins and leading to genetic differentiation between the Atlantic (type I) and the Mediterranean (type II) population. Similar phylogeographic breaks across the Atlantic-Mediterranean transition have been reported in a wide variety of species (see review in Patarnello et al. 11 ). These genetic discontinuities are frequently attributed to the effect of different oceanographic discontinuities acting in this area, i.e. the Almería-Orán Front, the Strait of Gibraltar, the Balearic Front or the Ibiza Channel 10 .
Back to P. elegans, genetic structure was not detected by Reuschel et al. 8 neither within the Atlantic population (type I) nor within the Mediterranean one (type II). As for P. elegans type III, it was only recorded from Mediterranean localities and showed high nucleotide divergence from the other two types 8 . This genetic pattern was confirmed in later phylogeographic studies in the Mediterranean Sea also using mtDNA markers 12,13 . Reuschel et al. 8 proposed that type III individuals could belong to a putative cryptic species existing within P. elegans. Cryptic species occur when two or more species are classified as a single nominal species because they are superficially morphologically indistinguishable 14 . The Mediterranean Sea is characterized by a high level of biodiversity and species richness with a proportion of endemism circa 25% 15 , largely related to its tortuous geological history. Consequently, cryptic speciation is not considered an uncommon phenomenon in the Mediterranean Sea and molecular markers have indicated the presence of cryptic species in several taxa including shrimps 16,17 . Given these findings in P. elegans, the need for further genetic studies using polymorphic nuclear markers, such as microsatellites, to clarify its phylogeography was highlighted.
Microsatellite markers have been used in a wide range of applications in population genetics, ecological, conservation and evolutionary studies 18 because of its codominant nature, biparental mode of inheritance and high levels of polymorphism 19 . Specifically in population genetics, microsatellites are extremely valuable tools as they might reveal the existence of genetically distinct populations even in fine-scale studies 20 . In fact, microsatellite loci have been crucial in understanding the population structure of several shrimp species as Penaeus vannamei 21 , Penaeus monodon 22 , Pandalus borealis 23 , Aristeus antennatus 24 or Palaemon serratus 25 . Thanks to the advent of high-throughput sequencing technologies, the development of microsatellite markers in non-model species has become a rapid and cost-effective process 26 . Precisely, we isolated and characterized 21 novel polymorphic microsatellite loci for Palaemon elegans using Illumina MiSeq sequencing 27 .
The aim of the present study was to address for the first time the genetic diversity and population structure and connectivity of the shrimp P. elegans along the European coastlines, with particular attention to the Atlantic-Mediterranean transition, using microsatellite loci previously developed 27 to shed light about the population biology of this species.

Results
Genetic diversity. A total of 400 individuals were genotyped for 21 microsatellite markers isolated and characterized in a previous work 27 . Thirteen coastal localities distributed along the natural distribution range of P. elegans were sampled. The number of individuals collected from each sampling site varied from 10 (Ebro Delta) to 50 (Collioure). No evidence of linkage disequilibrium between pairs of loci was detected across localities after the sequential Bonferroni correction. Hence, all microsatellites were considered as independent markers. Likewise, the Ewens-Watterson's neutrality test confirmed that all loci were selectively neutral. Estimators of genetic variability per locality, per locus and per locality and locus are available in Table 1, Supplementary  Tables S1 and S2, respectively. Across the microsatellite loci, the number of alleles ranged from 4 (loci Pe04, Pe13 and Pe20) to 25 (locus Pe06), with a grand mean of 5.3 alleles over localities and loci. Thirty-six out of the 217 identified alleles (16.59%) were private alleles. Every microsatellite locus showed private alleles from at least one sampling site (except loci Pe02, Pe03, Pe08 and Pe13 with no private alleles), and every locality showed private alleles in at least one locus (except Lanzarote) ( Table 1 and Supplementary Table S1). Average allelic richness (AR) across localities ranged from 2.032 (Santoña) to 2.453 (Livorno). The 21 microsatellites were highly polymorphic as the mean percentage of polymorphism was 98.53% (Table 1). The observed heterozygosity (Ho) over localities for each locus ranged from 0.134 (Pe09) to 0.842 (Pe14), and the expected heterozygosity (He) from 0.201 (Pe09) to 0.876 (Pe14) (Supplementary Table S1). Regarding the expected and observed heterozygosity for each sampling site, Ho ranged from 0.372 (Santoña) to 0.499 (Ebro Delta) and He from 0.468 (Santoña) to 0.607 (Livorno). Ho and He averaged 0.450 and 0.531, respectively (Table 1).
Genetic structure. The assessment of the modal value of DeltaK by the STRU CTU RE analysis (Supplementary Figure S1a) indicated that the most likely partition of the microsatellite dataset (400 sampled individuals) was two genetically distinct groups (K = 2, Fig. 1a). Most of the individuals (327 individuals) were assigned to one cluster, whereas the second cluster grouped 69 individuals. Accordingly, only 4 out the 400 sampled individuals could not be assigned to one specific cluster (probability of membership to a cluster < 0.75), most likely because they presented missing data at several loci. This partition in two groups did not follow any geographical or oceanographic pattern because the second cluster was formed by only some of the individuals collected in Granada, Almería, Collioure, Livorno and Mallorca, being the remaining individuals from these localities located in the main cluster. Factorial correspondence analysis (FCA) also showed genetic structuring in two main groups (Supplementary Figure S2) that coincided and supported the clustering inferred by STRU CTU RE. Taking into account that the smaller cluster was exclusively formed by Mediterranean individuals, coupled with the findings reported by Reuschel et al. 8 , we considered that this group is the putative cryptic species suggested to exist within P. elegans with mtDNA markers. Consequently, to check if there was such agreement between our nuclear markers and the previously analysed mtDNA markers with respect to the 69 individuals of the smaller cluster, we generated the cytochrome c oxidase subunit I (COI) sequence from them following the PCR conditions described in Reuschel et al. 8 . A COI fragment was successfully amplified for 64 out of the 69 individuals. BLASTn verified that the COI sequence of all these individuals corresponded to the COI sequence of those individuals named P. elegans type III (HE573177) and considered as a cryptic species in Reuschel et al. 8 . Therefore, the assignment of 69 individuals of our dataset to a putative cryptic species by the 20 microsatellite loci was also supported with a mtDNA marker. COI sequences from the individuals belonging to the putative cryptic species can be easily accessed in Supplementary File S1. In detail, distribution of the putative cryptic In accordance with the STRU CTU RE analysis, a peak of ΔK and a plateau in the plot of K versus ln Pr(X|K) (Supplementary Figure S1b) were also observed at K = 3, supporting a second partition of the data in three groups (Fig. 1b). One of these groups matched with the cluster corresponding to the putative cryptic species defined in K = 2. However, in this analysis, the individuals considered as P. elegans sensu stricto were divided into two clusters: individuals from Atlantic localities (Ré Island, Santoña, Ártabro Gulf, Cádiz, Tenerife and Lanzarote) When all sampling sites were analysed separately, a hierarchical locus-by-locus AMOVA (Table 2a) revealed that the proportion of global genetic variance was mostly attributed to differences among individuals within localities (83.60%), whereas the proportion of genetic variance among sampling sites was 16.40% (FST = 0.164). More AMOVA analyses were performed to test the grouping scenarios inferred by STRU CTU RE. Granada, Almería, Collioure, Livorno and Mallorca localities, with individuals belonging to the putative cryptic species, were subdivided by samples of each species (i.e. P. elegans sensu stricto and cryptic species) to conduct these analyses. Two-group AMOVA (P. elegans sensu stricto vs. cryptic species; Table 2b) maximized the FCT, revealing that the 26.23% of the total variance was explained by differences between P. elegans sensu stricto and the   Table 2c). Nevertheless, a genetic structuring of the whole genotype dataset in two groups yielded a higher FCT value.
Pairwise FST values were estimated between localities in which at least 10 individuals were sampled. Therefore, separating samples of Granada, Almería, Collioure, Livorno and Mallorca by species, FST values could not be calculated for some of these subdivided localities. Even so, pairwise FST values (Supplementary Table S3) suggested a marked and significant genetic differentiation in all the comparisons between localities of P. elegans sensu stricto and localities of cryptic species (FST ≥ 0.190), indicating the genetic divergence of these two species supported by the other analyses. Focussing on the population structure of P. elegans sensu stricto, pairwise FST values revealed significant genetic differentiation between the Atlantic and the Mediterranean localities (Table 3). Granada locality was connected with almost all the Atlantic localities and significantly differentiated from the remaining Mediterranean ones. Within the Mediterranean Sea, connectivity among Ebro Delta-Collioure-Marseille-Livorno was suggested by pairwise FST values. Along the Atlantic coast, connectivity was suggested among Ré Island-Santoña-Ártabro Gulf-Cádiz-Granada and some genetic differentiation was supported between each one of those and the Canarian localities. In other words, Tenerife and Lanzarote showed significant differentiation with all the localities regardless of the sampling source, while they were both connected to each other.
In order to increase the reliability of the study, the original microsatellite dataset (excluding locus Pe 11) was pared down to check the impact of the putative cryptic species on the genetic structure results of P. elegans sensu stricto. In fact, individuals belonging to the putative cryptic species presented a considerable amount of missing data in some loci, especially in loci Pe08 and Pe17. Thus, another STRU CTU RE run was carried out removing these 69 individuals and the 4 individuals that could not be assigned to any cluster, so the microsatellite dataset exclusively contained genotypes of individuals considered as P. elegans sensu stricto (327 individuals). DeltaK statistic distribution peaked at K = 2 (Supplementary Figure S1c), indicating again the partition of P. elegans sensu stricto samples into two clusters (Supplementary Figure S3a): one cluster joined Granada locality to the Atlantic localities (Ré Island, Santoña, Ártabro Gulf, Cádiz, Tenerife and Lanzarote), meanwhile the other cluster was formed by the Mediterranean localities (Almería, Ebro Delta, Collioure, Marseille, Livorno and Mallorca). The pattern of divergence between the sampling sites located on the Atlantic coast and those on the Mediterranean Sea was detected again for P. elegans sensu stricto. AMOVA was performed to test this two-group structure in this reduced dataset. When excluding the individuals of the putative cryptic species, differentiation among groups decreased (FCT = 0.094), explaining the 9.37% of the total variance (Supplementary Table S4).
By partitioning samples of P. elegans sensu stricto into two groups ( Fig. 1b and Supplementary Figure S3a) according to the Bayesian clustering analysis, most of the individuals from Tenerife and Lanzarote were assigned to the Atlantic Ocean cluster as expected, but a considerable number of Canarian individuals could not be assigned to either of the two clusters (probability of membership to a cluster < 0.75). Assessing the results of the STRU CTU RE analysis using the P. elegans sensu stricto dataset, a peak of ΔK and a plateau in the plot of K versus ln Pr(X|K) (Supplementary Figure S1d) were observed at K = 3, so the partition in three groups was plotted (Supplementary Figure S3b). With this grouping, the Mediterranean Sea cluster remained mainly invariable, whereas the Atlantic Ocean cluster was split into two subgroups, being then possible to clearly assign individuals from the Canary Islands localities to one of these two Atlantic clusters. Three-group AMOVA using this dataset unveiled that the 8.80% of the global variance is due to differences among groups, so the FCT value (FCT = 0.088) was lesser respecting to two-group AMOVA (Supplementary Table S4). Overall, the most likely genetic structure of P. elegans sensu stricto is two populations, one of them including the Atlantic localities and Granada, and the other encompassing the Mediterranean localities. Nevertheless, it is noteworthy to highlight that there is a substructure within the Atlantic population, with a subpopulation not exclusively but mainly formed by individuals from Table 3. Pairwise FST values (below the diagonal) between 11 sampling sites of Palaemon elegans sensu stricto (327 sampled individuals) as revealed by the microsatellite dataset (20 loci), and the corresponding P-values (above the diagonal). *Significance after FDR correction (P < 0.0033). **Significance after Bonferroni correction (P < 0.00091). www.nature.com/scientificreports/ the Canary Islands. This substructure is in line with the genetic differentiation between Tenerife and Lanzarote and the remaining Atlantic localities suggested by the FST pairwise values. Finally, isolation-by-distance (IBD) was retrieved in P. elegans sensu stricto by the Mantel test, as genetic and geographical distances resulted to be significantly correlated (R = 0.37, P = 0.001). Another two Mantel tests were performed to test IBD within each population of P. elegans sensu stricto. No IBD was detected within the Mediterranean basin (R = 0.00, P = 0.750), meanwhile a significant relationship between genetic and geographical distances was found along the Atlantic population (R = 0.40, P = 0.015).

Discussion
In this study, microsatellite markers were used for the analysis of the genetic diversity and population structure of the shrimp P. elegans. Levels of polymorphism were in line with those found in the counterpart common littoral shrimp P. serratus 25,29 . Sampling sites showed similar levels of genetic diversity (Na, A, AR, Ho or He). Thirty-one locality-locus combinations significantly deviated from HWE, mainly caused by heterozygote deficiency (29 out of the 31 cases) reflected in positive FIS values. Stutter-related scoring errors were discarded as the source of HWE departures because markers with ambiguous interpretation patterns were not included in the microsatellite set during the characterization work 27 30,31 . The 12 loci that deviated from HWE due to heterozygote deficit showed evidence of occurrence of null alleles, even though mainly in a low proportion. Only one locus, the locus Pe11, showed certain consistency in high null allele frequency estimates across localities and overall. For this reason, we decided to exclude locus Pe11 in population structure analyses as it was advised in González-Castellano et al. 27 .
Microsatellite data provided evidence for the existence of two clearly genetically distinct groups within P. elegans. Some individuals from Mediterranean localities showed a marked genetic divergence from the remaining specimens without following any geographical or oceanographic pattern. Since all shrimps from Granada, Almería, Collioure, Livorno and Mallorca were captured simultaneously in the same rockpools and unambiguously morphologically identified as P. elegans afterwards, we concluded that the explanation for the divergent individuals is the presence of a cryptic and sympatric species existing within P. elegans. Previous studies using mtDNA markers reported that some individuals from the Mediterranean Sea exhibited divergent haplotypes respecting the rest of the sampled collection, being those named P. elegans type III and considered as a cryptic species 8,12,13 . In fact, the individuals assigned with our microsatellite set to a cryptic species were confirmed to correspond with those P. elegans type III by their mtDNA COI sequences. The alternative hypothesis to explain this divergence in nuclear markers among our samples would be strong genetic structuring among individuals who live sympatrically and without exhibiting evident intrinsic reproductive barriers, what seems extremely unlikely.
Cryptic species are distributed among taxa, but they are particularly prone in the marine environment 32 , being detected within various species using microsatellite markers [33][34][35] . Here, the 400 sampled shrimps were assigned to one of the two unveiled genetic groups, i.e. P. elegans sensu stricto or cryptic species. Only two individuals from Granada and two from Mallorca could not be assigned to one specific cluster, most likely because they presented missing data at several loci. Consequently, no hybrids between P. elegans sensu stricto and the cryptic species were found. Hence, unknown reproductive isolating mechanisms between these two species should be acting to explain the absence of genetic hybrids. The main question to be solved is the origin of the cryptic species. Given that it was confined in the Mediterranean Sea, it should be originated there. The Mediterranean basin has a very complex geological history that largely contributed to its high level of biodiversity and endemism 15 . Reuschel et al. 8 postulated that P. elegans type III specimens might have their origin during the Messinian Salinity Crisis in late Miocene (around 5.6 Myr ago), when the Strait of Gibraltar was closed, isolating the Mediterranean Sea from the Atlantic Ocean and determining the extensive desiccation of the Mediterranean basin 36 . According to Reuschel et al. 8 , isolated P. elegans type III ancestors had to survive in local refugia in the Mediterranean Sea and genetically diverge from the Atlantic population of P. elegans, giving rise to the cryptic species detected within P. elegans nowadays. With the re-opening of the Strait of Gibraltar, the Mediterranean basin was flooded again with Atlantic water, so Atlantic P. elegans individuals should be re-introduced, but they should no longer be able to reproduce with the cryptic species ancestors, since none hybrids are currently detected. Similarly, the Messinian Salinity Crisis was considered the origin of two sister crab species of the genus Carcinus 37 and of two sister brittle star species of the genus Ophiothrix 38 . Our analyses using microsatellite loci support the existence of two species within P. elegans, providing nuclear evidence for the divergence signatures previously detected at mitochondrial level. Future efforts should be directed to clarify the origin of the cryptic species by phylogenetic or phylogenomic studies with precise calibration and to revise the taxonomy of P. elegans complex.
Regarding the population biology of P. elegans sensu stricto (P. elegans hereafter), a genetic structure was revealed in the whole study. Genetic differentiation analyses (Bayesian analysis, AMOVA and FST statistics) showed a pattern of two main populations: one composed by the Atlantic localities and Granada, and another composed by the Mediterranean localities. In other words, a clear geographic distinction between Atlantic and Scientific RepoRtS | (2020) 10:10784 | https://doi.org/10.1038/s41598-020-67824-7 www.nature.com/scientificreports/ Mediterranean sampling sites was detected, with an extension of the Atlantic genotypes into Mediterranean Sea waters. Granada locality was genetically closer to the Atlantic localities than to the other Mediterranean localities. Therefore, a biogeographical barrier restricting the gene flow between these two populations should be located between Granada (Atlantic population) and Almería (Mediterranean population). Precisely, this genetic break matched with the location of the Almería-Orán Front, a semi-permanent dynamic oceanographic front consisted in two anticyclonic gyres that entails abrupt changes of temperature and salinity 39 . The Almería-Orán Front has been considered as the most important hydrographic boundary to gene flow between the Atlantic Ocean and Mediterranean Sea surface waters 11 . Indeed, it was confirmed that the Almería-Orán Front acts on the population structure of several species as in P. elegans, e.g. the cuttlefish Sepia officinalis 40 , some fish species 41 , the barnacle Chthamalus montagui 42 or the wedge clam Donax trunculus 43 among others (see review in Patarnello et al. 11 ). Notwithstanding, it is known that population connectivity across the Atlantic-Mediterranean transition is reduced in other species by different oceanographic discontinuities, namely the Strait of Gibraltar 24,44 , the Balearic Front and the Ibiza Channel 45,46 . Most recently, in the congeneric common littoral shrimp Palaemon serratus, both mtDNA 47 and nuclear 25 markers have reported strong genetic differentiation among Atlantic and Mediterranean localities with an unusual phylogeographical break located west of the Strait of Gibraltar in the Gulf of Cádiz.
In the present study, microsatellite data unravelled that the Almería-Orán Front was the only oceanographic discontinuity affecting the population structure of P. elegans across the Atlantic-Mediterranean transition, so our results agreed with the phylogeographic pattern described in Reuschel et al. 8 using mitochondrial data. However, it is important to note that the Almería-Orán Front experiences seasonal fluctuations in intensity, size and position 39 that might mitigate its effects as oceanographic barrier between the populations of P. elegans.
Differentiation was more complex within each group. Within the Mediterranean population, all localities showed no genetic differentiation among them in the light of our data. Our shrimp samples only covered the northern of the Western Mediterranean basin, so it would be interesting to analyse more Mediterranean localities to study population connectivity throughout the entire Mediterranean Sea. As for the Atlantic group, genetic analyses indicated some substructure within this population. Specifically, Canary Islands localities (Tenerife and Lanzarote) showed low but significant genetic differentiation from the remaining Atlantic localities (Ré Island, Santoña, Ártabro Gulf, Cádiz and Granada) as suggested by the pairwise FST values and also observed in the Bayesian analysis. Thus, the above results allowed us to report two subgroups within the Atlantic population: one formed by the South of France and the Iberian Peninsula localities; and another group comprising the Canary Islands localities. The lack of connectivity between these subgroups is consisted with the isolation-by-distance pattern observed along the Atlantic population of P. elegans. In this regard, lack of gene flow between Macaronesian Islands (Azores, Madeira, Canary Islands and Cape Verde) and continental locations has been also reported for various marine invertebrates, such as limpets of the genus Patella 48 , starfishes of the genus Marthasterias spp 49,50 . or the sea cucumber Holothuria (Holothuria) mammata 51 , as well as for fishes as the triplefin Tripterygion delaisi 52 and the white seabream Diplodus sargus 53 . The remoteness of the Canary Islands from the Atlantic mainland sampling sites may explain the weak differentiation of Canarian P. elegans samples. The duration of the larval stage compared to the great distance that separates the Canary Islands and the continent may not be enough to prevent this differentiation. Related to this, present-day surface currents in Macaronesia, mainly the Azores Current and the Canary Current 54 , could also favour the relative isolation of the Canary Islands by disrupting larval dispersion. A more extensive sampling would be needed to reveal whether there is more genetic structuring along the Atlantic distribution of this shrimp.
Finally, the microsatellite loci used in this study were isolated and characterized from Atlantic individuals 27 , so they are specific markers for P. elegans sensu stricto. In order to study the genetic structure of the cryptic species found in the Mediterranean Sea, firstly we recommend revising its status as an independent new species, and then to develop novel specific microsatellite markers to carry out population analyses with them. De novo development of species-specific microsatellites is strongly recommended. Cross-amplification from congeneric species is not generally feasible because inherent problems like allele size homoplasy, polymorphisms biases, null alleles presence, broken repeats motifs or amplification of non-orthologous loci could arise [55][56][57][58] .
In conclusion, this study encompasses the first assessment of the genetic diversity and population structure of the shrimp P. elegans using microsatellite markers. Two clearly genetically distinct groups were detected within what is taxonomically considered as P. elegans, indicating the existence of a Mediterranean cryptic species within P. elegans. Furthermore, P. elegans sensu stricto is genetically structured in two populations: the Atlantic population and the Mediterranean population. The genetic break that separates these populations was found across the Atlantic-Mediterranean transition, matching with the location of the Almería-Orán Front. The microsatellite loci were sensitive enough to also reveal certain differentiation within the Atlantic population mainly involving individuals from the Canary Islands, which may be explained by an isolation-by-distance pattern.

Material and methods
Specimens of P. elegans used in this study were collected inshore between 2012 and 2019 from 13 European sampling sites: Ré Island (France), Santoña (Spain), Ártabro Gulf (Spain), Cádiz (Spain), Tenerife (Spain), Lanzarote (Spain), Granada (Spain), Almería (Spain), Ebro Delta (Spain), Collioure (France), Marseille (France), Livorno (Italy) and Mallorca (Spain) (Fig. 2). Animals were captured in intertidal rockpools with a fish trap and ethanol-stored at laboratory. Morphological identification of the 400 sampled shrimps was thoroughly carried out according to González-Ortegón & Cuesta 1 . Genomic DNA was extracted from 25 mg of abdominal muscle tissue using the NZY Tissue gDNA Isolation kit (NZYTech) following the manufacturer's instructions. Extracted DNA quality and concentration were determined with the NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific).  27 for P. elegans and following the wet-lab protocol described there. Fluorescently labelled polymerase chain reaction (PCR) products were run on a 3130xl Genetic Analysis System (Applied Biosystems) for fragment analysis in the Scientific Research Support Services (University of A Coruña), using GeneScan 500 (− 250) ROX internal size standard (Applied Biosystems). Geneious 8.0.5 59 was used for fluorescent profiles analysing and allele peaks calling.
Genetic variability was assessed with standard descriptive statistics. Number of alleles (Na), mean number of alleles (A), percentage of polymorphism (P), observed heterozygosity (Ho), expected heterozygosity (He) and deviations from Hardy-Weinberg equilibrium (HWE) for each locus and/or locality were estimated using GenAlEx 6.5 60 . Significance levels were adjusted applying the sequential Bonferroni correction whenever multiple tests were performed. Allelic richness per locality was calculated using FSTAT 2.9.4 61 . GENEPOP 4.2 62 was used to estimate inbreeding coefficient (FIS) and to test linkage disequilibrium between pairs of loci using Fisher's exact test. Significance was computed by a Markov chain method using 10,000 dememorizations, 5,000 batches and 5,000 iterations per batch. Departures from selective neutrality were evaluated for each locus applying the Ewens-Watterson test using the algorithm given in Manly 63 and implemented in POPGENE 1.32 64 after 1,000 simulations. Null allele frequency was estimated using FreeNA 28 and the EM algorithm 65 with a number of bootstrap replicates fixed to 10,000.
Genetic structure was assessed using the Bayesian model-based clustering method implemented in STRU CTU RE 2.3.4 66 . The admixture model with correlated allele frequencies between populations 67 was used. A burn-in of 100,000 iterations followed by 300,000 iterations for parameter estimation was set up. Each simulation was run 10 times, exploring values for K ranging from 1 to N + 2, where N was the number of localities. Results were processed using STRU CTU RE HARVESTER 68 . The best partition of the data was determined by examining both the log probability of the data (ln Pr(X|K)) and the ΔK statistic following Evanno, Regnaut, & Goudet 69 . CLUMPP 1.1.2 70 was used to permute the admixture coefficients for the several independent runs resulting for the chosen K-value. STRU CTU RE PLOT 2.0 71 was employed to visualize the output from CLUMPP. Genetic discrimination among samples was additionally assessed through a factorial correspondence analysis (FCA) using GENETIX 4.05 72 .
Genetic differentiation was estimated as pairwise FST values (10,000 permutations) between localities using ARLEQUIN 3.5.1.2 73 . Significance level was adjusted applying sequential Bonferroni correction for multiple comparisons and false discovery rate (FDR) correction (B-Y method as described in Narum 74 ). Several hierarchical locus-by-locus analyses of molecular variance (AMOVA) were carried out using ARLEQUIN 3.5.1.2 73 in order to estimate partitioning of genetic variance among groups (FCT), among localities within groups (FSC) and within localities (FST). Individuals were pooled according to a geographic/oceanographic criterion or based on the results previously obtained with STRU CTU RE to test different groupings.
Lastly, the relationship between genetic distance and geographic distance was assessed with a Mantel test 75 to test isolation-by-distance (IBD) using the program ISOLDE implemented in GENEPOP 4.2 62 with 10,000 permutations. Pairwise genetic differentiation among localities was estimated using linearized FST [FST/(1 − FST)] to construct a genetic distance matrix. Pairwise geographic distance between localities was measured in kilometres following the coastline in Google Earth (Google Inc.) to create a geographic distance matrix.

Data availability
Accession codes: Sequences containing the microsatellite loci were deposited in GenBank under Accession Numbers MH078079-MH078099. Genotype dataset generated during the current study is available from the corresponding authors on reasonable request. www.nature.com/scientificreports/