Spatial Distribution and Genetic Diversity of Turbot ( Scophthalmus maximus , Linnaeus, 1758) in Bulgarian Black Sea Waters Relative to Fishing Pressure and Their Abiotic Environment

: The present study examined the genetic diversity and spatial distribution of turbot ( Scoph-thalmus maximus ), an economically important species on the Bulgarian Black Sea coast. Maximum entropy (MaxEnt) modeling software Version 3.4.4. was utilized to develop a habitat suitability model for S. maximus in the Bulgarian Black Sea region. Data collected via demersal and pelagic surveys and genetic sampling from 2017 to 2021 were utilized to link species occurrence localities with selected abiotic factors. Our ﬁndings showed that the species’ habitat preferences are strongly inﬂuenced by temperature and dissolved oxygen, and projections based on simpliﬁed climatic scenarios indicated potential distribution shifts and a substantial reduction in reproduction habitats in the northern region. The assessment of genetic diversity was based on mtDNA COIII sequencing; MtDNA revealed a low level of polymorphism in all analyzed populations. The extensive ﬁshing pressure may have increased the likelihood of genetic and population bottlenecks and a consequent decline in genetic diversity in the Shabla, Nesebar, and Tsarevo populations. The Tajima’s D values for the latter indicated that turbot underwent a bottleneck followed by rapid population expansion. Our ﬁndings are essential for the conservation and effective management of S. maximus stocks in the region.


Introduction
In the marine environment, the greatest genetic threats appear to be the extinction of genetically unique subpopulations and the loss of genetic diversity through declines in abundance due to overfishing and climate change [1].The human overexploitation of natural resources has elevated conservation and management to one of modern society's most pressing concerns, particularly in the case of highly sensitive marine ecosystems [2].Fishing is one of the most significant causes of decline in populations of ocean species [3].In addition to depleting populations [4], overfishing can erode the age and size structure as well as the spatial distribution of stocks, making them more vulnerable to environmental fluctuations.The distribution of fish in the maritime environment follows distinct patterns.Fish are more prevalent in areas with sufficient resources to support higher densities.In general, the geographic dispersal of fish increases in proportion to abundance or as a density-dependent selection response, in which fish preferentially occupy a chosen location until maximum density is reached, at which point, they disperse into more peripheral regions [5].Fishing typically focusses on high-density core areas in order to maximize catch while minimizing effort.Hence, identifying core density areas and understanding the link between fish and their immediate environment are essential for an ecosystem-based approach to fisheries management [6].
Additionally, a growing number of studies in recent years indicate that substantial size-selective fishing can trigger evolutionary changes in the life-history traits of the harvested population, which can have negative implications for populations, ecosystems, and fisheries.Overfishing has the potential to generate three types of genetic change: changes in population structure and connectivity, selection-induced genetic changes, and genetic diversity loss [7].Despite the fact that data on stock abundance and age composition for economically important species have been collected for decades [8,9], fisheries statistics are based on short-term responses to exploitation.Genetic analysis, on the other hand, can provide information about long-term processes related to a species' adaptive capacity to respond to changes in environmental conditions.
By displaying the genetic level, patterns of dispersal (connectivity), and demographics (past and present), population genetic methods enable scientists to comprehend species' responses to ecological changes and anthropogenic stressors.Human activities have the potential to reduce the genetic diversity of a population in a very short period of time [10].The overexploitation of a wide variety of marine fish species has been linked to a decline in genetic diversity [11].
For economically important species, however, which are facing substantial fishery pressures throughout their entire distribution range, genetic diversity is likely declining globally [7].For many years, mitochondrial genetic markers have been used in phylogenetic and population genetics research.The loss of genetic diversity in mitochondrial genomes as a result of population declines can provide insights into the conservation status of fish species populations [7].
Fisheries management and conservation priorities, including fishing regulations, vary greatly around the world [8,9,12], as does the impact of fishing pressures, necessitating local studies of fish populations [7].Furthermore, the comprehension and effective conservation and management of ecosystems necessitate a solid understanding of population dynamics within a spatial and temporal context.In line with the latter, typically, dispersal and recruitment success are determined in the early development stages for many marine species.The final dispersal pattern, larval survival, and overall success of populations are influenced by a range of environmental conditions along with biological aspects such as ecophysiology, behavior, and reproductive strategy [13].
In a broader context, the application of genomic methods has the potential to enhance stock assessments through the facilitation, improvement, and estimation of crucial parameters.Additionally, these approaches offer valuable insights into stock structure, genetic diversity within populations, and connectivity, all of which are recognized as significant factors in the fields of conservation and resource management [14].
Turbot (Scophthalmus maximus, Linnaeus, 1758), as a transboundary stock, has played an important role in the fisheries of all Black Sea countries for decades, and is regarded as one of the most important key species in the Black Sea ecosystem.As with most benthic fishes, turbot display shifting reliance on coastal regions to complete their life-cycles, thereby exhibiting varying degrees of connectivity and environmental selection [15,16].Turbot is a top predator species, mainly feeding on fish, crustaceans, polychaetes, and mollusks.The distribution and migration of turbot across the continental shelf is evidently seasonal.The spring migration towards the coastal area typically starts in March.Spawning occurs in April to the middle of June with a peak in the first half of May, at depths ranging from 15 to 50 m.During the summer, turbot inhabit depths over 50 m and are spread out across a larger area.Migration towards the coastal feeding areas occurs in the autumn.In addition to the seasonal migrations offshore and towards the coast, adult turbot exhibit a tendency to move northwards.During all seasons, juvenile turbot are typically restricted to coastal areas.Considering the latter, it is evident that throughout its entire biological cycle, the turbot utilizes a variety of habitat types, including reproduction habitats (coastal waters, down to 15 m depth), growing habitats (near the seafloor, down to 50-60 m depth), and wintering habitats (near the seafloor, down to 60-80 m depth).To these depths, the turbot populates nearly all habitat types, but prefers sandy and muddy environments [17,18], and consequently, the environmental heterogeneity throughout the extensive distribution range of turbot allows for a variety of spatial differentiation and adaptation [15,16,19].Accordingly, it is anticipated that turbot populations will microevolve under dispersal limitation using models of isolation by distance and local adaptation [20].
Since 2007, Bulgaria has implemented the Total Allowable Catch (TAC) regime for turbot in an effort to prevent overfishing, and to maintain and restore the effective population size of the species inhabiting the Bulgarian waters of the Black Sea in order to provide high long-term yields that are consistent with the Maximum Sustainable Yield (MSY).Despite all conservation and technical measures adopted at the regional level to prevent overfishing and ensure sustainable exploitation, 75 percent of the assessed Black Sea stocks, including turbot, are reported to be fished at levels that are biologically unsustainable [8].
Although turbot has been the subject of extensive research, there is still a lack of understanding regarding its population structure throughout its entire distribution range.Current knowledge on wild populations is primarily limited to the Atlantic and Baltic regions.Furthermore, investigations into the genomic characteristics linked to environmental factors have only been conducted in a few populations [21,22] and with a limited number of markers [16,23].Novel findings about the adaptation of turbot have been documented [24], particularly in the Black Sea, highlighting the occurrence of parallel evolution in regions characterized by comparable environmental conditions, specifically temperature.Generally, the natural distribution of turbot exhibits the concurrent influence of robust neutral evolutionary factors and adaptive selection in geographically distinct populations [24].Accordingly, a more advanced conceptualization of management units can be developed by considering demographic variables and the adaptive capacity of turbot to respond to environmental changes throughout its whole distribution range.
There are few case studies addressing turbot stock assessments [25] along the Bulgarian Black Sea coast, and data on the species' spatiotemporal distribution and habitat suitability are limited.For a small number of populations, genetic diversity throughout the Bulgarian Black Sea coast was assessed using mitochondrial markers [26,27].However, no studies have been conducted to analyze the temporal decline in genetic diversity and the loss of genetic variation associated with harvesting in the past due to a lack of genetical monitoring.
In line with the latter, the primary objectives of our research were to provide an accurate and objective understanding of the spatial distribution and core density areas, the genetic diversity and variability of turbot populations along the Bulgarian Black Sea coast, as well as the environmental factors that govern the species' geographic dispersal.Providing new insights on turbot populations' genetic structure and spatial dispersal involving environmental, ecological, and genetic mechanisms can greatly support and adjust existing management strategies towards sustainable exploitation and conservation.

Study Area
The study area covered the Bulgarian Black Sea exclusive economic zone (EEZ) (Figure 1).Since the sampling data utilized in the spatial distribution modeling were only available for the Bulgarian coastal and shelf waters, the relative spatial extent was limited to the area where sufficiently large number of observations of S. maximus were available.

Spatial Modeling Software
Over the past two decades, significant advancements have occurred in the domain of species distribution modeling, resulting in the availability of numerous methodologies.One significant disparity between modeling approaches lies in the type of species data they utilize, the complexity of the algorithms, the underlying assumptions, and the simplicity of use and model interpretation [28,29].Regression methods such as generalized linear or additive models (GLMs or GAMs) and ensembles of regression trees are routinely utilized when species presence/absence and/or abundance have been systematically collected, typically via formal biological surveys [29,30].
In practice, however, the documenting of recurring observations of absence data can be challenging [31], and numerous spatial distribution modeling techniques, including regression models, various machine learning algorithms such as maximum entropy MaxEnt [32] models, presence and background learning (PBL) [33], and ensemble models, have been devised to address the issue of handling presence-only (PO) data.
The maximum entropy (MaxEnt) modeling package software Version 3.4.4.[32] was chosen to fit a correlative model for the studied species.The latter option was deemed appropriate for our case study due to several reasons.Generally, the predictive performance of Maxent has continuously been proven to be comparable to the most effective techniques [29], and it has been successfully employed to model the biogeographic distributions of species in various regions of the world [34][35][36][37][38][39][40][41][42][43][44][45].
Another advantage of MaxEnt is that its functionality performs variable transformation, maximum entropy model fitting, and lasso regularization to avoid model overfitting [29,32,[46][47][48].Moreover, the resulting output of the model is relatively straightforward to interpret.Another rationale for using a modeling technique based on PO (presence-only) data was the inclusion of auxiliary data, such as genetic samples and pelagic survey locations, where abundance data were unavailable.

Spatial Distribution Modeling 2.2.1. Spatial Modeling Software
Over the past two decades, significant advancements have occurred in the domain of species distribution modeling, resulting in the availability of numerous methodologies.One significant disparity between modeling approaches lies in the type of species data they utilize, the complexity of the algorithms, the underlying assumptions, and the simplicity of use and model interpretation [28,29].Regression methods such as generalized linear or additive models (GLMs or GAMs) and ensembles of regression trees are routinely utilized when species presence/absence and/or abundance have been systematically collected, typically via formal biological surveys [29,30].
In practice, however, the documenting of recurring observations of absence data can be challenging [31], and numerous spatial distribution modeling techniques, including regression models, various machine learning algorithms such as maximum entropy Max-Ent [32] models, presence and background learning (PBL) [33], and ensemble models, have been devised to address the issue of handling presence-only (PO) data.
The maximum entropy (MaxEnt) modeling package software Version 3.4.4[32] was chosen to fit a correlative model for the studied species.The latter option was deemed appropriate for our case study due to several reasons.Generally, the predictive performance of Maxent has continuously been proven to be comparable to the most effective techniques [29], and it has been successfully employed to model the biogeographic distributions of species in various regions of the world [34][35][36][37][38][39][40][41][42][43][44][45].
Another advantage of MaxEnt is that its functionality performs variable transformation, maximum entropy model fitting, and lasso regularization to avoid model overfitting [29,32,[46][47][48].Moreover, the resulting output of the model is relatively straightforward to interpret.Another rationale for using a modeling technique based on PO (presenceonly) data was the inclusion of auxiliary data, such as genetic samples and pelagic survey locations, where abundance data were unavailable.
MaxEnt model input data consist of a list of species' PO locations and a set of environmental predictors across a user-defined, grid-cell-divided spatial extent with an option for the manipulation of background/pseudoabsence data.MaxEnt extracts a sample of background locations, compares them to the presence locations, and generates an estimate of the probability of the presence of a species (or relative environmental suitability) ranging from 0 (least likely) to 1 (most likely) [32,48].

Species Presence Data
Data on species PO localities were compiled from demersal surveys records [49][50][51][52][53][54][55][56][57][58] aimed at the stock assessment of turbot, pelagic surveys where the accidental bycatch of turbot occurred, and genetics sampling.Survey expeditions were conducted in the period of 2017-2021 (Table 1).Both surveys' sampling designs, as well as genetic sampling, encompassed the Bulgarian Black Sea's coastal (0-30 m) and shelf waters (30-200 m), providing an accurate representation of S. maximus presence records within the area of interest.As a result, the species dataset was constructed by taking into consideration the PO localities (384 records in all).

Selection of Abiotic Variables
The abiotic variables chosen as predictors were based on previously documented impact extent and ecological concepts linked to species preferences form a particular habitat, as well as those influencing the biomass and recruitment variability of flatfishes [59][60][61][62][63][64][65][66][67][68][69].Accordingly, the water temperature, salinity, dissolved oxygen, and current speed were selected as the most relevant.
The Copernicus Marine Environmental Service (CMEMS) data portal (https://marine.copernicus.eu/accessed on 15 February 2023) was used for obtaining monthly mean environmental data layers for the chosen variables from April to December in the 2017to-2021 period (corresponding to the sampling expeditions).The latter are generated by numerical simulation models that integrate in situ and satellite data for the Black Sea profile (the hydrodynamic NEMO (Nucleus for European Modelling of the Ocean [70]) and the BAMHBI (Biogeochemical Model for Hypoxic and Benthic Influenced Areas [71][72][73])).The data were compiled using two products: the Black Sea Physics Reanalysis [74] for the period of April to December 2017 to 2020, with a spatial resolution of 0.037 • × 0.028 • , and the Black Sea Physics Analysis and Forecast [75] for the period of April to December 2021, with a spatial resolution of 0.025 The datasets were initially resampled to a finer resolution (0.025 • × 0.025 • ), and then, averaged in the water column and throughout the period of the study in MATLAB [76] over the entire Black Sea region.The resulting data were then imported into ArcGIS 10.2, where they were further reduced to match the spatial extent of the study area and extracted in ESRI (Environmental Systems Research Institute) ASCII grid format to be utilized in MaxEnt software to fit the model.

Model Evaluation and Validation
The species distribution model (SDM) was cross-validated with ten replicate model runs (using 10-fold cross-validation) for an adequate evaluation of its predictive performance, and 25% of PO data were set aside to be utilized as a randomly selected test sample for every model run.
The performance of the resulting SDM was assessed using the area under the curve (AUC) of a receiver operating characteristic (ROC) [77], which represents an overall measure of the model's performance across all thresholds and strengths of prediction, and is thus considered threshold-independent [78], the maximized Cohen's kappa [79], and true skill statistics (TSS) was introduced as a prevalence-independent alternative to Cohen's kappa to assess model reliability (the agreement between predictions and known occurrences, across various binary thresholds) [78].
Additionally, the mean relative probabilities of species occurrence (in terms of suitability) were cross-validated with the species' field-measured catch per unit effort (CPUE) to further support the models' applicability.2).Prior to the subsequent analyses, dorsal fin tissue samples were stored in 96% ethanol at 4 • C for DNA analysis.DNA was extracted using the DNeasy Blood & Tissue Kit (QIA-GEN).Prior to PCR amplification, all DNA extracts were analyzed via gel electrophoresis to determine the DNA's quality.

Statistical Analyses
The mtDNA sequence data were analyzed using MEGAX software (https://www.megasoftware.net/downloads/dload_win_gui,accessed on 15 February 2023) [81].DNA sequence polymorphisms comprising the largest number of haplotypes (h), number of segregating sites (S), haplotype diversity (Hd), nucleotide diversity (π), and neutrality test (Tajima's D and Fu's Fs) values for each population were estimated using DNASP 6.12.03 [82].Haplotype networks were visualized via PopART network analysis software (https://github.com/jessicawleigh/popart,accessed on 15 February 2023) [83] using the TCS network inference method [84].The MaxEnt grid output (maps) employ a gradient color scale to reflect the relative mean predictive probability (within the range 0-1) of the most favorable habitat for the species under consideration [32].The resulting model (Figure 2) clearly displayed a broad area where environmental conditions were deemed suitable, particularly within the 30-80 m depth range (clearly demonstrating the reproduction, growing, and wintering habitats due to the seasonal resolution of sampling), broadening across the entire region in front of the of the Bulgarian Black Sea coast.

Statistical Analyses
The mtDNA sequence data were analyzed using MEGAX software (https://www.megasoftware.net/downloads/dload_win_gui,accessed on 15 February 2023) [81].DNA sequence polymorphisms comprising the largest number of haplotypes (h), number of segregating sites (S), haplotype diversity (Hd), nucleotide diversity (π), and neutrality test (Tajima's D and Fu's Fs) values for each population were estimated using DNASP 6.12.03 [82].Haplotype networks were visualized via PopART network analysis software (https://github.com/jessicawleigh/popart,accessed on 15 February 2023) [83] using the TCS network inference method [84].The MaxEnt grid output (maps) employ a gradient color scale to reflect the relative mean predictive probability (within the range 0-1) of the most favorable habitat for the species under consideration [32].The resulting model (Figure 2) clearly displayed a broad area where environmental conditions were deemed suitable, particularly within the 30-80 m depth range (clearly demonstrating the reproduction, growing, and wintering habitats due to the seasonal resolution of sampling), broadening across the entire region in front of the of the Bulgarian Black Sea coast.For the entire Black Sea region, a strong negative spatial correlation between temperature and salinity values (averaged in the water column and over the period of the case study) was detected (r = −0.861,p~0); however, given that Maxent is robust to predictor collinearity when training the model and accounts for redundant variables [85], the exclusion of highly correlated variables has little to no effect on model improvement.Therefore, the original selection of predictor variables was retained.Furthermore, data collinearity is likely to have an effect when the model is applied to new environmental conditions if the variables' correlation and collinearity shift are not accounted for [85].For the entire Black Sea region, a strong negative spatial correlation between temperature and salinity values (averaged in the water column and over the period of the case study) was detected (r = −0.861,p~0); however, given that Maxent is robust to predictor collinearity when training the model and accounts for redundant variables [85], the exclusion of highly correlated variables has little to no effect on model improvement.Therefore, the original selection of predictor variables was retained.Furthermore, data collinearity is likely to have an effect when the model is applied to new environmental conditions if the variables' correlation and collinearity shift are not accounted for [85].

Species Distribution Modeling
In the context of variables' importance to the MaxEnt model, temperature showed the highest gain when used isolated, and therefore appears to have the most important information on its own (Table 3).Dissolved oxygen was the environmental variable that reduced the gain the most when it was omitted, and hence, seemed to contain the most information that was not present in the other variables.Salinity and current speed had little to no effect on fitting the training and test data models.The importance of temperature, identified as a major factor that governs turbot spatial dispersal, is of superior concern when accounting for the effects of climate change; hence, two distinct scenarios were examined to perceive whether rises of 0.5 • C and 1 • C in the annual mean temperatures could potentially lead to distribution shifts (Figure 3).
In the context of variables' importance to the MaxEnt model, temperature showed the highest gain when used isolated, and therefore appears to have the most important information on its own (Table 3).Dissolved oxygen was the environmental variable that reduced the gain the most when it was omitted, and hence, seemed to contain the most information that was not present in the other variables.Salinity and current speed had little to no effect on fitting the training and test data models.The importance of temperature, identified as a major factor that governs turbot spatial dispersal, is of superior concern when accounting for the effects of climate change; hence, two distinct scenarios were examined to perceive whether rises of 0.5 °C and 1 °C in the annual mean temperatures could potentially lead to distribution shifts (Figure 3).The projected results demonstrated a significant decrease in the relative probability of species occurrence in the coastal areas in front of Varna Bay and to the north, as well as a potential core area shift to deeper waters-a very concerning finding given that the Black Sea has hypoxic/anoxic conditions that prevent fish stock distribution beyond a 100 m depth.

Evaluation of SDMs Performance
The test omission rates and predicted areas as a function of the cumulative threshold for all replicate runs (including the averaged) were close to the predicted omission, according to the definition that a model that distinguishes well between presences and absences would have low omission (high sensitivity) and low commission (high specificity) The projected results demonstrated a significant decrease in the relative probability of species occurrence in the coastal areas in front of Varna Bay and to the north, as well as a potential core area shift to deeper waters-a very concerning finding given that the Black Sea has hypoxic/anoxic conditions that prevent fish stock distribution beyond a 100 m depth.

Evaluation of SDMs Performance
The test omission rates and predicted areas as a function of the cumulative threshold for all replicate runs (including the averaged) were close to the predicted omission, according to the definition that a model that distinguishes well between presences and absences would have low omission (high sensitivity) and low commission (high specificity) [84].The model performance metrics for the 10 replicate runs were calculated to be: AUC = 0.906 ± 0.016; maximized Cohen's kappa = 0.915 ± 0.006; and TSS = 0.687 ± 0.088.
In general, AUC values greater than 0.9 are very good [85,86], kappa > 0.75 is considered excellent, and TSS > 0.5 indicates good agreement between predictions and known occurrences [79].The MaxEnt predicted habitat suitability map was overlaid with the field-measured CPUE for all PO localities used in fitting the training data model to validate the expected association of the highest probabilities of occurrence of the studied species with areas where the species' density and performance were evidently high (Figure 4).
[84].The model performance metrics for the 10 replicate runs were calculated to be: AUC = 0.906 ± 0.016; maximized Cohen's kappa = 0.915 ± 0.006; and TSS = 0.687 ± 0.088.In general, AUC values greater than 0.9 are very good [85,86], kappa > 0.75 is considered excellent, and TSS > 0.5 indicates good agreement between predictions and known occurrences [79].The MaxEnt predicted habitat suitability map was overlaid with the field-measured CPUE for all PO localities used in fitting the training data model to validate the expected association of the highest probabilities of occurrence of the studied species with areas where the species' density and performance were evidently high (Figure 4).A good spatial overlap between the two grids-CPUE and the average relative predicted probabilities-was detected; nevertheless, the most abundant catches observed during the surveys were recorded in the southern region.

Mitochondrial DNA Analysis
Sequences from 133 individuals were obtained, and a total of 25 haplotypes for COIII (563 bp) were identified among the populations.
The number of haplotypes ranged from six in Shkorpilovtsi to twelve in the Nesebar population (Table 4).The level of genetic diversity was determined for S. maximus from each sampling location.The haplotype diversity (Hd) and nucleotide diversity (π) indices ranged from 0.363 ± 0.098 to 0.783 ± 0.057 and from 0.00128 ± 0.00044 to 0.00233 ± 0.00037, respectively.Overall, most locations exhibited moderate to high haplotype diversity (0.363-0.783) due to a large number of unique haplotypes.The mean haplotype diversity (Hd = 0.577) of turbot along the Black Sea coast was close to that reported by [26] for the A good spatial overlap between the two grids-CPUE and the average relative predicted probabilities-was detected; nevertheless, the most abundant catches observed during the surveys were recorded in the southern region.

Mitochondrial DNA Analysis
Sequences from 133 individuals were obtained, and a total of 25 haplotypes for COIII (563 bp) were identified among the populations.
The number of haplotypes ranged from six in Shkorpilovtsi to twelve in the Nesebar population (Table 4).The level of genetic diversity was determined for S. maximus from each sampling location.The haplotype diversity (Hd) and nucleotide diversity (π) indices ranged from 0.363 ± 0.098 to 0.783 ± 0.057 and from 0.00128 ± 0.00044 to 0.00233 ± 0.00037, respectively.Overall, most locations exhibited moderate to high haplotype diversity (0.363-0.783) due to a large number of unique haplotypes.The mean haplotype diversity (Hd = 0.577) of turbot along the Black Sea coast was close to that reported by [26] for the same species (Hd = 0.550), and higher than that (0.380) published by [27]; however, the nucleotide diversity was relatively low, ranging from 0.00128 to 0.00233.Due to its central position, the haplotype network analyses identified the haplotype Hap_1 as a potential maternal ancestral haplotype for the studied populations (Figure 5).Overall, five shared haplotypes were identified (Hap_1, Hap_4, Hap_7, Hap_9, and Hap_16) and 20 unique haplotypes (Shabla-5, Nesebar-7, Tsarevo-6, and Shkorpilovtsi-2) were found to derive from the above-mentioned haplotypes.
same species (Hd = 0.550), and higher than that (0.380) published by [27]; however, the nucleotide diversity was relatively low, ranging from 0.00128 to 0.00233.
Due to its central position, the haplotype network analyses identified the haplotype Hap_1 as a potential maternal ancestral haplotype for the studied populations (Figure 5).Overall, five shared haplotypes were identified (Hap_1, Hap_4, Hap_7, Hap_9, and Hap_16) and 20 unique haplotypes (Shabla-5, Nesebar-7, Tsarevo-6, and Shkorpilovtsi-2) were found to derive from the above-mentioned haplotypes.Both Tajima's D and Fu's Fs tests showed negative values for all populations (Table 4).Tajima's D values were statistically significant for three of the studied populations (Shabla, Nesebar, and Tsarevo).These results suggest that the lineages may have undergone a process of population expansion after a bottleneck event, which is evidenced by the high number of unique haplotypes.

Discussion
The key findings of our study may have been shaped by background events that have been unfolding over the last two decades.The unexpectedly low nucleotide diversity (π ≤ 0.002) observed across all studied populations may be attributable to the commercial overexploitation of wild stocks [87].The latter hypothesis was further strengthened by possible local depletion of the species under study identified in the northern region during the validation of the S. maximus spatial distribution model with the field-measured CPUEs recorded during the survey expeditions over the last five years.
The low haplotype and nucleotide diversity of the Tsarevo population (0.363 and 0.00139, respectively) was most likely the result of a founder event or population bottleneck [88].Such low nucleotide diversity of the COI mitochondrial gene has been discovered in heavily exploited fish stocks (Lophius budegassa, Merluccius merluccius, Mullus barbatus, and Mullus surmuletus) from the Balearic Islands [7], and is typically leading to population bottleneck and a subsequent genetic reduction in the total gene pool [3].
The negative and statistically significant Tajima's D values that were observed in three of the studied populations (Shabla, Nesebar, and Tsarevo) most likely indicate that several haplotypes in these populations were eliminated, resulting in a bottleneck effect [89].The latter has been reported to typically occur as a result of natural or human activities.Accordingly, overfishing can deplete a population by removing a large number of individuals, resulting in a population bottleneck, which may cause a loss of genetic diversity and lead to population decline.Consequently, the low level of nucleotide diversity observed in turbot populations along the Bulgarian Black Sea coast is likely due to population declines caused by highly selective and cost-effective fishing pressures, which target reproduction and growing habitats (30 to 50-60 m depths).The decline in genetic diversity is generally associated with the long-term adaptation and survival of a population and can further restrict the capacity of individuals to adapt to novel environmental conditions [90].
In line with the latter, species spatial distribution model validation with the fieldmeasured CPUE revealed persistently low catches (with regard to stock density) in northern coastal and shelf areas (up to 50-60 m depths) over the past five years, which could be attributed to overfishing and the applied physical disturbance in coastal benthic habitats while targeting other commercially important benthic species such as Rapana venosa and Mullus barbatus.Consequently, the locations with the lowest captures may indicate the partial effects of local depletion and habitat fragmentation, which are more evident in turbot growing habitat depths in the northern region, and the effective population size is likely to decrease as a direct result.
In addition to excessive fishing pressure and other stressors of anthropogenic and natural origin, the variables contributing most to the spatial distribution of the studied species are the temperature and dissolved oxygen, which raises serious concerns regarding how the species' stock's spatial integrity and dispersal would change as a result of global climate change.Both climatic scenarios (annual mean temperature rises of 0.5 • C and 1 • C) were intended to highlight prospective future habitat shifts of temperature-sensitive species, although species' tolerance to temperature rises and their overall influence on the abiotic environment, including species' biotic interactions, must still be addressed to yield more accurate outcomes in terms of species' responses to global warming.
Applying distribution modeling to marine benthic species to predict shifts in distribution patterns in response to climate change is rather challenging.Our simplified climatic scenario projections demonstrate a significant shift in the core areas and species' geographic dispersal, and furthermore, anoxic layers represent a natural barrier to spatial dispersion due to the Black Sea's specific vertical structure [91,92].
In addition to the annual mean abiotic variables that may be linked to the relative abundance and distribution of fish stocks, it is important to take into account other factors, such as biotic interactions (e.g., interactions with prey, predators, and competitor species) and time-lagged abiotic conditions on a finer temporal scale [93].In certain practical scenarios, it has been observed that ecological variables with time-lag effects, such as sea surface temperature and salinity, and proxies to biotic interactions (e.g., chlorophyll-a for planktivorous fishes), could facilitate a more accurate representation of fish spatial dispersal and productivity areas in response to environment alterations.The latter could potentially explain the spatio-temporal variability in the relative abundance and distribution of fish species more clearly; nevertheless, it also requires more frequent assessments and data collection for the species under consideration.
Processes such as species interactions, habitat change, and dispersal range and barriers usually remain unaccounted for in prognostic spatial distribution modeling as a response to climate change.Furthermore, changes in biotic interactions can significantly affect species distributions [94,95], which might be further altered by invasions of non-native species favored by climate change, resulting in novel combinations of species [96].
Furthermore, the examination of the physiological mechanisms underlying individual responses can facilitate the integration of various factors such as rising temperature, acidification, and hypoxia.This integration allows for a comprehensive assessment of the cumulative impact on organisms and the energetic costs associated with individual acclimation.Consequently, the inclusion of mechanistic models in predictive analyses may yield more accurate forecasts regarding species' spatial distribution in response to climate change and their ability to adapt to novel environmental conditions.

Conclusions
In conclusion, considering our research results, conducting mechanistic research to study species population dynamics and underlying development processes, including biotic interactions, is of fundamental importance, given the species preferences of habitat.Studying the sensitivity of S. maximus to temperature and oxygen availability can support the study of their response to climate change and provide more unbiased knowledge on the expected patterns in their spatial dispersal.Furthermore, given the findings of the effects of fishing pressure on S. maximus genetic diversity and stock spatial integrity, stock management strategies must focus on the development of an integrated ecosystem-based management approach to benthic fisheries, considering the impact of other benthic-species fisheries on turbot stock spatial dispersal, population dynamics, and species habitat integrity.

Figure 1 .
Figure 1.Sampling locations, with yellow dots displaying demersal and multispecies survey sites and red dots representing genetic marker sampling locations (annotated with S1-S4) for the study period (2017-2021).

Figure 1 .
Figure 1.Sampling locations, with yellow dots displaying demersal and multispecies survey sites and red dots representing genetic marker sampling locations (annotated with S1-S4) for the study period (2017-2021).

2. 3 .
Mitochondrial DNA Analysis 2.3.1.Sample Collection and DNA Extraction Tissue samples were obtained from 133 samples collected by a local fishing vessel in the Black Sea near Shabla, Shkorpilovtsi, Nesebar, and Tsarevo during the period of 2019 to 2021 (Figure 1, Table

Figure 2 .
Figure 2. MaxEnt Habitat suitability map of S. maximus in the Bulgarian Black Sea coastal and shelf waters (with dark red indicating the highest possible probability of suitable conditions, green indicating conditions equal to those where species were found, and lighter shades of purple corresponding to low predicted probability of suitable environment).

Figure 2 .
Figure 2. MaxEnt Habitat suitability map of S. maximus in the Bulgarian Black Sea coastal and shelf waters (with dark red indicating the highest possible probability of suitable conditions, green indicating conditions equal to those where species were found, and lighter shades of purple corresponding to low predicted probability of suitable environment).

Figure 3 .
Figure 3. MaxEnt Habitat suitability map projection under two different climate scenarios: (a) S. maximus habitat suitability in the Bulgarian Black Sea coastal and shelf waters if annual mean temperatures rise by 0.5 °C and (b) if annual mean temperatures rise by 1 °C (with dark red indicating the highest possible probability of suitable conditions, green indicating conditions equal to those where species were found, and lighter shades of purple corresponding to low predicted probability of suitable environment).

Figure 3 .
Figure 3. MaxEnt Habitat suitability map projection under two different climate scenarios: (a) S. maximus habitat suitability in the Bulgarian Black Sea coastal and shelf waters if annual mean temperatures rise by 0.5 • C and (b) if annual mean temperatures rise by 1 • C (with dark red indicating the highest possible probability of suitable conditions, green indicating conditions equal to those where species were found, and lighter shades of purple corresponding to low predicted probability of suitable environment).

Figure 4 .
Figure 4. MaxEnt model-predicted probabilities grid vs. field-measured CPUE kg/trawl grid (the red gradient dots) of S. maximus in the Bulgarian Black Sea coastal and shelf waters during the survey expeditions carried out in the 2017-2021 period [49-58] (with dark red indicating the highest possible probability of suitable conditions, green indicating conditions equal to those where species were found, and lighter shades of purple corresponding to low predicted probability of suitable environment).

Figure 4 .
Figure 4. MaxEnt model-predicted probabilities grid vs. field-measured CPUE kg/trawl grid (the red gradient dots) of S. maximus in the Bulgarian Black Sea coastal and shelf waters during the survey expeditions carried out in the 2017-2021 period [49-58] (with dark red indicating the highest possible probability of suitable conditions, green indicating conditions equal to those where species were found, and lighter shades of purple corresponding to low predicted probability of suitable environment).

Figure 5 .
Figure 5. TCS network of S. maximus based on COIII haplotypes.The sizes of the circles indicate the number of individuals with each given haplotype.Small lines illustrate the substitutions between the respective haplotypes.The small black circle indicates the intermediate haplotypes that are not present in the sample.

Figure 5 .
Figure 5. TCS network of S. maximus based on COIII haplotypes.The sizes of the circles indicate the number of individuals with each given haplotype.Small lines illustrate the substitutions between the respective haplotypes.The small black circle indicates the intermediate haplotypes that are not present in the sample.

Table 1 .
Survey expeditions from which the presence-only (PO) data were compiled.

Table 2 .
Information on sampling stations for genetic analyses-location, geographic coordinates in decimal degrees (DD), depth, and date.

Table 3 .
Variables' relative percentage contribution and permutation importance to training data model.Values shown are averages over replicate runs.

Table 3 .
Variables' relative percentage contribution and permutation importance to training data model.Values shown are averages over replicate runs.

Table 4 .
Genetic polymorphism and neutrality tests of S. maximus inferred from the mitochondrial DNA COIII sequences.