INTRODUCTION

Marine birds are known to be good indicators of environmental change (Piatt et al. 2007, Sydeman et al. 2017), yet rarely are the causes of long-term population variability tested directly. This is largely due to the prohibitively expensive, but indispensable, nature of systematic long-term ecological monitoring (Lindenmayer et al. 2012). Short-term environmental perturbations can produce annual population fluctuations around a fixed or changing carrying capacity (K) that are best described by inverse density dependence (Goyert et al. 2017). Explaining changes in K requires a quantifiable assessment of its relationships with climatic and habitat variation. We examined the relationships between these 2 factors and the carrying capacity of 5 Alaskan seabird species—Black-legged (Rissa tridactyla) and Red-legged kittiwakes (R. brevirostris), Common (Uria aalge) and Thick-billed murres (U. lomvia), and Tufted Puffins (Fratercula cirrhata)—along a 4-decade time series (1974–2014). These 5 species occupy an array of foraging niches (Byrd and Williams 1993, Gaston and Hipfner 2000, Ainley et al. 2002, Piatt and Kitaysky 2002, Hatch et al. 2009): They feed on a spectrum of diverse (e.g., murres) to narrowly specific prey (e.g., Red-legged Kittiwakes), exploited at the sea surface (e.g., both species of kittiwake) or at depth (e.g., murres and puffins), in offshore (e.g., murres and kittiwakes) and/or nearshore ranges (e.g., puffins). We assessed the effects of long-term variability in the quality of foraging habitat on the number of individuals supported by this habitat (i.e. carrying capacity or quasi-equilibrium abundance, K; Garton et al. 2016).

Our objective was to determine whether changes in habitat or climate explained why murres and kittiwakes have both increased in the Arctic, but have responded differently from each other to conditions in southern Alaska, USA (Goyert et al. 2017). We also aimed to describe why puffins have declined in the Gulf of Alaska, consistent with their neighboring populations in the North Pacific (which are under review for listing pursuant to the Endangered Species Act; Hodum 2013). Links between climate and the predictability or availability of prey have been shown to affect the productivity and population trends of seabirds breeding on the Pribilof Islands of Alaska (Byrd et al. 2008, Renner et al. 2012, 2014). We evaluated how the population dynamics of our study species responded to (1) climatic indices (e.g., Pacific Decadal Oscillation and Arctic Oscillation), and (2) habitat variability, as described by sea surface temperature (SST) and secondary productivity (e.g., zooplankton concentrations; Hermann et al. 2013, M. A. Smith personal communication). Our analysis is the first attempt to explain the population dynamics of marine birds relative to these climatic indices and habitat variability across the entire extent of Alaska and its 4 large marine ecosystems: Aleutian Islands, Arctic, Eastern Bering Sea, and Gulf of Alaska (Figure 1;  https://www.st.nmfs.noaa.gov/ecosystems/lme/index, M. A. Smith personal communication).

FIGURE 1.

Sites of monitored seabird colonies in Alaska, USA (yellow triangles), and their corresponding ecoregions. St. Lazaria Island is within the Gulf of Alaska, even though it is shown outside the scope of this large marine ecosystem (as delineated by  https://www.st.nmfs.noaa.gov/ecosystems/lme/index; M. A. Smith personal communication). White polygons represent the spatial extent of extracted covariates: sea surface temperature (rectangles) and euphausiids (circles).

We tested the hypotheses that population rates of change would be positively related to the density of euphausiids (i.e. krill) and negatively related to increasing SST associated with a changing climate (Irons et al. 2008, Descamps et al. 2017). Our 5 study species feed primarily on young, cold-water forage fish, particularly ammodytids (e.g., sand lance [Ammodytes hexapterus]), osmerids (e.g., capelin [Mallotus villosus]), and gadids (e.g., pollock [Gadus chalcogrammus]; Drummond 2016). Red-legged Kittiwakes additionally specialize on several species of myctophid (lanternfish). Many seabirds, Thick-billed Murres and Black-legged Kittiwakes in particular, also feed facultatively on krill (Gaston and Hipfner 2000, Williams et al. 2008, Ito et al. 2010, Drummond 2016). While long-term data on the distribution and abundance of forage fish remain limited in scope, SST and krill biomass serve as proxies of foraging habitat quality. Therefore, secondary productivity conventionally functions as an indirect mechanistic link between primary productivity and predator abundance, where cold, nutrient-rich upwelling water positively influences the availability of prey (e.g., zooplankton and/or forage fish) via bottom-up trophic forcing (Shealer 2001, Frank et al. 2006). We expected to find variability in each species' population response to habitat and climate as a result of disparities in their foraging ecology. Thus, we provide insight into the relative vulnerability of Alaskan seabirds to long-term, ecosystem-wide change.

METHODS

RESULTS

Across our 4-decade-long time series, SST increased significantly in the Aleutian Islands (y = 6.6 + 0.019t, calculated using a Theil-Sen estimator) by, on average, an annual increment (time t) of 0.019°C from the mean temperature of 6.6°C in 1974 to 7.3°C in 2014 (Figure 2); the Mann-Kendall test confirmed that this increase was significant (S = 248, P = 0.006), and the residuals of this regression were not autocorrelated. The mean temperature in the first 2 decades of the Aleutian Island time series (6.8°C, including 1994) was significantly different from that in the last 2 decades (7.2°C; t₃₈ = −2.9, P = 0.006). We found a significant positive correlation between SST from the ERSSTv4 and the derived CORE and CCCma models (Pearson's r = 0.39, P = 0.01), which supported our use of inferred zooplankton data from the latter source. This metric of euphausiid density at the sea surface decreased slightly in the Gulf of Alaska (y = 4.6 – 0.0072t; Mann-Kendall S = –174, P = 0.05; Figure 3); the residuals of this regression were not autocorrelated. Similarly, the mean euphausiid density in the first 2 decades in the Gulf of Alaska time series was significantly greater than that in the last 2 decades (t₃₆ = 2.0, P = 0.05).

FIGURE 2.

Annual April–October mean National Oceanic and Atmospheric Administration (NOAA) extended reconstructed sea surface temperature data (ERSST Version 4) across Alaska (AK), USA. The spatial extent covered 2° cells that overlapped seabird colonies in the Aleutian Islands (AI), Arctic (ARC), Eastern Bering Sea (EBS), and Gulf of Alaska (GOA) ecoregions. The regression line indicates significant increases in sea surface temperature (SST) across the 40-yr time series in the AI.

FIGURE 3.

Annual April–September mean euphausiid concentration at the sea surface across Alaska (AK), USA. The spatial extent covered cells within a 90-km radius of seabird colonies in the Aleutian Islands (AI), Arctic (ARC), Eastern Bering Sea (EBS), and Gulf of Akaska (GOA) ecoregions. Data were sourced from the Coordinated Ocean-Ice Reference Experiments model (1974–2004) and from the Canadian Centre for Climate Modelling coupled global climate model (2005–2014). The regression line indicates a decrease in euphausiid concentration across the 40-yr time series in the GOA.

Preliminary analyses of auto- and cross-correlation functions indicated that there were no strong time lags (P > 0.05 for cross-correlation coefficients) between population rates of change and environmental covariates, and that only the PDOa was autocorrelated, at a 1-yr lag (autocorrelation coefficient = 0.482, P = 0.002). Such temporal autocorrelation in the covariates, along with the structure of the model which allowed for density dependence and time lags, helped to accommodate nonindependence in the response variable (rates of change). This resulted in successful model diagnostics: The models showed no evidence of nonindependence, as diagnosed by the lack of significant autocorrelation in the top population dynamics models and in the residuals of the aforementioned covariates regressed against time (most notably those covariates that showed significant increases or decreases through time, e.g., SST in the Aleutian Islands and euphausiids in the Gulf of Alaska). Across the entire extent of Alaska, population rates of change were generally positively related to the AOIa (Common Murres), PDOa (Tufted Puffins, Red-legged Kittiwakes, and Thick-billed Murres), and euphausiids (Black-legged Kittiwakes and Tufted Puffins). They were negatively related to the NPIAss (Common and Thick-billed murres) and SST (Common Murres; Table 1, Appendix Figure 5). These relationships varied by species at the ecoregion level (Table 2).

TABLE 1.

Parameter effects of the top kittiwake, murre, and puffin population models across Alaska, USA. We evaluated density-independent and density-dependent Ricker (linear) and Gompertz (log-linear) models, using population reconstructions for Black-legged (BK) and Red-legged kittiwakes (RK), Common (CM) and Thick-billed murres (TM), and Tufted Puffins (TP). For each species (Spp), the relative importance of covariates is shown in the first row, weighted across the candidate set of 672 models. Below this, each row represents one of the most competitive models (difference from the top model in Akaike's information criterion corrected for small sample size [ΔAIC_c] < 2), with its Akaike weight (w_i), R² value, intercept (μ), and covariate effects. The covariates qualified as: (1) mutually exclusive density-dependent terms, where N_t is the population size at time t and negative effects suggest density dependence, while N_{t − 1} and N_{t − 2} indicate a time lag of dependence on population sizes up to 2 yr past; (2) climatic covariates, including the annual Arctic Oscillation (AOIa), spring North Pacific anomaly (NPIAss), and annual Pacific Decadal Oscillation (PDOa); (3) habitat covariates, including sea surface temperature (SST) and euphausiid concentration (Euph); and (4) linear and nonlinear (quadratic) effects of time (t) on the rates of change, which suggest an increasing (+) or decreasing (−) trend in carrying capacity (K) with year. For example, the top candidate model for BK contains only the terms ln(N_t) (negative Gompertz density-dependence), t, and t² (decreasing time trend), and the second most competitive model is the same but with a positive effect of euphausiids also included. The top null density-dependent models (with a constant time trend in K) are italicized; those with ΔAIC_c > 2 (i.e. significantly worse models) are listed below the dotted line. Models highlighted in bold font are depicted in Figure 4; refer to Supplementary Material Table S1 for candidate sets of models with cumulative Akaike weights < 0.95.

FIGURE 4.

Most competitive covariate models of population dynamics across monitored seabird colonies in Alaska, USA, for (A–B) kittiwakes, (C–D) murres, and (E) puffins. Reconstructed population sizes, N(t), are shown as black solid lines with 95% confidence intervals (gray shading). Solid red lines show model-estimated population sizes, N̂(t), and blue dashed lines indicate model-estimated carrying capacity, K̂(t). Note that 2 of the data points at either end of each time series were included in tests of time lags, and the last data point was included in rate of change calculations. The models shown here were the best of those that contained at least one environmental covariate (highlighted in bold font in Table 1), including: (A) euphausiids for Black-legged Kittiwakes, (B) PDOa for Red-legged Kittiwakes, (C) NPIAss and SST for Common Murres, (D) NPIAss for Thick-billed Murres, and (E) PDOa for Tufted Puffins.

TABLE 2.

Parameter effects of the top kittiwake, murre, and puffin population models across the Alaskan ecoregions (Reg) colonized by each species: Aleutian Islands (AI), Arctic (ARC), Eastern Bering Sea (EBS), and Gulf of Alaska (GOA). See Table 1 for more details; refer to Supplementary Material Table S2 for candidate sets of models with cumulative Akaike weights < 0.95.

Population models that assumed that the carrying capacity of Black-legged Kittiwakes varied on an annual basis with time, climate, and/or habitat performed significantly better than the null model (density dependence without any covariates, i.e. a constant trend in K; Tables 1 and 2). The time trend in K captured the most variation in their population dynamics, particularly in the Arctic, Eastern Bering Sea, and Gulf of Alaska ecoregions. K increased in the Arctic and decreased elsewhere. Population rates of change were negatively related to euphausiids in the Eastern Bering Sea and positively related elsewhere. In the Aleutian Islands, rates of change were positively related to the PDOa and negatively related to the AOIa and NPIAss.

Using covariates from the density-dependent models, including time where appropriate, explained 114% (Tufted Puffins) to 483% (Common Murres) more of the variation in growth rates for all of the species except Red-legged Kittiwakes, for which it increased the variation explained by only 13%, a negligible amount (see R² values in Table 1 for these calculations, compared between the null model and best covariate model). For Red-legged Kittiwakes across their range in Alaska, the most descriptive models estimated either a constant K or a varying K that changed with time or the PDOa. A dynamic K was significantly better than the null model in both ecoregions occupied by this species (Aleutian Islands, Eastern Bering Sea). K decreased with time, euphausiids, and NPIAss in the Eastern Bering Sea. Increasingly positive rates of change and higher carrying capacities were associated with higher abundance of euphausiids in the Aleutian Islands. However, rates of change in this latter ecoregion were negatively related to the AOIa, NPIAss, and SST.

For Common Murres across Alaska, K varied with time, but SST, the NPIAss, and/or the AOIa described significantly more variation than the null model (constant K; except in the Eastern Bering Sea ecoregion). K decreased linearly with time in the Eastern Bering Sea, but increased elsewhere. Rates of change were inversely related to changes in SST in the Arctic and Gulf of Alaska ecoregions, and positively related to the PDOa in the Gulf of Alaska.

In the density-dependent models for Thick-billed Murres across Alaska, climate covariates were most important, where they (NPIAss, PDOa) predicted significantly more variation than the time trend alone (i.e. those covariates likely explained the perceived time trend). A dynamic K was significantly better than the null model for Thick-billed Murres in the Aleutian Islands, Eastern Bering Sea, and Gulf of Alaska ecoregions. Higher population rates of change were associated with higher estimates of euphausiids and the PDOa in the Gulf of Alaska and Aleutian Islands, and with lower values of SST in the Gulf of Alaska and Arctic. Rates of change in the Aleutian Islands and Arctic ecoregions were positively related to the NPIAss, but in the Eastern Bering Sea and Gulf of Alaska were negatively related to this index. They showed an inverse relationship with the AOIa in the Arctic and Eastern Bering Sea.

For Tufted Puffins, the most descriptive models estimated either a constant K or a K that varied linearly with time, the PDOa, or euphausiids. A dynamic K was significantly better than the null model in the Gulf of Alaska. K increased with time and the PDOa in the Aleutian Islands and decreased with time and the PDOa in the Gulf of Alaska.

In summary, a carrying capacity that varied annually with climatic and/or environmental covariates was most descriptive of the population dynamics of all 5 species (Figure 4). The order of importance of each covariate, aggregated across all 5 species and 4 ecoregions, suggested that the Gompertz model with a time trend in K was most predictive (i.e. retained the most weight), followed by models with influences from euphausiids, the PDOa, AOIa, NPIAss, and SST. Similarly, across the entire state, the most important covariates were the NPIAss and PDOa, followed by euphausiids, SST, and the AOIa.

DISCUSSION

Consistent with our hypotheses that population rates of change would be negatively related to climatic changes in SST and positively related to densities of krill (euphausiids), rates of change generally declined in response to increasing SST, and decreased in response to decreasing secondary productivity (except in the Eastern Bering Sea ecoregion, where carrying capacities were inversely related to krill). Slight decreases in inferred krill density occurred in the Gulf of Alaska, which may explain declines of Black-legged Kittiwakes in this region. Along with the PDO, krill was positively related to rates of population change of Tufted Puffins, but did not appear to contribute to puffin declines in the Gulf of Alaska. Inconsistent with our hypothesis, negative relationships between krill and kittiwake rates of change occurred in the Eastern Bering Sea. The Aleutian Island chain was the only region to show a significant increase in SST and to contribute to negative effects in Red-legged Kittiwake populations. Thus, both species of kittiwake demonstrated the greatest vulnerability to significant increases in SST (Aleutian Islands) and decreases in secondary productivity (Gulf of Alaska). Warming waters over the past 4 decades seem likely to have affected the predictability or availability of their prey base.

Inverse relationships between Common Murre population rates of change and SST occurred across Alaska, particularly in 2 of the regions (Gulf of Alaska and Arctic) where their populations grew. These relationships are consistent with other findings on the circumpolar dynamics of Common Murre populations (Irons et al. 2008). Thick-billed Murres demonstrated similar relationships (Byrd et al. 2008, Renner et al. 2012, 2014), but showed higher sensitivity to the NPIAss. Their population rates of change were negatively related to the NPIAss (analogous to the strength of the Aleutian Low Index), consistent with published effects on the productivity of other alcids (auklets [Aethia spp.]; Bond et al. 2011). Given that the carrying capacity of murres has increased over the last 4 decades, despite their inverse relationship with ocean warming, these results suggest that they may be more resilient to the observed environmental changes than the other species (Goyert et al. 2017). The high, positive rates of change that we documented in Common Murres may have represented a steady increase in K, despite drops in their population sizes as a response to environmental perturbation (Southwood 1981). This ability of murres to adapt to a volatile environment is supported by evidence of stable breeding success across several decades (Dragoo et al. 2015, Sydeman et al. 2017).

Compared with murres, the increased sensitivity of kittiwakes to degraded habitat in southern Alaska (e.g., from reduced secondary productivity and warmer waters) likely can be attributed to differences in their foraging ecology. Black-legged Kittiwake rates of population change responded positively to habitat (e.g., euphausiids) and negatively to climate (e.g., the PDO) in the Arctic and Gulf of Alaska regions. In other words, their populations decreased in the Gulf of Alaska, where secondary productivity deteriorated, suggesting that kittiwakes rely on habitat that supports healthy euphausiid populations. Meanwhile, the carrying capacity of murre populations increased in these regions and elsewhere. The ability of murres to adapt to the environmental variables examined in our study likely had to do, at least in part, with their foraging flexibility. Murres are generalists with flexible diets that exploit the entire water column, adjust their foraging time budgets, and maintain high reproductive success when food supplies shift (Burger and Piatt 1990, Harding et al. 2007). In contrast, both kittiwake species face many constraints as surface feeders with rigid diets and limited energy budgets (Furness and Tasker 2000, Einoder 2009). Forage fish, in particular the young age classes depredated by seabirds, are highly sensitive to climatic cycles and shifts in environmental regimes due to their high demands for zooplankton (Batchelder et al. 2013). Unpredictable changes in the availability of forage fish can result in spatial or temporal (i.e. phenological) mismatches with seabird distributions (Burthe et al. 2012, Staudinger et al. 2013). The resulting prey limitation likely drove the population declines that we observed in kittiwakes. Our study suggests not only that kittiwakes are susceptible to degraded foraging habitat and climate change, but also that they are useful as marine indicators (Sydeman et al. 2017).

Given the vast spatial and temporal extent of our study, the data layers available from PMEL on the inferred historical distribution of zooplankton in Alaska were a valuable and unique source of information, despite their limitations. The marine ecosystem in Alaska is highly dynamic, and using modeled data as covariates adds a degree of uncertainty to our conclusions. Furthermore, although we minimized collinearity (i.e. nonindependence) among our model covariates, we acknowledge that climatic indices are defined by several interacting factors in the environment. Long-term studies that collect raw data on oceanographic and climatological factors throughout the North Pacific are crucial to avoid further uncertainty in predicting how wildlife populations will respond to climate change (Lindenmayer et al. 2012). Our attempt to explain the population dynamics of seabirds with respect to habitat and climatic variability provides preliminary evidence of effects of multidecadal processes. However, further analysis is required to document the influence of relationships (e.g., teleconnections) between the regional and global processes driving changes in the North Pacific Ocean (Di Lorenzo and Mantua 2016). This would help to identify the nature of the trophic forcing that may serve as the link between warm waters, low availability of prey, and low breeding success in seabirds, particularly as these factors influence declines in kittiwake populations (Frank et al. 2006, Sydeman et al. 2017).

By including environmental factors in our models, we were able to estimate and explain population change in association with variation in climate and habitat quality. Further studies should examine the demographic factors that influence population change in seabirds. For example, kittiwake breeding success has been shown to vary with SST (Frederiksen et al. 2007) and prey availability (Robert et al. 2002, Frederiksen et al. 2005). Their survival, on the other hand, has been shown to be relatively constant, and likely shapes the long-term population trends of these long-lived species (as in other larids; Monaghan et al. 1989, Oro et al. 1999). While moderate changes in food supply may influence breeding success, only severe changes in prey availability are likely to affect survival, by failing to replenish body reserves (Cairns 1987, Oro and Furness 2002). Such extremely poor conditions can influence the mortality of seabirds outside the breeding season (Orben et al. 2015). Thus, breeding success most likely contributes to the annual population dynamics of our 5 study species, whereas survival probably drives their long-term population trends.

Our results provide evidence of strong bottom-up forcing in Alaskan ecosystems, likely mediated by high primary and secondary productivity during negative (i.e. cold) PDO regimes (Pinchuk et al. 2008) and low productivity during warming events associated with anomalies in the Aleution Low Index and NPI (Springer 1998). The extensive spatiotemporal coverage of our study restricted our analysis to examining bottom-up influences on our focal species. Top-down forcing also occurs, particularly at high latitudes, in low-productivity marine environments (Frank et al. 2006) or in unfavorably warm years (Hunt et al. 2002), but its source is generally specific to each colony (e.g., from predation by Larus spp., bald eagles [Haliaeetus leucocephalus], or Arctic fox [Vulpes lagopus] introductions; Regehr and Montevecchi 1997, Byrd et al. 2005). Commercial fisheries may impose some pressure on prey availability but, because they don't generally target forage fish, their impacts are more likely to be related to seabird bycatch (Žydelis et al. 2013). Competition with other top predators (e.g., marine mammals; Sinclair et al. 2008, Estes et al. 2009) could additionally affect populations of Alaskan seabirds. Kittiwakes may be more vulnerable than murres to the effects of competitors such as pink salmon (Oncorhynchus gorbuscha; Sydeman et al. 2017), particularly in warm PDO years (Mantua et al. 1997, Springer and van Vliet 2014). Future research should consider the potential of multispecies models to gauge the sensitivity of Alaskan seabird population dynamics to interspecific competition and/or facilitation.

This study supports the use of a biologically relevant scale of analysis (Scales et al. 2017) to identify environmental influences on wildlife population dynamics. Our ecoregion-level analyses were essential to revealing euphausiid concentrations as a significant predictor of rates of seabird population change. Extracting covariate data at a resolution consistent with the foraging range of our study species (e.g., ∼90 km; Smith et al. 2012) avoided diluting the effect of functional relationships. Had we analyzed these data only at the state-wide level, the climatic covariates would have emerged as the most important predictors because they did not vary among ecoregions. Understanding the indirect links between climate change and forage availability is crucial to identify the ecological drivers of population change and inform monitoring at regional scales. Isolating these factors is necessary to develop effective conservation policies, such as protections for Important Bird Areas (M. A. Smith personal communication), to enhance the success of population management and habitat mitigation under intensified environmental volatility.