eDNA metabarcoding as a new surveillance approach for coastal Arctic biodiversity

Abstract Because significant global changes are currently underway in the Arctic, creating a large‐scale standardized database for Arctic marine biodiversity is particularly pressing. This study evaluates the potential of aquatic environmental DNA (eDNA) metabarcoding to detect Arctic coastal biodiversity changes and characterizes the local spatio‐temporal distribution of eDNA in two locations. We extracted and amplified eDNA using two COI primer pairs from ~80 water samples that were collected across two Canadian Arctic ports, Churchill and Iqaluit, based on optimized sampling and preservation methods for remote regions surveys. Results demonstrate that aquatic eDNA surveys have the potential to document large‐scale Arctic biodiversity change by providing a rapid overview of coastal metazoan biodiversity, detecting nonindigenous species, and allowing sampling in both open water and under the ice cover by local northern‐based communities. We show that DNA sequences of ~50% of known Canadian Arctic species and potential invaders are currently present in public databases. A similar proportion of operational taxonomic units was identified at the species level with eDNA metabarcoding, for a total of 181 species identified at both sites. Despite the cold and well‐mixed coastal environment, species composition was vertically heterogeneous, in part due to river inflow in the estuarine ecosystem, and differed between the water column and tide pools. Thus, COI‐based eDNA metabarcoding may quickly improve large‐scale Arctic biomonitoring using eDNA, but we caution that aquatic eDNA sampling needs to be standardized over space and time to accurately evaluate community structure changes.


| INTRODUC TI ON
In the Arctic, climate change and marine invasions are expected to result in over 60% species turnover from present biodiversity with substantial impacts on marine ecosystems (Cheung et al., 2009).
Climate change is opening new waterways in the Arctic Ocean, resulting in greater shipping traffic (ACIA 2004;Arctic Council 2009;Guy & Lasserre, 2016). Predicted increases in shipping frequency and routes (Eguíluz, Fernández-Gracia, Irigoien, & Duarte, 2016;Miller & Ruiz, 2014;Smith & Stephenson, 2013), increased infrastructure development in ports (Gavrilchuk & Lesage, 2014), and associated chemical/biological pollution will place other ecosystem services at risk. Furthermore, the introduction of nonindigenous species (NIS) may displace native species, alter habitat and community structure and increase aquaculture and fishing gear fouling in estuaries and coastal zones (Goldsmit et al., 2018;Grosholz, 2002;Parker et al., 1999). Currently, the continuous monitoring needed to evaluate large-scale changes in coastal biodiversity and faunal assemblages in the Canadian Arctic is limited (Archambault et al., 2010), hindering risk management and ecosystem sustainability planning (Larigauderie et al., 2012).
Recent advances in the collection and analysis of environmental DNA (eDNA) provide a new complementary approach that can help to fill gaps in regional species distribution data left by logistically difficult traditional methods (e.g., bottom trawl, SCUBA diving) (Deiner et al., 2017), particularly in remote and otherwise challenging locations. eDNA allows for the detection of traces of DNA in water from macro-organisms (Thomsen, Kielgast, Iversen, Wiuf, et al., 2012). Collecting water samples for eDNA surveys could allow rapid sample collection, reduce the cost associated with data collection/shipping, and is less destructive because it does not require the manipulation of organisms (Lodge et al., 2012;Taberlet, Coissac, Hajibabaei, & Rieseberg, 2012). eDNA metabarcoding (i.e., high-throughput eDNA sequencing) can enable the identification of millions of DNA fragments/sample, providing a powerful approach to survey aquatic biodiversity. Repeated eDNA surveys could potentially be used to evaluate long-term biodiversity changes such as detecting native species loss and declines, NIS introductions and range expansions, and community structure changes. However, the detection of species using eDNA varies as a function of the population densities (Lacoursière-Roussel, Côté, Leclerc, & Bernatchez, 2016;Lacoursière-Roussel, Dubois, & Bernatchez, 2016;Mahon et al., 2013), life history traits, shedding rates (Lacoursière-Roussel, Rosabal, & Bernatchez, 2016;Sassoubre, Yamahara, Gardner, Block, & Boehm, 2016) local environmental conditions and technical approaches such as sequencing efforts and primer biases (Freeland, 2017;Pawluczyk et al., 2015).
Moreover, major concerns with eDNA metabarcoding, including its ability to accurately identify sequences to species (Chain, Brown, MacIsaac, & Cristescu, 2016) and the unknown ecological dynamics of eDNA in coastal ecosystems, need to be studied before marine biodiversity can be compared across spatial and temporal scales using this method.
Little is currently known about the efficacy of eDNA metabarcoding in surveying long-term variation in marine coastal biodiversity (Lim et al., 2016;Port et al., 2016;Thomsen & Willerslev, 2015).
Relative to freshwater ecosystems where more studies have been conducted, eDNA in coastal marine ecosystems is diluted into a much larger volume of water and exposed to pronounced hydrodynamics (e.g., tides, currents) and variation in abiotic conditions (e.g., salinity, temperature), which is likely to affect eDNA transport and degradation (Foote et al., 2012;Thomsen, Kielgast, Iversen, Møller, et al., 2012). In spite of these challenges, a recent study of horizontal spatial eDNA distribution in the Puget Sound (Washington, USA;O'Donnell et al., 2017) was successful in revealing fine scale distribution of species in these communities. In Arctic ecosystems, higher eDNA transport and diffusion is expected due to slower DNA degradation in cold-water temperatures, but no study has yet characterized aquatic eDNA distribution in this environment. Improving our understanding of the ecology of eDNA-the myriad of interactions between extraorganismal genetic material and its environment (Barnes & Turner, 2016)-in various ecosystems is fundamental to determining how eDNA can and cannot improve biodiversity research.
Our objective is to explore the potential of eDNA as a biodiversity monitoring approach to assist in rapid detection of coastal biodiversity shifts on large spatial scale in two Arctic coastal areas: Churchill and Iqaluit. These two Arctic commercial ports are expected to be particularly prone to biodiversity changes because they are among the top three ports in the Canadian Arctic with respect to vessel arrivals and associated ballast and/or hull fouling invasions risk (Chan, Bailey, Wiley, & MacIsaac, 2013). More specifically, we estimate the proportion of the Arctic biodiversity that can be identified at the species level with eDNA, and we then characterize the spatio-temporal distribution of eDNA with respect to water column depths, tide pools, and seasons.

| MATERIAL S AND ME THODS
The spatio-temporal eDNA distribution was characterized at three different depths in the water column, in tide pools, and between summer and fall seasons. Specifically, water samples were collected in 13 subtidal sites at three different depths (surface, mid-depth and deep water (i.e., 50 cm from the bottom), 12 tide pool sites within three intertidal areas (N = 4 sites/area) and 20 samples were collected at a single site from the shore approximately 2 m spaced along a transect (Figure 1). For the summer period (without ice cover), Churchill and Iqaluit were surveyed in 2015 between August 11-14 and August 17-22, respectively (hereafter called S20). To evaluate seasonal effects (Iqaluit only), the 20 samples at a single site were collected during fall (November 18, 2015) near shore from water that rose between ice pans at high tide (hereafter called F20).
Each sample (250 ml water) was collected using a Niskin bottle and then rapidly filtered in the field through a 0.7 μm glass microfiber filter (Whatman GF/F, 25 mm) using syringes (BD 60 ml; Kranklin Lakes, NJ, USA). Field negative controls (i.e., 250 ml distilled water) were filtered for every 10 samples. Filters were preserved at 4°C in 700 μl of Longmire's lysis/preservation buffer within a 2 ml tube for up to 3 weeks (Wegleitner, Jerde, Tucker, Chadderton, & Mahon, 2015) and then frozen at −20°C until DNA extraction. To reduce risk of cross-contamination during sampling and the filtration process, individual sampling kits were used for each sample (bottles and filter housing sterilized with a 10% bleach solution and new sterilized gloves, syringes, and tweezers). Each sampling kit was exposed to UV for 30 min. To reduce the risk of laboratory cross-contamination, procedures for eDNA extraction, PCR preparation, and post-PCR steps were all performed in different rooms and PCR manipulations were performed in a decontaminated UV hood. Samples from a specific port were all treated together, and the bench space and laboratory tools were bleached and exposed to UV for 30 min prior to processing the next port. Sites within a port were processed in a randomized order.

| eDNA extraction, amplification and sequencing
DNA was extracted using a QIAshredder and phenol/chloroform protocol (see Supporting Information Appendix S1). Negative control extractions (950 μl distilled water) were performed for each sample batch (i.e., one for each 23 samples). Two pairs of universal metazoan mitochondrial cytochrome c oxidase subunit I (COI) primers that have been developed and tested on a broad array of marine species were used to amplify eDNA from as many metazoan taxa as possible: the forward mlCOIintF (Leray et al., 2013) and reverse jgHCO2198 (Geller, Meyer, Parker, & Hawk, 2013) amplifying 313 bp (hereafter called COI1) and the forward LCO1490 (Folmer, Black, Hoeh, Lutz, & Vrijenhoek, 1994) and reverse ill_C_R (Shokralla et al., 2015) amplifying 325 bp (COI2).
The performance of the two selected primer pairs used in this study was previously tested on 104 zooplankton species and was validated on mock metazoan communities collected in Canadian ports by Zhang (2017). Based on a total of 13 COI primer pairs selected from the literature and tested, Zhang (2017) showed the efficiency of using multiple COI primer pairs in a single Illumina run to recover species by metabarcoding and detected 32% of species using COI1 and 49% of species using COI2. Here, the DNA amplification protocols for both primer pairs were optimized in vitro using 12 Arctic specimens and 12 potential invaders (i.e., annealing temperature gradient using DNA extracted from tissue samples; Supporting Information Table S1). The primer sequences and sequence databases were also evaluated in silico for their ability to detect native and potential nonindigenous Arctic metazoans. A list of recorded coastal Arctic metazoans was obtained by pooling all Arctic species databases that we had access to (N total = 897 metazoan identified at the species level; Fisheries and Oceans Canada Arctic Marine Invertebrate Database (Supporting Information Appendix S2), Archambault unpublished data, Cusson, Archambault, and Aitken (2007), Goldsmit, 2016;Goldsmit, Howland, & Archambault, 2014;K. Howland, P. Archambault, N. Simard and R Young, unpublished data, Piepenburg et al., 2011;Link, Piepenburg, & Archambault, 2013;López, Olivier, Grant, & Archambault, 2016;Olivier, San Martín, & Archambault, 2013;Roy, Iken, & Archambault, 2015;Young, Abbott, Therriault, & Adamowicz, 2016). Potential NIS invaders (N = 130 species) were targeted based on (1) screening level risk assessments and predictive species distribution models indicating they were high risk (Goldsmit et al., 2017), (2) their presence in ports connected to the Canadian Arctic, and/or (3) their presence in ballast waters and hulls of ships based on monitoring at Canadian Arctic ports (Chan, MacIsaac, & Bailey, 2015;Chan et al., 2012). Historical data include many Arctic F I G U R E 1 Geographical locations of the sampling port in the Canadian Arctic (map a) and the site distribution within Churchill (map b) and Iqaluit (map c). Subtidal areas are shown in white and the intertidal areas in light gray. Circles depict the water column sites, triangles are the tide pools sites and the squares are the S20 and F20 shore sampling sites

| Taxonomic identification
Forward and reverse sequences for each sample were trimmed using Trimmomatic 0.30 (Bolger, Lohse, & Usadel, 2014). FastQC version v0.11.3 was used to confirm the quality of the trimmed reads (Andrews, 2010). The Fastq quality scores were all well above 20 for the trimmed reads. Reads were then merged using FLASH v1.2.11 with a minimum overlap of 30 bp (Magoč & Salzberg, 2011).
"Orphan" reads with <30 bp of overlap between forward and reverse reads were discarded and only merged reads were used in the analy- Supporting Information Appendix S3). OTUs represented by a single read (singletons) were excluded and the identity between the representative sequences and the BOLD database was performed using vsearch (Rognes, Flouri, Nichols, Quince, & Mahé, 2016). For each phylum, proportion of the biodiversity assigned to the species level was obtained from the number of OTUs between 97-100% (similar to threshold used to assign species for sequences in the BOLD database) relative to those between 80-97% (i.e., below species level).

| Statistical analyses
Sampling effort is an important factor to consider in both traditional and eDNA biodiversity surveys. Two levels of port-specific sampling effort were explored: number of unique species per read (a measure of sequencing effort) and the number of unique species per sample (a measure of eDNA collection effort). For water column (surface, mid-depth and deep), tide pool and shore (S20 and F20) sampling locations, we plotted both read and sample rarefied accumulation curves to visualize whether or when a plateau was reached (which would indicate adequate sequencing and sampling effort to characterize all species). We also inspected the relative position of the read curve compared to the sample curve, as read curves lying above sample curves typically indicate spatial aggregation of species (Gotelli & Colwell, 2010), or in this case eDNA sequences.
These sampling effort analyses were performed in R 3.4.1 using the specaccum function in the vegan package.
All further statistical analyses were performed using R 3.0.3.
The spatial distribution of eDNA and the seasonal variability in the community composition was represented using Principal component analysis (PCoA) and tested using PERMANOVA (Anderson, 2001) after Hellinger transformation. Hellinger transformation was appropriate to deal with the large proportion of zeros and reduces the importance of large abundances (Legendre & Legendre, 1998) that could be due to the eDNA origin (e.g., capture of cell or mitochondria vs. extracellular DNA) or the amplification process. Species that mostly contributed to the dissimilarity/similarity between the treatments (depths and tide pools vs. water column) were identified using SIMPER analysis using the simper() function of the vegan package. Shannon diversity indices were calculated with the R package vegan. Analyses of variance (ANOVAs) were used to test whether species diversity, richness and log 10 (reads abundance) varied as a function of sampling location (i.e., water column and tide pools; sites included as a random variable) and water depths for each port separately using the lme() function of the NLME package (Pinheiro, Bates, DebRoy, & Sarkar, 2017) with sites included as a random variable (interactions between sites and depths could not be tested due to unique values per depth). The seasonal effect on read abundance (i.e., metazoan reads, see section taxonomic identification), Shannon diversity and species richness was evaluated using a Student's t test comparing the S20 and F20 samples in Iqaluit. Sørensen and Jaccard nonparametric estimates were calculated for location pairs using the SimilarityPair function of the SpadeR package in R (Chao, Ma, Hsieh, & Chiu, 2016) to test for the level of similarity in species composition between sampling locations and seasons.

| RE SULTS
After bioinformatics filtering (Supporting Information   Figure 2). For both ports, the sampling effort could have been increased to reveal additional species as the sample and read accumulation curves did not plateau (Supporting Information Figure S1). However, there was little evidence for spatial eDNA aggregation within a location as sample-based curves fell only slightly below read curves, and within 95% confidence intervals, at all locations. At the species level, the primer sets detected a total of ten phyla; including nine phyla for the COI1 primer set (44 Annelida

| Spatial eDNA distribution
For both ports, the community structure differed significantly between the water column and the tide pools, but the proportion of The Shannon diversity index was significantly greater in the water column than tide pools in Churchill (ANOVA: p = 0.002), but there was no significant difference in Iqaluit (p = 0.2; Figure 5). In Churchill, despite a significantly greater number of reads in tide pools than the water column (averaging 23,276 and 11,623 reads in tide pools and water column samples, respectively; p = 0.06), there was no significant difference in species richness between water F I G U R E 3 Principal component analysis depicting the community structure at the species level among sampling locations: water column (surface, mid-depth and deep water), tide pools (i.e., intertidal zone) and surface water collected in a single site in summer (i.e., S20) and in fall (F20) for both Arctic sampling ports (Churchill and Iqaluit). Ports were analyzed separately because each port was treated on a separate sequencing run F I G U R E 4 eDNA community differences between sampling locations (i.e., water column (surface, mid-depth and deep), tide pools) and seasons (summer S20 and Fall F20 Ann elid a P e c t i n a r i a g r a n u l a t a 8 % ocu lata 7% pe la gi ca 3% N a is b r e ts c h e r i 3 % Deep A r t h r o p o d a P s e u d o c a l a n u s n e w m a n i 3 4 % A c a r t i a l o n g i r e m i s 1 7 % ba la nu s 5% 10 more A n n e li d a P e c ti n a r ia g r a n u la ta 2

| Seasonal variation
The community structure varied significantly between the summer and fall sampling (  ing. Here, we present evidence that eDNA may be used to assess

| D ISCUSS I ON
Arctic biodiversity and show that, despite the cold and well-mixed F I G U R E 6 Relationship between the species richness detected using eDNA metabarcoding and the salinity of the water collected for the surface layer (R 2 = 0.85, black; circles: sampling water column and triangles: S20) and mid-depth samples (R 2 = 0.44, gray squares) and deep water (gray cross) F I G U R E 5 Boxplots comparing Shannon indices, species richness, and read abundances detected using eDNA metabarcoding for each sampling location (i.e., water column (surface, mid-depth and deep), tide pools and S20 and Fall20) in Churchill and Iqaluit. The lines inside the boxes represents the median values, the top and bottom of the boxes represent the 75% and 25% quartiles and outliers are shown using empty circles (i.e., any data beyond 1.5*IQR) environment, standardized eDNA approaches to biodiversity monitoring will need to consider local spatio-temporal variation.

| Taxonomic assignment challenges
The high congruence between historical Arctic data and eDNA samples ( and Rotifera were less likely to be detected than other groups such as Annelida (Figure 2). The use of eDNA metabarcoding may thus become a powerful approach to guide reference database improvement (e.g., 97% Rotifera OTUs were not identified at the species level). Moreover, groups such as Bryozoans, Nemerteans and Rotifera are currently not included in the historical Arctic Canada species records that we compiled, but they are important to coastal ecosystems and could be good indicators of biodiversity shifts caused by ice cover changes.
The eDNA metabarcoding method might thus be a good practical approach to evaluate the community changes of such species groups, even when poorly identified at the species level.
The better our knowledge of local species richness, potential invaders, and their corresponding genetic information, the more accurate our eDNA biodiversity monitoring methods will become. However, even when not assigned to species, the eDNA sequences detected here provide a sequence reference baseline that can be used to evaluate future species loss, new invasions, or other changes in community structure.
Once a taxon has been firmly identified by taxonomic experts and its barcode sequence has been deposited in GenBank or BOLD, eDNA might eventually reduce the need for large teams of expert taxonomists to carry out routine biodiversity monitoring. Yet, the routine application of metabarcoding for Arctic monitoring requires overcoming various limitations. For example, here the eDNA metabarcoding identified Acartia tonsa, a potential invader that has been previously recorded in the ecoregions of ports connected to Churchill (Chan et al., 2012). However, the current available COI sequences for Acartia tonsa form several distinct clades, some of which cluster with Acartia hudsonica, raising the possibility that the eDNA sequences assigned to A. tonsa actually belong to the native A. hudsonica. Thus, taxonomic expertise remains crucial for reducing biases of species distributions related to increasing use of large-scale eDNA metabarcoding.
Using two COI primer pairs, we increased the level of genetic polymorphism recorded at the species level, thereby increasing the resolution of the method for biodiversity monitoring (Deagle, Jarman, Coissac, Pompanon, & Taberlet, 2014). In addition to increasing the number of species detected, combining multiple primers may also reduce bias of eDNA dominance among species groups (e.g., dominance shift between Arthropoda and Annelida; Figure 2). Despite the fact that the amplification of COI is often desirable to differentiate species using DNA barcoding procedures (Che et al., 2012), the degree of universality for COI primers is relatively low and so combining multiple COI primer pairs as we did enabled monitoring a greater proportion of the diversity.
Further studies are, however, needed to evaluate how the combination of the primer sets may depict local species diversity. On the other hand, targeting genes with lower taxonomic specificity (e.g., 18S) could improve the detection of biodiversity shifts at higher levels (e.g., phyla level; see Bik et al., 2012;Deagle et al., 2014;Elbrecht & Leese, 2015).

| Spatio-temporal eDNA variation
Our results clearly show that metazoan eDNA distribution in Arctic coastal environments has significant temporal and spatial variation. The transport of eDNA may be substantially higher Vertical eDNA distribution in the water column may vary as a function of the life cycle of species, transport and settling advection (Turner, Uy, & Everhart, 2015) and complex hydrodynamic processes. In addition to wave action on eDNA mixing Port et al., 2016), our data support the idea that in estuarine conditions, such as in Churchill, the freshwater flowing from the river over long distances may contribute to increasing the diversity in the surface water layer (e.g., Deiner & Altermatt, 2014;Jane et al., 2015). Community changes related to eDNA composition thus need to integrate information on temporal variation in river discharge. The variability in the eDNA capture zone should therefore combine complex interactions between community changes and hydrodynamic models.
The dominance of Mollusca reads in tide pools is consistent with the observed species composition in these habitats (e.g., Goldsmit, 2016

| Arctic conservation biology
As contributions of sequences from identified specimens increase to databases such as BOLD, so too will the ability to track biodiversity changes over time at the species level with powerful methods such as eDNA metabarcoding (Gibson et al., 2014;Ji, Ashton, & Pedley, 2013;Taylor & Harris, 2012). In the Arctic, the development of cost-effective monitoring methods is essential for better protecting the integrity of important natural environments and endangered species and to ensure sustainable subsistence harvesting by aboriginal people, as well as recreational and commercial harvest by non-Aboriginals. Applying eDNA metabarcoding to assess biodiversity in remote coastal regions offers several advantages toward increasing the speed and accuracy with which we can amass biodiversity data. As part of this research project, local community members and permanently stationed northern research staff were trained in eDNA sampling techniques with the goal of enabling a network of communitybased monitoring. In this context, we optimized eDNA strategies for remote regions. We first used a syringe method for filtering samples (Deiner & Altermatt, 2014), which allows for sampling at multiple sites simultaneously and limits cross-contamination between samples as each sample can be processed with independent equipment. Moreover, the simplicity of this approach allows inexperienced collaborators to collect more eDNA samples per unit of time relative to standard practices of using an electric pump. Second, as storing and shipping frozen samples in remote regions is risky and often not possible, we used methods that allowed for DNA preservation at room temperature (Renshaw, Olds, Jerde, McVeigh, & Lodge, 2014). Lastly, the cost-effective extraction method increases the ability to process large number of samples.
By overcoming methodological issues and improving knowledge about the ecology of eDNA in coastal area, this project creates the opportunity for future monitoring of metazoan coastal diversity in highly vulnerable ecosystems such as Arctic commercial ports. The combined benefits of being able to identify large numbers of species including local species and potential invaders, assess a large number of phyla, the local habitat variability and together with the effectiveness of the eDNA method under ice cover, are likely to make eDNA metabarcoding an efficient complementary approach to detect large-scale Arctic coastal biodiversity changes. As the eDNA method progresses, the use of eDNA is likely to expand and become a catalyst for increased research on coastal biodiversity, ecosystem services, and sustainability, particularly in remote regions of the world such as the Canadian Arctic.
However, spatio-temporal dimensions need to be considered in standardizing and optimizing the assessment of marine biodiversity using eDNA. Bonin and three anonymous referees for their constructive comments on a previous version of the manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
Anaïs Lacoursière-Roussel is a conservation ecologist evaluating anthropogenic impacts on the dynamic of aquatic communities. All authors of this manuscript are interested developing and calibrating the eDNA method in the application of aquatic species distri-