DOI: 10.5068/d18h47
title: eDNA captures depth partitioning in a kelp forest ecosystem

description: Environmental DNA (eDNA) metabarcoding is an increasingly important tool
 for surveying biodiversity in marine ecosystems. However, the scale of
 temporal and spatial variability in eDNA signatures, and how this
 variation may impact eDNA-based marine biodiversity assessments, remains
 uncertain. To address this question, we systematically examined variation
 in vertebrate eDNA signatures across depth (0 m to 10 m) and horizontal
 space (nearshore kelp forest and surf zone) over three successive days in
 Southern California. Across a broad range of teleost fish and
 elasmobranchs, results showed significant variation in species richness
 and community assemblages between surface and depth, reflecting
 microhabitat depth preferences of common Southern California nearshore
 rocky reef taxa. Community assemblages between nearshore and surf zone
 sampling stations at the same depth also differed significantly,
 consistent with known habitat preferences. Additionally, assemblages also
 varied across three sampling days, but 69% of habitat preferences remained
 consistent. Results highlight the sensitivity of eDNA in capturing
 fine-scale vertical, horizontal, and temporal variation in marine
 vertebrate communities, demonstrating the ability of eDNA to capture a
 highly localized snapshot of marine biodiversity in dynamic coastal
 environments.
description: We conducted our study at Leo Carrillo State Beach, Malibu,
 California, USA (34.0446° N, 118.9407° W). We sampled on three successive
 days (September 24 to September 26 in 2018) to test for temporal stability
 of spatial variation in eDNA signatures. On each day, we sampled at the
 highest tide of the mixed semidiurnal tide to minimize the impact of tidal
 variation on sampling. To ensure that results reflected
 variation in spatial sampling, rather than time, we synchronized watches
 and worked in multiple teams to simultaneously sample from five depths on
 SCUBA along a vertical transect in a kelp forest ~140 m from shore: 0 m
 (at the ocean surface), 1 m, 5 m, 9 m, and 10 m (just above the sea
 floor). We also sampled a sixth station along the shore in the surf zone,
 where we collected samples approximately 1 m below the water surface where
 depth was approximately 2 meters. At each location, we collected
 triplicate seawater samples using one-liter enteral feeding bags
 (Kendall-Covidien – 702500) following the methods of Curd et al. After
 returning to shore, we immediately gravity filtered all samples through
 0.22 μm Sterivex filters to isolate eDNA. We similarly filtered one liter
 of distilled water as a negative template control. Upon completion of
 filtration, we stored the dried filters at -20˚C until extraction 48 hours
 later. No permits were required for sampling seawater in state waters
 outside of designated protected areas. We extracted the
 DNA from the filters at UCLA using the Qiagen DNEasy Blood and Tissue kit
 (Qiagen, Valencia, CA, USA). To maximize eDNA recovery, we employed
 modifications made by Spens et al., adding proteinase K and ATL buffer
 directly to the filter cartridges before overnight incubation in a
 rotating incubator at 56˚C. We amplified the extracted eDNA using the
 <em>12S</em> MiFish Universal Teleost (MiFish-U) and MiFish
 Elasmobranch (MiFish-E) primers with linker modifications for Nextera
 indices. Though the primers target teleost fish and elasmobranchs, they
 can also amplify other vertebrate species such as birds and mammals. PCR
 amplification and library preparation was conducted following the methods
 of Curd et al.. After library preparation, we sequenced the library on a
 NextSeq at the Technology Center for Genomics & Bioinformatics
 (University of California, Los Angeles, CA, USA) using Reagent Kit V3 with
 30% PhiX added to the sequencing run. We processed the
 resulting sequencing data using the <em>Anacapa Toolkit</em>
 (version 1.0) for quality control, amplicon sequence variant parsing, and
 taxonomic assignment using standard parameters. The <em>Anacapa
 Toolkit</em> sequence quality control and amplicon sequence variant
 (ASV) parsing module relies on <em>cutadapt</em> (version
 1.16), <em>FastX-toolkit</em> (version 0.0.13), and
 <em>DADA2</em> (version 1.6) as dependencies and the
 <em>Anacapa classifier</em> module relies on
 <em>Bowtie2</em> (version 2.3.5) and a modified version of
 <em>BLCA</em> as dependencies. We processed sequences using
 the default parameters and assigned taxonomy using two
 <em>CRUX-</em>generated reference databases following the
 methods of Gold et al.. We first assigned taxonomy using the California
 Current Large Marine Ecosystem fish specific reference database. Second,
 we used the <em>CRUX</em>-generated <em>12S</em>
 reference database supplemented with California Current Large Marine
 Ecosystem fish specific references to assign taxonomy using all available
 <em>12S </em>reference barcodes to identify any non-fish taxa
 following the methods of Gold et al. using a Bayesian cutoff score of 60.
 Although <em>CRUX </em>relies on <em>ecoPCR
 </em>(version 1.0.1), <em>blastn</em> (version 2.6.0),
 and <em>Entrez-qiime</em> (version 2.0) as dependencies, we
 note that Bayesian cutoff scores are not directly analogous to percent
 identity from <em>blastn. </em>The <em>BLCA
 </em>classifier<em> </em>incorporates alignment metrics,
 including percent identity and percent overlap, into the underlying
 Bayesian model which then returns the Bayesian cutoff score metric as a
 measure of confidence for each taxonomic rank for a given ASV.
 The resulting <em>Anacapa</em>-generated taxonomic
 tables were transferred into R for further processing 
 (https://datadryad.org/stash/share/aMH1xTddyGgAhaWYoV3kmmmWgqCzv6Lt9YtU9s4F6NA). We then decontaminated the taxonomic tables using methods developed by Kelly et al. and McKnight et al. as implemented in Gold, which removes sequences from index hopping and negative controls and conducts a site occupancy model to identify true rare sequences. We also manually removed sequences for species with taxonomic assignments for non-marine taxa (e.g. terrestrial mammals) in R. We then merged ASVs by summing reads by assigned taxonomy (e.g. summed all sequences reads from the 7 ASVs that assigned to Garibaldi, <em>Hypsypops rubicundus</em>). Following decontamination, we converted the taxonomic tables into <em>phyloseq</em> objects (version 1.30.0) in <em>R.</em> We analyzed the eDNA signatures across depth and across nearshore vs. surf zone habitats. We analyzed differences in eDNA across depth using only the nearshore signatures, excluding the samples collected in the adjacent surf zone. To examine species richness across depths, we conducted ANOVA and post-hoc Tukey tests using eDNA read counts. We then transformed eDNA read counts to eDNA index scores, which better correlates to abundance, following the methods of Kelly et al.. The eDNA index calculation standardizes eDNA abundance across samples and across taxa. To calculate the eDNA index values, we first calculated the relative abundance of each taxa in each sample. We then divided the relative abundance of each taxa by it’s maximum observed abundance across all samples to standardize the read counts per species per sample. This results in an index that ranges from 0 to 1 for each species where a value of 1 corresponds to the sample with the greatest relative abundance observed for that species. Although these eDNA index values allow for direct comparisons of relative abundance within individual species, allowing for direct comparisons in relative abundance by depth, location and/or time, it cannot be used for comparisons among different species. See Kelly et al. for more detail. To analyze the importance of sampling depth on eDNA vertebrate community composition, we conducted a PERMANOVA test using the <em>vegan</em> (version 2.5-6) package in <em>R</em>. The PERMANOVA was run using Bray-Curtis dissimilarity and the model eDNA_Index ~ Depth + Day + Replicate. We also ran a multivariate homogeneity of group dispersions test using the betadisper function and Bray-Curtis dissimilarity using <em>vegan</em>. We then ran a Mantel Test and non-metric multi-dimensional scaling (NMDS) using <em>vegan </em>on Bray-Curtis dissimilarities to assess community composition differences across the depth gradient. We further analyzed vertical depth community composition by generating a gradient forest model using the <em>gradientForest</em> package (version 0.1-17) using 500 runs. The environmental variables in the vertical depth gradient forest model included sampling depth, sampling day, and replicate. We then extracted the taxa with the highest model performances and plotted their eDNA index values across depth. Using the broken stick method, we defined highest model performance as  ≥ 0.35 R<sup>2</sup> importance, as there was a steep drop-off in R<sup>2</sup> importance after 0.35. To analyze differences between the kelp forest and surf zone, we ran Welch t-tests, PERMANOVA and betadisper tests, only including samples taken at 1 m depths from nearshore and surf zone habitats. We used the Bray-Curtis dissimilarity for both tests and the model eDNA_Index ~ Habitat (nearshore vs. surf zone) + Depth + Day + Replicate for the PERMANOVA. We also ran an additional gradient forest model with the environmental variables sampling depth, nearshore vs. surf zone, sampling day, and replicate for eDNA index scores from all stations. We then extracted the top performing taxa and plotted their eDNA index distributions, defining highest model performance as ≥ 0.40 R<sup>2</sup> importance using the broken stick method, as there was a steep drop-off in R<sup>2</sup> importance after 0.40. Because of the additional nearshore vs. surf zone variable, the R<sup>2 </sup>importance values were higher in this model, resulting in a different threshold R<sup>2</sup> value than the one in the depth gradient forest model. To test whether vertical and horizontal variation in eDNA signatures were consistent over time, we compared species richness across sampling days in an ANOVA framework, looking at both total community diversity as well as the eDNA index abundances. The linear models used for the eDNA index ANOVA tests were eDNA_Index ~ Station + Day + Station:Day.