Samples and data analyzed in this study originated from previous efforts carried out in the framework of the EU (European Union) FishPopTrace Project (Martinsohn et al., 2009). Population samples were collected over the period from 1999 to 2009 from 13 geographical sites in the Mediterranean Sea ( Figure 1 and Table 1 ) as detailed in Nielsen et al. (2012). Genotype data at 426 SNP markers were obtained with the GoldenGate ™ (Illumina) high-throughput genotyping assay as described in Helyar et al. (2012) and Milano et al. (2011) and in the Supplementary Materials 1.1.1 . Different filtering criteria were applied to the Mediterranean samples, given their lower call rate compared to the Atlantic samples (Diopere et al., 2018). Individuals and loci with more than 20% and 10% of missing data, respectively, were removed from the dataset. Monomorphic markers in the Mediterranean dataset were also discarded prior to the analysis. Details of the evaluation of missing data and filtering procedure are presented in the Supplementary Materials 1.1.2 . The final Mediterranean dataset contained 352 individual genotypes at 380 SNPs.

Otolith data were also obtained for population samples from five geographical locations, maximizing the geographical coverage over the Mediterranean Sea ( Figure 1 and Table 1 ), full details of otolith processing and analysis are described in the Supplementary Materials 1.2 . In these five sites, the same individuals were processed for both molecular and otolith data. Based on geographical position of the 13 sampling sites, we identified nine centroid locations corresponding to a coastal slice defined as the portion of each GSA delimited by the continental shelf boundaries within the depth of 200 meters (m) ( Figure 1 ). The sampling sites were associated to the nearest centroids for the seascape analyses. We used the geographical coordinates (latitude and longitude) of the centroid for all analyses ( Figure 1 ). R version 3.6.3 was used for all analyses performed within the R environment in this study (R Core Team, 2022).

To assess genetic variation between samples, observed (Ho) and expected (He) heterozygosity, polymorphic rate, and F_IS estimate (Weir and Cockerham, 1984) were calculated using Genepop v.4.7.5 R packages (Rousset, 2008), respectively ( Supplementary Materials 1.1.3 ). Global and by locus Hardy-Weinberg (HW) tests were performed using R version of Genepop v.4.7.5 ( Supplementary Materials 1.1.3 ). We estimated global and pairwise values of F_ST index (Weir and Cockerham, 1984) and relative significance through 10,000 permutations to assess genetic differentiation among populations using strataG R package (Archer et al., 2017) ( Supplementary Materials 1.1.4 ). We then created a heatmap based on F_ST pairwise comparison to inspect population structure. A Discriminant Analysis of Principal Components (DAPC) was performed with adegenet R package (Jombart, 2008) to identify the individual-based population structure using both the complete and Mediterranean datasets. We ran the method both including as input the optimal number of clusters determined based on the Bayesian Information Criterion (BIC), as well as providing the sampling sites as prior. To describe the genetic structure, the first option allowed us to identify the optimal number of clusters, while the second one is used to assess the quality of the discrimination ( Supplementary Materials 1.1.5 ). Outlier loci were identified using four independent methods: (1) Fast Principal Component Analysis based on the pcadapt R package (Luu et al., 2017), (2) Fst-heterozygosity distribution analysis using the fsthet R package (Flanagan and Jones, 2017), (3) Bayesian analysis using Bayescan version 2.1 (Foll and Gaggiotti, 2008), and (4) the coalescent approach (FDIST) implemented in Arlequin v.3.5.2.2 (Excoffier and Lischer, 2010). Throughout, SNPs with minor allele frequency (MAF) lower than 0.05 were removed from the data set before applying the methods (Bekkevold et al., 2020), FDR was set at 5% and 95% CI was utilized to identify loci with fsthet R package (Flanagan and Jones, 2017). The method implemented in pcadapt R package allows the identification of genetic markers potentially involved in biological adaptation considering their relationship with the population structure, a detailed description of the method is included in Supplementary Materials 1.1.6 . We used the FDR correction method transforming the p-values into q-values and we choose α = 0.1 as the threshold for the selection of outlier SNPs (Benjamini and Hochberg, 1995). We used fsthet R package (Flanagan and Jones, 2017) to identify loci with extreme F_ST values relative to their heterozygosity H_t . Smoothed quantiles from empirical F_ST -H_t distribution were calculated. Loci lying outside of the quantiles for the 95% confidence level were regarded as outliers. Bayescan relies on differences in allele frequencies between subpopulations to identify candidate loci under selection. Subpopulation specific F_ST coefficients are divided into two components by logistic regression. The population-specific component β is shared by all loci, and the locus-specific component α is shared by all populations. Specifically, positive values of alpha indicate diversifying selection whereas negative values suggest balancing/purifying selection. For each locus, a posterior probability is computed ( Supplementary Materials 1.1.6 ) and SNPs with prior odds (PO) greater than a threshold are considered outliers. We set prior odds (PO) to 10 since the data set comprised less than 1,000 SNPs and we set the FDR to 5%. We applied FDIST approach implemented in Arlequin v.3.5.2.2 (Excoffier and Lischer, 2010) which uses simulations under an infinite island-model to generate the distribution of F_ST as a function of heterozygosity. SNPs found in the tails of the distribution are then identified as outliers. P-values were corrected through FDR approach (Benjamini and Hochberg, 1995), and SNPs with an associated p-value lower than 0.05 were considered as putative outliers. We then evaluated the overlap in loci using Venn diagram by visualizing SNPs in common to each method to create datasets of neutral and outlier SNPs that could suggest different patterns of genetic divergence. To create the outlier dataset, we considered the sum of all loci identified as outliers by each method, while the remaining markers constituted the neutral dataset. Genetic differentiation and multivariate analysis (PCA, DAPC) were repeated using the two resulting datasets to compare potential alternative patterns of divergence. To further investigate population structure, we ran STRUCTURE v.2.3.4 (Pritchard et al., 2000) using the datasets of neutral and outlier SNPs as separate inputs ( Supplementary Materials 1.1.7 ).

The variation in both otolith shape and trace element composition was investigated. Otolith contours were extracted from images and Fourier coefficients were calculated using Momocs R package (Bonhomme et al., 2014). The detailed procedures may be found in Supplementary Materials 1.2 . We calculated five shape indices: circularity, roundness, ellipticity, form factor and rectangularity (Tuset et al., 2003). In addition, the otolith shape was described through elliptic Fourier descriptors. We identified outliers with Interquartile Range (IQR) method to remove them before statistical analysis. Shape indices were tested for normality and homogeneity of residuals using Shapiro-Wilk Normality test and Levene’s test, respectively. Only normally distributed shape indices were retained for further analysis. To avoid redundancy of the data, Pearson correlations (Mérigot et al., 2007) were calculated using rcorr function in Hmisc R package (Harrell, 2019). We investigated whether the relationship between shape indices and total fish length differed significantly between sampling sites with analysis of covariance (ANCOVA) (Longmore et al., 2010; Keating et al., 2014; Rashidabadi et al., 2020) using the lm R function. If a significant interaction between sampling sites and fish Yc  =  Y  -  b  *  L where Y is the shape index, b is the common within group slope of the shape-size relationship and L is the measurement of total fish length (mm) (Keating et al., 2014). We used as slope the coefficient b resulting from the linear regression between total fish length and shape index. ANCOVA analysis was repeated to check that any co-variations were effectively treated (p > 0.01). We then performed a PERMANOVA analysis to uncover differences between geographical samples using adonis function in vegan R package (Oksanen et al., 2022). For the analysis, we computed a dissimilarity matrix using vegdist function with Mahalanobis method (Rashidabadi et al., 2020) and 1,000 permutations. The homogeneity of group dispersion was investigated with betadisper function. Then, we applied pairwise.adonis function in pairwiseAdonis R package to identify significantly different pairs of geographical sites, applying the Benjamini-Hochberg correction for multiple testing.

The first 40 elliptical Fourier harmonics were obtained from each otolith. Each harmonic n is described by 4 coefficients a_n , b_n , c_n and d_n , a total of 160 elliptical Fourier descriptors (EFDs) were available. We calculated the Fourier power (FP) for each otolith utilizing the calibrate_harmonicpower_efourier to estimate the number of harmonics required for the elliptical Fourier descriptors analysis. Using the efourier function, we computed EFDs, and we obtained a collection of coordinates (Coe object). Through a PERMANOVA analysis on the EFDs, we assessed differences between geographical sites using adonis function in vegan R package (Oksanen et al., 2022). For PERMANOVA analysis we followed the same procedure applied to the shape indices. Canonical analysis of principal coordinates (CAP) was performed on both types of morphometric data to visualize the spatial differences in classification accuracy obtained by leave-one-out cross-validation. We used CAPdiscrim function in BiodiversityR R package to perform the analysis based on Euclidean distances. To determinate the number of PCoA axes (m) to retain in the analysis we follow the general procedure used by Kindt and Kindt (2015). The argument m was set from 1 to 100 in the CAP analysis to explore the difference in the proportion of correct assignment of individuals to groups. Thus, we selected the m with the highest percentage, and we used it for the CAP analysis. After the PCoA, we checked for the percentage of total variance in the dissimilarity explained by the m PCoA axes. We checked that the axes chosen did not exceed 100% of the variance. Furthermore, the SIs and the EFDs were combined to evaluate their discriminatory power. We performed canonical analysis of principal coordinates (CAP) following the same procedure described above.

Trace elements composition data were obtained from core region and edge of sole otoliths using laser ablation inductively coupled mass spectrometry (LA-ICPMS) (Chang et al., 2012; Morales-Nin et al., 2014). After assessment of missing values, we replaced concentration of elements for which more than 25% of values were below the LOD with respective LOD value (Mercier et al., 2011). The final data set for both edge and core comprised ²³ Na, ²⁴Mg, ⁶³Cu, ⁶⁴Zn, ⁸⁵Rb, ⁸⁸Sr, ¹³⁸Ba, ²⁰⁸Pb. Trace element concentrations were log-transformed for further analysis. The same statistical analysis was performed for both core and edge elemental data. Before data analysis, we checked for normality and homogeneity of the residuals with Shapiro-Wilk Normality test (Shapiro and Wilk, 1965) and Levene’s test (Levene, 1960), respectively. The effect of otolith weight on each element was investigated by performing a non-parametric analysis of covariance (ANCOVA) (Moreira et al., 2018; Ferreira et al., 2019). To investigate any effect of otolith weight on trace element concentrations, we calculated Spearman correlations for each element (Geffen et al., 2003; Longmore et al., 2010), using the rcorr function in the Hmisc R package (Harrell, 2019). We performed a non-parametric Kruskal-Wallis statistics to investigate differences in elemental concentrations between geographical sites. We inspected pairwise comparisons between sites for each element through Dunn’s test using dunn.test function with Bonferroni correction for multiple testing. A Discriminant Analysis of Principal Components (DAPC) was performed with adegenet R package (Jombart, 2008) to predict the sampling area membership of each individual using the multivariate element composition assessing the quality of discrimination. As the optimal number of clusters identified using k-means algorithm did not provide evidence of correspondence with individual groups, we used the sampling sites as a prior. We performed a DAPC, following the procedure described in the Supplementary Materials 1.1.5 , on a combined dataset created merging genetic and otolith data when available for the same individual, resulting in 82 individuals. This combined dataset was employed to predict the sampling area membership of each individual. We ran two versions of the method, first including as input the optimal number of clusters determined based on the Bayesian Information Criterion (BIC), followed by providing the sampling sites as prior for clusters.

A decrease in the polymorphic rate of the population samples was found along the west-to-east geographical axis, ranging from 98.68% to 72.11%. Consistently, there was a significant decrease in average observed (Ho) and expected (He) heterozygosity, while no significant differences were found between Ho and He (α = 0.05). After FDR correction only a few SNPs per population sample were found significantly out of HW equilibrium. At the population samples level, no significant deviation from the HW equilibrium was found. Accordingly, the estimates of F_IS were low and consistent with levels of heterozygosity ( Table 2 ).

Pairwise F_ST values ranged from -0.003 to 0.172, with an overall significant F_ST of 0.069 (p = 0.001). After Benjamini–Yekutieli correction, all pairwise F_ST remained significant (α = 0.05) except for the comparisons involving samples from GSA6, GSA7, GSA10, GSA17 and GSA18 ( Supplementary Figure 1.1 ). A heatmap of pairwise F_ST pointed to three main areas with significant levels (α = 0.05) of genetic divergence corresponding to the Western Mediterranean, Adriatic Sea, and Levantine Sea ( Figure 2 ).

The clustering pattern of population samples resulting from DAPC analysis was concordant with the structure suggested by the genetic pairwise F_ST . Based on k-means algorithm and the lowest BIC score ( Supplementary Figure 1.2A ), three clusters were identified as corresponding to three main areas: Western Mediterranean, Adriatic Sea, and Eastern Mediterranean ( Figure 3 ). Western and eastern samples were separated along the first linear discriminant (LD1) axis. Along the second discriminant (LD2) axis a separation between samples from the Aegean and Levantine areas is shown when using prior information about the sampling sites. Moreover, the GSA18e_2000 sample separated from the rest of Adriatic samples and resulted in an intermediate position between samples from other Adriatic sites and Eastern Mediterranean area. In the DAPC analysis performed with prior knowledge of subgroups we retained 125 PCs and five discriminant functions. The proportion of conserved variance for the first two discriminant functions of whole dataset was 83.52% and the proportion of overall correct assignment was 84.38% ( Figure 3A ). Samples GSA10_2000, GSA10_2003, GSA18w_2000, GSA18e_2000 and GSA22_2009 showed a higher correct assignment rate with respect to the other population samples ( Supplementary Figure 1.3A ).

Methods for outlier detection identified different numbers of SNPs as outliers. The number of SNPs identified by pcadpat (eight SNPs) was not sufficient to reveal the population structure, thus its results were removed from further analysis. However, six of the eight loci resulting from pcadpat were also found by other methods and retained. As outlier, fsthet identified 17 loci, Bayescan identified 34 loci and Arlequin v.3.5.2.2 identified 20 loci ( Supplementary Figure 1.4 ). The sum of markers identified by all methods resulted in an outlier dataset of 43 SNP loci. Accordingly, the neutral dataset included the 337 markers which were not detected as outlier by any of the methods ( Supplementary Table 1.3 , Supplementary Figure 3.4 ). DAPC analysis performed on the neutral dataset corroborated the result obtained with the whole SNPs dataset, indicating a broad structure characterized by a separation of samples along the longitudinal axes of the Mediterranean. Two clusters were identified by the k-means algorithm and the lowest BIC score ( Supplementary Figure 1.2B ). In the DAPC analysis performed with prior knowledge of subgroups we retained 100 PCs and three discriminant functions. The proportion of conserved variance for the first two discriminant functions of neutral dataset was 77.51% and the proportion of overall correct assignment was 59.70% ( Figure 3B , Supplementary Figure 1.3B ). DAPC analysis performed on the outlier dataset revealed higher levels of structure detecting potential distinct subgroups within Western Mediterranean area in respect to whole and neutral datasets. Along the second discriminant (LD2) axis, a deeper separation between samples from Tyrrhenian and western basin areas is shown. Four clusters were identified by the k-means algorithm and the lowest BIC score ( Supplementary Figure 1.2C ). In the DAPC analysis performed with prior knowledge of subgroups we retained 25 PCs and three discriminant functions. The proportion of conserved variance for the first two discriminant functions of outlier dataset was 87.37% and the proportion of overall correct assignment was 47.44% ( Figure 3C , Supplementary Figure 1.3C ). In addition, we obtained concordant results from STRUCTURE, provided in Supplementary Figure 1.5, 1.6 , resulting barplots were highly consistent with the pattern revealed by DAPC analysis, supporting evidence of differentiation between GSA9 and GSA10 and other Western Mediterranean samples collected from the French and Spanish coasts (GSA6 and GSA7).

Eight individuals were detected as outlier data points and removed from the shape index dataset. Two of the five shape indices (circularity and form factor) were not normally distributed and removed from further analysis. Roundness, ellipticity and rectangularity were normally distributed (p < 0.01), thus retained for further analysis. Levene’s test for homogeneity of variance was significant (p < 2.2e-16}, α = 0.01). We found a significant negative correlation between roundness and ellipticity (r = − 0.802, p < 0.001), but none of the indices were correlated with more than one other index (p < 0.001). None of the shape indices was correlated with total length of fish, and only ellipticity showed a significant association (p = 0.008). The size effect on ellipticity was corrected by dividing values by the common ellipticity/fish length relationship since there was no significant interaction between total fish length and sampling sites (α = 0.01). Site was a significant explanatory variable at α = 0.01 (adonis test in the PERMANOVA analysis), with all pairwise site comparisons significantly different (p < 0.05), as reported in Supplementary Table 2.1 , except for GSA9_2009-vs-GSA17_2009, GSA17_2009-vs-GSA22_2009 and GSA17_2009-vs-GSA24_2009. The result of dispersion test was not significant (p = 0.159, α = 0.01) validating the adonis test results. CAP analysis highlighted a separation of samples from the Gulf of Lion with 75% correct reclassification compared to lower values for the other locations ( Figure 4A and Supplementary Table 2.2 ).

The otolith shape contour was summarized by 13 harmonics (thus 52 coefficients) as they reached 99% of the cumulative Fourier power ( Supplementary Figure 2.1 ). All pairwise comparisons among sampling sites resulted in significant differences in the PERMANOVA analysis (p < 0.05) ( Supplementary Table 2.1 ). In contrast, the dispersion test resulted significant with a p-value of 0.001 (α = 0.01) making the adonis test results less reliable. The CAP on Fourier descriptors using 13 harmonics identified Aegean and Turkish samples as separate groups ( Figure 4B ). In both cases Aegean samples (95.74%) showed a higher reclassification success compared to Turkish samples (68.97%) ( Supplementary Table 2.2 ). The CAP on integrated morphometric data identified the samples from the Aegean Sea and Gulf of Lion as separate groups ( Figure 4C ). Integrating the two data sets provided higher overall reclassification success, with the Aegean samples showing the highest value of correct reclassification ( Supplementary Table 2.2 ).

From the analysis of elemental data, only Zn (p = 0.02) fulfilled the assumption of normality and homogeneity (p < 0.01) for both core and edge data. ANCOVA analysis suggested the influence of otolith weight only on elemental concentration of Rb for α = 0.01 for core data, while for edge data the influence of otolith weight on elemental concentrations of Ba, Rb and Mg was significant (α = 0.01). None of the elements at either otolith core or edge were correlated with otolith weight (α = 0.01), based on Spearman’s correlation. Significant differences (α = 0.05) in the elemental concentration were found between sites for all the elements of both core and edge data ( Figure 5 ). The sampling sites were weakly identified by DAPC methods for core data, except for GSA22_2009, while the correct assignment of individual in sampling sites was higher for edge data ( Supplementary Figure 2.2 ). The elements that most contributed to the discrimination of geographical sites were Cu for otolith edge data and Ba and Zn for otolith core data ( Supplementary Figure 2.3 ).

DAPC analysis performed on the combined dataset confirmed the result obtained with the different dataset, indicating a broad structure characterized by a separation of samples along the longitudinal axes of the Mediterranean. Three clusters were identified by the k-means algorithm and the lowest BIC score ( Supplementary Figure 2.4 ). In the DAPC analysis performed with prior knowledge of subgroups we retained 33 PCs and four discriminant functions. The proportion of conserved variance for the first two discriminant functions of combined dataset was 70.02% and the proportion of overall correct assignment was higher (91.77%) compared to the different dataset. Samples of combined dataset showed a higher correct assignment rate ( Supplementary Figure 2.6 ) with respect to the other population samples considering only 380 SNPs ( Supplementary Figure 1.5A ). DAPC analysis performed on the combined dataset confirmed higher levels of structure detecting potential distinct subgroups within Western Mediterranean area with a separation of samples between GSA7_2009 and GSA9_2009. The analysis did not show a clear separation between samples from GSA22_2009 and GSA24_2009, probably due to the small number of individuals of GSA24_2009 for which we could retrieve both genetic and otolith data. However, the clustering obtained using this combined dataset groups several samples of Eastern Mediterranean origin with most of the Adriatic samples, as reflected by the third cluster displayed in Supplementary Figure 2.5 , probably due to the unbalance number of individuals.

We identified 13, 55 and 106 SNPs from Bayenv2, LFMM and Samβada respectively, which showed at least one significant association with an environmental variable. The intersection of the SNPs identified among the three GEA analysis (Bayenv2, LFMM and Samβada) resulted in six loci ( Supplementary Figure 3.2 , Supplementary Table 1.3 ), the intersection among outlier loci and GEA methods resulted in 22 loci ( Supplementary Figure 3.4 ) and the union of all SNPs identified by GEA methods resulted in 148 SNPs ( Supplementary Table 1.3 , Supplementary Figure 3.4 ). The regression analysis of locality scores from the sPCA performed on the GEA union dataset of 148 SNPs returned only dbMEM1, dbMEM3 and dbMEM5 as significant spatial variables, while dbMEM2 and dbMEM4 were not (p > 0.05) ( Table 3 ). All the environmental variables were significantly associated with the SNPs, sal_bottom_mean was the best predictor among environmental variables, followed by three weaker but still significant predictor variables (curM_Wclm_mean, temp_bottom_range and pco_Wclm_mean). Longitude was the best predictor of locality scores, followed by Latitude ( Table 3 ).

The global RDA model performed on both outlier and neutral SNPs identified significant patterns of variation (p < 0.05). The percentage of genetic variation explained was 38.13% for the outlier dataset of loci and 9.97% for the neutral ones, with the first two axes accounting respectively for 35.61% and 7.82% of the explained variation. The forward and backward procedure using ordistep function selected different spatial variables for the two datasets of SNPs, but temp_bottom_range, sal_Wclm_mean and pco_Wclm_mean were identified for both RDA on neutral and outlier SNPs ( Figure 6 , Supplementary Table 3.1 ). From marginal ANOVAs, spatial and environmental variables selected were all significant predictors of neutral and outlier genetic variation (p < 0.05). The results of variation partitioning based on neutral SNPs indicated that spatial variables (1.9%) explained approximately the same individual fraction of variance as the environmental variables (1.5%). Shared variance (6.3%) resulted in the highest amount of explained variance ( Supplementary Figure 3.5 ). In contrast, the spatial variables explained more of the individual fraction of variance in the outlier SNPs (7.0%), whereas environmental variables explained 5.5%. As for neutral SNPs, the highest amount of explained variance was represented by shared variance (25.6%) ( Supplementary Figure 3.6 ). The global RDA model performed on both western and eastern Mediterranean dataset separately identified significant patterns of variation (p < 0.05). The percentage of genetic variation explained was 1.56% for the western cluster of loci and 9.14% for the eastern ones, with the first two axes accounting respectively for 1.19% and 7.93% of the explained variation. The forward and backward procedure using ordistep function selected different spatial variables for the two sub-basin datasets, with dbMEM1 and dbMEM2 associated to the western dataset, while dbMEM3 to the eastern one. Concerning environmental variables, sal_bottom_mean was identified as significant for RDA on the western SNPs cluster, while sal_Wclm_mean, curM_Wclm_mean and pco_Wclm_mean were identified for RDA on the eastern dataset ( Supplementary Figure 3.9 ). From marginal ANOVAs, the spatial and environmental variables selected were all significant predictors of genetic variation for western and eastern clusters (p < 0.05).

The global RDA model performed on both otolith shape and elemental composition identified significant patterns of variation (p < 0.05). The percentage of variation explained was 19.66%, with the first axis accounting for all the explained variation. The forward and backward procedure using the ordistep function selected two significant spatial variables for the otolith datasets, and one significant environmental variable (sal_Wclm_mean). From marginal ANOVAs, spatial and environmental variables selected were all significant predictors of otolith variation (p < 0.05).

Based on the outcomes of transcript annotation and BLAST search, association with known genes was retrieved for 360 out of 380 loci. In details, 36 out of 43 outlier loci and 132 out of 148 loci associated with environmental variables were related to functional annotation. A total of 20 annotated markers were in common between the two datasets. Enrichment analysis revealed that the vast majority of significant results was due to GO terms connected to the outlier dataset, while only a negligible number of significantly enriched terms was ascribed to the environmental dataset. GO terms were significantly over-represented in all three aspects of gene ontology ( Supplementary Figure 4.1, 4.2 and 4.3 ). Enrichment analysis highlighted 93 biological processes, and the top three predominantly enriched GO terms were regulation of nucleobase-containing compound metabolic process, cell migration involved in gastrulation and embryo development. For the 19 molecular functions, the top three predominantly enriched GO terms were nucleic acid binding, RNA and DNA binding. Finally, 12 enriched GO terms were found for the cellular components. When retrieving SNPs associated to the enriched GO terms, a total of 22 loci were associated with GO terms significantly over-represented and belong to 20 successfully annotated genes (listed in Supplementary Table 4.1 ).

We used a large, wide-range dataset to examine the population structure of common sole in the Mediterranean Sea. We used different tracers to provide complementary information and identified factors contributing to genetic differentiation. Our results confirm that a multidisciplinary approach provides the opportunity to reveal fine population structure that is often missed with single markers. Molecular approaches were used to identify structure at the microevolutionary time scale, while otolith shape and elemental composition analyses for a subset of sampling sites described the structure at a life span scale (Tanner et al., 2016). The current work confirmed the presence of a high level of differentiation in common sole populations across the Mediterranean area. Understanding these patterns of population structure may contribute to marine fisheries conservation by identifying overlap or misalignment of current management units with biological population units (Kerr et al., 2017).

We acknowledge that the temporal and spatial scales of our sampling design could be a limitation for the identification of proper stock units for management purposes. However, common sole is a long living species (13+ years as from Carbonara et al., 2023), but reaches reproductive maturity relatively early. In fact, according to Scarcella et al., 2014, individuals start to mature after the second year of life, implying that only a few generations overlapped under the sampling scheme adopted by the present work. We observed a deep differentiation pattern between Mediterranean samples that was confirmed through the application and the agreement of different analysis methods. We found a strong signal of a well-defined spatial structure inside the Mediterranean basin, detected at a fine spatial scale for the first time. Genetic analysis revealed substantial differences and clustered Northern Mediterranean population samples of Solea solea into three main groups occupying three main areas (Western Mediterranean, Adriatic, Eastern Mediterranean), with higher levels of heterogeneity in the west compared to the east. We are aware that our study relied only on samples from the European areas, that are not fully representative of the species’ distribution in the Basin. In fact, Solea solea is the only Solea species widely distributed in the whole Mediterranean Sea, while the other two congeneric species currently occurring in the Basin have a different distribution following a longitudinal gradient, with S. senegalensis inhabiting the western Mediterranean Sea and S. aegyptiaca dwelling in the eastern Mediterranean Sea, even if the latter could be present in a wider area (Ofelio et al., 2020). The presence of three different Solea species in the Mediterranean Sea has been investigated, revealing evidence of hybridization among S. senegalensis and S. aegyptiaca (Souissi et al., 2017, 2018), but not with S. solea. In contrast, the latter share conserved external morphology with S. aegyptiaca, being currently considered as cryptic species, but reproductively independent (Sabatini et al., 2018). To overcome the limitations of our study and fully resolve the genetic structuring of the species, the collection of samples from geographical sites in the southern part of the Basin (i.e., the African coast and Sicily channel) would be beneficial to enhance the usefulness of the findings and the power to discriminate substructure at a finer geographical scale within the three main areas, even if sampling could be hampered by the aforementioned pitfall in species discrimination. In fact, attempts to collect S. solea specimens were carried out along the coast of Egypt, but only S. aegyptiaca could be retrieved, preventing further efforts at the time of the sampling campaign. A further differentiation pattern within the Adriatic Sea was confirmed, with Northern Adriatic and South-Western Adriatic population samples (GSA17_2009 and GSA18w_2000) separated from South-Eastern Adriatic population samples (GSA18e_2000). This pattern has already been suggested in previous studies based only on mitochondrial markers (Guarnieo et al., 2002; Sabatini et al., 2018). However, stock assessment has historically been carried out using only the GSA17 (Northern Adriatic Sea) as the management unit since the landings of common sole in the western part of the GSA18 are negligible (Sabatini et al., 2018). From the analysis on genetic diversity, common sole samples from the Eastern Mediterranean showed the lowest values of observed heterozygosity compared to those from other locations. A decreasing polymorphic rate could be indicative of demographic processes as bottleneck events, potentially linked to high mortality, although the causes have yet to be determined (Pinsky and Palumbi, 2014). As suggested with other species such as European hake (Merluccius merluccius) (Morales-Nin et al., 2014), the genetic differences observed between population samples of Solea solea into the three main areas in the Mediterranean basin are probably led by neutral drift. Limited connectivity across the Mediterranean is likely to be due to surface water-masses circulation with several cyclonic gyres as well as to biogeographic barriers (El-Geziry and Bryden, 2010). The life history traits of common sole also contribute to limiting its dispersal. Sole larvae start metamorphosis around 30 days after hatching and benthic settlement is completed within 2 months from hatching (Horwood, 1993). During the pelagic larval stage, larvae can travel considerable distances (Lacroix et al., 2013) but the retention in proximity of suitable nursery habitat, such as estuaries and lagoons (Primo et al., 2013; Morat et al., 2014), is a crucial factor for survival. Sole larvae have some ability for active swimming and use tidal transport to follow freshwater plumes from coastal rivers to reach their coastal nursery grounds (Lagardere et al., 1999). The connectivity pattern of sole is likely influenced by the coastline configuration and the availability of suitable environment (Tanner et al., 2017). For instance, the Strait of Sicily has been considered a geographical boundary separating the Eastern and Western Mediterranean at phylogenetic and population levels (Patarnello et al., 2007). In addition, the present-day genetic and demographic features of Mediterranean common sole could be also the signature of past processes of isolation into Pleistocene glacial refugia. These likely occurred in the Mediterranean linked to cycles of glacial and interglacial periods, with remarkable climatic fluctuations and severe sea-level changes (Sarnthein et al., 2003) that affected the life history of several marine species. Ecological niche models based on palaeoclimate reconstructions of sea surface temperature and bathymetry applied to zooarchaeological finds dating to the Last Glacial Maximum (ca.21ka ago) show strong evidence of important fish invasion and displacement through the Straits of Gibraltar in glacial periods - where they were exploited by Palaeolithic human populations around the western Mediterranean Sea - challenging an existing paradigm of marine glacial refugia (Kettle et al., 2011). The pattern of population structure was confirmed by all the genomic data analysis methods. DAPC analysis confirmed the separation of western from eastern population samples and, in the eastern population samples a separation from Aegean and Levantine areas. Furthermore, the complete overlap in the results of the temporal replicates available for almost all GSAs reinforces our findings. The use of different population differentiation methods helps to address both demographic and adaptive processes, indicated by variation at the neutral and outlier SNP loci, respectively. The presence of three main clusters was in fact revealed by using a complete genetic dataset as well as a selected dataset of only neutral markers. Genetic differentiation at a finer spatial scale was highlighted by the analysis conducted on a separate dataset of outlier SNPs. Consistency between neutral and outlier SNP loci signal was well supported, but the outlier SNP loci had greater resolution and revealed more fine scale segregation patterns compared to neutral markers. Outlier loci were able to distinguish the Tyrrhenian Sea population samples (GSA9 and 10) from the other Western Mediterranean samples. The subtle genetic clustering in the Tyrrhenian area may suggest that local selection could play a detectable role in shaping the genetic structure of common sole, since the analysis was based on outlier markers, putatively under selection. Diopere et al. (2018) used the same panel of genomic markers to study common sole populations in the North-Eastern Atlantic and obtained comparable results, indicating that the genetic structure of common sole might be shaped by a combination of neutral and adaptive processes. It appears that neutral genetic structure is due to dispersal limitation of this species throughout its distribution, while adaptive processes such as environmental and ecological features might play as key factors in local adaptation. In any case, we have demonstrated that outlier detection can reveal fine-scale population structure in areas in which connectivity and gene flow are limited as in our study. Notably, when applied in Atlantic area where gene flow due to geographical boundaries leads to stronger connectivity links, the same panel of loci supports the evidence for local adaptation (Diopere et al., 2018). Sixteen markers were identified as outliers in common in both Atlantic and Mediterranean studies, providing further insights of their involvement in adaptive processes (Diopere et al. (2018); Supplementary Table 3.1 ).

Another goal of the study was to investigate the influence of environmental variables on the population structure within a seascape analysis framework. Studies have shown that the use of constrained ordination methods such as RDA are effective tools in detecting the genetic basis of local adaptation (Forester et al., 2018). Seascape analysis on both neutral and outlier SNPs datasets identified clusters confirming the distinction between common sole population samples from the Western and Eastern Mediterranean areas, with the longitude as the best predictor. Furthermore, RDA models performed considering western and eastern Mediterranean dataset separately, confirmed the pattern of variation selecting as significant predictors both spatial and environmental variables. However, a larger proportion of the genetic variation was explained by clusters in the outlier loci dataset, confirming the higher power of resolution of outlier markers than neutral ones. In this respect, RDA models, which included spatial distribution and environmental factors, explained 38.13% of the outlier genetic structure while only 9.97% of the neutral one. Most of the variation remains unexplained in RDA models probably because it can be explained by additional factors. In fact, neutral genetic variation could be due to random processes such as genetic drift or mutation, but also other factors can contribute to demographic isolation (Gagnaire et al., 2015).

The analysis of environmental data and their association with genetic variation suggested that the population structure of Mediterranean common sole is likely to be partially due to environmental factors such as temperature or salinity and to spatial distribution. Salinity and temperature are factors shaping the survival and reproduction of common sole. Salinity is a key parameter for early stages which aggregates in estuarine areas and is a proxy for river input, which amount, and timing can have an influence on juvenile growth and survival driving therefore recruitment strength (Salen-Picard et al., 2002). Temperature is known for influencing the spawning seasonality (Lacroix et al., 2013). Studies shows that reproductive biology can differ within the basin, with a smaller length at first maturity observed in the eastern Mediterranean (19.6 cm; El-Aiatt et al., 2019) compared to the Adriatic Sea (25.9 cm; Carbonara et al., 2023). Suitable habitat for common sole is fragmented in the Mediterranean Sea and likely displays a longitudinal gradient in temperature and salinity conditions, with the eastern basin saltier and warmer. Statistical analysis suggests that the environmental variables indeed exhibit a spatial structure. Consequently, when arguing that spatial and environmental factors might significantly impact genetic structure, it is not possible to exclude that link between genetic differentiation and temperature or salinity is confounded with a broad geographical pattern. Especially for the outlier dataset, the spatial variables explained 7% of variance meaning that spatial structure related to other environmental, geomorphological, anthropic or connectivity features may act in shaping genetic differentiation among population samples at multiple scales. Among the anthropic pressures that were not accounted within this study, fishing exploitation may eventually play a role in shaping the observed variation. The common sole is a target species for commercial fisheries only in a few of the study areas (mainly the northern-central Adriatic Sea) and therefore the fishing pressure acting on the different sub-population identified is likely to be heterogeneous and spatially structured. High quality fishery-dependent data regarding catches and spatial distribution of trawlers are now common (e.g. DCF for catches, Global Fishing Watch for effort), but this was not the case during the time-period covered by our samples. Information for the beginning of 2000’s is eventually scarce and not coherent across the study area, therefore we decided to not include fishery data in the drivers to avoid biased sources. Our results show a combination of spatial and environmental influence on genetic structure, consistent with the life-cycle traits of this demersal flatfish. Environmental variables also explained 5.5% of the variance in the outlier dataset that was not explained by any spatial structure, even though there are significant salinity and temperature gradients across the Mediterranean Sea (Skliris, 2014; Spedicato et al., 2022). In fact, RDA models performed on otolith dataset underlined salinity as a significant predictor of the variation encountered between the populations of Solea solea in the Mediterranean Sea. Particularly, several outlier SNPs were also associated with environmental variables in the GEA analysis, suggesting the potential involvement of these loci in local adaptation.

Shape analysis of the otolith outlines showed statistically significant differences among the five population samples, hence otolith shape analysis is confirmed a valuable discriminating method in marine fishes. The analysis performed on shape indices showed significantly different dimensions for GSA7_2009 compared to the other samples, with a higher correct assignment to the original locations. This separation could suggest low levels of variability within samples from Gulf of Lion (GSA7_2009) and higher differentiation between GSA7_2009 samples and samples from the other locations. Individual variability might reduce the power of the classification method and may reflect local environmental conditions (Morales-Nin et al., 2022), although this is likely to vary between species. Even with correction for fish length effects, some otolith dimensions vary with age and growth conditions especially where there is a curvature in the otolith growth axis. This may contribute to variability within areas (Burke et al., 2008), but is not likely to be a major factor in the otolith shape in common sole. The elliptical Fourier descriptors (EFDs) were more successful than the shape indices in discriminating among individuals from different Mediterranean sampling locations. The results obtained from analysis of EFDs clearly separated Eastern Mediterranean population samples from the other sampling locations with a high percentage of accurate assignment to their original locations. Finally, overall reclassification success was higher when morphometric data were integrated, meaning that the combination of shape indices and EFDs improved the power of discrimination between different groups.

Although at a lower resolution because of the reduced sampling extent compared to genetic data, integrated morphometric data were effective in confirming significant divergence between Western and Eastern Mediterranean population samples separated from the rest of the common sole population samples, a pattern already suggested from genetics analysis. Furthermore, the analysis performed combining the genetic and the otolith dataset confirmed the agreement between the results of the two methods and underlined how the potential power of resolution could be increased by the combination of different data sources. Many factors can influence the trace element concentrations in otoliths, including environmental factors such as temperature, salinity and chemical composition of water masses, fish physiology and ontogenetic change in habitat occupation (Chang and Geffen, 2013; Hüssy et al., 2021). In relation to the otolith trace element composition, the correct population assignments depend on the otolith material deposited in a localized geographical area. The otolith core forms during embryo development and represents the early life history stage, including the spawning area, the pelagic or dispersal phase and possibly movement to nursery areas. The otolith edge represents the most recent resident area of that individual. In our case the discrimination methods detected stronger signals in the core trace element data, with population samples from the Aegean (GSA22_2009) more clearly differentiated. Whether this is limited by the number of samples in each population, or a reflection of the chemistry of the water masses, is difficult to say (Vasconcelos et al., 2007). In fact, there is an improvement in the discrimination power of core data compared to edge data, but not with perfect correspondence because of fish movement after their early life stages. Multiple factors might be limiting the power of discrimination methods, such as the small number of elements used in the analysis, low residence time in the location due to habitat shift in such a complex life cycle as the one of common sole. The core and edge mean concentrations measured in the single-element analysis ( Figure 5 ) were in line with the values assessed by other studies (Cuveliers et al., 2010; Chang and Geffen, 2013). The presence of barium (Ba) as one of the microelements that most contributes to the discrimination between groups in the core data is consistent with the fact that Ba is known as a proxy of coastal or freshwater influence occurring in nursery area ( Supplementary Figure 2.3 ) (Morat et al., 2014), as well as areas of high primary productivity (Chang and Geffen, 2013). Among other elements, copper (Cu), zinc (Zn) and lead (Pb) seem to be the trace elements that most contribute to the discrimination between groups both in core and edge data ( Supplementary Figure 2.3 ). This is consistent with the higher concentrations detected in fish from GSA7 and GSA17 which receive high levels of river input from industrial areas in the south of France and industry and agricultural areas in the northern Adriatic (Valbonesi et al., 2003; Zorita et al., 2008). Although coastal and shelf waters in many areas are naturally enriched with heavy metals, high concentrations of Cu, Zn and Pb can also indicate industrial, urban, or agricultural pollution detectable in otolith trace elements (Papadopoulou et al., 1978; Geffen et al., 2003; Chang and Geffen, 2013; Constenla et al., 2021). As reported in Constenla et al. (2021), common sole could be used also as a sentinel species for industrial pollution in coastal waters. The strontium (Sr) concentration in otoliths is sensitive to both salinity and water mass composition. This may explain why GSA22_2009 was so clearly differentiated from the other sampling sites based on the otolith core data. The salinity regime of the Eastern Mediterranean basin is substantially higher compared to the other parts of the Mediterranean ( Supplementary Figure 3.1 ) (Gačić et al., 2011). Moreover, this pattern is confirmed by RDA analysis in which, among environmental variables, salinity significantly contributed to the genetic differentiation of Eastern Mediterranean common sole population samples, although at a lower degree with respect to spatial variables. Also, among SNPs loci associated with environmental factors, salinity emerged as the best predictor among environmental variables with respect to the regression analysis on the loading score in sPCA analysis.

To validate putative SNPs under selection, research studies aim to link outlier markers to functional genes that can be involved in adaptive processes, by associating the function of known genes near to the identified locus to the variation of key environmental parameters. Several studies have successfully identified as outlier SNPs and SNPs associated with environmental variables found in genes involved in adaptation processes, using various methods of genome scan analysis (Benestan et al., 2016; S. Bernatchez et al., 2019; Selmoni et al., 2020). Different methods to detect outlier markers and loci associated to environmental factors, such as outlier detection methods and genetic-environment association (GEA), could be combined to obtain more information and enhance the possibility to identify loci potentially under selection (de Villemereuil et al., 2014; Rellstab et al., 2015). The use of additional information such as environmental variables improves the power of GEA analysis in identifying loci, compared to outlier detection analysis, and reveals signals of selection. As suggested by Forester et al. (2018), the overlapping approach which combines results from different methods can be effective in detecting adaptation while maintaining a low false-positive rate. As a result, both lists of putative loci under selection (43 SNPs) and/or influenced by environmental factors (148 SNPs) were used to perform Gene Ontology (GO) enrichment analysis. Our aim in performing the GO enrichment tests was to preliminarily explore whether the functions of known genes were under or over-represented in sets of genes associated with environmental variables or outlier SNPs. Most of the loci highlighted in the enrichment analysis, which are mainly the loci also included in the environmental dataset, were associated with different gene product functions and all of them were significantly over-represented in several aspects of gene ontology. However, the successfully annotated genes connected with these loci in most cases did not directly correspond to gene expression regulation processes generally connected with our environmental variables. Nevertheless, these loci could have an impact on overall gene expression because they relate to genes involved in RNA and DNA binding, molecular interaction, and metabolic processes, among others. The exception were a few genes that are involved in muscle contraction and encoded proteins that support mitochondrial metabolic demand that were expressed in response to high pCO2 (Cline et al., 2020). Outlier SNPs 1023 and 464 that are located into these annotated genes are associated with pCO ₂ environmental variables (pco_Wclm_mean) in at least one GEA method and they are also included in the loci that most contribute to the pattern identified by sPCA. Overall, based on our results, further work including more complete genomic information is necessary to expose the influence of environmental features on the microevolutionary process shaping the population structure of this species, with special emphasis on the assessment of putative functional local adaptation processes, which, ultimately, could lead to differences in life history traits.

This study has produced a major advance in the knowledge-based management of common sole in the Mediterranean Sea, defining and assessing the population structure of this relevant demersal resource with multidisciplinary analyses. In our study, the identified west-to-east differentiation pattern observed is fully coherent with the population structure described by Rolland et al. (2007). Using whole and neutral datasets, we identified three areas (Western Mediterranean, Adriatic, Eastern Mediterranean) that largely approximate to the FAO division of the Mediterranean Sea in the Major Fishing Areas (MFAs) 37.1, 37.2 and 37.3. In addition, the outlier dataset revealed a more subtle division within the western area (37.1). A large portion of MFA 37.2 is the Ionian Sea (including the coast of Libya), and our sampling did not cover the whole Mediterranean area. Even if more samples from the Southern and insular Mediterranean Sea are needed to achieve a complete pattern of common sole population units in the Basin, our findings could provide a new perspective for management decisions by reiterating the existence of a mismatch between management and biological units. In the case of the southern Adriatic Sea, the south-east population may be closer to the Ionian samples than to the main Adriatic Basin. The genetic spatial configuration does not match the statistical unit boundaries, which aggregates the shores of the southern Adriatic and divides the western coast. Following a multidisciplinary approach, we also tried to elucidate the selection forces that influence the observed pattern of population structure. The identified genetic clusters agreed with the subdivision in MFA areas, with population samples reassigned with high success rate by DAPC analysis on genomics data as well as CAP analysis on otolith data. Moreover, the power of outlier loci in detecting finer population structure helps to reveal metapopulations, and it should be considered in the assessment and management of fisheries plans. Both genomics and otolith data provided adequate discrimination power to delineate population structure in the present study, but the temporal and spatial scales of our sampling design could be a limitation for the identification of proper stock units for management purposes. Nevertheless, no major external forcing modifications were expected to occur during the sampling period, considering that the reduction of fishery exploitation started after the conclusion of the sampling. The multidisciplinary approach in our work is in line with the GFCM 2030 Strategy to focus on subregional level processes (FAO, 2021). In addition, comparison with recent samples would allow an evaluation of the genetic diversity of the species following two decades of management measures that reduced the fishing pressure on this stock (FAO-GFCM, 2021). In conclusion, our multidisciplinary approach was reliable in revealing the population structure of common sole, and providing new insights into the importance of local adaptations that produce finer scale structure in this highly valuable marine resource. These findings represent an asset to implement an effective integrated approach combining comprehensive quantitative results to delineate better aligned stock management units for this species.