This study builds on previous spatial comparisons of Arctic coastal metazoan communities in a Churchill port using both eDNA metabarcoding and conventional species collections (Lacoursière-Roussel et al., 2018; Leduc et al., 2019). Churchill is a key Arctic port in Canada and is situated on a large estuary, creating dynamic eDNA assemblages from marine water and freshwater communities (Lacoursière-Roussel et al., 2018; Sevellec et al., 2021). The estuary is located at the mouth of the Churchill River on Hudson Bay (Manitoba, Canada, 58°46′N, 94°11′W). We compared three different sampling approaches at different temporal scales each representing similar time commitments, which could be realistically implemented in long-term monitoring (i.e., sampling for 3–4 days distributed across consecutive days, consecutive months, consecutive years). The three temporal sampling strategies correspond to years (interannual: August 2015 and 2016; N=80), months (August 2015/2016, September 2017, October 2017, November 2017 and December 2017; N=158), and days (August 10^th, 11^th, 12^th 2016; N=18). There were 20 samples also collected in August 2017, as planned in our original design, but were compromised during the extraction process. All the samples were collected in the same geographic area, within 0.67 km of each other ( Supplementary Table S1 ). Although spatial distance influenced the eDNA community composition in this study, the relatively close proximity of the samples resulted in a much greater proportion of explained variance being attributed to temporal variation (PERMANOVA Temporal COI, R2 = 21.128, P < 0.001, PERMANOVA Spatial COI, R2 = 4.683, P < 0.001, PERMANOVA Temporal 18S, R2 = 35.158, P < 0.001, PERMANOVA Spatial 18S, R2 = 6.137, P < 0.001). All samples from the different temporal strategies were collected and processed following protocols detailed in Lacoursière-Roussel et al. (2018). In brief, we collected 250 mL of surface water (at ~1 m–2 m depth) which was filtered (0.7 μm, 25 mm diameter GFF) using a syringe. Filters were preserved in longmire buffer, placed on ice during field transport, and frozen at −20°C until DNA extraction.

DNA was extracted from filters using a QIAshredder and phenol/chloroform protocol as described in Lacoursière-Roussel et al. (2018). Technical variability of primer biases from high-throughput sequencing was reduced by using four different universal primer pairs (Goldberg et al., 2016; Kelly et al., 2019). Thus, eDNA was amplified by two pairs of mitochondrial cytochrome c oxidase subunit I (COI) primers (mlCOIintF/jgHCO2198 and LCO1490/ill_C_R) (Folmer et al., 1994; Geller et al., 2013; Leray et al., 2013; Shokralla et al., 2015) and two pairs of ribosomal gene 18S primers (F-574/R-952 and TAReuk454FWD1/TAReukREV3) (Hadziavdic et al., 2014; Stoeck et al., 2010). A one-step dual-indexed PCR approach with Illumina barcoded adapters was performed using 6 µl Qiagen Multiplex Mastermix, 4 µl diH₂0, 1 µl of each primer (10µM), and 3 µl of DNA. The PCR program consisted of an initial denaturation step at 95°C for 15 min, followed by 35 cycles at 94°C for 30 s, 54°C for 90 s (except for primer LCO1490/ill_C_R which was at 52°C for 90 s), and 72°C for 60 s, and a final elongation at 72°C for 10 min. Three PCR replicates were done for each sample and primer pair combination. After purification using Ultra AMPure beads and quantification by PicoGreen, the aliquots were pooled together in equal molar concentrations. Sequencing was carried out on three different flow cells (one per year) using an Illumina MiSeq at the Plateforme d’Analyses Génomiques (IBIS, Université Laval, Québec, Canada, www.ibis.ulaval.ca).

Meticulous care was taken and good practices implemented to avoid contamination risks. In the field, the sampling kits, composed of sterilized bottles, filters, gloves, syringes, and tweezers, were bagged and sealed and then exposed to UV for 30 min following assembly. Multiple field negative controls were performed for each sampling event, with 250 ml of sterilized distilled water, (three in 2015, four in 2016, and four in 2017), which were all treated as regular samples for the remaining manipulations until sequencing. In the laboratory, (i) eDNA extraction, PCR preparation, and post-PCR steps were performed in three separate rooms; (ii) all PCR manipulations were performed in a UV hood; (iii) all laboratory bench space and tools were treated with a 10% bleach solution and exposed to UV for 30 min before performing any manipulations. At the extraction step, additional negative controls were performed, with 950 µl of distilled sterile water (four in 2015, five in 2016, and eight in 2017), and were also integrated in the amplification and sequencing processes. At the amplification step, a negative PCR control was done for each sample and primer pair (i.e., 3 µl of sterilized distilled water) because barcodes were different for each sample. No positive amplification was visualized in the PCR negative control via 1.5% agarose gel electrophoresis. Finally, sequences from taxa deemed to be associated with sample contamination were removed ( Supplementary Table S2 ) following a decision process based on their relative abundance detected in the negative controls as detailed in Leduc et al. (2019) and Sevellec et al. (2021). In brief, contaminant taxa were removed if the total number of sequences detected in negative controls (field and extraction; N=28) was greater than 2% of the total number of sequences detected across all samples for a given genus or species. This contamination protocol aims to limit potential result bias by retaining dominant Arctic taxa that would be expected to contribute low levels of contamination (Drake et al., 2022).

Cleaning and classification of the sequences were performed with the MiSeq Control software v2.3 then with Barque version v1.5.2 (www.github.com/enormandeau/barque). In brief, raw forward and reverse reads were trimmed, merged, cleaned of chimeras, and denoised, and then clustered as operational taxonomic units (OTUs) with a threshold of 97% similarity. Finally, the COI sequences were annotated using BOLD (http://www.boldsystems.org) and the 18S sequences were annotated using SILVA (https://www.arb-silva.de/). Sequences from non-marine species (insects and most mammals), as well as those that could not be taxonomically assigned, were removed from our dataset. Because of the different filtration steps, and based on our contamination protocol, nine (COI) and seven (18S) samples collected in September and October 2017 were removed from the data set. For the year, month, and day sampling designs, we successfully blasted a per sample average of 587.4, 565.5, and 411.6 sequences for COI1, 386.6, 547.3, and 253.6 sequences for COI2, 1,575.1, 1,043.9, and 239.6 sequences for 18S1, and 632.0, 517.0, and 135.8 sequences for 18S2, respectively. Because the results obtained between the gene primers COI1 and COI2 data as well as 18S1 and 18S2 data were highly similar ( Supplementary Figure S1 ), reads were summed for each gene. For results below, COI includes both COI1 and COI2 and 18S includes both 18S1 and 18S2.

Coastal marine eDNA communities in Churchill Port were compared across the three temporal sampling strategies hereafter refer to as years (interannual: August 2015 and 2016), months (August 2015/2016, September 2017, October 2017, November 2017, and December 2017) and days (August 10^th, 11^th, and 12^th, 2016). To determine the effect of sampling effort, genus-level rarefaction curves were created with BiodiversityR and ggplot2 packages using R (R Core Team, 2013) using relative abundance without prior transformation for the three temporal communities. Since samples were available from two different years in the August period for the monthly sampling strategy, this analysis was done using two subgroups to avoid overestimating biodiversity: (i) August 2015 with the months of September, October, November, and December in 2017, (ii) August 2016 with the months of September, October, November, and December in 2017.

To investigate the differences among eDNA communities from the various temporal sampling strategies, two analyses were performed using a genus matrix after Hellinger transformation with the R vegan package: PERMANOVA analyses (number of permutations = 10,000) and principal coordinate analyses (PCoA) were performed on a Bray–Curtis distance matrix representing (dis)similarity among eDNA communities. Furthermore, to confirm that the significant results of the PERMANOVA and PCoA analyses were not driven by variations in group dispersions, multivariate homogeneity was evaluated using the betadisper function (R vegan package). No significant differences in dispersion homogeneity were observed across temporal strategy groups (betadisper P>0.01). To examine the taxonomic composition of each temporal sampling strategy, three tree maps were generated to represent yearly, monthly, and daily strategies at the genus level, using Hellinger-transformed data and the R treemap package.

To visualize similarities among eDNA communities, Venn diagrams were created using InteractiVenn (Heberle et al., 2015). These were complemented by Similarity Percentage (SIMPER) analyses performed with PAleontological STatistics software (PAST 4) ((Hammer et al., 2005) to evaluate dissimilarities of temporal sampling strategies on overall average dissimilarity index (OAD index). This index ranged from 0 (samples have identical taxa composition) to 100 (the samples do not share any taxa). The SIMPER analysis also identified the taxa that explain similarity (similarity percent) or dissimilarity (contribution percent) within the temporal strategies group. In this study, the SIMPER analyses were based on Hellinger transformed data. Finally, to determine the differences between taxonomy levels, temporal variation was contrasted for each phyla using time-series analysis (R vegan package) on the most abundant representative genus within each (i.e., annelida - Nais, chordata - Halocynthia, cnidaria - Aurelia, mollusca - Mytilus, porifera - Spongilla). More precisely, the distribution of Hellinger transformed data for each representative genus was plotted separately for each temporal sampling strategy. To compare if the relative abundance of shared genera varied among temporal sampling strategies, a non-parametric Kruskal–Wallis test was performed using Hellinger transformed data.

A total of 206,492 (COI) and 151,433 (18S) aquatic metazoan sequences were obtained after trimming the raw dataset composed of 148 (COI) and 150 (18S) eDNA samples from the Churchill port (for more details, see Supplementary Table S3 ). For both the genes, the accumulation curves for the year and month sampling strategies reached a plateau, whereas the eDNA community based on daily samples were at the end of the semi-logarithmic curves for each primer ( Figure 1 ). At the same given level of sampling effort, the monthly strategy yielded the highest number of taxa, followed closely by the yearly strategy, whereas daily sampling resulted in the lowest number of taxa.

The PCoA results were significantly different between years (i.e., between August 2015 and 2016; PERMANOVA year, P < 0.001 for both primers), but the samples from 2015 and 2016 clustered side by side, especially with the COI primers ( Figure 2 ). The monthly eDNA communities formed distinct clusters (PERMANOVA month, P < 0.001) that showed a gradual seasonal change ( Table 1 ; Figure 2 ). In the opposite way, no distinct communities were observed among the three consecutive days ( Table 1 ). The PCoA revealed scattered sample data points along axis-2 for the COI primers and along axis-1 for the 18S primers, suggesting high stochasticity and a lack of clear community structure at the daily temporal scale ( Figure 2 ).

Despite differences in the temporal sampling design, the eDNA communities shared similar dominant (most abundant) genera ( Figure 3 for COI and Supplementary Figure S3 for 18S). Using the COI gene, the annelida Pectinaria, Nais, Chaetogaster, and Harmothoe were highly represented in the three temporal communities. Other commonly and highly detected taxa were the arthropods Acartia and Pseudocalanus, the mollusks Mytilus and Macoma, the cnidarians Aurelia and Hydra, and the chordate Delphinapterus. Similar patterns were observed with the 18S gene, where dominant genera included the mollusks Macoma and Mytilus, the cnidarian Laomedea, and the bryozoan Alcyonidium, all consistently present across the eDNA communities.

Consistent with the above previous results, the highest proportion of genera were shared between the yearly eDNA communities (COI: 63.9%; 18S: 55.1%). In contrast, a much lower proportion of shared genera were detected among the 5 months (COI: 25.1%; 18S: 16.7%) and the three consecutive days (COI: 48.3%; 18S: 38.7%) with both genes ( Figure 4 ). Results from SIMPER were generally consistent with the Venn results. For COI, the overall average dissimilarity index was lower based on the yearly sampling strategy (63.49%) than for monthly (67.05%) and daily (64.71%) eDNA sampling strategies, indicating greater community similarity across years compared with months and days ( Table 2 ). For the 18S gene, the overall average dissimilarity index was lower for the daily communities (61.41%) compared with the yearly (74.24%) or monthly communities (75.24%). For both COI and 18S, the overall average dissimilarity index decreased from September 2017 to December 2017, indicating that winter months presented a more stable eDNA community than summer or fall eDNA communities.

For the COI primers, the three most dissimilar genera between years included the mollusk Mytilus (90.48% similarity), the copepod Pseudocalanus (93.17%), and the annelid Pectinaria (93.52%) ( Supplementary Table S4 ). Across months, the most dissimilar genera for COI were the annelid Nais (94.56%), the bivalve Mytilus (94.78%), and the annelid Chaetogaster (95.19%) ( Supplementary Table S5 ). For the daily sampling strategy, the percentage of similarity for the three most dissimilar taxa was 95% for the annelid Pectinaria, 95.84% for the copepod Acartia, and 96.14% for the annelid Galathowenia ( Supplementary Table S6 ). The most dissimilar genera detected with the 18S primers are presented in the Supplementary Tables S7–S9 .

For five highly detected genera representing different phyla, the relative abundance varied within and among the different temporal sampling strategies ( Figure 5 ). All the genera had significant variation in abundance between months (Kruskal–Wallis test, P < 0.01, Table 3 ). Three genera, Nais, Aurelia, and Mytilus, displayed variation in abundance among years for the both gene primers ( Table 3 ). Although no genera showed significant variation in abundance between days with the 18S primers, differences were observed for the genera Aurelia and Mytilus with the COI primers ( Table 3 ; Figure 5 ). Interestingly, the temporal variation in relative abundance of these five genera was similar with the COI and 18S combined gene primers ( Figure 5 ). For example, the cnidarian Aurelia was detected in August 2015 and 2016 as well as from September to November 2017, but not detected in December for either of the genes. Other genera with patterns of note included the tunicate Halocynthia which was characterized by detection in only a few samples forming a succession of eDNA peaks in August 2015 and 2016. The relative abundance of the freshwater sponge, Spongilla, tended to increase from August to December. Interestingly, the 1-week difference between the sampling of year temporal strategy (6 August 2016) and the daily strategy (10–12 August 2016) resulted in a rapid increase of annelid Nais relative abundance.

Identifying optimal eDNA sampling frequencies and time periods for a given system will increase the detection rates of local communities. Annual recurrence in eDNA communities combined with monthly sampling indicates potential for optimizing sampling during peak periods of marine life activity. For well-studied taxa, sampling at an annual scale during optimal periods could substantially increase knowledge of phenology and track how life history events and populations of individual taxa and/or biological communities may change over longer time scales with, for example, changing climatic conditions or increased industrial development. In order to gauge such optimal periods, monthly or even weekly eDNA surveys over potential periods of interest during the year would be recommended.

Efficient sampling strategies are essential to acquiring eDNA datasets that are as robust and informative as possible. Here, we observed that sampling at multiple time points provides a more holistic view of Arctic estuarine communities. This study suggests that, when using eDNA metabarcoding, monthly sampling is likely the most optimal strategy to effectively capture a comprehensive portrait of estuarine communities. Ideally, such a sampling strategy can be studied over multiple years to follow annual patterns of eDNA communities. However, if time and resources are limited, repeat sampling at a consistent time of year could offer a viable alternative to monitoring for long-term changes in coastal communities, based on our observations of annual stability in a large component of the eDNA community. A better understanding of the biological and physical factors affecting eDNA presence for taxa of interest is fundamental for optimizing eDNA sampling and resources, achieving more accurate predictions of eDNA communities.