Biocomplexity in Populations of European Anchovy in the Adriatic Sea

The sustained exploitation of marine populations requires an understanding of a species' adaptive seascape so that populations can track environmental changes from short- and long-term climate cycles and from human development. The analysis of the distributions of genetic markers among populations, together with correlates of life-history and environmental variability, can provide insights into the extent of adaptive variation. Here, we examined genetic variability among populations of mature European anchovies (n = 531) in the Adriatic (13 samples) and Tyrrhenian seas (2 samples) with neutral and putative non-neutral microsatellite loci. These genetic markers failed to confirm the occurrence of two anchovy species in the Adriatic Sea, as previously postulated. However, we found fine-scale population structure in the Adriatic, especially in northern areas, that was associated with four of the 13 environmental variables tested. Geographic gradients in sea temperature, salinity and dissolved oxygen appear to drive adaptive differences in spawning time and early larval development among populations. Resolving adaptive seascapes in Adriatic anchovies provides a means to understand mechanisms underpinning local adaptation and a basis for optimizing exploitation strategies for sustainable harvests.


Introduction
A growing body of evidence indicates that marine organisms with potentially dispersive planktonic larvae can show genetically patchy distributions that may be uncorrelated with spatial complexity [1,2]. This unexpected genetic population structure has been attributed to natural populations and iv) identify environmental variables influencing the adaptive seascape of Adriatic anchovies.

Materials and Methods
Ethical statement, Sample collection and DNA extraction A total of 531 mature European anchovies were collected at 13 localities in the Adriatic between July and September 2012 and at two localities in the Tyrrhenian Sea between May and June 2013 (Fig 1A; Table 1A). Since samples of European anchovies were collected during research cruises, special ethical permits for scientific sampling were not required. Samples from the Adriatic, except the Boka Kotorska sample (KOT), were collected by midwater pelagic trawls during annual acoustic surveys between July and September 2012 in southern (July) and northern-central (September) areas of the Adriatic Sea. Tyrrhenian samples were collected with pelagic trawls during the EVATIR 2013 oceanographic research cruise that took place May-June 2013. No area under protection or special management was included in our sampling design. Permission to carry out surveys in Italian waters was granted by the MIPAAF "Ministero per le Politiche Alimentari, Agricole e Forestali", Stato Maggiore della Marina (MARISTAT), Istituto Maridrografico, Comandi Militari Marittimi di La Spezia and Taranto and the local Adriatic/Tyrrhenian Port Authorities. Sampling in the waters of Montenegro was permitted by the Local Ministry of Agriculture and Rural Development, and sampling off Albania was permitted by the Ministria Mjedisit, Pyjeve dhe Administrimit Ujrave Drejtorise Politikave te Peshkimit. Permission to conduct the acoustic survey and related biological samplings of small pelagic fish in the eastern part of the Adriatic Sea during September 2012 on board the r/v BIOS DVA was granted by the Croatian Ministry of Agriculture and Ministry of science, education and sports. Whole animals were preserved in absolute or 95% ethanol at -20°C until DNA extraction. Total DNA was extracted from gill or muscle tissue using a standard phenolchloroform protocol (gill) [24], or with the NucleoSpin Tissue Kit (muscle) (Macherey-Nagel).

Molecular markers and PCR amplification
All samples were screened at 14 previously described microsatellite loci [25][26][27][28]. Eleven loci were labeled with fluorescent dyes and multiplexed in three polymerase chain reactions (PCR) (S1 Table, S1 File). These loci were genotyped using an ABI-PRISM 3130xl Genetic Analyzer (Applied Biosystems). The remaining loci, Enja-148, Ee2-508 and Ee2-165b, were amplified individually due to moderate levels of polymorphism and were screened with a silver staining protocol [29]. Additional information about PCR protocols, amplification profiles, fragment genotyping, internal standards, silver staining adopted for the non-labeled PCR products and genotyping accuracy are provided in S1 File.

Statistical treatment of microsatellite and mtDNA data
The incidences of null alleles, allelic dropout and stutter peaks were assessed with MICRO-CHECKER 2.2.1 [30], and loci affected by null alleles were corrected with the Brookfield algorithm [31] in MICROCHECKER. In addition, the Dempster algorithm [32] in FreeNa [33] was used to estimate null allele frequencies per locus and per sampling locality. Number of alleles per locus (N A ), allelic richness (R S ), observed (H O ) and expected (H E ) heterozygosity, and the inbreeding coefficient (F IS ) were calculated for each sample using FSTAT 2.9.3 [34].
Conformation to Hardy-Weinberg proportions was tested with the exact test implemented in Genepop 4.0.10 [35] using a Markov Chain Monte Carlo (MCMC) method with 100 batches of 10 000 iterations each, with the first 1000 iterations discarded before sampling. Genepop was also used to test for linkage disequilibrium between loci by a MCMC chain executed with 1000 batches of 2000 iterations each. A sequential Bonferroni adjustment of P-values was used to account for an increase in type-I error for multiple comparisons [36].
Genetic differentiation between samples was estimated with Ɵ ST [37] and compared with divergences (F ST ) estimated with the Exclusion of Null Allele Method (ENA) in FreeNa [33]. High heterozygosity and large number of alleles (20 or more) can inflate estimates of F ST [38]. To evaluate this effect, genetic differentiation between samples was also estimated with G ST [39], G' ST [40] and D EST [41]. G ST , G' ST and D EST were calculated with GenAlEx 6.5 [42]. Genetic population structure was also examined using a Bayesian approach implemented in STRUCTURE 2.3.2.1 [43,44] and by the multivariate ordination method in DAPC [45], implemented in the adegenet package [46] in R 2.15.3 (R Development Core Team, 2009; http:// www.r-project.org).
Runs in STRUCTURE were made assuming K = 1-10 and imposing an admixture model with correlated allele frequencies. Each K value was replicated with ten independent runs of 10 6 MCMC iterations, after a burn-in of 10 5 iterations. A set of STRUCTURE simulations were  [17,18]. (b) Map showing outcomes from the BARRIER [50] and Migrate-n [54]. The three lines (solid, dotted and dashed) represent major barriers to dispersal. The arrows represent the directionality of gene flow from each locality and a map (in the right up corner) of major seasurface currents allows a direct comparison among them.
doi:10.1371/journal.pone.0153061.g001 made using i) all loci, ii) only neutral loci and iii) only candidate outlier loci. These simulations were repeated with (LocPrior mode) and without Prior information. When Prior information was included, the following clustering criterion (based on geographical partition of data) was assumed: Cluster 1 (northwest Adriatic) [ [47] in STRUCTURE HARVESTER [48] was used to determine the most likely K value from all simulations and graphical representation of multiple runs per K values were produced by the CLUMPAK Server (http://clumpak.tau.ac.il) [49].
DAPC [45] is a multivariate approach using principal components of genetic variation (PCA) that maximizes differences among groups while minimizing differences within groups. A set of 150 principal components was retained as predictors for discriminant analysis among samples.
The geographic boundaries among sampled populations were identified with BARRIER 2.2 [50], which uses Monmonier's algorithm to identify boundaries characterized by an abrupt break in gene flow (genetic barriers) and plots these boundaries on a map with sample coordinates. The strengths of putative boundaries were evaluated using a bootstrapping procedure separately for each locus. Contemporary gene flow was evaluated with an individual-assignment test in GENECLASS 2.0 [51]. For each locality, the criterion of Rannala and Mountain [52] was applied using 10 000 individuals and the Paetkau et al. [53] algorithm to detect first generation migrants.
The second method used to detect gene flow and connectivity was based on the coalescent calculation of historical migration rates between sampling locations using MIGRATE-n 3.5.1 (http://popgen.sc.fsu.edu/Migrate/Migrate-n.html). A Bayesian method was applied [54] and F ST estimates among localities were used as baseline to calculate the mutation parameter Θ (Θ = 4N e μ, where N e is the effective population size and μ is the locus mutation rate) and the historical migration rate, M (M = m/μ, where m is the immigration rate per generation). A Brownian motion model was used and mutation rates were assumed to be constant over a set of 1000 neutral loci. The MCMC procedure consisted of one long chain of 500 000 genealogies for each locus, with the firs 10 000 genealogies discarded as burn-in. The sum of M values from each location was used to determine directionality of gene flow.
The relationship between genetic and geographic distances between samples (isolation by distance) was tested using Mantel's difference matrix test implemented in Genepop (ISOLDE) [35]. Genetic distances (Ɵ ST ) between samples were compared to linear shortest sea-distances between samples with a permutation test of 10 000 iterations.
Non-neutral loci were identified with four methods that are based on different statistical assumptions. First, the coalescent approach in LOSITAN [55] detects outlier loci from the joint distribution of F ST and expected heterozygosity (H E ) under the island model of migration [56]. Runs were made with both the Infinite Allele Mutation Model (IAM) and Stepwise Mutation Model (SMM) each with 100 000 simulations. A second, Bayesian approach implemented in BayeScan 2.0 [57], uses a Reversible Jump Markov Chain Monte Carlo (RJ-MCMC) algorithm to obtain posterior distributions and is useful for small sample sizes because it reduces bias. We used default parameters of 20 independent runs with 5000 iterations each, 50 000 iterations burn-in, sample sizes of 5 000 and thinning of 10. BayeScan results were inferred following the Jeffrey's interpretation, as described in the BayeScan User Manual. Third, a hierarchical finite island model [58] in Arlequin 3.5.1.3 [59] was used to identify potential candidate loci by partitioning the dataset into population clusters as inferred from the STRUCTURE results. This method is particularly efficient for species showing significant population structure, which reduces false positives in the hierarchical analysis of genetic differentiation. Replicate coalescent runs of 100 demes were performed with 20 000 iterations. Only loci identified as outliers in all three approaches were considered as true outliers. Forth, we used the LnHR method [60], which is specifically designed to detect outlier microsatellites by evaluating changes in H E between pair of samples under the following formula: We used a threshold of ±2.58 and assumed that LnRH values under neutrality were normally distributed after data standardization (mean = 0, SD = 1). The values of 99% of neutral loci are expected to fall within this range.
A set of 13 environmental variables (listed in S2 Table) was recovered from the Sea Data-Net Climatologies Pan-European Infrastructure for Ocean and Marine Data Management web site (http://gher-diva.phys.ulg.ac.be/web-vis/). In this database, chlorophyll and oxygen values were reported as averages over 1890 and 2008, whereas salinities and temperatures were reported as averages between 1900 and 2009. These values were used as independent variables in a redundancy discriminant analysis (RDA) implemented in CANOCO 4.5 [61]. These ordinations were compared to the allelic-frequencies distributions of outlier loci to test for alleleenvironment associations. Alleles with frequencies less than 5% were removed to reduce noise, and the environmental data series was tested for a neutral distribution by a Shapiro-Wilk test. Non-neutral variables were Log 10 or Ln transformed for seascape analyses. The forward selection of environmental variables at each locus was performed with a Monte Carlo permutation test of 999 iterations (P < 0.01). When a locus was associated with more than one variable, a permutation test as above was repeated to establish the conjunct significance of these variables on the distribution of allele frequencies. The associations between allele frequencies in samples and environmental variables were further investigated with a Canonical Correspondence Analysis (CCA) in PAST 2.17 [62] with 10 000 permutations to determine significance.

Detection of multiple anchovy species in the Adriatic Sea
Based on the assumptions of Borsa et al. [21], the most important feature distinguishing shallow-water silver anchovies (E. albidus sensu Borsa et al. [21]) from open-water blue anchovies (E. encrasicolus) is bathymetry. Hence, we collected samples from shallow coastal sites (<50m) and from open-water areas (>50m). We first performed a STRUCTURE simulation to test the effect of bathymetry. Two major groups of samples were defined based on the threshold of 50 m of the bottom depth and were used as a priori information. Cluster 1 included SLO+NAD +ANC+PEA+KOT +BAA+CDG, and Cluster 2 included the remaining samples. Second, we verified the consistency of clustering samples based on their sampling bathymetry (as performed in STRUCTURE) using an AMOVA test. Finally, estimates of genetic differentiation between samples were used to test for large genetic differences, as expected between divergent species.

Microsatellite variability
The PCR success rate was high for the 14 microsatellite loci with only 1.17% missing genotypes. Although allele dropout and stuttering were minimal, null alleles were detected in 50 of 210 global tests. Null-allele corrections with the Brookfield algorithm improved the results, but 26 of 210 tests were still significant. Ten of the significant tests were for Enja148, which we excluded from the following analyses.
Locus polymorphism varied from 11 alleles at Ee2-508 to 41 alleles at Ee2-507 (S3 Table). The mean number of alleles (N A ) per sample over all loci varied from 11.4 alleles in JAB to 15.8 alleles in KOT. Allelic richness (R S ), standardized with 19 individuals per sample, varied from 8.1 (DUG) to 11.2 (CDG) ( Table 1; S3 Table). The smallest values of H E and H O were in samples NAD and RIJ, whereas the largest values were in samples SLO and KOT (Table 1; S3  Table). The smallest values of H E and H O by locus were in Ee2-165b and Ee2-508, respectively, and the largest values were in loci Ej2 and Ee2-507, respectively (Table 2).
Significant departures from the Hardy-Weinberg (HW) proportions over all loci were found in 4 samples, MNA, MNB, BAB and ANC. F IS was significant in 4 of 210 tests, and three values of F IS indicated homozygote excess (S3 Table). Three loci deviated significantly from HW proportions over all populations (Ej41-1, Eja17 and Ee2-508) ( Table 1). Only 4 of 1170 tests (78 pairwise tests X 15 localities) showed significant linkage disequilibrium, which occurred between EJ35 x Ee2-507 and Enja83 x Ee2-507 in the NAD sample, EJ27.1 x EJ35 in the DUG sample and Ee2-507 x Eja17 in the RIJ sample.  Table). STRUCTURE simulations of neutral loci, with and without LocPrior information, yielded two major clusters (S1A,S1B and S1E Fig): i) DUG+JAB+RIJ (Northeast Adriatic samples) and ii) the remaining samples. However, a high probability of three major clusters emerged using only candidate outlier loci without the LocPrior information (S1C Fig) and using all loci with the LocPrior function (Fig 2A; S1D Fig). These three clusters were largely represented by i) DUG+JAB+RIJ (Northeast Adriatic samples), ii) NAD (Northern Adriatic sample) and iii) the remaining samples. Other possible substructures yielding 5 and 8 clusters appeared in the simulations, but with low q values (S1A-S1F Fig).

Microsatellite genetic population structure
In the DAPC analysis, the first 2 of the 150 Discriminant Axes (DA) were significant. These DA's defined three major genetic Clusters: i) DUG+JAB+RIJ (Northeast Adriatic samples), ii) NAD (Northern Adriatic sample) and iii) all remaining samples (Fig 2B). The first DA separated mainly the Northeast Adriatic Cluster from the NAD Cluster, while the second DA largely separated those two Clusters from a third, which included the remaining samples (Fig 2B).
A total of 38 of 531 individuals were identified as first generation migrants. The analysis of migration rates and first generation migrants between clusters revealed a higher intra-cluster migration than inter-cluster migration in the two or three genetic groups. The directionality of gene flow was largely from populations MNA and ANC (S5A Table). Estimates of Θ from the historical gene flow ranged between 0.003 in CDG to 0.046 in JAB (S5B Table). Migration rates (M) ranged between 18.50 (M JAB!NAD ) and 944.01 (M DUG!PEA ). The directionality of historical gene flow suggested a pattern that was similar to sea-surface circulation (Fig 1B). The main sources of migrants were DUG, ANC and CDG (S5C Table).
Additional approaches resolved other features of geographic structure. First, BARRIER identified at least three major allele-frequency discontinuities, indicating a heterogeneous pattern of partial isolation between locations. These results indicated two strong barriers, one between NAD and nearby sampling areas ( Fig 1B, barrier a, dotted line) and a longitudinal barrier across the Mid Adriatic Pit that divided Adriatic anchovies into east-west groups (Fig 1B,  barrier b, solid line). In addition, a north-south barrier that divide the Mid Adriatic Pit was detected (Fig 1B, barrier c, dashed line). This barrier was stronger in the east than in the west.
Even though differences appeared between populations, no isolation-by-distance over the length of the Adriatic was detected with Mantel's test between pairwise ϴ ST and geographical distances (P > 0.05; a = 0.0165172, b = 0.00000108) (data not shown).
The Redundancy Discriminant Analysis (RDA) identified significant correlations between allelic frequencies at four outlier loci (Ee2-507, Ee2-407, Ee10 and EJ2) and four environmental variables, including temperature at sampling depth (T samp ), salinity at sampling depth (S samp ) and at 10 m (S 10 ) and dissolved oxygen concentration at 10 m (Oxyg 10 ). At Ee2-507, alleles 279 and 291 were positively correlated and alleles 259 and 283 were inversely correlated with S 10 (both P < 0.019) (Fig 3A). Significant correlations appeared for Ee2-407 with both Oxyg 10 and S samp (P < 0.002); allele 169 varied positively and allele 147 inversely with S samp , and alleles 163 and 171 were positively and allele 159 was inversely correlated with Oxyg 10 (Fig 3B). At locus Ee10, alleles 224, 232 and 258 were positively and alleles 226 and 256 were inversely correlated with T samp (P < 0.003) (Fig 3C). Finally, at EJ2, alleles 196, 200 and 210 were positively and alleles 178 and 186 were inversely correlated with T samp (P < 0.033) (Fig 3D).  (Fig 2A). Graphical outcome of DAPC plot (Fig 2B). Both the first and second axes of the CCA ordination of samples with four environmental variables were significant (P < 0.05 and P < 0.001, respectively). The results indicated a positive association for allele frequencies in EMB, EMC and EMD and a negative association for allele-frequencies in sample EMA with salinity S samp and S 10 (Fig 4). Samples SLO, KOT and PEA were positively and samples SPE and CDG from the Tyrrhenian Sea were negatively associated with temperature at sampling depth (T samp ). Finally, SPE and CDG were directly associated with dissolved oxygen content in the first 10 m (Oxyg 10 ) and inversely associated with T samp .

Detection of multiple species of anchovies in the Adriatic Sea
The STRUCTURE analysis using the LocPrior function to cluster sampling sites by bathymetry showed three genetic clusters, as in the other STRUCTURE analyses (S1G Fig). In general, the pattern did not show a clear association with bathymetry as postulated by Borsa et al. [21]. Furthermore, the AMOVA test partitioned by ocean depth of 50 m, indicated that the largest variance was among individuals within populations (97.86% of total variance; S7 Table). Although variation among populations within groups (F SC = 0.0017) was significant (P < 0.05), it explained only a small amount of the total variation (1.79%). Divergence among groups (F CT = 0.0012) was not significantly different from zero (P > 0.05).

Discussion
The results of our survey of microsatellite DNA variability among populations show that while anchovies in the Adriatic are genetically subdivided among regions to some extent, none of these populations represent separate species, as previously proposed. The mechanisms producing divergence between populations include both demographic processes influenced by ocean currents and fronts and natural selection that leads to locally adapted populations. Before exploring these factors in more detail, we offer the following evaluation of the study. While the sample sizes used in this study were not large, they provided sufficient statistical power to detect population differences. Statistical power was enhanced by the use of a large number of microsatellite loci encoded in nuclear genes. These markers provided complementary insights into population structure and local adaptation. We additionally used several approaches to the analyses of the microsatellite data that illuminated various dimensions of population structure.
How many anchovy species in the Adriatic Sea?
Borsa [63] and Borsa et al. [21] postulated that two species of anchovy inhabit the Adriatic Sea on the basis of north-south allele-frequency gradients reported by Bembo et al. [22] and on morphological parallels with two putative anchovy species in the Gulf of Lyon [21,63]. The recently described estuarine species, E. albidus Borsa, Collet et Durand, 2004, resembled shallow-water anchovies in low-saline areas of the northern Adriatic [21]. The larger, bluish anchovies in the Adriatic resembled European anchovies elsewhere in open, deep waters with salinities greater than 38 psu, and were postulated to be E. encrasicolus [21]. However, our analysis of nuclear markers in anchovies from throughout the Adriatic did not show genetic differentiation between coastal and deep-water populations that could be interpreted to indicate multiple species. Nevertheless, it might be argued that the genetic differences in the northern Adriatic Sea support the occurrence of E. albidus in the low-saline plume of the Po River (NAD). Although these anchovies share some ecological and morphological traits with putative E. albidus in the Gulf of Lyon, Po River anchovies are genetically most similar to oceanic anchovies in the central Adriatic. The level of genetic differentiation between NAD and the other samples suggests a substantial genetic separation, but this level of divergence is not large enough to indicate the presence of two divergent species. The color and morphological variability observed among anchovy populations in the Adriatic Sea (and in the Gulf of Lyon) likely reflects plastic phenotypic responses to environmental variability in anchovies that are typical of anchovies in the Mediterranean Sea [64].

Genetic variability among populations
The small values of F ST and ϴ ST between Adriatic populations are typical of divergence between other population groups of anchovies or between populations of other pelagic fishes. The overall levels of genetic differentiation between populations in the Adriatic Sea were similar to those between anchovy populations in the eastern Atlantic [65,66]. As in many marine species, gene flow in Adriatic anchovies occurs by larval dispersal in currents, or by adult movement at low levels over long periods, and prevents the accumulation of genetic differences among populations [67]. Nevertheless, significant differences between some samples indicate that Adriatic anchovies are not entirely panmictic.
The analysis of neutral makers with STRUCTURE resolved two major population groups: populations in northeastern coastal areas (DUG, JAB and RIJ) and the remaining populations in the Adriatic and Tyrrhenian seas. The genetic uniqueness of northeastern anchovies may be due to adaptation to habitats in the numerous coves and embayments and around offshore isles [17,18,68]. The northeastern Adriatic appears to be conducive to spawning and early larvae growth and is an area where spawners converge [69].
A north-south genetic discontinuity in the northern Adriatic Sea was previously reported by Bembo et al. [22], who found an allozyme-frequency gradient between northern and southern populations. Bembo et al. [22] focused on the possible role of depth, defining the 40-50 m sill as the population boundary, rather than on other environmental variables. In our analyses, depth did not appear to be an important isolating variable, as we found no genetic partitioning of microsatellite variability by ocean depth in areas where samples were collected.
Instead, the oceanic front between the Istrian Peninsula and Mid Adriatic Pit may explain, in part, the north-south genetic discontinuity by limiting gene flow. The anticlockwise gyre off the Po River delta may also act as a barrier to gene flow [17,18]. Although both the front and the gyre are ephemeral, they are strongest in autumn [17,18] at the end of the anchovy spawning season when larvae and juveniles are most abundant [19,68]. Larval retention by ocean fronts and gyres is a well-known isolating mechanism in other pelagic fishes [70,71]. Hence, these oceanic mechanisms may limit the movement of larvae and juveniles and prevent mixing between larvae and juveniles from neighboring spawning sites. The effect of small coastal eddies limiting the dispersals of early life-history stages was also invoked by Borrell et al. [66] to explain genetic differences between nearby populations in the Bay of Biscay.
Individual assignment showed that about 10% of the fish in our samples represented firstgeneration migrants from other populations. However, patterns of dispersals revealed by these tests have to be interpreted with caution, because small samples sizes and finite geographical sampling introduce some randomness into the results. The results also represent only a snapshot of a dynamic process that may change over time. Nevertheless, the results indicate that dispersals between populations are asymmetrical: some populations act as 'sources' and others act as 'sinks'. The extent that populations receiving migrants depend on immigration to persist has to be evaluated with long-term abundance data. With these caveats, most migrants originated from only a few source locations (mostly MNA and ANC samples) and moved to several nearby sink locations following a pattern of dispersion that is similar to major sea-surface currents in the area. The evaluation of historical patterns of gene flow using Migrate-n suggests a similar scenario for both the magnitude and pattern of connectivity. This pattern of migration may also be related to the lack of divergence between Adriatic and Tyrrhenian anchovies because of reciprocal gene flow between the two seas. These results are in line with a recent analysis using SNPs markers showing genetic similarities between Adriatic and Western Mediterranean population of anchovies [72].

Local adaptation
While the small genetic distances between samples indicate high levels of gene flow, as expected for a pelagic species, the pattern of geographic differentiation and the associations of microsatellite loci with environmental variables indicate that gene flow may be countered by adaptation to local conditions. For example, the timing of spawning, relative to seasonal bursts of productivity can have a strong effect on year-class survival [73,74]. The persistence of genetic population structure in Adriatic anchovies may indicate differential reproductive success or reduced fitness of migrants reproducing outside their natal spawning areas. If so, adaptation to local habitats explains some of the genetic population structure observed in Adriatic anchovies.
Our inferences about natural selection are based on non-coding microsatellite loci, which are unlikely to be the direct focus of selection, except when microsatellite alleles modify the functions of genes in which they are embedded [75,76]. Alternatively, microsatellites may show signals of selection when they are linked to genes under selection. The strength of a signal from this 'hitchhiking' effect depends on the strength of selection, recombination rates and genetic drift [77,78]. Our analyses identified 5 of 13 loci (38.46%) as candidate outliers, and two loci under directional selection (Ee10 and Ee2-507) (15.38%), indicating the potential role of natural selection in driving differentiation among populations. The proportion of candidate outliers was large, but consistent with other studies that showed between 5-15% outlier loci under directional selection [79,80].
Several life-history traits can interact to facilitate local adaptation. For instance, it is well known that this species has a high skew in reproductive success, which determines the level of genetic diversity in each cohort. Sweepstakes reproduction can promote selection at local scales, especially when coupled with large demographic variation in census size. During population bottlenecks, the reduction in effective population size (N e ) can produce transient genetic sweeps that lead to locally adaptive genotypes [81]. A recent study showed that the effective sizes of anchovy populations in the Adriatic are several orders of magnitude smaller than census sizes [82]. These small effective sizes appear to be in part a legacy of a severe reduction in population abundances from overharvesting in the 1980s. Small effective population sizes in anchovies suggest that genetic drift, together with repeated selective sweeps, can mold the anchovy genome, resulting in local adaptation to major environmental drivers.
To further investigate the molecular footprints of selection, we tested candidate outlier loci for associations with environmental variables suspected to be important in the physiologies and ecologies of local anchovy populations. Four of the five candidate outliers were significantly correlated with environmental variables. Interestingly, these variables, especially salinity and temperature, influence larval development, survival and growth [19,68], as well as, adult maturation and reproductive success [83]. Associations with salinity and temperature variation explained patterns of genetic variation. The northeastern group of anchovies showed a positive association with higher salinities, whereas anchovies in the Po River Plume showed an inverse association with salinity. The waters in these adjoining areas showed the largest difference in salinity among the sites sampled in our study, with the Po River plume having the lowest salinity, and waters of the eastern Adriatic along the northern Croatian coast having the highest salinities in the whole Adriatic Sea [17,18]. Natural selection from strong differences in salinity among some areas may explain the genetic divergence between populations in the northern Adriatic Sea, despite substantial levels of gene flow between these areas.
In the CCA analysis, a set of samples from shallow Adriatic coastal waters (SLO, KOT and PEA) were positively associated with temperature variability, whereas samples from the Tyrrhenian Sea (SPE and CDG) were inversely associated with temperature variability. This is consistent with high seasonal fluctuations in temperature in northern and western Adriatic coastal areas, in contrast to greater thermal stability in the Tyrrhenian basin [68,84]. Temperature also appears to be an important selective factor among populations of anchovies in the Eastern Atlantic [85].
Adriatic anchovies are also influenced by variability in dissolved oxygen, which is inversely correlated with temperature, because oxygen saturation in water decreases at higher temperatures. Populations in the Tyrrhenian Sea (SPE and CDG) were positively associated with dissolved oxygen content, but coastal populations in shallow waters (SLO, KOT and PEA) were inversely associated with oxygen content. The dissolved oxygen content may interact with temperature to define local selective constraints. Higher water temperatures in summer generally lead to a drop in dissolved oxygen, producing hypoxic and anoxic conditions [86], especially in the Adriatic Sea.
Genetic discontinuities among populations of anchovies may be related to the effects of environmental variables on the growth and survival of early developmental stages, as well as on adult reproductive success. The complex genetic patterns among populations of marine organisms have only recently been correlated with the influence of environmental variables on the genomes of several marine species [11]. Numerous studies show that environmental variables affect adaptive responses in marine fishes on local (European hake [15]) and oceanic scales (Atlantic anchovies [85]; Bluefin Tuna [14]).

Conclusions
An important result of this study indicates that anchovies throughout the Adriatic belong to a single species, E. encrasicolus. Nevertheless, the individuals collected in the study area, which covered the entire Adriatic basin and part of the Tyrrhenian Sea, were partitioned into two major populations which were sometimes further sub-structured. Adriatic anchovy subpopulations show a patchy distribution of genetic variability that appears to reflect selective processes mediated by several environmental variables, including salinity, temperature and dissolved oxygen content. The relationship between environmental variables and the distribution of genetic differences among population groups is consistent with adaptation to local environmental conditions. Salinity, in particular, appears to be important in structuring northeastern Adriatic anchovy populations. Our geographical survey of microsatellite DNA variability reveals a more complex genetic structure among Adriatic anchovy populations than was previously detected, but does not confirm the presence of a second species of Engraulis. Further resolution of this biocomplexity will provide a basis for understanding adaptive mechanisms and a basis for optimizing management strategies for sustainable harvests.
Supporting Information S1 Fig. Graphical plotting of simulations carried out using STRUCTURE [43,44] under various settings. The simulations were run using i) all 13 loci, ii) only neutral loci and iii) only candidate outliers with and without the assistance of prior information (S1A-S1F Fig). Each simulation shows the mean logarithmic probability with standard deviation (SD) obtained for each K tested (1 < K < 10; figure on the left) and the delta K probability obtained after the use of Evanno method [48] (2 < K < 9; figure on the right). Additionally, a simulation was carried out dividing sampling sites based on bathymetry using a threshold of 50 m to detect population structure associate with potential occurrence of different anchovy species in the covered area