Predicting trawl catches using environmental DNA

Abstract

Quantifying the biomass, or number of individuals, diversity, and distribution of marine species is a critical aspect of understanding and managing marine ecosystems. In recent years, there has been growing interest in using environmental DNA (eDNA) for marine ecosystem management and biodiversity assessment. However, the main challenge hindering eDNA applicability has been the inability to infer absolute species abundances from multispecies analysis (eDNA metabarcoding). In this study, we demonstrate a way forward by estimating the abundance of commercially important fish species in a Norwegian fjord using a joint Bayesian statistical model of traditional trawl-catch data and molecular data derived from eDNA. Using this model, we accurately predict out-of-sample trawl catches using eDNA alone. Moreover, our model provides empirical estimates for key processes linking marine eDNA concentration to the fish population abundance estimated from trawl observations, including trawl catchability, DNA shedding, degradation, dilution, transport, recovery rate, and isolation efficiency. These processes, including amplification efficiencies correcting for Polymerase Chain Reaction (PCR) bias, are species-specific and enable the translation of eDNA metabarcoding data into abundances. These findings have broad implications for the use of eDNA in marine ecosystem management and conservation efforts.

Introduction

Quantifying the abundance, diversity, and distribution of species is a cornerstone of fisheries management and ecosystem conservation (Preston 1948, Jetz et al. 2012, Hill et al. 2014, Thomsen et al. 2016, Callaghan et al. 2021, Farr et al. 2022). Traditional survey methods, including scientific, and commercial trawl catches, close kin mark-recapture, telemetry, hydro-acoustics, and electrofishing surveys, have been widely used to gather ecological data and form the foundation for most fisheries management and conservation (Fraser et al. 2007, Neebling and Quist 2011, Crossin et al. 2017, Lees et al. 2021). However, these methods are invasive, and a few have a net negative ecological impact (Biju Kumar and Deepthi 2006, Eigaard et al. 2017). Furthermore, trawl and acoustic surveys are cost-intensive, requiring a large number of crew to sort the trawl catches, taxonomists onboard to identify species, and may be unsuitable for specific habitats, such as rocky bottoms (Lacoursière-Roussel et al. 2016).

Environmental DNA (eDNA), free genetic material shed by organisms which can be captured from environmental samples and sequenced, has emerged as a promising molecular approach that offers a less invasive method to identify and quantify species (Taberlet et al. 2012, Thomsen et al. 2012). Through “eDNA metabarcoding”—molecular identification of multiple species—researchers and resource managers can infer community composition, biodiversity changes, and ecological shifts due to anthropogenic impacts (Hansen et al. 2018, Jeunen et al. 2019, Atienza et al. 2020, Larson et al. 2022, Turon et al. 2022, Guri et al. 2023). Despite the recent technological advancements, eDNA methods have not yet reached the level of accuracy to fully replace traditional surveys for fish ecosystem monitoring or stock assessments. However, the integration of eDNA surveys with traditional methods can increase our understanding of marine ecosystems (Valdivia‐Carrillo et al. 2021, Pont et al. 2023), fisheries management (Stoeckle et al. 2021), and conservation (Hill et al. 2014, He et al. 2023). By combining these two survey types, we can obtain more comprehensive and accurate information on the species distribution and biological factors related to it (Veron et al. 2023). This becomes particularly pivotal when considering the ecological repercussions of trawl surveys, including adverse effects on non-target species and vulnerable benthic habitats (Lacoursière-Roussel et al. 2016). The integration of eDNA and trawl data can help to address some of the drawbacks of traditional survey methods, such as their high costs, invasive nature, inability to access restricted or vulnerable areas, and reliance on taxonomic expertise (Closek et al. 2019, Pont et al. 2023, He et al. 2023, Veron et al. 2023), thereby mitigating the dependence solely on the latter survey method (Valdivia‐Carrillo et al. 2021).

Despite the potential of eDNA as a method for understanding marine ecosystems, there are still several challenges that need to be addressed (Ramírez-Amaro et al. 2022). Metabarcoding methods do not reflect absolute concentrations of eDNA in environmental samples due to a combination of PCR stochasticity and species-specific amplification efficiencies, especially when using universal markers (12S, 16S, or COI). Furthermore, the compositional nature of metabarcoding data provides only inferable proportions (Gloor et al. 2017). Consequently, quantitatively assessing biomass or the abundance of targeted organisms through metabarcoding remains elusive (Rourke et al. 2022). Such constraints have hindered the broader adoption and integration of eDNA metabarcoding methods in cross-disciplinary fields such as fisheries, conservation biology, and ecosystem-based management (Ramírez-Amaro et al. 2022). Species-specific quantitative PCR methods (real-time PCR—qPCR or digital droplet PCR—ddPCR; Pont et al. 2023) offer an attractive alternative to metabarcoding because they are quantitative, hence reflect absolute eDNA concentrations. For example, Salter et al. (2019) showed a strong correlation (R² = 0.66; P = 0.008) between absolute eDNA quantities and the biomass of Atlantic cod, and Shelton et al. (2022) showed reliable hake distributions from eDNA on a continental scale. More recently, Maes et al. (2023) employed ddPCR and found eDNA concentrations can represent fish abundance and biomass for two commercial flatfishes (Solea solea and Pleuronectes platessa). In the context of ecosystem-based management of multiple species, however, species-specific quantitative methods can be more expensive and less efficient than metabarcoding multi-species, primarily due to the requirement for prior knowledge or isolation of the target species (Schneider et al. 2016).

eDNA metabarcoding has rapidly advanced as a tool for exploring fish communities. Several studies have shown positive relationships between eDNA metabarcoding data and trawl data (Thomsen et al. 2016, Afzali et al. 2021, Stoeckle et al. 2021, Yates et al. 2022). Stoeckle et al. (2021) compared log proportional metabarcoding reads to log proportional abundances from bottom trawls and observed a positive correlation when comparing across different months (mean R² = 0.565). Taken together, the available evidence indicates that metabarcoding reads without any type of transformation accounting for PCR biases cannot yield absolute abundance but may reveal information of rank, relative, or semiquantitative abundances if differences in amplification rates among species are small (Salter et al. 2019, Stoeckle et al. 2021, Guri et al. 2023, Veron et al. 2023). Pont et al. (2023) and Allan et al. (2023) addressed this issue by combining eDNA metabarcoding with qPCR (multiplying metabarcoding relative read-abundances with total qPCR fish DNA concentration) and found a strong correlation between total DNA concentrations and total fish biomass measured through traditional electrofishing in river systems. However, when comparing metabarcoding with the species-specific DNA concentration, they found roughly two orders of magnitude variation, which may have resulted from differences in species-specific amplification efficiencies. Alternatively, Shelton et al. (2023) developed a model to estimate initial DNA proportions in metabarcoding samples while accounting for species-specific amplification biases by sequencing mock community samples of known concentration of DNA extracts for a given list of taxa alongside environmental samples. To further improve the application of these models and methods, additional information is required to advance beyond proportions and provide absolute species concentrations from metabarcoding data. Addressing all aforementioned challenges will be critical to the effectively applying eDNA metabarcoding in ecological research and further integrating it into ecosystem-based management and fisheries stock assessments.

In this study, we developed a comprehensive and reliable framework for analyzing eDNA data in conjunction with other ecological data sources, such as trawl catch data, to improve the quantitative accuracy of metabarcoding analyses and provide a more holistic understanding of aquatic ecosystems. Our Bayesian joint model framework integrated eDNA concentrations from eDNA metabarcoding and ddPCR with fish density from trawl surveys (Fig. 1) to estimate of biological parameters of interest, such as abundance, and parameters describing the links between eDNA and traditional sampling data. Additionally, using these links, we show how we can predict out-of-sample trawl catches solely using eDNA observations.

Methods

We developed and joined five distinct datasets into a common analysis: trawl-catch counts, ddPCR observations of known concentration from standard samples for a reference species, ddPCR observations of unknown concentration from environmental samples for the reference species, metabarcoding of a known mixture of DNA extracts from multiple species of interest, and metabarcoding reads of environmental samples for target fish species. Below, we describe how each dataset was produced before providing detailed information on the joint statistical model and subsequent analyses.

In brief, the model uses species-specific ddPCR and a mock community to both calibrate metabarcoding results and to derive eDNA concentrations in environmental samples. We use this molecular dataset in combination with traditional trawl data to estimate the abundances of several species and rigorously test the resulting joint model against out-of-sample data to show its reliability.

Study area and samples

Balsfjord, northern Norway (Supplementary Fig. S1), is 40 km long and averages 150 m in depth with a sill at the fjord’s entrance, and the archipelago limits the water exchange between the fjord and the Norwegian Sea (open sea). Nearly all high-latitude Norwegian fjords, including Balsfjord, are ice-free and characterized by an Arctic light regime (Reigstad and Wassmann 1996). GPS coordinates and other metadata of sampling stations including GPS coordinates are provided in Supplementary Table S1.

Trawl surveys

Study samples were collected on research cruises in October 2019, 2020, and 2021 on the R/V Kristine Bonnevie as part of the Norwegian coastal annual surveys. The main aim of the survey was to collect data for abundance estimation of Norwegian coastal cod (Gadus morhua) and Northeast Arctic saithe (Pollachius virens) and catch count and weight of all fish species are recorded. Bottom trawl surveys were conducted annually in four distinct stations in Balsfjord alongside eDNA samples (Supplementary Fig. S1). Trawl surveys used a standard Campelen 1800 sampling trawl with an 80 mm (stretched) mesh size in the front section and 22 mm in the cod-end. The trawl sweeps were 40 m in length, and rockhopper gear was employed. Samples of the catch were sorted, weighed, and measured according to the recommended instructions (Mjanger et al. 2019) by the Institute of Marine Research crew.

eDNA water sampling, filtration, and extraction

Water was collected prior to trawling in three different 5-l Niskin bottles attached to a rosette which was raised and lowered using a deck-mounted winch at the established sampling station (Supplementary Fig. S1). At each station, we filtered triplicates of 2 l of seawater through 0.22 μm Sterivex filters (MerckMillipore) from three depths, surface (10 m), pycnocline (depth of highest density, ∼50 m) and bottom (10 m above bottom). Each filter was thereafter transferred to a sterile 50 ml Falcon centrifuge tubes and directly stored at −20°C. DNA extraction of water samples followed the manufacturer’s protocol using the DNeasy PowerWater Sterivex Kit (Qiagen GmbH, Hilden Germany), with minor adjustments (excluding PowerBead Tubes steps). For more details regarding these workflows see Guri et al. (2023). Negative controls such as field water and air blanks, laboratory blanks, and PCR blanks were taken throughout the workflow as described in Shu et al. (2020).

Droplet digital PCR

Quantitative molecular assays targeting the mitochondrial DNA (mtDNA), specifically 103-bp region of the ATPase gene of Atlantic cod (Gadus morhua, Table 1) and the cytochrome b sequence of Atlantic herring 86-bp (Clupea harengus, Table 1), were performed using ddPCR on seawater samples. Prior to running the environmental samples, assays were optimized using different primer/probe concentrations and temperature gradient to give the greatest fluorescence contrast between positive and negative droplets. Herring and cod ddPCR assays were run in duplex reactions using 6-FAM and VIC dedicated channels, respectively. All ddPCR runs were conducted using a QX200 ddPCR system (Bio-Rad) where ca. 20 000 droplets were generated in 20 µl reaction volumes. Each reaction used 5 µl of DNA template from environmental samples, 10 μl of ddPCR Supermix with no dUTP (Bio-Rad), 900 nM of final concentration of forward and reverse primers (Table 1), and 250 nM final concentration of TaqMan probe for each assay (cod and herring, respectively). The thermocycler reactions were run in C1000 Touch Thermal Cycler with 96-Deep Well Reaction Module (Bio-Rad) using the PCR program as follows: 10 min at 95°C for enzyme activation, followed by 44 cycles of denaturation for 1 min at 95°C and primer annealing and elongation for 2 min at 55.6°C, with a ramp rate of 2°C per s, followed by 10 min at 98°C and stored at 4°C. Alongside the environmental samples, we ran standard samples of known concentration from 10⁻²–10⁴ copies/µl for each species to construct the relationship between positive droplets and known concentration (Supplementary Table S2).

Mock communities

To calibrate metabarcoding observations and account for amplification bias (Gold et al. 2023, Shelton et al. 2023a), we selected ten fish species that were common in previous metabarcoding or trawl surveys in the study area (Guri et al. 2023). For the species Maurolicus muelleri (silvery lightfish), Gadus morhua (cod), Leptoclinus maculatus (daubed shanny), Hippoglossoides platessoides (long rough dab), Myoxocephalus scorpius (shorthorn sculpin), Cyclopterus lumpus (lumpsucker), Pholis gunnellus (rock gunnel), Brosme brosme (cusk), we used tissue samples collected by and DNA extracted at the University of Tromsø. The remaining species tissue samples, Mallotus villosus (capelin) and Pleuronectes platessa (European plaice), were purchased from the local fish market in Tromsø (Dragøy Fisk AS). Genomic DNA from bought fish muscle tissue samples was thereafter extracted using the DNeasy Blood & Tissue Kit (Qiagen GmbH, Hilden Germany) following the manufacturer’s protocol. A mock community of representative fish species was constructed from these genomic DNA samples by amplifying ca. 172 bp-long region of the 12S rRNA gene using MiFish-U universal primer set (Forward: 5′′-GTCGGTAAAACTCGTGCCAGC-3′′; Reverse: 5′′-CATAGTGGGGTATCTAATCCCAGTTTG-3′′; Miya et al. 2015). Reactions were run in 20 µl volume containing 12.5 µl of TaqMan Environmental Master Mix 2.0, 1.25 µl of each forward and reverse primers (10 nM concentration each), 8 µl of dH₂O and 2 µl of DNA template of above-listed species. The thermocycler program consisted of an initial polymerase activation for 10 min at 95°C, followed by 40 cycles of denaturation for 1 min at 95°C, annealing for 1 min at 60°C and primer extension for 1 min at 72°C, followed by a final extension for 10 min at 72°C and lastly storage at 4°C. PCR products were purified using the MinElute PCR Purification Kit (Qiagen GmbH, Hilden Germany) to ensure complete removal of primer-dimers. All samples were measured for concentration of dsDNA in purified PCR products using the High Sensitivity dsDNA assay with Qubit quantification system. PCR products were diluted in 10-fold steps to achieve a final concentration of ∼10⁵ copies/µl.

In total, six separate (but similar in concentration) mock community samples (250 µl each) were created by combining 25 µl of amplicon DNA from all 10 species in equimolar concentrations resulting in ca. 10⁴ copies/µl for each species (Supplementary Table S4). Mock community samples were amplified and sequenced alongside the environmental samples following the same workflow and protocol. Only six species in the mock community were also caught in trawl surveys (Supplementary Table S5) hence were selected for downstream analysis.

Library preparation, and NGS sequencing

In total, 109 samples divided into two libraries, including eDNA samples (n = 84), PCR blanks (n = 2), positive controls (n = 6), extraction blanks (n = 5), fieldwork water and air blank (n = 9 and 3, respectively), were amplified using the same primer set used for mock community samples (see above). We used fusion primers (primers containing adaptor and index), allowing for a one-step PCR protocol, and followed a thermocycler program and all following steps as described in Guri et al. (2023). Libraries were sequenced on GeneStudio S5 sequencer (Thermo Fisher Scientific) using the Ion 540-sequencing chip with the 200 bp protocol (Thermo Fisher Scientific).

Bioinformatics

Sequences were automatically demultiplexed and quality filtered after the sequencing process using Torrent Suite™ software inbuilt in the sequencer following their default settings. The sampled sequence dataset was thereafter filtered for chimeric sequences using a uchime-denovo algorithm in VSEARCH (Rognes et al. 2016). The sequences were clustered into Molecular Operational Taxonomic Units (MOTU) using SWARM 2 (Mahé et al. 2014) with a distance of d = 3, and subsequently the singletons were removed. MOTUs with a larger number of reads than 50 of the total library reads were retained and annotated by performing BLAST search against NCBI nucleotide (nt) database (23 August 2023) using the BLASTn algorithm. We used an E-value threshold of 1e⁻³⁰, maximum target sequences of 50 and a minimal percentage identity of 90. Lastly, we removed all taxa that were not assigned to fish (class: Actinopterygii or Chondrichthyes) and sequences whose read abundance was larger than 10% in negative controls (cumulatively) compared to the environmental samples.

Bayesian model

We developed a Bayesian joint model to synthesize information from trawl catches, ddPCR, and metabarcoding, to estimate fish quantities. We used mock community samples (box 3a in Fig. 1; only six species used) with known initial DNA concentrations to estimate species-specific amplification efficiencies. Knowing the amplification efficiencies enables the model to determine the species’ initial proportion and thus a species’ relative abundance in environmental samples of metabarcoding data (box 3b in Fig. 1; Shelton et al. 2023). Simultaneously, the model determines a link between known concentrations and positive droplets observed by the ddPCR system (box 2a in Fig. 1). Using this link-function, the ddPCR model estimates the absolute concentration of the reference (one) species arbitrarily chosen in ddPCR and mock community data in the environmental samples (box 2b in Fig. 1). Having both the relative abundance of species and the absolute concentration of one species in the environmental samples, the model expands proportion into absolute concentration for all species, or at least all species for which amplification efficiency can be estimated. Lastly, the model assumed a species-specific correlation between trawl catches and absolute DNA concentration which is a conversion parameter, and primary parameter of interest (θ) linking the two methods. See Fig. 1 for a simplified schematic representation of the joint Bayesian model.

In detail, in this model, we defined true fish density as N (units: fish/km²) and q for the catchability parameter on a scale of 0 to 1 (Zhang et al. 2020). Subscripts i, j, y, d, b, r, m, and p indicate species, station, year, depth category, biological replicate, technical replicate, observation in mock community aliquot samples, and observation in ddPCR standard aliquot samples, respectively.

The process (unmeasured biological parameters, Fig. 1)

Let X be the unobserved fish density detectable by trawling (fish/km² for species i in sample j in year y. Based on equations 1 and 3 in Mahévas et al. (2011), we model unobserved fish density (X) as a function of catchability (q) and true fish density (N):

Let C be the unobserved DNA concentration (copies/µl) from the same subscripts as X and let λ be the “integrated eDNA factor”—the conversion factor between the true fish density and unobserved DNA concentration. By analogy to catchability in the context of trawl data, we model the concentration of environmental samples as a function of the true fish density (N) and an integrated eDNA factor (λ) expressed as:

Note that we assume a simple proportional species-specific form for both catchability (q_i) and the integrated DNA factor (λ_i). These proportionalities are a common and reasonable assumption both in fisheries (Beare et al. 2005, Fraser et al. 2007, Mahévas et al. 2011) and environmental DNA (Stoeckle et al. 2021, Yates et al. 2022). However, different, and more complicated functional forms may certainly be reasonable and are worth further exploration (see more in discussion).

Biological parameters (Fig. 1)

Trawl model compartment (MC 1)

We model species-specific trawl catch count observations (Z_ijy) as negative binomial observations that are a function of annual mean fish density (V_iy) for each species (i) and year (y) across all the j stations (j = 1, 2, …, J_y), known fishing effort (E_jy, area-swept, km²; Fraser et al. 2007), and an overdispersion term φ_iy:

The annual geometric mean fish density (V) provides a link between the trawl observations (Z) and the unobserved fish density (X):

allowing X to be flexible at each station (j) but constrained for each year (across all the stations J_y; see more in discussion why this flexibility is beneficial).

ddPCR model compartment (MC 2)

Let S be the DNA concentration (units: copies/µl) from standards and C be the latent DNA concentration (Equation 1.2). The number of positive droplets (W) in ddPCR follows a binomial distribution (Guri et al. 2024, submitted), where the probability of a droplet being positive (ω) is a linear function of log_e DNA concentration (S or C) with an intercept (β0) and slope (β1) using clog-log (link-function), given the total number of droplets generated (U) expressed as:

Metabarcoding model compartment (MC 3)

Based on Shelton et al. (2023) we established that the number of observed reads (Y) from metabarcoding is a draw from a multinomial distribution given the proportions for each species ψ and the total number of reads per sample R:

Equations 4.2.1 and 4.2.2 link proportions ψ for each species i to observations, via a softmax transformation of the ratios of species abundances γ.

The known initial concentrations of all species i from mock community samples is converted into additive log ratios (alr_i) relative to a reference species of choice (G. morhua in this study, see Shelton et al. 2023), using the following equation:

where NPCR is the number of PCR cycles used to amplify the amplicons, α is the species-specific amplification efficiency (again, expressed as a log-ratio relative to that of a reference species; see Shelton et al. 2023), and η is the parameter that allows for overdispersion in the counts beyond the variability provided by the multinomial distribution, capturing the substantial variance among replicates often observed in metabarcoding data.

Having estimated the species-specific amplification efficiencies in Equation 4.3.1, we can derive the absolute DNA quantities for each species, station, and year (C_ijy), given the known concentration of the reference species (C_i=ref, G. morhua) estimated in Equation 3.2.2 for the same samples.

Note that because we have information about the absolute concentration of DNA from some species the term alr in Equation 4.3.1 can be replaced by the explicit ratio (difference in log-space) of DNA concentrations (eq. 4.3.2).

Joining of the model compartments

Using Equations 1.1 and 1.2, we can derive the relationship between unobserved fish density and the eDNA concentration as following:

Without additional information, we can estimate neither the catchability (q) nor the integrated eDNA factor in the marine environment (λ; these two parameters are unidentifiable). Therefore, we introduce a new parameter (θ) that encapsulates both the integrated eDNA factor and catchability,

For species detected by both trawl and eDNA we treat unobserved fish density (X) as a function of eDNA concentration (C) and the new parameter (θ), which is interpreted as the conversion factor between fish density observed by trawl toward eDNA concentrations (units: fish/km²/copies/µl):

By extension of Equation 5.2 and 5.3, where a species is observed in only one data stream, indicating no information on λ or q, θ inferences are uncertain (see the “Results” section, Fig. 3b and d). The integrated eDNA factor (λ) is a parameter that encapsulates multiple biological processes such as DNA shedding and degradation rate, DNA transport and dilution, and DNA isolation and extraction efficiency (θ, φ, ψ, and ξ parameters in Shelton et al. 2016).

The joint model (Supplementary Fig. S2) was implemented using the Stan language as implemented in R (package: Rstan) running four independent MCMC chains using 4000 warmup and 6000 sampling iterations (Table 2 for parameters and their priors). The posterior predictions were diagnosed using |$\hat{R}$| statistics (Gelman and Rubin 1992) and considered convergence for |$\hat{R}$| values less than 1.1 and effective sample size (ESS) greater than 1000 for all parameters.

To test the robustness of the parameters and the ability to use eDNA observations to estimate fish density, we conducted an out-of-sample analysis using two years of data as a training subset (part I) and the remaining year as a test subset (part II). In part I, the model estimated the parameters (θ, φ, β0, β1, α, η, ψ, γ, and τ) using observations (trawl, ddPCR, and metabarcoding) from only two years. Subsequently, for part II, we predicted fish density in the third year using the estimated parameters from part I in the trawl model (MC 1) and the eDNA model (MC 2 and 3). We defined delta (δ) as the difference between natural log fish density estimated by eDNA models and trawl models thus the log-fold difference for the out-of-sample year. Delta values closer to 0 indicate that trawl observations were well predicted using the eDNA model (or vice versa). This process was iterated three times separately with a different year data left out each time.

We included the seven species in the Bayesian model for which we had joint presence between trawl surveys or eDNA observations and ddPCR standard or mock community samples (Supplementary Table S5).

Results

The joint model successfully linked disparate datasets in a comprehensive and quantitative framework, yielding both meaningful estimates of fish abundances and derived parameters of substantial value. The posterior distribution of model parameters resulted in |$\hat{R}$| values <1.005 indicating that the chains converged to a common distribution, ensuring the stability of the outcomes. Additionally, the model had weak autocorrelation thus indicating efficient mixing as the effective sample sizes (ESS) were larger than ca. 2000. Posterior summaries of parameters together with |$\hat{R}$| and ESS can be found in Supplementary Table S6. The results regarding trawl catches, metabarcoding, ddPCR and mock community samples can be found in Supplementary Material S1.

Posterior predictive checks consistently indicated substantial agreement between predicted and observed data across all three observation methods: ddPCR, metabarcoding, and trawl surveys (Supplementary Fig. S3 a, b, and c, respectively). This result indicated the model’s capacity to accurately capture the characteristics of the observed data and their ability to extrapolate meaningful predictions forecasting future outcomes and deepening the understanding of the connections between these three observation methods.

The DNA concentration data generated by the eDNA models (ddPCR and metabarcoding) indicated a strong correlation with the unobserved fish density (X) on the trawl data (Fig. 2) given the species-specific conversion factor (θ). Such correlation indicated the reliability of the parameters (θ, β0, β1, α, η, and τ) in accurately representing the intricate relationships between the underlying processes and the observed data. This empirical evidence reinforced the dependency of the species-specific conversion factor (θ; Supplementary Table S7) for linking eDNA observation with other conventional methods (i.e. trawl catches) as the explanatory variable for underlying biological mechanisms between organisms and their shed and captured DNA.

The model output of the conversion factor (θ) ranged considerably from e^5.95 for C. harengus to e³⁰ for P. platessa (with θ fish/km² equalling 1 copy/µl observed in a sample extract; Table 3). The posterior distribution of the conversion factor differed significantly between the remaining species and varied between e⁸ and e¹⁴ fish/km² per copies/µl (Table 3). We noted that we observed no sequences of P. platessa in the metabarcoding data (Supplementary Table S7), thus the θ values for that species was estimated at the lowest boundary possible by the model. Similarly, due to extremely small catch counts, the conversion parameter of C. lumpus was hard for the model to estimate, thus a large standard deviation of the parameter's posterior distribution was observed. Rather than being reasonable estimates of θ, these two species instead represented lack-of-data scenarios in which insufficient information–in one case, no sequencing observations, and in the other, no trawl observations–results in unreasonable estimates. We included them specifically to illustrate some of the pitfalls and model prerequisites with respect to observation data.

The out-of-sample analysis showed that molecular data alone reliably predicted trawl catches (Fig. 3a and c). Five of seven species demonstrated a strong correlation between the predicted and observed fish densities, thus resulting in small δ-values (Fig. 3a and c). Among these, four species (C. harengus, G. morhua, L. maculatus, and M. villosus) resulted in δ-values of less than, or equal to, one for each predicted year, underscoring the robustness of the parameters governing the underlying relationships. Notably, H. platessoides displayed greater variability than the other four species, thus suggested slightly less robust and reproducible parameters among years. This variation may be attributed to an outlier year, although it was challenging to definitively ascertain this due to the limited dataset of n = 3 years in total (and therefore n = 2 years in training-data subsets). Increased availability of coupled data (trawl and eDNA) could enhance the robustness of θ and δ parameter estimation.

The remaining two species, C. lumpus and P. platessa, offered a useful illustration: when insufficient data were available, no meaningful inferences could be made. These species showed a notable average discrepancy of two orders of magnitude between the eDNA model-predicted fish densities and those derived from trawl observations (Fig. 3b and d). Furthermore, the mean within-year-out variability spanned over eight orders of magnitude, signifying the data insufficiency in either the trawl or the eDNA modules yielded challenges in estimating parameters.

Discussion

The current study has successfully bridged the gap between eDNA and trawl surveys, quantifying their bias alongside the biological relationships (jointly) between fish abundance and DNA concentration in the marine environment. Even more significantly, this bridge has been extended to multiple species, thereby enhancing the versatility and applicability of eDNA methods for stock assessments and marine biological monitoring. Extensive discussion regarding the outputs of metabarcoding, ddPCR and mock community can be found in Supplementary Material S2.

Biological interpretation of the model

Given the inextricable nature of DNA shedding and degradation, transport and dilution, and recovery and isolation processes, we unified them under the comprehensive term “integrated DNA factor–λ” in this study (framework in Shelton et al. 2016). We emphasize that the processes of catchability (q) and the integrated DNA factor (λ) are not assessed through this model but have been circumvented using the trawl to DNA conversion parameter (θ). Below, we discuss how these processes biologically interplay with the latent parameters estimated in this model.

The strong replicability of θ over years (for species with available data on both survey methods), as indicated by small values of δ in leave-one-out analysis, demonstrates that trawl metabarcoding and ddPCR observations yielded measurable and reliable results using the established model. Replicable results support the assumptions that trawl catchability (q) and the integrated DNA factor (λ) have remained constant throughout the years of survey period. Our rationale for assuming constant q is rooted in the standardized nature of the annual trawl surveys conducted by the IMR (Institute of Marine Research) and λ in similar environmental and biological factors (Zhang et al. 2020) during the sampling period.

The joint modelling of trawl and eDNA observations allows us to enhance the spatial resolution of fish densities. Following established research protocols, we model fish density as the annual average (V), based on four trawling catch count observations (assuming homogenous fish densities among stations), given the trawling effort. This approach becomes necessary because trawls at a specific location lack replication thus catches are commonly averaged over larger survey areas (Beare et al. 2005). Conversely, when we complement trawl observations using eDNA data, we gain the ability to develop models at the site level (X) and potentially future models at the depth level, rather than restricting eDNA analysis to the large aggregated, fjord level in our study. This flexibility offers a twofold advantage. First, it enables the ability to incorporate replicates of observation at different unit levels (site, area, station, and depth) thus increasing the accuracy of estimates. Second, it enables the model to establish the link between eDNA concentrations and fish densities without requiring a fixed relationship at every individual station, a factor that has been the downfall of previous studies attempting to establish this connection (Pont et al. 2023; Veron et al. 2023). This matches our understanding of the ecology of eDNA, in which DNA can be transported among stations (Rourke et al. 2022). This model advantage is thus particularly important as it mitigates issues related to lateral transport and dilution of eDNA (Andruszkiewicz et al. 2019), corresponding in maintaining the consistency of λ. Additionally, factors of eDNA recovery and isolation efficiency must have remained the same as the sampling technique and the laboratory workflow remained unchanged among the entire set of samples. Lastly, since all observations were conducted in October of different years, we posit that the behaviours of the fish population must have remained largely consistent, as did the relevant abiotic factors, both responsible for DNA shedding and degradation (Goldberg et al. 2011, Strickler et al. 2015, Jo et al. 2019, Mauvisseau et al. 2022). These results lead us to conclude that the integrated DNA factor (λ) has likely remained unchanged, particularly for species with small δ values.

For species with large δ values, several factors may explain variations of the θ parameter among years, which in turn could indicate either changes in q or λ between those years of observation. In the case of C. lumpus, the metabarcoding data suggest relatively small quantities, but trawl surveys observed minimal to no catches, suggesting small fish densities, or consequently pointing to an exceptionally small catchability term. This catchability term has been previously observed and indicated that bottom trawl surveys rarely catch C. lumpus, especially at depths greater than 60 m (Kennedy et al. 2016). Conversely pelagic trawls have been indicating larger catch rates of such species (Eriksen et al. 2014). Consequently, an unidentifiable catchability term can hinder the parametrization of θ. Conversely, flatfish species like H. platessoides and P. platessa exhibit strong catchability due to their preference for sandy and muddy bottoms. Trawl observations revealed relatively large fish densities for these species. However, in the case of P. platessa, no reads were observed through metabarcoding data, indicating an apparent absence of eDNA in the water, thus rendering the λ term ambiguous (hence also θ). An additional cause to this could be our filtering pipeline, as we removed amplicons with less than 50 reads, thus considered rare species (Jerde et al. 2016). Nevertheless, H. platessoides exhibited substantial read numbers and relatively large eDNA concentrations, but its δ values were substantial, indicating variation in the θ parameter among years of observation. This presents a more complex problem, and additional observations may be necessary to identify the source of this issue. We hypothesize that two potential errors could contribute to this problem. First, there may be an abnormality in the catchability terms in one of the years, making the outlier year unidentifiable given only three years of observations. Second, water samples instead of sediment samples could introduce biases for benthic species in the processes encompassed by the integrated DNA factor, such as DNA vertical transport and DNA recovery.

The variability among species highlights a general point that biological knowledge matters for the interpretation and use of eDNA observations, and that the use of eDNA for biomass estimation is likely species- and context-specific. Additional applications using eDNA information need to be developed to understand the kinds and magnitudes of uncertainty associated with using eDNA information in different contexts and with different species. Having a deeper understanding of these parameters and estimating their empirical value can enhance our biological understanding of the molecular sampling tools, thus boosting the applicability of eDNA to fisheries management and stock assessments (Jo et al. 2019). However, our work provides a roadmap for moving away from simple sample-to-sample correlations among sampling methods and toward a method that can inform and improve management and conservation actions. In this study, the specific conversion parameters are associated with the catchability of bottom trawling (the method’s inherent bias), and therefore, their applicability cannot be broadly extended to other methods or study areas. However, this statistical framework can be adapted to other sampling methods.

The future of eDNA and the trawl joint model

To our knowledge, this is the first empirical demonstration of species-specific conversion factors between eDNA metabarcoding and bottom trawl datasets. Future refinements should expand upon the existing model to include alternative sampling techniques and further generalize the framework and parameters we have presented here.

One potential area for improvement is the external estimation of trawl catchability (q; Fraser et al. 2007), which could lead to a more robust conversion parameter θ and also inform estimates of the integrated DNA factor (λ). Making λ a function of biological or environmental variables could allow for complex relationships between the trawl catchability and λ. However, exploring alternative forms is beyond the scope of this paper. One aspect of interest is understanding the effect sample depth, which could affect spatial- and species-specific DNA recovery and inform apparent differences between eDNA and trawl sampling (Jeunen et al. 2019, Canals et al. 2021). A further point of interest is accounting for the potential influence of allometric scaling on DNA shedding, metabolic rate, consumption, and excretion, and how potential interactions with environmental factors may affect the interpretation of eDNA survey results (Urban et al. 2023). An understanding of the effect of allometric scaling seems likely to be important for enabling molecular tools to estimate the size distribution of populations in marine environments (Yates et al. 2022). This could significantly further enhance the application of eDNA surveys into stock and population assessments.

The model can also be adapted to accommodate alternative monitoring techniques and data inputs. For example, catch-release, acoustic surveys, and trawl biomass data could either replace or be used in conjunction with the trawl model by adding another model compartment similarly to the presented model compartments (MC1, MC2, and MC3). Incorporating these additional methods would expand the range of available data collection techniques and provide further insight into the biological processes that underpin the molecular tools used in ecological assessments. Although this may alter the interpretation of parameters such as θ or λ, it would contribute to a more comprehensive understanding of their underlying biology. Our work shows how eDNA can inform and be included in the class of integrated analyses which are the foundation for modern stock assessments and fisheries management (Maunder and Punt 2013).

In summary, the future of this model holds significant potential for development and improvement. By meticulously documenting and mapping all parameters used in the model, along with the specific biological processes they represent, we can increase transparency and facilitate comparisons with other studies. This, in turn, enhances the replicability and reliability of surveys, a capability that conventional surveys alone cannot offer. Ultimately, these advances will contribute to the development of molecular sampling tools and their applicability in fisheries management and stock assessments.

Management applications

Trawl surveys provide consistent time-series data essential for stock and ecosystem assessments. However, the reliability and continuity of these time series can be influenced by various factors. Area accessibility (i.e. closed areas for conservation or energy transition), the necessity to explore new areas as fish stocks migrate, funding limitations, vessel availability, and adverse weather conditions can all impact the consistency of data collection. In response to these challenges, embracing innovative molecular tools and coupling them with conventional methods can provide more flexibility and resilience to changing circumstances while significantly reducing the ecological impact and maximizing the efficiency of surveys (Schneider et al. 2016, Di Muri et al. 2020, Veron et al. 2023). This strategic fusion can enhance the robustness and adaptability of survey programs, thereby strengthening our capacity for effective stock and ecosystem assessments in an ever-changing marine environment.

Conclusions

This study establishes the framework for species quantification and community composition of eDNA metabarcoding. We stress the need for species-specific conversion parameters to accurately estimate species abundances from eDNA metabarcoding. Furthermore, accounting for amplification efficiencies is necessary for ensuring the accuracy and replicability of metabarcoding data. These insights hold considerable implications for the application of eDNA in the management and conservation of marine ecosystems.

Acknowledgements

We thank all the technicians at IMR for the contribution during the field work and Agneta Hansen for the contribution during the field and lab work.

Author contributions

G.G., T.J., N.Y., J.I.W., R.K., J.R., A.S., T.H., J.F., and K.P. conceptualized the study; G.G., T.H., and J.I.W. conducted the fieldwork; G.G. and T.H. did the lab work; G.G., T.H., and J.I.W. did the bioinformatic analyses; G.G., A.S., and R.K. did the statistical analyses; G.G. with support from T.J., R.K., J.R., J.F., A.S., J.I.W., N.Y., and K.P. interpreted the statistical and ecological results; G.G. wrote the manuscript draft with contributions from all co-authors, and all authors approved the submitted version.

Conflict of interest

Authors declare no conflict of interest.

Funding

This research was supported by the project FISHDIV (NRC301691) funded by the Norwegian Research Council (NRC) and the program “Marine Ecosystem processes” at Institute of Marine Research.

Data availability

The data and the model’s script used in this article are available in Zenodo repository. (Dataset reference https://zenodo.org/records/12721619).

References