The French part of the Mediterranean coast constitutes the Gulf of Lyon, where both species were captured in 12 lagoons in 2004. This January-May sampling campaign was the shortest possible. Beach seine nets were used by scientists and fyke nets by professional fishermen. The bigger Thau lagoon (19×8 km) was visited at six localities in March (Table 1 and Fig. 1) and the salinity was measured in each place (Table 1). After laboratory observation and sometimes rapid microsatellite genotyping at one or several strictly diagnostic loci (Pmar-03, Pmar-08, Pmin-08, Pmin-10, Pmin-11, according to Berrebi et al. 2006), the species P. minutus was removed from the 17 collected samples. Finally, 512 gobies were clearly determined as belonging to the two sedentary species P. microps and P. marmoratus and constituted the sample for this study. The fish were individually preserved in 90% ethanol until dissection and DNA extraction. Table 1 gives detailed information about the samples and Figure 1 indicates the position of each sampling site.

Eleven lagoons and one small river mouth from the Spanish border to Marseille were sampled (Fig. 1). These sites are all more or less brackish but display high ecological diversity. A tentative classification of the 12 water bodies, highlighting ecological categories according to the area and salinity is given in Table 2. They can be divided into “big lagoons” exceeding 5000 ha in area (Berre, Thau, Vaccarès and Salses-Leucate) and with a mean depth of 2 to 6 m, and “small lagoons” of 45 to 3700 ha and a mean depth 0.3 to 1.3 m. Apart from size and water volume, salinity is a very important parameter for lagoon inhabitants, generally fluctuating according to marine or continental influences that depend on marine or continental water input (rivers, rains), essentially driven by the seasons. The water exchanges and salinities given in Table 2 are only indicative because they always fluctuate. Haloclines have been described (Poizat et al. 2004) but are not considered in this simplified description. Because only P. marmoratus cannot bear null salinity (Rigal et al. 2008), the minimum salinity is another important ecological parameter. Therefore, Canet, Mauguio and Vaccarès are three clearly unfavourable lagoons for P. marmoratus. In the Vaccarès-Impériaux lagoon complex, constituting a large part of the Camargue region or Rhône Delta, the Impériaux lagoon in the south, where P. marmoratus can be found (Berrebi et al. 2005), is saltier and more euryhaline than the Vaccarès lagoon in the north.

DNA was extracted from fin tissue samples by the Chelex/proteinase K protocol of Estoup et al. (1996). Among the 14 loci tested, six microsatellite markers allowing cross-priming for three Pomatoschistus species were retained. Details of each locus are given in Table 3. The 5’ end of one of the two primers was covalently linked to fluorescein, Cy3 or Cy5 labels. Polymerase chain reactions (PCR) were performed in an Eppendorf Mastercycler programmable thermocycler with a 10 µL reaction total volume containing 0.1 U of Taq polymerase (Sigma-Aldrich), 2.5 mM MgCl₂, 0.4 mM of each dNTP (Invitrogen), 1 µL of 10X reaction buffer and 0.75 µM of each primer (Eurofins MWG). The thermal cycle was composed of an initial denaturation (94°C, 5 min); then of 34 repetitions of denaturation (94°C, 45 s), annealing (45 s at the temperatures given in Table 3 for each locus) and extension (72°C, 45 s); and then a final extension (72°C, 5 min).

Amplified DNA fragments were separated according to their size in an 8% polyacrylamide denaturing gel (Bio-Rad). The PCR products were visualized on a Hitachi FMBIO-II fluorescent imaging system scanner. Allele size was determined by comparison with a fluorescently labelled ladder of known size (100–600 bp, Promega), with the help of an image analysis software FMBIO ANALYSIS 8.0 (Hitachi).

The genotype matrix was then constructed and used as the basis for all of the following statistical analyses.

Assignment tests were performed using the Structure software (Pritchard et al. 2000). The model for sorting individuals is based on minimizing deviations from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium between loci inside the randomly built subgroups (clusters). The admixture ancestry and correlated allele frequencies models were chosen. The calculation estimates the rate of admixture (Q) for each individual (here interpreted as an introgression rate). Each run consisted of 100000 burn-in steps and 200000 MCMC steps, repeated ten times. Assignments to K subgroups (K, the number of subgroups, was tested between 1 and 12 because of the 12 lagoons involved) were evaluated with the Evanno et al. (2005) delta K method through Structure Harvester (Earl and vonHoldt 2012). The possible bias should be the large HWE deviation. For the best K, the ten runs obtained on the 512 individuals/17 samples were processed by CLUMPP software (Version 1.1.2) in order to obtain a consensus histogram.

In order to consolidate the results and obtain information on the status of potential hybrids (F1, F2, backcrosses of various classes), the same data were analysed with NewHybrids software, version 1.1 beta (Anderson and Thompson 2002). According to Anderson (2008), NewHybrids is applicable to the present situation, where there are only two diploid species that seem to be hybridizing, with some samples containing both pure and hybrid individuals. In order to provide references for pure species, the whole sample of 512 individuals was analysed in the same run. Any possible bias should be the differentiation among populations illustrated by factorial correspondence analyses (FCAs).

FCAs (Benzécri 1973) were performed using the Genetix 4.04 program (Belkhir et al. 2004), providing the overall genetic structure of the samples. This method is well adapted to genotype data. The classes are alleles and the correspondences that are calculated are the co-occurrence of two given alleles in the same individual or in a group of individuals. Clusters (clouds) detected on the diagram correspond to homogeneous genetic lineages, independent of the fish geographical origin. The inertia values (i.e. the proportion of the total information contained by an axis, expressed as a percentage) along each axis were shown to be equivalent to linear combinations of the monolocus fixation index (Fst) values (Guinand 1996). More mathematical details of the method are given in She et al. (1987).

Population genetics parameters were calculated using Genetix software. These parameters were calculated for each sample. Some samples include the two species and some hybrids. Because hybrid detection depends on the software used and, for assignment, on the threshold chosen, the whole inhabitants of a given lagoon were included in the calculations. HWE can therefore be used to detect Wahlund effect. Genetic diversity was tested with observed heterozygosity, Ho, and with unbiased estimated heterozygosity, Hnb (Nei 1978), in order to limit the effect of small sample sizes. Parameter A (mean number of alleles per locus) was also calculated. Wright’s Fis and Fst indices were determined through the f and θ estimators of Weir and Cockerham (1984). Statistical significance of the estimated values was evaluated with 5000 permutations of alleles in each population (Fis) or 5000 permutations of individuals among compared samples (Fst). Bonferroni sequential corrections (Rice 1989) were applied to multiple tests.

The Mantel test was run with Genetix, using the genetic distance based on Fst/(1-Fst), as recommended by Rousset (1997) and the geographical distance (km) computed on the map as the coastline distances between the main openings to the sea of each sampled lagoon. The value of Z, the Mantel coefficient, between the two matrices of distances was calculated with the true data, then the significance of each test was assessed by comparison with the series of pseudo-values produced by 5000 permutations of the populations’ order of one of the two matrices of distances.

Finally, genetic relationships among populations based on Reynolds’ genetic distance were graphically represented in Treeview software (Page 1996) according to a neighbour-joining (NJ) construction with a bootstrap (100 reconstructed matrices) testing the robustness of the branching as proposed in the Phylip software (Felsenstein 2004).

One of the challenges of this analysis is the detection of hybrid individuals. Using assignment (Structure software), an individual is first considered to belong to a given species when the Q value is over 0.9, which means that more than 90% of its genome would belong to this species, though other thresholds are explored, mainly 0.8. When Q is below 0.9 (or 0.8), the individual is considered hybridized. NewHybrids gives probability for an individual to belong to a pure species or to a series of hybrid categories.

All genotypes are deposited in Dryad (DOI: http://doi.org/10.5061/dryad.70rxwdbv7).

The first information sought with the Structure software was the number of subgroups contained in the whole sampling of 512 gobies. The assignment tests were performed for K=1 to K=12. The procedure of Evanno et al. (2005) designated K=2 as the more informative partition. The following partitions (K>2) brought no better likelihood and no more information (Appendix 1). The ten runs performed for K=2 were strictly identical: the assignment and the goodness-of-fit values were very similar among runs, indicating that the structure was stable and the run length sufficient (confidence interval estimated to 1/1000). Figure 2 show the consensus histogram using CLUMPP software. Berrebi et al. (2005) and Rigal et al. (2008) described the unique presence of P. microps in Mauguio Lagoon and of P. marmoratus in most of Thau lagoon, allowing identification of the two assignment clusters to one species each. The alternate presence of each species along the 12 lagoons corresponds to what was expected (Fig. 2: P. microps in red and P. marmoratus in green). The NewHybrids analysis gave exactly the same distribution (Table 4).

The Structure and NewHybrids softwares retained 268 (with Q>0.9 but 275 with Q>0.8) and 279 pure P. microps individuals, respectively. They were present in 11 lagoons of the 12 sampled ones. Very small numbers of P. microps were present in the St Cyprien, Salses-Leucate, La Palme and Bages-Sigean lagoons. However, it dominated or was exclusive in the Canet, La Peyrade, Vic, Mauguio, Vaccarès, Impériaux and Berre lagoons (Table 4). P. marmoratus was designated as a pure species in 215 (222) and 226 individuals by the two softwares, respectively. It was observed in seven lagoons with both softwares and was dominant in five of the seven (St Cyprien, Salses-Leucate, La Palme, Bages-Sigean, Thau). A very small number of individuals of P. marmoratus were found in the lagoons of La Peyrade and Vic (Table 4). Until now, the lagoon of Thau was thought to harbour only P. marmoratus (Berrebi et al. 2005). The present investigation detected the presence of pure P. microps in a tiny part of the lagoon, very close to a pure population of P. marmoratus (locations Marseillan 1 and 2 respectively).

The genetic discrimination of the two cryptic species revealed that their geographic distributions exhibit opposite patterns: P. microps in the east, P. marmoratus in the west. Two instructive exceptions were recorded: P. microps in Canet lagoon and at the Marseillan 1 station of Thau lagoon (Fig. 3). There were also intruders: (i) isolated individuals of P. microps in the west included three individuals in the St Cyprien estuary (=12%), one in the Salse-Leucate lagoon (4%), two in Bage-Sigean lagoon (5%) and two at Bouzigue station in Thau lagoon (7%); (ii) P. marmoratus in the east included three individuals in the La Peyrade lagoon (4%) and one in the Vic lagoon (5%).

While false hybrid detection cannot be rejected, two programs were used to ascertain the presence of hybrids. Obviously the assignment method is more generous, detecting 16 to 30 hybrids in the whole sample (for Q<0.8 or 0.9 respectively), while NewHybrids detected just five of them, all in the F2 category. Using assignment, among hybrids, two categories were distinguished: F1 for 0.4<Q<0.6, and backcrosses for the remaining hybrids. All NewHybrids detected hybrid individuals were among the 16 or 30 assignment hybrids. NewHybrids recognized fewer hybrids: five hybrid individuals distributed in four lagoons (Table 4). These five hybrids were all assigned to the F2 category. Structure and NewHybrids were used together to search for convergence of results obtained by two radically different methods. This was the case for species distribution, but hybrid detection was very disparate, with Structure giving three to six times more hybrids than NewHybrids. However all hybrids recognized by NewHybrids were among those found with Structure.

Heterozygosity was slightly lower in P. microps (mean 0.63, 0.59<Hnb<0.68) than in P. marmoratus (mean 0.77, 0.69<Hnb<0.81). The interspecific difference in number of alleles was in the same direction: mean values were 8.9 for P. microps and 9.6 for P. marmoratus (Table 5). Nearly all populations showed a significant departure from HWE, even after removing intruders and hybrids. The Mantel test, used to determine the effect of isolation by distance, gave no significant results for either species. The NJ construction confirmed the populations’ differentiation, with a slight geographic organization for P. marmoratus but not for P. microps (Fig. 4).

The organization of populations inside each species is described by FCA. This method was applied separately for each species, only on pure species individuals. An individual is considered here to belong to a species when the Q value is over 0.9, which is a stringent selection removing from the calculation 13 individuals whose genome is dominated by P. microps and 17 whose genome is dominated by P. marmoratus. With P. microps individuals, there was considerable multidimensional overlapping of the lagoon populations (Fig. 5), with three differentiated populations (all others are in the common centre of the diagram). Thus, Canet, Mauguio and Impériaux populations seem to represent partly differentiated groups. With P. marmoratus, only a big cloud is constituted, gathering all populations except the Thau-ULM (sample 10) into one (projection not shown).

Intraspecific variability among the lagoons can be tested through Fst estimations (Table 6). Among the eight P. microps samples, only two comparisons showed non-significant Fst: La Peyrade/Vaccarès (samples 2/12) and Impériaux/Berre (16/17). Among the nine P. marmoratus samples, most lagoons hosted a differentiated population, except seven comparisons giving a non-significant Fst, involving consecutive lagoons (stations 3/4) or consecutive stations of the Thau lagoon (stations 7/8 and 8/9) and suggesting rare local exchanges, as already observed (Berrebi et al. 2009).

The Fis estimator used to test genotypes against HWE detected an imbalance for both species in almost all lagoons (Table 5). As revealed by Micro-Checker, null alleles can be present in all loci and explain the deviations. Other technical artefacts, such as large allele drop-outs were not expected. Several explanations have been proposed, in addition to artefacts such as scoring errors or null alleles: the Wahlund effect, inbreeding, family structure and natural selection. According to Pogson et al. (1995), the Wahlund effect may be a cause when there are significant Fis deviations at more than one locus. In contrast, natural selection is less probable, as microsatellite markers are assumed to be neutral and therefore little affected by natural selection (Côrte-Real et al. 1994). Therefore, the general HWE departure appears to have a biological cause rather than being the result of artefacts. The Wahlund effect is probably the first cause in lagoons where intruders and hybrids have been found, but this disequilibrium has also been detected in lagoons inhabited by one species only (Canet and Impériaux lagoons for P. microps; Marseillan 2, Mèze and ULM stations in Thau lagoon for P. marmoratus). These significant deviations from HWE due to heterozygote deficit were documented by González-Wangüemert and Vergara-Chen (2014) in all five populations and eight microsatellites in P. marmoratus in Mar Menor coastal lagoon in SE Spain. Large deviations have also been recorded in P. minutus in the Mediterranean (Boissin et al. 2011). This imbalance for these species has been observed by some authors (Stefanni et al 2003, Pampoulie et al. 2004) but not by others (Jones et al. 2001, Gysels et al. 2004). In fact, departure from HWE is not rare in marine organisms (García De León et al. 1997, De Innocentiis et al. 2001).

When we selected the number of hybrids according to different levels of Q value, we found for the whole sampling 30 hybrids with Q>0.1, 16 with Q>0.2, 9 with Q>0.3 and 3 with Q>0.4. By comparison, NewHybrids detected five hybrids corresponding to assignment with a threshold between 0.3 and 0.4. This discrepancy is due to different software methods: Structure calculation is based on the frequency of each allele in each species and individual, while NewHybrids reconstitutes the possible genealogy among individuals after assigning each allele to a species. Both techniques confirm the presence of hybrids, but their number depends on the statistical method. When hybrids were detected, the introgression rates observed in each individual generally differed from 0.50 (Fig. 2). In other words, backcrosses (according to Structure) or F2 (according to NewHybrids) had taken place. These findings allow us to conclude that hybrids are fertile.

The history of the present lagoons and therefore of their sedentary goby populations went through two successive phases. First the lagoons, at least those between Mauguio and Thau known as Palavasian lagoons near the small city of Palavas, were isolated from the sea by a sand bar (lido), giving the lagoonal zone its particular ecology with alternate continental and marine influences. Second, the formation of each lagoon happened later, due to both marine influence and anthropic actions. The principal process was isolation of the lagoon zone from the sea. Several sedimentological investigations have attempted to date this isolation. The more conclusive are those of Sabatier et al. (2010a) and Raynal et al. (2010), which both place closure of the Palavasian coastal lagoons by the sand barrier at around 1220±120 AD (after correction). This estimation was based on the transition of gastropod shells along the sediment core, from the marine sea snail Bittium reticulatum (Cerithiidae) to the lagoonal sea snail Hydrobia acuta (Hydrobiidae). Their density shifted at about 1.80 m depth in the present sediments. This stratum was dated by Sabatier et al. (2010b) taking into account a 14C reservoir age correction of 618±30 14C yr. for the Palavasian lagoons.

The next step was dated from ancient maps of the coast. Among others (in Lenthéric 1876), those of Johannes Solivet 1570, the Marcator map of 1585 and the Nolin geographer publication of 1692 all show the Palavasian lagoons from Thau to Mauguio as a continuous ria isolated from the sea by a sand bar that is well constituted but pierced by several passages. Fragmentation of this continuum to constitute the present small lagoons has not been dated, but according to the Vidal map, it happened after 1744.

According to phylogenetic studies, the two species belong to the same evolutionary branch (Huyse et al. 2004, Mejri et al. 2011, Vanhove et al. 2012). The date of their separation has been estimated to be from 1.37 (Vanhove et al. 2012) or 2.5 (Huyse et al. 2004) to 10 Mya according to Thacker et al. (2019). While very divergent, these dates are long before the isolation and formation of the Gulf of Lion lagoons.

The genetic structure of sedentary sand gobies can therefore be described as a macroevolution, separating P. microps and P. marmoratus millions of years ago and a microevolution beginning by settlement of the two species on the Mediterranean coasts of southern France. The ria structure of the coast lasted approximately from 1200 to 1750, i.e. 5.5 centuries, during the period when sedentary gobies settled in the Palavasian lagoon. Inter-lagoon differentiation has lasted two centuries only. The progressive fragmentation, first between the Palavasian lagoons and the more eastern and southwestern lagoons, then between the Palavasian lagoons themselves, explains the observed intraspecific interpopulation differentiation (Berrebi et al. 2009, present investigation).

P. microps, the common goby, is known to be very abundant in the French Mediterranean lagoons (Quignard et al. 1984, Bouchereau et al. 1989a). As confirmation, it was present at most sites sampled here, though differences in abundance are to be noted. Thus, in some lagoons it was almost the only sedentary species of the genus (with P. minutus): i.e. Canet, La Peyrade, Mauguio, Vaccarès, Impériaux, Berre and the Marseillan 1 station in Thau lagoon. It was predominant in Vic lagoon. Its presence was anecdotal at the other stations, except St Cyprien (12%). Salinity is an environmental factor that exerts selection pressure on aquatic animals and especially fish (Poizat et al. 2004). Bouchereau et al. (1989b) and De Casabianca and Kiener (1969) described the enormous resistance of P. microps to extreme environmental conditions and variations, as in the lagoons of Mauguio and Palo (Corsica). Rigal et al. (2008) demonstrated that P. microps is able to survive in absolutely fresh water, whereas P. marmoratus cannot.

Despite the phylogenetic tree which showed low bootstrap values for P. microps branches and no geographic logic (Fig. 4), and despite the absence of IBD according to the Mantel test, the multidimensional analysis, performed on pure P. microps individuals only (Q>0.9), was able to distinguish several populations (Fig. 5). A major genetic form, positioned at the centre of the graph, gathered populations from five main lagoons: La Peyrade, Vic, Vaccarès and Berre, together with the Marseillan 1 station of Thau lagoon. The P. microps populations of the Canet, Mauguio and Impériaux lagoons showed a slightly different genetic composition from that of all other lagoons. These results reflect a structure in meta-populations of the P. microps species in the Gulf of Lions, as a consequence of recent simultaneous lagoon isolation affecting the differentiation of gobies populations with no geographic logic, as observed on mitochondrial sequences (Berrebi et al. 2009). There is no specific information published about the pelagic larval duration (PLD) of sedentary gobies (González-Wangüemert and Vergara-Chen 2014), although congeneric species have a PLD of 30 to 39 days under laboratory conditions (Fonds 1973). In the few studied Pomatoschistus species, very limited migratory behaviour has been observed (Gysels et al. 2004, Berrebi et al. 2005). Their poor swimming abilities are possibly partly due to the shape of their pelvic fins, which constitute a suction disc typical of fixed animals (Miller 1986, Bardin and Pont 2002).

The marbled goby, P. marmoratus, showed a single large cluster in FCA analysis (diagram not shown), with only one sample differentiated from the remainder (the ULM station in Thau lagoon), with no clear explanation for this. Distribution of P. marmoratus on the French Mediterranean coast has been little studied (De Casabianca and Kiener 1969). The species has been reported in Thau lagoon and in the Vaccarès complex (Berrebi et al. 2005) and far in the southwest, in the Mar Menor lagoon in Spain (González-Wangüemert and Vergara-Chen 2014). The present study widens our knowledge of its range in the lagoons of the French Mediterranean coast. It was detected in seven lagoons and constituted the main sedentary Pomatoschistus species in five of them (St Cyprien, Salses-Leucate, La Palme, Bages-Sigean and Thau). In other lagoons, it was represented by a very small number of individuals, frequently between 1% and 5% and always below 10% (Table 4). The lagoons in which this species predominated are marinized (Table 2). The mouth of the Tech River (St Cyprien station) is strongly influenced by marine waters. At Bages-Sigean lagoon, salinity does not drop below 13.7‰ (Labat 1976). Rigal et al. (2008) demonstrated experimentally that the marbled goby cannot survive in a completely desalinated environment, as individuals die at the threshold of 1.4‰.

The differential distribution of the two sedentary Pomatoschistus species is far from presenting a random pattern with P. marmoratus in the west and P. microps in the east. This geographic distribution suggests a historical event, i.e. opposite migrations with contact at the middle. However, this process cannot wholly explain why the two hypothetic directional migrations ceased between the Thau and La Peyrade lagoons (Fig. 3) or some exceptions in the geographic picture, i.e. a P. microps population in the west (Canet lagoon) and a P. microps population (Marseillan 1 station) within a P. marmoratus lagoon (Thau). Concerning this last surprising phenomenon, the resurgence of freshwater springs in the Thau lagoon (Bakalowicz et al. 2003, Elbaz-Poulichet et al. 2005, Plus et al. 2006) might be an explanation because P. marmoratus does not survive in freshwater (Rigal et al. 2008). Other surprising features are the very limited intruders in both directions and relatively few fertile hybrids.

Absence of P. marmoratus can be explained in lagoons where freshwater crises occur due to intense rains. These include the Canet, Mauguio, Thau, La Peyrade, Vic, Vaccarès-Impériaux and Berre lagoons, i.e. the whole eastern part of the sampled region. The absence of P. marmoratus should be linked to its incapacity to survive in freshwater conditions. Its preference for marine water has been known for a long time (De Casabianca and Kiener 1969). The absence of the common goby, P. microps, needs some more hypotheses to be explained. This species has been detected in the west as hybrids in La Palme and as individual intruders in the St Cyprien, Salses-Leucate, Bages-Sigean lagoon and at the Thau-Bouzigues station. However, as far as we know, these are not inhospitable environments for the common goby (Rigal et al. 2008). A possible explanation could be that the marbled goby is a better competitor in these lagoons for an unknown ecological reason, which would exclude P. microps from these sites.