Investigating passive eDNA samplers and submergence times for marine surveillance

Field deployment

A total of 72 samples (excluding negative controls) were deployed from pontoon B65 (−35.836802, 174.468679) in Marsden Cove, Whangārei, Aotearoa New Zealand from the 19^th–20^th of June 2021 for a 12 h period. Cawthron holds a Special Permit with the Ministry for Primary Industries (SP822-2) that allows the taking of fish, aquatic life and seaweed for the purposes of education and investigative research. Four substrates were explored; nylon filters (5 µm pore size; Sterlitech, Auburn, WA, USA), positively charged nylon discs (0.45 µm pore size; Biodyne™; Thermo Fisher Scientific, Waltham, MA, USA), UV-treated nylon mesh cut in 50 × 50 mm pieces (20 µm pore size; NITEX™; Sefar Ltd., Thal, Switzerland), sterile Whirl-Pak^® Speci-Sponges^® (Nasco, Merck, Germany), and additionally the sponge water was squeezed from the Whirl-Pak^® Speci-Sponges^® and DNA extracted directly from the supernatant as in von Ammon et al. (2018a) (Fig. 1). Substrates were attached with cable ties to a bottle carrier and deployed at ~1 m below the water surface. Three replicates of each substrate were retrieved at 10, 30 min, 1, 3, 6 and 12 h after deployment. Upon retrieval, the membranes were cut in half and placed in 2 ml Eppendorf tubes prefilled with LifeGuard Soil Preservation Solution (Qiagen, Hilden, Germany), except Whirl-Pak^® Speci-Sponges^® which were placed into individual bags due to their large size and put on ice.

The water isolated from the Whirl-Pak® Speci-Sponges® was later analyzed as an additional individual sample type.

For reference (waterborne eDNA detection), triplicate 500 ml seawater samples were collected directly in front of the bottle carrier at the start (0 h) and end (12 h) of passive eDNA sampler deployment. These samples were vacuum filtered through sterile Whatman filters (pore size 0.45 µm; Merck KGaA, Darmstadt, Germany) placed in a Sterifil filtration system (Merck, Darmstadt, Germany) using a 12V Seaflo 21 Series Water Pressure Pump (SEAFLO, Fujian, China). Each filter was stored in a 2 ml Eppendorf tube containing 1.5 ml of LifeGuard Soil Preservation Solution (Qiagen, Hilden, Germany).

Negative field controls (blanks) consisted of (i) 500 ml ultrapure water filtered through the same 0.45 µm filter type as water field controls, and (ii) all materials used for passive filtration were taken out briefly in the field but then immediately placed into clean tubes or bags.

Laboratory processing

All samples were extracted using a modified protocol of the Qiagen DNeasy Blood & Tissue Kit (Jeunen et al., 2019). Materials were placed into new 2 ml Eppendorf tubes with 0.5 mm Zirconia/Silica Beads (Dnature, Aotearoa New Zealand) and 720 µL ATL buffer (Qiagen, Hilden, Germany), and 80 µL Proteinase K (Qiagen, Hilden, Germany) added. The original 2 ml Eppendorf tubes with 1.5 ml of LifeGuard Soil Preservation Solution (Qiagen, Hilden, Germany) were centrifuged at 10,000 g for 5 min, supernatant discarded, and the remaining pellet re-eluted with ATL buffer and transferred to the bead beating tube. Samples were homogenized at maximum speed (1,500 strokes/min) via bead beating for 2 min (1600 MiniG Spex SamplePrep) and incubated overnight at 56 °C. The following day samples were vortexed for 1 min, followed by a centrifugation step for 1 min at 6,000 g. 600 µL of supernatant was transferred to a new 2 ml Eppendorf tube, after which the standard Qiagen DNeasy Blood & Tissue Kit protocol was followed.

Species-specific assays

Droplet digital polymerase chain reactions (ddPCRs) were run on a QX200 Droplet Digital PCR System™ (Bio-Rad, Hercules, CA, USA) for four targeted assays amplifying the Mediterranean fanworm Sabella spallanzanii, the clubbed tunicate Styela clava, the Japanese kelp Undaria pinnatifida and the brown bryozoan Bugula neritina. Sabella spallanzanii specific Cytochrome c oxidase subunit 1 (COI) primers forward Sab3F: 5-GCT CTT ATT AGG CTC TGT GTT TG-3 and reverse Sab3R: 5′-CCT CTA TGT CCA ACT CCT CTT G-3′ and Taqman probe Sab3 5′-FAM/AAA TAGT TCA TCC CGT CCC TGC CC/BkFQ-3′ were used as in Wood et al. (2017). For Styela clava, we used the COI primers forward SC1F 5′-TCCGGCGGTAGTCCTTTTATT-3′ and reverse SC1R 5′-GAGATCCCCGCCAAATGTAA-3′ as well as and the TaqMan® probe SC1 5′-TTA GCTAGGAACTTGGCCCA-3′, as in Gillum (2014). For both assays, each reaction included 450 nM (1 µl) of each primer and probe, 1 × ddPCR Supermix for probes (10 µl, no UTP; Bio-Rad, Hercules, CA, USA) and 3 µl DNA and 5 µl sterile water for a total reaction volume of 21 µl. Amplification followed the cycling protocol: 95 °C for 10 min, 40 cycles of 94 °C for 30 s, 60 °C for 1 min, and 98 °C for 10 min. For Undaria pinnatifida, we used the COI forward primer U_pinnatifida2F 5′-TACAGCAATGTCTGTTTTTATCC-3′, reverse primer U_pinnatifida2R 5′-ACATTATACAACTGATGATTTCCC-3′ and probe U_pinnatifida2Probe_BHQ2 5′-FAM-ATTGCAATTAGCTAGCCCTG/3BHQ_2-3′ (Bott, Giblot-Ducray & Deveney, 2015). Each reaction included 450 nM of each primer and probe (1 µl), 1× ddPCR Supermix for probes (10 µl, no UTP; Bio-Rad, Hercules, CA, USA) and 3 µl of DNA and 5 µl sterile water for a total reaction volume of 21 µl. The amplification cycling protocol was as follow: 95 °C for 10 min, 40 cycles of 95 °C for 30 s, 57 °C for 30 s, and 98 °C for 10 min. For Bugula neritina, an Evagreen assay was used with the COI primers BuNe_SF 5′-GGTACATTATACTTTTTATTTGGAC-3′ and BuNe_SR 5′-CCCCCA ATTATAACTGGTATG-3′ primers as in Kim et al. (2018), and conditions adapted as in Scriver et al. (2024a) to reactions that included 450 nM of each primer, 1× ddPCR Evagreen Supermix (no UTP; Bio-Rad, Hercules, CA, USA) and 3 µl of DNA and 6 µl sterile water for a total reaction volume of 21 µl. The amplification protocol was as follow: 95 °C for 5 min, 40 cycles of 95 °C for 30 s, 57 °C for 1 min (2′ ramp), and 90 °C for 5 min. Each well of the plate was then individually analyzed on the QX200 instrument to establish the threshold value separating negative and positive droplets and perform absolute quantification of target DNA (see Chiu et al., 2017). A negative control (DNA-free water) and positive control (consisting of extracted DNA from tissue of the targeted species) were included on each plate.

18S rRNA metabarcoding

For metabarcoding analyses, polymerase chain reactions (PCRs) were performed to amplify the V4 region of the nuclear small ribosomal subunit (18S rRNA) gene. The primers for the 18S rRNA gene were Uni18SF: 5′-AGG GCA AKY CTG GTG CCA GC-3′ and Uni18SR: 5′-GRC GGT ATC TRA TCG YCT T-3′ (Zhan et al., 2013) with Illumina overhang adaptors (Kozich et al., 2013). Amplifications were undertaken on an Eppendorf Mastercycler (Eppendorf, Hamburg, Germany) in a total volume of 50 μl using 25 μl of MyFi Mix (Bioline, Meridian Bioscience, Cincinnati, OH, USA), 1 μl of each primer, 20 μl of DNA-free water, and 3 μl of template DNA. Thermocycling conditions were as follows: 95 °C for 2 min, followed by 40 cycles of 95 °C for 20 s, 52 °C for 20 s, 72 °C for 20 s, and a final extension of 72 °C for 10 min. Negative PCR controls (3 µL of DNA-free water) were run alongside the samples in the PCR. Amplicons were cleaned and normalized using the SequalPrep Normalization kit (Thermo Fisher Scientific, Waltham, MA, USA) resulting in a concentration of ~1 ng/μl. Amplicons were sent to Sequench Limited (Nelson, Aotearoa New Zealand) for indexing with Nextera XT (Illumina, San Diego, CA, USA) kit and high-throughput sequencing on an Illumina MiSeq platform. The library pool was diluted to a final loading concentration of 6 pM with a 15% PhiX spike and sequenced using MiSeq Reagent Kit v3-600 cycle (2 × 301 bp) (Illumina, San Diego, CA, USA). A sequencing control (DNA-free water) was added to each sequencing run.

Bioinformatics

FASTQ files were demultiplexed and primers removed using Cutadapt (version 3.7; Martin, 2011), using a minimum overlap of 15 bp. Forward and reverse sequences were truncated on their 3′ end at 225 and 216 bp for forward and reverse reads, respectively, to remove low quality sections of the reads. Sequences were quality filtered and denoised using the default parameters of the ‘DADA2’ R package (version 1.21; Callahan et al., 2016). Forward and reverse reads were merged using a minimum overlap of 10 bp, and chimeric sequences removed using the consensus method of ‘DADA2’. Remaining sequences were taxonomically assigned with the SILVA database (version 132) using a Naive-Bayesian classifier and with blastn and megablast against the GenBank nt database using default settings of blastn_taxo_assignment from the biohelper R package (https://github.com/olar785/biohelper, following Laroche et al. (2020). Taxonomic assignments for the three approaches were merged using the taxo_merge function of biohelper, where taxonomy is first normalized using NCBI’s curated taxonomy database, using the highest taxonomic assignment (resolution) across methods if there is consensus among them (>50%), across assigned taxonomic ranks. In case of discrepancy, the last common ancestor is assigned. Sequences found in blanks were investigated and potential contamination removed using the MicroDecon R package (version 1.0.2, McKnight et al., 2019). For removal of potential contamination, microDecon compares the taxa (ASVs) present in the negative controls to those in the study samples. Those present in both the negative controls and the study samples are flagged as potential contaminants. The microDecon algorithm assumes that the proportion of each putative contaminant in the study samples can be approximated based on its relative abundance in the negative controls. For each shared ASVs (present in both controls and samples), it calculates the expected contaminant proportion in the sample which is then subtracted. The advantage of this method over simply removing all ASVs found in blanks is that it enables keeping relevant and often abundant/important taxa that may have been detected in blanks. In this study, we have used the microDecon’s default parameters and have included ASVs identified in negative controls in the Supplemental Material (see Table S2).

Raw sequence reads were deposited in the NCBI short read archive under the accession number PRJNA1154297.

Statistical analysis

Bar plots were used to visualize species-specific ddPCR detections (copies/µl) using ggplot2 package (v3.5.0.; Wickham, 2016). The mean value of active filtration was overlayed as a horizontal line including standard deviation. The relationship between deployment time and detected DNA (copies/µl) via ddPCR for the four investigated non-indigenous species (NIS) was tested per species per matrix using a generalized linear model (GLM) with an exponential link function and visualized using the ggpubr R package (version 0.6.0; Kassambara, 2018).

For the 18S rRNA metabarcoding dataset, sequences identified as non-eukaryotic or non-marine taxa were discarded, and rarefaction curves performed on the remainder of the data to investigate sequencing depth per sample using the ‘rarefy_even_depth’ function (Fig. S1, phyloseq v.1.40–1 R package). Samples with read sum <1,000 reads were removed. The effect of deployment time and matrix on overall amplicon sequence variants (ASV) richness and those solely associated to NIS was investigated and visualized in scatter plots including R² and p-values with ggpubr. Eukaryotic community composition was visualized with bar plots using the biohelper R package (https://github.com/olar785/biohelper), following Laroche et al. (2020), for overall biodiversity at phylum and family level. For subsampled NIS taxa at genus and species level, the non-rarefied fasta reads for 18S rRNA were screened with the Pest Alert Tool (https://gitlab.com/cawthron-public/marine-biosecurity-toolbox/pest-alert-tool, Zaiko et al., 2023) using Minimum % sequence identity match = 99.5% and Minimum sequence length (300 bp). Beta-diversity was initially assessed with a principal component analysis (PCA) on the 18S rRNA dataset with deployment time displayed with an ordisurf layer. In addition, the correlation of deployment time with biological assemblages is showed with an arrow, with Pearson (r) and p-value displayed. The effect of deployment time and matrix was tested with a permutational analysis of variance using the vegan R package (version 2.6.9; Oksanen et al., 2018). The progression of beta-diversity dissimilarities for each matrix and along submersion time was plotted in a line diagram using phyloseq_group_dissimilarity function of the metagMisc R package (v 0.5.0).

Species-specific detections of key non-indigenous species

From the four species-specific ddPCR assays, only Undaria pinnatifida was not detected, neither with any of the passive samplers nor with the vacuum filtering, the later undertaken for reference purposes. However, Sabella spallanzanii, Styela clava and Bugula neritina were detected throughout submersion times from the different substrates but displayed variable detection signals and patterns (Figs. 2 and 3). Sabella spallanzanii showed the strongest amplification signals of the tested assays, ranging from a maximum value of 66 copies/µl to a median value of 0.1185 copies/µl. Interestingly, vacuum filtration for S. spallanzanii, which was only undertaken at time point 0 and 720 min (12 h), showed consistent copy numbers (mean = 0.9725 ± 0.29 copies/µl). All passive samplers but not sponge water showed a significant positive temporal trend of copy numbers/µl with time (p ≤ 0.05, Fig. 3), with nylon discs reaching the reference mean value of vacuum filtered samples after 60 min (Fig. 2). After 720 min, nylon disc exceeded the reference value with 2.33 copies/µl and nylon mesh with a maximum value of 66.1 copies/µl (median 0.3 ± 14.74 copies/µl).

For Styela clava, species-specific detection was more sporadic across sampling events with 0 detections for all samplers at 180 min (Fig. 2). This species was detected only in 50% of the water samples, with relatively low signal (mean = 0.027 ± 0.03 copies/µl). The passive sampler nylon mesh still displayed a significant increase in copy numbers over time (p ≤ 3e−04, Fig. 3). Nylon disc and nylon mesh exceeded the reference values from active filtration after 30 min and 1 h, respectively. At the final sampling time point (720 min), all passive sampling matrices showed detections (Fig. 2), with nylon disc (0.08 copies/µl), nylon mesh (0.19 copies/µl) and sponge water (0.02 copies/µl) reaching or exceeding the reference waterborne values.

The B. neritina assay yielded positive detections across the different sample types and time points (Fig. 2) but remaining at the lower end of the overall copy numbers (0.05–0.2 copies/µl) and showing no significant pattern of eDNA signal accumulation over time except a decreasing trend for nylon, highly influenced by the high initial copy numbers (p ≤ 0.0019, Fig. 3). The reference value from water samples (mean = 0.06 ± 0.08) was exceeded with nylon, nylon disc, nylon mesh and sponge within 10 to 60 min. All passive samplers except for sponge water showed detection signals after 720 min (Fig. 2).

Metabarcoding data

High-throughput sequencing yielded a total of 12,318,746 18S rRNA reads, of which 85% remained after primer trimming, 44.1% after denoising, 40% after merging, and 36.1% after chimera removal (Table S1). Removing potential contamination from negative controls (see Table S2) and unidentified amplicon sequence variants (ASVs) resulted in a total of 4,076,470 sequence reads (mean of 42,463 reads per sample) and a total of 4,046 ASVs. Rarefaction curves plateaued before 10,000 reads, therefore samples were rarefied to the minimum number of reads per sample (19,470 reads) prior to diversity analysis (see Fig. S1).

The overall ASV richness displayed significant decreasing temporal trends for all passive samplers, including sponge water (p ≤ 0.05, Fig. 4A). Similarly, while not significant, a negative trend was observed for the NIS ASVs except for the sponge water samples, which indicated some increase over time, also not significant (Fig. 4B).

Looking at the taxonomic composition across sampling times and sample types, a clear difference could be observed between filtered water samples and the passive sampling (Fig. 5A). The majority of taxa sampled by water filtering was assigned to the copepods Calanoida, Oithonidae and Styelidae (Chordates), with a clear shift to a different family of copepods, Tachidiidae, observed at 720 min compared to the beginning of the experiment.

For the passive sampling, a successive shift through time with a clear dominance in relative abundance of Styelidae and unidentified Eukaryotes was reported (Fig. 5A). In addition, the relative abundance of the polychaete family Syllidae substantially increased from 10 to 30 min, before decreasing at later sampling time points. After 720 min, Podocopida (Ostracoda) became the most abundant taxonomic group, particularly for all nylon samplers.

Looking at the taxonomic composition of exclusively non-indigenous species, filtration and passive sampling approaches appeared very similar (Fig. 5B). Earlier samples (10 min–6 h) were dominated by Ciona savignyi and Styela plicata. After 720 min, Ascidiella aspersa clearly dominated across all samples. In filtered seawater samples, relative abundance of NIS decreased between the first and last sampling point, while for nylon, nylon disc, nylon mesh and sponge water, relative NIS abundance increased over time. The relative abundance values varied more drastically and appeared less consistent among sampling time points for the sponge sampler compared to the other passive samplers (Fig. 5B).

A permutational analysis of variance (PERMANOVA) performed to assess the effect of submersion time and sampling method among passive samplers confirmed significant effects of both factors and their interaction (p ≤ 0.001, Table 1). When tested separately, the effect of submersion time remained significant for all sample types (p ≤ 0.001, Table S3, Fig. S2).

In Fig. 6, the progression of beta-diversity dissimilarity shows an increasing trend over time for all matrices (10–720 min), although for nylon it nearly reaches a plateau at 360 min. While all nylon matrices very consistently increase in beta diversity dissimilarities, the sponge matrix shows a drop of dissimilarity at 60 and 360 min, while the sponge water steeply increases in dissimilarity until 60 min, drops slightly at 180 min and then increases again until 720 min. The beta diversity dissimilarity for the nylon matrix increased gradually from 30 to 360 min, with a more significant increase observed between 360 and 720 min. Overall, dissimilarity in community composition at the end of the experiment for passive samples was substantially higher than that of the active filtration samples (mean = 0.63, sd = 0.04; Table S4).