eDNA captures depth partitioning in a kelp forest ecosystem

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.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 12S 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 Anacapa Toolkit (version 1.0) for quality control, amplicon sequence variant parsing, and taxonomic assignment using standard parameters. The Anacapa Toolkit sequence quality control and amplicon sequence variant (ASV) parsing module relies on cutadapt (version 1.16), FastX-toolkit (version 0.0.13), and DADA2 (version 1.6) as dependencies and the Anacapa classifier module relies on Bowtie2 (version 2.3.5) and a modified version of BLCA as dependencies. We processed sequences using the default parameters and assigned taxonomy using two CRUX-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 CRUX-generated 12S reference database supplemented with California Current Large Marine Ecosystem fish specific references to assign taxonomy using all available 12S reference barcodes to identify any non-fish taxa following the methods of Gold et al. using a Bayesian cutoff score of 60. Although CRUX relies on ecoPCR (version 1.0.1), blastn (version 2.6.0), and Entrez-qiime (version 2.0) as dependencies, we note that Bayesian cutoff scores are not directly analogous to percent identity from blastn. The BLCA classifier 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 Anacapa-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, Hypsypops rubicundus). Following decontamination, we converted the taxonomic tables into phyloseq objects (version 1.30.0) in R. 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 vegan (version 2.5-6) package in R. 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 vegan. We then ran a Mantel Test and non-metric multi-dimensional scaling (NMDS) using vegan 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 gradientForest 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 R2 importance, as there was a steep drop-off in R2 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 R2 importance using the broken stick method, as there was a steep drop-off in R2 importance after 0.40. Because of the additional nearshore vs. surf zone variable, the R2 importance values were higher in this model, resulting in a different threshold R2 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