A study of feasibility and detection sensitivity of environmental DNA for fish biodiversity monitoring and stock assessment in the Black Sea

Ivanova, Petya; Zlateva, Ivelina; Dzhembekova, Nina; Popov, Ivan; Raykov, Violin; Dimitrov, Dimitar; Mihova, Svetlana; Raev, Yordan; Stefanova, Kremena

doi:10.3389/fmars.2025.1648741

ORIGINAL RESEARCH article

Front. Mar. Sci., 26 November 2025

Sec. Marine Biology

Volume 12 - 2025 | https://doi.org/10.3389/fmars.2025.1648741

A study of feasibility and detection sensitivity of environmental DNA for fish biodiversity monitoring and stock assessment in the Black Sea

Institute of Oceanology, Bulgarian Academy of Sciences (IO-BAS), Varna, Bulgaria

Introduction: Human activities and their effects on ecosystems are threatening biodiversity and causing a decline in species, population, and genetic diversity. Environmental DNA (eDNA) metabarcoding represents a powerful instrument to assess marine biodiversity, enabling an efficient, cost-effective, and non-invasive approach. eDNA metabarcoding technology has significantly progressed in recent years and has been widely employed to assess marine fish populations.

Methods: This study assessed its potential, detection sensitivity, and effectiveness for an initial evaluation of fish biodiversity in the Black Sea. The efficacy of eDNA was evaluated by contrasting metabarcoding results with trawling data obtained from two sampling campaigns (summer and autumn 2022) across 16 selected locations along the Bulgarian Black Sea coast. Mitochondrial genes 12S and 16S were utilized as standard markers for metabarcoding given their high specificity and sensitivity in amplifying fish DNA, and the availability of large reference databases. A multi-model analytical approach utilizing Bayesian regression and Generalized Additive Models (GAMs) was implemented to assess the relationships between environmental DNA (eDNA) 12S metabarcoding data and abundance indices obtained from catch-per-unit-effort (CPUE). The method addresses significant challenges in predicting abundance from eDNA, including scarce, noisy datasets, excessive zeros, and nonlinear predictor-response relationships. Both frameworks reliably detected biologically meaningful associations between eDNA signal strength, environmental gradients, and trawl-derived abundance.

Results: Metabarcoding, employing the 12S mitochondrial gene with MiFish-U primers, exhibited superior sensitivity in identifying a broader array of fish species, highlighting its efficacy in detecting eDNA traces of rare and migratory species. During the autumn survey, environmental DNA (eDNA) analysis identified 23 fish species, compared to 15 species recorded through trawl sampling. In the summer expedition, 12 species were detected using molecular methods, while trawl surveys recorded 9 species. The Bayesian framework provided robust uncertainty quantification and reliable inference with few samples, whereas the GAM framework effectively captured nonlinear environmental and spatial effects. The results indicate that eDNA data can replicate patterns observed in conventional trawl surveys. This supports the assumption that eDNA data may serve as a low-impact, scalable method for monitoring fish populations. Model uncertainty informs the need for more sampling or a redesign of covariates to improve the accuracy of eDNA-based abundance predictions.

Discussion: Our findings demonstrate that the eDNA 12S primer combination exhibited superior taxonomic diversity compared to trawling, while necessitating reduced sampling effort. Furthermore, the eDNA-based approach can identify species across diverse life stages and sizes, unlike traditional trawling, which predominantly captures adult specimens within particular size categories. Although eDNA metabarcoding has limits in providing estimates for absolute biomass and size distributions, it provides significant insights on the presence of numerous fish species within the studied ecosystem. The consistency of results from multiple statistical frameworks supports the validity of eDNA counts as an abundance measure, provided that analytical models effectively account for dispersion and environmental variability. The Bayesian regression approach enabled thorough uncertainty quantification and dependable parameter inference with constrained sample sizes, whereas Generalized Additive Models (GAMs) provided an advanced representation of potentially nonlinear relationships among environmental DNA, environmental variables, and abundance.

Introduction

Marine biodiversity is rapidly declining as a result of human activities, including climate change, overfishing, the introduction of invasive species, and pollution. These factors expose a significant portion of marine life at risk of extinction (Brito-Morales et al., 2020; Yang et al., 2024). The Black Sea, a unique and highly vulnerable marine environment, is strongly affected by these challenges (Erdogan et al., 2009; Oguz, 2017). In recent years, selective commercial overfishing, amplified by illegal fishing practices and the use of non-selective fishing gears, has led to a rapid decline in the populations of economically important fish species, surpassing their regenerative capacity and resulting in a dramatic reduction in their overall abundance (Dai et al., 2022; Radulescu, 2023). This tendency represents significant threats to regional biodiversity economic stability, and human well-being. Hence, assessments of the marine living resources are essential for sustainable fisheries and ecosystem-based management (Trenkel et al., 2019).

The assessment of marine fish populations has traditionally depended on landings and catches, or on utilizing highly selective and invasive techniques that target commercially important species in specific geographic areas with favorable conditions. Consequently, it is essential to reevaluate current fisheries research methodologies and integrate complementary, minimally invasive techniques alongside traditional trawl surveys. This integrative strategy would reduce ecological consequences from sampling while improving spatial resolution for monitoring fish migration patterns and developmental stages and broaden community-level observation (Hammerl et al., 2024; Westgaard et al., 2024). Furthermore, the risk of species misidentification is significantly magnifying these challenges (Thomsen et al., 2012). Traditional species identification techniques are becoming increasingly limited due to dependence on the shrinking pool of skilled taxonomists having specialized phenotaxonomic knowledge (Hopkins and Freckleton, 2002). This highlights the pressing necessity for more precise, efficient, and cost-effective fish survey methodologies such as environmental DNA (eDNA) that will address current gaps in biodiversity assessment and improve fisheries management effectiveness.

Environmental DNA (eDNA) metabarcoding has emerged as a powerful tool for assessing marine biodiversity, allowing for the identification of many taxa within a single sample (Rourke et al., 2022; Tsuji et al., 2023; Haderlé et al., 2024). The data collected through eDNA metabarcoding contain complex and comprehensive biodiversity information. This method has the potential to become a cutting-edge instrument for biodiversity study and aquatic ecosystem conservation (Sahu et al., 2023). Recent advancements in eDNA technology have led to its increasing application for monitoring marine fish populations (Miya et al., 2015; Zhang et al., 2020; Afzali et al., 2021; Stoeckle et al., 2021; Valsecchi et al., 2021; Dai et al., 2022; Maiello et al., 2022; Zhang et al., 2022; Capurso et al., 2023; Cote et al., 2023; He et al., 2023; Osathanunkul and Suwannapoom, 2023; Veron et al., 2023; Westgaard et al., 2024), Despite its growing application globally, eDNA metabarcoding remains underutilized for assessing fish diversity and community structure in the Black Sea.

The application of eDNA to assess abundance in aquatic ecosystems is presently most advanced in freshwater and controlled environments (Rourke et al., 2022). Nonetheless, the complex interplay of external biological, physical, and environmental factors makes direct quantitative assessments in marine ecosystems highly challenging (Iversen et al., 2015; Hansen et al., 2018; Gilbey et al., 2021).

Using environmental DNA (eDNA) to quantitatively assess fish abundance in marine environments presents several dificulties. Key issues include establishing a robust relationship between eDNA concentrations and actual fish abundance, mitigating false positives and false negatives, elucidating the processes of eDNA degradation and dispersal in marine systems, and quantifying uncertainty in abundance estimates. Furthermore, the development of comprehensive DNA reference databases and the adoption of standardized protocols remain essential (Baetscher et al., 2025; Nneji et al., 2025). Despite this, significant evidence supports the use of eDNA as a supplementary tool for assessing fish population abundance and/or biomass across various spatio-temporal scales, dependent on preliminary studies to identify species- and context-specific factors influencing the eDNA abundance/biomass relationship (Salter et al., 2019; Rourke et al., 2022; Baetscher et al., 2025).

Clarifying the efficacy of a quantitative eDNA methodology and understanding its potential as a complementary instrument in traditional stock assessments would signify substantial progress in developing a cost-effective strategy for monitoring regional fisheries of both local and migratory species (Gilbey et al., 2021; Rourke et al., 2022). Further investigation is necessary to improve this prospective method as a viable tool for monitoring catch composition, assessing stock spatial distribution, and systematically recording changes in fish ecosystems across marine ecosystems (Rourke et al., 2022; Maiello et al., 2024).

The present study aims to advance the ongoing discourse on the application of eDNA methodologies in the Black Sea, with a particular focus on evaluating their effectiveness in:

a. improving fish biodiversity monitoring;

b. detecting and quantifying species presence;

c. providing robust scientific data to inform the conservation and sustainable management of marine fishery resources;

d. evaluating the feasibility of integrating eDNA monitoring into existing fisheries management frameworks in the Black Sea;

e. assessing the potential of eDNA metabarcoding counts to serve as a viable alternative to conventional trawl surveys;

f. support decision-making by engaging relevant stakeholders (e.g., policymakers, conservationists, and fisheries managers) in the adoption of eDNA as a complementary or alternative monitoring tool.

Materials and methods

Study area

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).

Figure 1

Figure 1. Map of sampling stations along the Bulgarian Black Sea coastline. The stations indicate the sites for eDNA sample that align with the endpoints of the corresponding trawling.

Trawl surveys

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 $L_{50} \approx 9.4 c m$ 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 $L_{m 50}$ (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.

Table 1

Table 1. Estimation of retention lengths for key commercial species using single trawl codend diamond mesh, size 63mm stretched, targeting $L_{50} \approx 9.4 cm$ and $S F \approx 3$ .

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.

eDNA sampling

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.

eDNA metabarcoding analyses

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).

Statistical analyses

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:

A b_{h} = β_{0} + \sum_{k = 1}^{p} β_{k} S e a s o n_{i k} + ϵ_{i j},

where: $A b_{h}$ is the Hellinger-transformed scaled abundance or sequence counts of species $j$ at site $i$ ; $E n v_{i k}$ represents the Season as categorical predictor throughout the sampling period and $ϵ_{i j}$ 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:

\log (\frac{p_{i, s}}{1 - p_{i, s}}) = β_{0} + β_{1} M_{i} + β_{2} S p_{i} + β_{3} (M_{i} * S p_{i}) + β_{4} T + β_{5} p H + β_{6} S + u_{S p_{i}} + u_{S i t e_{i}},

where: $p_{i, s}$ is the probability of detection of species $S p_{i}$ by sampling method $M_{i}$ (eDNA12S, eDNA 16s or Trawl survey), assuming that the probability of observing species detection by $M_{i}$ is a binary outcome with Bernoulli distribution: $y_{i} ~ B e r n o u l l i (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 $S p_{i}$ across different methods and covariates; $u_{S p_{i}}$ and $u_{S i t e_{i}}$ 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:

\begin{array}{l} B r i e r s c o r e = 1 / N \sum_{i = 1}^{N} {(p_{i} - y_{i})}^{2}, & () \end{array}

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

A U C = \int_{0}^{1} T P R (t) d (F P R (t))

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

A c c u r a c y = \frac{T P + T N}{T P + T N + F N + F P}

where: $T P$ and $T N$ are respectively the true positive and the true negative predictions/classifications, whereas $F N$ and $F P$ represent the false positive and false negative predictions/classifications;

T S S = T P R + T N R - 1 = T P R - F P R,

where: $T P R$ is the true positive rate (sensitivity) and $F P R$ 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.

Generalized additive models multi-model framework for best model selection

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.

Table 2

Table 2. Generalized additive models (GAMs) multi-model framework - summary.

Bayesian regression multi-model framework for best model selection

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).

Table 3

Table 3. Bayesian regression multi-model framework - summary.

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 $R_{a d j}^{2}$ (GAMs) and Bayesian $R_{a d j}^{2}$ 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.

Results

Fish diversity via eDNA metabarcoding and trawling

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).

Figure 2

Venn diagram with three circles labeled eDNA, 12S (blue), eDNA, 16S (red), and trawl (green). The overlapping regions contain numbers: 5 in 12S only, 1 in 16S only, 3 in trawl only, 3 shared between 12S and trawl, 0 between 16S and trawl, and 3 in the intersection of all three. Venn diagram showing shared and unique items among three categories: eDNA, 12S; eDNA, 16S; and trawl. The eDNA, 12S category has nine unique items, eDNA, 16S has none, and trawl has three. Two items are shared between eDNA, 12S and eDNA, 16S, five between eDNA, 12S and trawl, and none between eDNA, 16S and trawl. All three categories share seven items.

Figure 2. Venn diagram illustrating comparative analysis of species identified using eDNA and trawling in July (A) and November (B).

Table 4

Table 4. Species identified via eDNA (12S and 16S markers) and parallel trawl surveys throughout seasonal samplings in 2022.

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).

Biodiversity Indices

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.

Table 5

Table 5. Biodiversity indices for each sampling station computed using both the trawl surveys and eDNA (12S marker) data.

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

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.

Figure 3

Species accumulation curve graph showing the number of species on the vertical axis and sampling sites on the horizontal axis. The curve increases steeply and then levels off, with a shaded area indicating variability. Species accumulation curve for trawl data. The graph shows the number of species on the y-axis and sampling sites on the x-axis. The curve rises steeply, then levels off, with a shaded area representing variability.

Figure 3. Species accumulation curves. (A) computed for eDNA 12S marker and (B) for trawl survey as sampling effort increases.

Seasonality and species preferences

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.

Figure 4

Scatter plot illustrating multivariate data points grouped by season, with red circles for autumn and blue triangles for summer. The plot shows two distinctclusters along axes RDA1 and PC1, indicating seasonal variation.

Figure 4. Redundancy analysis (RDA) ordination plot showing the variation in species composition constrained by season.

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 ( $R^{2} = 0.7; R_{a d j}^{2} = 0.51$ ), demonstrating that the model effectively captured general relationships while accounting for the significant contribution of temperature and pH.

Sampling method detection sensitivity and species-specific detection sensitivity

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.

Figure 5

Graph titled “Posterior Predictive Check,” showing simulated density against observed data. The plot features curves representing observed data $y$ and replicated data $y_{rep}$, highlighting deviations. Most density lies between 0.00 and 0.25, and increases again beyond 0.75 on the x-axis.

Figure 5. Posterior predictive check for the detection model. Comparison of observed data (bold black line) with simulated posterior predictive distributions (grey density curves) to assess model fit.

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).

Figure 6

Estimated probability of detection by method, comparing eDNA and trawl surveys. Bars represent the mean estimated detection probability, with error bars indicating the 95% credible intervals.

Figure 6. Estimated probability of detection by method, comparing eDNA and trawl surveys. Bars represent the mean estimated detection probability, with error bars indicating the 95% credible intervals.

Figure 7

Estimated detection probability for each species, comparing eDNA 12S and 16S markers and trawl sampling methods. The points represent mean estimated detection probabilities, with error bars indicating 95% confidence intervals.

Figure 7. Estimated detection probability for each species, comparing eDNA 12S and 16S markers and trawl sampling methods. The points represent mean estimated detection probabilities, with error bars indicating 95% confidence intervals.

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 $p_{t r a w l d e t e c t i o n} > p_{e D N A 12 S d e t e c t i o n}$ equals 0, signifying that in all posterior samples the detection probability of eDNA 12S was greater than detection probability of trawl, and $p_{e D N A 16 S d e t e c t i o n} > p_{e D N A 12 S d e t e c t i o n}$ 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.

Table 6

Table 6. Effects of environmental covariates on sampling method detection sensitivity.

The reliability and accuracy of the model were assessed using the Brier score, AUC, TSS, and accuracy metrics, providing the following results: $B r i e r s c o r e = 0.105$ ; $A U C - R O C S c o r e = 0.901$ ; $a c c u r a c y = 0.862$ and $T S S = 0.6375$ . A $B r i e r s c o r e$ 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) $A U C - R O C S c o r e$ 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.

Modelling species abundance based on eDNA 12S sequence counts

Best models identified for the test case species through the application of the GAMs multi-model framework

Increasing the sample size enabled the GAM smoothers to identify functional relationships instead of being susceptible to noise overfitting (Tables 7, 8). A higher number of samples led to improved model stability, an increased proportion of deviance explained by the models, and demonstrated distinct, biologically credible relationships between eDNA 12S sequence counts and CPUE-derived abundance for multiple species. The results, indicated that selecting optimal distribution families is essential for accurately capturing the underlying data dispersion and zero prevalence among species. Cross-validation metrics further verify that predictive generalization, particularly rank-based consistency, may be achieved with moderate sample size, while small-n models were predominantly influenced by noise and sampling variance.

Table 7

Table 7. Best model, diagnostics (diagnostic plots – Supplementary Figures S4, S6, S8, S10; Observed vs. predicted - Supplementary Figures S5, S7, S9, S11) and model metrics for small sample size – $n = 16$ .

Table 8

Table 8. Best model, diagnostics (diagnostic plots – Supplementary Figures S12, S14, S16, S18; Observed vs. predicted - Supplementary Figures S13, S15, S17, S19) and model metrics for moderate sample size – n=100 (synthetic dataset).

Best models identified for the test case species through the application of the Bayesian regression multi-model framework

The results clearly indicated that the Bayesian regression multi-model framework effectively identified the relationships between eDNA 12S and CPUE-derived abundance for both small ( $n = 16$ ) and intermediate ( $n = 100$ ) sample sizes (Table 9). The Negative Binomial family with spline extensions demonstrated higher performance across all tested species, effectively capturing nonlinear eDNA-abundance trends while addressing overdispersion and zero inflation. M. merlangus exhibited robust correlations (r > 0.9) and moderate $R_{B a y e s i a n}^{2}$ (≈ 0.5), so confirming a reliable eDNA - CPUE-derived abundance relationship. The remaining species demonstrated greater variability in fits, perhaps due to biological and sampling limitations in the detectability of the eDNA 12S marker and/or trawl survey experiment design. The Bayesian regression approach demonstrated effectiveness with small, noisy datasets, yielding consistent parameter estimation and robust uncertainty quantification.

Table 9

Table 9. Best model, diagnostics (posterior predictive check plots – Supplementary Figures S20–S27) and model metrics for moderate sample size – $n = 100$ (synthetic dataset).

Discussion

DNA metabarcoding is a powerful tool for assessing biodiversity in environments that are challenging to explore with conventional approaches. This research employed simultaneous trawling and eDNA metabarcoding to compare trawl catches with eDNA data from corresponding water samples along the Bulgarian Black Sea coast. Consistent with other studies (Díaz-Ferguson and Moyer, 2014; Gilbey et al., 2021; Rishan et al., 2023; Suren et al., 2024; Çevik and Çevik, 2025; Sahu et al., 2025) our findings demonstrate that eDNA represents and excellent complementary technique, especially for improving detections capacities and providing a more comprehensive understanding of marine fish diversity. Our results demonstrated that eDNA reliably identifies a greater number of species across multiple habitat types. Collectively, our findings underline the value of eDNA as an additional approach that improves detection capabilities, strengthens biodiversity monitoring, and presents promising prospects for its incorporation into routine stock assessment frameworks.

Biodiversity

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.

Species-specific detection sensitivity and methods detection sensitivity

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).

Estimation of CPUE-derived abundance from the sequence read counts

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.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Ethics statement

Ethical approval was not required for the study involving animals in accordance with the local legislation and institutional requirements because Sampling studies were carried out within the scope of the permission given by the Executive Agency for Fisheries and Aquaculture, Ministry of Agriculture and Food, Bulgaria, according to Contract No D-156 of 13.03.2023. This study does not require Local Ethics Committee Approval as experimental animals were not used in this study.

Author contributions

PI: Writing – review & editing, Writing – original draft. IZ: Writing – review & editing, Writing – original draft. ND: Writing – review & editing, Writing – original draft. IP: Writing – review & editing, Visualization. VR: Writing – review & editing, Resources. DD: Writing – review & editing, Resources. SM: Resources, Writing – review & editing. YR: Writing – review & editing, Resources. KS: Writing – review & editing, Project administration.

Funding

The author(s) declare financial support was received for the research and/or publication of this article. This research was supported by the grants: European Union’s Horizon 2020 BRIDGE-BS project (Grant A. No. 101000240) and MASRI (project of the National roadmap for scientific Infrastructure (2020–2027) of Republic of Bulgaria).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2025.1648741/full#supplementary-material

References

Adamo M., Voyron S., Chialva M., Marmeisse R., and Girlanda M. (2020). Metabarcoding on both environmental DNA and RNA highlights differences betweenfungal communities sampled in different habitats. PloS One 15, e0244682. doi: 10.1371/journal.pone.0244682

PubMed Abstract | Crossref Full Text | Google Scholar

Afzali S. F., Bourdages H., Laporte M., Mérot C., Normandeau E., Audet C., et al. (2021). Comparing environmental metabarcoding and trawling survey of demersal fish communities in the Gulf of St. Lawrence, Canada. Environ. DNA 3, 22–42. doi: 10.1002/edn3.111