Environmental variables, habitat discontinuity and life history shaping the genetic structure of Pomatoschistus marmoratus

Abstract

Coastal lagoons are semi-isolated ecosystems exposed to wide fluctuations of environmental conditions and showing habitat fragmentation. These features may play an important role in separating species into different populations, even at small spatial scales. In this study, we evaluate the concordance between mitochondrial (previous published data) and nuclear data analyzing the genetic variability of Pomatoschistus marmoratus in five localities, inside and outside the Mar Menor coastal lagoon (SE Spain) using eight microsatellites. High genetic diversity and similar levels of allele richness were observed across all loci and localities, although significant genic and genotypic differentiation was found between populations inside and outside the lagoon. In contrast to the F _ST values obtained from previous mitochondrial DNA analyses (control region), the microsatellite data exhibited significant differentiation among samples inside the Mar Menor and between lagoonal and marine samples. This pattern was corroborated using Cavalli-Sforza genetic distances. The habitat fragmentation inside the coastal lagoon and among lagoon and marine localities could be acting as a barrier to gene flow and contributing to the observed genetic structure. Our results from generalized additive models point a significant link between extreme lagoonal environmental conditions (mainly maximum salinity) and P. marmoratus genetic composition. Thereby, these environmental features could be also acting on genetic structure of coastal lagoon populations of P. marmoratus favoring their genetic divergence. The mating strategy of P. marmoratus could be also influencing our results obtained from mitochondrial and nuclear DNA. Therefore, a special consideration must be done in the selection of the DNA markers depending on the reproductive strategy of the species.

Introduction

Studies of coastal lagoon fauna are extremely important from ecological and evolutionary perspectives. Coastal lagoons are semi-isolated ecosystems separated from the sea by both physical barriers and ecophysiological boundaries. These ecosystems suffer frequent environmental disturbances, being exposed to wide fluctuations on their physical–chemical parameters (Pérez-Ruzafa et al. 2005, 2007). Such environmental variation may subject the local fauna to severe adaptive challenges that could have a direct effect on the genetic composition of some species such as Aphanius fasciatus, Diplodus sargus, Elysia timida, Cerastoderma glaucum, Ostrea edulis, Ostreola stentina and Bursatella leachii (Cognetti and Maltagliati 2000; Maltagliati 2002; González-Wangüemert et al. 2004, 2006, 2009; Giménez-Casalduero et al. 2011; González-Wangüemert and Pérez-Ruzafa 2012) or separate species into different populations (Bilton et al. 2002; Trabelsi et al. 2004; González-Wangüemert et al. 2009; Richards et al. 2010; Fluker et al. 2011; Vergara-Chen et al. 2013; González-Wangüemert et al. 2014). The isolation and habitat fragmentation have been also supposed to contribute to the genetic structure of marine species with populations inhabiting coastal lagoons, enhancing the effects of microevolutionary processes (Watts and Johnson 2004).

Several physical and biological factors can reduce gene flow between populations and hence produce the appearance of spatial genetic structure in marine species. Patterns of gene flow in benthic shore fishes appear to stem in part from aspects of life history, larval dispersal, adult migration, habitat discontinuity and coastal water currents (Marques et al. 2006; Giovannotti et al. 2009; Earl et al. 2010; Riginos et al. 2011; Hirase et al. 2012; Durand et al. 2013). Moreover, differences in reproductive strategy may have an important impact on genetic diversity and population structure (González-Wangüemert and Pérez-Ruzafa 2012; Portnoy et al. 2013). In addition, there is evidence of genetic adaptation as result of natural selection capable of sustaining adaptive divergence on contemporary time scales (Conover et al. 2006; Nielsen et al. 2009; Yoboué et al. 2012; Wang et al. 2013). All these factors could have a profound effect on the arising of population differentiation.

Population differentiation in marine fishes through drift is hoped to be weak because of their populations have relatively shallow histories and are often very large (Hauser and Carvalho 2008). However, as it was indicated above, other factors as a restricted gene flow, selection events or some mating systems such as polygyny could favor the genetic structuring of fishes. In fact, evidence of genetic differentiation across large and small spatial scales has been reported in a great number of marine fishes using microsatellites (Pampoulie et al. 2004; Larsson et al. 2007; Bradbury et al. 2009; Earl et al. 2010; González-Wangüemert et al. 2010, 2012; Larmuseau et al. 2010a; Horne et al. 2011; González-Wangüemert and Pérez-Ruzafa 2012). These nuclear markers are especially useful to study the fine geographical variation of marine populations because they tend to be highly variable and can discern even subtle genetic differences (Waples 1998; Hedgecock et al. 2007), resolving population structures that are not detected by mitochondrial DNA and allozyme markers (De Innocentiis et al. 2001; Pampoulie et al. 2004; Bisol et al. 2007; Canino et al. 2010). In this sense, several studies based on microsatellite markers have suggested that estuaries and coastal lagoons may offer special opportunities (because of their environmental features) to find local genetic structuring of fish populations (Beheregaray and Sunnucks, 2001; Bisol et al. 2007; Roberts et al. 2010; McCraney et al. 2010; González-Wangüemert and Pérez-Ruzafa 2012).

The marbled goby Pomatoschistus marmoratus (Risso 1810; Teleostei: Gobiidae) is a small benthic fish, with lagoonal, estuarine and marine populations, inhabiting sandy, inshore, shallow waters from the eastern Atlantic, Mediterranean, Black and Azov seas and Suez Canal (Miller 1986) usually in high densities (Verdiell Cubedo et al. 2008). During the reproductive season, marbled goby males build a nest by cleaning the inside of empty bivalve shells and covering the outside with sand (Mazzoldi and Rasotto 2001). They defend the nest and take care of the eggs deposited on the upper valve by one or more females, cleaning and fanning the eggs until larvae hatch (Mazzoldi and Rasotto 2001). There is no information about the pelagic larval duration (PLD) of P. marmoratus, although congeneric species, such as Pomatochistus minutus (Pallas 1770) and Pomatochistus lozanoi (de Buen 1923), have a PLD of 30–39 days under laboratory conditions (Fonds 1970). Several species of the genus Pomatoschistus show little or no migration behavior (Gysels et al. 2004; Berrebi et al. 2005) and have a limited swimming ability because the pelvic fins have been fused into a suction disk (Miller 1986; Bardin and Pont 2002).

Therefore, considering all these biological and ecological features, P. marmoratus could be a good target species to study small-scale genetic structure and adaptation under highly variable environmental conditions in coastal lagoons. However, in the last years, only works considering great spatial scales have been published using different species from Pomatochistus genus and no study considering simultaneously mitochondrial and nuclear markers. Mejri et al. (2009) using mtDNA markers proposed that the complex genetic structure of Pomatochistus tortonesei was shaped by recurring shifts in sea level and sea surface temperatures of Mediterranean Sea which caused the desiccation of shallower lagoons and therefore colonization and re-colonization events of the brackish populations. Moreover, Berrebi et al. (2009) studied the genetic structure of two sedentary gobiid fishes (Pomatoschistus microps and P. marmoratus) along several French Mediterranean coastal lagoons using mitochondrial DNA RFLP markers; they found that populations of P. microps inhabiting neighboring lagoons showed a high level of isolation, which was not shown by P. marmoratus. In general, previous data published using mtDNA markers and medium/great spatial scales pointed that the Pomatoschistus populations inhabiting different brackish habitats could be characterized by unique mitochondrial haplotypes which are well defined in relation to the limited gene flow and restricted dispersal abilities of the different species (Mejri et al. 2011, 2012). Therefore, as concluded by Mejri et al. (2011), the role of the Mediterranean lagoon habitat should be related to how much it amplifies the effects of historical (e.g., past sea level changes) and environmental (e.g., present-day hydrographic regime) processes in regard to the genetic structure of the inhabiting Pomatochistus species.

Recently, Vergara-Chen et al. (2010a) carried out a study about the genetic structure of P. marmoratus considering small spatial scales (<25 km) into Mar Menor coastal lagoon (SE Spain) and between adjacent localities in Mediterranean Sea (<100 km) using mitochondrial DNA sequences (control region). This study found very low population differentiation considering F _ST values and high gene flow rates between lagoon and nearby open sea locations. However, a higher genetic diversity and occurrence of exclusive haplotypes into the different lagoon localities were observed, suggesting a possible adaptive potential in this species.

Considering these interesting data, we analyze the same five localities using eight polymorphic microsatellite loci, which could show us the genetic gradient of P. marmoratus at small spatial scale with higher accuracy. We also look for factors explaining the observed genetic patterns, including environmental features and life history. Finally, we assess the concordance of genetic patterns obtained from our previous mitochondrial results (control region; Vergara-Chen et al. 2010a) and current nuclear data (eight polymorphic microsatellite loci).

Materials and methods

Study sites

The Mar Menor is a hypersaline coastal lagoon located in a semiarid region in the southeastern Spain (Fig. 1). It is one of the largest coastal lagoons in the Mediterranean basin, with a surface area of 135 km² and an average depth of 3–4 m. It is separated from the Mediterranean Sea by a 22-km-long sand bar called La Manga with three narrow inlets connecting the lagoon to the open sea: Las Encañizadas, El Estacio and Marchamalo. The lagoonal salinity ranges from 38 to 51 psu, and temperature varies from 10 °C in winter to 32 °C in summer (Pérez-Ruzafa et al. 2005; González-Wangüemert et al. 2009). Its bottom is principally covered by dense meadows of the invasive macroalga Caulerpa prolifera, although shallowest areas are covered by scarce patches of the phanerogams Cymodocea nodosa and Ruppia cirrhosa (García-Sánchez et al. 2012; Martinez-Garrido et al. 2014). The lagoon maintains a diverse fish community due to its environmental heterogeneity: sandy without vegetation and muddy bottoms, rocky bottoms and seagrass beds composed of monospecific and multispecific meadows of C. nodosa, R. cirrhosa and C. prolifera (Pérez-Ruzafa et al. 2006). The Mar Menor coastal lagoon hydrodynamic is driven by the wind (mainly northeast), generating a circulation pattern comprising three main gyres and dominant currents from north to south along the internal coast of the lagoon (González-Wangüemert et al. 2009; de Pascalis et al. 2012). According to its hydrographical characteristics, three main basins have been identified in the lagoon (Fig. 1) (Pérez-Ruzafa et al. 2007; González-Wangüemert et al. 2009). The northern basin has the greatest influence of Mediterranean waters through the entrance channels at the north of La Manga and El Estacio with the lowest mean salinity values. The southern basin is a more confined area with the most saline water, and the central basin has intermediate salinity values due to the mixing of Mediterranean and lagoonal waters.

Fig. 1

Study area: the Mar Menor lagoon-SW Mediterranean transition (SE Spain) showing sampling locations for Pomatoschistus marmoratus. The Mar Menor map is showing the interpolated surface of mean salinity and the current pattern representing by the arrows (González-Wangüemert et al. 2009). White circles indicated lagoonal sites: Lo Pagán (LP), Los Urrutias (LU), Playa Honda (PH). Black triangles indicated marine locations: Veneziola (VZ), Mazarrón (MZ)

Otherwise, the southwestern Mediterranean Sea is characterized by lower extreme values of temperature and salinity than Mar Menor lagoon, oscillating its salinity between 36.84 and 37.41 psu and its temperature between 13.83 and 26.44 °C (available information at: www.puertos.es/oceanografia_y_meteorologia/).

Thus, the environmental heterogeneity observed along this ecological gradient may be considered as an outdoor laboratory to assess the relationships between the environmental conditions and the spatial distribution of the genetic variation in coastal transitional waters. Three sampling sites were selected inside the Mar Menor coastal lagoon: Lo Pagán, located on the northern coast of the Mar Menor; Los Urrutias, located on the west side; and Playa Honda, located on the south. Two sampling locations were selected in the Mediterranean Sea: Veneziola (seaward side of La Manga) and Mazarrón (southern Murcia) (Fig. 1). We highlight the lack of appropriate habitats and absence of P. marmoratus populations along the Murcia coast from Mar Menor to Mazarrón.

Field sampling procedures

Fish samples were collected using a beach seine over shallow sandy bottoms. The sample size consisted of 35–40 individuals per collecting site. The marbled gobies were identified on the basis of external morphological characters, killed by freezing and preserved in 100 % ethanol until tissue dissection. Tissue samples of muscle were removed from each specimen and preserved in 100 % ethanol.

DNA extraction and microsatellite genotyping

The tissue was dissolved using lysis buffer solution with proteinase K, and total genomic DNA was isolated by protein precipitation and final cleanliness with ethanol (Sambrook and Russell 2001). PCRs were carried out for eight microsatellite loci (Pmar03, Pmar05, Pmar08, Pmin04, Pmin09, Pmin29, Pmin35 and Pmin38) (Berrebi et al. 2006; Larmuseau et al. 2007) (Table 1). The loci Pmin04, Pmin09, Pmin29, Pmin35 and Pmin38 designed on P. minutus and tested only by cross-species amplification for P. microps, P. lozanoi and P. pictus (Malm 1865) (Larmuseau et al. 2007) amplified well for P. marmoratus. The PCR was a volume of 10 μl containing 300 ng of genomic DNA, 10× buffer, 1.75 mM MgCl2, 0.4 mM of each dNTPs, 3.75 μm of each primer [using forward primers labeled with FAM (SIGMA), HEX (SIGMA) and NED (Applied Biosystems)] and 0.3 U of Taq polymerase (Ecogen). Amplification was performed in a 2,720 thermal cycler (Applied Biosystems) programmed for an initial denaturation at 95 °C for 2 min, followed by 35 cycles composed of denaturation at 95 °C for 1 min, annealing at an appropriate temperature (Table 1) for 1 min, extension at 72 °C for 1 min and a final extension at 72 °C for 3 min. Allele sizes were scored according to the protocols from Molecular Biology Service (Center of Marine Science, Faro, Portugal) using an ABI Prism 3130 automated genetic analyzer (Applied Biosystems).

Table 1 Characterization of the eight microsatellite loci used to genotype five populations of Pomatoschistus marmoratus from the Mar Menor coastal lagoon and adjacent marine sites (SE Spain)

Microsatellite data analysis

The genotype data were scored using the STRand software (v. 2.4.59 http://www.vgl.ucdavis.edu/STRand). The mean numbers of alleles per locus, allele frequencies, observed (H _O) and unbiased expected heterozygosity (H _E) were calculated in GENETIX version 4.05 (Belkhir et al. 2004). Deviations from Hardy–Weinberg equilibrium (HWE) were characterized by F _IS and tested using exact test in the software GENEPOP version 4.0.10. Also, we registered the F _IS values per locality and locus, estimating the exact P values by the Markov chain method (GENETIX version 4.05). In instances where the observed genotype frequencies deviated significantly from HWE, the program Micro-Checker version 2.2.3 (van Oosterhout et al. 2004) was used to infer one of the most probable causes of such HWE departures. Some authors have described that the presence of null alleles led to the overestimation of both F _ST and genetic distance in cases of significant population differentiation (Chapuis and Estoup 2007).

We used the software GENEPOP version 4.0.10 to estimate the genotypic (G-based) and genic differentiation for all pairs of populations (Raymond and Rousset 1995; Rousset 2008). Genetic linkage disequilibrium between locus pairs was estimated according to Weir and Cockerham (1979) and tested on contingency tables under the null hypothesis of independence (P < 0.05). To determine whether a population exhibits a significant number of loci with heterozygosity excess, we used the Sign and Wilcoxon tests implemented in the program BOTTLENECK version 1.2.2 (Cornuet and Luikart 1996). Computations were based only on the infinite allele model (IAM) because the two-phase mutation (TPM) model still assumes that loci mutate within repeat number, but some our loci were off repeat.

The Wright’s single locus F-statistics (Wright 1969) were calculated from allele frequencies for all loci examined for each population according to Weir and Cockerham (1984) in ARLEQUIN version 3.11 (Excoffier et al. 2005). Significance of F _ST for all loci and pairwise population comparisons was assessed by permutation of the values 10,000 times. Standard deviations of single-locus F _ST values were obtained by jackknifing over all populations. The Bonferroni correction for multiple comparisons (Rice 1989) was applied to all P values from F _ST estimates to compensate for possible type I errors resulting from multiple pairwise comparisons. Genetic distances (D _CE) (Cavalli-Sforza and Edwards 1967) were computed between pairwise samples in GENETIX. We used Cavalli-Sforza distance because it is less affected by null alleles (Chapuis and Estoup 2007). Probabilities of random departure from zero for distance values, according to the null hypothesis, were read directly from the distribution of 10,000 randomized matrices computed by permutation. The Bonferroni correction was applied to all P values.

Populations were spatially clustered using correspondence analysis (CA) implemented with BiodiversityR package in R software (R Development Core Team 2007), which utilizes the genotype frequencies of populations as variables in order to visualize similarities among locations without assuming tree-like relationships. We used as genotype data, the allele frequencies obtained from the microsatellite data and the haplotype frequencies from our previous study with mitochondrial DNA control region sequences (Vergara-Chen et al. 2010a). An analysis of molecular variance (AMOVA) (Excoffier et al. 1992) was carried out in ARLEQUIN to assess the hierarchical partitioning of genetic variability within and among populations and among lagoon and marine groups.

Population structure was inferred using Structure version 2.3 software by the method of Pritchard et al. (2000) from multilocus genotype data. Each K was replicated 20 times for 100,000 iterations after a burn-in period of 5,000, without any prior information on the population of origin of each sampled individual. The height of the modal value of distribution for the posterior probability of the data for a given K was used as an indicator of the strength of the signal detected by Structure (Evanno et al. 2005).

To investigate gene flow linking to the lagoonal localities and between lagoonal and marine sites, migration rates based on maximum likelihood were obtained with the program MIGRATE v. 3.2.7 (Beerli and Felsenstein 2001; http://popgen.scs.fsu.edu). MIGRATE uses a Markov Chain Monte Carlo-based (MCMC) approach to explore all possible gene genealogies to provide maximum likelihood estimates of the population size and migration rates compatible with the data. These estimates are computed as θ = N _e μ as population size and M = m/μ as the mutation-scaled migration rate for migration (nomenclature according MIGRATE-n software), where N _e is the effective population size, m is the fraction of the new immigrants in the population per generation, and μ is the mutation rate of the gene. The MCMC run consisted of 10 short chains (sampling 10,000 trees) and one long chain (sampling 10,000 trees) with a burn-in period of 10,000 trees.

To assess the presence of microsatellites loci under selection, we used LOSITAN version 1.0.0 (http://popgen.eu/soft/lositan/) (Antao et al. 2008). This program evaluates the relationship between F _ST and expected heterozygosity (H _E) to identify outlier loci. We ran 75,000 simulations with “neutral mean F _ST” and “force mean F _ST,” to increase the reliability of the mean F _ST and the entire microsatellite dataset under the infinite allele model (IAM). The number of populations was four according to F _ST values and Cavalli-Sforza genetic distances. We chose the confidence intervals of 99 % to carry out a more conservative test for selection.

To test the possible relationships between environmental variables and genetic diversity, we developed generalized additive models (GAMs) for the first principal component of the CA (which usually explains the higher values of variance). The coordinates of each locality derived of the CA were included in GAMs as dependent variables. GAMs are a nonparametric extension of generalized linear models (GLM) that fit a wide variety of forms of stochastic variation in the response. GAMs represent the relationship between the response variable and the predictors by smooth functions, which can take virtually any form (Hastie and Tibshirani 1990). These models have been applied on genetic and environmental data previously (Snäll et al. 2004; Parisod and Bonvin 2008; González-Wangüemert et al. 2009; Vergara-Chen et al. 2010b). The GAMs were evaluated by examining the proportion of explained deviance and minimizing the generalized cross-validation (GCV; Wood 2000; Wood and Augustin 2002) and the Akaike information criterion (AIC; Venables and Ripley 2004) scores. As independent variables, we used maximum, minimum and mean values of temperature and salinity (data from: http://www.noaa.gov; Research Group “Ecología y Ordenación de Ecosistemas Marinos Costeros”). GAMs were performed using “ade4” (Chessel 1992) and “mgcv” (Wood 2006) packages from R statistical software (R Development Core Team 2007).

Results

Genetic diversity

A total of 177 fishes were analyzed from three lagoonal and two marine locations. Mean expected heterozygosities (H _E) were high, ranging from 0.669 (Los Urrutias) to 0.748 (Lo Pagán). A similar level of allelic richness was observed across all samples with a mean of 9.025 alleles across all loci (ranging from 8.375 to 10.0 alleles). Regarding to the fixation index (F _IS), the Mazarrón sample showed the highest F _IS value (0.161) while the lowest value (0.065) was detected in the Veneziola sample (Table 2). Significant deviations from Hardy–Weinberg equilibrium (heterozygote deficit) were detected at all five populations. Four loci, Pmar03 and Pmin38 (3 localities), Pmin04 (4 localities) and Pmin29 (5 localities), deviated significantly from Hardy–Weinberg expectations (“Appendix 1”). Analysis of all samples indicated that some of the heterozygote deficits could be due to the presence of null alleles, but they were not attributable to genotyping errors. The recalculated frequency data using MICROCHECKER software to consider null alleles were not significant different from the detected first values, and therefore, we continued using them to all posterior analyses. Also, inclusion or exclusion of these loci did not quantitatively affect either population comparisons (F _ST values and genetic distances, data not shown).

Table 2 Estimates of genetic diversity of the five samples of Pomatoschistus marmoratus from the Mar Menor coastal lagoon and adjacent marine sites (SE Spain) based on eight microsatellite markers

Significant linkage disequilibrium was detected in four loci pairs from some localities (Pmar03/Pmin29 in Lo Pagán, Pmar03/Pmin38 and Pmar05/Pmin38 in Lo Pagán and Mazarrón, Pmar08/Pmin35 in Playa Honda). The IAM model was applied to find possible bottlenecks; however, they were not detected in any populations.

Population genetic differentiation

Pattern of population differentiation based on F _ST statistics pointed to significant differences between six sample comparisons, while we detected lack of differentiation between four sample comparisons (Table 3). The highest F _ST value (0.0303) was observed between Lo Pagán (lagoon) and Mazarrón (marine) samples, whereas the lowest significant value (0.0137) was found between Los Urrutias (lagoon) and Veneziola (marine) samples. Pairwise F _ST estimates among lagoonal samples indicated that Los Urrutias and Playa Honda were genetically different (F _ST = 0.0157) while Lo Pagán–Playa Honda and Lo Pagán–Los Urrutias showed absence of genetic differentiation. The marine locality of Veneziola showed genetic similarity with two lagoonal localities, Lo Pagán and Playa Honda. All lagoonal localities and the Veneziola marine sample were genetically differentiated from the Mazarrón marine sample. These results were corroborated by Cavalli-Sforza genetic distances (Table 3).

Table 3 Pairwise fixation indices (F _ST, below) and Cavalli-Sforza genetic distances (above) between five sampling localities of Pomatoschistus marmoratus from the Mar Menor coastal lagoon and adjacent marine sites (SE Spain) based on eight microsatellite markers

Furthermore, genotypic differentiation (G-based) for all pairs of populations showed significant differences between all comparisons across all loci except for some lagoon localities, Lo Pagán–Playa Honda (P = 0.3198) and Lo Pagán–Los Urrutias (P = 0.0863). The genic differentiation test showed the same pattern without significant differences between Lo Pagán and Playa Honda (P = 0.1087) (“Appendix 2”).

This structure was confirmed using a correspondence analysis (CA) based on allele frequencies. The surface plot of the two first ordination axes scores 64.58 % of the total variance in the data revealing a small-scale geographic structure (Fig. 2) and showing two groups (Fig. 2a): (1) Mazarrón sample on the positive side of Axis I; (2) a group composed of all remaining localities on negative side of Axis I. The Axis II discriminated the Veneziola sample from the lagoonal localities. This pattern indicating the differentiation of Mazarrón sample was reconfirmed with the CA plot using our previous data (frequencies of haplotypes) from mitochondrial DNA control region (Vergara-Chen et al. 2010a) pointed a score of 67.47 % of the total variance for the two first ordination axes (Fig. 2b).

Fig. 2

Correspondence analysis of genotype frequencies of five populations of Pomatoschistus marmoratus along the Mar Menor lagoon-SW Mediterranean transition (SE Spain). a Plot based on allele frequencies from eight microsatellite loci. b Plot based on haplotype frequencies from mitochondrial DNA control region sequences obtained by Vergara-Chen et al. (2010a). Populations are identified as follows: Lo Pagán (LP), Los Urrutias (LU), Playa Honda (PH), Veneziola (VZ), Mazarrón (MZ)

The analysis of molecular variance (AMOVA) pointed to nonsignificant differences among groups: coastal lagoon (Lo Pagán, Los Urrutias and Playa Honda) and marine (Veneziola and Mazarrón). A low percentage of the variance was attributed to differences among groups (0.12 %; P > 0.05). However, the analysis revealed significant differences among populations within groups and within populations (1.23 and 98.65 %, respectively; P < 0.05).

The P. marmoratus populations did not show a definite structure such as was obtained by the Bayesian approach. This analysis suggests that K = 1 best describes the situation among samples (Fig. 3).

Fig. 3

Estimated population structure of Pomatochistus marmoratus and surface population using STRUCTURE (Pritchard et al. 2000) for K = 1. Each individual is represented by a thin vertical line, which is partitioned into K segments that represent its estimated population group membership fractions

The MIGRATE analysis pointed the lowest value of effective population size (Θ) in Los Urrutias (2.04) and the highest value in Lo Pagán (3.06). Migration rates showed a predominant gene flow from Playa Honda to Lo Pagán (3.78) and the lowest value of genetic exchange from Lo Pagán to Veneziola (0.87) (Table 4).

Table 4 Effective population size (Θ = N _e μ) and migration rates (m/μ = Nem/Θ) based on eight microsatellite loci of Pomatoschistus marmoratus from the Mar Menor coastal lagoon and adjacent marine sites (SE Spain)

LOSITAN tests did not identify any microsatellite as an outlier (Fig. 4), so no loci can be consider under positive or balancing selection.

Fig. 4

Comparison of F _ST and H _e in polymorphic loci to identify outliers and potential candidates for selection using LOSITAN software. Graphical output shows the simulated confidence area for neutral loci (pale gray shading)

Relationships among genetic structure and environmental variables

Our analyses, except for AMOVA and assignation tests, have detected a small but significant genetic differentiation between coastal lagoon and marine populations of P. marmoratus, which could be influenced (among others factors) by environmental variables.

The application of GAMs to our previous published data from control region in P. marmoratus (Vergara-Chen et al. 2010a) supports this relationship among genetic structure and environmental variables, obtaining a linear response to maximum salinity. The deviance explained by the GAM was high (Dev_{maximum salinity} = 98.2 %, Pr (intercept) = 0.0014, Pr_{maximum salinity} = 0.0010). The GCV and AIC scores were significant (GVC_{maximum salinity} = 0.0329, AIC_{maximum salinity} = −1.9807). Similar results are obtained when we apply the GAMS to our data from microsatellites (Dev_{maximum salinity} = 98.5 %, Pr (intercept) = 0.0017, Pr_{maximum salinity} = 0.0016; GVC_{maximum salinity} = 0.0685; AIC_{maximum salinity} = 1.6816).

Comparison between mitochondrial and microsatellite data

In contrast to the F _ST values from mitochondrial DNA (Vergara-Chen et al. 2010a), the microsatellite data exhibited significant evidence for population genetic structure between coastal lagoon and marine samples. This was also reflected in overall F _ST values that were higher for nuclear microsatellites (global F _ST = 0.014) than for mitochondrial DNA (global F _ST = 0.006). These differences were more important between Mazarrón marine sample and all remaining localities.

This pattern was confirmed when we compared the correspondence analyses from allele and haplotype frequencies (microsatellites and mtDNA, respectively); the results pointed the same grouping with Mazarrón isolated from the remaining sampling sites. Otherwise, both types of molecular markers indicated variable rates of genetic exchange among lagoonal samples and between coastal lagoon and marine localities although the overall of the pattern was very similar. The parameter Θ (a function of H _e) based on mitochondrial DNA data suggested values between 0.01 (Playa Honda) and 0.05 (Los Urrutias), while the migration rates showed significant gene flow from Playa Honda to all remaining localities, and restricted flow from Mazarrón to Veneziola (3.52e−11) and from Lo Pagán to Playa Honda (9.36e−14) (Vergara-Chen et al. 2010a). Microsatellite data pointed values of effective population size higher than those obtained with the mitochondrial marker, with the lowest (Θ) value for Los Urrutias and the highest value for Lo Pagán. Migration rates showed a predominant gene flow from Playa Honda to the Lo Pagán (Table 4). The AMOVA for both mitochondrial DNA and microsatellite loci revealed no significant evidence for genetic subdivision among groups (coastal lagoon vs. marine).

Discussion

Microsatellite loci showed high levels of polymorphism in P. marmoratus which is consistent with the mean number of alleles observed in other Pomatoschistus species (Berrebi et al. 2006; Larmuseau et al. 2007; Marques et al. 2012). The studied populations showed high genetic diversity, with values of expected heterozygosity ranged from 0.67 to 0.79. These levels of genetic variability were also similar to those found in a congeneric species, P. minutus, analyzed with microsatellites markers: 0.74–0.86 (Pampoulie et al. 2004), 0.74–0.78 (Larmuseau et al. 2010a) and 0.58–0.63 (Boissin et al. 2011).

Departures from Hardy–Weinberg equilibrium (HWE) were found in all Pomatochistus samples, such as were found in P. minutus and P. microps populations by other authors (Pampoulie et al. 2004; Larmuseau et al. 2010a; Marques et al. 2012). The most probable causes of HWE departures may be attributed to Wahlund effect due to the recruitment of genetically variable cohorts of larvae (González-Wangüemert et al. 2007), inbreeding, null alleles (Horne et al. 2011; González-Wangüemert and Pérez-Ruzafa 2012; Milana et al. 2012), cryptic species, groupings of relatives, scoring errors or selection against heterozygotes (Pampoulie et al. 2004). The Wahlund effect, the most common explanation of heterozygote deficiency, should result in significant F _IS values at more than one locus, as drift causing population structuring should affect all polymorphic loci similarly (Pogson et al. 1995), but we did not find this pattern in our results (“Appendix 1”). We think that the high F _IS values found in the studied localities are mainly due to structuring and mixing of sand goby populations. Inbreeding seems an unlikely explanation in fish with large populations such as gobies that are not subject to drastic reduction in their effective population size by fishery. However, we cannot exclude the possibility of inbreeding derived by the polygyny in P. marmoratus species. We do not favor the hypothesis of null alleles because all F _IS estimates were positive and significant; it seems highly improbable that all loci exhibit null alleles with such a constant frequency. The heterozygotes deficiency detected in all populations can be boosted by the use of locus primers designed for other species (cross-amplification) (Rungis et al. 2004; Pashley et al. 2006; Keever et al. 2008). Finally, the selection against heterozygotes cannot be demonstrated from our results, although we have detected a significant link among genetic structure and maximum salinity (implying that selection could be acting on P. marmoratus). LOSITAN did not detect any locus under balancing or positive selection. Further research will be necessary to test this hypothesis.

In spite of the lack of a geographic structure at total scale showed by the STRUCTURE results, signs of a genetic substructure were found. The correspondence analyses, Cavalli-Sforza genetic distances and the F _ST values showed significant genetic differentiation between some localities after Bonferroni correction; however, AMOVA did not provide evidence of differences between sample groups (coastal lagoon vs. marine), being Veneziola included in the marine group. Our genetic data seem to demonstrate that this locality could be considered more lagoonal than marine; also, the Veneziola’ habitat is more similar to lagoonal localities than Mazarrón, but it is not correct to carry out an AMOVA considering four samples inside “lagoon group” and one sample as “marine group.” On the other hand, the lack of significance in the AMOVA comparisons between groups could be also due to low power associate with a low number of populations per group such as was described by Fitzpatrick (2009). Therefore, according to our results from F _ST values, Cavalli-Sforza genetic distances, correspondence analysis and significant genic and genotypic differentiation, we suggested that the highest genetic divergence observed between samples (Mazarrón and the rest of localities) could be consequence of the isolation of lagoonal populations from open sea which has been already reported in other coastal fishes such as Pomatoschistus minutus (Pampoulie et al. 2004; Berrebi et al. 2009), Dicentrarchus labrax (Lemaire et al. 2000), Atherina lagunae (Trabelsi et al. 2004), Atherina boyeri (Kraitsek et al. 2008; Milana et al. 2012), Diplodus sargus (González-Wangüemert and Pérez-Ruzafa 2012) and several invertebrate species (Camilli et al. 2001; González-Wangüemert et al. 2006; Tarnowska et al. 2010; Vergara-Chen et al. 2010b).

The genetic differentiation among our P. marmoratus samples from lagoonal and marine sites could be also influenced by environmental discontinuities along the coastal lagoon–marine transition and even inside coastal lagoon. The habitat discontinuities could cause the spatial isolation of the populations mainly on benthic marine species (Riginos and Nachman 2001; Johansson et al. 2008; Riginos et al. 2011). Inside Mar Menor lagoon, there are important habitat discontinuities related to spatial distribution of substrate type and submerged marine vegetation (González-Wangüemert et al. 2009; Quintino et al. 2010) which can be influencing the distribution and density of P. marmoratus populations. In fact, some studies about the distribution of fish communities in Mar Menor showed significant differences among densities and standing stocks of P. marmoratus depending on bottom granulometry and seagrass and algae cover (Verdiell Cubedo et al. 2008). This discontinuity of available habitats for P. marmoratus is also found in the Mediterranean coast from Veneziola to Mazarrón, explaining the genetic isolation and the important genetic break detected between coastal lagoon localities (including Veneziola) and Mazarrón.

On the other hand, the reproductive strategy of P. marmoratus plus a highly variable environmental conditions create the potential for sweepstakes reproductive success (Hedgecock 1994; Hedgecock et al. 2007) which could be also explaining the genetic differentiation detected on our target species such as was registered on other fish species previously (Planes and Lenfant 2002; González-Wangüemert et al. 2007; Christie et al. 2010; Hogan et al. 2010).

The presence of exclusive alleles and differences on the allele frequencies existing between localities could also suggest that some environmental variables are influencing the genetic differentiation among populations; this hypothesis is supported by GAM results. As P. marmoratus is associated with brackish and hyper-saline habitats, which are systems suffering frequent environmental fluctuations and disturbances, this species could have been exposed to extreme conditions more frequently than other marine species inhabiting open sea habitats. Our results point a significant link between extreme lagoonal environmental conditions (mainly maximum salinity) and genetic composition, such as it was already demonstrated on cockles and sea cucumbers populations from the Mar Menor lagoon (González-Wangüemert et al. 2009; Vergara-Chen et al. 2010b). This find is very interesting because previously several authors had defined different ecological requirements for the various Pomatochistus species; for example, P. marmoratus survives to higher salinities than P. microps, especially in habitats whose temperature and salinity change through the year, but do not vary abruptly (Berrebi et al. 2005). Thereby, extreme environmental conditions could be acting on genetic structure of coastal lagoon populations of P. marmoratus causing genetic divergence (Conover et al. 2006; Marshall et al. 2010; Sanford and Morgan 2011).

The genetic changes detected on P. marmoratus are not explicitly linked to functionality, but really the signal recovered might be a signature of adaptive selection. This intraspecific pattern might be the base underlying evolution at the interspecific level. Some studies focused on several Pomatochistus species using a well-defined gene under selection (rhodopsin gene) have demonstrated that even the local light availability acts as selection pressure, sand gobies seem to be genetically adapted on the rhodopsin gene (RH1) to the differences in light between seas, lagoons and rivers (Larmuseau et al. 2009, 2010b).

Finally, the detected differences using two molecular markers (mitochondrial and nuclear) can be explained by their intrinsic characteristics and mutation rates. Nevertheless, the reproductive strategy of P. marmoratus could be also influencing our results. Mitochondrial DNA is a good marker for studying matrilineal movements because of its haploid nature, maternal inheritance and lack of recombination (Avise 2000). Therefore, the mitochondrial markers are influenced by mode of reproduction of species (Consuegra and de Leaniz 2007; Cano et al. 2008). The absence of significant genetic structure and high gene flow detected in our previous genetic study using control region (mitochondrial DNA) (Vergara-Chen et al. 2010a) could be considered indicative of extensive female dispersal related with P. marmoratus breeding behavior. This reproductive behavior shows a female spawning in different nests of different males; each male defends its nest and takes care of the eggs until larvae hatch (Mazzoldi and Rasotto 2001). Hence, the females are expected to be the most mobile sex during the breeding season (Lindström et al. 2006) favoring genetic population homogeneity detected specially using mitochondrial DNA markers. However, nuclear DNA markers (microsatellites) are diploids, with biparental inheritance. This leads to the expectation that more genetic variability and population differentiation would be detected with nuclear markers (Ruzzante et al. 1998). In fact, our microsatellite results supported this expectation, showing higher values of effective population size and lower rates of gene flow in contrast to our previous mitochondrial data (Vergara-Chen et al. 2010a). Therefore, these results corroborate the conclusions of Daemen et al. (2001) affirming that cytoplasmic markers provide a record of the historic and recent gene flow in females, while the nuclear markers estimate recent gene flow in both sexes.

Appendices

Appendix 1

See Table 5.

Table 5 F _IS values per locality and locus

Appendix 2

See Table 6.

Table 6 Total genic and genotypic differentiation between populations