Seawater samples were collected aboard the R/V Miguel Oliver during the PELACUS 0324 survey (27 March–16 April 2024), which is part of a historic series of annual surveys dedicated to the assessment of fish stocks in continental shelf waters off the northern and northwestern Spanish coast. During this survey, seawater samples were collected at three sampling stations in the northwest of the Iberian Peninsula (Fig. 1): Northern Galicia (Station 1; 43°48.08'N, 7°52.61'W), Costa da Morte (Station 2; 43°05.29'N, 9°20.86'W), and Ria de Arousa (Station 3; 42°30.98'N, 8°57.45'W).

A total of 8 L of seawater was collected at each sampling station using a continuous circuit intake at a depth of 4–5 m. During sample collection, seawater was pre-filtered using a 60 µm mesh to remove larger sediments that could clog the filters. The water collected was divided into two replicates of 4 L each. Subsequently, the entire volume of each replicate was filtered through a cartridge filter (Sterivex-HV Durapore PVDF, 0.45 µm pore size) using a peristaltic pump (Pegasus Alexis XL V4.0). Filtration control samples were obtained by filtering 2 L of Milli-Q water for each sampling station. All materials used were previously cleaned by soaking them for 15 min in a 10% bleach solution, followed by several rinses with Milli-Q water. After filtration, the samples were stored at –20 °C until extraction in the laboratory. Then, eDNA was extracted using two different methods. The DNeasy PowerWater Sterivex Kit (Qiagen) was used according to the manufacturer’s instructions for Method 1. For Method 2, eDNA was extracted using a version of the DNeasy PowerWater Kit (Qiagen) protocol modified in this study to include an initial homogenization step to test whether breaking filtered cells prior to the employment of the standard DNeasy PowerWater protocol could result in greater uptake of vertebrate DNA. This involved removing the filter from the cartridge, cutting it into small fragments with a sterilized scalpel, placing the pieces into an additional Eppendorf tube along with a metal bead, and adding 1 mL of PW1 solution. The samples were then homogenized using a Precellys Evolution homogenizer (Bertin Technologies) according to the following protocol: two cycles of 30 s at 6,800 rotations per minute. The liquid was then transferred to PowerBead Tubes, and the manufacturer’s protocol was followed as indicated from that step onward. The final elution volume of extracted DNA was set at 100 µL for all samples.

Recovered eDNA was quantified using the Qubit™ dsDNA High Sensitivity (HS) Assay Kit (Thermo Fisher Scientific). One blank extraction control was used for each tested method (Method 1 and Method 2).

Controlled water samples containing DNA from common dolphin (Delphinus delphis) skin were prepared at 100, 1,000, and 10,000 ng/L in sterilized seawater. The DNeasy Blood & Tissue Kit (Qiagen) was used for DNA extraction. Ozone irradiation, 0.22 μm filtration, and autoclaving were used to provide the sterilized seawater. Three 1 L replicates of each sample were filtered and extracted according to the two methods described in Test 1. All spiked samples were eluted in 100 μL of DNase- and RNase-free water. For each method, filtration controls were obtained by filtering 1 L of sterilized seawater, and blank extraction controls were included. The total amount of extracellular DNA recovered was quantified using a Qubit™ dsDNA Broad Range (BR) Assay Kit.

Two DNA polymerase enzymes with different degrees of fidelity were evaluated on mock community (MC) and eDNA samples: a high-fidelity (>300 × vs. Taq enzymes) proofreading polymerase (Platinum SuperFi II DNA Polymerase, Invitrogen) and a non-proofreading polymerase (Supreme NZYTaq 2 × Green Master Mix, NZYTech).

To prepare MC samples, DNA was first extracted from tissue samples using the DNeasy Blood & Tissue Kit (Qiagen), following the manufacturer’s instructions with minor adjustments. DNA quantification was performed using a Qubit™ dsDNA Broad Range (BR) Assay Kit. Then, the extracted genetic material of the 16 vertebrate species was pooled together in equal concentrations to form an equimolar MC. The MC included DNA extracted from skin samples of nine cetacean species, tissue samples (muscle or fin) of six fish, and human blood (Suppl. material 1: table SS1). PCR amplification (following the first amplification protocol for environmental samples detailed in the Library preparation section), product cleanup using ExoSAP-IT, and subsequent Sanger sequencing were performed on the extracted DNA samples to confirm the identification of the species and obtain the amplicon sequences of the individuals pooled in the MC (amplicon sequences are provided in Suppl. material 2). Diluted MC (DMC) samples (1/10 in a mixture of different eDNA samples) were also used to evaluate the effectiveness of both polymerases under possible PCR inhibitors commonly present in environmental samples (Schrader et al. 2012). In addition to MC and DMC samples, eDNA samples extracted by Method 1 from Test 1 were also used for this experiment to further verify polymerase performance in regular environmental samples.

The performance of both polymerases in MC and DMC samples was compared statistically using paired t-tests. Each MC and DMC sample was amplified in triplicate with both polymerases from the same DNA preparation, yielding six paired comparisons (three MC and three DMC), and the paired test was therefore applied to account for the non-independence between samples sharing a common template. Reads assigned to fish and cetaceans were tested separately, and exact two-sided p-values were calculated. Normality of the within-pair differences was assessed with the Shapiro-Wilk test and was not rejected for either taxonomic group (fish: W = 0.956, p = 0.792; cetaceans: W = 0.937, p = 0.634). Reads assigned to fish and to cetaceans were tested separately, using two-sided tests with a significance level of α = 0.05, and effect sizes were expressed as Cohen’s dz. Statistical testing was not applied to comparisons involving environmental samples. At each of the three sampling stations, a single sample was processed per method, so no technical replicates were available within each station and method combination, and the resulting three observations per method were insufficient to assess distributional assumptions or support reliable inference. The resulting three paired observations cannot yield statistical significance regardless of effect size.

All environmental samples from Tests 1 and 3, plus the respective extraction and filtration controls, were used for library preparation and sequencing. A short region (approximately 245 base pairs) of the vertebrate mitochondrial 16S rRNA gene was amplified and reamplified using MarVer3(A) primers (Valsecchi et al. 2020). These primers were modified from the original set (swapping the purine base for an adenine at position 16 in the 5'–3' orientation of the forward sequence) to increase specificity to target taxa, such as marine mammals and fish (Afonso et al. 2024). PCR protocols were employed as described in Table 1.

A negative control sample containing no DNA (BPCRP) was included to check for contamination during library preparation. Library size was verified by running the libraries on 3% agarose gels stained with GreenSafe (NZYTech) and imaging them under UV light. After amplification, the libraries were purified using Mag-Bind RXNPure Plus magnetic beads (Omega Bio-tek), following the manufacturer’s instructions. The finished libraries were pooled in equimolar amounts according to the results of Qubit dsDNA HS Assay (Thermo Fisher Scientific) quantification. Finally, the pool was sequenced in a fraction of a NovaSeq PE250 flow cell (Illumina), aiming for a total output of 1 gigabase.

For Test 1, all eDNA samples were processed using the described methods. For Test 3, MC and DMC samples were used. In addition, because sequencing samples in different runs could introduce unwanted non-biological variation due to technical biases, all environmental DNA samples collected and extracted using Method 1 were resequenced simultaneously in the same sequencing run using both the Platinum SuperFi II and Supreme NZYTaq protocols.

All sample codes used in this study can be found in Suppl. material 1: table SS2. All raw data were made available in the European Nucleotide Archive under project PRJEB92411.

The bioinformatics analysis followed the Marine Vertebrate eDNA Metabarcoding bioinformatics pipeline (hereafter MVeM, available at https://github.com/pascoalf/MVeM). First, any potential traces of adapter dimers were removed using Cutadapt v3.5 (Martin 2011), and quality was assessed by verifying FASTQ files using FastQC software (Andrews 2010). Reads were then trimmed using Cutadapt v3.5 to remove primer sequences.

To infer amplicon sequence variants (ASVs), raw reads were processed in R (v4.4.1) using the DADA2 R package (v1.22; Callahan et al. 2016) with default parameters. Other packages required for raw data processing included ShortRead (v1.64.0; Morgan et al. 2009), stringr (v1.5.2; Wickham 2010), seqinr (v4.2.36; Charif et al. 2017), and dplyr (v1.1.4; Wickham et al. 2014). Based on quality profiles, reads were trimmed to optimize overall read quality. Specifically, forward and reverse reads were truncated at positions 180 and 150 for mock communities and 170 and 150 for eDNA samples, respectively. After ASV inference, the chimera filter was applied to reduce the chance of capturing artifact sequences. The raw data filtering process was tracked (Suppl. material 1: table S3). The output of the DADA2 pipeline was an ASV abundance table for all sequenced samples (Suppl. material 1: table S4). Rarefaction curves were generated by progressively subsampling sequence reads from each sample and plotting the number of observed ASVs against sequencing depth. Curves reaching an asymptotic plateau were considered indicative of sufficient sequencing depth to capture the diversity present in the samples and support subsequent diversity comparisons.

Taxonomic assignment of each ASV was performed by querying its representative sequence against the NCBI Nucleotide database (last accessed on 21 October 2025) using the BLASTn algorithm (v2.13.0+; Camacho et al. 2009). The resulting hits were filtered based on the following criteria: a minimum percent identity of 99%, an E-value threshold of 1 × 10^–5, and a minimum query coverage percentage of 80%. To retrieve full taxonomic information for the hits, the NCBI Taxonomy Toolkit—TaxonKit (Shen and Ren 2021)—was used. The BLASTn results were then processed using the following R packages: tidyr (v1.3.1; Wickham 2017), purrr (v1.1.0; Wickham et al. 2020), and ulrb (v0.1.8; Pascoal et al. 2025b). Sequences were filtered with a minimum alignment length of 190 bp to avoid misassignments to shorter non-targeted sequences. Additionally, matches labeled as ‘environmental samples’ were excluded to restrict the reference set to sequences derived from individual, tissue-based specimens, ensuring that taxonomic assignments relied on vouchered material rather than bulk eDNA sources. A geographic ban list containing species known to be outside the study’s geographic or ecological scope was also used to automatically exclude potential misassignments (e.g., Tursiops australis, Sousa chinensis, Stenella attenuata), with the aim of improving the resolution of target species that may occur in the study region. To further minimize erroneous matches to nuclear mitochondrial DNA segments (NUMTS), only ASVs with accession numbers corresponding to Chordata mitochondrial and 16S rRNA sequences (based on a curated list from the NCBI Nucleotide database, last updated on 16 October 2025) were retained for downstream analysis.

To refine taxonomic assignments, only the best hit per ASV was retained based on the highest bit score, a metric representing sequence similarity (Madden 2013). In cases where multiple hits shared identical top scores, further tiebreakers were applied based on percent identity and alignment length. If ties persisted among different species with equivalent bit score, percent identity, and alignment length, the Lowest Common Ancestor (LCA) approach was applied, assigning the lowest common taxonomic level (e.g., genus or family) to the respective ASV. The taxonomic assignments were filtered by target taxa (i.e., the classes Chondrichthyes, Actinopteri, and Mammalia) and merged with the ASV abundance table. The resulting ASV abundance-taxonomy table (Suppl. material 1: table S5), containing read counts and assigned taxonomy, was then imported for additional filtering steps. To address potential mis-tagging as a consequence of potentially incorrect index assignment during library preparation, sequencing, or demultiplexing steps (Illumina 2017), ASVs occurring at a frequency below 0.01% in each sample or with a single read were removed. For eDNA samples, an additional contamination filtering step was added to check for any ASV found in the filtration, extraction, and PCR amplification blanks processed alongside each sample batch. Some ASVs found in control samples may still represent true biological signals in environmental samples if they occur at high read abundance. To avoid bias from manually setting read thresholds, an automatic, unsupervised machine-learning method (via the ulrb R package) was used to classify ASVs within each sample based on abundance in relation to the number of reads in the control sample, ensuring consistent and sample-specific contaminant detection despite varying sequencing depths. From the filtered ASV abundance-taxonomy table (Suppl. material 1: table S6), only marine vertebrate assignments were retained for statistical analysis and data visualization using the ggplot2 (v3.4.3; Wickham 2011) and vegan (v2.7.2; Oksanen et al. 2007) R packages. Human reads were included in the mock community sample analysis, as humans were one of the species included in that mix. Downstream analyses included computation of α-diversity metrics (observed ASVs and the Shannon index at the species level). The bioinformatics pipeline was first validated using mock community samples.

The concentration of eDNA and diversity of marine vertebrate species obtained by both eDNA extraction methods (Methods 1 and 2) were compared using environmental samples collected at three different sampling stations. The eDNA concentration per sample obtained using each method is shown in Table 2.

The DNeasy PowerWater Sterivex standard protocol (Method 1) yielded higher eDNA concentrations per sample than the modified protocol (Method 2), which included an initial filter bead-beating step outside the Sterivex cartridge unit. The observed differences in eDNA concentrations corresponded to approximately six times higher yields for samples collected at Station 1, 17 times higher at Station 2, and 12 times higher at Station 3. Both methods yielded the lowest eDNA concentrations for samples collected at Station 1. eDNA concentrations in the filtration and extraction controls were near the detection limit of the Qubit HS Kit, being ≤ 0.01 ng/µL. The sequencing output resulted in a total of 218,833 raw reads obtained from the environmental samples using Method 1, whereas 323,218 raw reads were obtained using Method 2. After processing, 145,813 (Method 1) and 205,143 (Method 2) reads were retained for subsequent analyses. After applying a taxonomic filter to remove non-target ASVs, 38,470 reads (Method 1) and 90,750 reads (Method 2) were obtained. Across all stations, Method 1 recovered slightly more ASVs and species than Method 2 (Table 2). At Stations 1 and 2, both methods detected similar species diversity, whereas at Station 3, Method 1 detected a higher number of species than Method 2 (Table 2, Fig. 2A). The complete list of species detected at each station and by each method is also provided in Suppl. material 1: table S7. Moreover, variability in species composition was observed among replicate samples collected from the same sampling station and processed using Methods 1 and 2. Regarding α-diversity (Fig. 2B, C), the average richness (number of observed ASVs) was higher with Method 1, whereas the Shannon index showed slightly higher median values with Method 1, with wider dispersion for Method 2, largely driven by sample M2.1 (collected at Station 1).

For Test 2, results showed a positive trend between the amount of spiked DNA in the seawater samples and the concentration of DNA recovered using both extraction methods (Suppl. material 1: table S8). However, recovery efficiencies were low, with lower average values (0.28% with Method 1 and 0.09% with Method 2) in samples with a high f-exDNA content (100 ng) compared with those with medium spiked DNA amounts (10 ng/µL; 1.28% with Method 1 and 0.33% with Method 2). For Methods 1 and 2, samples with a low spiked DNA quantity (1 ng) exhibited an average DNA concentration (0.024 ng/µL and 0.007 ng/µL, respectively) similar to that observed in the filtration controls (0.017 ng/µL and 0.010 ng/µL, respectively). This observation suggests that minor environmental contamination may have occurred during the filtration process or that DNA quantification variability may have affected the results, as these values were very close to the detection limit.

A total of 516,283 reads were obtained for MC and DMC samples using Platinum SuperFi II, and quality filtering for analysis yielded a total of 348,984 reads. The sequencing output for the same samples using Supreme NZYTaq was 523,766 reads, of which 354,416 were retained after read quality filtering. Both enzymes successfully detected the same species present in the MCs, with dilution in other environmental samples not interfering with the results (Fig. 3A). Nevertheless, although the DNA from the different species in the MC was at equimolar concentrations, some species showed lower ASV read relative abundance than others (e.g., Phocoena phocoena, Engraulis encrasicolus, and Merluccius merluccius). Across all samples, amplification with Platinum SuperFi II generated 331,610 ASV reads assigned to the target taxa, comprising 234,848 cetacean reads and 96,762 fish reads. In comparison, amplification with Supreme NZYTaq produced 337,083 ASV reads assigned to the target taxa, including 302,201 cetacean reads and 34,882 fish reads. Thus, while Platinum SuperFi II yielded a higher proportion of fish-derived reads, Supreme NZYTaq consistently produced more cetacean-derived reads (Fig. 3B). Paired t-tests on the six paired comparisons confirmed that Platinum SuperFi II generated significantly more fish-derived reads than Supreme NZYTaq (mean difference = 10,313 reads, 95% CI [8,372, 12,255]; t(5) = 13.65, p < 0.001; dz = 5.57). Conversely, Supreme NZYTaq generated significantly more cetacean-derived reads (mean difference = 11,226 reads, 95% CI [2,835, 19,617]; t(5) = –3.44, p = 0.018; dz = 1.40), although the effect was considerably weaker and more variable across pairs than that observed for fish.

The same pattern was not observed when eDNA samples from Method 1 were rerun using Platinum SuperFi II or Supreme NZYTaq. Here, the species detected with each polymerase varied (Fig. 4). In this test, Platinum SuperFi II generated 451,631 raw reads, of which 192,164 remained after filtering. Under the same conditions, Supreme NZYTaq produced 337,607 raw reads, with 44,760 retained after filtering. Using Platinum SuperFi II to sequence the eDNA samples obtained by Method 1, 13 species were detected, whereas 11 taxa were detected using Supreme NZYTaq. The listed species are provided in Suppl. material 1: table S7. Both enzymes yielded distinct species compositions, with seven species detected in both cases, six detected only using Platinum SuperFi II, and four detected only using Supreme NZYTaq (Fig. 4A). By examining α-diversity metrics (Fig. 4B), samples amplified with Platinum SuperFi II exhibited a higher median number of ASVs compared with Supreme NZYTaq. Shannon diversity showed slightly higher median values for Platinum SuperFi II, whereas Supreme NZYTaq exhibited greater variability (Fig. 4B).

The rerun of environmental samples with Platinum SuperFi II also revealed differences between Tests 1 and 3. For the same environmental samples extracted using Method 1, both runs of Tests 1 and 3 resulted in the identification of 13 species in total, including seven fish species common to both runs (Sardina pilchardus, Boops boops, Engraulis encrasicolus, Zeugopterus punctatus, Diplodus sargus, Atherina presbyter, and Belone belone). The runs diverged in the remaining six species identified: in Test 1, five additional fish species (Oncorhynchus mykiss, Arnoglossus laterna, Parablennius gattorugine, Conger conger, and Spondyliosoma cantharus) and one cetacean (Balaenoptera physalus) were detected. In Test 3, five other fish species (Scomber scombrus, Merluccius merluccius, Ctenolabrus rupestris, Gaidropsarus mediterraneus, and Helicolenus sp.) and one different cetacean (Tursiops truncatus) were identified using the same amplification protocol.

In Test 1, two eDNA extraction strategies from Sterivex filters were compared: the standard DNeasy PowerWater Sterivex protocol (Method 1) and a modified PowerWater protocol based on a preliminary bead-beating homogenization step performed on the filters (Method 2). This second approach has previously been shown to enhance metabarcoding results in sediment-derived benthic eDNA (Hestetun et al. 2021). Cruaud et al. (2017) suggested that the limited available space inside the cartridge in Sterivex filters may reduce lysis efficiency during the first steps of the standard extraction protocol, which are carried out within the Sterivex unit.

The results showed that both extraction methods were effective for eDNA analysis, yielding similar overall diversity of marine vertebrates. However, the species identified by both methods varied at Station 3, probably owing to the intrinsic heterogeneity of the sample rather than differences in the extraction method. Although Method 1 consistently yielded higher eDNA concentrations across all three stations, the lack of technical replication prevents formal statistical validation of these differences, and they should therefore be interpreted as exploratory trends. The higher DNA yield with Method 1 suggests that removing the filter from the cartridge may not be necessary for effective lysis of eukaryotic cells in seawater samples during homogenization. This preliminary step added in Method 2 involves direct handling of the filter, which, in addition to increasing the risk of contamination, could even lead to a loss of genetic material during this process instead of ultimately increasing it. Nevertheless, the possibility that the modified protocol enhances DNA recovery from other target groups, such as prokaryotic organisms, or from other types of environmental samples, such as sediments, cannot be excluded, although this was beyond the scope of the present study.

Test 2 was based on the assumption that a portion of the f-exDNA could still be retained by the PES filter membranes because of specific chemical properties of the material (Liang and Keeley 2013). The results suggest a limited contribution of the f-exDNA fraction to the eDNA extracted using the two methods tested. Moreover, analysis of the recovery efficiency revealed a negative trend with increasing DNA concentrations in the spiked samples, which is likely linked to the saturation of available binding sites on the filter membrane. As in Test 1, the limited number of independent replicates per treatment precluded formal statistical testing, and the trends described here are therefore descriptive. Recent studies have suggested that extracellular DNA represents a substantial component of the total eDNA pool in aquatic environments (Torti et al. 2015; Nagler et al. 2022), accounting for up to 90% in marine systems (Torti et al. 2015). In addition to the f-exDNA subfraction, exDNA can also occur as weakly bound exDNA (wb-exDNA) and tightly bound exDNA (tb-exDNA), which mainly differ in the strength of their association with particles of the extracellular matrix, cell membrane proteins, and other organic or mineral matter. Previous eDNA studies on fish have shown that particle-associated exDNA in the size range of 1–10 µm is the most abundant fraction in aquatic environments (Turner et al. 2014; Wilcox et al. 2015). Although f-exDNA was not highly recovered in the eDNA extracted under the experimental conditions of this study, it is reasonable to assume that exDNA bound to organic and inorganic particles (wb-exDNA and tb-exDNA) in natural waters may occur at larger effective particle sizes than f-exDNA and therefore be more efficiently retained by filters with the pore sizes used here. Further research is needed to determine whether particle-bound exDNA subfractions represent substantial contributions to the extracted eDNA.

Analysis of the different polymerases (Test 3), which varied in their proofreading capacity, made it possible to assess their effectiveness in different taxonomic groups. The results for the MC samples showed that Platinum SuperFi II consistently amplified fish DNA more successfully, whereas for cetaceans, Supreme NZYTaq recovered significantly more reads, although with a considerably weaker and more variable effect. Oliver et al. (2015) reported that the use of a non-proofreading enzyme can overestimate species richness in environmental fungal community datasets, yielding a greater number of operational taxonomic units. However, this trend was not evident in this study, either in the MC or in the environmental samples, where the proofreading enzyme produced a greater number of ASVs at the stations tested. This suggests that, in the overall context of marine vertebrate eDNA, particularly from low-biomass or degraded samples, the enhanced fidelity and sensitivity of the proofreading polymerase may improve detection efficiency, even if it reduces amplification bias that favors certain taxa. Nevertheless, Supreme NZYTaq performed well for cetacean detection, suggesting that, in certain taxonomic groups, even non-proofreading enzymes may offer advantages in target recovery. The use of non-proofreading enzymes in library preparation offers a more cost-effective method for metabarcoding analyses of marine vertebrate communities, providing comparable species detection profiles. The accessibility of tools such as eDNA analysis is a current concern, and there is a need to reduce workflow costs to allow these methodologies to be applied more broadly (Hirsch et al. 2024). It is also important to consider that polymerase selection can introduce technical biases in metabarcoding results, particularly through sequence composition preferences, such as guanine-cytosine (GC) content. While PCR remains essential for metabarcoding, especially when the target DNA comprises only a small fraction of the total extracted material, differences in polymerase performance can affect the accuracy of relative abundance estimates. For example, taxa with a larger biomass may already be overrepresented in environmental samples, and amplification biases can further distort their apparent abundance. Recent studies have shown that many commonly used polymerases introduce GC-related biases that alter their relative proportions during PCR cycling. However, there is no evidence that occurrence data are similarly affected (Nichols et al. 2018). Although GC bias was not explicitly tested in this study, the observed variation in taxon detection patterns may partially reflect such effects.

Statistical comparison of polymerase performance was possible for the MC and DMC samples, where each sample was amplified with both enzymes from the same DNA preparation. For the environmental samples, the absence of replicates within each station and method combination did not allow such testing. Observed variability among eDNA samples in Test 3 likely reflects inherent environmental heterogeneity rather than methodological bias, although the limited sampling effort precludes definitive conclusions. However, the rerun of Method 1 samples with Platinum SuperFi II further illustrates this point, as distinct marine vertebrate diversity profiles were obtained compared with Test 1 despite the use of the same amplification protocol. In this case, such discrepancies are consistent with the methodological stochasticity inherent in metabarcoding workflows, encompassing random template dropout, variable amplification efficiencies, and sequencing-level variation in cluster formation and depth distribution (Leray and Knowlton 2016; Shirazi et al. 2021). These sources of variability can compound across the workflow, and no single strategy can fully resolve them (Gold et al. 2023). Approaches such as sequencing PCR replicates independently or increasing sequencing depth may help characterize this variability, but their effectiveness in marine vertebrate eDNA contexts remains to be evaluated and was beyond the scope of the present study. Overall, while proofreading polymerases offer some advantages, their performance gains appear modest and taxon-dependent, suggesting that cost-effective alternatives can be suitable for eDNA metabarcoding workflows, provided that replicate amplifications are included and results across sequencing runs are interpreted with appropriate caution.

MVeM automates the post-BLAST filtering steps, prioritizing reproducibility, adaptability to diverse case studies, and biological realism in taxonomic assignment. Sequence matches were sorted by descending bit score—a widely accepted metric for alignment quality that accounts for the number and quality of matched residues while penalizing gaps or mismatches. This method has been adopted in various metabarcoding studies targeting various taxa (de Vere et al. 2017; Kim et al. 2025). However, BLAST-based assignments can be confounded by closely related species with minimal sequence divergence, especially within families such as Delphinidae, where markers such as 12S and 16S rRNA often lack sufficient resolution to distinguish species, resulting in undetermined or incorrectly assigned species within this taxonomic group (Valsecchi et al. 2020; Agersnap et al. 2022). To mitigate potential misassignments in such cases, the pipeline includes a tie-resolution step that assigns the Lowest Common Ancestor (LCA) when the top hits share identical bit scores, percent identities, and alignment lengths. This conservative strategy prevents false-positive detections and overinterpretation of ambiguous results and enables classification at more reliable taxonomic levels. To further reduce uncertain matches, a customizable ban list was implemented to filter known non-target or geographically non-native species. This flexible exclusion tool is especially helpful in marine eDNA studies, where it is often difficult to distinguish closely related species inhabiting different areas. However, this step should be used with caution and in alignment with the specific objectives of the study, particularly when the goal is to detect shifts in species distribution patterns or the presence of non-indigenous taxa in aquatic ecosystems. In such cases, overly stringent filtering may obscure ecologically relevant signals, particularly when early detection is critical. Given these considerations, users should compile exclusion lists that reflect the ecological context and management goals of their work. A NUMTs filtering step was also included, as some mitochondrial sequences can be assigned to nuclear gene regions. NUMTs are fragments of mitochondrial DNA that have been transferred to the nuclear genome and can be co-amplified during PCR, potentially leading to false positives in metabarcoding studies (Song et al. 2008; Andújar et al. 2021). This step minimizes the risk of overestimating biodiversity due to nuclear pseudogenes and improves the accuracy of community composition assessments.

Given the sensitivity of eDNA workflows to contamination (Sepulveda et al. 2020), the pipeline also incorporates a contamination control filter to remove ASVs detected in filtration and extraction blanks while retaining ASVs that were substantially more abundant in samples than in the corresponding blanks. Selecting a relative abundance threshold to define abundant or rare taxa is subject to bias (Pascoal et al. 2025a). Therefore, an option was included to automatically identify which ASVs were abundant within each environmental sample using unsupervised machine learning (Pascoal et al. 2025b). It is important to note, however, that this automatic method does not allow users to manually adjust which ASVs are classified as potential contaminants. Importantly, this framework also helps preserve ecologically meaningful detections that might otherwise be discarded, as in cases in which a true biological signal is overwhelmingly abundant in environmental samples but appears at very low levels in controls. Such patterns highlight how data-driven classification can reduce the risk of losing key ecological information while maintaining a conservative approach to contamination filtering, thereby supporting traceable mitigation strategies that enhance both the precision and interpretability of eDNA detection during extensive or long-term sampling efforts.