Bowhead whale faeces link increasing algal toxins in the Arctic to ocean warming

Sample and data collection

To define the relationship between warming climate conditions and algal toxin prevalence and concentration in Arctic food webs, we used the following data: algal toxin concentrations in bowhead whale faeces, A. catenella cell and cyst densities, presence of Pseudo-nitzschia frustules in bowhead whale faeces, ocean heat flux, wind velocity, SLP, open water area, SSTs and annual minimum sea ice extent.

Bowhead whale faecal sample collection

During 2004–2022, faecal samples were collected from 205 bowhead whales harvested for subsistence purposes and landed at Utqiaġvik (formerly Barrow), Alaska, during autumn harvest seasons (August–October; Supplementary Table 1). The whales are known to feed near Utqiaġvik at depths ranging from shallow continental shelf (45 m) to the deeper waters (more than 300 m) of Barrow Canyon¹⁸. Depths at which Alexandrium cells mainly occur are in the upper 25 m (ref. ⁴) and Pseudo-nitzschia particulate DA measurements have been do***ented from the surface (about 2 m) to chl-a maximum depths (20–40 m)²⁷. Whales are typically harvested within a 30-mile radius of Utqiaġvik³⁷. Sections of colon from each whale were dissected and faecal matter was removed using plastic spoons. Samples were stored frozen in Whirl-Pak bags at −20 °C until blockysed for algal toxins.

Quantification of algal toxins in faecal samples

Algal toxins were extracted from frozen bowhead whale faecal samples (n = 205; Supplementary Table 1) by slowly thawing and stirring cold faecal material followed by subsampling into approximately 1 g for blockysis. To each aliquot, 50% methanol was added at a volume of three times the aliquot weight for a one-in-four dilution. Samples were vortexed briefly on high (Analogue Vortex Mixer, sn 060223013, VWR) and homogenized using a generator probe (GLH 850, Omni-International) for 1 min at 2,100 rpm. Homogenates were then centrifuged at 4,100 rpm for 20 min at 4 °C (CR3i centrifuge, Jouan) and supernatants were poured off and stored in 4-ml amber glblock vials in a dark refrigerator (about 1 °C). Directly before blockysis, 200 μl of each sample were filtered through 0.22-μm Ultra-Free Centrifugal filters (UFC30GVNB, Millipore Sigma) in a tabletop centrifuge (AccuSpin Micro 17, Fisher Scientific) at 12,000 rpm. Methanol (50%) is the standard extraction solvent for DA enzyme-linked immunosorbent blockay (ELISA) blockyses³⁸ and has also been validated as an effective extraction solvent for STX ELISA blockyses³⁹; hence, these 50% methanol extracts were used for both DA and STX quantifications.

Algal toxins were quantified in bowhead whale faecal samples by means of commercially available direct-competition ELISA kits. DA was quantified using Biosense ASP ELISA kits (PN A31300401, Biosense Laboratories) for samples collected in 2004–2021 and using the comparable ABRAXIS Domoic Acid ASP ELISA kits (PN 520505, Gold Standard Diagnostics) for samples collected in 2022. STX was quantified using ABRAXIS Saxitoxins PSP ELISA (PN 52255B, Gold Standard Diagnostics) for all samples. Although these kits are designed to blockyse shellfish and water samples, previous studies have determined appropriate dilutions to avoid matrix effects from marine mammal matrices and have validated ELISA results compared to other blockytical methods^12,40,41. Kits were used according to manufacturer’s instructions with dilution modifications from ref. ⁴¹ and ref. ¹² for DA and STX, respectively. Sample extracts were diluted 1:100 for DA and 1:50 for STX (sample to kit-provided sample diluent). Standards, controls, blanks, samples and kit-provided reagents were then added to ELISA plates in duplicate and processed following the kits instructions. Toxin quantifications were obtained using a BioTek Epoch plate reader (sn 257814) and concentrations (ng of toxin per g of sample) were determined using a four-parameter logistic curve model based on the known standards concentrations. Samples with concentrations exceeding the detection range of the kit (defined as 20–80% of the standards range) were diluted further and reblockysed until concentrations fell within the detection range. Minimum blockay detection limits were 4 ng g⁻¹ for DA and 4.7 ng g⁻¹ for STX. All faecal samples collected from 2010 to 2022 were blockysed within the year they were collected. Faecal samples from 2004 to 2009 were blockysed within 5 years of collection. To rule out any potential toxin degradation issues, studies were performed with bowhead faecal material stored over 4 years under various storage conditions. Results from both STX and DA storage studies confirmed that long-term frozen storage did not impact toxin concentrations over time^39,42. In several years of bloom sampling across the region during this project (2019, 2022 and 2023), the suite of toxins produced by Pacific Arctic A. catenella has been consistently dominated by gonyautoxin-4, neosaxitoxin, gonyautoxin-3 and STX¹³. Unfortunately, both gonyautoxin-4 and neosaxitoxin have low cross-reactivities with the ELISA test (less than 2%), but STX is picked up at 100% and gonyautoxin-3 at 23%. So, although the ELISA is probably underestimating the total amount of toxin in these faecal samples, the overall consistency observed in toxin profiles of regional A. catenella strains across years indicates that ELISA data are appropriate for blockessing relative temporal trends in toxicity, such as the results reported in this study. These STX quantifications are representative of prevalence and relative concentrations equally over the two decades of sampling.

A. catenella cyst sample collections

To map A. catenella cyst abundance in the study region, surface sediment samples were collected during 12 different cruises over 5 years (2018–2022; Supplementary Table 4). Sediments were collected using a Van Veen or Smith-McIntyre grab, and a cut syringe was used to collect a plug from the 0–3-cm layer of each grab; in some cases, several plugs were collected and pooled together. Each subsample was homogenized, sealed in an airtight container and maintained in the dark at 0–4 °C.

Cyst microscopy and mapping

All sediment samples were processed using a primulin stain (Extended Data Fig. 1), allowing A. catenella cysts to be enumerated following established methods⁴³. Briefly, an aliquot of each homogenized sediment sample was diluted (1:5) in filtered seawater and sonicated (Brandon Sonifier 250, 40% amplitude, 60 s) in an ice bath. The resulting slurry was sieved to isolate the 20–80-μm size fraction, which was then resuspended in filtered seawater and preserved with 5% formalin. The formalin-preserved samples were chilled (4 °C) for 1–3 h, after which they were centrifuged (3,000g, 10 min), formalin–seawater supernatant was aspirated and sediment pellet was resuspended in chilled methanol. After refrigeration (4 °C) in methanol for at least 72 h, the samples went through a series of centrifugation and aspiration steps to transfer from methanol to deionized water and from deionized water to primulin stain (2 mg ml⁻¹, 2 ml per sample). Samples were rotated on a Labquake for 1 h at 4 °C, after which more centrifugation steps were used to wash the sample in deionized water and resuspend each sample in a final volume of 10–15 ml of deionized water.

A. catenella cysts were enumerated using a Zeiss Axioscope epifluorescence microscope equipped with a FITC filter set (Zeiss 09, excitation 450–490 nm band pblock; emission 515 nm long pblock) under a ×10 objective. Samples that were too dense to count were diluted 1:10; all raw counts were normalized to A. catenella cysts per cubic centimetre. Results from all cruises were compiled together and a map of A. catenella cyst abundance (Fig. 1d) was produced in Matlab (R2024a) using the m_map package and an interpolation method which weights along isobaths⁴⁴ and which has been previously used to map cyst distributions in the region⁴.

A. catenella vegetative cell sample collections

During the HLY1901 cruise (Supplementary Table 4), vegetative cell concentrations of A. catenella were quantified in discrete water samples collected and preserved at process stations. These samples were used to detect and characterize a bloom of A. catenella on the Barrow Canyon East transect line on 20 August 2019 (Supplementary Table 2). At each station, conductivity–temperature–depth (CTD) data were collected using a Sea-Bird 911 plus CTD mounted to a 24-position rosette⁴⁵. Discrete 2-l water samples were collected from Niskin bottles representing surface (about 3 m) and 10-m, chlorophyll maximum and bottom depths. These samples were immediately sieved through a 15-μm Nitex mesh; all captured particles were backwashed with filtered seawater into a 15-ml conical tube and fixed with formalin (5% final concentration). Water samples were stored at 1 °C for up to 72 h, at which point they were centrifuged (3,000g, 10 min) and overlying seawater–formalin mixture was aspirated. The phytoplankton pellet was resuspended in chilled methanol and all samples were stored at −20 °C.

Vegetative cell microscopy and enumeration

A. catenella were enumerated in preserved seawater samples using a fluorescence in situ hybridization method following previously published procedures⁴⁶. Briefly, an aliquot of each methanol-preserved sample was transferred to a filtration manifold column fit with a 5.0-μm pore size, 25-mm diameter Cyclopore membrane filter. Vacuum suction was used to remove the methanol from each manifold column and replace it with 1 ml of prehybridization buffer (5× SET (750 mM NaCl, 5 mM EDTA, 100 mM Tris-HCl, pH 7.8), 0.1 µg ml⁻¹ of polyadenylic acid, 0.1% IGEPAL CA-630, 10% formamide). After a 5-min room-temperature incubation period, the prehybridization buffer was removed and replaced with 1 ml of hybridization buffer (prehybridization buffer augmented with 4.8 µg ml⁻¹ of Cy3 NA-1 probe). The NA-1 oligonucleotide probe (5′ Cy3-AGT GCA ACA CTC CCA CCA-3′) was selected to label A. catenella large subunit ribosomal RNA. The samples were then incubated in the dark (50 °C, 1 h), after which the hybridization buffer was removed and replaced with 1 ml of wash buffer (0.2× SET) for an extra 5-min room-temperature incubation. All remaining buffer was removed by means of vacuum filtration and filters were mounted on slides with a small volume (20–40 µl) 80% glycerol 25× SET solution. Samples were stored in the dark at 4 °C for up to 3 days before enumeration. In each sample, all A. catenella vegetative cells were enumerated at ×10 on a Zeiss Axioscope M1 using a Cy3 filter set (Chroma no. 49016/TRITC long pblock); these cell counts were then normalized to cells per litre to determine in situ concentrations.

Pseudo-nitzschia frustules in faecal samples

Subsamples of bowhead whale faecal samples (about 0.05–0.2 g) were prepared for SEM using published methods⁴⁷. Briefly, faecal samples were rinsed three times with 1 ml of distilled water in 1.5-ml microcentrifuge tubes. Each rinse step included vortexing and centrifugation of the pellet. Pellets were then oxidized for 2 h with four or five drops of saturated potblockium permanganate solution, cleared with three rinses of concentrated HCl, rinsed again three times with distilled water and then filtered onto 13-mm diameter, 1.2-µm pore size, polycarbonate filters (Millipore Corp). Filters were then glued to aluminium stubs, coated with gold-palladium and viewed in a JEOL 6360LV SEM. Species determinations were made using published morphological characteristics⁴⁸ (Extended Data Fig. 1).

Mooring data and heat flux calculation

We used the hydrographic and velocity data from a mooring close to Barrow Canyon, maintained since 2002, as part of Arctic Observation Network^49,50. To compute the lateral heat flux in the upper layer of the Beaufort shelfbreak jet, we used the temperature and velocity data measured at the top float of the mooring (approximately 35-m depth). The ac***ulated lateral heat flux (H) is computed as:

where ρ is potential density, C_ρ is the specific heat of seawater, θ(t) is the time-dependent potential temperature, θ_r = −1.91 °C is the reference temperature and u(t) is the time-dependent alongstream component (125° clockwise from north) of the velocity. The heat flux was averaged for each composite toxin concentration group of bowhead whales (low DA n = 90, STX n = 80; moderate DA n = 73, STX n = 63; high DA n = 42, STX n = 62).

Wind velocity and SLP reblockysis

We used the hourly wind field and SLP from the ERA5 reblockysis, with a spatial resolution of 0.25°, provided by the European Centre for Medium-Range Weather Forecasts⁵¹. The ERA5 data have shown good agreement with observations in the western Arctic Ocean⁵².

Open water anomalies and SST baseline departure

Environmental data from the Beaufort Sea including SST (°C) and open water area (km²) data were retrieved from the NOAA OI SST v.2 (ref. ⁵³) and the National Snow and Ice Data Center (NSIDC)³⁶ databases, respectively, for years when whales were harvested for subsistence purposes (2004–2022) and to compare with environmental baselines (1982–2011) (SST data, Supplementary Table 5; open water data, Supplementary Table 6).

July SSTs departures from baseline (that is, z-score correction) in the Beaufort Sea were calculated using the equation: July SST (°C) − mean July baseline (1982–2011) SST (°C))/standard deviation July baseline SST (°C) (Supplementary Table 5). Anomalies of open water in the Beaufort Sea were calculated by dividing summer monthly (June, July, August and September) open water areas (km²) by their respective monthly mean baseline (1982–2011) values (Supplementary Table 6).

Years were categorized into two DA prevalence groups; years with 100% of whales testing positive (more than 0 ng g⁻¹) for DA (n = 7 (years) with n = 56 total whales, ‘100% DA prevalence’) and years with less than 100% DA prevalence (n = 12 (years) with n = 149 total whales, ‘less than 100% DA prevalence’) and compared to open water anomalies in the Beaufort Sea during summer months (June, July, August and September). Linear models were constructed to test whether the open water anomalies of each summer month (June–September) were significantly different during years of 100% DA prevalence and less than 100% DA prevalence in whales while weighting each model by the number of whales tested for DA per year. Weighted estimated marginal means of monthly open water anomalies for each DA prevalence group were tested using unpaired t-tests (Fig. 4a). All months (June, July, August and September) had significantly higher open water anomalies in the Beaufort Sea in years when DA prevalence was 100% compared to years when DA was present in less than 100% of whales harvested from the Beaufort Sea (Fig. 4a and Supplementary Table 3). Weighted estimated marginal mean comparisons were repeated for two more prevalence categories: (1) more than 90% versus less than 90% prevalence and (2) more than 75% versus less than 75% (Extended Data Fig. 2 and Supplementary Table 7). Although similar statistical relationships exist for the 90% prevalence comparison for all months, except September, they are not as significant as for the 100% prevalence comparisons (Fig. 4a, Extended Data Fig. 2a and Supplementary Tables 3 and  7a). For the greater than 75% prevalence comparison, the relationship remains for June, July and August as for 100% prevalence comparisons, but anomalies are similar among groups in September (Extended Data Fig. 2b and Supplementary Table 7b). For STX, all months (June, July, August and September) had significantly higher open water anomalies in the Beaufort Sea in years when STX prevalence was 100% compared to years when STX was present in less than 100% of whales harvested from the Beaufort Sea (Extended Data Fig. 3 and Supplementary Table 8).

Pearson correlations were performed among June open water anomalies and July SST anomalies for years blockociated with bowhead toxin blockyses (n = 19 years) (Fig. 4b). Analysis was done using software programs R⁵⁴ and R studio⁵⁵ and R packages emmeans⁵⁶, lme4 (ref. ⁵⁷) and ggpubr⁵⁸.

Average SSTs (1900–2023) for the different Seas

Average SST (°C) during May–September for years 1900–2023 in the Bering, Chukchi and Beaufort Seas (Fig. 5a–c) were obtained from the NOAA Extended Reconstructed SST v.5 data provided by the NOAA PSL, Boulder, Colorado, USA, from their website at (ref. ³⁵).

Annual minimum sea ice (1979–2024) for the Seas

Annual minimum sea ice extent (km²) data for the Bering, Chukchi and Beaufort Seas during 1979–2024 were acquired from the NSIDC³⁶. Daily sea ice extent for the Bering, Chukchi and Beaufort Seas were summed and the minimum daily extent for each year was plotted in Fig. 5d.

Inclusion and ethics statement

This study was a mutually beneficial collaboration between NOAA, NWFSC, WARRN-West and the whaling communities of the NSB. The stakeholders’ needs and concerns were the top priority for all aspects of this work. The project stemmed from a 15-year collaboration in service to the NSB Department of Wildlife Management (DWM), the Alaska Eskimo Whaling Commission and Whaling Captains’ Associations, for health blockessments of harvested bowhead whales.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.