A power analysis (Zar 1999) was used to determine the required number of experimental replicates per treatment to test the above null hypothesis using the following equation:

The power analysis calculates sample size (n) or the number of 100 ml experimental replicates, for each response variable; including filter material, pore size and extraction type. We performed the test assuming a 90% confidence that the average DNA concentration after amplification or count of OTUs deviates (d), with a confidence interval (CV) no more than 10% from the true mean (µ) post-amplification DNA concentration or count of OTUs. The true mean is assumed equivalent between replicates. We used the Critical Values calculated from the t distribution assuming a one-tailed test (α=0.05, (α=0.05, Zar 1999) and iteratively solved for sample size until the initial and predicted sample sizes were equivalent. Our statistical power analysis with the above assumptions indicated that we needed at least five replicates per experimental treatment. We therefore strove for ten replicates per experimental treatment (Figure 1) to increase our confidence and compensate for any event during the experiment that would lead to the loss or omission of any replicate.

The Port of Southampton in the United Kingdom is amongst the busiest and oldest commercial ports in Europe (McCutcheon 2008). The port is protected from wave exposure by the natural harbour created by the Solent Estuary and the Isle of Wight. A modelling approach to identify invasion risk categorised the Solent as one of the highest risk areas for marine species invasions (Pearce et al. 2012). Furthermore, surveys of Southampton marinas have found many of the most common NIS in the UK, making it an ideal location to test whether eDNA metabarcoding methods are sensitive enough to detect NIS in commercial ports (Bishop et al. 2015).

To allow for broad taxonomic richness in the water used for the experiment, thirteen 1 litre seawater samples were collected from Empress Dock at two sites: seven from the harbour master’s pontoon (50°53’24.7128”N, 1°23’47.2704”W) in the southwest corner and six from the National Oceanography Centre pontoon (50°53’28.2732”N, 1°23’38.0724”W) in the northeast corner. The 1 litre Nalgene™ plastic bottles were decontaminated with a 4% bleach solution (0.58% NaOCl) and rinsed using sterile Milli-Q water. Plastic bottles were then used to collect water from the top 0.5 m of the water column, first being rinsed out twice in the water and, finally, collecting the 1 litre water sample with the third fill. After collection, samples were immediately transported to the laboratory directly near the dock where they were refrigerated at 4 °C. Samples were refrigerated for 12.5 hours until filtration could begin.

Prior to filtration, the thirteen 1 litre seawater samples were poured into a single container. Experimental replicates were created by taking 100 ml volumes of seawater from the single container (Figure 1). The single container of seawater was continuously mixed to minimise sedimentation of particulates and to ensure homogeneity across experimental replicates. Each 100 ml experimental replicate was randomly assigned to a filter type (Figure 1). To compare pore size and filter material, two pore sizes (0.7 µm and 1.2 µm) of Glass Microfibre Filters (GF) and three pore sizes (0.45 µm, 0.65 µm and 1.2 µm) of Cellulose Nitrate Filters (CN) were tested: GF 0.7 µm (Whatman), Grade GF 1.2 µm (Whatman), CN 0.45 µm (GE Healthcare), CN 0.65 µm (GE Healthcare) and CN 1.2 µm (GE Healthcare). All filter types had 47 mm diameters. Filters where chosen because they represent commonly used material and pore sizes. Pore sizes between GF 0.7 µm and CN 0.65 µm were considered equivalent pore sizes for statistical analysis of the experiment. Filtration occurred on a three-piece stainless-steel manifold with the capacity to processes three experimental replicates in parallel.

After each filtration, filters was folded and placed into a 2 ml microcentrifuge tube filled with 700 µl of Longmire’s tissue lysis buffer (Longmire et al. 1997), ensuring that the filter was completely submerged and stored at -20 °C. The filters in Longmire’s tissue lysis buffer were shipped to the University of Notre Dame at ambient temperature for processing. Longmire’s buffer, as it is commonly called, has been shown to be effective for storage of filtered eDNA samples at room temperature (Renshaw et al. 2015). Use of a lysis buffer for preservation of eDNA samples is favourable because it has a pH suitable for stabilising DNA and additives capable of deactivating nucleases.

A rigorous cleaning protocol, before and during the procedure, ensured minimal cross-contamination between filters. All equipment was decontaminated before the procedure using a 4% bleach solution (0.58% NaOCl). The three-piece stainless-steel manifold and working area were first wiped down with the 4% bleach solution and then rinsed with distilled water and dried. New gloves and forceps were used for every experimental replicate and the decontamination procedure was repeated between replicates.

Experimental replicates were randomly assigned by filter type to an extraction method as indicated in Figure 1. Three DNA extraction protocols were tested. Two consisted of phenol-chloroform-isoamyl protocols (PCI-1, PCI-2) and the third protocol used a chloroform-isoamyl protocol (CI). The exact step-by-step protocols that were followed are provided in Suppl. material 1. Protocol PCI-1 is a modified protocol that was developed by combining aspects of both the PCI-2 and the CI protocols. The PCI-1 protocol was previously used by Lacoursière-Roussel et al. (2018) and the main difference from either of the other two protocols is that it utilised QIA shredder filters (Qiagen, Germany), which have been reported to assist with performance of extracting DNA from environmental samples (Goldberg et al. 2013). Protocol PCI-2, which was previously used in comparison to commercial extraction kits and was shown to perform well when targeting both eukaryotes and bacteria (Deiner et al. 2015), was used here unmodified except to account for the difference in filter size used by increasing the volume of tissue lysis buffer added. Protocol CI was also previously used in comparison to kit extractions and is an inexpensive alternative that performed as well as other extraction methods (Renshaw et al. 2015). The CI protocol used in this study was unmodified. The CI protocol omits the chemical phenol, which is a toxic substance requiring extra safety precautions and added processing time and thereby considered a safer and less time-consuming option. DNA extractions were eluted in 75 µl of water or buffer as indicated in Suppl. material 1.

To remove potential inhibitors, extracted DNA was treated with the OneStep™ PCR Inhibitor Removal Kit (Zymo Research, Irvine, California, USA). Each replicate was quantified with 2 µl of sample using Qubit dsDNA High-Sensitivity Kit and Qubit 2.0 Fluorometer (Life Technologies, Grand Island, NY).

Our 18S ribosomal RNA gene primer set contained two sequences: Illumina adapter sequence and gene-specific target sequence. The Illumina adapter sequence provided the requisite flanking region for subsequent dual-indexing PCR and sequencing platform compatibility. The gene-specific sequence targeted a 378 bp fragment of the V4 variable region of the Eukaryote 18S ribosomal RNA gene. Hadziavdic et al. (2014) demonstrated in silico that the 18S_574F and 18S_952R primer set was the most optimal primer pair to cover eukaryote taxa while efficiently excluding prokaryotes. Together, these primers allowed for the generation of a eukaryotic 18S ribosomal RNA gene library compatible with Illumina sequencing platforms.

The library preparation workflow is illustrated in Figure 2. The eukaryote 18S ribosomal RNA gene V4 variable region (378 bp) was amplified from extracted DNA by high fidelity PCR using iProof PCR reagents (Bio Rad Laboratories, Hercules, CA), and our 18S ribosomal RNA gene primer set in a 0.2 ml PCR tubes with individually attached cap. The 50 µl-reaction contained 5.0 µl of extracted DNA (regardless of concentration), 27.0 µl of PCR-grade water, 10.0 µl of 5× HiFi Buffer (Bio Rad Laboratories, Hercules, CA), 1.5 µl of 50 mM magnesium chloride (Bio Rad Laboratories, Hercules, CA), 1.0 µl of 10 mM dNTP mix, (G-Bioscience, St. Louis, MO), 2.5 µl of 10 µM 18S_574F primer (Integrated DNA Technologies, Coralville, IA), 2.5 µl of 10µM 18S_952R primer (Integrated DNA Technologies, Coralville, IA) and 0.5 µl of 2 U/µl iProof HiFi DNA Polymerase (Bio Rad Laboratories, Hercules, CA). PCR reactions were cycled in a Master Pro S (Eppendorf, Westbury, NY) with the following conditions: initial denaturation at 98.0 °C for 2 minutes, followed by 25 cycles of denaturation at 98.0 °C for 10 seconds, annealing at 55.0 °C for 20 seconds and primer extension at 72.0 °C for 30 seconds, with a final extension at 72 °C for 10 minutes. Reactions were purified with a 0.8× bead:sample volume ratio (40.0 µl) of Agencourt Ampure XP Beads (Beckman Coulter, Webster, TX) and separated with a DynaMag-2 magnet rack (Life Technologies, Grand Island, NY). Ampure XP Beads were resuspended in 52.5 µl PCR-grade water to elute PCR products. After clearing the solution of Ampure XP Beads with the DynaMag-2 magnet, only 50.0 µl was recovered to avoid Ampure XP Bead carry-over. PCR product concentration was determined with 2.0 µl of sample using Qubit dsDNA High-Sensitivity Kit and Qubit 2.0 Fluorometer (Life Technologies, Grand Island, NY).

Figure 2.

Library preparation workflow illustrating pooling by volume rather than by normalising to avoid removing the effect of the experimental treatments.

The Illumina adapter and dual-index barcodes (i7 and i5 index primers) were added to the PCR-1 products by high fidelity PCR using iProof PCR reagents (Bio Rad Laboratories, Hercules, CA) in a 0.2 ml PCR tube with individually attached cap. The 50 µl-reaction contained 5.0 µl of PCR-1 Template (regardless of concentration), 22.0 µl of PCR-grade water, 10.0 µl of 5× HiFi Buffer (Bio Rad Laboratories, Hercules, CA), 1.5 µl of 50 mM magnesium chloride (Bio Rad Laboratories, Hercules, CA), 1.0 µl of 10 mM dNTP mix, (G-Bioscience, St. Louis, MO), 5.0 µl of 10 µM i7 Index primer (Integrated DNA Technologies, Coralville, IA), 5.0 µl of 10 µM i5 Index primer (Integrated DNA Technologies, Coralville, IA) and 0.5 µl of 2 U/µl iProof HiFi DNA Polymerase (Bio Rad Laboratories, Hercules, CA). PCR reactions were cycled in a Master Pro S (Eppendorf, Westbury, NY) with the following conditions: initial denaturation at 98.0 °C for 2 minutes, followed by 8 cycles of denaturation at 98.0 °C for 10 seconds, annealing at 60.0 °C for 20 seconds and extension at 72.0 °C for 30 seconds, with a final extension at 72.0 °C for 10 minutes. Reactions were purified with a 0.8× bead:sample volume ratio (40.0 µl) of Agencourt Ampure XP Beads (Beckman Coulter, Webster, TX) and separated with a DynaMag-2 magnet (Life Technologies, Grand Island, NY). Ampure XP Beads were resuspended in 32.5 µl PCR-grade water to elute the dual-indexed library. After clearing the solution of Ampure XP Beads with the DynaMag-2 magnet, only 30.0 µl was recovered to avoid Ampure XP Bead carryover. Library concentration was determined with 2.0 µl of sample using Qubit dsDNA High-Sensitivity Kit and Qubit 2.0 Fluorometer (Life Technologies, Grand Island, NY).

No-template controls were introduced at multiple steps in the workflow (DNA extraction, eukaryotic 18S ribosomal RNA gene PCR-1 and PCR-2). No-template controls followed the same procedure as experimental replicates. In place of extracted DNA, 5.0 µl of PCR-grade water was used as the template for PCR-1 and PCR-2. A total of ten no-template controls were generated to detect contamination during seawater sampling and preparation of libraries.

Due to the project’s scale, 31 of 129 dual-indexed libraries (~24%), derived from experimental replicates, were randomly chosen for Bioanalyzer analysis to control overall cost of validation. To verify the final library fragment size and diagnose primer contamination, 1.0 µl of dual-indexed library was analysed using Bioanalyzer 2100 instrument and Bioanalyzer DNA 7500 kit (Agilent Technologies, Santa Clara, CA). All no-template controls (1.0 µl) were analysed using Bioanalyzer DNA High-Sensitivity kit (Agilent Technologies, Santa Clara, CA).

Due to the number of samples in the study, the libraries where split across two MiSeq runs. To avoid a run effect, libraries were randomly assigned to one of two pools by sub-group (explained below), extraction method and filter material. All ten no-template control libraries were added to both pools to test for consistency between runs. As PCR-2 DNA concentrations from each filter material were the response variables, the volume of template added to the eukaryotic 18S ribosomal RNA gene PCR and dual-index PCR was held constant without normalising the template concentration as this would remove the effect of our experiment (Fig. 2). Furthermore, we multiplexed equal volumes of library or no-template control by group and pool assignment (Fig. 2). Specifically, 5.0 µl of library (regardless of concentration) and 5.0 µl of no-template control were combined by group. This created six (6) sub-pools: PCI-1 Pool-1, PCI-1 Pool-2, PCI-2 Pool-1, PCI-2 Pool-2, CI Pool-1 and CI Pool-2. To remove primer contamination, each sub-pool was purified with a 0.8× bead:sample volume ratio of Agencourt Ampure XP Beads (Beckman Coulter, Webster, TX) and separated with a DynaMag-2 magnet (Life Technologies, Grand Island, NY). Ampure XP Bead volume was calculated by multiplying 0.8 and the total volume of the sub-pool. Ampure XP Beads were resuspended with 1.089 µl PCR-grade water per library and no-template control in the sub-pool. This volume was chosen based on the desire for an elution volume of about 32 µl to effectively concentrate the pooled libraries by 4 to 5 times. In practice, 32 µl would be divided by the sum of the number of libraries and no-template controls in the pool to determine the elution volume per library (Fig. 2). (Example: If a sub-pool with 25 libraries and 4 no-template controls had a total volume of 145 µl, then116 µl of Ampure XP Beads would be added. Ampure XP Beads would be resuspended with 31.58 µl of PCR-grade water and then 29.0 µl recovered after bead separation.). After clearing the solution of Ampure XP Beads with the DynaMag-2 magnet, 2.5 µl less than the elution volume was recovered to avoid Ampure XP Bead carry-over. Finally, 0.5 µl per library and no-template control of the concentrated sub-pools were combine by final pool assignment. (Example: Subpool Group M with 25 libraries and 4 no-template controls contributed 14.5 µl to the final pool while Subpool Group K with 20 libraries and 3 no-template controls contributed 11.5 µl to the final pool).

The Genomics Core Facility (University of Notre Dame, Notre Dame, IN) validated the final library pools using a combination of Qubit dsDNA High-Sensitivity Kit and Qubit 2.0 Fluorometer (Life Technologies, Grand Island, NY), Bioanalyzer 2100 instrument and Bioanalyzer DNA 7500 Kit or Bioanalyzer DNA High-Sensitivity Kit (Agilent Technologies, Santa Clara, CA) and Kapa Illumina Library Quantification qPCR assay (Kapa Biosystems, Wilmington, MA). Final molar concentration of the final library pool was based on the average of the values determined by Qubit and qPCR analysis. For each of the final library pools, a solution containing 6 pM denatured library and 3 pM denatured PhiX Control v3 (Illumina, Inc., San Diego, CA) was sequenced on Illumina MiSeq Sequencer operating MiSeq Control Software v2.5 with MiSeq flowcell (v3) and MiSeq Reagent Kit (v3). Sequencing format was 300 cycles for read 1, 8 cycles for index 1 read, 8 cycles for index 2 read and 300 cycles for read 2. Each library pool was run on a different MiSeq flow cell. MiSeq Reporter v2.5 Real Time Analysis (RTA) v1.18.54 performed base calling. Illumina Bcl2fastq v2.18 demultiplexed the RTA output and converted the data from bcl to fastq format. Finally, BaseSpace Broker v2.1 reported the files to BaseSpace Sequencing Hub.

Lastly, one set of samples was unfortunately mixed during the second stage of DNA extraction. While they were fully processed and contributed to the total raw reads observed, we removed them from any statistical analysis due to not being able to differentiate which replicate they were. This error reduced the number of experimental replicates for the treatment of cellulose nitrate filters with a pore size of 0.65 µm (CN 0.65 µm) extracted with the PCI-2 protocol from ten to four. We also lost one sample during DNA extraction from the treatment with the filter material of glass fibre with a pore size of 0.7 (GF 0.7 µm) extracted with PCI-1. The replicate could not be processed for sequencing and thus reduced the number of replicates for this treatment from ten to nine.

Sequencing adapters were removed and reads were quality filtered using Trimmomatic v0.32 (Bolger et al. 2014). An in-house Perl script was used to check for 100% match to both forward and reverse primers (Github: https://github.com/pfrender-laboratory/epps) and non-matching reads were discarded. Paired end reads were then merged using the -fastq_mergepairs command in USEARCH v8.1.1861(Edgar 2013) and singleton merged sequences (defined as a merged sequence observed only once across all replicates) were discarded. All replicates were pooled for OTU generation and OTU clustering was performed using the -cluster_otus command in USEARCH with two clustering thresholds of 97% and 95% similarity. Potential chimera sequences were also filtered by -cluster_otus command during the OTU clustering process as suggested by USEARCH. Centroid sequences from each OTU were taxonomically annotated using the Species Assignment Package v1.9.4 (Munch et al. 2008) against the NCBI nr database. OTUs were counted as metazoan if they received an assignment to the phylum level or lower. Specific parameters for each filtration step are shown in Suppl. material 2.

To explore the effect of method choice on a globally relevant NIS, we determined the detection of Styela clava in different technical replicates based on species assignments of OTUs. Styela clava, a solitary ascidian, is native to the NW Pacific. Over the last 100 years, it has become common in fouling communities in harbours and ports across the globe, including Europe, Australasia and North America (Goldstien et al. 2011). In the United Kingdom, it was first discovered in 1953 in Plymouth and then described as Styela mammiculata by Carlisle (1954). This introduction was attributed to the transport via warships after the end of the Korean war (Holmes 1968). It was first recorded in Southampton water in 1958 reaching population densities of 200 adults per square metre (Holmes 1968). Since then, the contemporary presence of Styela clava in Southampton water is well documented (Bishop et al. 2015). In each replicate, detection was considered positive if a single read or greater was assigned to Styela clava.

All statistical analyses were performed using the programme R, version 3.3.1.1 (R Core Team 2016). R scripts and required data files are provided on Figshare (https://doi.org/10.6084/m9.figshare.c.4191167.v1). Least squares regression models were used to test the additive effect of filter material (GF or CN), extraction type (CI, PCI-1 and PCI-2), pore size (0.5 μm, 0.7 μm and 1.2 μm) and the nested effect of filter material on pore size on PCR-2 DNA concentrations and OTU richness (both eukaryotic and metazoans only). As mentioned above, GF 0.7 µm and CN 0.65 µm were considered equivalent pore sizes for statistical analysis of the experiment. Richness of an experimental replicate was calculated by summing all OTUs with a single read or more. Due to the possible violation of homogeneity with pore size, resulting from the imbalance in factor levels (0.5 μm = 10, 0.7 μm = 53, 1.2 μm = 60), we also tested the models after removing the smallest pore size (0.5 μm) from the dataset. In addition to the two response variables of PCR-2 DNA concentrations and OTU richness, the variables of DNA concentration after PCR-1, the number of raw merged reads, number of merged reads after quality filtering, the number of OTUs clustered at 95% and the subset of metazoan OTUs clustered at 95% were tested for a correlation to establish an association between each step in the laboratory workflow and the outcome of methods’ choices. To evaluate the detectability of the NIS Styela clava across the experimental treatments, reads from the five OTUs assigned to Styela clava were pooled (Suppl. material 3) and detection was indicated by a single read or more. The effect of experimental treatments on detection was modelled using the same least squares regression model as above.

Raw fastq files have been deposited in the Short Read Archive (PRJNA395904). The two MiSeq runs produced a total of 45,600,412 raw reads for the experimental replicates and 41,522 reads for the no-template controls. After merging paired reads and quality filtering, a total of 16,144,713 reads remained for the experimental replicates and 837 reads for the no-template controls. Experimental replicates had an average of 125,153 ± 111,631 reads and no-template controls had an average of 42 ± 25 reads.

The number of unique reads summed across both runs was 996,242. After removal of singletons, 260,711 unique reads remained. Clustering at 97% and 95% resulted in 4,554 and 3,634 OTUs, respectively (Suppl. material 3). Based on taxonomic assignment of OTUs clustered at 97% and 95%, 2,112 and 1,519 were identified to at least the taxonomic level of phylum and represented 123 and 102 metazoan eukaryotes, respectively. Using clustering at 97% as the example, the no-template controls had reads assigned to 61 of the OTUs (Suppl. material 4). Of these, 28 OTUs had only one read assigned to a no-template control inclusive of both MiSeq runs. When comparing the other 33 OTUs, only two were consistently detected between the two MiSeq runs. As so few OTUs were consistently detected between runs and nearly half were represented by one read, we did not consider this evidence of contamination due to preparation error in the sampling or library preparation workflow, but rather evidence of the known tag-jumping or cross-talk errors generated on Illumina instruments (Schnell et al. 2015; Edgar 2016). As this error results in so few reads, we considered this error random and did not apply a correction factor to the data.

DNA concentration varied as a function of filter material, pore size and extraction method (Fig. 3a). Using PCR-2 DNA concentrations as the response variable for an estimate of DNA concentrations of eukaryotic environmental DNA (Fig. 3b), a significant effect of extraction type (P <0.001), pore size (P <0.001), filter material (P <0.001) and a nest effect of pore size on filter material (P = 0.035) were found (Table 1, Fig. 4a). Specifically, DNA concentration after PCR-2 ranged from 0.08 to 45.70 ng/μl (mean = 18.427, s.d. = 15.770) (Fig. 4a). Chloroform extractions yielded PCR-2 DNA concentrations from 0.08 to 30.6 ng/μl (mean = 4.314, s.d. = 6.429), with PCI-1 extractions producing 0.18 to 45.7 ng/μl (mean = 26.414, s.d. = 12.861) and PCI-2 extractions producing 0.16 to 45.5 ng/μl (mean = 30.022, s.d. = 12.323). DNA concentrations, after PCR-2, varied with corresponding pores sizes as well as with samples filtered with 0.5 μm filters showing 2.78 to 30.6 ng/μl (mean = 11.173, s.d. = 9.780), 0.7 μm resulting in 0.08 to 45.7 μm (mean = 15.849, s.d. = 15.453) and 1.2 μm filters resulting in 0.15 to 45.5 ng/μl (mean = 21.914, s.d. = 16.184). Samples processed with CN filters had PCR-2 DNA concentrations of 0.08 to 45.7 ng/μl (mean = 22.021, s.d. = 16.682) and samples processed with GF filters had PCR-2 DNA concentrations of 0.09 to 42.6 ng/μl (mean = 14.400, s.d. = 13.733). Models, from which we removed the smallest pore size to assess possible type I errors, showed significant differences across extraction method (P <0.001), pore size (P = 0.004) and filter material (P <0.001) and between the interaction of pore size nested within filter material (P = 0.033) (Suppl. material 5).

Figure 3.

Observed DNA concentrations or number of Operational Taxonomic Units by experimental treatment. Observed DNA concentrations after DNA extraction (a) indexed PCR (PCR-2) (b) and the number of estimated OTUs per experimental treatment clustered at 97% for all eukaryotes and (c) and metazoans (d). The upper and lower whiskers indicate the minimum and maximum point within 1.5 times the Interquartile Range extended from the 25^th and 75^th percentile, respectively. Colours indicate filter material where red is Cellulose Nitrate (CN) and blue is Glass Fibre (GF). The three extraction methods are abbreviated as in Figure 1.

Table 1.Download as CSV 

Statistical results from the least squares regression models with response variables of PCR-2 DNA concentrations, OTUs with 97% sequence similarity and Styela clava detection. Provided for each explanatory variable (extraction type, filter material, pore size and pore size nest in filter material) are the degrees of freedom (DF), Sum of squares (Sum sq), Mean squares (Mean sq), F value and p-value (P-value).

Response: PCR-2 DNA Concentration
	DF	Sum Sq	Mean Sq	F value	P value
Extraction type	2	17018.2	8509.1	111	< 0.001
Pore size	2	1231.1	615.5	8.06	< 0.001
Filter material	1	2883.3	2883.3	37.8	< 0.001
Pore size / Filter material	1	347.8	347.8	4.55	0.035
Residuals	116	8859.8	76.4
Response: OTU 97
Extraction type	2	1980701	990351	55	< 0.001
Pore size	2	378264	189132	10.5	< 0.001
Filter material	1	1152720	1152720	64	< 0.001
Pore size / Filter material	1	66101	66101	3.67	0.058
Residuals	116	2090261	18019
Response: Styela clava detection
Extraction type	2	4.148	2.074	10.595	<0.001
Pore size	1	0.177	0.177	0.905	0.343
Filter material	1	0.587	0.587	2.998	0.086
Pore size / Filter material	1	0.007	0.007	0.033	0.855
Residuals	123	24.074	0.196

OTU detection for eukaryotes and the subset of metazoan eukaryotes varied with method choice (Figure 3c, d). DNA extraction concentrations and PCR-2 DNA concentrations were correlated with OTU richness estimated for all eukaryotes and the subset of metazoan eukaryotes (Fig. 5). Significant effects were observed for the number of OTUs (clustered at 97% similarity) detected across extraction types (P <0.001), pore size (P = 0.010) and filter material (P <0.001) (Fig. 4b, c). However, a non-significant interaction effect of pore size nested in filter material (P = 0.058) was detected (Table 1, Fig. 4b, c). Models, from which the smallest pore size was removed to assess possible type I errors, showed comparable results in that OTU detection and significant differences amongst extraction types (P <0.001) and filter material (P <0.001), but there was a non-significant effect of pore size nested in filter material (P = 0.056) and a non-significant effect of pore size (P=0.061) (Fig. 4b) (Suppl. material 5). Results were similar for the subset of OTUs that were assigned to the metazoan phylum and clustered at 95% similarity (Fig. 4c, Fig. 5). Furthermore, each step of the workflow from DNA extraction to library preparation to bioinformatic choices for clustering threshold showed strong positive correlations (Fig. 5).

Based on the taxonomic assignment of OTUs, Styela clava was detected in 45 of 73 (62%) PCI replicates compared to 10 of 50 (20%) CI extraction replicates (Suppl. material 3). A significant effect (P = 0.001) of extraction method was found; the CI method was less effective than both phenol-based extractions (Table 1, Fig. 6).

Figure 4.

Observed DNA concentrations or number of Operational Taxonomic Units by experimental factor. Mean PCR-2 DNA concentrations (indexed libraries) (y-axis) for each set of explanatory variables including extraction type (A), filter pore size (B) and filter material (C). Extraction types include chloroform (Cl), phenol chloroform 1 (PCl 1) and phenol chloroform 2 (PCl 2) extractions. Filter materials included Cellulose Nitrate (CN) and Glass Fibre (GF). The upper and lower whiskers indicate the minimum and maximum point within 1.5 times the Interquartile Range extended from the 25^th and 75^th percentile, respectively.

Figure 5.

Pairwise correlations between non-independent response variables corresponding to quantified values taken during the analyses. Lower triangle consists of scatterplots with individual sample comparisons shown. The upper triangle shows the R correlation value corresponding with the lower triangle. The diagonal shows the distribution curve for the variable indicated in the corresponding column.

Figure 6.

Detection rate of Styela clava by extraction method. The proportion of positive detections for Styela clava across experimental treatments (left). Photographic image of Styela clava individual from Chichester Harbour, United Kingdom (Right). Colours are indicating the three different extraction methods which are abbreviated the same as in Figure 1.