Mass mortality events are increasing in frequency and magnitude, potentially linked with ongoing climate change. In October 2016 through January 2017, St. Paul Island, Bering Sea, Alaska, experienced a mortality event of alcids (family: Alcidae), with over 350 carcasses recovered. Almost three-quarters of the carcasses were unscavenged, a rate much higher than in baseline surveys (17%), suggesting ongoing deposition and elevated mortality around St Paul over a 2–3 month period. Based on the observation that carcasses were not observed on the neighboring island of St. George, we bounded the at-sea distribution of moribund birds, and estimated all species mortality at 3,150 to 8,800 birds. The event was particularly anomalous given the late fall/winter timing when low numbers of beached birds are typical. In addition, the predominance of Tufted puffins (Fratercula cirrhata, 79% of carcass finds) and Crested auklets (Aethia cristatella, 11% of carcass finds) was unusual, as these species are nearly absent from long-term baseline surveys. Collected specimens were severely emaciated, suggesting starvation as the ultimate cause of mortality. The majority (95%, N = 245) of Tufted puffins were adults regrowing flight feathers, indicating a potential contribution of molt stress. Immediately prior to this event, shifts in zooplankton community composition and in forage fish distribution and energy density were documented in the eastern Bering Sea following a period of elevated sea surface temperatures, evidence cumulatively suggestive of a bottom-up shift in seabird prey availability. We posit that shifts in prey composition and/or distribution, combined with the onset of molt, resulted in this mortality event.

Funding: These analyses were supported by the National Science Foundation (NSF) EHR/DRL award 1322820 and Washington Department of Fish and Wildlife award 13-1435, to JKP. Support for algal toxin analyses was provided by NSF OCE-1314088 to KAL and Dr. David Marcinek at the University of Washington. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability: All relevant data are within the manuscript and its Supporting Information files, with the exception of wind data that is openly available from the Earth System Research Laboratory (ESRL) website. Wind data used in this study came from the NCEP North American Regional Reanalysis (NARR) Pressure level dataset available from ESRL (URL: https://www.esrl.noaa.gov/psd/data/gridded/data.narr.pressure.html ), listed as 3-hourly zonal (uwnd) and meridional (vwnd) wind speeds.

This is an open access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.

In this paper, we document a MME of marine birds, primarily Tufted puffins (Fratercula cirrhata), on St. Paul Island, Pribilof Islands, Alaska, in the eastern Bering Sea ( Fig 1 ) during the late fall/winter of 2016/2017. The Pribilof Islands, located near the edge of the Bering Sea continental shelf, support one of the largest concentrations of breeding seabirds (>2 million) in the North Pacific [ 42 , 43 ]. The islands have also been hunting and harvesting grounds to Unangan (or Aleut) for millennia, with permanent settlements on both islands established in the late 1700s. Several species of seabirds are important cultural and subsistence resources [ 44 ], and as such seabird mortality events are both an ecological and societal concern for island residents. Using a combination of long-term standardized beached bird surveys and intensive surveys during the event, we characterize this MME in terms of timing, abundance, species composition and carcass condition as compared to baseline measures. We use wind forcing and carcass count phenology to model daily deposition and provide estimates of total mortality. Our results add to the growing body of literature documenting marine bird MMEs in the northeast Pacific associated with recent and persistent warming [ 10 , 38 , 45 ].

As abundant, visible, upper-trophic organisms, seabirds have been proposed as indicators of marine ecosystem shifts due to climate, with documented effects of climate variability on both reproduction [ 28 – 30 ] and adult survival [ 31 – 33 ]. Large-scale shifts in climate have been punctuated by large mortality events of marine birds [ 34 – 38 ]. These “massive mortality events” (MME)—defined as catastrophic, but often short-lived, periods of elevated mortality—can affect substantial proportions of a population, occasionally with long-term consequences to population size [ 39 ]. Seabird MMEs are perhaps one of the most frequently occurring and widely reported types of MME in the literature [ 40 ], potentially due to their perceived and absolute (mortality often exceeding 10,000s-100,000s birds; [ 35 , 38 , 40 , 41 ]) magnitude.

The Bering Sea is a high latitude, semi-enclosed sea between the north Pacific and Arctic Oceans [ 12 ], notable for having an extensive continental shelf and seasonal ice-cover that varies in extent on interannual to multi-decadal time scales [ 2 , 13 ]. Ecosystem structure, including the timing and composition of primary production and primary/secondary consumers, varies markedly among early and late ice-retreat years [ 14 – 19 ]. The eastern Bering Sea supports some of the most economically important fisheries in the world [ 1 , 20 ], hosts large populations of marine mammals [ 21 , 22 ], and is the breeding and/or summering ground for ~30–40 million marine birds [ 23 – 25 ]. Bering Sea food webs are particularly sensitive to bottom-up climate effects, as changes in atmospheric forcing impacts sea ice, as well as the extent of the ‘cold pool’, a lens of cold (< 2°C) near-bottom seawater that acts as a refuge for cold-water associated species, and also promotes primary production through summer/fall [ 18 , 25 , 26 ]. Over the last two decades, several multi-year stanzas of warm (2000–2005, 2014-present) and cold (2007–2013) conditions have been observed in the southern Bering Sea, which have been linked to variability in phytoplankton biomass (lower in warm years), copepod species composition (reduced abundance of large lipid-rich species in warm years) and forage fish energy density (lower in warm years) [ 25 – 27 ].

Climate change has been increasingly linked with shifts in marine ecosystem processes and structure [ 1 – 3 ]. In addition to long-term global warming [ 4 ] and large-scale modes of climatic variation (i.e. Pacific Decadal Oscillation [ 5 ]; North Atlantic Oscillation and El-Niño Southern Oscillation [ 6 ]), marine heatwaves (MHW)—prolonged periods of elevated sea surface temperatures (SST)–have emerged as a phenomena of ocean climate variability [ 7 – 8 ] that can significantly affect marine ecosystems [ 9 , 10 ]. Although climate change is predicted to alter marine ecosystems globally, the effects of global warming are predicted to be the most extreme at higher latitudes [ 11 ].

Methods

Data collection Beached bird data on St. Paul Island (SPI) including the date, location, taxonomic identification, condition and effort-controlled count were collected by the Aleut Community of St. Paul Island Ecosystem Conservation Office (ACSPI-ECO) in collaboration with the Coastal Observation and Seabird Survey Team (COASST). COASST is a citizen science program in which trained participants conduct monthly standardized surveys, recording all new and previously observed carcasses within prescribed beaches. Field identifications are made from recorded morphological evidence (foot type; standardized body measurements) and consultation with a bird identification guide [46]. Carcasses are individually marked, photographed and subsequently verified by experts using morphological and photographic evidence. Bi-weekly to monthly surveys have been carried out on four 1 km beaches on SPI (Fig 1) since 2006. Additional COASST surveys from nearby St. George Island (~80 km distant; SGI; 2 beaches), as well as throughout the Aleutian Islands (14 beaches on 5 islands), provide a baseline (inclusive of surveys conducted 2006–2015) of effort-standardized carcass abundance and taxonomic composition. During the 2016 MME, extremely high numbers of carcasses and difficult weather conditions necessitated the development of a streamlined protocol (COASST Die-Off Alert). Created and tested by ACSPI-ECO and COASST, the Die-Off Alert protocol requires collection and removal of all carcasses on a set length of beach to a safer location off the beach where they are sorted by species, age class, and carcass condition (i.e. intactness), and photographed in groups (Fig 2). Although primary evidence (e.g. measurements, body condition) is not recorded, the Die-Off Alert protocol facilitates collection over a larger beach area, and gross anatomical features (e.g. intactness, molt) are visible from photographs, enabling post-collection verification. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 2. Photo of carcasses found on North Beach, St. Paul Island, Alaska, on 17 October 2016. Birds pictured are 2 murres (Uria spp.—top row left corner), 8 Horned puffins (Fratercula corniculata—top row center), 2 juvenile Tufted puffins (middle row right side) and 27 adult Tufted puffins (middle and bottom rows). Scale bar in photo is 15 cm total length. https://doi.org/10.1371/journal.pone.0216532.g002

Necropsies Eight intact carcasses (6 Tufted puffins—5 adults, 1 juvenile, and 2 adult Horned puffins–Fratercula corniculata) collected on SPI in October 2016 were sent to the National Wildlife Health Center (USGS). Tissues collected for histopathology (5 birds) were fixed in 10% neutral buffered formalin, processed routinely, embedded in paraffin and sectioned at approximately 5 μm. Routine bacterial cultures of the liver were conducted on five birds. Tracheal and cloacal swabs from all birds were tested for avian influenza [47], and feather pulp was tested for avian paramyxovirus-1 (2 birds) [48] and for West Nile virus (1 bird) [49]. Cloacal contents (4 birds) and stomach content samples (2 birds) were sent to the Wildlife Algal-toxins Research and Response Network (WARRN-West) to analyze for harmful algal toxins domoic acid and saxitoxin using an ELISA (Abraxis, Inc.). Remaining birds could not be tested for algal toxins due to insufficient stomach contents for diagnostic analyses, or decomposition state [50].

Quantifying event versus baseline We conducted two quantitative comparisons of event versus baseline data: lineal carcass encounter rate (ER—carcasses per km of beach surveyed) and taxonomic composition. Monthly baseline ER was calculated as the average for all month-years of data available for the Pribilof Islands (SPI and SGI, 6 beaches). Ranges in baseline estimates were calculated as the bootstrapped 95% confidence interval (95% CI) of the mean for each calendar month. To create a broader geographic comparison, we extended taxonomic comparisons to data collected on the Aleutian Islands (14 beaches on 5 islands).

Particle trajectory modeling To determine the likely origin of carcasses at sea (e.g. catchment area) and to estimate proportional beaching rates, we ran a series of wind-forced particle simulation experiments. Daily releases of 10,000 surface-trapped particles (i.e. replicating a floating bird carcass) were simulated according to wind conditions observed from 1 October 2016 to 24 November 2016, capturing the period where the majority of carcasses were deposited on SPI. Because we had no a priori knowledge of the at-sea distribution of birds prior to the mortality event other than the occurrence of beached birds on SPI, but not on SGI, we randomly generated starting particle locations at a uniform density around SPI, with initial distances, d 0 , up to 100 km. A maximum distance of 100km from SPI was chosen as it results in simulated deposition on SGI, allowing us to identify maximal and closer distributions where deposition on SGI would have been minimal. Previous studies have found strong correlations between carcass deposition and wind-speed and direction [51–54]. We obtained wind velocity fields from the North American Regional Reanalysis (NARR) dataset [55], which consists of 3-hourly averaged grids (32 km resolution). Particle movement from one time-step to another (i.e. 3 hours) was modeled using 4th order Runge-Kutta numerical integration, assuming particle windage equal to 2.5% of the location (linearly interpolated from the nearest 4 NARR grid points), and time-specific wind velocity [51, 52, 56, 57]. Although local surface currents likely contribute to carcass dispersal, we are unaware of surface current information resolved at a suitable temporal (3-hourly) scale to capture nearshore current dynamics. Particle trajectories were simulated from release until intersection occurred with the coastline of SPI or SGI, or 14 days [54], whichever came first. To account for sinking, each particle intersecting either island was assigned a probability of reaching shore, modeled as a logistic function of float duration [54]: (Eq 1) where p(f) is the proportion of carcasses remaining afloat as a function of float duration, f, in hours. Parameters η 1−3 control the shape (i.e. rate and mid-point) of the float function, and determine the rate at which simulated carcass sinking occurs (modeled after [54]). Values for η 1−3 were specified to match observations of carcass float duration from Alaska, where cooler temperatures may delay decomposition, allowing carcasses to remain afloat longer [54, 58]. Ford et al. [58] found median float durations of 7 and 9 days in Prince William Sound, Alaska, with nearly all carcasses sunk by 14 days. We specified two alternate sink functions, with median float durations of 7 and 9 days (S1 Fig), and tested among them to determine how alternate representations of float duration affected our mortality results. For particles remaining at sea for the entire 14 days, p(f) was set to zero.

Catchment analyses To identify catchment area, or the area of ocean within which carcasses could have originated based on observed deposition, we focused on three temporal windows of carcass deposition: (1) 17 to 21 October 2016, (2) 27 October to 1 November 2016, and (3) 15 to 23 November 2016, as these three periods had consistent survey effort, differential patterns of carcass deposition, and differential wind direction (north versus south; S2 Fig). For each window, we took release sets ranging from three days prior to the first MME date (70% of simulated deposition occurs within 3 days–S3 Fig) to the end of the MME window (inclusive) and calculated a grid of proportional deposition (5×5 km grid cells arrayed from SPI to 100 km offshore) by summing p(f) for particles originating in each grid cell that ‘beached’ on SPI and SGI, respectively. In the absence of at-sea distribution data for birds prior to mortality, we assumed that birds were uniformly distributed around SPI out to a maximum distance, d max , which we varied (i.e. by subsetting all particle releases) to investigate alternate spatial distributions. For a given d max we estimated proportional deposition on SPI and on SGI for all particles with initial distance, d 0 ≤ d max , for each of the three observed deposition time windows. To explore the relationship between offshore distribution and island-specific proportional deposition, we varied d max from 2 km to 100 km. Calculating the ratio of proportional deposition on SGI relative to SPI as a function of d max allowed us to bound at-sea distributions by identifying the value of d max at which deposition on SGI would have been expected given deposition on SPI. We also calculated proportional deposition on day i from release j, P j (i), by summing p(f) for particles deposited on day i, from release j. Alternate time-series of expected deposition on SPI and SGI per day (i.e. ∑ j<i P j (i)) were then calculated for alternate values of d max . We explored the likelihood of alternate at-sea distributions (proxied by d max ) by calculating the ratio of simulated deposition on SGI relative to SPI for different values of d max . This allowed us to identify those distributions that would have resulted in minimal expected deposition on SGI, versus those where expected deposition rates would have been comparable among islands. For daily proportional deposition rates we restricted this analysis to time windows when deposition occurred on SPI as these were time intervals when carcasses were known to be afloat, and therefore could have been deposited on SGI.