Population Genetic Diversity of Two Marine Gobies (Gobiiformes: Gobiidae) from the North-Eastern Atlantic and the Mediterranean Sea

: Gobies (Gobiiformes: Gobiidae) are the most species-rich family of ﬁshes in general, and the most abundant ﬁsh group in the European seas. Nonetheless, our knowledge on many aspects of their biology, including the population genetic diversity, is poor. Although barriers to gene ﬂow are less apparent in the marine environment, the ocean is not a continuous habitat, as has been shown by studies on population genetics of various marine biota. For the ﬁrst time, European marine goby species which cannot be collected by common ﬁshery techniques were studied. The population genetic structure of two epibenthic species, Gobius geniporus and Gobius cruentatus , from seven localities across their distribution ranges was assessed, using one mitochondrial (cytochrome b) and one nuclear gene (ﬁrst intron of ribosomal protein gene S7). Our results showed that there is a great diversity of haplotypes of mitochondrial gene cytochrome b in both species at all localities. Global ﬁxation indices (F ST ) indicated a great di ﬀ erentiation of populations in both studied gobies. Our results did not show a geographic subdivision to individual populations. Instead, the data correspond with the model of migration which allow divergence and recurrent migration from the ancestral population. The estimated migration routes coincide with the main currents in the studied area. This matches well the biology of the studied species, with adults exhibiting only short-distance movements and planktonic larval stages.


Introduction
Despite the fact that the seas and oceans are interconnected, marine organisms can show a strong genetic differentiation [1]. Several phylogeographical studies have shown that even the ocean recolonizing the Atlantic when more favourable temperatures were re-established during interglacial phases like the present one [29].
The studied area concerns the Mediterranean Sea province and the South European Atlantic Shelf ecoregion. The Mediterranean Sea is classified as a specific, well defined biogeographic province within the Temperate Northern Atlantic realm [30]. It is further divided into seven ecoregions (Adriatic Sea, Aegean Sea, Levantine Sea, Tunisian Plateau/Gulf of Sidra, Ionian Sea, Western Mediterranean, and Alboran Sea). Its neighbouring marine provinces are the Black Sea in the east and the Lusitanian province in the Atlantic Ocean in the west. The Mediterranean Sea is connected with the South European Atlantic Shelf and Saharan Upwelling ecoregions (both belonging to the Lusitanian province and having the boundary in the Gibraltar Strait area).
Circulation of water in the studied area, and especially in the Mediterranean Sea, is very complex. The main currents and gyres are depicted in Figure 1. Within the studied region, multiple biogeographical barriers have been identified, of which the Almeria-Oran front, the Strait of Sicily and the Otranto Strait (see Figure 1) are considered to be the major ones influencing genetic diversity of various marine organisms [1,13]. However, dissimilar influence of biogeographic barriers has been found even for the closely related taxa with the same biology and ecology [31][32][33][34]. The studied area concerns the Mediterranean Sea province and the South European Atlantic Shelf ecoregion. The Mediterranean Sea is classified as a specific, well defined biogeographic province within the Temperate Northern Atlantic realm [30]. It is further divided into seven ecoregions (Adriatic Sea, Aegean Sea, Levantine Sea, Tunisian Plateau/Gulf of Sidra, Ionian Sea, Western Mediterranean, and Alboran Sea). Its neighbouring marine provinces are the Black Sea in the east and the Lusitanian province in the Atlantic Ocean in the west. The Mediterranean Sea is connected with the South European Atlantic Shelf and Saharan Upwelling ecoregions (both belonging to the Lusitanian province and having the boundary in the Gibraltar Strait area).
Circulation of water in the studied area, and especially in the Mediterranean Sea, is very complex. The main currents and gyres are depicted in Figure 1. Within the studied region, multiple biogeographical barriers have been identified, of which the Almeria-Oran front, the Strait of Sicily and the Otranto Strait (see Figure 1) are considered to be the major ones influencing genetic diversity of various marine organisms [1,13]. However, dissimilar influence of biogeographic barriers has been found even for the closely related taxa with the same biology and ecology [31][32][33][34]. Gobies (Gobiiformes: Gobiidae), with currently recognized 1915 valid species, are the most species-rich family of fishes in the world [35]. It is also the most speciose family of fishes of the European seas and, in particular, of the Mediterranean Sea. The majority of gobies are benthic fishes with planktonic larval stages. Marine gobies predominantly occupy shallow shelf bottoms, with limited number of species extending to deeper shelf, and only a few of them reaching bathyal depths [36,37]. In total, there are more than 90 marine species of gobies in European seas listed in the last review [38]. However, new species are still being discovered, e.g., [36,37,39,40], and the knowledge about the distribution of many species is still quite limited e.g., [41][42][43]. The information about the population genetic diversity and phylogeography of European marine gobies is very poor too. Population genetic studies have been done so far in only a few species, which are easy to collect by the common fishery techniques. Nevertheless, on the basis of limited information, it is evident that for some studied goby species there is a clear genetic differentiation between distant populations [44][45][46][47][48][49][50], but not for the others [51]. Within the genus Gobius, the population genetic structure was studied so far in the only species, Gobius niger, where the mitochondrial marker nicotinamide adenine dinucleotide dehydrogenase (NADH) was analysed by restriction fragment length polymorphism (RFLP) [51]. The results suggested the existence of population subdivision in this species. Gobies (Gobiiformes: Gobiidae), with currently recognized 1915 valid species, are the most species-rich family of fishes in the world [35]. It is also the most speciose family of fishes of the European seas and, in particular, of the Mediterranean Sea. The majority of gobies are benthic fishes with planktonic larval stages. Marine gobies predominantly occupy shallow shelf bottoms, with limited number of species extending to deeper shelf, and only a few of them reaching bathyal depths [36,37]. In total, there are more than 90 marine species of gobies in European seas listed in the last review [38]. However, new species are still being discovered, e.g., [36,37,39,40], and the knowledge about the distribution of many species is still quite limited e.g., [41][42][43]. The information about the population genetic diversity and phylogeography of European marine gobies is very poor too. Population genetic studies have been done so far in only a few species, which are easy to collect by the common fishery techniques. Nevertheless, on the basis of limited information, it is evident that for some studied goby species there is a clear genetic differentiation between distant populations [44][45][46][47][48][49][50], but not for the others [51]. Within the genus Gobius, the population genetic structure was studied so far in the only species, Gobius niger, where the mitochondrial marker nicotinamide adenine dinucleotide dehydrogenase (NADH) was analysed by restriction fragment length polymorphism (RFLP) [51]. The results suggested the existence of population subdivision in this species.
In this work we studied the population genetic diversity of Gobius geniporus and Gobius cruentatus (Gobius-lineage, gobiine-like clade sensu Agorreta et al. [52]), two bottom dwelling (epibenthic) species. Gobius cruentatus occurs in the north-eastern Atlantic Ocean, from south-western Ireland to the coasts of Senegal, and in the Mediterranean and the Black Seas [53,54], while G. geniporus is a Mediterranean endemic [53]. Gobius cruentatus is typically detected on mixed bottom habitats dominated by stones, boulders or seagrass [55]. It can grow up to 18 cm and occurs in depths between 1.5 and 40 m [53,54]. Gobius geniporus prefers sandy bottoms mixed with gravel, cobbles and boulders, with at least some amount of rocky formations present scattered over sand. It reaches a size of 16 cm and the depth at which it was observed ranges from 1 to 30 m [53,56,57]. These species are restricted to shallow shelf bottoms, so their real area of occupancy is only a narrow stripe of bottom along the coastline. Being small, epibenthic, territorial and non-migratory in adulthood, these species are expected to have only a limited dispersal capability [53].
The aim of this study was to assess the genetic diversity of seven geographically distant populations of two species of European marine gobies, G. geniporus and G. cruentatus, across their distribution ranges, using mitochondrial and nuclear markers, in order to reveal a possible population subdivision and a potential existence of biogeographical barriers which would affect the connectivity of the populations.

Samples
It is very difficult to collect benthic marine gobies in general, unless they live on sandy or muddy bottoms. They are usually not commercially used, so it is not possible to purchase them or to obtain them from fishermen even as a bycatch. As both species investigated in this study typically occupy mixed bottom habitats it is not possible to use common fishery techniques, such as trawl nets or push-nets, to collect them, even though these two species belong to the largest ones among gobies from the north-eastern Atlantic and the Mediterranean Sea. For research purposes, these two species are being collected individually, during scuba-diving and using hand nets and anaesthetic, which is very time consuming. They do not occur in shoals and it can be difficult to spot them. A total of 74 specimens of G. geniporus from seven localities in the Mediterranean Sea and 41 specimens of G. cruentatus from two localities in the Atlantic Ocean and five localities in the Mediterranean Sea were included in this study (see Figure 1 and Table 1). Tissue samples and voucher specimens are deposited in the ichthyological collection of the National Museum, Prague, Czech Republic.

DNA Extraction, Amplification, Sequencing
Total genomic DNA was extracted from the finclips using Geneaid ® DNA Isolation Kit following the manufacturer's protocol. The samples were analysed for two genes, one mitochondrial, cytochrome b (cyt b), and one nuclear, first intron of the ribosomal protein gene S7 (S7). For the amplification of cyt b, the primers GluF and ThrR [58] were used. S7 was amplified with the primers S7RPEX1F and S7RPEX2R [59]. The polymerase chain reaction (PCR) was performed in 25 µL total volume containing 12.5 µL of PPP Master Mix (TopBio), 9.7 µL of Ultrapure H 2 O, 0.65 µL of each primer and 2 µL of DNA isolate. Amplification of cyt b followed the protocol described in Šanda et al. [60]. For S7, a specific touch-down protocol was used with the following steps: initial denaturation at 95 • C for 5 min, followed by 5 cycles of denaturation, annealing, and elongation: 94 • C for 40 s, 60 • C for 1 min, 72 • C for 2 min, followed by 35 cycles of denaturation, annealing, and elongation: 95 • C for 30 s, 56 • C for 1 min, 72 • C for 2 min, and the final elongation at 72 • C for 20 min. PCR products were purified with the use of ExoSAP-IT and sequenced at Macrogen Europe. For the sequencing of cyt b, the specific internal primers were designed: GcruF1 (5 -GGT GCA ACC GTC ATC ACT AA-3 ) and GcruR1 (5 -AGT GGG TTG GCA GGA ATG-3 ) for G. cruentatus and GgenF1 (5 -GTA GGC TAT GTC CTG CCC TG AG-3 ) and GgenR1 (5 -TTG GAG CCT GTC TCG TG GA-3 ) for G. geniporus. Nuclear gene S7 was sequenced using the amplification primers. Sequences were deposited in GenBank under accession numbers MT774412-MT774485 (cyt b) and MT893746-MT893891 (S7) for G. geniporus and MT684467-MT684507 (cyt b) and MT684508-MT684585 (S7) for G. cruentatus.

Data Analyses
Obtained cyt b and S7 sequences were checked manually in Chromas v2.6.4 and aligned in Bioedit v7.2.6.1 [61]. The appropriate model of nucleotide substitution was determined using jModelTest v2.1.9 [62], based on Akaike Information Criterion (AIC) [63]. DnaSP v6.11.01 [64] was used to assess the haplotype diversity (Hd) and nucleotide diversity (π), as well as to perform Fu and Li's F and Tajima's D neutrality tests. The results of these tests can point to possible selection or a change in population demography. To evaluate the amount of genetic variance within and between populations, analysis of molecular variance (AMOVA) was performed using ARLEQUIN v3.5 [65]. It estimates population differentiation with the use of individual haplotypes and their frequency in the studied populations. This further enables calculation of fixation indices, a global F ST and pairwise F STs . F ST expresses a degree of genetic differentiation between the individual populations. The two populations from Cyprus were grouped together, as well as the two populations from the Adriatic Sea (Montenegrin and Croatian); the grouping was made based on location, proximity and the water circulation. The remaining populations represented individual groups. The statistical significance of the F ST values was tested by executing 16,000 permutations. For pairwise F STs , a Bonferroni correction was subsequently applied to correct for multiple tests. Further, genetic distances (uncorrected p-distances) between and within populations of each species were calculated in MEGA 6 [66]. The datasets of S7 were phased by the program PHASE v2.1.1 [67]. All sequences were phased with a probability of 0.9 and the final datasets with inferred phased sequences consisted of 146 sequences for G. geniporus and 78 for G. cruentatus. The phased S7 data were then used for calculating diversity measures and constructing haplotype networks. The rest of the analyses were not performed on S7 due to a very low polymorphism of S7 datasets. A detailed reconstruction of relationships of the haplotypes of populations was performed by a statistical parsimony method under a 95% connection limit [68], using PopART [69].
Isolation by distance hypothesis was tested by Mantel test [70] using R v3.5 software (package adegenet), executing 1000 permutations. Mantel test compares genetic distances estimated by pairwise F STs with geographical distances between locations. The matrices of geographical distances were derived from the coordinates of the individual localities. Another approach was applied using ARLEQUIN v3.5 [65], where the shortest marine paths between each pair of localities, estimated from the GoogleEarth, were used in matrices of geographical distances. Statistical significance of the Mantel test was estimated by executing 1000 permutations.
We estimated migration routes between pairs of populations using Migrate-n software v4.4. [71]. Five demographic scenarios were tested for each population pair: (1) model assuming full migration between populations, (2) model assuming migration from the population A to the population B, (3) model assuming migration from B to A, (4) model allowing divergence, where A splits off from B, with migration from B to A after the split, (5) model allowing divergence, where B splits off from A, with migration from A to B after the split. For each scenario, the migration rate and population sizes were estimated; for the scenarios (4) and (5) also the time of divergence between the two populations. Pairs of populations are listed in the Table S1. Migrate-n analyses were conducted using a static heating strategy with four short chains with temperature values of 1.0, 1.5, 3.0, and 1.0 × 10 6 and a single long chain. 1,000,000 steps were recorded every 100 generations with 200,000 steps discarded as burn-in to ensure the convergence of the analyses. Appropriate mutation model was assessed using jModelTest [62] resulting in Hasegawa-Kishino-Yano (HKY). Priors were set as follows: Bayes-priors = THETA * * UNIFORMPRIOR: 0.001, 0.000 0.0100, Bayes-priors = MIG * * UNIFORMPRIOR: 0.000, 100000.000, 10000.000.
Past population demography of each species was inferred using the linear Bayesian skyline plot model [72], implemented in BEAST v1.8.4 [73]. It allows observing fluctuations of effective population sizes from the present, backwards in time, to the coalescence in the most recent common ancestor, and is expressed graphically. Analyses were conducted under the Bayesian coalescent method, with corresponding nucleotide substitution model for each species and using a strict molecular clock. The x-axis of the plot shows the time in mutation units per nucleotide position and y-axis scaled effective population size. Simulations ran for 100 million Markov chain Monte Carlo (MCMC) steps with sampling every 10,000th generation. Results from three independent runs were combined using LogCombiner and burn-in was set to 20 million iterations in each run. Finally, TRACER v1.7.0 [74] was used to check the parameter estimates and visualize Bayesian skyline plots.

Gobius geniporus
In G. geniporus, 74 specimens were analysed ( Table 1). The alignment of cyt b had a length of 1113 bp and contained 56 polymorphic sites, while there were only two variable sites in the 594 bp long alignment of 146 sequences of S7 (Table 2). A total of 45 haplotypes were found for cyt b and only three haplotypes for S7 within seven Mediterranean populations. The best-fit substitution model selected for cyt b was general time reversible with proportion of invariable sites (GTR+I). Haplotype diversity of cyt b was high, while nucleotide diversity low (Hd = 0.969; π = 0.004) and for S7 both haplotype and nucleotide diversity were extremely low (Hd = 0.054; π = 0.0001) ( Table 2). Diversity measures calculated per each locality are listed in Table S2. The values of neutrality tests (Tajima's D, Fu and Li's F) for cyt b were negative and significant, indicating a recent population expansion or purifying selection ( Table 2). The Bayesian skyline plot of G. geniporus depicts a gradual population size growth since the coalescence and its stabilization in the present (Figure 2a).  In the haplotype network of G. geniporus based on cyt b, there was an indication of a certain geographical pattern: there were two major haplotype groups. In one of them there was one more frequent haplotype shared between the central Mediterranean populations (Montenegrin, Croatian and Sicilian), while in the other one, there were two more frequent haplotypes, shared between Cypriot and Greek populations (Figure 3a and Figure S1a). Most haplotypes from Cypriot populations grouped together, as well as the majority of the haplotypes from the Sicilian one. Unique haplotypes prevailed in the network. Practically all sequences of S7 in G. geniporus were of the same haplotype (142 out of 146, Figure 3a) and due to this low polymorphism, all following analyses were performed only on the cyt b dataset. AMOVA for G. geniporus based on the cyt b showed that most of the genetic variance is distributed within populations (76%). FST index indicated a high level of genetic differentiation (FST = 0.237, p < 0.01). Pairwise FSTs showed in most cases a pronounced or high level of genetic differentiation between the pairs of populations, but several values were low (Table  3); however, only a half of values were significant. Statistically significant values indicating high or pronounced differentiation were for most comparisons for Sicilian and both Cypriot populations. Mean p-distances between the populations were low and ranged between 0.2 and 0.6% (Table 3). Mean p-distances within populations were of a similar range (0.1-0.5%), while the maximum intraspecific p-distance for G. geniporus was 1.08%.  In the haplotype network of G. geniporus based on cyt b, there was an indication of a certain geographical pattern: there were two major haplotype groups. In one of them there was one more frequent haplotype shared between the central Mediterranean populations (Montenegrin, Croatian and Sicilian), while in the other one, there were two more frequent haplotypes, shared between Cypriot and Greek populations (Figure 3a and Figure S1a). Most haplotypes from Cypriot populations grouped together, as well as the majority of the haplotypes from the Sicilian one. Unique haplotypes prevailed in the network. Practically all sequences of S7 in G. geniporus were of the same haplotype (142 out of 146, Figure 3a) and due to this low polymorphism, all following analyses were performed only on the cyt b dataset. AMOVA for G. geniporus based on the cyt b showed that most of the genetic variance is distributed within populations (76%). F ST index indicated a high level of genetic differentiation (F ST = 0.237, p < 0.01). Pairwise F STs showed in most cases a pronounced or high level of genetic differentiation between the pairs of populations, but several values were low (Table 3); however, only a half of values were significant. Statistically significant values indicating high or pronounced differentiation were for most comparisons for Sicilian and both Cypriot populations. Mean p-distances between the populations were low and ranged between 0.2 and 0.6% (Table 3). Mean p-distances within populations were of a similar range (0.1-0.5%), while the maximum intraspecific p-distance for G. geniporus was 1.08%. Table 3. Mean genetic distances between Gobius geniporus populations for cytochrome b (uncorrected p-distances, in %, above the diagonal), intrapopulation distances (on diagonal), and pairwise F STs (below diagonal). Significant values of F STs (at α = 0.05/number of pairs) indicated by asterisk.  Mantel test was significant using both approaches (adegenet: observation 0.299, expectation 10 −5 , p-value < 0.001; ARLEQUIN: correlation coefficient 0.69, p-value < 0.001), indicating a possible pattern of isolation by distance (see Figure 4a). Among the modelled migration scenarios, for each pair of populations, the model which allows divergence and the recurrent immigration from the ancestral population after the split was the one with the highest probability. The divergence directions and the migration routes are schematically depicted in the Figure 5a, while the estimates of immigration rates, divergence times and population sizes are listed in the Table S3. The system of migration routes is rather circular, anticlockwise, with a large circle between Sicily, western Cyprus, eastern Cyprus, Greece, Montenegro, Croatia and Sicily, and two smaller ones: Sicily, Greece, Montenegro, Croatia, Sicily, and Sicily, Montenegro, Croatia, Sicily. All the routes, with the exception of the one between Sicily and Greece eastwards, can be well explained by the prevailing currents (see Figure 5a). The highest rate of migration among the modelled pairs of populations was estimated between the western and eastern Cyprus, correspondingly with their proximity and the prevailing eastward current. Mantel test was significant using both approaches (adegenet: observation 0.299, expectation 10 −5 , p-value < 0.001; ARLEQUIN: correlation coefficient 0.69, p-value < 0.001), indicating a possible pattern of isolation by distance (see Figure 4a). Among the modelled migration scenarios, for each pair of populations, the model which allows divergence and the recurrent immigration from the ancestral population after the split was the one with the highest probability. The divergence directions and the migration routes are schematically depicted in the Figure 5a, while the estimates of immigration rates, divergence times and population sizes are listed in the Table S3. The system of migration routes is rather circular, anticlockwise, with a large circle between Sicily, western Cyprus, eastern Cyprus, Greece, Montenegro, Croatia and Sicily, and two smaller ones: Sicily, Greece, Montenegro, Croatia, Sicily, and Sicily, Montenegro, Croatia, Sicily. All the routes, with the exception of the one between Sicily and Greece eastwards, can be well explained by the prevailing currents (see Figure 5a). The highest rate of migration among the modelled pairs of populations was estimated between the western and eastern Cyprus, correspondingly with their proximity and the prevailing eastward current.

Gobius cruentatus
Mitochondrial marker cyt b was analysed in 41 specimens of G. cruentatus, and nuclear marker S7 in 39 specimens (Table 1). The alignment length of cyt b was 1117 bp and contained 47 polymorphic sites, while that of S7 had a length of 555 bp and contained only three segregating sites (Table 2). There were 32 haplotypes of cyt b and only four of S7. The best-fit substitution model selected for cyt b was general time reversible with proportion of invariable sites (GTR + I). An overall haplotype diversity of cyt b was high, while nucleotide diversity low (Hd = 0.985; π = 0.006), which can indicate a recent population expansion. This was also suggested by negative values of neutrality tests (see Table 2). On the other hand, for nuclear gene S7, values of haplotype and nucleotide diversity were markedly lower (Hd = 0.212; π = 0.0004). Diversity measures calculated per each locality are listed in Table S2. The Bayesian skyline plot showed a constant population size in the past, followed by a gradual population expansion and a stable population in the present (Figure 2b).
In G. cruentatus, cyt b haplotype network reconstruction did not reveal any well-defined spatial structure (Figure 3b). The network consisted mostly of unique haplotypes, very diverse for each locality ( Figure S1b). Three haplotypes were shared between two or three distant populations (Spain/Sicily, Sicily/Croatia/Cyprus, Cyprus/Croatia). The network displayed two highly variable haplogroups separated by six mutational steps. Interestingly, there was a particular geographic pattern: while haplotypes from the central Mediterranean populations (Croatia, Sicily and Montenegro) were in both haplogroups, the haplotypes from the westernmost populations (Portugal, Spain and France) were placed only in one haplogroup, and haplotypes from the easternmost population, Cyprus, occurred only in the other haplogroup. On the contrary, there was nearly no polymorphism in S7. The network was formed by one dominant haplotype which included 69 alleles (out of 78) and was shared among all populations, one less frequent shared haplotype, which included only 7 alleles, and two unique ones. All haplotypes were very similar (Figure 3b). The remaining analyses were performed only on the cyt b dataset due to the low polymorphism in S7.
Similarly to the results of haplotype network, AMOVA performed on the cyt b showed that most of the genetic variance is distributed within the populations, with a ratio of approximately 1:4 between the variability among and within populations. Computed F ST index indicated a high level of genetic differentiation (F ST = 0.216, p < 0.01). Some values of pairwise F STs performed on G. cruentatus indicated a pronounced differentiation (the highest values were between the easternmost population, Cyprus, and three westernmost populations-Spain, Portugal and France), but the majority were low or moderate; however, none of the values was significant ( Table 4). The mean p-distances between the populations ranged between 0.3 and 1%, and the range of intrapopulation p-distances was similar (0.2-1.1%) ( Table 4). The highest interpopulation divergences were between the westernmost and easternmost populations (0.9-1%, see Table 4). The highest overall intraspecific p-distance for G. cruentatus was 1.52%. Table 4. Mean genetic distances between Gobius cruentatus populations for cytochrome b (uncorrected p-distances, in %, above diagonal), intrapopulation distances (on diagonal), and pairwise F STs (below diagonal). None of the values of F STs was significant at α = 0.05/number of pairs. The result of the Mantel test was significant using both approaches (adegenet: observation 0.085, expectation 0.002, p-value 0.03; ARLEQUIN: correlation coefficient 0.81, p-value < 0.01), indicating a possible pattern of isolation by distance (see Figure 4b). According to our results, the most probable scenario of migration of G. cruentatus between the studied localities was also the one allowing divergence and the migration from the ancestral population after the split. The divergence and migration routes are schematically depicted in the Figure 5b, while the estimates of immigration rates, divergence times and population sizes are listed in the Table S3. The results were not unequivocal, as can be seen from the Figure 5b, in several cases models proposing the opposite direction of divergence and migration between two populations had almost the same probability. This might be due to the low number of samples.

Discussion
Many previous phylogeographic studies have shown the existence of the geographical structure in populations of zoobiota within the north-eastern Atlantic and Mediterranean (i.e., Northern European Seas, Lusitanian and Mediterranean Sea provinces of Temperate Northern Atlantic realm sensu Spalding et al. [30]), e.g., [7,8,[10][11][12]. In this term, however, the European marine gobies, despite being the most speciose and abundant fish family in this area, have been little studied so far.
Our results on two epibenthic goby species (G. geniporus and G. cruentatus) showed that the most plausible model which can explain the genetic structure of populations of both species is a model of divergence and recurrent migration from the ancestral population after the split. In the case of G. geniporus, the direction of divergence and the migration routes match well the prevailing currents between the studied localities ( Figure 5a). The main feature of the migration route is anticlockwise circulation from Sicily towards east and then turning westwardly back to Sicily, and making a smaller circle from Sicily to the Adriatic Sea, following the Montenegrin and Croatian coast and subsequently the Italian coast, and back to Sicily. More dense sampling would be useful to confirm these findings, as also smaller gyres can substantially influence the genetic structure of epibenthic fishes [75]. The directions of divergence and migration between the pairs of populations of G. cruentatus were ambiguous. The observed pattern may be an outcome of a low number of individuals used to infer the migration routes in this species. Alternatively, it might be a consequence of higher complexity of water circulation in the species range, and/or biology of this species (see later discussion on hyperbenthic juveniles).
The lifestyle of the two studied species matches the model of divergence and recurrent migration. Being epibenthic and territorial, G. geniporus and G. cruentatus most probably exhibit only short-distance movements in adulthood, which allow the divergence between populations. Their main dispersal route is thus via a transport of planktonic larval stages, which can be dispersed by currents. The distance which a larva can reach depends mainly on the hydrodynamics and on the duration and behaviour of the larval stage, but the dispersal of planktonic larvae is much more complex and still not well understood [76]. The high multiscale variability of topography, temperature and salinity in the Mediterranean Sea generates free and boundary currents, bifurcating jets, meander and ring vortices, permanent or temporary cyclonic and anticyclonic gyres and eddies [77]. Recently, computer simulations that integrate a high number of biological and marine physical information have been successfully used in several works focused on the role of marine currents on the dispersion and genetic structure of marine organisms [75,[78][79][80].
The influence of currents on genetic structure of the populations of epibenthic marine fish species was found for Tripterygion tripteronotum [75], where the population structure matched well the gyres in the Adriatic Sea, and also for other marine organisms [78][79][80][81][82].
Where known, the planktonic life stages in different European goby species have a variable duration, with a minimum of 13 days in Zosterisessor ophiocephalus to 51 days in Gobius paganellus [83]. However, in many Mediterranean gobies, nothing is known about their larvae, the duration of this stage, nor about their dispersion routes or distances. A similar range of planktonic larval duration (PLD) was observed in other Mediterranean fish species. In epibenthic Mediterranean littoral fish species of the genus Tripterygion, the PLD is estimated to be two to three weeks [84], while blennies (Blenniidae) have a PLD between 22 and 71 days. In species of both these fish groups a population genetic subdivision was observed [75,[85][86][87].
The dispersal capability of fish larvae can broadly differ, while it is around 120 km during 80-170 days of PLD in Sebastes melanops [76,88], it is only 100-500 m during 30-50 days in Chaetodon vagabundus [76,89]. This underlines the complexity of the dispersion process of fish larval stages.
Apart from having mobile larval stages, G. cruentatus has hyperbenthic juveniles, swimming in shoals within 1 m above the sea bottom. It is not known which distances this stage can cover and whether the dispersion during this stage has any influence on the gene flow. In the aquarium, this stage lasted two months [90]. It is not known whether other European gobies have hyperbenthic juveniles. Also other biological traits, such as reproduction strategy (European benthic gobies are iteroparous [53]) and timing of spawning, can influence the genetic structure of populations [75].
Our results showed that there was a high diversity of haplotypes of cyt b at each sampled locality. As discussed above, no clear population subdivision was found in two studied species, as it was disturbed by the recurrent migrations between the populations. There was a certain structuring in both species, as two haplogroups are observable in the networks (Figure 3). In G. geniporus, in the most frequent haplotypes of each haplogroup, specimens from different areas dominate: in one, the specimens from the eastern (Cyprus and Greece), while in the other, the specimens from the central Mediterranean Sea (Italy, Montenegro and Croatia). However, haplotypes of specimens from the Sicilian population are prevailing in the haplogroup with the eastern Mediterranean Sea haplotypes. In G. cruentatus, one haplogroup includes all specimens from the western part of the species range, from the Atlantic coast of Spain and Portugal, as well as from the western Mediterranean French coast, while the other haplogroup includes all samples from the eastern Mediterranean Sea (Cyprus). However, the central Mediterranean samples (Sicily, Montenegro and Croatia) are present in both haplogroups. Similar situation, where haplotypes from different haplogroups were found at the same geographic locality, with no clear geographical pattern, was observed also in other fish species in the Mediterranean Sea [91,92]. It was attributed to the secondary contact between the isolated populations which diverged in allopatry and came to a contact again after the removal of the migration barrier [91]. Additionally, our migration scheme for G. cruentatus shows convergence of the routes from the eastern and western Mediterranean Sea and the Adriatic Sea near Sicily, corresponding to the situation in the haplotype network.
Most of the research on population genetic structure of marine gobies from Europe have been conducted on epibenthic species of the genus Pomatoschistus (gobionelline-like gobies [52]), usually inhabiting lagoons and shallow coastal waters with fine substrates. Population genetic differentiation was observed in all four studied Pomatoschistus species [44][45][46][47][48][49][50]93]. Population genetic diversity of species from the gobiine-like gobies [52] has been studied in only two European marine species [51], epibenthic goby G. niger, living on the muddy substrates, and Aphia minuta, a pelagic shoal species. Giovannotti et al. [51] found a spatial genetic structure in epibenthic G. niger, while no structure in the pelagic A. minuta.
There are several recognised biogeographic breaks in the Mediterranean Sea and the north-eastern Atlantic Ocean. Our data did not point to the existence of any biogeographic boundary preventing a gene flow between the studied populations for neither of the two species. However, the effect of a small sample size cannot be excluded. The Strait of Gibraltar, or rather the Almeria-Oran front, which is an important biogeographic barrier for some marine organisms [4,8,10,33,94], did not have any influence on the gene flow between Atlantic and western Mediterranean populations of G. cruentatus. Similarly, this break does not present a barrier to gene flow of the various fish species, neither pelagic, e.g., Sardina pilchardus (nDNA microsatellite loci) [95], Thunnus thynnus (mtDNA d-loop) [34], Scomber colias (mtDNA d-loop) [31], Diplodus sargus (mtDNA d-loop, nDNA S7 first intron) [96], nor benthic ones, ranging from widespread eurybathic Lophius piscatorius (mtDNA d-loop), able to reach depths down to 500 m [33] to Parablennius sanguinolentus (mtDNA d-loop, nDNA S7 first intron), which is restricted to very shallow littoral of 0-1 m depth [86,87]. Neither did the Sicily Channel influence the genetic structure of the two studied goby species, unlike is the case of some other fish species, e.g., Dicentrarchus labrax (nDNA microsatellite loci) [97], Sprattus sprattus (mtDNA d-loop) [11], and P. tortonesei (mtDNA 16S, COI) [48], where the Sicily Channel presents an important breakpoint. Although many studies showed a genetic differentiation between populations of the biota of the Adriatic and the Mediterranean Seas, separated by the Otranto strait, e.g., in P. minutus (mtDNA d-loop, cyt b, allozymes) [44][45][46]50], Platichthys flesus (allozymes) [98], Gouania willdenowi (mtDNA COI and 9 nDNA markers) [99] and Sparus aurata (allozymes) [100], neither was the Otranto Strait a biogeographic barrier for G. cruentatus and G. geniporus.

Conclusions
Our data revealed that the population genetic structure of the two studied epibenthic goby species (G. geniporus and G. cruentatus) can be well explained by the model of migration, allowing divergence between each pair of populations, with the ongoing migration from the ancestral population. This corresponds well with the biology of these gobies, having poorly mobile adults on one hand, and planktonic larval stages, which can be dispersed by currents, on the other hand. The population genetic structure of G. geniporus is influenced by currents: the estimated migration routes between the studied populations follow the main current directions in the study area.
Supplementary Materials: The following are available online at http://www.mdpi.com/2077-1312/8/10/792/s1, Figure S1: List of pairs of populations modelled in Migrate-n, Table S2: Diversity measures for Gobius geniporus and G. cruentatus calculated per each locality based on cytochrome b and S7 sequences, Table S3: Posterior distribution  table of