Benthic 0.1m² Van Veen grab samples were collected in July from the Canadian Coast Guard ship Sir Wilfrid Laurier during 2018 and 2019 cruises in the Bering Sea and Chukchi Sea ( Figure 1 ). M. calcarea were obtained from 11 sampling stations with 7–10 clams selected per station per year (n = 190) for dissection ( Table 1 ). Clams were stored in collection bags by station and frozen at -20°C prior to dissection. Dissections occurred in 2020 between the months of January and October. The whole wet weight of collected clams (including shell) and length, height, and width were recorded. Clams were placed on an individual weighing boat that was changed for every sample. Stomach and digestive gland (hereafter “stomach”) tissues were excised from adductor muscles, mantle, foot, and gills (outer surface was not rinsed prior to excision). Dissection tools were cleaned with 10% bleach solution and rinsed with water between samples. Precautions were taken throughout to prevent cross-contamination of PCR products. Only stomach tissues were used to ascertain diatom prey of Macoma clams. DNA was obtained from clam stomachs (and their contents) using a DNeasy PowerSoil Pro Kit (Qiagen) and eluted in 50µL of Solution C6. The DNeasy PowerSoil Pro Kit protocol contains a vortex step with beads that homogenize the stomach contents prior to DNA extraction. A single extraction protocol was performed on each stomach. Most clams (n = 122 of 190) had stomach samples that were < 250 mg of tissue, the maximum amount suggested in the DNeasy PowerSoil Pro Kit protocol. Of the remaining 68 clams, most of the stomach tissue and associated content was subsampled and included in the DNA extraction protocol. The portion of the stomach that was not included in the DNA extraction protocol was largely extraneous tissue. Extraction blanks were included every 23 samples.

We performed DNA metabarcoding on two genes, the 18S V9 SSU ribosomal RNA gene (18S rDNA) and the large subunit of ribulose bisphosphate carboxylase (rbcL) gene located in the chloroplast genome. We selected the 18S gene because it is broadly amplified in eukaryotes and its wide use in the literature allows for cross-study comparisons as well as a richer reference database from public repositories. We also amplified the rbcL gene to focus on organisms with chlorophyll, which are the presumed primary diet of Macoma clams (Lovvorn et al., 2005).

We followed the 18S Illumina Amplicon Protocol adapted from the Earth Microbiome Project (Amaral-Zettler et al., 2009) that amplifies ~260 base pairs (bp; +/- 50 bp). A Macoma clam blocking primer was designed and included in reactions to limit amplification of host DNA during PCR (Vestheim and Jarman, 2008). Macoma balthica, M. nasuta, and M. secta 18S sequences were obtained from the National Center for Biotechnology Information (NCBI) GenBank and the blocking primer was constructed at a region with conserved sites among species with a C3 spacer at the 3’ end (ClamBlock_1391F: GCCCGTCGCTACGACCGATTGT[I][I][I][I][I]TTAATGAGCGAT[SpC3]). No sequence data was available for M. calcarea. The Earth Microbiome Project protocol amplifies 18S in a single PCR step. Eur1391F and EurBr primers contain the Illumina adapter sequence, 18S locus specific primer, with the EurBr_R primer fitted with a unique 12 bp Golay barcode (Thompson et al., 2017). PCR amplifications were performed in 25µL volumes containing 4µL of genomic DNA, 1X Phusion High-Fidelity Master Mix with HF buffer (Thermo Scientific), 4% bovine serum albumin (BSA), 0.4µM Euk1391F primer and 0.6µM ClamBlock_1391F with 0.4µM EurBr primer added separately to each well. A negative control was included in each PCR amplification. Thermocycler conditions were 98°C for 30 s, 35 cycles of 98°C for 10 s, 65°C for 15 s, 57°C for 30 s, 72°C for 40 s, and concluded with 72°C for 5 min. PCR amplifications were performed in triplicate. Products were visually inspected on an 1.5% agarose gel to verify successful amplification and pooled by sample.

We adapted the protocol by Vasselon et al. (2017) for amplification of 312 bp of rbcL. Specifically, we adapted the primers to conform to Illumina adapters for library preparation on an Illumina platform using a two-step PCR protocol. Briefly, locus-specific primers Diat_rbcL_708F_1, Diat_rbcL_708F_2, Diat_rbcL_708F_3, rbcL_R3_1 and rbcL_R3_2 were fitted with Illumina Nextera overhangs on the 5’ end (F: TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG; R: GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG). Initial PCR reactions were carried out in 25µL volumes containing 4µL of genomic DNA from excised stomachs, 1X Phusion High-Fidelity Master Mix with HF buffer (Thermo Scientific), 4% Bovine serum albumin (BSA), 0.4µM primer pool Diat_rbcL_708F and 0.4µM primer pool Diat_rbcL_R3. Thermocycler conditions were 98°C for 30 s, 30 cycles of 98°C for 10 s, 55°C for 30 s, 72°C for 40 s, and concluded with 72°C for 5 min. PCR amplifications were performed in triplicate. Products were visually inspected on an 1.5% agarose gel to verify successful amplification. PCR products were pooled by sample, and excess dNTPs and primers removed with ExoSAP-IT (USB Corporation, Cleveland, OH). The second PCR reactions incorporated Nextera indices and adaptors and were carried out in 15µL volumes containing 1µL pooled PCR product, 1X Phusion High-Fidelity Master Mix with HF buffer, and 1.0µM each of P5 and P7 adaptors. Thermocycler conditions were the same as above for rbcL except the number of cycles was reduced to 20 and the annealing temperature was increased to 60°C.

PCR products for 18S and rbcL marker were quantified using a Broad Range Quant-iT dsDNA Assay Kit (Invitrogen, Carlsbad, California, USA). 18S and rbcL sample libraries with unique barcodes or index pairs were pooled in equimolar concentration (sets of 48 samples for 18S and 96 samples for rbcL) by locus, purified with a 1.8x concentration of AMPure XP beads (Beckman Coulter, Inc.) to remove excess dNTPs and primer, then quantified using a High Sensitivity Quant-iT dsDNA Assay Kit (Invitrogen, Carlsbad, CA). Pooled libraries (8–10pM library and 30% PhiX) for each locus were sequenced on an Illumina MiSeq (250 bp pair-end reads) at the U.S. Geological Survey Alaska Science Center. Extraction blanks and PCR negatives failed to generate libraries that could be sequenced. Sequences and metadata are deposited in the NCBI database under BioProject (PRJNA1190063; BioSample accessions SAMN45024226-45024419). Detailed sample information is available in Sonsthagen et al. (2024).

Reads were trimmed with BBDuk (Bushnell, 2025) using a modulus trim of 5 and a Phred-scaled minimum quality of 10. Read pairs were then merged with vsearch v. 2.14.1 (Rognes et al., 2016), requiring an overlap of 100 (rbcL) or 50 (18S) and a maximum difference in the overlap region set to 6%. For the 18S locus, the “fastq_allowmergestagger” option was set since many amplicons are expected to be shorter than the cycle number. For 18S, both merged read pairs and the first read of unmerged read pairs were analyzed, whereas unmerged rbcL read pairs were discarded. Primers were not trimmed because the primers were not degenerate and thus did not affect the distribution of alignment scores. Each locus was then clustered within each sample with Swarm v. 2 (Mahé et al., 2015) with the “fastidious” option invoked to link adjacent low-abundance swarms. Small clusters of fewer than ten reads in a given sample were set to zero to reduce potential contamination or tag jumping errors.

The resulting operational taxonomic units (OTUs) were assigned a taxonomy using the least common ancestor (LCA) approach. OTUs were first aligned to the nucleotide database of the NCBI with BLAST+ v. 2.10.0 (Camacho et al., 2009) with low complexity masking disabled and a maximum of 25 matches reported. The bit score was used as the scoring metric, and the assigned taxon was the lowest common ancestor of all taxa within 3% of the highest bit score, provided that the bit score was at least 400 for rbcL or 150 for 18S. Alignments to accessions with no established taxonomy (e.g., uncultured, environmental, or metagenomic sequences) were excluded. The LCA algorithm was then evaluated for each set of matched accessions at the standard ranks of species, genus, family, order, and class, as listed in the NCBI taxonomy database (accessed April 1, 2021). The cluster size was then considered to be the abundance of the OTU in the sample, and the taxon abundance was summed across all OTUs with the same assignment.

The LCA method is useful when well-scoped, discrete reference databases are not available, as is typically the case for these markers in marine environments (see Table 1 in Weigand et al., 2019, for example). LCA uses relative alignment-score distributions independently for each OTU, rather than probabilistic training on a discrete reference database that allows estimation of assignment confidence at species and genus level; otherwise these OTUs were demoted to the next taxonomic rank for re-evaluation. Specifically, for rbcL, species-level assignments required the average percent identity of matches to be 97% or greater and the maximum bit score of matches to be 500 or greater, whereas genus-level assignments required an average percent identity of 93% and a maximum bit score of at least 450. For 18S, species-level assignments required an average percent identity of 97% and genus-level assignments required an average percent identity of 93% (no minimum bit score was required due to the variation in amplicon size). Given that we are not aware of studies quantifying intraspecific and intergeneric variation at barcode loci in these taxa, these thresholds are necessarily somewhat arbitrary but were informed by the overall distribution of best alignment scores for OTUs at each locus (See Supplementary Figure S1 ).

OTUs that were unlikely to be considered prey items (e.g., mammals) and known parasitic dinoflagellates and apicomplexans were filtered from the data set. Read counts were normalized to counts per million (cpm). We applied a 1% threshold to taxa within samples to account for inherent errors associated with amplification sequencing and alignment methods that make it difficult to distinguish among low-abundance prey, environmental background, amplification artifacts, and assignment errors (Deagle et al., 2019).

To assess the relative frequency of OTUs with the same taxonomic assignments across all sampling stations, we utilized heatmaps as a visualization tool. For the 18S gene data, we filtered the dataset to include only the top 15 taxa based on normalized read counts, due to the many low-abundance taxa assigned, the detection of which is expected to be highly stochastic, and which contribute minimal dietary information. This filtering allowed us to focus on the most abundant taxa, providing a clearer picture of the dominant prey items. The heatmaps were generated using normalized read counts, with color gradients representing the relative abundance of each OTU with similar taxonomic assignments across each station. This approach facilitated the identification of spatial patterns in the distribution of taxa and highlighted key differences in community composition among the sampling sites and years.

Initial trimming and filtering of nontarget taxa (see Methods) for the 18S marker resulted in a final read count of 235,064,859 (mean 1,237,184 per sample) prior to normalizing read counts and applying read count thresholds to samples. After normalizing samples, the lowest taxonomic rank that could be assigned to each OTU was as follows: 11% to class, 15% to order, 30% to family, 25% to genus, and 19% to species. The rbcL marker analysis resulted in a final read count of 221,898,210. The lowest taxonomic rank that could be assigned to each OTU was as follows: 3% to class, 3% to order, 21% to family, 27% to genus, and 46% to species. The read counts were grouped by station and means calculated for each OTU. A list of OTU and read count assignments are available at the Arctic Data Center (Koch et al., 2025).

To allow for the analysis of all available taxonomic ranks across the 18S gene data, read counts were then assigned to primary diet items by eukaryotic groups ( Figure 2 ). These diet categories included cercozoans, prasinophytes, bacillariophytes (henceforth, referred to as diatoms), dinophytes (henceforth referred to as dinoflagellates), silicoflagellates, and other, which included raphidophytes, prymnesiophytes and ciliates. The relative 18S read count of each eukaryotic group by year and location indicated a dominance of diatoms in both years across nearly all stations. In 2018, SLIP1 (85%), DBO4.3 (89%) and DBO4.5 (86%) had the highest relative read count of diatoms. In 2019, UTN3 had ~92% relative read count of diatoms. In 2018, the proportional 18S read abundances of dinoflagellates in the clams’ diet at UTN3 and UTN4 increased by 43% and 41% respectively; these stations are located immediately north of the Bering Strait ( Figure 1 ). In 2019, relative dinoflagellate 18S read counts increased moving north of the Bering Strait into the Chukchi Sea. In the southeast Chukchi Sea (SEC3), proportional read counts of dinoflagellates in the clam diets exceeded those for diatoms (67% versus 25%; Figure 2 ).

The proportions of prasinophytes, cercozoans, and silicoflagellates decreased within the 18S read counts from 2018 to 2019 ( Figure 2 ). Of these groups, prasinophytes comprised a substantial proportion of all normalized 18S reads, which included the species Micromonas polaris and Pyramimonas parkea ( Figures 3A, B ). In 2018, the station UTN2 had the highest relative proportion of prasinophytes at 27% ( Figure 2A ). In 2019, this group was highest at SLIP1 and SLIP2 (~19% proportion of normalized read counts at both stations; Figure 2B ).

Dinoflagellates comprised a substantial proportion of the normalized 18S mean read counts in both years. In 2018, dinoflagellates ranged from 6% at DBO4.3 to 43% at UTN3 ( Figure 2 ). In 2019, the range of proportional normalized read counts ranged from 3% at UTN3 to 67% at SEC3. We also detected frequent normalized read counts of the dinoflagellate Alexandrium in the 18S results throughout both years and stations ( Figure 4 ). In 2018, the highest normalized read counts for Alexandrium occurred at station UTN2, which is north of Bering Strait (6803 counts) and these dinoflagellates were not detected in clams collected from the St. Lawrence Island polynya region, in one station immediately north of the Bering Strait (UTN3), or in the southeast Chukchi Sea (SEC stations). Normalized read counts ranged from 1046 to 3486 at the other UTN stations in 2018. In 2019, normalized read counts of Alexandrium were identified across all stations (noting that DBO4 was not sampled in 2019). Normalized read counts were greatest at UTN3 (3139 counts) and UTN4 (3849 counts). Detection of Alexandrium reads in the southeastern Chukchi Sea notably increased in 2019 compared to 2018.

To address our initial objective of identifying ice-associated and pelagic diatoms, and to explore the full extent of our dataset, we categorized the 18S and rbcL read counts for diatoms by family including: Chaetocerotaceae, Thalassiosiraceae, Melosiraceae, Paraliaceae, and two general categories, Other (centrics), and Other (pennates; Figure 5 ). These families were largely represented by the following genera: Chaetoceros, Thalassiosira, Melosira, and Paralia. The ‘other–centric’ diatom category included Biddulphia, Leptocylindrus, Minutocellus, Rhizosolenia, and Stephanopyxis. The ‘other–pennate’ diatom category included the genera Fragilariopsis and Navicula. Approximately 89% of the OTUs within the Chaetocerotaceae family were assigned to species level, and Chaetoceros is the probable genus represented. Similarly, ~80% of the Thalassiosiraceae OTU assignments were to genus or species-level for Thalassiosira and are the probable genus represented.

In 2018 and 2019 across both 18S and rbcL markers, Chaetoceros sp. dominated clam diets, with dietary proportions greater than 75% in all stations in 2018 and greater than 80% at all stations in 2019, with the exception of SEC3 ( Figure 3 ). Normalized read counts at SEC3 were composed of ~30% other (centrics), which could only be assigned to the order Thalassiosirales ( Figure 5 ). Otherwise, the genera Thalassiosira was the second-most frequently occurring in read count abundance at most other stations. For the normalized read counts assigned to species, the 18S results indicated a dominance of Chaetoceros gelidus across all stations, while Chaetoceros radicans was the most abundant species assignment for the rbcL marker ( Figure 3 ). Melosira sp., a centric, ice-associated species, increased in relative frequency (>1%) at the northernmost stations, DBO4.3 and DBO4.5 in 2018 ( Figures 5A, C ). The pennate ice-associated diatom, Nitszchia frigida, was not detected and nor was the pelagic diatom, Cylindrotheca closterium, major components of the diatom flux in 2018 and 2019 (Lalande et al., 2021; Fukai et al., 2020).

In this study, we explored the primary food sources for the clam, M. calcarea, in the Pacific Arctic during extreme low sea ice years, using DNA metabarcoding techniques to provide increased understanding of bivalve diets at a higher taxonomic resolution than possible through gut content and microscopy counts (Reid and Reid, 1969; North et al., 2014). Our findings are consistent with a diet consisting primarily of pelagic diatoms, specifically Chaetoceros (C. gelidus or C. radicans, Figure 3 ) as the dominant food source for these clams in July of 2018 and 2019. C. gelidus occasionally dominates the Arctic spring bloom (von Quillfeldt, 2000; Sergeeva et al., 2010). In the summer of 2018, C. gelidus was reported to be a major constituent of the diatom community in the Alaska Coastal Water (ACW) mass south of St. Lawrence Island (Fukai et al., 2020), which corroborates the high proportional reads of this species in our samples. While these two markers suggest different species, it is more likely that these markers represent the same species and exemplify how these two gene markers can produce differing results (Dermastia et al., 2023).

Linkages between the pelagic and benthic components of the ecosystem can further be interpreted using diatom fluxes from sediment traps located within the DBO regions 2 (closest to UTN2) and 3 (nearest SEC2), as reported by Lalande et al. (2021). The flux data from these sediment traps are not fully consistent with our data reported here but generally align with the proportion of normalized read counts for Chaetoceros and Thalassiosira that we observed. For example, in Lalande et al. (2021) the particle flux in 2018 following the spring bloom at DBO2 had a large proportion of Chaetoceros diatoms followed by Thalassiosira. In 2019, the spring bloom had a large proportion of Thalassiosira and Fragilariopsis. At DBO 3, large fluxes of Thalassiosira and Fragilariopsis also occurred from April through June in 2018 and 2019, but less so for Chaetoceros sp. Notably, Fragilariopsis, a pennate diatom, appeared minimally in the clam diets in both years, perhaps indicative of a preference for centric diatoms (Ward et al., 2004). A higher proportion of centric diatoms had also been observed for Macoma sp. in the northern Bering Sea based on gut content and microscopy analysis (North et al., 2014). The lack of Cylindrotheca closterium in our samples was notable, given the large fluxes of this pennate diatom in nearby sediment traps at DBO 2 and 3 (Lalande et al., 2021). Comparisons of 18S, rbcL, and microscopy methods in other Arctic studies also reported a lack of expected reads for C. closterium based on microscopy results (Dermastia et al., 2023). They concluded that these diatoms were previously misidentified by microscopy and were probably Pseudo-nitzschia sp., which look similar and were well-represented in their read assignments (Dermastia et al., 2023). Despite these specific inconsistencies, linkages between the diet of the clams and the sediment trap data underscores the importance of pelagic, centric diatoms as a significant resource for the benthic food web in this region and demonstrates the role of benthic-pelagic coupling.

Diet studies using biogeochemical approaches such as fatty acid biomarkers, stable isotope analysis, and (or) direct examination of gut contents of bivalves in the Pacific Arctic region previously indicated dietary plasticity and a strong reliance on pelagic phytoplankton (e.g., Budge et al., 2007; Iken et al., 2010; Oxtoby et al., 2016; Kędra et al., 2019; North et al., 2014). DNA metabarcoding is a promising addition to complement these methodologies, because it potentially facilitates study of a greater diversity and number of species that are not practical for other analytical techniques. Shifts in biodiversity are poorly understood in the Arctic (Feng et al., 2021), which underscores the complexity of Arctic marine food webs and the importance of adopting a multifaceted approach to understand the ecological dynamics of the changing Arctic environment. DNA metabarcoding therefore offers a comprehensive and nuanced perspective on the complex interactions within marine food webs, enhancing our understanding of ecosystem dynamics and biodiversity. This approach not only enables tracing of the intricate trophic links that sustain marine life but also provides the knowledge to predict and respond to the ecological shifts anticipated in a rapidly changing Arctic (Geraldi et al., 2024).

DNA metabarcoding is not without its limitations (Coissac et al., 2012). One significant constraint for our study was the differentiation between DNA derived from ingested food and DNA from environmental contamination or from the digestion of prey, which can complicate interpretations of diet (Garrison et al., 2022). More specifically, samples could potentially be contaminated with phytoplankton assemblages not actually consumed, but rather introduced into the interior clam tissues when the siphons retract into the shells and/or during digging with its foot (Reid and Reid, 1969). Bacterial breakdown of large diatoms may be necessary for incorporation into bivalve diets in this region (North et al., 2014). Another limitation in characterizing diet is the incomplete representation of potential prey items in reference databases, which can lead to the underestimation or misidentification of components within a diet (Dermastia et al., 2023). This gap in reference data particularly affects rare or less-studied species, whose DNA sequences might not yet be cataloged. Furthermore, DNA metabarcoding is sensitive to the relative abundance of DNA, meaning that it may overemphasize the presence of species that shed more DNA or whose DNA degrades more slowly, potentially skewing the perceived dietary composition (Murray et al., 2011; Garrison et al., 2023).

The winter of 2018 was remarkable for the lack of sea ice formation in the northern Bering Sea, which corresponded to more typical winter conditions of the southeastern Bering Sea (Stabeno and Bell, 2019; Thoman et al., 2020). As a result of a marine heat wave, sea ice was thinner and formed later in the Chukchi Sea in 2017, leading to the late arrival of sea ice in the Bering Sea in November (Stabeno and Bell, 2019). Warm, southerly winds in February blew the sea ice present in the Bering Strait region out by mid-February 2018 (Kikuchi et al., 2020). There was no sea ice present south of St. Lawrence Island in the spring of 2018, an event not identified in observed or reconstructed conditions since 1850 (Walsh et al., 2017). As a result, no ice-associated spring bloom occurred, but a phytoplankton bloom occurred ~1 month later than usual and only once the surface waters thermally stratified, in contrast to the salinity stratified surface waters typically resulting from ice melt (Duffy-Anderson et al., 2019; Kikcuhi et al., 2020). Accordingly, the summer phytoplankton community in the northern Bering Sea had transitioned from cold-water species in 2017 to more cosmopolitan species in 2018, including Chaetoceros sp (Fukai et al., 2020). The marine heat wave continued into 2019, the second lowest winter sea ice in the northern Bering Sea on record (Stabeno and Bell, 2019). In these years, the sea ice retreat in mid-winter provided unfavorable conditions for a spring ice-associated bloom due to lack of sunlight and atypical winter conditions, leading to another unusual spring phytoplankton bloom (Kikuchi et al., 2020). Therefore, the clams in our study (with the exception of DBO4.3 and DBO4.5 in 2018) were sampled approximately 5 months after sea ice retreat, from nutrient depleted water in the Bering Sea, with more cosmopolitan species of diatoms (Fukai et al., 2020).

The appearance of Melosira sp. from the DBO 4 region in 2018 suggests a nuanced aspect of food availability relative to sea ice persistence and timing. While we were not able to assign this OTU to species, Melosira arctica is the most probable species to be present based on taxonomic studies (Poulin et al., 2011, 2014). M. arctica are a centric, ice-associated colonial diatom typically associated with sub-ice habitat that are released as the ice retreats in July in high latitude locations (Booth and Horner, 1997; Lalande et al., 2019). In 2018, the sea ice had retreated from the Chukchi Sea (DBO4 region) in mid-June, approximately 3 weeks prior to sampling as opposed to the 5+ months in the northern Bering Sea (Perovich et al., 2018). Therefore, the presence of Melosira suggests that there was a recent ice algae bloom in the Chukchi Sea. While not represented in the 18S or rbcL marker reads, the presence of ice-associated species such as Nitzschia frigida have been observed in greater abundances within particle fluxes in the Chukchi Sea, highlighting the spatial variability in food sources (Koch et al., 2020b; Lalande et al., 2021). It is probable that ice algae and ice-associated blooms are an important component in the diets of Macoma clams under specific sea ice conditions and at certain times of the year.

Our findings were also consistent with lipid biomarkers analyzed from Macoma clams collected from the same grab samples as the 2018 sampling effort reported here (Koch et al., 2020a). The biomarker results also indicated a minimal or absent ice algal signature in the SLIP region and in the Bering Strait/southern Chukchi Sea, which was attributed in that study to clams suspension feeding rather than surface deposit feeding in the summer. Sub-surface deposit feeders in the northern Bering Sea have access to benthic reserves of ice-associated organic matter (North et al., 2014; Pirtle-Levy et al., 2009). With the early pulse of food from ice-associated blooms in the spring, benthic activity is stimulated, and bioturbation increases, leading to rapid utilization and burial of ice-associated organic matter (Grebmeier et al., 2006; Clough et al., 1997; McMahon et al., 2006). With this species capable of engaging in both suspension and surface deposit feeding, those results were likely reflective of the prevailing conditions in the water column following the most recent bloom, whereas subsurface deposit feeding bivalves in the sediments can access ice-associated organic carbon buried in the sediments long after the ice algal bloom has concluded (Koch et al., 2020a). Given the timing of sampling in 2018 relative to sea ice retreat and the overall lack of ice-associated production, ice-associated diatoms are unlikely to be a significant component of the summer diet of Macoma clams in the northern Bering Sea.

The contributions of the prasinophytes, Micromonas and Pyramimonas, to the diet of these clams were also potentially indicators of nutrient availability and bloom progression. Micromonas polaris had abundant read counts in both 2018 and 2019 throughout the region and Pyramimonas parkaeae had elevated normalized read counts in the St. Lawrence Island polynya region in 2019 ( Figure 3 ). These species are indicative of shifting phytoplankton assemblage structures in a changing Arctic, where prasinophytes thrive under very low nitrogen (N) conditions while diatoms require high N and Si (Ardyna and Arrigo, 2020). Accordingly, in 2018, the spring bloom was delayed by ~1 month (Stabeno and Bell, 2019) and nutrients were depleted in the upper mixed layer resulting in non-diatom phytoplankton blooms (Fukai et al., 2020).

The temporal mismatch between our DNA metabarcoding of benthic consumers and the window of ice-associated production highlights the potential for additional work. Conducting metabarcoding studies throughout various seasonal transitions—from ice breakup through to resuspension and the autumn blooms—could provide a more comprehensive understanding of the variation in clam diets and resource utilization across different bloom phases. A similar approach examining the diet of the arctic copepod, Calanus glacialis, revealed similar complexities where higher proportions of ice-associated diatoms were observed at certain times of the year but also found higher abundances of the pelagic diatom, Thalassiosira, when these organisms were observed foraging directly under the sea ice (Cleary et al., 2017). There are gaps remaining in our understanding of ice-associated diatom community composition, species that span these niches (sympagic diatoms seeding phytoplankton blooms, and vice versa), and increasingly the role of under ice algal blooms (Ardyna et al., 2020).

Over the last few years, the Pacific Arctic has had increasing attention focused on harmful algal blooms (HABs), which may pose an escalating threat to coastal communities and marine ecosystems (Anderson et al., 2022; Lefebvre et al., 2022). Alexandrium catenella stands out for its toxicity and potential to disrupt marine food webs, particularly within bivalves, seabirds, and marine mammals, as well as impact human health (Lefebvre et al., 2022). Our study revealed the presence of Alexandrium sp. across all stations sampled during our study period, with the highest detection observed in 2019 at station UTN2. Like other organisms of interest in this study, we were only able to identify this OTU to genus, but A. catenella is widely distributed and abundant in this region as confirmed by other studies (Anderson et al., 2022; Brosnahan et al., 2020; Lefebvre et al., 2022). Additionally, over half of the 30 known species of Alexandrium produce paralytic shellfish toxins (Anderson et al., 2012). The widespread detection of Alexandrium in clams aligns with recent findings from Anderson et al. (2022) and Lefebvre et al. (2022), who reported large dinoflagellate cyst beds in the sediments at the same stations we sampled. Furthermore, these studies noted the presence of saxitoxins—a potent neurotoxin produced by Alexandrium—in clams and walruses in 2019 within the Bering Strait and southeast Chukchi Sea regions, areas that correspond to where we detected signatures of Alexandrium.

In addition to Alexandrium, we also detected the presence of other potentially harmful algal species, such as the raphidophyte Chattonella subsalsa (Hallegraeff and Hara, 2003), albeit less frequently than A. catenella. While these taxonomic assignments were rarer, they nonetheless signify the potential for early detection and monitoring of future HAB events in Alaska waters. The presence of these species, along with A. catenella, highlights the potential threat that HABs present to marine ecosystems and human health in the Arctic (Anderson et al., 2022; Pućko et al., 2023). Our findings contribute to the growing body of evidence indicating that monitoring and research into emerging HAB dynamics in a changing Pacific Arctic should be prioritized alongside efforts to manage their impacts on marine ecosystems and the health of coastal communities in Alaska (Anderson et al., 2022). The ability to use environmental DNA (eDNA) metabarcoding to rapidly identify harmful algal species has profound implications for public health, particularly in the context of monitoring and managing risks associated with paralytic shellfish toxins (Jacobs-Palmer et al., 2021). The detection of harmful species through this method not only highlights the direct links between marine biodiversity and human health but also the potential of eDNA metabarcoding as an accessible tool for supporting community-based monitoring programs (Larson et al., 2020).