Almost one thousand wild sea breams consistently genotyped at 1159 high quality SNP

A total of 767.1 M read pairs were obtained, with an average of 804,463 ± 505,000 reads per individual (range 97–3,104,802 reads). After filtering, an average of 72.2% ± 6% of reads was retained. The initial number of called SNPs was 11,662, included in a total of 216,713 tags. After filtering out low-quality markers, 1165 SNPs (10.6%) were retained. After filtering, the level of data missing per sample ranged from 0% to 19.9%, with an average of 5.4%. Of 30,290 tests for departure from H-W equilibrium carried out, after sequential Bonferroni correction, two loci showed significant deviation from equilibrium (both for an excess of heterozygosity) in more than half of the natural populations and were excluded. A total of 675,703 tests for LD were carried out and four loci pairs showed r² values higher than 0.7 and for each pair, the locus with lower missing data was retained. The remaining 1159 SNPs were used for subsequent analysis (S2 Table). The mismatch rate between replicated samples at the 1159 filtered loci ranged from 3.4% to 5.8%, with an average of 4.0%. Mean number of alleles per locus (Na) was 1.782, ranging from 1.598 (SPA-1) to 1.931 (TUN-1). Mean expected and observed heterozygosity were 0.153 (range 0.141–0.157) and 0.138 (range 0.120–0.149), respectively (Table 1). All these parameters tended to be lower in the Atlantic populations. Minor Allele Frequency (MAF) ranged from 0.112 to 0.168, with an average of 0.135. Effective Population Size (Ne) for all Atlantic samples were very high ('Infinite'), whereas Mediterranean samples showed values ranging from 106.4 (ITA-6) to 'Infinite' (ten samples). Nevertheless, ranges of confidence varied substantially especially for smaller samples (Table 1), which suggests caution when interpreting these data. The inbreeding index varied from 0.033 to 0.105 and averaged 0.074 (Table 1). No general trend was found for this parameter when comparing Atlantic and Mediterranean samples. Relatedness analysis gave clues for the presence of potential full-sib pairs within the SPA-2 sample (Wang relatedness > 0.375). As a consequence, four individuals that were part of highly related pairs were eliminated from the dataset to avoid biases in the clustering analysis. Genetic differentiation between time replicates was not significant. Heterozygosity was also similar for the time replicates, but for Basova Kavalas and Thermaikos Gulf recent samples showed significantly higher Ne (Table 1).

Four genetic clusters of sea bream populations were identified

Comprehensively, the results obtained with the full SNPs dataset (representing the combination of demographic and selective factors) suggested a subdivision of sea bream into four major genetic clusters: Atlantic (ATL), West Mediterranean (WMED), Ionian/Adriatic seas (ION) and the Aegean Sea (AEG), the strongest differentiation being between Atlantic and Mediterranean samples (Fig 2 and Table 2). The Evanno’s method identified four clusters as being most likely and all runs at this number of ancestral groupings showed the same cluster pattern. While samples showed a high level of admixture, those from different basins differed in the proportions of each component. Within-basin differences were much lower (Fig 2). Samples from the Northern and Southern Adriatic Sea (i.e. ION01/ITA7 and WMED10/ITA6, respectively) showed no genetic differentiation (F_ST = 0.0000), however under Structure’s Bayesian analysis they clustered separately: with ION and WMED populations, respectively. AMOVA analysis further confirmed this subdivision, with the proportion of genetic variability among basins being statistically significant, while the variability within groups was not (Table 2).

Fig 2

Structure clustering analysis.

Result for the most likely number of clusters (k = 4). Labels under the graph indicate the sample groups. Labels above the graph indicate the genetic clusters suggested by the analysis and supported by AMOVA.

Table 2

Results of the AMOVA analysis with sample subdivided according to the clustering analysis.

Source of Variation	Degrees of freedom	Sum of squares	Variance components	Fixation indices	Percentage of Variation	P-value
Among groups	3	433.523	0.29020	0.00363	0.36	<0.01
Among populations within groups	19	732.743	-0.66080	-0.00830	-0.83	Ns
Among individuals within populations	801	68161.210	4.82579	0.06012	6.04	<0.01
Within individuals	824	62165.500	75.44357	0.05576	94.42	<0.01
Total	1647	131492.976	79.89876		-	-

DAPC was based on 150 PCs, after cross-validation analysis. Differentiation signal was weak, with the first axis roughly matching the geographical West to the East distribution of the samples. Along the second axis of the DAPC ATL and AEG clusters showed some level of genetic similarity when compared to the two remaining clusters within Mediterranean Sea (S1 Fig).

In general, F_ST values were low (average F_ST = 0.0031) with a trend of increased values in the comparisons between Atlantic and Mediterranean (S3 Table). Within the Mediterranean, values tended to increase and be more significant in the comparisons between samples from the Western and Eastern Mediterranean (Ionian and Aegean basins). When samples were subdivided into the proposed four groups, F_ST values were low but highly significant (p<0.001) and ranging from 0.0015 in the comparison between WMED and AEG to 0.0070 in the comparison between the ION and ATL clusters. No population specific private alleles were observed. However, when analysed at group level, an allele at locus 7704_13 (S2 Table) was observed only within the WMED group, with nine heterozygous individuals being observed (allele freq = 0.024)

Dissecting genetic differentiation: Atlantic-Mediterranean differentiation at neutral loci and signs of convergent selection at outlier loci

Bayescan analysis detected 20 loci with log10(PO) > 1.5; while Arlequin detected 28 potential diverging outliers at p < 0.05. Of these, 16 loci (1.3% of the entire dataset) shared by the two methods were selected to create the ‘outlier dataset’ (OL, S2 Table). Two additional loci were included in the OL dataset after the analysis with samples divided into four clusters (see above). The OL dataset was subsequently used for analysis focused on exploring the “functional” divergence between populations. The putative neutral loci dataset consisted of 1126 remaining loci, after excluding potential outliers identified by at least one of the approaches used.

Single SNP allele-specific assays were designed and performed for the 16 identified outlier loci set (data not shown). Expected polymorphisms were validated for all but one locus (8727_39), which was removed from subsequent analyses.

BayEnv detected 24 loci, of which only six were among those identified by the joint analysis carried out with Bayescan and Arlequin (S2 Table).

OL allele frequencies did not show a constant behavior when moving from West to East populations (S2 Fig): some followed a gradual variation (e.g. 2689_65); some others showed more abrupt changes at basin boundaries (e.g. 12615_64). finally, OL identified by BayEnv showed patterns of allele frequency that overlapped with the correlated environmental variables variation (e.g. 13310_71).

The genetic structure detected from analysis of the panel of putative neutral loci alone was weaker than that when using the entire dataset, especially within the Mediterranean. Pairwise F_ST values were rarely significant, occurring mainly in comparisons between Atlantic samples and samples from Ionian and Aegean seas (S4 Table). Structure clustering analysis, using the most likely value of k = 2 (according to Evanno's method), showed a high degree of admixture, though differentiation was apparent in cluster representations between ATL, WMED and East Mediterranean (ION + AEG) samples (S3 Fig).

The presence of within basin differentiation emerged (from analyses involving the OL panel (S4 Table). The Northern Atlantic sample ATL1/FRA-1 was differentiated from more southern samples ATL2/SPA-1 and ATL3/SPA-2. In general, ATL1/FRA-1 is the most differentiated sample with pairwise F_ST values ranging up to 0.1751 (with ION3/GRE-2 sample). WMED was the only group found to be highly homogeneous at outlier loci. Astriking difference uncovered by the OL panel analyses compared to the neutral loci data set was the strong differentiation between Ionian and Aegean samples, reflected in both F_ST values and Structure results at the most probable k = 3 (S4 Fig).

Genetic differentiation at outlier loci reflected environmental differentiation between localities, as confirmed by SEASCAPE analysis, which indicated longitude and winter surface salinity (S_0_W) as significant explanatory variables (Fig 3). Climate explained 54% of the variance in the data, while geography explained 14%. The joint effect of the two sets of variables explained 32% of the total variance. The PCA plot indicated that longitude had the strongest effect on the first axis, while winter surface salinity was correlated with the second axis. Along this axis, ATL samples are roughly grouped in the lower half of the graph and closer to AEG and ION-1/ITA-7 samples, reflecting differences in salinity values. ION and ATL samples showed the strongest pairwise differences both in terms of environmental variables and also in terms of differentiation at OL loci.

Fig 3

PCA of the redundancy analysis performed using outlier SNP dataset.

Sampling groups are labeled according to the genetic clusters identified with Structure and AMOVA and plotted according to the most explanatory axes. Longitude and surface winter salinity (S_0_WIN) were the environmental variables most significantly related with genetic data.

Almost all of the outlier loci (28 of 34, including those identified by Bayenv) were detected in the genome assembly. Mapped loci were scattered among 16 of the 24 chromosomes of the species. Eighteen of them were located within genes, 11 of which were among the markers indicated by BayEnv as correlated with environmental variables. Fourteen of these tags were located in a Coding DNA Sequence (CDS) region (S2 Table). The enrichment analysis suggested a significant over-representation of genes involved in the biological processes “kinetochore assembly” and “kinetochore organization” (GO:0051382 and GO:0051383 respectively) and “endoplasmic reticulum calcium ion homeostasis” (GO:0032469). In addition, pathways related to the expression of cellular membrane components and many “transferase” molecular functions were highlighted (S5 Table).

Long term demographic factors explain the neutral genetic structure

F_ST levels, calculated from the full SNP dataset, were lower than those usually found in marine fish (average F_ST = 0.062 [52]) and agree with prior sea bream studies [18, 19]. A previous broad geographical range analysis of gilthead seabream, using allozyme and microsatellite markers (putatively neutral markers) by Alarcón et al. [18] concluded that the population structuring was not associated with known geographic or oceanographic factors. In the current study, a much larger marker panel was employed in an attempt to identify and resolve more subtle genetic structure, if present. Despite pairwise genetic differentiation levels being statistically significant in only a few comparisons, other approaches (i.e. Structure, DAPC and AMOVA) pointed to a significant subdivision of the species in three main basins: the Atlantic, the West Mediterranean and the East Mediterranean, including the Ionian and Aegean seas.

The presence of two major water mass boundaries separating these basins, namely the Almeria-Oran front and the Sicily Channel front, likely explain this pattern. In addition, the Otranto Strait front effect may also explain the differentiation between WMED10/ITA-6 and the Ionian samples. The population structure of another coastal fish species, the East Atlantic peacock wrasse, has previously been associated with these fronts [53]. The authors also found a connection between Ionian and Adriatic samples along the East coast of these basins. The same oceanographic features could explain the similarity between Ionian and North Adriatic (ION1/ITA-7) samples. These patterns could also be influenced by a stronger dispersal along the coast rather than across the opposite stretch of the sea, considering that sea bream prefers shallower water [54].

Evidence of gilthead sea bream genetic differentiation between Atlantic and Mediterranean basins has been previously reported using other typically neutral markers (i.e. microsatellites in [19]), and significant differentiation has also been found between Atlantic samples from north and south coasts of Spain [18]. In the current study, genetic structuring was reflected by loci under selection and thus may be the result of evolutionary responses to different environments (see further below). Atlantic samples collected closer to the Mediterranean (ATL3/SPA-2) were more strongly differentiated from Mediterranean samples than the more distant ATL1/FRA-1 sample. A similar result was found by Alarcón et al. [18], using allozymes and microsatellites. Genetic differentiation within the Mediterranean has been previously reported by other authors [22, 55, 56]. These studies focused on specific areas of the Mediterranean and agreed with our study in detecting the presence of a genetic differentiated cluster in the Adriatic basin [22, 56] and reduced gene flow through Strait of Sicily [55]. They also suggested the presence of differentiation at finer level between samples from West Mediterranean (e.g. between Tyrrhenian and Sardinian seas) based on F-statistics, but this result was not confirmed by their Bayesian analysis [22, 56]. Our results confirm the absence of differentiation, which is not detectable also using the more sensitive OL panel. Low levels of differentiation within basins, coupled with limited differentiation across fronts are not unexpected for this species. Indeed, given the long period of larval dispersion of seabream (>30 days) [57] and its pelagic/benthic vagile lifestyle, retention rates within basins are expected to be low, with fronts and currents having a dominant effect on shaping the neutral genetic structure. In a previous study, based on microsatellite markers [58], some degree of differentiation between samples from the north and south west Mediterranean Sea was reported. This was not apparent in the current study, with the single sample (Tunis) from the south west Mediterranean Sea area grouping with the other West Mediterranean samples.

The lack of previous genetic analyses of sea bream wild populations based on SNPs does not allow for direct comparisons of genetic diversity within the species, but the value reported here (He = 0.153) was slightly lower than that reported for other species using SNP panels (0.198–0.220 in East Atlantic Peacock wrasse [53]; 0.233–0.290 in European sea bass [59]; 0.270–0.310 in Atlantic herring [60].

Outlier allele frequency analysis: Convergent adaptation in distant sea bream populations?

Analyses based on outlier SNPs revealed a) a strong subdivision of East Mediterranean wild populations into the two “sub-basins” ION and AEG, and b) similarity between Atlantic and Aegean samples, despite the geographic separation. This is reflected by F_ST values showing the highest values between ATL and ION samples and highly significant values also between ION and AEG. Within-basin, significant differences were identified in the Atlantic between the environmentally different groups ATL1/FRA-1 and ATL3/SPA-2, which were not apparent from analysis of the entire SNP dataset. The role of environmental variables is also highlighted by seascape analysis. Separating loci influenced by varying environmental factors from those whose variation is solely due to demographic factors is challenging [61]. The patterns of some OL loci across populations showed allele frequency variations more likely related to geographical distance (e.g. 10524_58), but other OLs seem to reflect the effect of barriers between basins (e.g. 7513_19 and 12615_64). This distinction is important for identifying genes involved in adaptation; in this case, a test based on the correlation between allele frequency and environmental factors is expected to work better [49]. Bayenv analysis added 18 potential outlier loci to those previously detected with Bayescan and Arlequin, many of which showed correspondence with annotated genes and identified two biological processes likely related to adaptation to the environment. While it is acknowledged that the functional analysis conducted was limited, the results point toward a real functional role of the regions sourrounding these SNPs in the adaptation to the environment and deserve further exploration. Therefore, Bayenv, used in combination with redundancy analysis can be useful to detect hidden signature of adaptation despite the reported rate of false positives (around 20–50% [62–64]).

The results presented in this paper are a pivotal step towards the development of practical traceability and conservation tools for the gilthead sea bream. In addition, in a context of massive aquaculture production of this species, an accurate assessment of baseline wild genetic variability and structure is fundamental to the monitoring of potential aquaculture impacts going forward. Sea cage fattening is indeed a source of potential escapees of cultured individuals (often selected for productive traits and/or coming from different areas) into the wild, which could alter the genetic characteristics of wild counterparts [65]. The data from Greek temporal replicates, despite being limited in space, suggest that during the last ten to 15 years the genetic makeup of local wild populations has not been significantly impacted by aquaculture, despite many reported escape events, as also found in other studies of the same area [66].

Assessing population specific neutral and potentially adaptive genetic traits will allow a proper monitoring of the changes that escapees can induce if they successfully reproduce in the wild. The results of our study, based on samples from the entire gilthead sea bream distribution, described the presence of a genetic population structure of the species in the Atlantic and Mediterranean basins and answered biological questions that will support the management of wild stocks of this economically important species and can be used as a baseline for the assessment of the genetic impact of sea cage aquaculture.