Can the Population Structure of Three Greek Marine Species ( Sardina pilchardus, Penaeus kerathurus, Mullus barbatus ) Become a Tool for their Future Characterization as PGI Products?

.


Introduction
In order to satisfy the increasing number of consumers demanding for edible products of certified origin, the European Commission established three protection labels: the "Protected Designation of Origin (PDO)" assigned to a product or a foodstuff produced, processed and prepared in a given area; the Protected Geographical Indication (PGI) that indicates a link of the product with the area where at least one of the stages of production, processing or preparation possesses a specific quality, reputation or other characteristics attributable to that geographical origin; the Traditional Specialities Guaranteed (TSG) for traditional products with specific characteristics Sardine (Sardina pilchardus, Walbaum 1792) is a small pelagic fish species of great interest to fisheries in central-eastern and north-eastern Atlantic Ocean. This cupleiform is found from the North Sea to Senegal, as well as in the Mediterranean Sea, Sea of Marmara, and Black Sea [1]. Northern and southern limits seem to be related to the average water temperature, being located within 10 and 20°C isotherm [2]. Nevertheless, several authors have hypothesized that sardine distribution and abundance are dependent on the hydrological regime [3]. Like other marine pelagic fishes, sardines show schooling and migratory behaviour, as well as great dispersal capabilities both at the larval and adult stages. Sardine supports important fisheries in the northeast Atlantic and in the Mediterranean Sea, with approximately 130,000 tons fished on the European coast, 660,000 tons fished on the African coast and 80,000 tons fished in the Mediterranean area [4][5][6]. In particular, Spain and Morocco are the countries with the largest captures (representing about the 77% of the total annual catch of sardines). In Greece it is one of the two most commercially important fish, with a total catch of 10,890.5t in 2016, representing 14.6% of the total catch of marine species [7].
The caramote prawn Penaeus kerathurus (Forskäl 1775) is an ecologically and economically important penaeid species. It is widespread in the Mediterranean and it also ranges from the south coast of England to Angola in the Eastern Atlantic. The benthic adults inhabit nearshore and offshore waters, to a depth of about 80 m, and prefer muddy or sandy-mud flats. In summer, adults migrate to reproduce in coastal areas and spawn in offshore waters. After a planktonic larval stage (about 4 weeks) post-larvae move into shallow waters, where they enter the juvenile stage until they reach 5-8 cm in length, and then join the adult population [8]. This species is extensively fished; the annual global capture is around 6,000 tons [5]. In Greece the total catch for 2016 was 1,404.1t, representing almost 2% of the total catch of marine species [7].
The red mullet Mullus barbatus Linnaeus, 1758 belongs to the Mullidae family and it is distributed in the eastern Atlantic, along the European and African coasts from the British Isles to Senegal, as well as in the Mediterranean Sea [9]. It is a gregarious benthic species living on sandy and muddy bottoms, usually between 10 and 100 m, occasionally in deeper waters up to 270 m [10]. It feeds on small benthic invertebrates such as crustaceans, worms, mollusks. The reproduction period extends from April to August at depths between 10 and 55 m on sandy or muddy bottoms and the postlarvae is pelagic until 28 mm SL [9]. Two subspecies are present, Mullus barbatus barbatus and Mullus barbatus ponticus which occur in the Black Sea and Azov Sea [9]. It is caught mainly by trawling fleet [11] and thus subject to intense fishing pressure; red mullet stocks are composed mainly of young fish [11]. The red mullet is one of the most commercial fish resources in the Mediterranean. In Greece it ranks among the most important demersal fish, with a total catch of 1,758.7t for 2016 (2.4% of the total catch of marine species) [7].
DNA-based methods like microsatellite genotyping [12,13], Internal Transcribed Spacer analysis (ITSs) [14,15], Single Nucleotide Polymorphisms (SNPs) [16] and High-Resolution Melting (HRM) analysis of mitochondrial markers [17] offer the possibility to identify individuals, breeds or populations along the food chain, in order to define them as PGI products. Once DNA is extracted from a biological sample (blood, muscle or even processed food), it is analyzed by molecular markers to assess variations among individuals or among populations [18].
The mitochondrial genome has a higher rate of mutation compared to the nuclear genome, is maternally inherited, has less hybridization and has a high copy number, which facilitates PCR amplification and sequence recovery from degraded tissue [19,20]. Furthermore, the mitochondrial genome lacks introns, pseudogenes and repetitive sequences, which makes sequence alignments of the amplified genes easier. Finally, complete mitochondrial DNA (mtDNA) genome sequences are available on public databases; primers can therefore be designed to amplify and sequence any species that has a published mtDNA genome [21-24].
Population genetic structure has been described for numerous species, using parts of the mitochondrial genome. A primary advantage of mtDNA is the inheritance pattern: clonal inheritance through the maternal line allows tracing of oceanographic events lasting multiple generations [25]. The . 03 .
Ocimum Scientific Publishers inheritance pattern of mtDNA involves that male dispersal does not participate in the homogenization of the population. Gene flow caused by male dispersal will not affect spatial patterns of variation in mtDNA [26], so, in general, mtDNA genes show more population structure than do nuclear genes [27]. This makes them particularly appropriate indicators of the population genetic differentiation of marine organisms, which are generally high gene-flow species.
Most population genetic studies of marine organisms have focused on coastal species, including sessile species with planktonic larvae. Significant population genetic structure at scales of hundreds of kilometres and smaller has been observed for a number of marine fish [28-33] and marine invertebrate species [34][35][36][37][38][39]. Two general principles emerge from many studies: first, that both marine fish and invertebrates are quite variable at the molecular level, and second, that this variability may determine genetically distinguishable, geographical populations in some of the species [40,25].
The aim of this research was to assess the potential of using the population differentiation of three Greek marine species (S. pilchardus, P. kerathurus, M. barbatus) as a tool for their future PGI definition. For this purpose, the PCR and further Sanger sequencing analysis of three mitochondrial segments (Cytochrome Oxidase subunit I -COI, cytochrome b -cytb, control region -D-loop) were applied. These mitochondrial markers have been repeatedly used for population differentiation analyses in all three species.  . 04 .

Submit Manuscript
-20°C until DNA preparation. In some cases, tissue samples were immediately extracted and preserved in absolute ethanol (Applichem Panreac, Germany).

DNA extraction
Total genomic DNA was extracted from muscle, with the UltraClean™ Tissue and Cells DNA Isolation Kit (MoBio, Germany), following the manufacturer's recommendations, and with the CTAB methodology [41]. The DNA concentration was determined by using the NanoDrop Spectrophotometer (ND-2000 Thermo Fisher Scientific, USA).

PCR amplification and sequencing analysis
Three mitochondrial segments were screened as potential markers for population differentiation analysis in this study: Cytochrome Oxidase subunit I (COI), cytochrome b (cytb) and the control region (D-loop). Double stranded DNA was amplified in a total reaction volume of 25μl containing 1 unit of Taq polymerase (Biolabs, USA), 2.5μl of 10 × reaction buffer (Biolabs, USA), 2.5 mM MgSO 4 (Biolabs, USA), 0.24 mM of each dNTP (Biolabs, USA), 1 pmole/μl of each primer (IDT, USA) and approximately 50-100 ng of DNA.
Electrophoresis of 3μl of the PCR product was performed in 1 × TBE buffer (Applichem Panreac, Germany) for 1h at 150V, in 1.5% agarose gel (Applichem Panreac, Germany) containing 0.5μg/ml Midori Green (Neppon Genetics, Germany). The size of the PCR products was checked against a 100bp DNA ladder (GeneOn, Germany). After the end of the electrophoresis, the resulting DNA fragments were visualized by UV transillumination and photographed.
PCR products were purified using MOBIO Ultra Clean PCR Clean-Up Kit (MoBio, Germany), following the manufacturer's instructions. In order to avoid false positive polymorphisms and to double check mutations, sequencing was carried out in both directions on a ΑΒΙ 3500 Genetic Analyzer (Applied Biosystems, Thermo Fisher Scientific, USA). Furthermore, for all the ambiguous sites the PCR and sequencing reactions were repeated. Both strands of each sequence were aligned with CLUSTAL W algorithm [42] as implemented in the software BioEdit v.7.0.5.3 [43] and visually confirmed. All sequences were confirmed as those of the corresponding genes by BLAST searches on GenBank except for P. kerathurus, the cytb and D-loop sequences of which were the first DNA sequences of the species submitted to any Genetic Database.

Statistical analysis
Different haplotypes for each mitochondrial segment were detected in all three populations of each species with the DAMBE6 software package [44]. Subsequently, haplotype frequencies were estimated for each population. In order to examine relationships among individuals, the maximum likelihood approach was applied, using MEGA7 software [45] and the complete nucleotide data set (COI, cytb, D-loop) of each analyzed sample. The best-fit substitution model was provided by the MEGA7 software [46]. Non-uniformity of evolutionary rates among sites were modelled by using a discrete Gamma distribution (+G) with 5 rate categories and by assuming that a certain fraction of sites are evolutionarily invariable (+I). All positions with less than 95% site coverage were eliminated. That is, fewer than 5% alignment gaps, missing data, and ambiguous bases were allowed at any position. Branch support was assessed by 1,000 replicates and sites with missing data were removed only when needed.
Genetic distances among the three populations per species were estimated based on both the nucleotide sequences of each segment and the total nucleotide sequences (COI, cytb, D-loop), using the 2-parameter Kimura distance model [47]. Phylogenetic relationships among them were estimated with

S. pilchardus
The size of the PCR products was approximately 700 base pairs (bp) for COI gene, 340 bp for cytb gene and 1600 bp for the D-loop region. In total, 641 bp at the 5΄ end of the COI gene, 294 bp at the 5΄ end of the cytb gene and 473 bp at the 5΄ end of the D-loop region were sequenced for the majority of the individuals.
Fifty-eight individuals were sequenced for COI gene; namely, 19 individuals from Amvrakikos Gulf population, 20 individuals from Ionian Sea and 19 individuals from Kalloni Bay population. Nine samples were replicated and the good quality sequences were successfully obtained with the second sequencing effort. Thirty-three haplotypes were detected for COI gene in the three S. pilchardus populations (Table 2SM). The most common Haplotype (H1) was found in 20 individuals of all three populations and two haplotypes were found in two populations. The high frequency and wide geographical distribution of haplotype 1 indicate that this is probably the ancestral haplotype. All the other haplotypes were detected in only one population.
Forty-nine individuals were sequenced for cytb gene; namely, 16 individuals from Amvrakikos Gulf population, 17 individuals from Ionian Sea and 16 individuals from Kalloni Bay population. Twelve samples had to be replicated. Twentyfive haplotypes were detected for cytb gene and the most common Haplotype (H1) was found in 23 individuals of all three populations (ancestral haplotype) (Table 3SM). All the other haplotypes were detected in only one population.
Finally, 55 individuals were sequenced for D-loop region; namely, 16 individuals from Amvrakikos Gulf population, 20 individuals from Ionian Sea and 19 individuals from Kalloni Bay population. Only five samples were repeated and the good quality sequences were successfully obtained with the second sequencing effort. Fifty-three haplotypes were found for the D-loop region (Table 4SM). All these haplotypes were detected in only one population, with the exception of H2, which was detected in two out of three populations. All haplotypes were deposited to GeneBank under the accession numbers MH141137-MH141169 (COI haplotypes), MH127862-MH127879 (cytb haplotypes) and MH141170-MH141222 (D-loop haplotypes).
The complete nucleotide dataset included 46 sequences with a total length of 1408 bp. The maximum likelihood topology estimated with the complete nucleotide dataset of each individual ( Figure 2) revealed no clustering of the samples in relation to the population they are coming from. Besides, most of the bootstrap values were lower than 40%.
The Neighbour Joining phylogenetic trees, constructed with the nucleotide sequences of each mitochondrial segment and with the complete nucleotide dataset (Figure 3), revealed the same structure: populations from Amvrakikos Gulf and Kalloni Bay were grouped together, whereas the population from Ionian Sea formed a separate clade. The most significant outcome of the phylogenetic analysis was that both the UPGMA and the Minimum Evolution trees (based on the nucleotide sequences of each segment and on the total sequences) had identical topology: the first cluster included populations from Amvrakikos Gulf and Kalloni Bay while the second one was formed by the Ionian Sea population.
For population differentiation estimates, the total nucleotide sequences (COI, cytb, D-loop) for each individual were used.
The AMOVA analysis showed that most of the genetic variation was present within samples (98.77%) and 1.23% of genetic heterogeneity was apportioned among them ( Table 1). The total F st value was F st = 0.01232. F st and statistical P values among all the studied population samples are shown in Table  2. The highest F st value occurred between the populations from Ionian Sea and Amvrakikos Gulf. It is also obvious that the Ionian Sea population was clearly differentiated from the Amvrakikos Gulf population (P < 0.05), while there was no significant differentiation among the other pairs of populations.

P. kerathurus
The size of the PCR products was approximately 700 base pairs (bp) for COI gene, 350 bp for cytb gene and 600 bp for the D-loop region. In total, 632 bp at the 5΄ end of the COI gene, 316 bp at the 5΄ end of the cytb gene and 552 bp at the 5΄ end of the D-loop region were sequenced for the majority of the individuals.
Fifty-nine individuals were sequenced for COI gene; namely, 20 individuals from Amvrakikos Gulf population, 20 individuals from Ionian Sea and 19 individuals from Thracian Sea population; only two samples had to be replicated. Thirty haplotypes were detected for COI gene in all three populations of P. kerathurus, and only Haplotype 1 (H1) was found in 14 individuals of both Amvrakikos Gulf and Ionian Sea populations (Table 5SM). All the other haplotypes were detected in only one population.
Fifty-seven individuals were sequenced for cytb gene; namely, 20 individuals from Amvrakikos Gulf population, 19 individuals from Ionian Sea and 18 individuals from Thracian Sea population. Overall, 15 samples were replicated and the good quality sequences were successfully obtained with the second sequencing effort. Forty-one haplotypes were determined for cytb gene while the most common Haplotype (H1) was found in 12 individuals of all three populations (ancestral haplotype) ( Table 6SM). All the other haplotypes were detected in only one population.
Fifty-one individuals were sequenced for D-loop region; namely, 18 individuals from Amvrakikos Gulf population, 19 individuals from Ionian Sea and 14 individuals from Thracian Sea population. Overall, five samples were re-sequenced and fifty-one haplotypes were found for the D-loop region, none of which was common among the populations: 18 haplotypes appeared in Amvrakikos Gulf population, 19 haplotypes were detected in Ionian Sea population and 14 haplotypes were found in the population from Thracian Sea (Table 7SM). All haplotypes were deposited to Gene Bank under the accession numbers MF939109-MF939138 (COI The complete nucleotide dataset included 49 sequences with a total length of 1500 bp. The maximum likelihood topology estimated with the complete nucleotide dataset of each sample ( Figure 4) showed that most of the individuals from Thracian Sea population were grouped separately, with high bootstrap values (66%, 96% and 98%). The Neighbour Joining phylogenetic trees, constructed with the nucleotide sequences of each mitochondrial segment and with the complete nucleotide dataset ( Figure 5), revealed the same structure: populations from Amvrakikos Gulf and Ionian Sea were grouped together whereas the population from Thracian Sea formed a separate clade. The most remarkable outcome of the phylogenetic analysis was that both the UPGMA and the Minimum Evolution trees (based on the nucleotide sequences of each segment and on the total sequences) had the same structure: the first cluster included populations from Amvrakikos Gulf and Ionian Sea while the second one consisted of the Thracian Sea population.  For population differentiation estimates total nucleotide sequences (COI, cytb, D-loop) for each individual were used.
The AMOVA analysis showed that 87.90% of the genetic variation was present within samples and 12.10% of genetic heterogeneity was found among them (

M. barbatus
The size of the PCR products was approximately 700 base pairs (bp) for COI gene and 620 bp for cytb gene. For the amplification of the D-loop region, previous universal primers were tested [52]. The resulted PCR product had many bands but there was no band with the expected size. Therefore, the mitochondrial D-loop segment was not used for further analysis. In total, 637 bp at the 5΄ end of the COI gene and 576 bp at the 5΄ end of the cytb gene were sequenced for all the individuals.
Fifty-five individuals were sequenced for COI gene; namely, 20 individuals from Amvrakikos Gulf population, 21 individuals from Ionian Sea and 14 individuals from Kalloni Bay population. Nine samples were re-sequenced and thirty-two haplotypes were detected for COI gene in the three M. barbatus populations (Table 8SM). The most common Haplotype (H1) was found in 14 individuals of all three populations. The high frequency and wide geographical distribution of haplotype 1 indicate that this is probably the ancestral haplotype. Two other Haplotypes (H3 and H4) were common for individuals of the three populations but in lower frequencies. Two haplotypes were detected in two out of three populations and all the other haplotypes were found in only one population (Table 8SM).
Altogether, fifty-five individuals were sequenced for cytb gene; namely, 20 individuals from Amvrakikos Gulf population, 21 individuals from Ionian Sea and 14 individuals from Kalloni Bay population. Overall, ten samples were replicated and the good quality sequences were successfully obtained with the second sequencing effort. Nineteen haplotypes were found for cytb gene (Table 9SM). Three haplotypes were common among the three populations and two haplotypes were detected Submit Manuscript . 11 .
in two out of three populations. All the other haplotypes appeared in only one population. All haplotypes were deposited to Gene Bank under the accession numbers MG969169-MG969200 (COI haplotypes) and MG969201-MG969219 (cytb haplotypes).
The complete nucleotide dataset included 55 sequences with a total length of 1213 bp. The maximum likelihood topology estimated with the complete nucleotide dataset of each individual ( Figure 6) revealed no clustering of the samples in relation to the population they are coming from. In addition, most of the bootstrap values were lower than 40%.  For population differentiation estimates total nucleotide sequences (COI, cytb) for each individual were used. The AMOVA analysis showed that most of the genetic variation was present within samples (96.85%) and 3.15% of genetic heterogeneity was apportioned among them (

Discussion
Populations constitute interbreeding units with more or less autonomous dynamics. Population genetic structure is determined by the level of connectivity or exchange between individuals, and dispersal potential is a key factor in shaping the structure of populations [53]. In terrestrial and freshwater environments, populations are often well defined and distinct from each other, often physically separated by barriers to mixing and interbreeding [54]. Fish populations can be described from a biological perspective, which implies some level of reproductive isolation, or from a fishery perspective, which concerns a practical description of a group of fish exploited in a specific area [55].
The genetic stocks are intraspecific groups of randomly mating individuals with temporal and spatial integrity that are at some extent reproductively isolated from each other, and thus will react independently to exploitation [56]. Within marine populations it can be difficult to apply the concept of structural subdivisions, as many marine ecosystems lack obvious barriers [57]. The continuous water represents a potential means for dispersal, favouring intermixing of individuals over the species range. In practice, it is not easy to identify stocks, as the delimitation of adjacent populations involves many issues, especially in the sea where there are no clear geographical barriers. Tides and ocean currents may further act to mix passively drifting organisms, primarily eggs and larvae, over appreciable distances [58]. In marine species it is generally assumed that a high capacity for dispersal in the early life history stages results in reduced intra-specific differentiation over smaller spatial scales [59] and less genetic structure. Moreover, the effect of genetic drift in promoting differentiation is diminishing due to the very large population sizes [60]. In fact, marine populations often consist of localised sub-populations that are relatively independent and have distinct ecological and genetic properties [61].
Despite the general observation of relatively little geographic variation among populations of high dispersal marine species, many cases of divergence have been observed. Several evolutionary forces affect the amount and distribution of genetic variation among marine populations and thereby population differentiation [62]. Geographic distance and physical barriers enhance reproductive isolation by limiting the migration, and increase genetic differentiation between populations [63]. In cases where divergence among populations has been observed in an otherwise seemingly highly dispersal species, natural selection has usually been invoked to explain the pattern of genetic variation. The rationale for this argument is straightforward: when larval dispersal is large, there is virtually no capacity for populations to diverge by random processes, so significant genetic divergence among populations must be driven by selection [64]. The use of genetic markers for population description requires a detectable level of genetic differentiation, and this has presented problems in studies of many marine organisms [65]. In the marine environment many studies have failed to detect statistically significant population structuring because of low differentiation, especially over small geographical distances [66].
The different haplotypes inside a population represent the various nucleotide sequences having its individuals. Genetic distance among individuals or among populations is a measure of the genetic divergence between them [67], and arises primarily from the separate nucleotide haplotypes. The genetic differentiation of populations is the result of uneven spatial distribution of genetic variation in a species. All these genetic parameters constitute the population structure which is mainly based on the individuality of specimens, an item that is necessary for the PGI designation. For these reasons the population structure of a species could be used as a reliable marker for its future PGI definition.

S. pilchardus
Morphological studies based on gill raker counts and head length [1] reported enough phenotypic variation to differentiate two subspecies, specifically the S. p. pilchardus (Eastern Atlantic Ocimum Scientific Publishers . 14 .

Submit Manuscript
Journal of Nutrition, Food and Lipid Science 10.33513/NFLS/1901-08 Ocean, from the North Sea to Southern Portugal) and the S. p. sardina (Mediterranean Sea and Northwest African coast) [1]. Although no private mitochondrial control region sequence haplotype could be found for each proposed subspecies, they were suggested to be genetically distinct based on significant pairwise haplotype frequency differences [68].
There have been many studies addressing the genetic structure of sardine in different parts of its global distribution, using different molecular markers such as allozymes [69][70][71][72][73][74], mtDNA [68,[75][76][77][78], microsatellites [78][79][80] and intron polymorphism [81]. The results of these studies, though not completely congruent, suggest a very weak genetic structure and low levels of genetic differentiation. In populations of pelagic fish, little genetic structure is expected, because adults often disperse over large distances. Populations of pelagic fish are therefore unfragmented over large areas, because groups of fish can migrate thousands of kilometres [73], usually leading to weak genetic differentiation [65].
In the present study the common haplotypes revealed in all three sardine populations for COI and cytb segments, resulted in similar findings with low levels of genetic structuring among populations. Furthermore, most of the genetic variation detected within samples and the total F st value was F st = 0.01232, which means that only 1.23% of genetic heterogeneity was apportioned among populations. The mean F st revealed from the present work has a higher/lower value compared with previous findings (Table 10SM).
According to a previous study [69], the genetic distances based on electrophoretic variation and multivariate analysis of several morphometric and meristic characters, suggest that populations of sardine from the Aegean and Ionian seas do not form one panmictic population. The comparison of the two basins is highly significant, accounting for 90% of the total variance. These findings do not support the idea of the pure or discrete stock concept, and suggest a fair amount of gene exchange between seas, but not enough to homogenize the populations completely. Instead, they provide evidence for the dynamic population structure model according to which physical (e.g., hydrographic) or biological (e.g., predation or behaviour) factors impose a population structure. The existence of well-defined geographic and hydrographic boundaries in the Greek peninsula and related hydrographic conditions, can partially structure populations of this pelagic fish inhabiting adjacent areas [69].
Our findings do not seem to support completely the hypothesis of the population genetic structuring among the two basins, as there is no genetic differentiation between the Ionian Sea and the Kalloni Bay populations (F st = -0.00126; P = 0.45045).
Nevertheless, the Ionian Sea population formed a separate cluster in all four dendrograms. In addition, the highest F st value was revealed between the populations from Ionian Sea and Amvrakikos Gulf (F st = 0.03119) and it was significant (P = 0.01802). A value of F st > 0.02 for marine fish means that there is a genetic population structure [65]. These findings imply a subtle but significant genetic structuring of the Ionian Sea sample. Such a pattern could be explained by the passive dispersion of larvae with the marine currents, rather than active migration of adult individuals [28], or by recent colonization events. Furthermore, the weak but significant differentiation between the two populations may be the result of genetic drift, which is likely maintained by the deep waters of the Ionian Sea (the Ionian Sea is the deepest sea of the Mediterranean basin) that separate this population from the shallow (maximum depth: 60 m), semi-enclosed and protected Amvrakikos Gulf, with different hydro ecological characteristics [82]. These findings also support the hypothesis that the physicochemical characteristics (temperature, salinity, water column) [82,83] and oceanographic conditions (current speeds) [84] prevailing in the Amvrakikos Gulf may reduce gene flow between locations within and outside the Gulf. Further mtDNA or microsatellite analysis of sardine samples from more locations in the Ionian Sea could facilitate the discrimination of these populations, and give more data for a possible characterization of this stock as a PGI product.
The apparent lack of a spatial pattern to the genetic differentiation observed in the present study, especially the lack of a correlation between differentiation and geographical distance, indicates a mechanism of gene flow that is independent of distance. The observed structuring was unlikely to be caused by isolation due to geographical distance. Passive transport of egg or larvae seems the most likely explanation, affected by a combination of various processes, such as hydrodynamics, geographical structure, and environmental conditions. Furthermore, any differentiation at such a restricted geographical scale may be interpreted as an innate tendency for population substructuring and not just a consequence of external forces.

P. kerathurus
Multiple studies of genetic diversity in natural and hatchery populations of several penaeid species have been conducted, most of which have been of allozyme variation [35]. Many authors suggested population genetic structuring over large or short-geographic distances in several penaeid species [85][86][87][88] and this may affect the management of the species. Levels of worldwide genetic population structure vary widely in different penaeid species, so that genetic structure can be strong over very short distances and weaker across large geographic scales [35, [88][89][90][91]. Despite its economic and ecological importance,  [92][93][94][95] but there are also few mtDNA studies [39,96,97].
Our results revealed a significant population differentiation in the total data set and a strong genetic structure among all the pairs of populations. The total F st value found in the present study indicated that about 12% of the total gene diversity observed was due to population differentiation, and almost 88% was due to variation among individuals within populations. The mean F st revealed from the present work has a higher/ similar value compared with previous findings (Table 11SM).
Unique haplotypes were found in COI gene and D-loop region for individuals from Thracian Sea population. Additionally, as for cytb gene, the common H1 haplotype was found in just a single individual from Thracian Sea population, and all the other haplotypes were unique for this population. As a result, the highest F st values were revealed for prawn population from Thracian Sea, while populations from Amvrakikos Gulf and Ionian Sea were less differentiated. The study of the genetic variability of nine Mediterranean and two Atlantic samples of caramote prawn with the mtDNA COI region revealed similar significant genetic structure between populations from Kavala (Thracian Sea) and Amvrakikos Gulf (with a higher value of F st = 0.532) [39]. Generally, geographic patterns of genetic variation in invertebrates species suggest that isolation by distance occurs, but only over large geographic scales [40]. In the present study the distance between the Thracian Sea population and the Amvrakikos Gulf population was approximately 1100 km, whereas the estimated distance between the Thracian Sea population and the Ionian Sea population was 1200 km.
The genetic structure observed here for P. kerathurus populations reflects the two distinct biogeographical regions of Aegean and Ionian Sea. The Mediterranean has sometimes been referred as a "sea of seas" because of its division into different sub basins, each having its own distinct characteristics, which includes partially enclosed current systems and likely different ecological conditions. Also, the Aegean Sea has a complex archipelago, highly irregular coastline, and it is in fact a combination of semi-isolated deep basins [98]. Thus, distinct hydrographic and ecological conditions may be sufficient to reduce gene flow [31]. It is possible that the three localities sampled in the present study belong to distinct reproductive units (stocks) that have apparently been genetically isolated from one another. Natural selection through local adaptation, vicariance events, and/or current dispersal could explain this pattern of genetic differentiation.
Furthermore, this pattern of genetic structure may be explained by hydrographic incidents during the Pleistocene. Throughout this period the sea level has dropped and modified coast lines, splitting apart the eastern and western basins. During this period, climate fluctuations produced episodes of habitat fragmentation and promoted genetic discontinuities across geographical ranges. However, other contemporary factors that limit effective genetic dispersal, including oceanographic currents (in benthic species with long pelagic larval stages like P. kerathurus, water currents are assumed to play an important role in shaping the structure of genetic polymorphism) and larval behaviour may also be of critical importance [99,100]. However, the predominant mechanisms leading to population differentiation are not always clear [40].
In our case, this model of heterogeneity is probably linked to the particular biological cycle of P. kerathurus. This cycle strictly depends on low water salinity, thus limiting the geographical dispersion of the species to a few favourable places, where the larva-adult-larva cycle can be completed [101]. Hence, the distribution of the reproductive populations of this species can be very patchy. Genetic changes are thought to alter substantially the genetic architecture of such populations, allowing rapid accumulation of many genetic differences that can lead to reproductive isolation. The well-known genetic processes of mutation and selection may be the most powerful forces creating reproductive isolation [40].
The P. kerathurus population from Thracian Sea (Northern Aegean Sea) revealed a significant genetic differentiation and a higher value of F st , compared with the other samples. Based on our finding, it could be suggested that this population constitute a different genetic stock and a possible PGI product. Future studies could reinforce the hypothesis of its characterization as a PGI product.
In the present study three common haplotypes among all red mullet populations were found for both COI and cytb genes. This finding is an evidence of important gene flow between collecting sites. There was no clear discrimination of Amvrakikos Gulf population from Ionian Sea and Kalloni Bay populations in the NJ tree. Only population from Kalloni Bay was significantly differentiated from the Ionian Sea population (P < 0.05), while there was no significant differentiation among the other pairs of populations.

Submit Manuscript
Low levels of differentiation in marine organisms are most likely due to extensive gene flow [54,57] and exchange of pelagic forms is assumed to be the major mechanism uniting spatially discrete populations. Theoretically, gene flow of few individuals per generation would be sufficient to prevent the accumulation of significant genetic drift between locations.
On the other hand, subtle population genetic differentiation does not necessarily imply that structuring does not exist, but that more powerful tools are required to detect it [57]. Marine organisms, even weakly differentiated on a small geographical scale, often show evidence of differentiation over larger distances, probably because the long distance acts as an isolation mechanism. Indeed, because marine species can often experience very high gene flow, the differentiation of populations can be detected on very large spatial and longtime scales [64].
Furthermore, M. barbatus is a demersal species, displaying a close relationship between habitat characteristics (depth, salinity, temperature, etc.) and life history. Larvae and juveniles are pelagic and when they move they follow currents at the surface. On the contrary, the adults inhabit deeper water (10-300 m) and they probably move over large areas covering long distances. Consequently, the distribution pattern of red mullet involves continuous admixture.
The genetic polymorphism of Greek M. barbatus populations has been studied with different techniques: allozymes [104,120], RAPDs [108,120] and mtDNA-RFLPs [111,120]. In both allozymes [104] and RARDs [108] analyses of Greek red mullet samples the populations from Ionian Sea formed a separate cluster in the UPGMA dendrogram and were differentiated from the Aegean Sea populations. On the contrary, samples of M. barbatus showed no genetic structuring with mtDNA-RFLPs analysis, and Corfu/Amvrakikos populations were clustered separately with no obvious differentiation [111]. A non-significant differentiation between the populations from Amvrakikos Gulf and Ionian Sea and zero differentiation between Amvrakikos Gulf and Kalloni Bay populations were revealed in the present work. The assumption of high gene flow has been discussed above. In addition, the only pelagic stage that the species goes through is the post larval stage, during which it exhibits intentional movement. Consequently, the distribution pattern of M. barbatus must be less random and more fixed, directed toward species biotopes.
In a previous study [108] the Aegean Sea samples seemed to diverge from the Ionian Sea samples and it was concluded that gene flow levels were insufficient to homogenise completely the red mullet populations. In the present study a significant genetic differentiation among Ionian Sea and Kalloni Bay populations (F st = 0.04768, P < 0.05) was revealed. These results strengthen the hypothesis that there is a differentiation between the two basins (i.e., Ionian basin, Aegean basin). Bathymetric constraints may act as physical barriers, limiting migration of red mullet adults between Aegean Sea and Ionian Sea. Mass transportation of pelagic eggs and larvae between the two basins may also be limited by oceanographic conditions [108]. Both of these assumptions could support an explanation in our case. A similar significant differentiation between an Ionian Sea and an Aegean Sea population was revealed (F st = 0.021, P < 0.01) from microsatellite data [119]. Additionally, red mullet populations from Ionian Sea were found to be significantly differentiated from other eastern and western Mediterranean Sea populations (F st values ranged from 0.021 to 0.082) [115], or to have subtle but significant differences from Adriatic and western Mediterranean Sea populations [116]. The mean F st found in the present study (0.03146) has a higher/similar/lower value compared with previous findings (  [116] 0.005 Adriatic Sea microsatellites [113] 0.007 Western Mediterranean Sea microsatellites [118] 0.009 Mediterranean Sea allozymes [105] 0.027 Adriatic Sea, Mediterranean Sea microsatellites [119] 0.043 Mediterranean Sea allozymes [104] 0.047 Gulf of Pagasitikos allozymes [120] 0.057 Gulf of Pagasitikos RAPDs [120] 0.5 Mediterranean Sea mtDNA-RFLPs [111] 0.5 Gulf of Pagasitikos mtDNA-RFLPs [120]  Data from the present work as well as from previous studies [108,115,116,119] indicates a significant genetic heterogeneity of the Ionian Sea red mullet stock. These preliminary findings if combined with further genetic analyses from more sampling sites in the Ionian basin, could contribute to a future definition of this stock as a PGI product.
The significant population differentiation and the strong genetic structure among all pairs of P. kerathurus populations are probably connected with the particular biological cycle of the species: the nursery grounds are located in shallow waters near river estuaries areas and the benthic adults inhabit nearshore and offshore waters. These conditions limit the geographical dispersion of the species to a few favourable places. Mutations modify substantially the genetic architecture of the reproductive populations, allowing rapid accumulation of genetic differences that can lead to reproductive isolation. On the other hand, the weak genetic structure and low levels of genetic differentiation found for S. pilchardus and M. barbatus populations could be justified by their wide dispersal range. Specifically, S. pilchardus is a pelagic species with pelagic eggs and M. barbatus life cycle includes several pelagic phases involving eggs, larvae, postlarvae and juveniles. As a consequence, populations of the two species can widen their dispersal range to a scale of hundreds of kilometres, usually leading to weak genetic structure.

Conclusion
In this work the PCR and Sanger sequencing analysis of three mitochondrial segments were used to investigate the possibility of characterizing certain populations of three