Between 20 and 23 individuals of Arbacia lixula were sampled by scuba diving, from autumn 2014 to summer 2015, from 11 locations spanning the Mediterranean basin (for a total of 240 individuals, Table 1 and Figure 1). Western Mediterranean samples included three populations from the transitional zone of the Alboran Sea: Tarifa (TA), Torremuelle (TO), and La Herradura (LH), plus Carboneras (CA), Palos (PA), Xàbia (XA) and Colera (CO) from the remaining Iberian coastline and Cap Bon (CB) at the north-eastern margin of the Gulf of Tunis. Eastern Mediterranean localities were Murter (MU) and Otranto (OT) from the Adriatic Sea, and Dalyan (DA) from the eastern side of the Aegean Sea (Figure 1). This sampling spans ca. 3,100 km (shortest distance by sea) from TA to DA. The layout allows testing the influence of the main oceanographic barriers in the area (Pascual et al., 2017): The Almería-Oran front between Atlantic and Mediterranean, the Ibiza Channel separating North and South of the Iberian Peninsula, the Siculo-Tunisian divide between western and eastern Mediterranean, and the Otranto Strait separating the Adriatic Sea.

Sea urchins were dissected in situ to extract Aristotle’s lantern and were stored in absolute ethanol. Between 10 and 15 mg of muscle tissue of the Aristotle’s lanterns was used to extract the DNA using QIAamp^® DNA Mini Kit (QIAGEN) following manufacturer’s instructions.

Between 0.58 and 3.02 μg (mean ± SD = 1.96 ± 0.43) of genomic DNA (with ratio 260/280 ranging from 1.61–2, and ratio 260/230 ranging from 1.23–2.77) were used to perform a Genotyping by Sequencing protocol (GBS, Elshire et al., 2011) at the “Centre Nacional d’Analisis Genòmica” (CNAG-CRG, Barcelona, Spain). We selected the restriction enzyme ApekI based on digestion profile sizes of the genome of A. lixula (data not shown). The resulting fragments were ligated to individual barcodes and to a common adapter. A PCR was performed for fragment enrichment before paired-end sequencing of 2 bp × 100 bp fragments in an Illumina HiSeq 2000 platform using 96 multiplexed individuals in each lane of a flow cell. Each individual was sequenced in different lanes to obtain a mean number of reads of ca. 4 million.

The GIbPSs toolkit (Hapke and Thiele, 2016) was run to analyze the raw Illumina sequences. This program has already been used in previous studies of marine invertebrates (Casso et al., 2019; Carreras et al., 2020) and has proven its suitability to analyze GBS paired-end sequences of non-model organisms. Pre-processing of the sequences (which include trimming, filtering and assembling of paired-end sequences), locus identification (formation of the svars, definition of loci and alleles), and analysis of inferred loci (checking for indels and discarding deeply sequenced loci) were carried out following methods in Carreras et al. (2020).

In short, all raw sequences were trimmed to 80 bp. Filtering values were set to discard sequences shorter than 32 bp, with N’s or with a phred score quality lower than 22 in a sliding window of 5 bp. A minimum overlap of 5 bp was required to assemble forward and reverse sequences of a paired-end read. Duplicate sequences originated in loci shorter than the read length were identified and only the forward read was kept for further analyses (Carreras et al., 2020). Svars were identified for each individual separately and grouped into loci using the default distance settings, requiring at least five sequences for a locus to be deemed as valid. Note that all the polymorphic sites in each locus were used for identifying alleles. We thus used sequence variants (multiallelic), where each allele at a given locus is defined as the whole sequence with all polymorphic sites, not just one SNP per locus (biallelic). Multiallelic markers have been shown to provide more information and resolution power than biallelic markers (Ryman et al., 2006), and empirical tests on population genomics of different marine organisms have demonstrated this superior resolution on GBS data (Casso et al., 2019; Carreras et al., 2020). A global loci catalog was assembled with all the alleles found in all loci when combining all the individuals. Loci flagged by the program as potentially including an indel or with more than two alleles per individual were discarded, and deeply sequenced loci (those with a median depth percentile above 0.2, likely corresponding to paralogous regions with multiple copies) were likewise removed. A final locus selection was done to retain only loci present in at least 170 individuals (70% of the total dataset).

The “pegas” R package (Paradis, 2010) was used to evaluate Hardy-Weinberg Equilibrium (HWE) for each locus in each sampled site, and significant departures were identified. A Bonferroni correction was not applied since it dramatically increases the probability for type II errors (Cabin and Mitchell, 2000; Moran, 2003). Instead, we applied a Benjamini and Yekutieli (B-Y) false discovery rate (FDR) correction for multiple comparisons (Benjamini and Yekutieli, 2001) as detailed in White et al. (2019) at an overall corrected 0.05 α-level. Any locus with significant departures from HWE in seven or more sampling sites (>60%) was removed from the final dataset.

Pairwise genetic distances among sampling sites were assessed using the Weir and Cockerham (1984) estimator of F_ST using the “hierfstat” R package, with 999 permutations to obtain the respective p-values (Goudet, 2005) and significance thresholds were set after B-Y FDR correction as detailed above. The results were plotted with heatmaps and dendrograms using the “gplots” R package (Warnes et al., 2016). Genetic diversity values (observed and expected heterozygosity and allelic richness) and F_IS values per sampling site were calculated using the R packages “genetics” and “hierfstat” (Goudet, 2005; Warnes et al., 2019). The significance of the inbreeding coefficient was tested with function boot.ppfis and 999 permutations. If the generated values did not include 0, the corresponding F_IS was deemed significant.

A discriminant analysis of principal components (DAPC) was carried out with the “adegenet” R package (Jombart, 2008) using localities as prior groupings. We used the function “xvalDapc” to determine the number of PCs to retain and three discriminant functions using localities as prior grouping information. In order to ascertain whether there is genetic differentiation between the two major groups of sampled sites, western (TA, TO, LH, CA, PA, CB, XA, and CO) and eastern (OT, MU, and DA) Mediterranean (see section “Results”), an analysis of the molecular variance (AMOVA) was performed, grouping populations in these regions, with the program Arlequin version 3.5.2.2 (Excoffier and Lischer, 2010). Two further AMOVAs were performed with the western Mediterranean populations to test for differentiation between the Alboran Sea and the remaining localities, and between the two northernmost western localities (XA and CO, see section “Results”) and the remaining populations.

We also performed a STRUCTURE analysis, which implements a Bayesian clustering method to identify the most likely number of genetic groups (K) even in the presence of linked loci (Falush et al., 2003). Following Evanno et al. (2005), we carried out 10 runs per each value of K ranging from 1 to 11. We used the model of correlated allele frequencies and a burn-in of 50,000 steps followed by 500,000 Markov Chain Monte Carlo steps. We estimated the statistic ΔK to infer the most likely number of groups using STRUCTURE HARVESTER (Earl and vonHoldt, 2012). The 10 runs of STRUCTURE for the most probable K were averaged using CLUMPP version 1.1.2 (Jakobsson and Rosenberg, 2007).

A Mantel test was done to test isolation by distance in our dataset using the “vegan” R package (Dixon, 2003; Oksanen et al., 2019) considering the shortest distance by sea among sampling localities, measured with Google Earth. In order to avoid the influence of any possible inter-basin genetic structure, the Mantel test was repeated considering only sampling sites of the western Mediterranean (TA, TO, LH, CA, PA, CB, XA, and CO).

A network analysis was done using the F_ST matrix with the program EDENetworks (Kivelä et al., 2015). We used the F_ST values among populations as the dissimilarity matrix and then transformed them into a network. We chose the automatic thresholding option, that sets the maximal distance values required to keep an edge just above the percolation threshold (at which the network breaks down into its main components). Finally, geographic coordinates of sampling populations were introduced to adjust the network to the actual geographic location.

Both monthly mean sea surface temperature and salinity data were obtained from Copernicus Marine Environment Monitoring Service database (CMEMS¹, product identifier: MULTIOBS_GLO_PHY_REP_015_002) using the R package “ncdf4” version 1.16 (Pierce, 2017) from January 1987 to December 2015 and averaged.

We used the envfit function of the vegan R package (Oksanen et al., 2019) to fit the environmental variables into the DAPC ordinations as vectors, whose direction indicated the main gradient of change in the environmental variable and whose length was scaled to reflect the correlation between the spatial configuration and the environmental variable. The significance of these correlations was assessed with 999 permutations of environmental variables and using the squared correlation coefficient (r²) as a goodness of fit statistic.

In addition, for each abiotic variable, we computed distance matrices using the Euclidean distances for the mean, maximum, minimum and range of the variable. These distance matrices were compared individually with the F_ST matrices using partial Mantel tests as implemented in “vegan.” The geographic distance matrix was used as a controlling variable, as both temperature and salinity have an east-west pattern of variation in the Mediterranean (Coll et al., 2010) that can co-vary with geographic distance among localities (as it is mostly due to the longitudinal span of the samples).

Candidate outlier loci were identified using two different programs: Arlequin (Excoffier and Lischer, 2010) and BayeScan (Foll and Gaggiotti, 2008). The analysis with Arlequin was performed using 100,000 coalescent simulations and 100 demes in a finite Island model. Benjamini-Yekutieli FDR corrections were applied to the p-values as mentioned above. BayeScan search was implemented with a total of 150,000 iterations, with a burn-in period of 100,000 and 40 pilot runs of 10,000 iterations each. Outlier loci in BayeScan were identified using a qval of 0.05. In order to avoid false positives only loci detected with the two programs were kept as candidate outliers.

Two different sample groupings were used to detect outlier loci: “All populations,” in which all samples were included and grouped by sampling site to obtain a general overview of the outlier loci across the sampled geographical range; and “West vs. East” grouping, in which all samples were included in two groups according to their geographical distribution, western or eastern Mediterranean. We created two datasets using the results of these analyses. First, we defined as the “outlier” dataset the one comprising the outlier loci found by any of the two groupings (“All populations” or “West-East”). Second, we defined a “neutral” dataset in a conservative way by removing from the general dataset the loci of the “outlier” dataset plus any outlier found in any of the two groupings by any of the two programs using a conservative α < 0.05 for the p-value provided by Arlequin and the q-value provided by BayeScan. In this way we made sure that no potentially selected locus was included in the “neutral” dataset.

The analyses of population structure carried out with the complete loci dataset (DAPC, F_ST, comparisons with environmental variables) were repeated using only the “outlier loci” and the “neutral loci” datasets.

We used a Redundancy Analysis (RDA) (Forester et al., 2018) in order to identify candidate loci significantly associated to environmental variables (association loci). As this type of analysis needs a SNP biallelic dataset, we selected from our multiallelic variant dataset only the first polymorphic position (SNP) of each locus, using the option –uS 2 of GIbPSs (Hapke and Thiele, 2016). Given the restriction of the analysis that does not allow missing data, we imputed missing genotypes by replacing them with the most common genotype across all individuals. Finally, we identified significantly associated loci with the “rda” function of the R package “vegan” (Oksanen et al., 2019) using as response variable the matrix of SNP genotypes of all loci and the environmental variables as predictors.

A mean of 4,043,922 reads per individual were obtained. A dataset of 5,241 loci was retained after all filtering steps and used for the analyses. A genepop file with the allele composition of all individuals is available in Supplementary File 1. A list of all alleles with their sequences is presented in Supplementary File 2. In order to link loci names in both files an equivalence table is provided as Supplementary File 3. The raw sequence data was uploaded to the SRA repository (Bioproject PRJNA746276).

The length of the selected loci ranged from 35 bp to a maximum of 152 bp (mean ± SD = 124.44 ± 35.26 bp). Thus, there was a wide range of alleles per locus (from 2 to 190; mean ± SD = 25.77 ± 21.59), with 27 loci having only two alleles and 68 having more than 100 alleles (Supplementary Figure 1).

Table 1 shows the main features of the sampled populations. The mean number of alleles per locus per sampling site was similar (mean ± SD = 7.20 ± 0.14), ranging from 6.84 (Colera) to 7.33 (Xabia). The number of private alleles per locus (mean ± SD = 1.20 ± 0.048) ranged from 1.127 (Colera) to 1.275 (Murter).

There was a deficit of heterozygotes at all sampling sites (Ho = 0.379 ± 0.002 and He = 0.441 ± 0.002, mean ± SD). F_IS values were in all cases positive and low (mean ± SD = 0.141 ± 0.006), being nevertheless significant at all sites.

Pairwise F_ST values (Supplementary Table 1) were very low in general with values near zero or even negative (mean = 0.00228). Nonetheless, comparisons among western (TA, TO, LH, CA, PA, CB, XA, and CO) and eastern (OT, MU, and DA) Mediterranean sampling sites showed the highest genetic distances (mean = 0.00464, Supplementary Table 1 and Figure 2), and were all significant. On the other hand, pairwise F_ST comparisons within the two major regions had lower values (among western localities, mean = 0.0005; among eastern localities, mean = −0.0001) and were mostly not significant, albeit XA and CO, the two northernmost sites of the western Mediterranean, showed slightly higher genetic distances with the rest of the western Mediterranean sites (mean = 0.0012). The transitional localities found in the Gibraltar area and the Alboran Sea (TA, TO, and LH) did not show any significant differentiation with the remaining western Mediterranean populations, except for the comparison XA-TO. Likewise, the heatmap and the dendrogram based on F_ST population pairwise distances showed two major groups, western and eastern Mediterranean, with XA and CO forming a clade slightly separated from the rest of the western sampling sites (Figure 2).

The two major groups are also evident in DAPC results (Figure 3). Western sites appear clustered in one extreme and the eastern localities in the other extreme of the first axis. Within the western group, the Alboran Sea populations appear concentrated in one extreme of the first axis and XA and CO are placed on the other extreme, while CB (the easternmost population in the western basin) occupied a central position in the western Mediterranean cluster. Further, AMOVA results (Supplementary Table 2) showed significant genetic differences between the western and eastern groups of populations, although the percentage of variation explained was low (0.41%, p < 0.001). Likewise, an AMOVA of the western populations comparing the Alboran Sea with the remaining localities showed no significant differentiation between groups (p = 0.085), while a comparison between XA and CO with the other western Mediterranean populations explained a small but significant percentage of variation (0.06%, p < 0.001). In all AMOVAs, no differentiation was found among populations within groups, and most variance was concentrated within individuals (Supplementary Table 2).

The change in ΔK of the STRUCTURE analyses (Supplementary Figure 2) pointed to the existence of three genetic clusters in our dataset. The plot of the assignment probabilities (Figure 4) showed that one of the clusters was mostly restricted to the eastern Mediterranean, with some presence in the northern populations of the western basin. The other two groups were present throughout the western Mediterranean.

The Mantel test between geographic and genetic distances was highly significant (r = 0.816, p = 0.003, Supplementary Figure 3). This effect, however, was due to the separation between eastern and western basin populations, which had both wide intervening distances and genetic differentiation. When the analysis was repeated with only the western populations, no relationship was found (r = −0.144, p = 0.260).

The results of the network analysis, superimposed to the geographic layout (Figure 5) again reinforced the existence of two main groups of localities, linked by a weak edge between XA and OT. Accordingly, these two localities had the highest betweenness centrality, which reflects the importance of a node in connecting other nodes of the network.

The environmental vectors plotted in the DAPC ordination in Figure 3 show that the mean temperature had the shortest vector length, while mean, minimal, and maximal salinity were the longest vectors, indicating higher correlation of salinity variables with the obtained ordination. The fit of the environmental variables to the ordination scores is presented in Table 2. All variables considered had a significant fit with the ordination based on genetic distances. However, the r² coefficients were much higher for salinity measures (above 0.78 for mean, maximal, and minimal salinity).

As expected, when using partial Mantel tests to compare distance matrices correcting for the effect of geographic distance (Table 2) the correlations between genetic and environmental distance matrices were lower, but they were still significant for the mean and maximal salinity measures, thus confirming that salinity had an effect on the genetic make-up of the populations, over and above the effects of geographic distance per se. No significant effect of temperature variables was detected after controlling for distance.

A total of 25 loci were identified as outliers (putatively under selection): 13 loci were found when considering all populations, these same loci plus another 12 were found when joining populations in western and eastern groups. Additionally, we ended up with a set of 4,540 neutral loci after removing 701 loci, including the 25 outlier loci mentioned above and 676 additional loci found by any of the outlier detection methods applied with a non-corrected threshold of 0.05. Hereafter we will call them the “neutral” loci and the “outlier” loci datasets. Supplementary File 3 provides the IDs of the loci belonging to both datasets. The F_ST values among populations obtained using these datasets are presented in Supplementary Table 3. As expected, the mean F_ST was an order of magnitude higher when using outlier loci than when using the neutral dataset (0.0669 vs. 0.0022), although the pairwise comparisons found to be significant with the neutral and outlier loci are mostly the same as those found with all the loci (involving basically comparisons of the eastern Mediterranean with the other localities).

The population structure revealed by the neutral and outlier loci was similar to that obtained with the complete loci dataset, as shown by heatmaps (Supplementary Figure 4) and DAPCs (Figure 6), with a major break between the populations of the west and east of the Mediterranean (with XA and CO appearing as a separate group in the western Mediterranean cluster, Supplementary Figure 4). Both ordinations showed a similar pattern of correlation with the environmental variables as the complete dataset. The envfit analyses revealed a significant fit of all environmental vectors with both loci datasets (Supplementary Table 4). As before, the r² values were much higher for salinity (with the exception of salinity range) than for temperature. When partial Mantel tests were run to filter out the effect of geographic distance on the correlation between genetic and environmental distance matrices, only salinity (mean and maximal) showed a significant correlation with the genetic differentiation (F_ST) measures (Supplementary Table 4), as already detected with the complete dataset (Table 2).

The first two axes of the RDA including all the samples explained an accumulated 38.43% of the genetic variability found (Figure 7). The first axis explained most of this variability and sorted the samples along a longitudinal gradient related mainly to salinity-related variables, with eastern localities at one extreme and western localities at the other. The second axis was related to temperature variables. A total of 243 loci were found to be significantly associated to the environmental variables tested, being salinity variables the ones with the higher number of associated loci (Supplementary Figure 5). Eight of these loci were also detected as F_ST outliers in the previous analyses.