Distribution and hybridization of two sedentary gobies (Pomatoschistus microps and Pomatoschistus marmoratus) in the lagoons of southern France ; Distribución e hibridación de dos gobios sedentarios (<em>Pomatoschistus microps</em> y <em>Pomatoschistus marmoratus</em>) en las lagunas del sur de Francia

Pomatoschistus marmoratus and Pomatoschistus microps are small sedentary gobies inhabiting the lagoons of European Mediterranean and Atlantic coasts. Along the French Mediterranean coast their respective geographical distribution is not precisely known, in part because they are cryptic species. In this study, 512 gobies of both species were caught as 17 samples in 12 lagoons of the Gulf of Lion on the French Mediterranean coast. They were genotyped at six microsatellite loci and investigated statistically using multidimensional analyses, Bayesian assignment (Structure) and NewHybrids classification. This allowed the contrasted distribution of each species (P. microps in the east, P. marmoratus in the west) to be described, with several exceptions. Neither geographic structure nor isolation by distance was detected among differentiated populations of each species. The suggested mechanism is a deep sedentary behaviour associated with foundations following extinctions. The two species are sympatric or even in syntopy in five or six sampled lagoons producing rare fertile hybrids.


INTRODUCTION
Gobies are teleost fish of the order Perciforms and the family Gobiidae (Miller 1986). Nearly 2000 species of this family have been identified worldwide (Miller 1984, Nelson 2006 in all aquatic environments: continental, brackish and marine. In the Mediterranean, Miller (1986) considered 44 species described in marine and lagoon environments. This number climbed to 60 species in 25 genera (Quignard and Tomasini 2000) and then to 62, 26 of them endemic to the Mediterranean (Ahnelt andDorda 2004, Engin andSeyhand 2017).
The monophyletic "sand gobies" clade, morphologically recognizable by the head canals and vertebrae characteristics, includes the genus Pomatoschistus Gill, 1863, among four paraphyletic genera. Thirteen Pomatoschistus species have been described, all less than 10 cm in length, the smallest being P. nanus (Engin and Seyhand 2017). Based on the pairwise distance matrix estimated with nuclear DNA (ITS1 locus) and mitochondrial DNA (12S and 16S fragments), Huyse et al. (2004) showed that P. microps, P. marmoratus and P. knerii clustered separately and in parallel with the P. minutus complex. P. microps and P. marmoratus are considered to have appeared in the Mediterranean ). However, recent molecular data suggest another story for P. microps. Namely, this species could have originated from Portuguese ancestral populations in glacial refuges. The species might have colonized northern Europe first and then the Mediterranean Sea (Tougard et al. 2014).
P. microps is present from the Moroccan Atlantic coast to the Baltic and to the Mediterranean North West (Miller 1986). Gene flow is almost non-existent between the Atlantic and the Mediterranean populations, which are clearly differentiated into two lineages  with different life-history characters (Healey 1972, Bouchereau andGuélorget 1997). Locally, there is evidence for slight genetic differentiation between close lagoons (Berrebi et al. 2009). It is the most abundant goby in Mediterranean lagoons with P. minutus (Quignard et al. 1984, Bouchereau et al. 1989a. P. marmoratus is distributed on the Atlantic west coast of the Iberian Peninsula, the north coasts of the western and eastern Mediterranean, the Black Sea, the Azov Sea and Lake Qarun, Egypt, where it was introduced via the Suez Canal (Miller 1986), as well as along the southern coasts of the Mediterranean, but limited to the Tunisian lagoons (Mejri et al. 2009). Overall, a west/east genetic structure is hinged on either side of the Sicilian-Tunisian strait (Mejri et al. 2011). An intra-lagoon genetic structure has been observed in the Mar Menor coastal lagoon, SE Spain (González-Wangüemert and Vergara-Chen 2014), together with a high gene flow along the coast (Vergara-Chen et al. 2010). It is a sedentary species of Mediterranean lagoons little known because of the difficulty of discriminating it from the well-known P. microps.
Both species live primarily on sandy substrates. P. microps is a consumer of meiofauna, i.e. mainly planktonic and benthic crustaceans such as copepods and amphipods (De Casabianca andKiener 1969, Leclerc et al. 2014). They live in very shallow lagoon zones, mainly on the shores, thus constituting a source of food for many species of the macrofauna. This gives them a key role in energy transfer from meiofauna to fish and birds (Pampoulie 1999, Ray 2005. Each species is a batch spawner with one or two spawning seasons (Bouchereau et al. 1989b). Because it is strictly sedentary, the fertility and growth of these species of gobies are greatly affected when environmental conditions become unfavourable (Pampoulie et al. 2000).
Based on discriminating allozymic markers, P. marmoratus has been shown to represent over 90% of gobies of the Thau lagoon, while Mauguio Lagoon seems to be exclusively inhabited by P. microps (Berrebi et al. 2005). These two lagoons will be considered as references for each species in the present study.
Both species are called cryptic because it is very difficult to determine which is which using morphological criteria (Webb 1980). The safest test able to distinguish the two species morphologically is based on the presence (in P. marmoratus) or absence (in P. microps) of the mucosal occulo-scapular channel, which necessitates unstable staining protocols to be observed (Huyse et al. 2004, Rigal et al. 2008. The relatively close positions of P. microps, Kroyer, 1838 (common goby) and P. marmoratus, Risso, 1810 (marbled goby) in the phylogeny and the fact that the two species coexist in sympatry in lagoons of southern France explains why they hybridize occasionally in a southern France lagoon complex: the Vaccarès and Impériaux lagoons (Berrebi et al. 2005), corresponding to stations 15 and 16 of the present survey (Table 1). This discrete hybridization between small and cryptic species is neither really understood nor explained and deserves a detailed regional investigation on the distribution of the two species, their hybrids and the ecological conditions of their hybridization. Because southern France is the only zone where P. microps and P. marmoratus are in sympatry (together with P. minutus), and because this is the first zone where hybridization between these two species has been recorded, the present investigation aimed to determine precisely the reciprocal distribution of the two species and of their hybrids on the French coasts of the Mediterranean Gulf of Lyon. For this study 12 lagoons were selected along the French Mediterranean coast, from the Spanish border to Marseille (Fig. 1), for a total of 512 fishes. Gobies were genotyped at six microsatellite loci and the data were analysed using several statistical methods.

Sampling
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.

The lagoon ecology
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.

Microsatellite genotyping
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 2 , 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.

Statistical methods
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 link-age 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 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 Table 2. -Short description of the ecological characteristics of the lagoons (according to Labat 1976, Quignard et al. 1983, Pampoulie 2001, Ifremer 2002, 2003 Jones et al. 2001 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.
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   Table 1).
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%).

Hybrid detection
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.

Organization of species populations
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).  (Nei 1978); Ho, observed heterozygosity; A, mean number of allele by locus; Fis, f estimator of Fis (Weir and Cockerham 1984); signif., significance of departure of the Fis estimator from zero (5000 permutations and Bonferroni correction). The p value is p<0.01 for a highly significant departure (**) and p<0.001 for a very highly significant departure (***).

Map
Lagoon

DISCUSSION
The efficiency of the six microsatellite markers used was high enough to assign each individual to one of the two pure species or to different classes of hybrids, as seen in Figure 2 and Table 4. The presence of two species (P. microps and P. marmoratus) in the sampled lagoons has been demonstrated earlier. They were discovered and described using microsatellites (Berrebi et al 2005), mtDNA (Berrebi et al 2009, Tougard et al 2014, adaptation (Rigal et al 2008) and morphology (Travers et al 2011). Therefore, the detection by assignment and NewHybrids of two taxa in the samples corresponds necessarily to these two species, despite the partially overlapping allele sizes recorded between them (Table 3).

Hardy-Weinberg equilibrium departure
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 Berre -P. microps 0 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 but not by others (Jones et al. 2001). In fact, departure from HWE is not rare in marine organisms (García De León et al. 1997, De Innocentiis et al. 2001).

Hybridization
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.

Macro-and microevolution in lagoonal Pomatoschistus gobies
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 Bit-tium 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 14 C reservoir age correction of 618±30 14 C 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).

Distribution, differentiation and adaptation of P. microps
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 , 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 andPont 2002).

Distribution, differentiation and adaptation of P. marmoratus
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‰.

Differential distributions, historic events and competitive exclusion
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.

CONCLUSIONS
A new picture of Pomatoschistus inhabitants of southern French lagoons is now available. Each Mediterranean coastal lagoon is a site of exclusion, competition or hybridization between two cryptic species of sand gobies. While composed of differentiated populations, neither species shows a geographic structure in the NJ tree (Figs 4 and 5) or IBD. This reinforces the idea that the species exhibit a very sedentary behaviour. This status makes them dependent on extinction/ foundation processes. For example, a major dystrophic crisis in August 1988 resulted in the disappearance of resident fish populations in the Mauguio Lagoon, including P. microps. Recovery of the P. microps population took several years and was carried out from a few lagoon survivors or from external intakes from healthy biotopes (Bouchereau et al. 1990). These random recolonizations would explain the phylogeographic "disorder" observed for both species. APPENDIX Appendix 1. -DeltaK calculations according to the Evanno et al. (2005) method, through Structure Harvester online software (Earl vonHoldt 2012). Here K=2 (in bold type) is the uppermost hierarchical level of structure.