The study was conducted in the western Black Sea, within Bulgarian coastal and shelf waters between Cape Shabla (north) and Tsarevo (south). This region is ecologically and socio-economically important, influenced by seasonal hydrographic variability (summer stratification, autumn mixing) that shapes fish community composition. It supports diverse pelagic and demersal assemblages, and several locally rare or threatened species. The area also overlaps with national fisheries management zones and monitoring frameworks (Figure 1).

The trawl surveys were conducted under the National Fisheries Data Collection Program, which aims towards assessing the biomass of various species, including sprat, whiting, red mullet, and horse mackerel, with some species obtained as bycatch (Raykov et al., 2020). Sixteen randomly selected sites were surveyed to collect fisheries data, biological samples, and eDNA, with sampling carried out in July 2022 (summer) and November 2022 (autumn) as part of the National Data Collection Program (project № BG14MFOP001-3.003-0002-C02) (Figure 1; Supplementary Table S1).

A small-mesh bottom otter trawl designed for coastal and shelf surveys featured a four-panel body and a single codend of diamond mesh (63 mm stretched – bar ≈ 31.3 mm), and selection factor (SF) of approximately 3, targeting a L50≈9.4 cm was used during the surveys. This configuration ensured high retention for multi-species assemblages, while maintaining consistent selectivity baseline across species. Each tow lasted 30 to 40 minutes, with an average trawling speed of 2.7 knots, conducted solely during daylight hours according to the standard protocol, when target species were abundant in the deeper strata reaching the thermocline zone.

Each tow was treated as independent haul. The total catch was categorized to species through morphological identification, catch weight (kg)/numbers were recorded, length frequencies were processed to the nearest 0.5 cm (total length), and biological data (sex, maturity) were collected in laboratory conditions following standard protocols. The length-frequency data (LFD) obtained during the surveys were used to compute the implied retention thresholds for the four key commercial species (Table 1). Given the higher retention rate of the fishing gear the immature fractions as a percentage were computed based on the species-specifc logistic maturity ogives, parameterized by the species-specific Lm50 (the length at wich 50% of the population has reached sexual maturity). For each length class, the immature proportion [1−M(L)] ( M(L) represents the proportion of mature individuals) was applied to adjust the expected catch-at-length resulting in a catch-weighted immature fraction. To produce mature-only indices, we calculated the catch-per-unit-effort CPUE based on the standardized biomass (kg) and abundance (the latter derived from biomass using mean weight-at-length values for Black Sea fish species where available; Supplementary Table S2), and subsequently multiplied each by the maturity proportion M(L). The present study utilized abundance derived from maturity-adjusted CPUE.

Environmental covariates (temperature, salinity and pH) were logged at haul start and end.

Ethical approval was not required for this study, as no experimental animals were used.

Seawater was collected at each sampling station from three distinct depths within the water column: surface/near-surface (0–5 m), mid-depth (10–20 m), and deeper layers (30–50 m), using 5 L Limnos water samplers (KC Denmark A/S), with the exception of station 13S, where only two depth layers were sampled due shallow water depth. Environmental parameters (temperature, salinity and pH) were logged at each depth for all stations using internal thermometer of the water sampler, a digital refractometer (MA886, Milwaukee Electronics Kft), and a pH meter (HI83141, Hanna Instruments Ltd). To avoid contamination, the sampling device was thoroughly washed before each use. Two liters of seawater from each sampling depth were filtered using a 0.22 μm Millipore Sterivex filter with a peristaltic pump, labeled, and initially stored at −20°C onboard the research vessel, then transferred to −80°C in the laboratory until DNA extraction.

Genomic DNA was extracted from a total of 47 seawater samples (Sterivex Filter Units) using MO BIO PowerSoil^® DNA Isolation Kit (QIAGEN). DNA concentration was quantified with a Qubit 2.0 Fluorometer (Thermo Fisher Scientific). A set of 12S (MiFish-U) and 16S primers targeting fish eDNA was used for PCR amplification (Berry et al., 2019; Stoeckle et al., 2021). The preparation of paired-end libraries and high-throughput sequencing on the Illumina Miseq 300 PE platform (Illumina, San Diego, CA, USA) were performed by Macrogen Inc. (Seoul, South Korea).

Raw demultiplexed sequencing reads were quality-checked using FASTQC v0.12.1 (Andrews, 2010), and primers were removed with Cutadapt v4.0 (Martin, 2011). Sequence processing was carried out with DADA2 pipeline v1.32.0 (Callahan et al., 2016), including quality trimming, error rate modeling, merging of paired-end reads, chimera removal, and generation of Amplicon Sequence Variants (ASVs). In the filterAndTrim step, reads were truncated at 160 nucleotides (12S), and 200 nucleotides (16S). Sequences with exepected error rate (EE) higher than 2, or containing ambiguous bases were removed. For error rate modeling the default parameters were used with standard pooled processing for sample inference. Paired-end reads with fewer than 12 nucleotide overlaps were removed, with no mismatches allowed in the overlap region. Taxonomic assignment was performed with the BLAST + tool blastn: v2.16.0, using a Top Hit approach against custom-made 12S and 16S reference databases. These custom databases consist of complete or partial sequences obtained from the NCBI nucleotide database (https://www.ncbi.nlm.nih.gov), comprising 637 sequences (12S) and 769 sequences (16S), representing 117 and 121 fish species, respectively, known to occur in the Black Sea (Yankova et al., 2014).

Biodiversity indices, including species richness, Pielou evenness, Shannon index, Simpsons’ and its inverse and unbiased forms (Margalef, 1958; Moreno and Halffter, 2000; Spellerberg and Fedor, 2003; Colwell, 2009; Bowman and Hacker, 2020), were calculated to quanitify diversity patterns. Species accumulation curves were generated to illustrate the relationship between sampling effort and cumulative species detection (Moreno and Halffter, 2000; Bowman and Hacker, 2020).

To ascertain statistically significant correlations between eDNA sequence counts and trawl data, Spearman rank correlation (Simpson, 1949; Lehman, 2005) was applied. Additionally, Moran’s I was calculated to assess spatial autocorrelation prior to data modeling; nevertheless, no significant values were identified, hence incorporation of spatial lag was not considered feasible.

Community-environment relationships were analyzed using Redundancy Analysis (RDA) (Legendre and Gallagher, 2001; Borcard et al., 2011). A representative model for examining community dynamics in relation to seasonal trends was:

where: Abh is the Hellinger-transformed scaled abundance or sequence counts of species j at site i; Envik represents the Season as categorical predictor throughout the sampling period and ϵij indicates the residual error. Equivalent model formulations were employed for the examination of all continuous (including salinity, temperature, and pH) and categorical (sampling depth and season) environmental variables.

A Hierarchical Bayesian Logistic regression model was implemented using Stan, incorporating site and species-level random effects to account for uncertainty in species and site variation, along with environmental covariates, to evaluate the detection sensitivity of the sampling methods (eDNA and trawl survey) and the species-specific detection sensitivity for each method. The model can be expressed as follows:

where: pi,s is the probability of detection of species Spi by sampling method Mi (eDNA12S, eDNA 16s or Trawl survey), assuming that the probability of observing species detection by Mi is a binary outcome with Bernoulli distribution: yi~Bernoulli(p); β0 is the intercept; β1,2,3,4,5,6 are fixed effects which capture the differences in the probability of detection of species Spi across different methods and covariates; uSpi and uSitei represents random effect which accounts for species- and site-specific variability in detection.

Model performance was evaluated by partitioning the dataset into training (80%) and testing (20%) subsets. Prediction accuracy was assessed using Brier score (Brier, 1950; Stone, 1974; Yeung et al., 2005), overall model accuracy, the area under the ROC curve (AUC) (Fielding and Bell, 1997; Junge and Dettori, 2018), and true skill statistics (TSS) (Allouche et al., 2006), which are the most commonly used metrics for assessing classification model performance. The Brier score is a metric that evaluates the precision of probabilistic forecasts, ranging from 0 to 1, where 0 signifies perfect predictive accuracy and 1 denotes complete inaccuracy (Brier, 1950). AUC is a threshold-independent statistic that indicates a model’s ability to differentiate between positive and negative occurrences, whereas TSS is a threshold-dependent and prevalence-independent metric, providing it particularly useful for evaluating classifier performance in imbalanced datasets. The latter were computed with the subsequent formulas:

where: N denotes the total number of observations, pi represents the predicted probability for observation i and yi is the actual classification outcome (1 – detected; 0 – not detected);

where: TPR is the true positive rate (sensitivity), which represents the proportion of positive instances that were correctly classified by the model and FPR is the false positive rate, which is the proportion of positive instances that were incorrectly classified by the model. The ROC curve plots TPR against FPR across various probability thresholds t;

where: TP and TN are respectively the true positive and the true negative predictions/classifications, whereas FN and FP represent the false positive and false negative predictions/classifications;

where: TPR is the true positive rate (sensitivity) and FPR is the false positive rate.

Four commonly detected species in both eDNA and trawl survey expeditions (M. merlangus, M. barbatus, S. sprattus and S. maximus) were chosen as a test case to model the relationship between the maturity-adjusted CPUE-derived abundance and the eDNA 12s counts along with environmental covariates.

Given the limited size of the empirical data (only 16 observations per species (Supplementary Table S3), which is incufficient for robust modeling using hierarchical or multivariate approaches, intended for methodological testing, sensitivity analyses and comparability across various model families, a synthetic dataset ( n=100) was generated. The synthetic data was construced using the actual raw eDNA 12s sequence counts and CPUE-derived species abundance within a simulation framework, specifficaly designed to maintain the marginal distributions of each variable trough empirical or parametric models, while preserving the dependence structure among the variables and between species. The framework addressed the overdispersion and zero-inflation of eDNA data resulting from known causes, including stochasticity, insufficient sampling, and detection thresholds (Deiner et al., 2017; Peixoto et al., 2023). Consequently, eDNA observations were analyzed using a Zero-Inflated Negative Binomial (ZINB) distribution (Warton, 2005; Hilbe, 2011), which accounts for both detection failure and variability in biological counts. Species abundances were modeled conditionally on simulated eDNA sequence counts and environmental covariates using Generalized Additive Models (GAMs) (Wood, 2017), to capture nonlinear ecological responses. This hybrid ZINB-GAM model was intended to resemble the hierarchical relationship among eDNA generation, persistence, and organismal abundance as outlined in empirical case studies (Deiner et al., 2017; Yates et al., 2019; Harrison et al., 2019).

Environmental variables which have been reported to influence eDNA persistence and fish dispersal, such as temperature, pH, salinity, and site locations, were included in the analysis (Barnes and Turner, 2016; Goldberg et al., 2016; Harrison et al., 2019). Each covariate was treated as a continuous random variable, with empirical marginal distributions derived from observed data and subsequently resampled using quantile functions. The latter maintains the natural variability and environmental covariate ranges commonly observed in aquatic ecosystems. A Gaussian copula was employed to model the combined dependency between species’ eDNA and environmental variables, providing an adaptive multivariate framework that separates marginal distributions from their dependence structure (Nelsen, 2006). This is especially significant in ecology, where interactions among species as well as between species and their environment are nonlinear and non-Gaussian (Clark et al., 2017; Krupskii and Joe, 2013).

The result was a synthetic dataset including 100 observations, ensuring that all simulated variables maintain realistic interrelationships and conform to the actual covariance structure of the original sample data (Supplementary Figures S1-S3).

Various modeling methodologies were employed to determine their effectiveness in modeling CPUE-derived fish species abundance through eDNA 12S counts, alongside environmental covariates. These included Generalized Additive Models (GAM), Bayesian Generalized Linear/Aditive models (GLM/GAM), Bayesian hierarchical regression models, and ensemble machine learning (ML) models, the latter of which notably failed to achieve cross-validation due to insufficient data. The model selection inference focused on two robust multi-model approaches: a Generalized Additive Models (GAMs) Multi-model framework utilizing bootstrap blocks as a reliable alternative to cross-validation methods like Leave-One-Out Cross-Validation (LOOCV) due to the limited sample size, and Bayesian Regression Multi-model framework employing species-specific configurations.

A sampling size adaptive modeling approach, adequate for small and moderate datasets ( n<25 and n>25) was adopted, featuring 11 Generalized Additive Models (GAMs) with penalized regression splines, recognized for their ecological validity in analyzing overdispersed and continuous zero-rich data, such as abundance counts (Table 2). The latter was utilized to predict maturity-adjusted CPUE abundance for the test case species, integrating predictors including eDNA 12S sequence counts, environmental covariates (temperature, salinity, and pH), and spatial structure (geographic coordinates of the sampling sites), while applying a variety of response distributions. The framework was intended to facilitate the comparison of results and the selection of the most suitable model - “Best model”, given the expectation that the species may display diverse statistical distributions, including skewness, zero-inflation, over- and aunderdispersion, and nonlinear relationships with the response variable.

Furthermore, a Bayesian regression multi-model framework was employed, beginning with a simple hierarchical pipeline that included only eDNA sequence counts as predictors, focusing on a small dataset, and then expanding to include temperature in an attempt to conform with ecological plausibility. The latter evolved to the final version which is well suited to working with both small and moderate datasets, and provides distinct per-species models, allowing for species-specific flexibility and model comparison. A similar model selection approach to the GAM framework was adopted, featuring four possible models (Table 3).

A set of metrics was used to evaluate the resulting models, including Out-of-fold (OOF) Pearson correlation to assess the relationship between observed and predicted values (in the context of models’ ability to capture overall variability and accurately reproduce the linear trend in abundance), and OOF Spearman rank correlation, which provided an insight into how the model preserves sample ordering. Furthermore, the Root Mean Squared (RMSE) and Root Mean Squared Logarithmic Errors (RMSLE), as well as the Mean Absolute Error, were calculated to allow for effective comparison of predictive error magnitudes among models and to provide a metric for model predictive accuracy. The best model was chosen based on LOOIC (Leave-one-out Information Criterion), WAIC (Watanabe-Akaike Information Criterion) for the multi-model Bayesian framework, and RMSE or RMSLE values for best model selection in the multi-model GAM framework. ANOVA (analysis of deviance) was used to examine the significance of models’ terms and quantitative evaluation of deviance explained by the “Best model”, as well as the Radj2 (GAMs) and Bayesian Radj2 value for measuring Bayesian regression models’ overall signifficance.

All analyses and the accompanying graphical outputs were carried out in the R programming environment (R Core Team, 2021) utilizing R packages ‘brms’ (Bürkner, 2017a, 2017) with Stan (Carpenter et al., 2017), ‘rstan’ (Guo et al., 2020); ‘ggplot2’ (Wickham, 2016); ‘vegan’ (Oksanen et al., 2013), ‘mgcv’ (Wood and Wood, 2015; Wood, 2001), ‘copula’ (Hofert et al., 2014) and all related dependencies.

The eDNA analysis using the 12S marker demonstrated seasonal fluctuations in the number of detected fish species (demersal and pelagic), totaling 23 species. A higher number of species was detected in eDNA samples during autumn (23 species) compared to summer (12 species). Correspondingly, trawl catches revealed a greater diversity of fish species in autumn (15 species) than in summer (9 species), reflecting expected seasonal variations in species assemblages (Figure 2; Table 4).

Species diversity was more than twofold lower in both seasons when using the 16S primers, yielding only 9 species. The species identified using the 16S marker coincided with those detected using the 12S marker (Figure 2; Table 4). Furthermore, six species were identified in summer and twelve in autumn using both methods (eDNA with 12S and trawling) (Table 4). Species such as turbot (Scophthalmus maximus Linnaeus, 1758), whiting (Merlangius merlangus, Linnaeus, 1758), S. sprattus, M. barbatus, flounder (Platichthys flesus, Linnaeus, 1758), round goby (Neogobius melanostomus Pallas, 1814), and Alosa sp. were documented in trawl surveys during both summer and autumn, while migratory species like T. mediterraneus, anchovy (Engraulis encrasicolus Linnaeus, 1758), and bluefish (Pomatomus saltatrix Linnaeus, 1766) were found only in the specific seasons corresponding to their migratory patterns. The disparity in the autumn sampling survey was mainly attributable to eleven species (Alosa sp., transparent goby (Aphia minuta Risso, 1810), garfish Belone belone Linnaeus, 1761), leaping mullet (Chelon saliens Risso, 1810), little tunny (Euthynnus alletteratus Rafinesque, 1810), black goby (Gobius niger Linnaeus, 1758), sand goby (Pomatoschistus minutus Pallas, 1770), thornback ray (Raja clavata Linnaeus, 1758), chub mackerel (Scomber japonicus Houttuyn, 1782), Atlantic bluefin tuna (Thunnus thynnus Linnaeus, 1758) and Atlantic horse mackerel (Trachurus trachurus Linnaeus, 1758) that were not caught by trawling but were detected through eDNA. The last four species were likely missed by the trawl due to gear selectivity and mesh size. Conversely, four species: stellate sturgeon (Acipenser stellatus Pallas, 1771), pontic shad (Alosa immaculata Bennett, 1835), tub gurnard (Chelidonichthys lucerna Linnaeus, 1758) and Atlantic stargazer (Uranoscopus scaber Linnaeus, 1758) were captured by trawling but not detected via eDNA (Table 2). Additionally, the number of species detected in eDNA samples during summer differed by two species compared to the trawl data (Table 2). Rare species were detected exclusively in eDNA samples during the autumn survey at stations located in the southern region of the Bulgarian Black Sea coast: E. alletteratus (33A), S. japonicus (27A), T. thynnus (27A), and T. trachurus (16A, 18A, 25A, 27A, 32A).

The highest species richness values (Table 5) were predominantly recorded at most stations sampled by eDNA12S during the autumn expedition, with the exception of station 16A, which nonetheless exhibited a species count exceeding the average determined through trawling. The greatest number of species (12), obtained through trawling, was observed at station 18A during the autumn sampling season.

Pielou’s evenness ranging from 0 to 1, with values nearing 0 was noted at three stations: trawl survey hauls gathered at stations 13S, 16A, 32A attributed to the significantly large captures of S. sprattus (13S) and T. mediterraneus (16A and 32A), which predominated the catch composition.

The Shannon index reached its maximum values during the trawl survey at station 18A and in the eDNA (12S marker) assessment at stations 25S, 31S, 33S, 16A, and 33A. The first three eDNA samples were obtained during the summer expedition, whereas stations 16A and 33A were sampled in the autumn expedition, revealing heightened biodiversity indices in both seasons, likely due to particular developmental stages (in summer) and/or migratory species (in autumn) whose pathways include the coastal and shelf waters of the Bulgarian Black Sea.

Species accumulation curves were used to illustrate the increase in species count relative to sampling effort. The eDNA species accumulation curve (Figure 3A) indicated that roughly 10 more species could be detected using eDNA (12S marker) compared to those identified through trawl surveys considering the sampling effort (Figure 3B). Thus, the extensive reach and non-invasive characteristics of eDNA analysis enable the identification of a greater variety of species within the environment.

Seasonal patterns in species communities, detected by both sampling techniques, were analyzed using constrained ordination (Redundancy Analysis - RDA). Seasonal clustering, as indicated by confidence ellipses, was distinctly observable for both seasons (Figure 4), presumably attributable to species habitat suitability and/or particular biotic interactions relevant to each season. Only a limited number of species within the detected community structure during the autumn season showed preferences for an abiotic environment comparable to that of the summer season, possibly due to species interactions or early recruitment developmental stages.

The seasonal pattern explained over 50% of the variance in community structure (based on abundance/eDNA 12S sequence count data). Identical clustering patterns were identified for species associations constrained to environmental covariates (temperature, salinity, and pH), which also show seasonal cycles. The RDA, that includes all continuous and categorical variables, accounted for 65% of the variance in species CPUE (kg)/sequence counts ( R2=0.7; Radj2=0.51), demonstrating that the model effectively captured general relationships while accounting for the significant contribution of temperature and pH.

The posterior predictive check for the Bayesian Hierarchical model, utilized to evaluate method detection sensitivity and species-specific detection sensitivity per method (Figure 5), demonstrated a strong agreement between the actual data and the simulated posterior predictive distribution. The posterior predictive check was conducted utilizing new observations unknown to the model (train-test split 80–20%), indicating that the model sufficiently captured the data variability.

The model results highlight the influence of sampling methods and environmental covariates on method detection sensitivity (Figure 6) and determine the probability of detection for each species, utilizing the sampling method as a predictor while considering the effects of environmental covariates (Figure 7).

The detection probability for each method was obtained from the posterior distribution, based on the null hypothesis that trawling will identify more species than eDNA. The results indicated the chance of ptrawl detection>peDNA 12S detection equals 0, signifying that in all posterior samples the detection probability of eDNA 12S was greater than detection probability of trawl, and peDNA 16S detection>peDNA 12S detection equals 0, signifying that in all posterior samples the detection probability of eDNA 12S was greater than detection probability of eDNA 16S, hence validating the advantages of 12S primers over the eDNA 16S primers and trawl sampling for species detection.

The estimated species-specific detection probability evidently differs among various species; those that are naturally more abundant and prevalent in the environment (such as M. merlangus, S. maximus, and S. sprattus) and those with seasonal occurrences (E. encrasicolus, P. saltatrix, and T. mediterraneus) have a greater probability of detection. Furthermore, these species exhibit elevated detection probability across all sampling methodologies, despite the natural selectivity of trawling, which targets particular species based on trawl design and fishery type.

Rare species exhibit low detection probabilities, evident at the detection threshold, and are likely to be identified solely through eDNA (12S), given the constraints of trawling as a sampling technique and the lower sensitivity of 16S primers, which captures only the most abundant species and exhibited inability of detection of rare species and such present in low abundance.

The model findings revealed a statistically significant negative effect of temperature (Table 6) on method sensitivity, suggesting that high temperatures may adversely affect eDNA detection sensitivity, with no discernible influence from salinity and negligible impact from the pH.

The reliability and accuracy of the model were assessed using the Brier score, AUC, TSS, and accuracy metrics, providing the following results: Brier score=0.105; AUC−ROC Score = 0.901; accuracy=0.862 and TSS=0.6375. A Brier score of 0.105 signifies that the model’s probability estimates are well-calibrated, indicating accurate reflection of species detection likelihood (Brier, 1950; Trick et al., 2023) AUC−ROC Score  of 0.8912 suggests that the model has very good discriminatory capacity, meaning that it correctly ranks both detected and undetected instances; An accuracy of 86% is strong, indicating that 86% of the model’s predictions are correct.

Our findings indicate that the eDNA 12S primer pair revealed greater taxonomic diversity than trawling, while requiring less sampling effort. Moreover, the eDNA-based method may detect organisms at various life stages and sizes, in contrast to conventional trawling, which selectively capture mature individuals within specific size ranges (Afzali et al., 2021). Although eDNA metabarcoding cannot provide estimates of absolute biomass, developmental stage, and size distribution, it offers substantial insights into the presence of numerous fish species within the studied ecosystem.

Another advantage of eDNA metabarcoding is its efficiency in a broader functional spectrum compared to bottom trawling and is capable of detecting both demersal and pelagic species present in the water column (Veron et al., 2023). Along the Bulgarian Black Sea coast, eDNA metabarcoding (12S) detected 13 pelagic species, including four rare species (T. trachurus, S. japonicas, E.s alleterattus, and T. thynnus), eight migratory species, and twelve demersal species (Table 4). Nevertheless, some case studies indicate that pelagic species are captured during demersal surveys in the Black Sea, as summer stratification influences pelagic species such as S. sprattus and T. mediterraneus which remain beneath the seasonal thermocline, and only ascend during the summer-autumn mixing (Tserkova, 2013; Panayotova and Todorova, 2015; European Commission, Joint Research Centre (JRC), 2025). The latter is not relevant for juvenile fish wich tipically inhabits shallower waters. Accordingly, it is reasonable for the region to expect co-occurrence of both pelagic and demersal species in the deeper layers (25-100m). The rare species, although not captured by trawling, are known to occur in the Black Sea (Yankova et al., 2014). Overall, eDNA metabarcoding significantly enhanced detection rates of rare and migratory species across diverse habitats without requiring physical capture or direct observation. This capability is especially beneficial for monitoring species that are difficult to detect using traditional methods.

Biodiversity indices revealed substantial seasonal fluctuations in species communities throughout the Bulgarian Black Sea coastline. Both methods revealed seasonal clustering patterns in fish communities (Table 4), with autumn exhibiting greater biodiversity and species richness compared to summer. eDNA detected higher number of species compared to trawling, including rare and those that are less prevalent during specific seasons, such as E. encrasicolus and M. barbatus in July, and S. sprattus in November. Accordingly, eDNA metabarcoding may represent a viable instrument in monitoring of shifts in species distribution and behavior associated with climate change, and assist conservation efforts by identifying locations of the core-density areas for key-species, such as S. maximus, S. sprattus, and M. barbatus for targeted conservation management. These new discoveries improve our understanding of migration timing and mechanisms, aiding the formulation of successful conservation measures.

Only six species were identified to overlap when employing both sampling approaches (trawl and eDNA 12S and 16S). Eleven species were detected only with eDNA, while four species (Table 4) were exclusively identified through trawling. Additionaly one species was identified on the genus level (Alosa sp.). The efficacy of eDNA-based biodiversity assessments depends on the thoroughness and precision of reference databases relevant to local biological populations, as emphasized in the study of Zhang et al. (2020). The lack of reference data might hinder taxonomic assignment, resulting in ambiguous species identification or, in severe instances, misidentification due to dependence on ‘best match’ algorithms (Gold et al., 2022). To address this, our study used custom reference databases, that included species known to inhabit the Black Sea, however, sequences for many Black Sea species are absent from NCBI and other public repositories, underscoring the need for greater efforts to incorporate additional species-specific mitochondrial DNA sequences into open access databases (Kasmi et al., 2023; Kumar et al., 2022; Xiong et al., 2022). Another potential explanation for the non-detection of some species via eDNA could be their rare occurrence or low abundance in the survey area, resulting in insufficient detectable eDNA. Species detected exclusively through trawling (Alosa immaculata, Chelidonichthys lucerna, Uranoscopus scaber, and Acipencer stellatus), exhibited minimal biomass in the catches (individual catch ≤ 4 kg). Furthermore, low-abundance taxa are often under-dected by traditional survey gear because encounter rates and capture probability decline with density (Williams et al., 2019; Adamo et al., 2020).

Scientific trawling remains the conventional method for assessing marine ecosystems and fish populations along the Bulgarian Black Sea coast, while eDNA-based species identification signifies an emerging scientific field and one of the most rapidly evolving tools in fisheries science and management (Lacoursière-Roussel et al., 2016; Casas Castaño and Saborido-Rey, 2022). The efficacy of fisheries management and conservation strategies depends heavily on the accurate identification of exploited species. This issue is especially challenging when dealing with cryptic species, larval stages, or morphologically similar bycatch (Casas Castaño and Saborido-Rey, 2022; Hansen, 2019). Accurate species-level identification is crucial for ensuring accurate catch records while supporting sustainable fisheries management. An illustrative instance of a challenging scenario involves the Ponto-Caspian species of the genus Alosa, which are assigned as vulnerable on the IUCN Red List. Morphological and genetic variations illustrate the taxonomic complexity of Alosa species in the Black Sea (Dobrovolov et al., 2012). The inconsistencies may result from misidentification in catches, as shads exhibit significant morphological plasticity, with overlapping morphological traits among species (Orlova et al., 2024). Genetic markers are essential for acquiring accurate information regarding species status and enhancing conservation efforts. Given these identification issues, some Alosa specimens gathered during our concurrent trawl research were identified as Alosa immaculata, while others could not be identified to the species level and were reported as Alosa sp. Similarly, environmental DNA (eDNA) metabarcoding at the same sampling locations (33S, 18A, and 25A), as well as an additional site (27A), using solely molecular approaches, also failed to resolve determination to species level. The absence of genetic differentiation within the genus Alosa, based on various mitochondrial markes (Antognazza et al., 2019; Orlova et al., 2024; Bhendarkar et al., 2025), including those used in the present study, hampered species-level identification. Thus, regardless of the fact the ASVs were affiliated with A. falax (the only Alosa species represented in the reference library used) in the present study, the results could only be interpreted at the genus level, and intragenus differentiation remains highly challenging.

Our findings indicated that species that are naturally more common and prevalent in the ecosystem have a higher likelihood of detection. Moreover, these species demonstrate increased detection probability across all sampling methods, including the trawling, which is selective and targets specific species according to trawl design and fisheries type. Moreover, species with high abundance are anticipated to exhibit a greater probability of correlating eDNA sequence counts with biomass or abundance, as their eDNA signal is expected to be more pronounced due to their naturally higher densities in the environment, thus increasing the likelihood of detection by both conventional and eDNA methodologies. Higher abundance implies higher DNA shedding rates (Tillotson et al., 2018), and effective eDNA sampling with species-specific mitochondrial markers (12S rRNA, COI) can provide quality data to infer actual abundance and spatial distribution.

In contrast, rare species demonstrated low detection rates, apparent at the detection threshold, and are likely to be discovered only via the 12S eDNA marker, due to the limitations of trawling as a sampling method and the reduced sensitivity of the 16S primers.

The concentration of eDNA in marine ecosystems may vary due to numerous abiotic factors influencing its dispersal and degradation rates, including currents, temperature, and pH (Collins et al., 2018; Andruszkiewicz et al., 2019; Fukaya et al., 2021; McCartin et al., 2022). Our data revealed statistically significant adverse impacts of temperature on method sensitivity. Elevated temperatures, particularly in shallow coastal waters, may induce behavioral avoidance in fish, compelling them to seek cooler environments (Pörtner and Peck, 2010; Perry et al., 2005), and may also lead to reduced eDNA persistence due to accelerated degradation (Strickler et al., 2015; Barnes and Turner, 2016). The latter emphasizes the importance of considering environmental fluctuations when evaluating methods sensitivity, and particularly when inferring abundance from eDNA sequence counts (Parsley et al., 2024; Jo, 2025).

There is increasing interest in employing eDNA sequence counts to assess actual fish abundance and biomass (Dukan et al., 2024). Many studies have indicated substantial correlations between eDNA data and fish abundance/biomass in marine ecosystems (Thomsen et al., 2016; Salter et al., 2019; Afzali et al., 2021; Stoeckle et al., 2021; Rourke et al., 2022; Russo et al., 2021). Several research studies have found that using specific statistical methodologies, including Generalized Linear Models (GLM), Generalized Additive Models (GAM), non-linear models, or machine learning, can enhance the precise modeling of species abundance based on sequence counts (Collins et al., 2022; Fowler et al., 2024; Kasmi et al., 2025; Ip et al., 2025). However, some studies have demonstrated a lack of quantitative association between biomass/abundance and metabarcoding sequences, coupled with significant uncertainty (Lamb et al., 2019; Knudsen et al., 2019; Fraija-Fernández et al., 2020).

This research presents a robust multi-model analytical approach providing two multi-model frameworks implementing Bayesian hierarchical regression and frequentist Generalized Additive Models (GAMs) to establish relationship between eDNA 12S metabarcoding data and CPUE-derived abundance indices. The method addresses the primary challenges of eDNA–abundance modeling including constrained and noisy datasets, excessive zeros, and the potential for nonlinear environmental effects. In all frameworks, the models consistently identified biologically significant correlations between the abundance counts environment, spatial structure and/or eDNA signal strength. For Merlangius merlangus, both the Tweedie GAM and Bayesian Negative Binomial-spline models exhibited significant explanatory power. However, the GAM-derived model failed to account for any explanatory power of eDNA counts in relation to the predicted abundance, and the fit was significantly affected by the spatial structure.

Nonetheless, the uniformity of results from diverse statistical frameworks enhances the credibility of eDNA counts as a measure of abundance, assuming that analytical models adequately address dispersion and environmental variability. The Bayesian regression method facilitated comprehensive uncertainty quantification and reliable parameter inference in limited sample sizes, while Generalized Additive Models (GAMs) offered a sophisticated depiction of possibly nonlinear relationships among environmental DNA, environmental variables, and abundance. The latter illustrates that this complex approach enables eDNA data to replicate the spatial and environmental trends identified in trawl surveys, so facilitating the monitoring of fish populations with less environmental impact and increased frequency.

Our findings revealed that, depending on the strength of species-specific abundance and eDNA signals, many taxa (such as M. barbatus) had limited or absent relationships between observed and predicted abundance. This underscores the need for integrating additional molecular, ecological, and physical variables in abundance estimation to improve model explanatory power. A significant problem is the necessity for a carefully designed targeted sampling approach in both eDNA and trawl surveys to accurately capture actual species-specific signals, particularly when species coexist, leading to mixed catches. Concerning the latter, it is essential to highlight that the outcomes provided by the Bayesian regression framework can enable the direct translation of model uncertainty into recommendations for future eDNA survey design. The posterior variance and cross-model agreement identify species for which specific variables strongly affect predictive uncertainty, emphasizing areas where more sampling would improve model accuracy. Species displaying broad credible intervals or unstable correlations may require further spatial replication or more accurate temporal coverage, particularly for species that are seasonally present in the sampling area and have distinct migratory patterns. In contrast, strong relationships (such as those evident in M. merlangus model) suggest that eDNA signal strength reliably reflects abundance based on the present effort.

Another observation concerns the integration of environmental variables, including temperature and salinity, which consistently improved fit stability. This highlights the need of jointly predicting eDNA concentration and habitat context. The latter demonstrate that a probabilistic, multi-model Bayesian framework can improve adaptive sampling and concentrate sampling efforts, hence strengthening confidence in eDNA-based monitoring especially through application of quantitative markers as an alternative or supplement to traditional trawl surveys.

However, as with many emerging technologies, eDNA faces several limitations. The accuracy of reference databases, environmental heterogeneity, and DNA degradation rates are critical factors that can compromise the reliability of eDNA-based analyses. Establishing standardized protocols for data acquisition, sampling, and interpretation remains essential (Çevik and Çevik, 2025).

Improvements remain feasible through the application of internal DNA standards in 12S metabarcoding (Stoeckle et al., 2021) and qMiSeq approach (Tsuji et al., 2022), combined with an understanding of species-specific DNA shedding rates, eDNA transport dynamics, and signal degradation. The latter may undoubtedly provide a more robust basis for the precise assessment of biomass in less abundant species, particularly uncommon, newly introduced, and invasive ones. Moreover, several challenges persist, including the need for reliable data on DNA shedding and degradation rates, species mobility in the water column, hydrodynamic conditions influencing DNA transport and decay, and the necessity for precise ecological and abundance data to validate the model outcomes. Furthermore, the integration of eDNA data with ecological models and traditional approaches is essential to improve precision and deepen the understanding of ecological dynamics. By overcoming these constraints, eDNA could emerge as an essential instrument for biodiversity evaluation and conservation efforts (Çevik and Çevik, 2025). In conclusion, current advancements in molecular biomonitoring have expanded opportunities for aquatic ecosystem assessment (Blancher et al., 2022), and its application is becoming increasingly widespread (Ruppert et al., 2019). eDNA metabarcoding is an evolving method that is likely to continue developing as technology advances and standardization progresses. We conclude that it is a promising instrument that can complement and substitute the scientific trawling for biodiversity monitoring and that it can be successfully adapted for abundance estimation, in agreement with Veron et al. (2023). Further advancements are needed before this method can be widely used for fisheries assessments in the Black Sea basin, with a reasonable initial emphasis on the most abundant species. Future biomonitoring programs in the region should integrate eDNA metabarcoding with conventional surveys to yield a more comprehensive and precise evaluation of ichthyological diversity and fish stock, hence facilitating evidence-based policy formulation.