Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Spatial Distribution and Temporal Patterns of Cassin’s Auklet Foraging and Their Euphausiid Prey in a Variable Ocean Environment

  • Suzanne Manugian,

    Affiliation California Current Group, Point Blue Conservation Science, Petaluma, California, United States of America

  • Meredith L. Elliott,

    Affiliation California Current Group, Point Blue Conservation Science, Petaluma, California, United States of America

  • Russ Bradley,

    Affiliation California Current Group, Point Blue Conservation Science, Petaluma, California, United States of America

  • Julie Howar,

    Affiliation California Current Group, Point Blue Conservation Science, Petaluma, California, United States of America

  • Nina Karnovsky,

    Affiliation Department of Biology, Pomona College, Claremont, California, United States of America

  • Benjamin Saenz,

    Affiliation Department of Biology, University of South Florida, Tampa, Florida, United States of America

  • Anna Studwell,

    Affiliations California Current Group, Point Blue Conservation Science, Petaluma, California, United States of America, Department of Geography, San Francisco State University, San Francisco, California, United States of America

  • Pete Warzybok,

    Affiliation California Current Group, Point Blue Conservation Science, Petaluma, California, United States of America

  • Nadav Nur,

    Affiliation California Current Group, Point Blue Conservation Science, Petaluma, California, United States of America

  • Jaime Jahncke

    jjahncke@pointblue.org

    Affiliation California Current Group, Point Blue Conservation Science, Petaluma, California, United States of America

Abstract

Krill (Euphausiids) play a vital ecosystem role in many of the world’s most productive marine regions, providing an important trophic linkage. We introduce a robust modeling approach to link Cassin’s auklet (Ptychoramphus aleuticus) abundance and distribution to large-scale and local oceanic and atmospheric conditions and relate these patterns to similarly modeled distributions of an important prey resource, krill. We carried out at-sea strip transect bird surveys and hydroacoustic assessments of euphausiids (2004–2013). Data informed separate, spatially-explicit predictive models of Cassin’s auklet abundance (zero-inflated negative binomial regression) and krill biomass (two-part model) based on these surveys. We established the type of prey responsible for acoustic backscatter by conducting net tows of the upper 50 m during surveys. We determined the types of prey fed to Cassin’s auklet chicks by collecting diet samples from provisioning adults. Using time-depth-recorders, we found Cassin’s auklets utilized consistent areas in the upper water column, less than 30 m, where krill could be found (99.5% of dives were less than 30 m). Birds primarily preyed upon two species of euphausiids, Euphausia pacifica and Thysanoessa spinifera, which were available in the upper water column. Cassin’s auklet abundance was best predicted by both large scale and localized oceanic processes (upwelling) while krill biomass was best predicted by local factors (temperature, salinity, and fluorescence) and both large scale and localized oceanic processes (upwelling). Models predicted varying krill and bird distribution by month and year. Our work informs the use of Cassin’s auklet as a valuable indicator or krill abundance and distribution and strengthens our understanding of the link between Cassin’s auklet and its primary prey. We expect future increases in frequency and magnitude of anomalous ocean conditions will result in decreased availability of krill leading to declines in the Farallon Islands population of Cassin’s auklets.

Introduction

In the world’s oceans, climate change has led to widespread shifts in ocean conditions including increases in temperature, decreases in pH, loss of ice and shifts in winds, currents and strength of upwelling events [1]. These changes have led to widespread impacts to the distributions, phenology, and demography of marine species [1]. Under increasingly variable conditions, understanding marine species’ responses is essential to the managers of our ocean resources, allowing them to prioritize areas and species for protection. The objective of this study was to understand how variability in the California Current System (CCS), a highly productive eastern boundary current, impacts the distribution and relationship of an important prey resource, krill (Euphausiids), and one of its predators, the Cassin’s auklet (Ptychoramphus aleuticus).

In the CCS, recent decades have demonstrated increased variability in ocean conditions along central California (e.g., increased incidence of warm water events) and increased frequency of El Niño-like events [2,3]. Also, the nature of El Niño-like events has changed. Sea surface temperatures have increased since 1995 [4,5,6,7,8] though both warming and cooling are occurring in different locations at different scales. Over the past 30 years, local cooling has been documented along the continental shelf in central California [9,10]. While this upwelling-driven cooling may help maintain healthy food-web dynamics, there has been an increase in the frequency and intensity of extreme events, which has reverberated up to upper trophic predators [2,3].

The Cassin’s auklet is a small, zooplanktivorous, central-place foraging seabird in the Northeast Pacific (Baja California, Mexico to the Alaskan Aleutians) that is capable of diving to depth of up to 40 m [11]. Their reproductive success has been monitored on Southeast Farallon Island (SEFI) since 1972. SEFI is the main breeding colony within the central CCS surrounded by pelagic waters that fall within two National Marine Sanctuaries: Gulf of the Farallones National Marine Sanctuary (GFNMS) and Cordell Bank National Marine Sanctuary (CBNMS). In the past decade extreme variation in breeding occurred: successful double-brooding was recorded in some years while complete colony abandonment was observed in 2005–2006 on SEFI [12,13,14,15,16]. To gain insight into the link between ocean conditions and reproductive success we examined the distribution and abundance of prey in relation to the distribution and abundance of Cassin’s auklets foraging at sea. Establishing the link between foraging seabirds and their environment is critical to understanding how marine food webs are affected by climate change. Seabirds are highly sensitive to localized changes in oceanographic conditions, prey availability, and ecosystem dynamics [17,18,19,20].

The CCS is a highly productive eastern boundary current system dominated by upwelling in the spring and summer months and is home to a large array of marine birds and mammals [14,21,22,23,24]. This upwelling system is replete with mid-trophic level organisms, linking primary producers and upper-trophic level organisms [25]. Euphausiids (family Euphausiidae, known commonly as krill) are marine crustaceans and important mid-trophic level organisms. The euphausiid assemblage in the central CCS is dominated by two species, Euphausia pacifica and Thysanoessa spinifera, both of which show sensitivity to the strength and duration of upwelling [26]. Cassin’s auklets rely on these mid-trophic level prey and therefore are indirectly influenced by ocean forcing [12,27,28,29,30,31]. Both krill species exhibit diel vertical migration [22,32,33,34], although T. spinifera has been documented to swarm at the surface during daylight hours within the Sanctuaries [35].

Within the CCS, zooplankton are often aggregated in hotspots, defined as persistent locations where birds and mammals prey on them [36,37,38]. In the central CCS, the greater Gulf of the Farallones region is one of these persistent hotspots [36]. Here, seabirds are associated with specific bathymetric and oceanographic features such as islands and the shelf-break that play a role in aggregating krill [37,38,39, 40,41,42,43]. However, birds’ associations with many of these hotspots vary interannually [38, 43,44] and are often related to the onset, timing, and strength of upwelling.

We hypothesized (1) that the distributions of both krill and Cassin’s auklets are influenced by local bathymetry and local and basin-scale oceanographic conditions, (2) Cassin’s auklet at-sea distribution is correlated with measures of prey biomass in the upper water column, and (3) Cassin’s auklets utilize consistent hotspots where their krill prey are aggregated.

Methods

To evaluate our hypotheses, we (1) modelled how Cassin’s auklets’ at-sea distributions were related to in situ oceanographic, large-scale climate, and bathymetric variables (2004–2013) and (2) examined how their at-sea spatial distributions were related to the distribution of krill detected acoustically (2004–2013). We also (3) examined the zooplankton that Cassin’s auklets collected for their chicks to establish whether what they were taking was the same prey species as what was caught in nets (on SEFI, 1985–2013), (4) used Time-Depth-Recorders (TDRs, 2008–2013) to establish the diving depths of foraging Cassin’s auklets, and (5) collected zooplankton from hoop nets (2008–2011) to determine what was the composition of the acoustically detected prey biomass. Together, these analyses characterize the ecological linkages between Cassin’s auklets and their prey.

Ethics Statement

All TDR and diet collection work was conducted on the Farallon National Wildlife Refuge with permission and approval of the U.S. Fish and Wildlife Service after thorough review of the project proposal. Auklet diet collection has occurred annually since the 1970s as part of the regular monitoring on the Farallones while TDR deployments have been done annually since 2008. The study was conducted under the terms of Cooperative Agreement #81640-5-J046 (May 2005 –September 2009) and Cooperative Agreement #81640AJ008 (October 2009 –September 2014) and was in agreement with all required state and federal collecting permits. Point Blue Conservation Science’s Federal Banding Permit (#09316) specifically authorizes application of TDRs to Cassin’s auklets and Rhinoceros auklets though field studies only involved Cassin’s auklets which do not have any special protected status. All sampling procedures and manipulations (including deployment/retrieval of TDRs and diet sampling) have been specifically reviewed and approved by U.S. Fish and Wildlife Service and the Farallon National Wildlife Refuge manager. The federal Bird Banding Laboratory (BBL) also reviewed and approved a proposal for the TDR work and subsequently amended Point Blue Conservation Science’s banding permit to specifically include the use of auxiliary markers (including TDR, GPS and GLS tags) on Cassin’s auklets. Personnel concluded that a compatibility determination was not necessary because of the low probability of impact (distress) to the sampled Cassin’s auklets, the small number of birds sampled, the project’s short duration, lack of impact to other island natural resources, and negligible impacts on accomplishing other island duties, personnel resources, and FWS resources.

Study area and survey design

This study used ten years of survey data collected in GFNMS and CBNMS (Fig 1) as part of the Applied California Current Ecosystem Studies (ACCESS) program (www.accessoceans.org), a partnership of GFNMS, CBNMS, and Point Blue Conservation Science. Data were collected between April and October 2004–2013 on three vessels: the R/V John H. Martin (Moss Landing Marine Laboratories, Moss Landing, California, USA), the R/V McArthur II (National Oceanic and Atmospheric Administration (NOAA), Seattle, Washington, USA), and the R/V Fulmar (NOAA, Office of National Marine Sanctuaries, West Coast Region, Monterey, California, USA).

thumbnail
Fig 1. Study area within the Gulf of the Farallones and Cordell Bank National Marine Sanctuaries, California, USA, showing ACCESS program transect lines and oceanographic sampling stations visited for CTD casts and zooplankton tows.

https://doi.org/10.1371/journal.pone.0144232.g001

ACCESS conducts multidisciplinary research and monitoring to understand changes in marine wildlife populations and how wildlife responds to changing ocean conditions in the Sanctuaries. ACCESS cruise transect lines primarily ran east-west through the Sanctuaries to sample the area between southern Bodega Bay (38.3°N) and San Pedro Rock (37.6°N); the amount of data collected varied by cruise and year as a function of time constraints and inclement weather. Transect lines 1–7 were used in these analyses and each transect was subdivided into 3 km bins by cruise using methods from prior studies [36,38,42]. Each bin’s centroid was subsequently assigned in-situ oceanographic and remotely-collected environmental data, biological data (seabird and krill), and location-related data (e.g., distance to certain features, depth, and contour index). Finally, we discarded end-of-line bins (< 1 km in length) from analyses to avoid introducing errors [38,45,46].

We created a prediction matrix of 1 km2 cells based on the Sanctuaries’ boundaries buffered externally by 5 km. We used this matrix to make predictions using selected final models for krill and Cassin’s auklets incorporating environmental information from the 3 km bins. The scale of the matrix allowed for highly detailed visualization of model results.

Environmental variables: oceanographic variables and oceanographic surfaces

Sea surface temperature (SST), sea surface salinity (SSS), and sea surface fluorescence (SSF) were recorded in situ by a thermosalinograph (TSG) located on the hull of each vessel. SSF is a measure of chlorophyll a, a proxy for the phytoplankton used as forage by zooplankton species, such as krill. We predicted missing sea surface values, a result of instrument malfunction, using data collected in situ at pre-determined, permanent stations using a Conductivity-Temperature-Depth (CTD) profiler (SBE 19plus SeaCAT Profiler CTD with WETStar fluorometer). The closest permanent station to each missing TSG value’s associated 3 km bin was used. Cruises missing all TSG values for a variable were predicted from the CTD station values corresponding to that cruise. To fill in missing TSG values, we used a regression equation incorporating observed TSG values from all ACCESS cruises and the associated CTD values while controlling for month and year. We generated separate models for predicting average fluorescence (R2 = 0.63), salinity (R2 = 0.88), and temperature (R2 = 0.94), all p < 0.001, and used those values to complete the dataset.

Using ordinary kriging with optimization, we created interpolated surfaces in ArcGIS (10.2, ESRI Redlands, CA) for the oceanographic variables (SST, SSS, SSF) to predict oceanographic conditions in locations not sampled within and extending 5 km outside the Sanctuaries’ boundaries (n = 106). Kriging fits a spatial-dependence model to a set of data and generates a landscape of values of the measured sample points around a predicted location providing a measure of accuracy to the predictive surface [47,48,49]. We employed parameter optimization to standardize the kriging process in order to minimize mean square error for landscape predictions and we applied a default search radius to weight measured locations. Kriging can be sensitive to anisotropy (i.e., can vary with latitude and/or longitude) but there was a spatial trend in some of the oceanographic measurements given the large extent of our study area. To check for raster surface outliers, we compared each variable’s set of surfaces’ average standard errors to ensure all values were within three standard deviations of the mean. We excluded one cruise because of lack of data (n = 3 surfaces) and only used values extracted from the study region south of ACCESS transect line 3 to the southern boundary of the region for two others (n = 6 surfaces). We then sampled all generated surfaces to extract predicted oceanic variable values for each bin and for each 1 km2 prediction cell within the Sanctuaries.

Environmental variables: climate indices, bathymetric variables, and distance variables

Climate indices.

In our modeling, we included four variables that fluctuate on multiple timescales and which are known to influence ocean conditions in the CCS: the (1) Pacific Decadal Oscillation (PDO) which is the first dominant mode of variability in sea surface temperature north of 20°N [50]; (2) North Pacific Gyre Oscillation (NPGO) which is the second dominant mode of variability in sea surface height in the northeast Pacific [51]; (3) Southern Oscillation Index (SOI) which is a standardized index based on the observed sea level pressure differences between Tahiti and Darwin, Australia [52]; and (4) Upwelling Index (UI) which is the amount of upwelled water from the base of the Ekman layer per 100 m of coastline (Table 1). We used as a single averaged monthly UI value from two locations along the California coastline (36°N 122°W and 39°N 125°W) to best match the study region. Monthly values of all climate indices were assigned to cruise by month. We included climate lags in monthly increments in modeling to examine the possible influence of ocean conditions up to three months prior to each cruise month.

thumbnail
Table 1. Description and ranges of variables used to model krill distribution and Cassin’s auklet abundance using 3 km bins (n = 4313).

https://doi.org/10.1371/journal.pone.0144232.t001

Bathymetric and distance variables.

For each observation point along ACCESS cruise transects and each 1 km2 prediction cell midpoint, we calculated distance to important locations and depths. Important locations included nearest mainland feature, the island breeding colony SEFI, Cordell Bank, and the 200 m isobath (or shelf-break). Distances were calculated in ArcGIS from the midpoint of each 3 km bin or 1 km2 prediction cell. For each bin and cell, we also calculated average depth using a bathymetric surface (available from California Department of Fish and Wildlife) and generated a contour index value, a metric of ocean floor roughness: (maximum depth—minimum depth) / maximum depth [42].

Krill hydroacoustics surveys.

Zooplankton was sampled acoustically every 2 s along transects using a Simrad EK60 echosounder. Raw volume backscatter data were post-processed and integrated into 200 m horizontal by 5 m depth bins using Echoview software (Sonardata, Pty. Ltd.). To avoid surface interference, backscatter from 0–5 m were excluded. The level of noise reduction was adjusted by visually comparing echograms until noise at depth appeared equal [53] and time-varied noise was subtracted from each echogram before integration into bins.

Watkins & Brierley [54] found that by selecting for a decibel difference between 120 kHz and 38 kHz of 2–16 dB, portions of the echogram could be classified as Antarctic krill (Euphausia superba). We evaluated the Stochastic Distorted-Wave Bourne Approximation (SDWBA [55,56]) for models of Euphausia pacifica generated from images of body shape and published animal densities [57], and found no significant difference between the E. pacifica and E. superba models. We therefore used the simplified SDWBA [58] to generate krill acoustic target strengths for 8–30 mm length krill, roughly the size-range of juvenile, immature, and adult stages of the euphausiid species of interest. Based on these target strengths, we classified echogram bins as krill using a decibel difference of 11–19 dB. Echogram bins that fell within the 11–19 dB range were then used as a positive mask over the original 120 kHz echogram data for further processing. This method of classification functions well in single species aggregations; however, in the CCS, there are many more types of zooplankton that may be similar acoustically to krill, adding to classification uncertainty. Furthermore, the sampled area may contain not only other zooplankton but also fish, altering the average target strength sampled by the echosounder. Despite these limitations, we believe the dB difference method of krill classification is more accurate than manually classifying aggregations in the echograms.

The dB differencing operation resulted in the inclusion of some layers we judged not to be krill aggregations, due to either weak scattering (e.g. < -75 dB) or the position and shape of the scattering layer in the water column (e.g. too shallow and/or diffuse). Tucker trawls were deployed at the shelf-break down to 200 m depth on the 3 designated CTD lines to sample the zooplankton and ground-truth the acoustic record. Additionally, we manually included several bins excluded from analysis during the dB difference operation. These bins appeared at the center of dense krill swarms, but were excluded due to a large difference between backscatter at 120 and 38 kHz. We theorize that the most dense krill swarms may not exhibit pure Raleigh scattering at 120 kHz, causing backscatter values that fell outside of the 11–19 dB differencing classification. These manual classification adjustments accounted for less than 0.5% of volumetric acoustic returns ultimately attributed to krill.

Krill biomass (g m-2) was calculated by first apportioning the acoustic backscatter to a length-frequency distribution of individual krill determined from trawl surveys concurrent to acoustic sampling. The backscatter was divided by the corresponding individual backscattering cross sections and then multiplied by a length-weight relationship to obtain the mass of krill [53]. The length-weight relationship was derived from sampled individual krill ML = 2.0 x 10-4L2–1.7 x 10-3L, where L is an individual krill length (mm). Acoustic data were integrated throughout the water column to 30 m depth and binned to 3 km bins to match other datasets (n = 3046; maximum krill biomass: 24,628 g m-2).

Cassin’s auklet surveys

Seabird observations were systematically collected during each ACCESS cruise from the vessel flying bridge using standardized strip-survey methods while on effort, defined as transiting a transect line with an observer actively recording seabird observations during daylight hours with vessel speed of ~10 knots [14,59]. Data on factors that affect detection (Table 1) were noted and a single observer recorded encounters of seabirds within a 90° arc from the bow to the vessel’s starboard side and within 50–300 m of the vessel (following methods in Jahncke et al. [14] and McGowan et al. [38]). Seabirds’ behaviors were recorded as flying, foraging, ship-attracted, and sitting. Only actively foraging birds or birds sitting on the water (assuming birds fed or were about to feed in the vicinity of where they were resting) were included in analyses [19]. The total number of Cassin’s auklets counted along each transect line was summed for each 3 km bin and assigned to that bin’s midpoint (n = 3879; 3047 zero-count bins, 832 non-zero count bins).

Modeling: krill and Cassin’s auklets

We used Generalized Linear Modeling (GLM) to examine how krill distribution and abundance were related to Cassin’s auklets at-sea distribution and abundance. We treated krill biomass and seabird counts as separate dependent variables and developed two distinctive predictive models (one for each organism) in the Sanctuaries (STATA version 13.0, StataCorp 2013, Statistical Software, College Station, TX). To predict Cassin’s auklet abundance, we used a zero-inflated negative binomial regression model [60]. This approach is well-suited for count data, and allows for the presence of “false zeroes” [60]. For modeling krill biomass at depths of up to 30m (due to lack of auklet foraging below this depth) we fit a two-part model, using the procedure twopm in STATA 13.0. The first part was a logistic regression model estimated for probability of observing positive value (presence of krill) vs. zero. The second part was a negative binomial model that was conditional on a positive value in the first part of the model. Thus it is a type of hurdle model, allowing one to separately model presence vs. absence from abundance, where present [60]. We used twopm because it allows for a positive outcome that is continuous, since the dependent variable is krill biomass per unit volume, rather than a count.

Univariate testing and initial modeling.

For krill modeling, we conducted univariate testing of variables of interest (Table 1) for both parts of the two-part model to determine whether relationships were linear or not. We generated a full multivariable model including, for each modeled variable, the most significant functional relationships in each part of the model (determined using Likelihood Ratio Statistics (LRS)). Depending on the variable, we used one of several strategies. Significant quadratic terms (p < 0.05) were always included; if the quadratic term was non-significant, the linear form of the variable (provided p < 0.20) was included in this initial phase. For distance variables, we chose the most significant of four possible relationships to krill based on a priori knowledge: linear, quadratic, logarithmic, or inverse logarithmic. For climate indices, we selected either a no-, 1-, 2-, or 3-month lag index by testing each variable for each of the lag periods and selecting the lag with the greatest significance. For average depth, we chose the most significant of eight possible relationships: linear, quadratic, cubic, logarithmic, quadratic logarithmic, inverse logarithmic, quadratic inverse logarithmic, or square root.

We then included the selected form of all variables in a full multivariate model and used manual backwards stepwise removal until all remaining variables were significant (p < 0.05), either in the first or second part of the final two-part model. The exception to this rule was with respect to average depth. For Cassin’s auklets, the cubic transformation was preferred (p < 0.0001 for the cubic term). For krill abundance (part 2 of the two-part model), we forced in the cubic term (p > 0.7), to allow for an asymmetric polynomial fit (quadratic transformation alone imposes a symmetric pattern); dropping the cubic term would have yielded an artifactual spatial pattern. Thus, both krill and Cassin’s auklet models allowed for a cubic term for depth.

Zero-inflated modeling.

For Cassin’s auklet modeling, we first used negative binomial regression to test all variables of interest using the same modeling guidelines outlined above. This regression is recommended for data that have a larger variance than expected from a Poisson-distributed dataset [61,62,63]. However, to account for an excess of zeroes compared with a negative binomial distribution, we used zero-inflated negative binomial regression [60].

Zeroes in count data may reflect observer detection biases and ignoring them can lead to poor modeling of a species’ distribution. Using a model that accounts for both “true zero counts” and “false zero counts” is needed to obtain accurate predictions [61,63]. There are four ways zeroes can arise in ecological data: the species is not present because of an unsuitable habitat; the species is not seen because it does not fill the entire suitable habitat; the species readily occurs in the habitat but is not present in the surveyed location at the time of survey; and the species is present during the survey but individuals are not detected due to observer error [64]. In an example of this last case, a failure to detect an individual might be caused by poor observation conditions that prevent the observer from seeing the bird even though it was present.

For both the krill and Cassin’s auklet datasets, variables were included to standardize sampling efforts. As a measure of effort in krill modeling exercises, we used the total number of acoustically-sampled vertical cells within each bin. An offset of the logarithm of the area surveyed within each bin was used for the Cassin’s auklet dataset [38,45]. We tested variables affecting detection probability (i.e., variables used to account for excess of recorded zeroes) in the Cassin’s auklet modeling (Table 1). Bias variables of statistical significance (p < 0.05) were left in the model using previously mentioned manual backwards step-wise removal. Finally, we used the Vuong test statistic (STATA version 13.0, [65]) to test if using the zero-inflated negative binomial regression was superior to the standard negative binomial regression [60].

Once we generated final models, we examined potential interactions between year and specific variables based on a priori knowledge. We tested for interactions of year with UI, salinity, and temperature for krill and with UI and salinity for Cassin’s auklets and built models incorporating significant interactions, if any, as determined by the likelihood ratio test. We also considered allowing for multiple interactions in a final model and report the significance of those interactions with respect to the likelihood ratio statistic. Variance inflation factor (VIF) testing for multicollinearity [63] was conducted among variables in the final models.

Model validation.

We validated model fit using k-fold cross validation (k = 10, with 20 replicate runs) based on the predictions to the 3 km bins. One at a time, each one of the 10 subsets was used as a “test” dataset and the other nine were combined and used as a “training” dataset. For each run, we calculated a median result; we repeated the procedure 20 times and calculated a median of the 20 medians. We compared the final model pseudo R2 value with the median pseudo R2 value and reported the comparison. A model validation pseudo R2 value within 40% of the final model pseudo R2 value was considered a minimally adequate model fit [36].

Mapping predators, prey, and their overlap

We predicted predator (Cassin’s auklet) and prey (krill spp.) distributions in each 1 km2 cell using model results and inputs that included observed oceanographic variables (collected during surveys), climate indices (collected remotely), and bathymetric- and distance-related data from 3 km bins (Fig 2). We averaged spatial occurrence predictions separately by month (across all years 2004–2013) and by year (May through September) for each cell of the prediction matrix. We excluded April because it was not sampled as frequently as other months during the study period. We mapped each month or year independently and scales were chosen using percentiles of predicted values from the full predicted monthly or annual dataset to allow for same-scale comparison for each species’ maps. We identified areas of high use habitat for both individual models and combined species in all-inclusive hotspot maps. We mapped areas of co-occurrence, examining the locations of top 5-, 10-, 20-, and 30-percentile habitat use by both krill and Cassin’s auklets.

thumbnail
Fig 2. Schematic study modeling methodology showing binned and environmental data, predicted to 1 km2 prediction cells, linked with Cassin’s auklet diet data, Cassin’s auklet TDR dive data, and zooplankton tow samples from the upper water column.

https://doi.org/10.1371/journal.pone.0144232.g002

Cassin’s auklet diet data

We determined the diet of Cassin’s auklets feeding chicks by capturing adults returning to SEFI at night and collecting regurgitations intended to be fed to their young. Cassin’s auklets carry partially digested prey to their chicks in a sublingual pouch [66]. Adult Cassin’s auklets were captured by hand as they returned to the colony at night. 10 random diet samples were collected each week throughout the chick-rearing period (May to July / August) for 1985 to 2013, excluding 1992. Immediately upon capture, the Cassin’s auklet was inverted with its head over a sterile Whirl-Pak specimen bag and collected samples were frozen until the diet composition analysis could be performed. In most cases, diet samples were obtained when the birds spontaneously regurgitated the contents of their sublingual pouch. If spontaneous regurgitation did not occur, we attempted to induce regurgitation by gently massaging the sublingual pouch for a few seconds. Birds were subsequently released on the ground where they were initially captured immediately after a sample was obtained. See Sydeman et al. [12,18] and Abraham and Sydeman [31,67] for details on diet collections and analysis. For this paper, we show prey composition as the annual proportion of the taxa in the diet, as determined by number. While proportion by number inflates the importance of smaller prey items, we used this measure as we had counts of prey items for all years, while we only had wet weights beginning in 1994 (Point Blue, unpublished data). Results using wet weights were similar but are not reported as they require assumptions about prey masses. We used one-way ANOVA to test for differences in proportion of euphausiids vs. all other taxa (arc-sin transformed) found in the diet in each sample amongst years.

Cassin’s auklet dive records

We used TDR records to provide justification for use of an acoustic record in the upper 30 m of the water column. Cassin’s auklets at SEFI in the Farallon National Wildlife Refuge were outfitted with TDRs between 2008 and 2013; procedures followed those of Karnovsky et al. [68]. TDR deployment was conducted concurrently with ACCESS cruises to fully explore the link between birds, their prey, and ocean conditions. The TDRs (CEFAS Technology Limited, United Kingdom) were approximately 8 mm diameter x 31 mm length in size, weighed 2.7 g in the air (or 1.5% of the average initial mass of instrumented birds), and had a cross sectional area of 0.50 cm2 (or between 1 and 2% of instrumented birds’ cross sectional area). They were programmed to record temperature and pressure every 5 s while the bird was on the surface or at the colony and every 0.5 s whenever the bird dove below 1.5 m. Instruments were cylindrical in shape and had a rounded tip to minimize water resistance [69]. The instrument was attached to the ventral side of the bird by gluing it to the feathers of the belly using a few drops of Loctite adhesive. We deployed TDRs (n = 73) on adult Cassin’s auklets which were brooding young chicks in artificial nest boxes. We attempted to balance the sex ratio to account for any potential biases. We used one-way ANOVA to test for differences in dive depths amongst years.

To evaluate potential negative impacts on Cassin’s auklet breeding biology of TDRs, we compared reproductive success (RS; chicks fledged) of long-term study birds (non-instrumented; Point Blue, unpublished) with RS of pairs where one adult was fitted with a TDR on SEFI between 2008 and 2013 using logistic regression controlling for year as a categorical variable. Additionally, we compared the growth rates of chicks of long-term study birds and pairs where one adult was fitted with a TDR on SEFI for 2008 using a 2-tailed t-test.

Zooplankton tow data

In order to evaluate the potential available prey for Cassin’s auklets, we examined zooplankton data from daytime surface tows. We sampled zooplankton during ACCESS cruises July 2004 through September 2011 in daytime hours (0600–1800 hours) using a 1 m diameter hoop net fitted with 333 μm mesh. Oblique tows (50 m to the surface) were performed at each CTD station, towing the net at 1–2 knots for about 10 min. Volume filtered was estimated from a mechanical flowmeter (General Oceanics model 2030RC) attached to the opening of the hoop net frame. The contents of each tow were preserved in 10% formalin solution and stored in 1 L plastic jars.

Preserved samples were identified and enumerated at the Institute of Ocean Sciences, British Columbia, Canada. Zooplankton identification protocols were consistent with previous zooplankton analyses [70]. In brief, entire samples were scanned for abundance of large and/or rare taxa. Using a Folsom splitter, the sample was quantitatively sub-sampled for counts of small abundant zooplankton taxa. Total counting effort was ≥ 400–500 individuals per sample, sufficient to give an expected subsampling error of ≤ 20% for the dominant species [70,71]. We used a chi-square test (χ2, Pearson’s) to determine differences of adult and juvenile krill counts in hoop net samples amongst years.

Results

Acoustically-sampled krill: habitat association and predicted distribution

We found that all tested variables were significant predictors for the amount of krill in one or both parts of the two-part model. SST, SSS, and SSF all significantly predicted presence for krill. All climate indices (SOI, PDO, NPGO, and UI) were significant, with varying month lags. Bathymetric and all distance variables were also significant predictors of the amount of krill observed (Table 2, S1 Table). No significant multicollinearity was found between variables in the final model (VIF < 10). Cross validation (k = 10) pseudo R2 was within 40% of krill final model pseudo R2 (Table 2).

thumbnail
Table 2. Krill and Cassin’s auklet model results with transformation (L = linear, Q = quadratic, C = cubic) and coefficient sign for significant variables.

Predictive models used a two-part model combining logistic and negative binomial regressions for krill biomass (g m-2; n = 3046) and a zero-inflated negative binomial regression for Cassin’s auklet counts (birds km-2; n = 3879). The bold variables represent significant (p < 0.01, Likelihood Ratio Test) interactions of that variable with year and coefficient sign refers to the main effect without any interactions.

https://doi.org/10.1371/journal.pone.0144232.t002

The influence of oceanic processes varied among years with respect to the amount of krill present in the upper 30 m as shown by the inclusion of two significant year-interaction terms for both parts of the two-part species model: 3-month lag UI (LRS = 437.17, df = 26) and SST (LRS = 414.63, df = 35). The influence of processes driven by changes in climate on predicted krill biomass varied among indices: the localized influence of upwelling (UI (LRS = 547.19, df = 29)) was an order of magnitude greater than climate variables that operate on a basin-wide scale [PDO (LRS = 26.0, df = 1), SOI (LRS = 50.12, df = 1), and NPGO (LRS = 74.38, df = 2)]. Krill habitat use varied across months (Fig 3) and across years (Figs 4 and 5) with June and 2007–2008 exhibiting the least predicted krill biomass within the Sanctuaries. Krill biomass varied with climate processes (Fig 5). Krill hotspots were consistently predicted across all months and years in the areas surrounding SEFI, southeast of SEFI, and Cordell Bank (Fig 6).

thumbnail
Fig 3. Predicted krill biomass (g m-2), 2004–2013: variations in modeled habitat use for May (a, predictions averaged over n = 8 years), June (b, n = 7), July (c, n = 7), and September (d, n = 8); each gradation represents a decile, monthly mean values with all years pooled.

The 200 m isobath and boundaries of National Marine Sanctuaries included in blue and black, respectively.

https://doi.org/10.1371/journal.pone.0144232.g003

thumbnail
Fig 4. Predicted krill biomass (g m-2), 2004–2013 (a-j, sequentially): variations in annual modeled habitat use (predictions averaged over n = 3 months except 2010, n = 4, and 2013, n = 2); each gradation represents a decile, annual mean values with all months pooled per year.

The 200 m isobath and boundaries of National Marine Sanctuaries included in blue and black, respectively.

https://doi.org/10.1371/journal.pone.0144232.g004

thumbnail
Fig 5.

Annual variation in predicted (a) Cassin’s auklet abundance (birds km-2) and (b) krill biomass (g m-2), yearly averaged values for (c) the Southern Oscillation Index (SOI), (d) the Pacific Decadal Oscillation (PDO), (e) the North Pacific Gyre Oscillation, and (f) Upwelling Index within the Gulf of the Farallones and Cordell Bank National Marine Sanctuaries by year (2004–2013) with all months pooled.

https://doi.org/10.1371/journal.pone.0144232.g005

thumbnail
Fig 6. Predicted (a) krill biomass (g m-2) and (b) Cassin’s auklet abundance (birds km-2), 2004–2013: top habitat use model applied across all months and years; each gradation represents a decile.

The 200 m isobath and boundaries of National Marine Sanctuaries included in blue and black, respectively.

https://doi.org/10.1371/journal.pone.0144232.g006

Cassin’s auklets: habitat association and predicted distribution

Sea surface salinity (SSS) was significantly associated with Cassin’s auklet distribution along with three climate indices with varying time lags: SOI, PDO, and UI. The influence of SSS on individuals also varied annually; we included the interaction of that oceanographic variable with year (LRS = 102.43, df = 9). Overall, Cassin’s auklets were found in regions of higher salinity. Additional significant predictive variables were distances to SEFI and the 200 m isobath, contour index, and average depth (Table 2, S2 Table). Birds were found closer to the shelf break in areas of higher bathymetric relief and in association with shallower depths near SEFI. From all tested detection variables, only sea state and visibility significantly contributed to the Cassin’s auklet zero-inflation model (p < 0.01). The influence of climate on predicted Cassin’s auklet counts varied among indices: those that operate on a basin-wide scale [PDO (LRS = 89.94, df = 2) and SOI (LRS = 74.50, df = 2)] were an order of magnitude greater than that of a localized process [UI (LRS = 5.69, df = 1)]. No significant multicollinearity was found between variables in the final model (VIF < 10). Cross validation (k = 10) pseudo R2 was within 40% of Cassin’s auklet final model pseudo R2 (Table 2).

We found that hotspot locations of Cassin’s auklets varied across months (Fig 7) and across years (Figs 5 and 8). Like krill, predicted values of Cassin’s auklet counts varied with climate processes (Fig 5). Auklets were present in the study area during all surveys reflecting use of the Sanctuaries across all surveyed months. Auklet occurrence was predicted to be high in 2004 and 2010, but low in 2007–2008. Cordell Bank and the area to the east of the 200 m isobath along the Farallones Escarpment were highlighted as areas where birds consistently foraged (Fig 6).

thumbnail
Fig 7. Predicted Cassin’s auklet abundance (birds km-2), 2004–2013: variations in modeled habitat use for May (a), June (b), July (c), and September (d); each gradation represents a decile, monthly mean values with all years pooled.

The 200 m isobath and boundaries of National Marine Sanctuaries included in blue and black, respectively.

https://doi.org/10.1371/journal.pone.0144232.g007

thumbnail
Fig 8. Predicted Cassin’s auklet abundance (birds km-2) by year, 2004–2013 (a-j, sequentially): variations in annual modeled habitat use; each gradation represents a decile, annual mean values with all months pooled per year.

The 200 m isobath and boundaries of National Marine Sanctuaries included in blue and black, respectively.

https://doi.org/10.1371/journal.pone.0144232.g008

Predator and prey co-occurrence

We examined areas of high predicted abundance for both krill and Cassin’s auklets across all months and years, combined. Areas of most likely spatial co-occurrence (as indicated by top 20% use by both species) were the Farallones Escarpment surrounding SEFI and over Cordell Bank (Fig 9); however, these areas varied both annually and by month (Figs 10 and 11).

thumbnail
Fig 9. Intersection of high predicted spatial co-occurrence in red for krill and Cassin’s auklets combined over May-July and September 2004–2013 at top 20% habitat use for each predictive model.

Consistent high use areas (Cordell Bank and the Farallones Escarpment) are noted. 200 m isobath and boundaries of National Marine Sanctuaries included in blue and black, respectively.

https://doi.org/10.1371/journal.pone.0144232.g009

thumbnail
Fig 10. Intersection of high predicted spatial co-occurrence in red for krill and Cassin’s auklets for May (a), June (b), July (c), and September (d) at top 20% habitat use for each predictive model.

200 m isobath and boundaries of National Marine Sanctuaries included in blue and black, respectively.

https://doi.org/10.1371/journal.pone.0144232.g010

thumbnail
Fig 11. Intersection of high predicted spatial co-occurrence in red for krill and Cassin’s auklets 2004–2013 (a-j, sequentially) at top 20% habitat use for each predictive model.

200 m isobath and boundaries of National Marine Sanctuaries included in blue and black, respectively.

https://doi.org/10.1371/journal.pone.0144232.g011

Cassin’s auklet diet

Euphausiids dominated the diet of the Cassin’s auklet in 22 of the 27 years of prey composition data from SEFI (Fig 12). We found significant differences in euphausiid consumption among all years analyzed (ANOVA F(27,2351) = 20.40, p = 0.000). While sample size varied among years (Fig 12), the few samples collected in 2005–2006 (n = 9 and n = 6, respectively) contained almost entirely mysids.

thumbnail
Fig 12. Annual proportion of major zooplankton taxa, by number, in the chick provisioning diet of Cassin’s auklets on Southeast Farallon Island, 1985–2013.

Numbers above the bars indicate sample sizes.

https://doi.org/10.1371/journal.pone.0144232.g012

Cassin’s auklet diving behavior

The maximum dive depth recorded by a Cassin’s auklet was 39.5 m, although most dives were shallower with an overall median dive depth of 7.9 m (Table 3). We found significant differences in dive depths among the 6 years of data collection (ANOVA F(5,57510) = 1371.64, p = 0.000). While 55% of tagged birds’ records (n = 40) had at least one dive greater than 30 m, 99.5% of all dives across years (n = 57,516) were shallower than 30 m (n = 57,216).

thumbnail
Table 3. Summary of TDR records (n = 73) overall and by year (2008–2013).

https://doi.org/10.1371/journal.pone.0144232.t003

Attaching TDRs to Cassin’s auklets had no significant influence on breeding success (chicks fledged per breeding pair) or growth rates of chicks where one adult was fitted with a TDR (chick weight change g day-1). Between 2008 and 2013, breeding success of TDR-outfitted birds compared to long term study birds was higher for 5 of the 6 years (Point Blue, unpublished data; TDR birds: success = 0.887, SD = 0.318, n = 106; Non TDR birds: success = 0.780, SD = 0.318, n = 255). In 2008, there was no statistical difference between chick growth rates for the two groups (p = 0.86); chicks from nest sites with non-instrumented parents grew at an average of 5.25 g day-1 (SE = 0.10; n = 40) while chicks from nest sites with one adult fitted with a TDR grew at an average of 5.37 g day-1 (SE = 0.074; n = 10).

Zooplankton tows

We found that adult krill (combined E. pacifica and T. spinifera) were present in 20.38% of the surface hoop net samples, and juveniles were present in 42.68% of samples in the eight year time series (Fig 13). We found a significant difference in krill occurrence in hoop net samples for both adults and juveniles among the eight years of data (adults: χ2 = 42.3451, df = 7; juveniles: χ2 = 62.9727, df = 7). Abundance (shown through percent frequency of occurrence) varied among years, especially for juveniles; adults ranged from 5.00% (in 2004) to 37.33% (in 2010), and juveniles ranged from 19.12% (in 2004) to 70.73% (in 2011). These tows suggest krill were frequently available in the upper 50 m of the water column.

thumbnail
Fig 13. Percent frequency of occurrence of E. pacifica and T. spinifera adults and juveniles in surface hoop net samples (2004–2011).

https://doi.org/10.1371/journal.pone.0144232.g013

Discussion

Within the years of this study we found variability in oceanic conditions, krill distribution and abundance, and responses of the Cassin’s auklet. Both krill and Cassin’s auklets aggregated along the shelf break in areas of productive upwelled waters; their distribution was influenced by local ocean processes and basin-wide climate conditions. Though hotspots varied temporally, we generally found persistent predicted krill habitat in the Cordell Bank region and around SEFI over the study period (Figs 3, 4 and 6). Krill were concentrated away from the mainland, close to the Farallones Escarpment and Cordell Bank in areas of higher bathymetric complexity. Our results corroborate previous studies that found a large center of abundance of E. pacifica off central California [72], often in the form of long narrow bands along the outer shelf [43]. Additionally, other research found large densities of auklets within 25 km of SEFI and along the shelf break [43]. These areas have been proposed by other studies to be krill hotspots important to seabirds [38,43,73] and new studies (including the present one) have corroborated the high availability of krill in these areas [37].

Foraging Cassin’s auklets, like krill, were associated with Cordell Bank, the shelf break, and the Farallones Escarpment, and we found seasonal and interannual variability associated with these features (Figs 6, 7 and 8). Seasonal variation showed that the shelf break west of SEFI was consistently important during all surveyed months, while Cordell Bank was most important during the post-breeding season. Auklets’ annual use of the shelf-break region west of SEFI and the Farallones Escarpment also varied. Cordell Bank was important to auklets in half of the surveyed years, notably 2005–2006 and 2010–2011 (2005 displayed particularly low krill density while 2010 and 2011 displayed high density).

Localized and basin-scale processes driven by changes in climate affected the abundance of Cassin’s auklets and krill biomass within the Sanctuaries. Our models indicated that predicted krill abundance was associated with waters upwelled 3 months prior to the prediction month. Previous studies on euphausiids within the CCS have found significant associations between krill and UI at different lags (30-, 60-, or 365-day), revealing the dependence of adult krill survival and development upon both short- and long-term ocean climate conditions [31,67]. Additionally, Cassin’s auklet abundance was predicted to be higher with greater sea surface salinity, signifying the importance of upwelled waters.

Cassin’s auklets are sensitive to changes in prey availability which is directly driven by changes in climate [74]. Previous modeling efforts [38] suggested Cassin’s auklets are more strongly influenced by larger-scale ocean forcing compared with local conditions. Conversely, our models show krill biomass to be more strongly influenced by local conditions than basin-scale climate. This finding suggests that predators and prey respond to ocean conditions at differing spatial scales likely associated with disparate individual ranges. Warm El Niño years have been intensifying [4] and a major reduction in zooplankton biomass has been recorded [75]. Recent research has shown how the Cassin’s auklet’s response to different climate variables has changed over time [14,15,67,76]; this may be due, in part, to the increase in frequency of a different type of Pacific warming event (termed El Niño modoki), which is caused by a disruption in the normal atmospheric circulation [77]. An increase in frequency and magnitude of these unproductive conditions will negatively impact marine species on multiple trophic levels [2,3].

While Cassin’s auklet diets from 1985 through 2013 were generally comprised of krill, non-euphausiid food items were present during anomalous years. For example, when combining diet results from 2005–2008, mysids made up >50% of chick diet composition, yet in all other years mysids comprised less than 50% of the chick diet. While mysids and decapods are the dominant shrimp-like crustaceans at depth (>500 m below the surface), krill are one of the most abundant crustaceans in the pelagic environment [78]. Additionally, although few diet samples were collected in 2005 and 2006, composition was almost entirely mysids. The lack of euphausiids in the diet during the unprecedented breeding failures of 2005 and 2006 points to the importance of krill to this species’ breeding success. From 2007 through 2013, the proportion of euphausiids in the auklet diet has increased and breeding success has been high.

Past studies on the Farallon Islands identified Euphausia pacifica and Thysanoessa spinifera as the main prey consumed by Cassin’s auklets, making up ~80% of the prey items identified [12,18,29,67,79]. When considering the krill portion of the auklet diet, these two species comprise ~96% of the krill observed in the diet in all years (Point Blue, unpublished data). Between 1985 and 2002, the cumulative adult stages of both euphausiid species comprised 82% of euphausiids eaten by Cassin’s auklets. Only in two years (1995 and 2014) did younger age classes (juvenile and immature stages) outnumber adults in the euphausiid portion of the diet, and in only two other years (2000 and 2009) did adult stages make up less than 70% of the euphausiid portion of Cassin’s auklet diet. In the remaining 23 years, adult euphausiids made up between 72% and 100% of euphausiids consumed by Cassin’s auklets. The larger mature stages of these species are the dominant age class consumed (adult lengths: 11.0–25.0 mm for E. pacifica, 16.0–25.0 mm for T. spinifera), while smaller juveniles and immature stages of both species are eaten to a lesser degree (juvenile/immature lengths: 5.0–11.0 mm for E. pacifica, 6.0–16.0 mm for T. spinifera; [78]).

A study on the Farallones revealed that Cassin’s auklets generally feed on less energy-rich E. pacifica during pre-breeding and the early part of the breeding season after which there may be a switch to more energy-rich T. spinifera [31,80]. Stomach samples from Cassin’s auklets on SEFI revealed different krill species were consumed during different upwelling scenarios; E. pacifica was consumed during times of upwelling relaxation (e.g., warmer surface water), while T. spinifera was the main prey during times of increased upwelling (e.g., cooler surface water [30]). Abraham and Sydeman [67] documented that Cassin’s auklet productivity correlated with chick growth rate, and growth rate was correlated with proportion of T. spinifera in the Cassin’s auklet diet between 1996–2001, suggesting the importance of this krill species to auklet reproductive success and chick development. A warming period in the CCS from the 1970s to the late 1990s, marked by more intense warm-water events, was associated with a decreased overall zooplankton biomass and a decline in many seabird populations [4,75,81].

Dive records from TDR-instrumented Cassin’s auklets revealed significant use of the upper water column. Burger and Powell [11] showed similar results, with auklets expending most of their foraging efforts in shallower portions of the upper water column (surface to 30 m) targeting euphausiid prey. Dives of different depths are a signature of alcid foraging behaviors associated with searching for prey [82,83,84]. We found adult krill at or near the surface in daytime towed hoop net samples across sampling years and Cassin’s auklet dive depths showing use of the upper water column by foraging, both interesting findings given later life stages of krill are known to exhibit diel vertical migration as a predator-avoidance mechanism [85].

Cassin’s auklets foraging for krill in the upper water column varied on a year-by-year basis. From 2010–2011, auklets were capturing abundant krill within the Sanctuaries as indicated by hoop net sampling; predicted krill density was also high. Conversely, in 2008, krill was less available for auklets as indicated by predicted densities and hoop net sampling. The story is less clear in 2009, when krill availability in hoop net samples and diet samples was high (just as in 2010–2011), but predicted krill density showed average levels.

Two years of particular interest were 2005 and 2010: 2005, when the first complete abandonment by Cassin’s auklets of the SEFI breeding colony was observed, and 2010, the year of the highest recorded breeding success to date. Low breeding success has been documented in strong El Niño years (e.g., 1983, 1992, and 1997), but in such years complete abandonment was not observed, as observed in 2005 and 2006 [31,67]. In 2005 and 2006, chick diet was highly anomalous: these are the only years exhibiting an absence of euphausiids in the chick diet; in the previous 19 years, euphausiids comprised at least 50% of items fed to the chick. In 2005, modeled krill abundance was low (Fig 5), though not extremely so, suggesting that the nature of the timing or spatial location of krill may have limited the ability of Cassin’s auklets to provide krill to their chick. Low krill abundance in turn was associated with strongly negative NPGO in that year. Other studies in the CCS have demonstrated the effect of a chronological mismatch between auklets’ breeding demands and plankton productivity [86,87].

In spring 2010, central CCS waters demonstrated a quick recovery to normal upwelling conditions following a mild El Niño event [16]; predicted Cassin’s auklet and krill biomass data reflected productive conditions. Krill biomass was higher still in 2011 than in 2010, yet reproductive success was lower in 2011 than in 2010. Euphausiid abundance in the southern CCS, off Baja California, though not directly spatially comparable, was relatively high towards summer and fall 2010, possibly a consequence of a productive cool-water La Niña event [16]. Again, 2011 was a stronger La Niña year than 2010, yet, as noted, reproductive success, while good that year, was not as high as 2010. The reproductive strategy of Cassin’s auklets, and attendant success, reflects the complexities of timing and location of prey availability as well as the condition of parents during the course of the breeding season [88].

The broader implications of this study are to highlight the value of hotspots along shelf breaks for krill predators in all upwelling systems, and that daytime feeding on diel migrating krill may still play a major role in the diet of these predators. In addition, although persistent krill habitat was identified, characteristics of predator prey relationships will likely vary with local and wide scale oceanographic conditions. For example, a shift was documented in trophic level patterns when seabird density in the southern California CCS changed, corresponding to the variability in krill and forage fish abundance [89]. Elements like the optimal foraging theory and predator response to prey shift and availability should be considered when managers are preparing for the impacts of climate changes across varying geographic scales.

Our study highlights persistent areas of high use by foraging Cassin’s auklets around Cordell Bank and the Farallones Escarpment, areas in which krill consistently occur. Changes in the local oceanography and basin-scale climate can lead to weakened or delayed upwelling and warmer ocean temperatures, thus decreasing prey availability. Diet samples confirmed a largely euphausiid-dominated diet among Cassin’s auklets that foraged in the upper water column. Despite diurnal migrations, some adult and late juvenile krill remain in the upper water column during daylight hours, when and where they are preyed upon by Cassin’s auklets. A projected increase in frequency and magnitude of anomalous ocean conditions may affect the long-term stability of the SEFI Cassin’s auklet population, as krill becomes less available, causing reduced breeding success and adult survival [13]. More broadly, krill predators in other systems may face similar impacts with oceanographic changes.

Supporting Information

S1 Table. Coefficients, standard errors, z values, p values, and associated 95% confidence intervals for all quantitative variables (including year and interactions with year) for the two-part model combining logistic regression and negative binomial regression for krill.

Model output assumes base year is 2004.

https://doi.org/10.1371/journal.pone.0144232.s001

(DOCX)

S2 Table. Coefficients, standard errors, z values, p values, and associated 95% confidence intervals for all quantitative variables (including years and interactions with years) for the zero-inflated negative binomial regression model for Cassin’s auklets.

Model output assumes base year is 2004.

https://doi.org/10.1371/journal.pone.0144232.s002

(DOCX)

Acknowledgments

We thank the many partners and contributors whose support has made this initiative and paper possible. This work was conducted by the Applied California Current Ecosystem Studies (ACCESS, www.accessoceans.org) Partnership, an ongoing collaboration between Point Blue Conservation Science and the Cordell Bank and Gulf of the Farallones National Marine Sanctuaries to support marine wildlife conservation and healthy ecosystems in northern and central California. We also thank Maria Brown and Dan Howard for their support, as well as Jan Roletto, Danielle Lipski, Kaitlin Graiff and Lisa Etherington and volunteers from Cordell Bank and the Gulf of the Farallones National Marine Sanctuaries for coordinating and facilitating the ACCESS cruises. Thanks to Pomona College students Zachary Brown, Eleanor Caves, Kristina McOmber, Nicole McDuffie, and professor Andre Cavalcanti (Pomona College) for help with TDR field work and/or analyses. Thanks to Danielle Lipski and Cotton Rockwood for comments on various stages of the draft which were helpful in improving the manuscript. Additional thanks to Dr. Alan Burger, Dr. Elliot Hazen, and an anonymous reviewer for their constructive reviews of this manuscript. This is Point Blue Conservation Science contribution number 2022.

Author Contributions

Conceived and designed the experiments: SM NN JJ. Performed the experiments: SM ME RB JH NK BS PW NN JJ. Analyzed the data: SM ME JH NK BS AS NN JJ. Contributed reagents/materials/analysis tools: SM ME JH NK JJ. Wrote the paper: SM.

References

  1. 1. Doney SC, Ruckelshaus M, Duffy JE, Barry JP, Chan F, English CA, et al. Climate Change Impacts on Marine Ecosystems. Ann Rev Mar Sci. 2012; 4: 11–37. pmid:22457967
  2. 2. Rhein M, Rintoul SR, Aoki S, Campos E, Chambers D, Feely RA, et al. Observations: Ocean. In: Stocker TF, Qin D, Plattner GK, Tignor M, Allen SK, Boschung J, et al. editors. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. New York: Cambridge University Press; 2013. pp. 255–316.
  3. 3. OEHHA. Indicators of Climate Change in California. Office of Environmental Health Hazard Assessment, California Environmental Protection Agency; 2013. Available: http://www.oehha.ca.gov/multimedia/epic/2013EnvIndicatorReport.html
  4. 4. McGowan JA, Cayan DR, Dorman LM. Climate-Ocean Variability and Ecosystem Response in the Northeast Pacific. Science. 1998; 281: 210–217. pmid:9660743
  5. 5. Enfield DB, Mestaz-Nunez AM. Multiscale variability in global sea surface temperatures and their relationships with tropospheric climate patterns. J Clim. 1999; 12: 2719–2733.
  6. 6. Sagarin RD, Barry JP, Gilman SE, Baxter CH. Climate-related change in an intertidal community over short and long time scales. Ecol Monogr. 1999; 69: 465–490.
  7. 7. Roemmich D, Church J, Gilson J, Monselesan D, Sutton P, Wijffels S. Unabated planetary warming and its ocean structure since 2006. Nat Clim Chang. 2015; 5: 240–245.
  8. 8. Palacios DM, Bograd SJ, Mendelssohn R, Schwing FB. Long-term and seasonal trends in stratification in the California Current, 1950–1993. J Geophys Res. 2004; 109: C10016.
  9. 9. Mendelssohn R, Schwing FB. Common and uncommon trends in SST and wind stress in the California and Peru—Chile current systems. Prog Oceanogr. 2002; 53: 141–162.
  10. 10. García-Reyes M, Largier J. Observations of increased wind-driven coastal upwelling off Central California. J Geophys Res. 2010; 115: C04011.
  11. 11. Burger AE, Powell DW. Diving depths and diet of Cassin’s auklet at Reef Island, British Columbia. Can J Zool. 1990; 68: 1572–1577.
  12. 12. Sydeman WJ, Hester MM, Thayer JA, Gress F, Martin P, Buffa J. Climate change, reproductive performance and diet composition of marine birds in the southern California Current system, 1969–1997. Prog Oceanogr. 2001; 49: 309–329.
  13. 13. Lee DE, Nur N, Sydeman WJ. Climate and demography of the planktivorous Cassin’s auklet Ptychoramphus aleuticus off northern California: implications for population change. J Anim Ecol. 2007; 76: 337–347. pmid:17302841
  14. 14. Jahncke J, Saenz BL, Abraham CL, Rintoul C, Bradley RW, Sydeman WJ. Ecosystem responses to short-term climate variability in the Gulf of the Farallones, California. Prog Oceanogr. 2008; 77: 182–193.
  15. 15. Sydeman WJ, Bradley RW, Warzybok P, Abraham CL, Jahncke J, Hyrenbach KD, et al. Planktivorous auklet Ptychoramphus aleuticus responses to ocean climate, 2005: Unusual atmospheric blocking? Geophys Res Lett. 2006; 33: L22S09.
  16. 16. Bjorkstedt EP, Gomez-Valdes J, Gaxiola-Castro G, Lavaniegos B, Warzybok P, Bradley R, et al. State of the California Current 2010–2011: Regionally variable responses to a strong (but fleeting?) La Nina. Calif Coop Ocean Fish Investig. 2011; 52: 36–68.
  17. 17. Ribic CA, Ainley DG, Spear LB. Effects of El Niño and La Niña on seabird assemblages in the Equatorial Pacific. Mar Ecol Prog Ser. 1992; 80: 109–124.
  18. 18. Sydeman WJ, Hobson KA, Pyle P, McLaren EB. Trophic relationships among seabirds in central California: Combined stable isotope and conventional dietary approach. Condor. 1997; 99: 327–336.
  19. 19. Hunt GL, Russell RW, Coyle KO, Weingartner T. Comparative foraging ecology of planktivorous auklets in relation to ocean physics and prey availability. Mar Ecol Prog Ser. 1998; 167: 241–259.
  20. 20. Ainley DG, Hyrenbach KD. Top-down and bottom-up factors affecting seabird population trends in the California current system (1985–2006). Prog Oceanogr. 2010; 84: 242–254.
  21. 21. Hickey BM. The California Current System—hypotheses and facts. Prog Oceanogr. 1979; 8: 191–279.
  22. 22. Chavez FP. Forcing and biological impact of onset of the 1992 El Niño in central California. Geophys Res Lett. 1996; 23: 265–268.
  23. 23. Tynan CT, Ainley DG, Barth JA, Cowles TJ, Pierce SD, Spear LB. Cetacean distributions relative to ocean processes in the northern California Current System. Deep Sea Res Part II. 2005; 52: 145–167.
  24. 24. Stevick PT, Incze LS, Kraus SD, Rosen S, Wolff N, Baukus A. Trophic relationships and oceanography on and around a small offshore bank. Mar Ecol Prog Ser. 2008; 363: 15–28.
  25. 25. Field JC, Francis RC, Aydin K. Top-down modeling and bottom-up dynamics: Linking a fisheries-based ecosystem model with climate hypotheses in the Northern California Current. Prog Oceanogr. 2006; 68: 238–270.
  26. 26. Dorman JG, Bollens SM, Slaughter AM. Population biology of euphausiids off northern California and effects of short time-scale wind events on Euphausia pacifica. Mar Ecol Prog Ser. 2005; 288: 183–198.
  27. 27. Ainley DG, Sydeman WJ, Norton J. Upper trophic level predators indicate interannual negative and positive anomalies in the California Current food web. Mar Ecol Prog Ser. 1995; 118: 69–79.
  28. 28. Hyrenbach KD, Veit RR. Ocean warming and seabird communities of the southern California Current System (1987–98): response at multiple temporal scales. Deep Sea Res Part II Top Stud Oceanogr. 2003; 50: 2537–2565.
  29. 29. Ainley DG, Boekelheide RJ. Seabirds of the Farallon Islands. Stanford, CA: Stanford University Press; 1990.
  30. 30. Ainley DG, Spear LB, Allen SG. Variation in the diet of Cassin’s auklet reveals spatial, seasonal, and decadal occurrence patterns of euphausiids off California, USA. Mar Ecol Prog Ser. 1996; 137: 1–10.
  31. 31. Abraham CL, Sydeman WJ. Prey-switching by Cassin’ s auklet Ptychoramphus aleuticus reveals seasonal climate-related cycles of Euphausia pacifica and Thysanoessa spinifera. Mar Ecol Prog Ser. 2006; 313: 271–283.
  32. 32. Simard Y, Mackas DL. Mesoscale aggregations of euphausiid sound scattering layers on the continental shelf of Vancouver Island. Can J Fish Aquat Sci. 1989; 46: 1238–1249.
  33. 33. Croll DA, Tershy BR, Hewitt RP, Demer DA, Fiedler PC, Smith SE, et al. An integrated approach to the foraging ecology of marine birds and mammals. Deep Res II. 1998; 45: 1353–1371.
  34. 34. Allen SE, Vindeirinho C, Thomson RE, Foreman MGG, Mackas DL. Physical and biological processes over a submarine canyon during an upwelling event. Can J Fish Aquat Sci. 2001; 58: 671–684.
  35. 35. Smith SE, Adams PB. Daytime surface swarms of Thysanoessa spinifera (Euphausiacea) in the Gulf of the Farallones, California. Bull Mar Sci. 1988; 42: 76–84.
  36. 36. Nur N, Jahncke J, Herzog MP, Howar J, Hyrenbach KD, Zamon JE, et al. Where the wild things are: predicting hotspots of seabird aggregations in the California Current System. Ecol Appl. 2011; 21: 2241–2257. pmid:21939058
  37. 37. Santora JA, Sydeman WJ, Schroeder ID, Wells BK, Field JC. Mesoscale structure and oceanographic determinants of krill hotspots in the California Current: Implications for trophic transfer and conservation. Prog Oceanogr. Elsevier Ltd; 2011; 91: 397–409.
  38. 38. McGowan J, Hines E, Elliott M, Howar J, Dransfield A, Nur N, et al. Using seabird habitat modeling to inform marine spatial planning in central California’s National Marine Sanctuaries. PLoS One. 2013; 8: e71406. pmid:23967206
  39. 39. Hunt . Foraging ecolgy as related to the distribution of planktivorous auklets in the Bering Sea. In: Vermeer K, Briggs KT, Morgan KH, Siegel-Causey D, editors. The status, ecology, and conservation of marine birds of the North Pacific. Ottawa: Canadian Wildlife Service Publications; 1993. pp. 18–26.
  40. 40. Ainley DG, Spear LB, Tynan CT, Barth JA, Pierce SD, Ford RG, Cowles TJ. Physical and biological variables affecting seabird distributions during the upwelling season of the northern California Current. Deep Sea Res Part II. 2005; 52: 123–143.
  41. 41. Genin A. Bio-physical coupling in the formation of zooplankton and fish aggregations over abrupt topographies. J Mar Syst. 2004; 50: 3–20.
  42. 42. Yen PPW, Sydeman WJ, Hyrenbach KD. Marine bird and cetacean associations with bathymetric habitats and shallow-water topographies: implications for trophic transfer and conservation. J Mar Syst. 2004; 50: 79–99.
  43. 43. Briggs KT, Tyler WMB, Lewis DB, Carlson DR. Bird communities at sea off California: 1975 to 1983. Studies in Avian Biology Number 11; Cooper Ornithological Society; 1987.
  44. 44. Briggs KT, Ainley DG, Spear LB, Adams PB, Smith SE. Distribution and diet of Cassin’s auklet and common murre in relation to central California upwellings. Proc Int Ornithol Congr. 1988; 19: 983–990.
  45. 45. Dransfield A, Hines E, McGowan J, Holzman B, Nur N, Elliott M, et al. Where the whales are: using habitat modeling to support changes in shipping regulations within National Marine Sanctuaries in Central California. Endanger Species Res. 2014; 26: 39–57.
  46. 46. Michael PE, Jahncke J, Hyrenbach KD. Relative influence of static and dynamic features on black-footed albatross (Phoebastria nigripes) habitat use in central California Sanctuaries. Fish Oceanogr. 2014; 23: 18–31.
  47. 47. Johnston K, Ver Hoef JM, Krivoruchko K, Lucas N. Using ArcGIS Geostatistical Analyst. Redlands, CA; 2001.
  48. 48. Barlow J, Ferguson MC, Becker EA, Redfern J V, Forney KA, Vilchis IL, et al. Predictive modeling of cetacean densities in the Eastern Pacific Ocean. SERDP Project S1-1391. National Oceanic and Atmosphereic Administration, Washington, DC; 2009.
  49. 49. Forney KA, Ferguson MC, Becker EA, Fiedler PC, Redfern J V, Barlow J, et al. Habitat-based spatial models of cetacean density in the eastern Pacific Ocean. Endanger Species Res. 2012; 16: 113–133.
  50. 50. Mantua NJ, Hare SR. The Pacific Decadal Oscillation. J Oceanogr. 2002; 58: 35–44.
  51. 51. Di Lorenzo E, Schneider N, Cobb KM, Franks PJS, Chhak K, Miller AJ, et al. North Pacific Gyre Oscillation links ocean climate and ecosystem change. Geophys Res Lett. 2008; 35: L08607.
  52. 52. Trenberth KE. Some Effects of Finite Sample Size and Persistence on Meteorological Statistics. Part I: Autocorrelations. Mon Weather Rev. 1984; 112: 2359–2368.
  53. 53. Hewitt RP, Watkins J, Naganobu M, Sushin V, Brierley AS, Demer D, et al. Biomass of Antarctic krill in the Scotia Sea in January/February 2000 and its use in revising an estimate of precautionary yield. Deep Res Part II. 2004; 51: 1215–1236.
  54. 54. Watkins JL, Brierley AS. Verification of the acoustic techniques used to identify Antarctic krill. ICES J Mar Sci. 2002; 59: 1326–1336.
  55. 55. Demer DA, Conti SG. Reconciling theoretical versus empirical target strengths of krill: effects of phase variability on the distorted-wave Born approximation. ICES J Mar Sci. 2003; 60: 429–434.
  56. 56. Demer DA, Conti SG. Validation of the stochastic distorted-wave Born approximation model with broad bandwidth total target strength measurements of Antarctic krill. ICES J Mar Sci. 2003; 60: 625–635.
  57. 57. Mukai T, Iida K, Ando Yy, Mikami H, Maki Y, Matsukura R. Measurements of swimming angles, density, and sound speed of the krill Euphausia pacifica for target strength estimation. Ocean ’04 MTTS/IEEE Techno-Ocean '04. 2004; 4: 383–388.
  58. 58. Conti SG, Demer DA. Improved parameterization of the SDWBA for estimating krill target strength. ICES J Mar Sci. 2006; 63: 928–935.
  59. 59. Tasker ML, Jones PH, Dixon T, Blake BF. Counting seabirds at sea from ships: A review of methods employed and a suggestion for a standardized aproach. Auk. 1984; 101: 567–577. Available: papers2://publication/uuid/4FC94E90-EA06-435F-9240-99CF460693E6
  60. 60. Zuur AF, Ieno EN, Walker NJ, Saveliev AA, Smith GM. Mixed effects models and extensions in ecology with R. New York: Springer; 2009.
  61. 61. Barry SC, Welsh AH. Generalized additive modelling and zero inflated count data. Ecol Modell. 2002; 157: 179–188.
  62. 62. Elith J, Leathwick JR. Species Distribution Models: Ecological Explanation and Prediction Across Space and Time. Annu Rev Ecol Evol Syst. 2009; 40: 677–697.
  63. 63. Zuur AF, Savaliev AA, Ieno EN. Zero inflated models and generalized linear mixed models with R. Newburgh, United Kington: Highland Statistics Ltd; 2013.
  64. 64. Martin TG, Wintle BA, Rhodes JR, Kuhnert PM, Field SA, Low-Choy SJ, et al. Zero tolerance ecology: improving ecological inference by modelling the source of zero observations. Ecol Lett. 2005; 8: 1235–1246. pmid:21352447
  65. 65. Long JS. Regression models for categorical and limited dependent variables. Thousand Oaks, CA: SAGE Publications, Inc.; 1997.
  66. 66. Speich S, Manuwal DA. Gular pouch development and population structure of Cassin’s auklet. Auk. 1974; 91: 291–306.
  67. 67. Abraham CL, Sydeman WJ. Ocean climate, euphausiids and auklet nesting: inter-annual trends and variation in phenology, diet and growth of a planktivorous seabird, Ptychoramphus aleuticus. Mar Ecol Prog Ser. 2004; 274: 235–250.
  68. 68. Karnovsky N, Harding A, Walkusz W, Kwasniewski S, Goszczko I, Wiktor J, et al. Foraging distributions of little auks Alle alle across the Greenland Sea: implications of present and future Arctic climate change. Mar Ecol Prog Ser. 2010; 415: 283–293.
  69. 69. Wilson RP, Culik BM, Peters G, Bannasch R. Diving behaviour of gentoo penguins, Pygoscelis papua; factors keeping dive profiles in shape. Mar Biol. 1996; 126: 153–162.
  70. 70. Mackas DL, Galbraith M. Zooplankton community composition along the inner portion of Line P during the 1997–1998 El Niño event. Prog Oceanogr. 2002; 54: 423–437.
  71. 71. Mackas DL, Tsurumi M, Galbraith MD, Yelland DR. Zooplankton distribution and dynamics in a North Pacific Eddy of coastal origin: II. Mechanisms of eddy colonization by and retention of offshore species. Deep Res Part II. 2005; 52: 1011–1035.
  72. 72. Brinton E. The distribution of Pacific euphausiids. Bull Scripps Inst Oceanogr. 1962; 8: 51–269.
  73. 73. Fontana RE. Physical process and their influence on zooplankton composition and foraging opportunities for seabirds in a strong, persistent upwelling region: Gulf of the Farallones-Cordell Bank (California). PhD dissertation, The University of California, Davis. 2013.
  74. 74. Wolf SG, Sydeman WJ, Hipfner JM, Abraham CL, Tershy BR, Croll DA. Range-wide reproductive consequences of ocean climate variability for the seabird Cassin’s auklet. Ecology. 2009; 90: 742–753. pmid:19341144
  75. 75. Roemmich D, McGowan J. Climatic warming and the decline of zooplankton in the California Current. Science. 1995; 267: 1324–1326. pmid:17812604
  76. 76. Wells BK, Field JC, Thayer JA, Grimes CB, Bograd SJ, Sydeman WJ, et al. Untangling the relationships among climate, prey and top predators in an ocean ecosystem. Mar Ecol Prog Ser. 2008; 364: 15–29.
  77. 77. Schmidt AE, Botsford LW, Eadie JM, Bradley RW, Di Lorenzo E, Jahncke J. Non-stationary seabird responses reveal shifting ENSO dynamics in the northeast Pacific. Mar Ecol Prog Ser. 2014; 499: 249–258.
  78. 78. Brinton E, Ohman M, Townsend A, Knight M, Bridgeman A. Euphausiids of the world ocean, Series: World Biodiversity Database CD-ROM Series Window version 1.0. Expert Center for Taxonomic Identification, Amsterdam, Netherlands. 2000.
  79. 79. Manuwal D, Thoreson AC. Cassin's auklet (Ptychoramphus aueuticus). In: Poole A, Gill F, editors. The Birds of North America, No. 50. Philadelpha: The Academy of Natural Sciences and Washington, D.C: The American Ornithologists' Union; 1993.
  80. 80. Ju S-J, Kang H-K, Kim WS, Harvey HR. Comparative lipid dynamics of euphausiids from the Antarctic and Northeast Pacific Oceans. Mar Biol. 2009; 156: 1459–1473.
  81. 81. Marinovic BB, Croll DA, Gong N, Benson SR, Chavez FP. Effects of the 1997–1999 El Niño and La Niña events on zooplankton abundance and euphausiid community composition within the Monterey Bay coastal upwelling system. Prog Oceanogr. 2002; 54: 265–277.
  82. 82. Falk K, Pedersen CE, Kampp K. Measurements of diving depth in dovekies (Alle alle). Auk. 2000; 117: 522–525.
  83. 83. Harding AMA, Egevang C, Walkusz W, Merkel F, Blanc S, Grémillet D. Estimating prey capture rates of a planktivorous seabird, the little auk (Alle alle), using diet, diving behaviour, and energy consumption. Polar Biol. 2009; 32: 785–796.
  84. 84. Karnovsky NJ, Brown ZW, Welcker J, Harding AMA, Walkusz W, Cavalcanti A, et al. Inter-colony comparison of diving behavior of an Arctic top predator: implications for warming in the Greenland Sea. Mar Ecol Prog Ser. 2011; 440: 229–240.
  85. 85. Bollens SM, Frost BW, Lin TS. Recruitment, growth, and diel vertical migration of Euphausia pacifica in a temperate fjord. Mar Biol. 1992; 114: 219–228.
  86. 86. Bertram DF, Mackas DL, McKinnell SM. The seasonal cycle revisited: interannual variation and ecosystem consequences. Prog Oceanogr. 2001; 49: 283–307.
  87. 87. Bertram DF, Harfenist A, Smith BD. Ocean climate and El Niño impacts on survival of Cassin’s Auklets from upwelling and downwelling domains of British Columbia. Can J Fish Aquat Sci. 2005; 62: 2841–2853.
  88. 88. Jacobs SR, Elliott KH, Gaston AJ. Parents are a drag: long-lived birds share the cost of increased foraging effort with their offspring, but males pass on more of the costs than females. PLoS ONE 2013; e54594. pmid:23382921
  89. 89. Sydeman WJ, Thompson SA, Santora JA, Koslow JA, Goericke R, Ohman MD. Climate-ecosystem change off southern California: time-dependent seabird predator-prey numerical responses. Deep Sea Res Part II 2015; 158–170.