Genetic diversity and morphological characterisation of three turbot (Scophthalmus maximus L., 1758) populations along the Bulgarian Black Sea coast

Turbot (Scophthalmus maximus L., 1758) is a valuable commercial fish species classified as endangered. The conservation and sustainability of the turbot populations require knowledge of the population’s genetic structure and constant monitoring of its biodiversity. The present study was performed to evaluate the population structure of turbot along the Bulgarian Black Sea coast using seven pairs of microsatellites, two mitochondrial DNA (COIII and CR) and 23 morphological (15 morphometric and 8 meristic) markers. A total of 72 specimens at three locations were genotyped and 59 alleles were identified. The observed number of alleles of microsatellites was more than the effective number of alleles. The overall mean values of observed (Ho) and expected heterogeneity (He) were 0.638 and 0.685. A high rate of migration between turbot populations (overall mean of Nm = 17.484), with the maximum value (19.498) between Shabla and Nesebar locations, was observed. This result corresponded to the low level of genetic differentiation amongst these populations (overall mean Fst = 0.014), but there was no correlation between genetic and geographical distance. A high level of genetic diversity in the populations was also observed. The average Garza-Williamson M index value for all populations was low (0.359), suggesting a reduction in genetic variation due to a founder effect or a genetic bottleneck. Concerning mitochondrial DNA, a Nature Conservation 43: 123–146 (2021) doi: 10.3897/natureconservation.43.64195 https://natureconservation.pensoft.net Copyright Petya Ivanova et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. RESEARCH ARTICLE Launched to accelerate biodiversity conservation A peer-reviewed open-access journal Petya Ivanova et al. / Nature Conservation 43: 123–146 (2021) 124 total number of 17 haplotypes for COIII and 41 haplotypes for CR were identified. The mitochondrial DNA control region showed patterns with high haplotype diversity and very low nucleotide diversity, indicating a significant number of closely-related haplotypes and suggesting that this population may have undergone a recent expansion. Tajima’s D test and Fu’s FS test suggested recent population growth. Pairwise Fst values were very low. The admixture and lack of genetic structuring found pointed to the populations analysed probably belonging to the same genetic unit. Therefore, a proper understanding and a sound knowledge of the level and distribution of genetic diversity in turbot is an important prerequisite for successful sustainable development and conservation strategies to preserve their evolutionary potential.


Introduction
Over the past few decades, human impacts on wild fish populations have increased drastically worldwide as a result of extensive aquaculture, exploitation of fish stocks for global consumption and human-induced climate changes (Olsson et al. 2007). Excessive fishing and large-scale environmental changes not only affect the spatial distribution and structure of populations, but also introduce extensive modifications into their genetic diversity (Madduppa et al. 2018). The genetic diversity is a fundamental estimation for a population genetic study and crucial to the fishery management and resource utilisation (Xu et al. 2019). Marine fishes generally show low levels of genetic differentiation amongst geographic regions because of higher dispersal potential during planktonic egg, larval or adult history stages, coupled with an absence of physical barriers to the movement between ocean basins or adjacent continental margins (Grant and Bowen 1998;Hewitt 2000). Loss of genetic diversity can reduce the adaptability, lessen the population persistence and lead to a decrease in the productivity of the target species. Such genetic diversity reductions in some of the world's most abundant species may add to the growing long-term impact of fishing on their evolutionary potential, particularly with the abundance staying low and diversity continuing to erode (Pinsky and Palumbi 2014).
The turbot (Scophthalmus maximus) is a marine flatfish, with a high commercial value living on the European continental shelf and drawing remarkable attention with respect to fisheries and aquaculture (Iyengar et al. 2000). Despite its economic importance, the turbot is considered vulnerable under the current IUCN Red List Criteria (IUCN 2019). The wild populations of turbot are exposed to a strong anthropogenic pressure as it is considered one of the most valuable commercial species subjected to intensive fishing and is characterised as exploited unsustainably and at a risk of extinction in the Black Sea (Nikolov et al. 2015).
Information on the genetic structure of commercially important fish species is crucial to prevent ecological damage and to ensure sustainable and effective management of exploited stocks and systems (Liu and Cordes 2004;Olsson et al. 2007). Molecular markers (microsatellites and mtDNA) were applied only to closely-related turbot species in different marine regions in order to assess the genetic diversity (Bouza et al. 2002;Suzuki et al. 2004;Vandamme et al. 2014Vandamme et al. , 2020do Prado et al. 2018). Thus, using these molecular markers, limited data on the population structure of S. maximus in the Black Sea were obtained (Atanassov et al. 2011;Nikolov et al. 2015;Bessonova and Nebesikhina 2019;Turan et al. 2019b;Firidin et al. 2020;Ivanova et al. 2020).
Accordingly, to avoid depletion of the genetic diversity, fisheries management should be based on a more comprehensive knowledge of the genetic integrity of the populations. The modern molecular methods developed over the past few years offer unique opportunities to rate the genetic population structures; moreover, the subsequent evaluations should be further used to smooth the process of local management and to promote increased harvest under a sustainable fisheries regime.
The aim of the present study is to evaluate the population genetic diversity of three turbot populations along the Bulgarian Black Sea coast and its applicability for the purposes of monitoring and conservation of genetic diversity in terms of sustainable management and rational exploitation of stocks.

Sample collection and DNA extraction
Seventy-two turbot samples were caught by a local fishing vessel in March 2019 and April 2020 at two locations in the Black Sea (Shabla, Shkorpilovtsi and Nesebar) (Fig. 1). For DNA analysis, tissue samples were taken from dorsal fin and stored in 96% ethanol at 4 °C prior to the following analyses. The genomic DNA was extracted by DNeasy Blood & Tissue Kit (QIAGEN). All the DNA extracts were analysed by gel electrophoresis to monitor the DNA quality prior the PCR amplification.

Microsatellite genotyping
Seven microsatellite loci (Sma1-125INRA, Smax-02, Sma3-12INRA, 3/9CA15, B12-IGT14, Sma-E79 and Sma-USC26) were amplified and analysed (Table 1). The polymerase chain reaction (PCR) was performed in a reaction volume of 50 µl containing 2 µl of each primer, 25 µl of Mastermix (MyTaqTM HS Mix) and 2 µl of the target DNA. The reverse microsatellite primers were fluorescently labelled at the 5' end with the 6-FAM dye ( Table 1). The PCR amplification was performed under the following conditions: 35 cycles (95 °C for 1 min, 95 °C for 45 sec, 56-60 °C for 50 sec depending on the primer type, 72 °C for 1 min) and 72 °C for 10 min. Applied Biosystems 3130 Genetic Analyzer (Thermo Fisher Scientific) was used to carry out the fragment analysis. The size of the fragments was determined with GeneMapper 4.0 (Thermo Fisher Scientific) software package.

PCR and sequence analysis of mitochondrial DNA (COIII)
The polymerase chain reaction (PCR) using mitochondrial primers (COIII) was carried out in a reaction volume of 50 µl containing 2 µl of each primer, 25 µl of the Mastermix (MyTaqTM HS Mix) and 2 µl of the target DNA. The mitochondrial cytochrome c oxidase subunit III (COIII) was amplified using universal primers (F: AGCCCATGACCTTTAACAGG and R: GACTACATCAACAAAATGTCAGTAT-CA, according to Valles-Jiménez (2005). The conditions for PCR amplification included the following parameters: 94 °C for 5 min, 95 °C for 30 sec, 49 °C for 1min, 72 °C for 1 min (35 cycles) and 72 °C for 5 min. The mtDNA sequencing was performed by Macrogen Europe B.V. The obtained COIII sequences were deposited in the GenBank under accession numbers MN556886-MN556913 and MW446249-MW446288.

PCR and sequence analysis of mitochondrial DNA (CR)
The polymerase chain reaction (PCR) using mitochondrial primers (CR) was carried out in a reaction volume of 50 µl, containing 2 µl of each primer, 25 µl of Mastermix (MyTaqTM HS Mix) and 2 µl of target DNA. The mitochondrial control region (CR) was amplified using universal primers (L15924: 5'AGCTCAGCGCCAGAGCGC-CGGTCTTGTAAA and H16498-5'-CCTGAAGTAGGAACCAGATG, according to Atanassov et al. (2011). The conditions of PCR amplification included the following parameters: 94 °C for 5 min, 95 °C for 45 sec, 50 °C for 1 min, 72 °C for 1 min (35 cycles) and 72 °C for 5 min. PCR product quality control was performed by electrophoresis on 2% agarose gel. The sequencing was performed by Macrogen Europe B.V. The obtained CR sequences were deposited in the GenBank under the accession numbers MN556856-MN556885 and MW446289-MW446330. The new sequences were compiled and analysed with previously deposited sequences extracted from the GenBank database covering all available CR sequences of different populations of S. maximus along the Bulgarian Black Sea coast.

Statistical analyses
The Hardy-Weinberg equilibrium (HWE) exact tests and loci combinations for linkage disequilibrium with the Markov Chain methods were conducted using GenAlEx 6.5 (Peakall and Smouse 2012). The indices for genetic diversity as Shannon's information index, observed, expected unbiased, expected heterozygosities and inbreeding coefficient were calculated, against that background also using GenAlEx 6.5 (Peakall and Smouse 2012). Polymorphic information content (PIC) and Garza and Williamson index were computed using Arlequin v.3.5.2.2 (Excoffier and Lischer 2010). The analyses of private alleles (alleles observed in only one population), the allelic richness and Mantel test were accomplished by the programme GenAlEx 6.5 (Peakall and Smouse 2012).
The inbreeding coefficient (Fis) was thereupon calculated using Genepop 4.7 (Rousset 2008) and mtDNA sequence data were further analysed by applying MEGA7 (Kumar et al. 2018). The sequence alignment was established, the number of haplotypes, haplotype network and molecular phylogenetic neighbour joining (NJ) tree (Saitou and Nei 1987) and TCS Networks were easily constructed by means of PopArt (Clement et al. 2002). The bootstrap consensus tree was inferred from 1000 replicates. The evolutionary distances were computed using the Tajima-Nei method (Tajima and Nei 1984) and the haplotype (H) and nucleotide (π) diversities were estimated with DnaSP 5.10.01 (Librado and Rozas 2009). The fixation index Fst calculated in Arlequin v. 3.5.2.2 (Excoffier and Lischer 2010) was used to infer the genetic structure. The significance test was obtained from 1023 permutations. Bayesian clustering analysis was performed in STRUCTURE v.2.3.4 (Pritchard et al. 2000) with 10 independent iterations at each K = 1-10 with 20 5 Markov Chain Monte Carlo iterations after a burning period of 20 5 repetitions. The optimum K was estimated using the Structure Harvester (Earl and Von Holdt 2012) following the Evanno method (Evanno et al. 2005) for BIC comparison. Structure graphical results were plotted with CLUMPAK (Kopelman et al. 2015).Correlation analysis was appropriately implemented to identify the specific pattern of linear relationships between morphometric and meristic characters to be further analysed and validated amongst the studied populations, aiming at assessing the similarity/ dissimilarity in morphology. The non-parametric statistical test -similarity analysis (ANOSIM) was conducted in order to identify statistically significant differences between samples and hierarchical clustering and NDMS were used to visualise the distribution of samples with respect to their morphology similarities. All statistical analyses and graphic representations were performed using the statistical and programming software R 4.3.8 (R Core Team 2020), packages: 'PerformanceAnalytics' (Peterson and Carl 2014), 'vegan' (Oksanen et al. 2019), 'tidyverse' (Wickham 2017) for data manipulation, 'cluster' (Maechler et al. 2018) for selection and implementation of the proper clustering algorithm and 'factoextra' (Lê et al. 2008) for clustering visualisation, as well as 'dendextend' (Galili 2015), available through the CRAN repository (www.r-project.org).

Genetic variability of microsatellite loci
All analysed microsatellite markers proved to be polymorphic and a total of 59 alleles was discerned within the seven loci. The length of the identified alleles in the investigated loci ranged between 85 and 320 base pairs (bp). The number of alleles per locus ranged between 2 and 13 ( Table 2). As a result, Sma-USC26 (13 alleles) was the most polymorphic locus and Smax-02 (two alleles) was the least polymorphic locus. Analysis of allelic frequencies by loci showed that the lowest allelic frequency (0.017) was found at locus Sma-E79 and the highest one (0.717) at locus Smax-02. Reasonable amount of polymorphism in turbot was evident from the allele frequency data, with the mean number of alleles (MNA) being 6.330 ± 0.634.
The Shabla, Nesebar and Shkorpilovtsi populations had the smallest number of alleles per Smax-02 locus in comparison with all Black Sea populations previously investigated (Table 2). For the other loci, the data obtained were similar (Turan et al. 2019b;Bessonova and Nebesikhina 2019;Firidin et al. 2020). The most polymorphic locus in the current study was Sma-USC26 in the Shabla population (Table 2).  The expected number of alleles varied from 1.529 (Smax-02) to 6.374 (Sma-USC26) (Table 3). Heterozygosity is a major indicator of the level of genetic diversity in a population. With the values reflecting the level of observed heterozygosity (Ho) ranging from 0.107 at locus Smax-02 to 0.929 at locus B12-IGT14 (Table 3), the expected heterozygosity (He) of genetic diversity would vary from 0.406 at locus Smax-02 to 0.843 at locus Sma-USC26, respectively. The total mean of Ho was 0.638 ± 0.470, with maximum and minimum values of 0.684 and 0.582 in Shkorpilovtsi population and Nesebar population, respectively. The overall mean value of He was 0.685 ± 0.026.
The polymorphic information content (PIC) and the values of Shannon-Wiener index (I) provided a relatively broad value range of 0.324-0.825 and 0.596-2.015, respectively. The fixation index (Fis) varied between 0.006 and 0.769, with an average value of 0.094 to testify to a slight excess of heterozygotes in the fish group. The Garza-Williamson index was lowest at the B12-I GT14 locus in the Shkorpilovtsi population and highest at the Smax-02 locus in all populations analysed.
All investigated loci differed in terms of the Garza-Williamson index (M), within the range 0.10959 to 0.66667. The M average value in the investigated populations was 0.361 for Shabla, 0.335 for Nesebar and 0.351 for Shkorpilovtsi (Table 3).
The mean values of unbiased expected heterozygosity (uHe) in the analysed populations were similar (0.686±0.049 for Shabla population, 0.701±0.052 for Nesebar and 0.719±0.043 for Shkorpilovtsi), which marked the similar genetic diversity (Bessonova and Nebesikhina 2019).
Pairwise Fst comparisons showed low genetic differentiation between sampling sites (Table 4) and was always not significant. The presence of migrants per generation varied between 15 and 19% ( Table 4). The highest rate of gene flow (Nm = 19.498) was observed between Shabla and Nesebar populations.
A Mantel test revealed positive, but not significant relationships between the genetic and geographic distances (R 2 = 0.8273, P = 0.336).

Genetic diversity of mitochondrial DNA
A total of 17 haplotypes for COIII (563 bp) and 41 haplotypes for CR (432 bp) were identified. The sequence analyses of COIII recovered seven haplotypes for Shabla, two haplotypes for Shkorpilovtsi and eight haplotypes for the Nesebar population (Table 5). The majority of the identified COIII haplotypes originated from one prevalent set of haplotypes (Hap 2) following a single nucleotide substitution (Fig. 2). Four haplotypes were unique for the Shabla and Nesebar populations (Table 5)  The COIII haplotype diversity ranged from 0.389 to 0.766 in the three populations, with the highest value presented in the Nesebar population (Table 5). Nucleotide diversity ranged from 0.001 to 0.002. The result of mean COIII haplotype diversity between populations was 0.550, higher than the data presented in Turan et al. (2019b) for five populations in the Black and Marmara Seas (0.380). For the control region (CR), haplotypes (H1, 2 and 4) were common for the three populations analysed (Fig. 3) Unique CR haplotypes were observed in all populations, i.e. Shabla (13), Shkorpilovtsi (7) and Nesrbar (8) (Table 5, Fig. 3).
All three turbot populations analysed showed high levels of haplotype diversity (0.892-0.954), but low nucleotide diversity (0.004-0.006) in the mtDNA control   Table 5). The population pairwise Fst revealed an overall low level of genetic structure between the turbot populations. Non-significant Fst values were observed except for the comparison between Nesebar and Shabla populations (Table 6).
Genetic structure and phylogenetic trees were unable to detect genetic differentiation between sampling sites due to the low genetic differentiation for the COIII markers (0.002-0.003) and for CR (0.005-0.006) between haplotypes ( Table 7). The lowest genetic distance (0.002) was detected between Shkorpilovtsi and Shabla populations and the highest value (0.003) was found between Shabla and Nesebar populations, based on COIII analyses. Using CR, the highest genetic distance (0.006) was found between Shabla and Shkorpilovtsi populations (Table 7).
A Mantel test revealed no significant relationship between genetic and geographic distances (R 2 = 0.868 and P = 0.358 for COIII; R 2 = 0.881 and P = 0.324 for CR).
Tajima's D values were negative for all populations (Table 5), but statistically nonsignificant, with the exception of the Shabla population, indicating an excess of rare nucleotide site variants. The results of Fu's FS test, which is based on the distribution of haplotypes, showed negative values for all regions except Shkorpilovtsi (COIII marker), indicating an excess of rare haplotypes. The overall negative values resulting  Table 6. Scophthalmus maximus population pairwise F ST (based on COIII below diagonal) and based on CR (above diagonal), (p-values < 0.05*, < 0.01**, < 0.001***, ns not significant).

Nesebar
Shkorpilovtsi Shabla Nesebar -0.004 ns -0.012 ns Shkorpilovtsi -0.018 ns --0.017 ns Shabla 0.117** 0.122*** -from both tests indicate that there is an excess of rare mutations in the populations, but the excess is statistically non-significant. Based on the mean log likelihood values LnP(K), Bayesian clustering analysis, implemented in STRUCTURE, indicated K = 3 (LnP = -1695.64) as the most likely number of clusters ( Fig. 4 and Suppl. material 1: Fig. S7).

Morphological variability
Meristic data analyses provided strong statistically significant correlations only between TL, SL and W (Suppl. material 1: Figs S1, S3, S5), which was rather expected with respect to individual growth theory basics and moderate positive correlation between VFR and BPFR with very close values of correlation coefficient values for Shabla, Nesebar and Shkorpilovtsi populations (r = 0.55, p = 0.00164; r = 0.58, p = 0.00126; r = 0.057, p = 0.0265), (Table 8). Only the Nesebar population showed moderate positive correlation between AFR and VFR (r = 0.58, p = 0.01216). As the age of samples was not determined, to provide a closer look at growth model parameters, the length-weight relationship (LWR) was studied for comparison. In the Shabla population samples, LWR was best approximated with a 5 th order polynomial model (R 2 = 0.6732 -non-linear harmonic regression would provide even better approximation); in Nesebar samples, LWR was modelled with the power model: W = 5e -06*L 3.3006 (R 2 = 0.0.9608) and in the Shkorpilovtsi population, LWR was best approximated with 2 nd order polynomial model (R 2 = 0.9035). These results might highlight either certain differences in food availability, accounting for the weight being a highly variable parameter and/or regional environmental specifics or some background noises as seasonal processes or certain physiological stages, which also may have an impact on the genetic diversity of the studied populations.  Correlation analysis applied on morphometric characters revealed a certain linear functional pattern evident with slight differences amongst the samples, taken from the populations under investigation. The morphometric characters' correlation pattern in Shkorpilovtsi samples was slightly different from the pattern identified for Shabla and Nesebar samples (although it has to be noted that the correlation coefficient p values are sensitive to the number of samples n).
Analysis of similarity (ANOSIM) was carried out to identify statistically significant differences between the samples. The results showed that there was no statistically significant differences between Shkorpilovtsi and Shabla and Nesebar and Shabla samples (Significance Shabla-Shkorpilovtsi = 0.359; Significance Shabla-Nesebar = 0.8869, with an even distribution of high and low dissimilarity ranks in and between populations). However, there were statistically significant differences between the samples taken from Shkorpilovtsi and Nesebar populations (Significance Shkorpilovtsi-Nesebar = 0.063). In addition, the hierarchical clustering outcome (Suppl. material 1: Figs S2, S4, S6) showed that the clusters formed are mixed and specimens of all three populations are represented evenly and, furthermore, no specific outliers were identified. The latter is also visible from the non-metric multidimensional scaling (NMDS) plot (Fig. 5). No specific clusters were formed with respect to morphology and only single specimens are dispersed close to the major cluster formed due to specific morphometric characters deviation (TL, SL and W), which are generally considered to vary in more broad ranges.

Discussion
The results obtained in this study correspond to the creation of a database of the genetic diversity for turbot populations along the Bulgarian Black Sea coast. A comprehensive analysis of the acquired morphological and molecular data will enable a subsequent assessment of the impact of fishing on the structure of turbot populations. By knowing the genetic characteristics of valuable populations, we can monitor relatively easy changes in heritable traits and in the level of average genetic diversity (Olsson et al. 2007). Apparently, the selection of microsatellites with a range of polymorphism has led to a reduction in the risk of overestimating genetic variability, which might occur with the selective use of highly polymorphic loci. On the whole, the allelic diversity (mean number of observed alleles per locus) for populations along the Bulgarian Black Sea coast (3.559) is higher than that reported by Bessonova and Nebesikhina (2019) for the same loci in turbot populations from the Crimean Black Sea (3.430) and Azov Sea (3.420) and lower than species from the Caucasian Black Sea coast (3.623). The mean values for the obtained and expected heterozygosity (Ho = 0.638 and He = 0.685) are close to those indicated by previous studies for the same loci in turbot populations (Estoup et al. 1998;Iyengar et al. 2000;Bouza et al. 2002;Bessonova and Nebesikhina 2019;Turan et al. 2019b;Firidin et al. 2020). The relatively high heterozygosity and allelic diversity of these populations suggest that local gene flow takes place amongst them. The observed heterozygosity for Smax-02 and Sma-E79 loci (Shabla population), for Smax-02, Sma-USC26 and Sma-E79 (Nesebar) and Sma-E79 (Shkorpilovtsi) are significantly lower than the expected value in turbot populations (P < 0.05) (Table 3). There are several reasons that can explain this significant reduction in Ho compared with He, as high rate of migration, errors in reading alleles and inbreeding (Skaala et al. 2004;Li et al. 2009). The unbiased expected heterozygosity (uHe) showed a mean value of 0.702 ± 0.027, similar for those described by Bessonova and Nebesikhina (2019) for the Black Sea turbot (0.622 ± 0.086), which marked an equal level of the genetic diversity in the three populations investigated.
Nevertheless, a significant deviation was detected, from the HWE at a 0.05 α-level at all of the investigated loci with the exception of Smax-02 locus (Shabla population), Smax-02, B12-I GT14, Sma-USC26 and Sma-E79 loci (Nesebar population) and B12-I GT14 and Sma-E79 loci (Shkorpilovtsi population) (P < 0.05). Departures from HWE at other loci may be the result of founder and/or bottleneck effects followed by a high rate of inbreeding (Kaczmarczyk and Wolnicki 2016). A marker with a low Fis value will show the genetic diversity of a group more accurately, which means that this specific marker has a high discrimination power (Han et al. 2018). The positive and high value for Fis for the three populations was recovered for Smax-02 locus, based on the small number of the recorded alleles (Table 3) and indicates heterozygote deficiency that could be caused by inbreeding, presence of null alleles and admixture of distinct populations (Guo et al. 2013). The Smax-02 locus is also the locus with the least information (PIC is 0.324) and the lowest Shannon index (0.596).
The studies of the turbot populations, characterised by genetic diversity parameters (PIC and I), indicate that these values were high for six of the loci analysed (PIC and I higher than 0.6 and 1.3, respectively) indicating, thereby, high genetic diversity (Froufe et al. 2004;Fopp-Bayat 2010) (Table 3). The PIC values from each loci varies between populations analysed (0.622-0.825) and were comparable with the data for turbot population from the Black Sea (0.713-0.853) (Firidin et al. 2020) and from Turkey, France and South Korea (PIC 0.3-0.9) (Han et al. 2018) indicating the random loci were under independent selection pressure and seem to be ideal to be used as a divergence index. For detection of the bottleneck events, the Garza-Williamson index (M) was successfully employed. The data acquired (ranging from 0.335 to 0.361) across the studied populations are comparable with the data found for the Varna population (western Black Sea), based on the same loci (G-W index 0.4 on average) (Turan et al. 2019b), which suggests a sustainable reduction of genetic variation in this population as a result of founder and/or bottleneck effects (Kaczmarczyk and Wolnicki 2016). It should be pointed out that, in all studied populations, M indices were lower than 0.68 indicating a recent reduction in population size due to overfishing. The results of the Mantel test did not detect a clear correlation between geographic distance and the size of genetic differences existing between populations' pairs.
Generally, genetic differentiation based on microsatellites is considered as very weak. The mean value of Fst was 0.014 and the AMOVA results revealed that it was mostly related to a within-population variation (99.4%) rather than variation amongst populations (0.6%). The highest Fst value was observed between Shabla and Shkorpilovtsi = 0.016). The population pairwise Fst revealed an overall low level of genetic structuring between the turbot populations (Table 4). Based on the Fst values, it can be concluded that the rate of differentiation between populations is low and a high rate of migration (Nm =17.484, average) can be the main factor underlying this low differentiation. The Fst results were lower when compared to other Black Sea and Marmara populations (Fst = 0.249 according to Turan et al. 2019 b;Fst = 0.138 according to Firidin et al. 2020) and evidence of genetic sub-structuring along the Bulgarian Black Sea coast was not detected.
Mitochondrial DNA polymorphism is widely used to determine population structure, species differences and evolutionary relationships (Avise et al. 1987). The COIII gene of mtDNA is a genetic marker used for species identification in the genus Scophthalmus in the Black Sea (Turan et al. 2019a). The S. maximus populations displayed a genetic pattern typical of a recent population expansion due to its one common COIII haplotype (Hap 2) present across the range and the high number of unique haplotypes (Fig. 2). The Shabla population was characterised by the low levels of haplotype diversity (0.389) and nucleotide diversity (0.001) which is probably a result from a founder event or population bottleneck followed by rapid population growth (Xu et al. 2019). Population expansion in Shabla was also supported by highly significant negative values for both Tajima's D (-1.857, P < 0.05) and Fu's FS (-5.235 P < 0.05). The values of genetic diversity, as well as the distribution of private alleles are comparable amongst the populations, in spite of the low number of samples from Shkorpilovtsi.
The mean haplotype diversity and nucleotide diversity indices calculated as 0.550 and 0.001 for COIII and 0.930 and 0.005 for CR are similar to the data for COI and Cyt-b for the Black Sea turbot populations presented by Firidin et al. (2020). The genetic distances calculated, based only on COIII, present clear structure/correlation between populations according to the geographic distances (Table 7), but they are not significant according to the Mantel to test.
Studying more variable regions such as mtDNA CR to investigate the genetic variability across the Black Sea would give more valuable information about population structuring, based on mtDNA analyses (Firidin et al. 2020;Ivanova et al. 2020). The high levels of genetic diversity and low levels of nucleotide diversity observed in S. maximus, based on the mtDNA CR sequences, are in agreement with earlier studies of turbot populations from the Black Sea (Suzuki et al. 2004;Atanassov et al. 2011;Ivanova et al. 2020). A similar spatial pattern of distribution of genetic variability in mtDNA control region was also reported for other fishes (Xiao et al. 2008(Xiao et al. , 2014. High haplotype diversity in CR suggests large, stable, effective population sizes over time in the continental shelf fishes (Stepien 1999) and, in concurrence with low nucleotide diversity, it has been linked to population growth after a period of low effective population size (Grant and Bowen 1998). The inference of population expansion is further supported by the starlike patterns in Figure 3. The low nucleotide diversity in CR indicate a high number of closely-related haplotypes and suggest that this population may have undergone a recent expansion (Slatkin 1993;Mendez-Harclerode et al. 2007), supported also by non-significant and negative Tajima's D and Fu's Fs (Alcaraz and Gholami 2020).
Pairwise Fst values were very low, with the highest value for COIII observed between Shabla and Shkorpilovtsi populations and for CR between Shkorpilovtsi and Nesebar populations. Therefore, no significant genetic differentiation was evident between any populations and showed that S. maximus within the examined range constitutes a non-differential mtDNA gene pool. Results from STRUCTURE analysis (Fig. 4) also revealed some degree of admixture in the studied populations. In this regard, the observed heterozygote deficiency could also be caused by inbreeding depression, as a result of wild populations experiencing considerable reduction mainly because of overfishing and increased level of pollution.
The present mtDNA and microsatellite analyses of turbot populations along the Bulgarian Black Sea coast using seven microsatellites and two mtDNA markers give no evidence for genetic subdivision of this species in comparison with population genetic structuring found between the north and south Black Sea populations using microsatellites (Firidin et al. 2020). Based on both microsatellite and mtDNA analyses, Turan et al. (2019b) found one genetic unit of turbot in the Black Sea (Trabzon, Duzdge, Varna and Sevastopol populations). In the present study, S. maximus populations were considered to be a single stock, which probably could be attributed to this genetic unit.
A lack of correlation between genetic and geographic distances along the Bulgarian Black Sea coast was recorded. This result may be evidence of the hydrodynamic factors that have an effect on the dispersal potential of larvae phases and subsequently affect genetic differentiation (Vandamme et al. 2014(Vandamme et al. , 2020. S. maximus spawn with pelagic eggs (Vasil'eva 2007) and perform compensatory migration to the north against taking the caviar to the south of the currents along the Bulgarian coast (Karapetkova and Zhivkov 2006). Considering the lack of long active migration of species, the ocean currents might be responsible for this homogeneity in examined groups. In the North Sea, the turbot larvae drift on average 102 km (Barbut et al. 2019). Seawater temperature, food availability and coastal currents explain a significant component of geographically distributed genetic variation, suggesting that these factors act as drivers of local adaptation (Ruggeri et al. 2016;Diopere et al. 2017). The strength of local adaptation is dependent on the interaction between connectivity, population size and environmentally and human induced pressures (Bernatchez 2016).
Fishing pressure on turbot stock affects the size of the catches; however, there is no evidence of the impact of fishing on population genetic diversity. Genetic monitoring, in addition to the stock assessment, should also be carried out to track and identify all the potential population-genetic changes. Placing a restriction on the maximum catch size of 45 mm prevents the loss of rare alleles from older and larger individuals. It could also be an effective tool for protecting the genetic diversity. The regular monitoring and quota determination of the S. maximus populations are necessary to control the turbot population status. Close collaboration between molecular geneticists and fisheries biologists would be required to undertake extensive research into the recruiting processes of the marine populations and their possible implications for fisheries management and conservation (Hauser et al. 2002).

Conclusions
The DNA barcoding has proved to be an invaluable tool for tracking and monitoring of endangered populations, thus giving a sharper focus on the strategic conservation of distinct genetic stocks and mitigation on human impacts along their range. The molecular characterisation and analysis of the genetic structure in turbot populations along the Bulgarian Black Sea coast contributes to the considerable knowledge about the levels and genetic diversity distribution patterns. Microsatellite and mitochondrial markers both indicated close genetic relationships between populations. It should, however, be pointed out that, in the examined populations, a high level of genetic diversity was observed. The lack of strong population genetic structure was probably due to the small geographic distances and high gene flow between them. The derived low values of the G-W index are specific for reduced populations and may indicate a dramatic decrease in the population size in the past. The applied molecular approach proves critical for any rigorous monitoring of the impacts of overexploitation and genetic management of threatened fish species along the Bulgarian Black Sea coast. The results confirmed the high effectiveness of the use of different types of markers for performing genetic analysis and relevant provision of reliable information with regard to the genetic diversity within the turbot populations. The genetic characteristics of turbot populations, revealed in this study, will provide useful information for continuous and effective resource management. Moreover, constant monitoring may be needed to maintain the high level of genetic diversity in natural populations.
Finally, the underlying rationale behind the adoption of a more integrated approach (genetic and morphological) to the study of S. maximus populations in the Bulgarian Black Sea coast will provide more accurate assessment of the population structure as well as it will facilitate detection of any probable changes in gene pools of the wild populations in connection with their more effective management. Therefore, a proper understanding and a sound knowledge of the level and distribution of genetic diversity in turbot is an important prerequisite for successful sustainable development and conservation strategies to preserve their evolutionary potential.