Spatial connectivity pattern of expanding gilthead seabream populations and its interactions with aquaculture sites: a combined population genetic and physical modelling approach

In gilthead seabream the number of domesticated individuals increased annually, and escape events occur regularly in the Adriatic Sea. Still there is a lack of population genetic characteristics and evidence of the extent and geographic scale of interbreeding resulting from fish-farm escapees. We screened 1586 individuals using a panel of 21 neutral microsatellite loci in several consecutive years and here report on the medium-scale detection of hybrid and farmed seabream in the natural environment. Wild adults showed a lack of genetic structure within basin and sampling years and reduced connectivity with wild offspring collection, suggesting their temporal residency within the Adriatic. On the contrary, by linking the results of multiannual genetic analyses with the results of coupled hydrodynamic and individual based models (IBM-Ichthyop), we observed a strong connection of wild seabream associated with tuna-aquaculture sites and offspring from the nursery grounds, indicating that the surroundings of tuna sea-cage farms can function as a spawning grounds. The study results present the genetic baseline of wild and farmed strains from the eastern Adriatic Sea, as a first step toward development of a mitigation strategy for fish escapees aimed at controlling further erosion of genetic integrity.

by farmed fish inside sea-cages, i.e. escapes through spawning 4,[9][10][11][12][13] . Escaped fish have been shown to have detrimental effects on native fish stocks, due to competition for resources, spread of disease, and alteration of genetic diversity due to hybridisation 14 . Hybridization events may break down locally adapted gene complexes, reduce competitive ability and overall fitness of wild populations 15 . The impact may be more enhanced in cases when farmed individuals are of a different geographical origin, as is the case in Croatian aquaculture, and/or have experienced strong domestication 9,16 .
On the other hand, ocean warming may also be behind the increasing gilthead seabream abundance. This is currently one of the main driving forces causing changes in species abundance and distribution and, thus, in species biodiversity and functioning of marine ecosystems [17][18][19][20] . Annual trend distribution of the Mediterranean sea surface temperature (SST) indicates a warming trend throughout the basin with average values of 0.035 ± 0.007 °C per year 21 . Thus, it could be argued that wild gilthead seabream, as a subtropical Sparid, is taking advantage of the present climate change in terms of increased larval survival and subsequent recruitment success, but also in establishment at the northern limits of distribution areas (e.g., Brittany and Denmark coast) 3,22,23 .
To date, molecular markers have been successfully employed in the discrimination of fish origins and in identifying hybridization processes among wild and farmed fish 4,9,[24][25][26][27][28] . In the case of gilthead seabream, the application of molecular markers is a challenge in detecting escapees and hybrids, since (i) gilthead seabream have been domesticated relative recently, i.e. have a weak selection footprint 29 ; (ii) the high dispersal capacity of this species produces a high exchange of migrants among subpopulations, allowing a large effective subpopulation size and low structuring, and (iii) historical samples are not available for comparison analyses.
To explore the genetic population structure and footprints of selection in gilthead seabream populations of different origin (wild vs. farmed) and ontogenetic state (juveniles vs. adults) from the Adriatic Sea, we genetically assayed 1586 individuals with 21 putatively neutral microsatellites 30 . Thus, this study aimed to: (i) examine the possible changes in genetic variation occurring over the spatial and short-term temporal scale of wild seabream populations from the eastern Adriatic, with special attention on populations sampled around fish farms and their potential impact on recruitments from surrounding nursery grounds, (ii) investigate the genetic structure of farmed populations originating from different hatchery sources and (iii) quantify the presence and magnitude of hybridization between wild and escaped farmed individuals.

Results
Genetic diversity. A total of 1586 individuals of Sparus aurata were genotyped at 21 neutral microsatellite loci ( Fig. 1; Tables 1 and S1) where the proportion of missing data per locus ranged from 0 to 1.5%, with an average of 0.3%, respectively.  Table 1. The figure has been created using MATLAB 2014a (www.mathworks.com) and GIMP 2.8.16 (www.gimp.org) software.
Several populations showed significant deviation from Hardy-Weinberg equilibrium, with tendencies towards heterozygote deficiency at the loci M3, H8, G3 and L7 (Supplementary Table S2), as revealed by Fisher's exact test. MICROCHECKER detected null alleles for these loci at estimated frequencies <0.3 for M3 and H8. Those two loci were removed from further analyses. Lower null allele frequencies <8% were detected at G3 and L7 and these loci were retained, since the estimation of F ST both using and without using the ENA correction method gave comparable results; 0.022 vs 0.023 for F ST with vs without the ENA, with overlapping 95% CI. There was no consistent evidence of linkage disequilibrium between any pair of loci following strict Bonferroni correction for multiple tests.
Among the 19 loci examined, all were polymorphic, with the number of alleles per locus ranging from 7 to 14. Significant positive correlations in allele frequencies between adjacent temporal wild populations sampled in 2009, 2015, and 2016 were recorded, suggesting stability in allele frequencies (Supplementary Table S3). Genetic diversity revealed varying degrees of genetic diversity among populations, ranging from 0.74 to 0.81 in expected heterozygosity (He) ( Table 2). Farmed gilthead seabream exhibited a significantly lower effective number of alleles per locus (6.5 vs 7.3), allelic richness (7.3 vs 8.6), observed (0.74 vs 0.76) and expected heterozygosity (0.79 vs 0.81) (ANOVA test, p < 0.05) in comparison to the wild counterparts. Wild young-of-the-year and farm associated populations followed the genetic diversity pattern observed for the wild group ( Table 2). The inbreeding coefficient, F IS , ranged from −0.06 to 0.10 and was significantly higher than zero in 33% populations of the dataset ( Table 2). Analyses of multilocus genotype data indicated strong disparity of contemporaneous effective population sizes (Ne) in respect to fish origin and ontogenetic state. On average, estimates of Ne were 3-fold, 7-fold and 17-fold smaller in the farmed group (164) in comparison to the wild group (549), farm-associated (1139) and wild YOY (2799), respectively.     www.nature.com/scientificreports www.nature.com/scientificreports/ nurseries for both sampled years (−0.002 < F ST < 0.009). Only the YOY group sampled in the area under aquaculture impact (16JB) showed gene flow towards France-origin farmed fish (15FF). All farmed populations from different origins showed significant pair-wise differentiations within the group (0.009 < F ST < 0.076) and between the rest of populations.
Initial Structure analysis of the dataset suggested K = 6 after application of the Evanno et al. 31 approach. Wild populations sampled within the Adriatic Sea were assigned to five different clusters, of which three clusters were exclusively associated with wild fish origin and two were associated with farmed origins (Fig. 2). Wild historical samples from 2009 formed the first cluster (grey), while recent temporal replicates of wild adults formed the second cluster (orange). In accordance with F ST pair-wise results, no gene barrier was observed between 2015 and 2016 temporal replicates of farm-associated adults and temporal replicates of YoYs sampled from different nursery grounds, grouping all populations into the third cluster (green). The YoY population sampled in the area near the finfish aquaculture site (16JB) was assigned to the fourth cluster (red), which also comprised farmed populations of Atlantic origin. Farm-associated adults from 2017 appeared to be admixed with individuals assigned partially to the third wild cluster (orange) and also to the fifth cluster (blue) that comprise farmed populations While no further structure was detected within the adjacent temporal wild populations (orange cluster) and within the adjacent temporal farm-associated populations and YoYs (green cluster), subsequent analysis of each of the cluster of farmed origin (red and blue) resulted in further structure in both subsets, with K = 3 for the Atlantic origin subset, grouping 16JB and 15FF populations into one cluster, and K = 3 for the Mediterranean farmed origin subset with included farm-associated adults 2017, grouping 17AK and 17AB with the eastern Adriatic broodstock (Croatia) (Fig. 2). No further substructure was detected, resulting in a total of K = 10 overall, with clear conformation of farmed origin genotypes in wild populations. These relationships between sampled populations, based on pairwise Nei's genetic distance, are depicted in the principal coordinates analysis (Fig. 3). An AMOVA arrangement based on the initial Structure analysis when K = 6 revealed that 1.9% of the total genetic variation is explained by differences between major clusters, with evident gene flow among populations within groups (0.8%). Within populations variance accounted for 97.3%. DAPC similarly identified a hierarchical structure in the gilthead seabream dataset and clustered groups following the defined structure observed in Structure software. In successive K-means clustering, the elbow curve of Bayesian Information Criterion (BIC) values indicated the optimal number of clusters of K = 10 ( Supplementary Fig. S1).
In respect to the high genetic similarity of YoYs and farm-associated populations from 2015 and 2016, using the 95% critical LOD value given by CERVUS, a total of 6 parent-offspring pairs were identified, connecting YoYs sampled from coastal nursery areas (Neretva, Raša and Pantana) with the parents sampled in the close vicinity of tuna farms. Moreover, 3 YOYs sampled in 2015 from the Neretva River were assigned to the parents sampled around the Kali tuna farm on Ugljan Island (15AK), 2 from the Raša River (2016) and one from Pantana (2015) were assigned to the parents sampled around the Brač tuna farm (15AB), respectively. No parent-offspring pairs were identified among parents from wild populations and YoYs.
Loci detection power and hybrid identification. Analysis of three simulated data sets, each having a different percentage of hybrids (15%, 33% and 66%), indicated 0.85 and 0.90 as the best thresholds to be applied to assign individuals as non-admixed or admixed in an empirical dataset, balancing among efficiency and accuracy of measurements for testing overall performance of the software Structure in hybrid detection (Supplementary Table S6). The assignment efficiency for all genotype classes was highest in the set with 10% hybrids, except for backcross classes due to a higher overlap with the parental classes that decreased with increasing hybrids in the data set ( Supplementary Fig. S2). Still, backcrossed individuals decreased both the accuracy and the efficiency of hybrid detection, as they were difficult to detect, affecting the levels of recommended thresholds (Supplementary  Table S6, Hybrid all vs Hybrid F1 F2). In comparing the Structure and NewHybrids results, a relatively similar pattern in individual assignment was observed, except for the class of 66% hybrids where high misclassification of purebred individuals as hybrids decreased the overall performance of NewHybrids.
By analysing the real dataset, the Structure results with a threshold q = 0.90 indicated that a total of 75% wild individuals belonged to one of the two recognized clusters with exclusively wild fish origin. Reciprocal hybridization of individuals from wild clusters (orange and green) was noted, where 14% of wild adults and 17.4% farm-associated and YoYs populations showed to be mutually admixed. Moreover, 2.5% individuals from the green cluster were assigned as non-admixed ones to the orange cluster, and vice versa of 4.5% individuals (Fig. 2).
Among 25% of wild individuals that were not assigned to pure wild origin, 10% displayed pure aquaculture ancestry and 15% displayed hybrid-aquaculture ancestry, including hybrids of different farmed origins (Atlantic, Eastern and Western Adriatic) ( Table 3). The assignment scores were similar to those obtained by NewHybrids using a threshold q > 0.5, where all hybrids were assigned to the F1 genotype class (Table 3). Different proportions of hybrids were observed in wild populations, with frequencies ranging from 0 to 65%. In general, the smallest hybrid frequencies were recorded in YoY populations (3-13%), whereas the largest frequencies of hybrids were  Table 1. Historical wild samples (09HW) were excluded from the analysis. (2019) 9:14718 | https://doi.org/10.1038/s41598-019-51256-z www.nature.com/scientificreports www.nature.com/scientificreports/ observed in the most recently sampled wild populations in areas impacted by tuna farms (17AK, 65%; 17AB, 37%). Moreover, all juveniles from the YoY group (16JB) sampled within the concession of the Brač farm were assigned to Atlantic farm origin (86.7%) or were identified as hybrids (13.3%). Individuals of pure farm origin, and likely recent escapees, were also found in farm-associated populations (17AK and 17AB), and in wild adult populations (16WR, 16WV and 16WK), with frequencies ranging from 6 to 43% (Table 3).
connectivity of spawning and nursery areas. Spatial particle distributions on 15 May 2016 obtained in the Ichthyop -Adriatic ROMS simulation showed that there was hydrodynamic connectivity between the aquaculture area near Brač Island, as a potential spawning ground for wild gilthead sea breams, and natural nursery ground areas near Pantana and Raša River (Fig. 4a). The same model setup showed connectivity between Ugljan Island and the Raša River (Fig. 4b). The Eastern Adriatic Current (EAC), which is prominent and stronger in southern and central parts of the basin during winter, spring and autumn, transported particles from Brač and Ugljan Islands to Raša River (see Supplementary Methods for details). The Ichthyop -ASHELF2 ROMS simulation for the same date clearly showed connectivity between the aquaculture area near Brač Island and Pantana as a natural nursery ground (Fig. 4c). Averaged modelled current fields produced by the ASHELF2 ROMS simulation indicate interplay between oppositely directed currents (see Supplementary Methods for details). Due to the impact of different wind forces during the first five and half months of 2016 superimposed on thermohaline flow, current direction varied between atypical alignment with the wind and typical inflow direction 32 .

Discussion
Significant increase of wild gilthead seabream populations in coastal areas of the Mediterranean Sea and its northward expansion has been documented recently, affecting the structure and productivity of ecosystems 3,5,23 . In this study, 19 neutral microsatellites loci were used to explore geographically fine-scale population processes of gilthead seabream within a short-temporal window in coastal areas of the eastern Adriatic Sea, to gain a deeper understanding of the factors shaping genetic connectivity and structure, in comparing 1586 sampled individuals grouped by fish origin (wild vs farmed vs farmed-associated) and ontogenetic state (juveniles vs adults).
The standardized panel of two microsatellites multiplex PCRs for gilthead seabream 30 applied in this study showed robustness and effectiveness for pedigree analysis and characterisation of populations. Allelic polymorphisms and expected heterozygosities were comparable with those observed in wild and farmed populations of   33,34 and were slightly increased in comparison to previous studies in the eastern Adriatic 4,9 where only ten loci were used.
The majority of farmed populations showed departure from HWE, due likely to non-random mating and population mixing [35][36][37] . Similar reports are known 4,9 , including HWE departure of farmed seabass in the Adriatic 24 . Adding wild fish into farmed stocks and swapping breeders between hatcheries is common practice 38 , particularly in case of gilthead seabream where wild males are continuously introduced into aquaculture since most individuals turn into females after their second year of reproduction 39 . Still, lower than average number of alleles per locus, significantly lower allelic richness, heterozygosity and effective size in cultured versus wild populations was recorded within the Adriatic (Table 2). Drops in genetic diversity as a reflection of small effective breeding population size occurs in hatcheries due to high variance in family size and fewer males than females contributing to each mass spawning event 35 . Even though estimates of contemporaneous Ne were significantly smaller in farmed populations than in other wild groups, the majority of farmed populations had F IS values close to zero, suggesting that genetic breeding programmes are implemented in commercial hatcheries of Mediterranean countries in order to control the risk of inbreeding 40 . The heterozygosity deficit and significant departures from HWE observed in Adriatic wild populations at some loci but not systematically for each location can be explained due to the temporary Wahlund effect 41 caused by several breeding subunits in a sample 4,9,36 or due to selection. Namely, human-mediated selection due to fisheries may induce evolution toward slow growth, early maturation at small size and higher reproductive investment where the latter two characteristics are already observed in certain gilthead seabream populations from eastern Adriatic 7 .
Even though gilthead seabream is an important commercial species for the aquaculture and fisheries sectors, population genetic characteristics throughout the range of species distribution are still not clearly resolved. Different types of markers produce incoherent and puzzling patterns of gene flow discontinuity among regions. Some authors indicated either a weak population structure in the Atlantic and Mediterranean Sea 42-45 or a strong genetic subdivision at short distances along the basins 9,46 . On the other hand, farmed fish generally form genetically distinct groups compared to their proximal wild counterparts 36,37,42 . This pattern is confirmed in the present study where the main source of genetic differentiation was found among farmed populations of various origins (F ST = 0.041) and between farmed and wild populations, except for interactions with wild populations (17AB, 17AK, 16JB) in which a high number of escapees/hybrids were recorded. A similar pattern of inter-and www.nature.com/scientificreports www.nature.com/scientificreports/ intra-farmed strain divergence, as a result of the selection in breeding programs, the founder effect and genetic drift, was observed for other farmed species 24,47 .
Another source of genetic divergence was found among wild genotype collections (F ST = 0.012). While wild adults from natural locations showed homogenous genetic pattern along the Croatian coast, wild adults sampled around tuna farms together with wild young-of-the-years sampled from three coastal nursery grounds showed multiannual genetic connectivity and a clear break in gene flow toward the wild adults from natural locations. At a small spatial scale, such gene flow discontinuity observed among wild fish collections adds new insight into population structure dynamic, considering that previous studies using fewer loci indicated a lack of any population subdivision within the Adriatic 9,33 . Low values of F ST and high gene flow in marine species do not necessary imply an absence of population structure 48,49 and the observed differences could be due to the different resolution of molecular markers used.
The extent of movements of adult S. aurata is still little known 50 , though the long pelagic larval duration up to 50 days at 17-18 °C is likely to contribute to the exchange of migrants among populations 51 . Considering that wild adults that appear in large school formations during late autumn and winter showed genetic divergence from wild offspring in temporal replicates, it could be argued that those populations function as independent units and conduct an exclusively trophic north-to-south migration along the eastern Adriatic coastline, and spawning occurs elsewhere. This is further supported by the parental assignment analysis where no parent-offspring pairs were identified. Thus, we can characterize wild gilthead seabream populations from Croatian waters into two contingents with non-overlapping spawning grounds, (i) migratory and (ii) residential adults, i.e. ones observed around tuna farms. Due to the high dispersal capacity of the studied species, these two contingents cannot be considered as completely isolated units, since the Structure analysis identified a bidirectional pattern of hybridization at lower frequencies (<17%). In Croatian waters, migratory adults first appear during autumn in the north-eastern parts and subsequently move southwards, and therefore, we can speculate that they originate from Italian natural lagoons, i.e. immigrant donor areas that are very important nurseries and feeding grounds for seabream along the western coast 52 . However, more sampling is required along the entire Adriatic coast to obtain a representative pattern of genetic connectivity.
Estimates of effective population size recorded from residential adults suggest the existence of large breeding populations, which is in line with several recent studies that portrayed Adriatic tuna farms areas as a permanent residency of highly abundant wild seabream 6,7 . The aggregative behaviour has been primarily attributed to the permanent and increased availability of tuna baitfish losses, composed of small pelagic species rich in omega-3 and 6 fatty acids. In addition to the observed morphological differences, farm-associated seabream had a higher reproductive potential than other wild populations, due to the fish-enriched diet 7 . It seems that reliance on this trophic source has a positive influence on successful spawning and development of fish egg and larvae 53 considering the strong genetic connectivity observed between resident adults from three tuna-farm locations with offspring collections in the three coastal nursery grounds, suggesting these sites are spawning areas. This is further supported by parental test analysis where six parent-offspring pairs were identified, connecting parents from the farm surroundings of Brač and Ugljan Islands with offspring from all three nurseries (Raša, Pantana and Neretva) situated along the coastline. In addition, the coupled numerical system with ROMS and Ichthyop models successfully utilized in tuna offspring distribution studies within the Adriatic 54 , has confirmed a possible passive transport between tuna-farms areas off the islands of Brač and Ugljan and the Raša and Pantana nursery grounds in 2016.
Both the multivariate and Bayesian based clustering analysis identified similar structure architecture within the basin, supporting the populations F ST pair-wise estimates. Hierarchical Structure and DEPC analysis gave 10 distinct clusters, three of which were exclusively associated with the wild fish origin and seven with farmed ones, where each unit corresponded to a different hatchery source of fry. Namely, due to the limited number of national hatcheries, the majority of fry is still imported from France, Italy and Greece. While all historical samples were completely assigned to a single cluster, several wild populations showed to be heavily admixed with different farmed origins. This is especially true for the recently sampled populations around both tuna farms (17AK, 17AB) where a complete shift in population structure was recorded in 2017 in comparison to 2015 and 2016. Both tuna-farms are located few kilometres from seabream-cages farms, and it is possible that seabream escapees from neighbouring farms were attracted by tuna waste fish feed at the time of fish sampling. Namely, the presence of other farms has influenced the movement patterns of escaped fish, by moving repeatedly from the release cages to other adjacent farms 55,56 . The overall number of escapes detected in the wild samples was close to 10%, which is consistent with previous studies 9,24,57 ; however, more than half of the escapees in this study were captured at the aquaculture sites. Moreover, offspring 16JB collected from the natural bay near a seabream farm in middle Adriatic showed pure Atlantic-farmed origin, and they clustered together with farmed adults sampled in 2015 from the same area, suggesting that they may have originated from eggs spawned by farmed fish inside the cage. Escape through spawning has already been recognized in Greece and authors suggested that this mechanism may be responsible for recruitment peaks in sea bream populations in certain coastal areas of Greece 5 . According to the risk assessment of the environmental impact of farming 58 , the observed number of escapees (<10%) in this study presents a moderate risk for further genetic changes in wild populations. Interestingly, genetic temporal changes in population structure has already been recorded here, considering that the wild historical samples collected during an earlier stage of farming (in 2009) showed high genetic divergence toward wild contemporary samples. Long-term and cumulative introgression of farmed fish and the density of the native populations have been recognized as the main drivers of the temporal genetic change in wild salmonids populations 25,26 .
Overall performance of simulated hybrid genotypes using both Structure and NewHybrids indicated that our dataset provided the necessary resolution to identify recent hybrids, but failed to distinguish backcross hybrids with high efficiency, as also seen elsewhere 24 www.nature.com/scientificreports www.nature.com/scientificreports/ About 15% of the hybrids in this study were found to have an aquaculture ancestry, with a similar non-random spatial distribution to those of escapees, demonstrating successful genetic hybridization and introgression of farmed escaped gilthead seabream in recipient wild populations. These wild populations showed slightly reduced allelic diversity (16WK, 16WD) or reduced Ne (17AK) in comparison to non-admixed ones. In the long run, repeated escapes of farmed fish in a wild population may cause a decline in the fitness of the recipient population due to reduced genetic variability that might affect the ability of native populations to cope with a changing environment 61,62 .
Studies of farmed seabream introgression in other regions are scarce in contrast to the extensive work conducted in Norway salmon rivers 25,26,[63][64][65] . The latest research 27 indicates significant introgression in half of the wild populations studied, with levels of introgression in excess of 10% in 27 of 109 rivers. Overall, the authors reported decreased genetic differentiation among populations over time [25][26][27] . Similar patterns of introgression level and temporal population instability have been recorded here. Due to the high larval dispersal potential and large natural variation in fish survival, impacted also by climate change and other anthropogenic factors, the degree to which these processes contribute in shaping contemporary regional population genetic patterns is more challenging than introgression estimation per se. Despite the continuous introgression of farmed seabream in the past decade 4,9 where wild populations near aquaculture facilities were most affected, it should be noted that this warm-water species shows a highly adaptive behavioural response and evolutionary capacity to prevailing conditions by increasing its dominance in coastal ecosystems and capacity for range changes towards the northern edge of its distribution. While the solid genetic base of wild and farmed strains from the eastern Adriatic is established here, as a first step toward development of mitigation strategies to prevent further erosion of genetic integrity, abovementioned adaptations require further attention.  (Table 1). All farms included in this study are active and annually produce more than 100 tons of fish. In addition, historical samples of wild gilthead seabream from Šegvić-Bubić et al. 9 , collected from the Brač Channel in 2009, were included in the analysis to test temporal stability in allele frequencies. Fish sampling was conducted in collaboration with local fishermen and farmers. Fish were measured and weighted, and the distal portion of fin-clips was clipped and stored in 96% ethanol for subsequent genetic analysis. Depending on location, sample sizes ranges from 19 to 124 individuals (Table 1).

Material and Methods
Microsatellite genotyping. Total genomic DNA from fins was extracted by proteinase K digestion, followed by a simplified DNA isolation procedure 66 . DNA quality and quantity were assessed by spectrophotometry (IMPLEN N50, Germany), and each sample was diluted to 15 ng μL −1 in DNAse/RNAse free water. Total data set containing 1586 individuals from 20 locations was characterized using two multiplex PCRs, SMsa-1 and SMsa2 (SuperMultiplex Sparus aurata), developed by Lee-Montero et al. 30  www.nature.com/scientificreports www.nature.com/scientificreports/ Genetic diversity and effective population size. The number of alleles (N) and mean effective number of alleles across loci (Ae) were calculated using POPGENE v.1.32 71 , while the inbreeding coefficient (F IS ) and allelic richness (Ar) were calculated using FSTAT v.2.3 72 . ARLEQUIN v.3.5 73 , was used to calculate observed (Ho) and expected (He) heterozygosity. Contemporary effective population size (Ne) was estimated in the program NeEstimator V2 74 using the single-sample linkage disequilibrium method for populations with a sample size over 30 individuals. Low frequency alleles ≤0.02 were excluded from the analysis to minimize potential bias caused by rare alleles. Using putatively neutral loci, parent-offspring matches were analysed using the software CERVUS 75 . The genotyping error rate was set to 1%, and three different estimates (0.01, 0.1 and 0.2) of the proportions of candidate parents sampled were evaluated. The critical LOD value was set at strict 95% confidence, where final LOD value cut-off for identifying parent-offspring pairs was always above 2. The tested dataset included allele frequencies of wild, farm-associated YOY populations.
Genetic differentiation and population structuring. POWSIM 76 was used to assess the statistical power of tests for genetic homogeneity on the applied set of markers and sample sizes. Each simulation was run 1000 times and power was determined as the proportion of simulations that Fisher's exact test detected as significant at the 0.05 level. Pair-wise and global F ST values calculation were performed in ARLEQUIN v.3.5 73 , where confidence levels were estimated by 1000 permutations of the datasets.
To infer genetic structure in the neutral microsatellite dataset and the extent of admixture between wild and farmed individuals, different statistical tests were performed, including multivariate and Bayesian based clustering analysis. Based on pairwise Nei's genetic distance matrix, a principal coordinates analysis (PCoA) was conducted using the GENALEX 6. Software 77 . Discriminant Analysis of Principal Components (DAPC) 78 from adegenetv1.4-1 79 in R was used to further characterize the potential differences between historical, present-day wild, and farmed fish. The function find.clusters and 'k-means' algorithm determined the most likely number of genetic clusters in the data and function xvalDapc was used sequentially using 1000 replicates, to determine the optimal number of principal components to retain in the discriminant analysis.
In addition to nonparametric analysis, ParallelStructure 80 an R-based implementation of the program Structure 28 on the CIPRES portal 81 was used to identify genetic clusters and to evaluate the extent of admixture among them. An admixture model was used to determine the number of population clusters (K) with a burn-in of 200,000 followed by 1,000,000 iterations with default parameters and no prior structuring information. K values from 1 to the maximal number of sampled groups were assessed, with 10 replicates per K value. The number of genetically homologous clusters (if greater than two) was determined using the ad hoc statistics ΔK 31 implemented in STRUCTURE HARVESTER 82 and visualized using CLUMPP 83 and DISTRUCT 84 . To determine the presence of hierarchical structuring of population where one genetic group can be made of several subgroups, the above procedure was repeated in successive steps for all subgroups identified with the deltaK or lnP (D) method until the lowest level of differentiation was recorded in all subgroups or when a subgroup corresponded to a single locality. Analysis of molecular variance (AMOVA) was conducted with Arlequin 3.5. using 20000 permutations on population groupings identified by the Structure analyses.
Hybrid detection. The optimal threshold values to distinguish an individual as admixed in the Structure analyses, or to assign it to a particular hybrid class in NewHybrids 85 , were identified following the approach of Vähä and Primmer (2006) 86 . In brief, we first selected 50 pure wild and pure farmed individuals whose Q-values were >0.95 in the preliminary Structure analysis using the neutral microsatellites panel. Next, we created simulated data sets of each of six genotype classes (pure parents, first and second-generation hybrids and F1 backcross with a pure parent) using the R packages hybriddetective 87 . To compare the overall performance of each software for hybrid detection, three sets containing different proportions of hybrids (15%, 33%, 66%) were constructed and analysed separately. The simulated data sets were then run through Structure using k = 2 and the above settings for the real dataset, and NewHybrids to determine the efficacy (the proportion of individuals in a group that were correctly identified), accuracy (the proportion of an identified group that truly belongs to that category) and performance (the product of efficiency and accuracy) of these programs in detecting hybrid individuals. For NewHybrids, we used a posterior probability of 50% as the threshold, while for Structure we used threshold values ranging from 0.60 to 0.95 to analyse the assignment of an individual to a specific class. Finally, the threshold value having the most suitable performance was later used to analyse the hybrid percentage in the real dataset.
We implemented NewHybrids v.1.1 85 , which provides the posterior probability that an individual belongs to one of six possible classes that differ in the extent of admixture, in our case farm, wild and hybrids (F1, F2 and backcrosses), by using Jeffreys prior probabilities and default genotype proportions, following a burn-in period of 20,000 generations and 200,000 MCMC. The simulated pure populations were included in the real population data sets using the 'z' and 's' designations indicating pure wild or farmed genotypes.
Hydrodynamic and Lagrangian modelling. Due to high gene flow observed among tuna farm-associated and YoY populations (see Results), a coupled modelling system with hydrodynamic Regional Ocean Modeling System (ROMS, www.myroms.org) [88][89][90][91][92][93] and Lagrangian individual based model Ichthyop v.3.3 (IBM, www.ichthyop.org) [94][95][96][97] was set up to infer larval supply in 2016 of the gilthead seabream adults aggregated around two main tuna aquaculture areas (16AB -Brač Island and 16AK -Ugljan Island) with the natural nursery ground areas (16JR -Raša; 16JP -Pantana; Fig. 1) and to test the correlation between genetic and transport connectivity. Two different ROMS setups were used to calculate input fields for IBM (hourly averaged current, temperature and salinity fields) for the period from 1 January, as seabream spawning regularly occurs during first two months of the year, to 15 May 2016 when juvenile sampling was conducted. The first ROMS setup (Adriatic ROMS) 54 covered the entire Adriatic Sea with a resolution of 2.5 km, while the second ROMS setup (ASHELF2 ROMS) 32 had a smaller domain encompassing the eastern coastal area of the middle Adriatic, with a resolution of 1 km