Environmental heterogeneity promotes microgeographic genetic divergence in the Mediterranean killifish Aphanius fasciatus (Cyprinodontidae)

Environmental stress can promote evolutionary forces and genetic divergence, but at which microgeographic scale divergence may arise, even in the presence of gene flow, remains poorly known. We studied the effects of eutrophication in a saltwork over a period of 6 years on the gene pool of a local population of the Mediterranean killifish Aphanius fasciatus, a species resistant to extreme conditions. We hypothesised that during eutrophication, the environmental stress may have acted differently along a gradient of salinity and oxygen concentration in promoting evolutionary forces, generating divergence over a small spatial scale (2 km), despite gene flow. We analysed 24 allozymes in three temporal samples each composed of four spatial sub-samples, collected along the gradient, during eutrophication (2003 and 2005) and after a recovery project (2008). The results suggest that eutrophication promoted natural selection, originating a genetic cline on one locus (adenosine deaminase) significantly linked to salinity and oxygen concentration. Together with selection, both genetic drift and gene flow contributed to shaping the genetic structure under stress by further promoting the genetic heterogeneity and giving rise to deficits of heterozygotes as a secondary effect of the divergence. Environmental stress thus increased diversification, with the effects of selection and drift prevailing on gene flow. When environmental stress was relaxed (2008), allele and genotype frequencies became homogeneous, likely because under less extreme conditions the gene flow prevailed. These results improve our understanding of microgeographic divergence, and highlight the role of environmental stress in moulding microevolutionary dynamics and genetic patterns of animal populations.


INTRODUCTION
Environmental stress can exert adverse effects on biological systems, but could also promote rapid evolutionary changes (Parsons 2005;Nevo 2011). Stressful environments can help to overcome evolutionary stasis by imposing natural selection and microgeographic adaptation (Hoffmann & Hercus 2000;Dvornyk & Jahan 2012;Le Corre & Kremer 2012;Chapman et al. 2013;Richardson et al. 2014). Periodic stressful conditions can promote the persistence of genetic variation through heterozygous advantage or fluctuating selection (Gillespie & Turelli 1989), and could even increase mutation and recombination rates in certain organisms (Yauk & Quinn 1996;Lamb et al. 1998;Somers et al. 2002). In contrast, if extreme conditions persist, the demographic and effective population size may be reduced, thus lowering genetic variation via genetic drift (Hoffmann & Hercus 2000). Environmental stress can thus enhance evolutionary forces by exerting contrasting effects on a population's gene pool.
The existence of antagonistic and overlapping interactions makes it difficult to disentangle the effects of environmental pressures on single evolutionary forces (Joseph et al. 2012), with the consequence that the genetic patterns observed in natural populations may be hard to interpret. Besides the old debate between neutralism and selectionism, recent empirical research has depicted a framework in which both deterministic (natural selection) and stochastic forces (mutation, genetic drift and gene flow) play their own roles in molecular evolution and in promoting genetic divergence, and the debate has thus shifted to evaluate the relative influence of each factor in shaping natural genetic patterns (Riginos et al. 2002;Wagner 2008). Divergence and local genetic adaptation can occur when drift or selection overcomes the homogenising effect of gene flow, which is assumed to be high, thus preventing differentiation, at a small spatial scale (Hendry et al. 2001;Björklund et al. 2007;Adams et al. 2016). Research has only recently challenged this assumption, so the scale and environmental conditions under which genetic divergence may arise despite gene flow largely remain to be explored (Maltagliati et al. 2003;Whitehead et al. 2011;Richardson et al. 2014).
Temporal studies of microevolutionary dynamics in populations that live in confined, heterogeneous, and unstable environments can be particularly suited to address this issue, and they may help to clarify both the relative contributions of evolutionary forces and their possible interactions in shaping the genetic structure of populations (Kawecki & Ebert 2004). Indeed, in isolated local populations, the disruption of gene flow favours the maintenance of different acquired genetic and phenotypic traits (Hoffmann & Hercus 2000) making microevolutionary dynamics easier to interpret with respect to local environmental events. Furthermore, in heterogeneous habitats, local sub-populations or demes will experience limiting conditions, and environmental changes over space and time are likely to be reflected by parallel variations at the genetic level (Vergara-Chen et al. 2013;Rodrigues et al. 2015).
Here, we studied the effects of environmental changes on the genetic structure of a local population of the Mediterranean killifish Aphanius fasciatus (Valenciennes 1821) over a period of 6 years. We conducted our study in the Tarquinia saltworks (Italy), a small area characterised by salinity and hypoxia gradients (increasing southward), with habitat conditions that underwent temporal changes from deep eutrophication to a recovery of typical conditions (Cimmaruta et al. 2010a;Bellisario et al. 2013). A. fasciatus has a short generation time, high reproductive rate, rapid population turnover and highly isolated local populations (Leonardos & Sinis 1998;Maltagliati & Camilli 2000;Pappalardo et al. 2015). Although resistant to extreme conditions, A. fasciatus may be affected by highly stressful conditions, such as high temperature and salinity (above 24°C and 70 ppt, respectively) that are regularly reached in the Tarquina saltworks, as they may cause the discharge of oocytes and affect its fecundity (Leonardos & Sinis 1998). Accordingly, our field observations suggest that A. fasciatus became progressively less frequent in the southern area during summer, where the high salinity and hypoxia levels can lead to the death of several individuals. Thus, the study area and the species possess the features potentially allowing divergence to arise over a small scale, because of the strong and heterogeneous environmental stress. In addition, previous studies showed that: (i) the population underwent significant genetic erosion associated with both a lower effective population size (N e ) and genetic drift during eutrophication, whose effect almost ceased when conditions improved ); (ii) allozymic analysis suggested a role for natural selection in shaping the gene pool under stressful conditions, particularly in relation to adenosine deaminase (ADA), as this locus presented a significant deficit of heterozygotes during eutrophication ); (iii) allelic and genotypic frequencies at the ADA locus, correlated with dissolved oxygen of the sampling points in the southern, most stressed area of the saltwork (Angeletti et al. 2009).
Since the gradients of increasing salinity and hypoxia became steeper and more sudden during eutrophication (Bellisario et al. 2010(Bellisario et al. , 2013, we predicted that such environmental stress could overcome the homogenising effect of gene flow, causing shifts in allele frequencies as a result of different intensities of genetic drift and selection along the gradient. In particular, we expected: (i) a general genetic divergence between groups of fishes collected along the gradient, due to drift and selection; (ii) a clinal, condition-dependent distribution of frequencies at the ADA locus following the environmental gradient of salinity and oxygen concentration, due to selection; (iii) the disappearance of genetic divergence after the environmental recovery, as a result of stress mitigation, combined with gene flow. To test these hypotheses, we analysed allozymes in three temporal samples of A. fasciatus, each composed of four spatial sub-samples, collected during eutrophication (2003 and 2005) and after the recovery (2008).

Study area and environmental conditions
The Tarquinia saltworks (Fig. 1) are a coastal habitat, mostly hyperhaline, located along the Tyrrhenian coast of northern Latium, central Italy (42°12ʹN, 11°43ʹE). It is a semi-artificial, patchy ecosystem composed of nearly 100 shallow ponds extending along the coast for about 2500 m and covering an area of 135 ha. The inlet of water is provided by water pumps through a single connection with the sea located north of the area; the ponds are connected with each other either directly or by a system of drainage channels, and become progressively shallower further away from the seawater channel. The water temperature varies seasonally, ranging from 5 to 28°C, but it is rather homogeneous between the ponds. Conversely, the salinity gradient can range from 35‰ in the northern ponds up to 250‰ in the southern ones, which influences both the abiotic parameters and the biological communities of plankton and benthos (Cimmaruta et al. 2010a). Species richness decreases from north to south, and the benthic community is progressively dominated by members more resistant to hyperhaline and hypoxic conditions, such as chironomid larvae (Bellisario et al. 2010(Bellisario et al. , 2013.
Studies of macrobenthic communities, which are well-known indicators of water quality in coastal lagoon ecosystems (Lardicci et al. 2001), have shown significant alterations in the aquatic environment since 1998 as a result of changes in the human management of the site. After cessation of salt production in 1997, the habitat suffered progressively severe degradation because of the lack of maintenance activities, and scarce water flow, together with the increasing sedimentation of organic material, causing eutrophication and hypoxic conditions between 2003 and 2005 (Cimmaruta et al. 2010a). Notably, macrobenthic and environmental data showed that the pre-existing gradients of increasing salinity and hypoxia became steeper and more sudden during eutrophication, thus limiting the presence of less-tolerant species to the northern ponds (Bellisario et al. 2010(Bellisario et al. , 2013.
This situation changed in 2006 when interventions were performed to prevent further degradation of the ponds, within the framework of a European habitat restoration project (LIFE02NAT/IT/8523). The project aimed to restore water volume and circulation by partially removing organic matter and banked up mud, and by repairing the water pumps and seawater- inlet channel. These actions led to a complete recovery of the richness, structure and distribution of the macrobenthic community through mitigation of the general environmental conditions and increases in the mean oxygen levels in the water (Bellisario et al. 2010(Bellisario et al. , 2013.

Study species and sampling
A. fasciatus is a cyprinodont fish living in coastal brackish waters, lagoons and salt marshes throughout the central and eastern coastal zones of the Mediterranean Sea (Cimmaruta et al. 2010b). Although it is a mobile species within this habitat, local populations are strongly isolated and characterised by low migration rates and gene flow as a result of their discontinuous habitat, because the adults do not actively disperse in seawater, and because the species lays benthic eggs and does not have larval dispersal stages (Leonardos & Sinis 1998;Pappalardo et al. 2015). Reproduction takes place from early spring to the end of summer with high reproductive rates (Maltagliati & Camilli 2000). In the Tarquinia saltworks, A. fasciatus lives in all the ponds with salinities up to 180‰, excluding a few ponds previously used for salt collection.
We collected three temporal samples of A. fasciatus (total n = 1270) using fish traps. Traps were monitored once a day and no deaths occurred. Immediately after capture, the individuals were euthanised by immersion in an anaesthetic solution (10 mg/L MS-222 for light anaesthesia followed by 500 mg/L MS-222 for euthanasia) and then transported to the lab and frozen at -80°C for subsequent analyses. Specimens were collected in 2003 (n = 321), when abandonment of the saltwork in 1997 had caused a substantial shift towards extreme conditions and eutrophication; in 2005 (n = 469), just before the start of the recovery activities; and in 2008 (n = 480), 2 years after the end of the restoration project. Each sample consisted of four sub-samples collected across the whole north-south gradient, from less extreme to more extreme conditions of high salinity and hypoxia. The sampling scheme was thus designed to take account of environmental heterogeneity and its changes over time. The details of the 12 subsamples are shown in Fig. 1 and Table 1.

Data analysis
Preliminary analysis and general statistics. A preliminary analysis was performed to check for the influence of uneven sample sizes using the program HP-Rare (Kalinowski 2005), which compares the number of observed alleles with the number of alleles computed by the rarefaction approach, standardised for the smallest sample size (n = 44). BIOSYS-2 software (Swofford & Selander 1999) was used to obtain allelic and genotype frequencies, and thus to check the distribution of frequencies during eutrophication (2003 and 2005) and after recovery (2008). Genotypic linkage disequilibrium and within-population significant departures from Hardy-Weinberg equilibrium (HWE) were tested using Fisher's exact test implemented in GENEPOP 3.4 (Raymond & Rousset 1995; http://www.cefe.cnrsmop.fr).
Genetic divergence and clinal distribution at the ADA locus. To verify the hypothesis that the local population is structured in sub-populations only during eutrophication as an effect of stress, the partitioning of genetic diversity and its significance across the four samples for each year (2003, 2005, and 2008) were analysed for variable loci using F-statistics (Weir & Cockerham 1984) in the computer package ARLEQUIN 3.5 (Excoffier & Lisher 2010). In addition, to test for allelic and genotypic differentiation across the four samples within each sampling year, the distributions of allelic and genotypic frequencies were checked using exact probability tests for all samples and all pairs of samples, implemented in GENEPOP 3.4 (Raymond & Rousset 1995; http://www.cefe.cnrsmop.fr). An unbiased estimate of the P value for Fisher's exact test was performed for each locus and over all loci, applying the Bonferroni correction to the significance threshold for multiple comparisons. The same test allowed us to check the significance of the divergence and of the clinal distribution of allelic frequencies at the ADA locus, expected during eutrophication.
Exploring the causes of the divergence. To explore the causes of the possible divergence among the four spatial samples within each sampling year, we performed the following analyses.
To verify whether the divergence depends on the geographic distances among the spatial samples, according to an isolation by distance model (IBD ;Wright 1943;Kimura & Weiss 1964;Slatkin 1993), we tested the correlation between the genetic F ST distance matrix and the geographic distance matrix of the sampling points, for each sampling year. To this end, we performed a Mantel test (10,000 permutations) using the ISOLDE program, as implemented in GENEPOP 3.4 (Raymond & Rousset 1995; http://www.cefe.cnrs-mop.fr).
To check whether natural selection may have generated the divergence at the ADA locus or at other polymorphic loci, the selection detection workbench LOSITAN (Antao et al. 2008) based on the F DIST F ST outlier methods of Beaumont and Nichols (1996) was performed for the four spatial samples from each year to detect candidate loci for directional or balancing selection. For all runs, 75,000 simulations were generated using the 'neutral mean F ST ' and 'force mean F ST ' options to increase the reliability of the mean F ST . Simulations were performed under the infinite alleles mutation model (IAM; Kimura & Crow 1964), which is considered to be the model by which allozymes evolve.
To verify the possible contribution of genetic drift in generating the divergence, through a drastic reduction of the population size, the software BOTTLENECK (Cornuet & Luikart 1996;Piry et al. 1999) was used to infer recent bottlenecks in the four groups of each temporal sample. We applied the sign test (Cornuet & Luikart 1996) and Wilcoxon's signed-rank test (one-tailed for heterozygosity excess; Luikart et al. 1998) under the IAM model.
Exploring the causes of heterozygote deficits. In addition to selection against heterozygotes, the deficits of heterozygotes can be due to both inbreeding (i.e., the mating of close relatives; Wright 1921) and the Wahlund effect (i.e., samples composed by a mix of individuals with heterogeneous allelic frequencies; Castric et al. 2002).
To test for possible inbreeding, we calculated the relatedness coefficient (COANCESTRY software; Wang 2002Wang , 2011, whose values range from 1 to − 1, indicating the percentage of alleles shared among individuals (Karamanlidis et al. 2012). Positive values increase with the relativeness of the breeding individuals, so finding values approaching 1 would strongly support that inbreeding occurred in our samples, causing the deficit of heterozygotes. We tested the significance of pairwise differences in average relatedness between the four spatial samples of each year, using the bootstrapping method implemented in the software COANCESTRY.
If the heterozygote deficits were due to a Wahlund effect, we would expect the presence of a subtle sub-structure within the saltwork. To test for this hypothesis, we used an individual-based assignment test to calculate the probability that the studied individuals come from different 'source sub-populations'. The software PartitionML (http://www.genetix.univ-montp2.fr/partitionml.htm) was used to calculate the assignment probability for each individual based on the genotypic data, then allocating each specimen to the population with the higher likelihood (Paetkau et al. 1995;Rannala & Mountain 1997).
Relationships between the genetic structure and environmental variables. We developed generalised additive models (GAMs) as described by González-Wangüemert et al. (2009) and González-Wangüemert and Vergara-Chen (2014), to test the possible relationship between salinity and oxygen concentration, whose values showed a gradient from the northern to the southern part of the saltwork, and the genetic structure of the population. We first considered all of the samples, then we repeated the analysis considering only the samples collected during eutrophication (2003 and 2008), as we assumed that the correlations might emerge more clearly during stress. Samples were spatially clustered using principal component analysis (PCA), using the allele frequencies of samples as variables. The GAM models were then developed using the first principal component of the PCA as a dependent variable and the minimum, maximum and mean values of salinity and oxygen concentration as independent variables. We assumed respectively the dependence (i.e. salinity influences oxygen concentration) and independence (i.e. salinity and oxygen concentration are independent) of those variables. Values of salinity and oxygen concentration (Table 2) were recorded monthly in the water column during the 3 years, using a refractometer and a field probe (Oxyguard, Denmark). PCA analysis and GAM models were performed using 'ade4' (Chessel et al. 2004) and 'mgcv' (Wood 2006) packages from R statistical software (R Development Core Team 2013).

Preliminary analysis and general statistics
The number of alleles computed by the rarefaction method was consistent with the observed number with no significant differences between the number of observed and computed alleles, either between all loci or at the single locus level (χ 2 test, both P ≥ 0.99). Thus, samples can be considered representative and comparable to each other, despite their unequal sizes.
Eleven of the 24 analysed loci were monomorphic in each sampling zone for the three temporal samples. The remaining 13 polymorphic loci had two to five alleles Table 2.
Maximum, minimum and mean values (± SD) of oxygen concentration (mg/L) and salinity (‰), recorded in the sampling ponds during the 3 years of the study.  Table A1), with the most common allele (*100) occurring at frequencies of 0.50-0.80 in most samples (LDH-1, ADA, MPI, GPI-1 and PGM-1) or > 0.80 (G3PDH-1, MDH-1, G6PDH, EST-3, PEP-B, ACON, GPI-3, PGM-2). There was no evidence of linkage disequilibrium following Fisher's exact test, and the genotypes at each locus could therefore be considered independent from genotypes at other loci. The spatial samples from 2003 and 2005 showed a total of five heterozygote deficits and one heterozygote excess, at the following loci: ADA (three deficits), GPI-1 (one excess), GPI-3 (one deficit), and PGM-1 (one deficit). Moreover, multilocus tests detected a significant heterozygote deficit in the 2003 north sample (Table 3). Notably, the deficits occurred mainly in the northern samples, and no deviations from HWE emerged in 2008.  Table 4 shows the allelic and genotypic differentiation among the four spatial samples for each year (GENEPOP 3.4 exact probability test; Raymond & Rousset 1995). Similarly to the previous analysis, the test revealed signs of divergence only for the samples collected during eutrophication. Significant differences between allelic and genotypic frequencies at the single locus level emerged from 2003 (loci ADA and GPI-3) and 2005 (loci ADA, PGM-1) spatial samples. Furthermore, the analysis computed over all loci indicated that the samples from 2005 were structured in genetically differentiated groups. Conversely, the samples from 2008 did not show any allelic or genotypic differentiation. Fig. 2 shows the spatial distribution of allele frequencies at the ADA locus during the 3 sampling years. As expected, this locus, which presented the highest divergence during eutrophication (2003 and 2005; see Table 4), in the same years also showed a clinal distribution of frequencies, as ADA*100 increased from the northern to the southern area (from 0.62 to 0.83). This increase is significant, according to the previous  Exploring the causes of the divergence

Genetic divergence and clinal distribution at the ADA locus
The genetic structuring of the population did not conform to an IBD model, both during and after eutrophication, as no significant correlation between the genetic F ST distance matrix and the geographic distance matrix of the sample points emerged, for any year.
With regard to selection, coalescent simulations (LOSITAN; Beaumont & Nichols 1996;Antao et al. 2008)  Finally, the sign test (Cornuet & Luikart 1996) performed by the software BOTTLENECK (Cornuet & Luikart 1996;Piry et al. 1999) under the IAM model revealed signs of recent bottlenecks for the 2003 southernmost sample (P = 0.049). No signs of recent bottlenecks emerged for the other samples. Table 5 shows the mean values of the relatedness coefficient (Wang 2002) and the proportion of individuals of each spatial sample assigned to a hypothetical second source population. Relatedness values range from -0.111 to 0.149 and tended to be lower in the northern samples, both during (2003 and 2005) and after eutrophication Table 5.

Exploring the causes of heterozygote deficits
Relatedness coefficients (RC; Wang 2002) and proportion of individuals assigned to a different source population (DSP%) (PartitionML; http://www.genetix.univ-montp2.fr/partitionml.htm), for the four spatial samples within each sampling year. Different letters in parentheses indicate significant pairwise differences in average relatedness of spatial samples, within each year; numbers in parentheses indicate heterozygote deficits observed in the samples (see Table 3).
Year 2003  (2008). These values, well below 1, indicate that no inbreeding is occurring in any of our samples. The lowest relatedness values are always found in the northern samples, which thus included significantly less-related specimens than the other groups did (Table 5). The assignment test supported the allocation of the studied individuals to more than one sub-population: a high percentage of the specimens was assigned to a different source sub-population, and in the samples from the northern area this percentage reached 79% in 2003, in coincidence with the higher number of heterozygote deficits.

Relationships between the genetic structure and environmental variables
Regarding the relationship between environmental variables and the genetic structure of the population, when considering the entire period (2003, 2005 and 2008) the first two components of the PCA explained 60% of total ordination for the 3 years. The first component explained 33% of total ordination, and was associated with the variation of allele frequencies in ADA (R = 0.67), but any model among those performed by the GAM analysis explained the relationship between dependent (first axis of PCA) and independent variables (salinity and oxygen concentration).
Considering the data set collected during eutrophication (2003 and 2005), in this case the first two components of the PCA explained 65% of total ordination. The first component explained 47%, and was strongly associated with the variation of allele frequencies in ADA (R = 0.87). Among all models, the only one explaining the relationship between dependent (first axis of PCA) and independent variables (salinity and oxygen concentration) was the one that considered the dependence between maximum salinity and minimum oxygen concentration, while maximum salinity and minimum oxygen concentration alone did not yield significant results (Table 6). 3.206 0.032* df = degrees of freedom; t = Student's t critical value; * = P ≤ 0.05.

DISCUSSION
As hypothesised, the genetic structure of the local population of A. fasciatus showed significant genetic heterogeneity and HW disequilibria at the time of severe environmental conditions, during eutrophication. This supports the prominent action of diversifying evolutionary forces under stressful conditions, which were strong enough to overwhelm the homogenising power of gene flow, even at the microscale (about 2 km) and in the absence of landscape barriers. Among the many mechanisms promoting divergence at the microgeographic scale (Richardson et al. 2014), the analyses highlighted the roles of natural selection and genetic drift. The intensity of action of these evolutionary forces appeared to follow the environmental gradient of the saltwork, and to decrease when environmental stress relaxed.
As postulated, a genetic structuring among fish groups was recorded following the environmental gradient during the eutrophication period. The F ST and the exact probability test analysis showed a relatively weak but significant divergence during eutrophication, which concerns mostly three loci (ADA, GPI-3 in 2003; ADA and PGM-1 in 2005), besides the overall bulk of loci in 2005 (north vs south; see Table 4). In particular, the spatial and temporal patterns of the ADA locus were consistent with our hypothesis: during the eutrophication period (2003 and in 2005) ADA showed significant F ST values and ADA*100 allele frequencies presented a clinal distribution, increasing significantly from the northern to the southern hyperhaline and hypoxic ponds. The low, although significant, F ST values at this locus may depend on its allelic frequencies, as Jacobsson et al. (2013) recently demonstrated that the upper bound of F ST is rapidly restricted to small values when the frequency of the most common allele and homozygosity are higher or lower than 0.5, as in this case.
The analysis excluded an isolation by distance model (IBD) as the cause of the genetic structuring during eutrophication. Conversely, the current data analysis supports the hypothesis of a role for natural selection in generating the divergence and the ADA cline during stress, whereas the other polymorphic loci turned out to be neutral. It has already been demonstrated that the genetic structure of populations inhabiting coastal lagoons may be influenced by natural selection, promoted by highly variable environmental conditions (mainly salinity) (González-Wangüemert et al. 2009;González-Wangüemert & Vergara-Chen 2014;Rodrigues et al. 2015). Further, geographical clines within natural populations have been considered to provide evidence for selection (Gysels et al. 2004;Umina et al. 2005), although they usually emerged on a much wider spatial scale (Lenfant & Planes 2002;Cimmaruta et al. 2005). For example, selection at the ADA locus has already been invoked to explain the clinal distribution of its allelic frequencies in the marine dusky grouper fish (De Innocentiis et al. 2001).
Concerning the divergence observed at the other loci (GPI-3 and PGM-1), and over all loci in 2005, the signs of a recent bottleneck that emerged in the 2003 southern sample suggest that the previously recorded genetic erosion, associated with the reduction of the effective population size (N e ) , may have derived from bottlenecks localised in the southern area, which promoted faster and random shifts of the frequencies. Additional explanations for divergence and genetic erosion may include variation-reducing selection (Leffler et al. 2012): if selection acts on one locus as ADA, it may also enhance the loss of genetic variation and shifts of allelic frequencies at nearby or linked sites. Alternatively, selection at the ADA locus and drift at other loci, due to the reduction of N e , may be simultaneous and independent effects of environmental stress.
A number of processes could explain the deficits of ADA heterozygotes recorded during the stressful periods. Selection against heterozygotes may generate a similar pattern, as it leads to heterozygote deficits and to increasing frequency of the most common allele (here, ADA*100), but it is not supported by the fact that the deficits occurred mainly in the less-stressed northern area, while the frequencies of the ADA*100 allele increased southwards. Inbreeding too may result in heterozygote deficits, but the analysis of the relatedness coefficients, which were always close to 0 and significantly lower in northern samples, does not support this hypothesis either. A Wahlund effect seems to be the most likely hypothesis, since the northern samples showed a high proportion of individuals assigned to a discrete sub-population, besides the significant F IS value at this locus in 2003. Moreover, the allozyme loci showing heterozygote deficits are the same, having differentiated allele frequencies through the saltwork, which is a signature of the Wahlund effect (Allendorf & Luikart 2007). In other words, this pattern could arise from the active displacement and mixing of groups of fish within the study area, whose allelic frequencies at the ADA locus become heterogeneous during stress, because of selection. González-Wangüemert and Vergara-Chen (2014) proposed a similar hypothesis to explain the pattern of genetic divergence and heterozygote deficits recorded in a lagoon population of Potamoschistus marmoratus (Gobiidae), which would result from the combined action of selection, random reproductive success linked to the variability of environmental conditions, and displacement of individuals. The Wahlund effect could explain the other heterozygote deficits recorded in 2003 (GPI-3, PGM-1 and Multilocus), since also in this case the deficits, low relatedness coefficients and high proportion of individuals assigned to discrete sub-populations match, though additional data may strengthen this conclusion.
The relationships between the genetic structure and environmental variables emerge only if we consider the data set collected in the 2 years when eutrophication occurred. The PCA/GAM indicated that the divergence at the ADA locus was relevant to defining the genetic structure of the population, which was in turn linked to higher values of salinity and low oxygen concentration, giving further support to the hypothesis of selection for this locus. ADA plays a regulatory function in the production of reactive oxygen species that can be linked to cellular oxygen levels and oxidative stress (Clanton 2007;Kälvegren et al. 2010). A connection between different functionalities of ADA genotypes in hypoxia could thus be assumed, although selection may involve a linked locus.
Considering also the post-recovery data set (2008), no relationship was detected between genetic and environmental data, likely because this sample did not exhibit any significant genetic structuring, and the environmental data showed less steep gradients. Indeed, according to our expectations, no deviations from HWE were detected in 2008, allelic and genotypic frequencies became homogeneous, and ADA showed similar frequencies at the four sampling sites and was no longer a candidate for directional selection. The strength of the diverging forces appeared to relax after habitat restoration, allowing gene flow re-establishing genetic homogeneity and HWE over a few generations.

CONCLUSIONS
The pattern of divergence observed, although concerning few loci, is relevant, given the relative scale of its occurrence, and evidenced the action of diverging forces that have previously been considered to be rare in nature at the microscale level (Kingsolver et al. 2001;Richardson et al. 2014). Theoretically, divergence and HW disequilibria may have little chance to arise in a small lagoon and in a mobile species such as A. fasciatus, and genetic structuring usually needs a wider geographical scale to emerge in the face of gene flow, as shown for many fish species (Berg et al. 2015;Guo et al. 2015;Konijnendijk et al. 2015). As hypothesised in our case, even the divergence at a single locus might be important for the conservation of local populations, if it has an adaptive value to persisting extreme conditions. Further, our study showed the potential for evolutionary forces to interact, causing secondary effects. Overall, the intensity of the environmental pressure appeared to define the relative strength of forces, with those promoting genetic divergence emerging only when stress levels were raised. APPENDIX Table A1.