Climate presses and pulses mediate the decline of a migratory predator

The climate in southern Argentina surrounding Punta Tombo is becoming warmer and wetter (34, 35). Regional average monthly precipitation extremes, defined as months with total precipitation greater than the 95th percentile, have increased by 41.6% in the last 30 y compared to the preceding 30 y, and average monthly maximum air temperatures have increased by 2.4% (Fig. 2C). In addition, mean sea surface temperature anomalies (SSTA) in the ocean surrounding the penguin’s foraging area during breeding warmed over the past 38 y ($\overline{\mathrm{x}}$ = −0.15 to 0.11) (Fig. 2C), which correlates with reduced chlorophyll-a and declining ocean productivity (36). Lastly, average plume strength has weakened over the past 38 y ($\overline{\mathrm{x}}$ = −0.31 to 0.29; positive values correspond to weaker plumes) (Fig. 2C), consistent with climate-driven changes in regional and local wind conditions (26, 28, 29). The increasing presses of climate and ocean variables drive more frequent and higher magnitude pulses of these variables (Figs. 1 and 2). Higher means and SDs in climate presses resulted in greater frequency and severity of climate pulses (SI Appendix, Fig. S1).

We collected survival data from 53,959 Magellanic penguins at Punta Tombo, as well as reproductive and population count data, between 1982 and 2019. To estimate the fecundity, juvenile and adult survival, and population size (number of breeding pairs) of the colony, we developed an IPM to fit the data. Adult survival has become more variable in the past 10 y than previously (ranging from 0.65 to 0.96 and SD from 2008 to 2019 = 0.13; SD from 1982 to 2007 = 0.051), whereas juvenile survival has become less variable during that time (ranging from 0.01 to 0.11 and SD from 2008 to 2019 = 0.037; SD from 1982 to 2007 = 0.13) (SI Appendix, Fig. S2), consistent with previous findings on survival trends in the colony (37). The penguin colony at Punta Tombo has fallen by approximately 50% in total size since the 1980s, with a population growth rate of 0.98 (95% CI: 0.977, 0.985) (SI Appendix, Fig. S2). Population growth was most sensitive to adult survival (r = 0.66, 95% CI: 0.434, 0.813; P < 0.001), followed by juvenile survival (r = 0.39, 95% CI: 0.079, 0.636; P = 0.01) and immigration (r = 0.33, 95% CI: 0.021, 0.686; P = 0.01), but not significantly sensitive to fecundity (r = −0.07, 95% CI: −0.383, 0.262; P = 0.69) (SI Appendix, Fig. S3). This indicates that strong fluctuations in adult survival and low juvenile survival may contribute to the population decline of penguins at the colony.

In order to estimate the effects of presses and pulses on penguin population persistence, we need to first calculate the relationship between climate variables and demographic rates. To estimate the general effect of climate variables on demographic rates, we used regressions in the IPM to relate climatic (rain and air temperature in the breeding grounds) and oceanographic (SSTA in the breeding grounds and plume strength in the wintering grounds) variables to fecundity, juvenile, and adult survival. Climate variables differentially impacted different life stages and demographic rates (Fig. 3). Rain negatively impacted fecundity ($\overline{\mathrm{x}}$ = −0.34, 95% credible intervals (CRI): −0.391, −0.286). Air temperature had a small negative effect on adult survival ($\overline{\mathrm{x}}$ = −0.09, 95% CRI: −0.129, −0.058), but had no influence on fecundity. SSTA in the breeding grounds negatively affected juvenile survival ($\overline{\mathrm{x}}$ = −0.40, 95% CRI: −0.440, −0.360) but had only very weak effects on fecundity and adult survival. The plume strength in the wintering grounds strongly influenced adult and juvenile survival, with higher juvenile and adult survival in years with weak plumes ($\overline{\mathrm{x}}$ = 0.61, 95% CRI: 0.555, 0.674; and $\overline{\mathrm{x}}$ = 0.29, 95% CRI: 0.222, 0.345; respectively). However, the prior winter’s plume did not impact penguin fecundity.

To determine how presses and pulses differentially impacted penguin population growth rates, we performed a sensitivity analysis of the IPM by varying the means and SDs (long-term presses) and frequency and magnitude of extreme events (pulses). We first varied means and SDs to examine effects if only presses, but not pulses, were considered. Under this scenario, mean plume strength had the greatest impact on population growth ($\overline{\mathrm{x}}$ = 0.96, 95% CI: 0.954, 0.963), followed by the shape and rate of rain ($\overline{\mathrm{x}}$ = −0.14, 95% CI: −0.202, −0.065; and $\overline{\mathrm{x}}$ = 0.38, 95% CI: 0.333, 0.443; respectively) and mean SSTA ($\overline{\mathrm{x}}$ = 0.20, 95% CI: 0.130, 0.269) (Fig. 4A). Importantly, however, incorporating pulses changed the results of our sensitivity analysis. When both presses and pulses were considered, the means of SSTA and air temperature (presses) no longer affected population growth (Fig. 4 A and B), whereas the frequency of air temperature extremes (pulse) had a strong negative effect ($\overline{\mathrm{x}}$ = −0.15, 95% CI: −0.225, −0.090) (Fig. 4B). In addition, years with higher frequency and magnitudes of weak plumes significantly improved population growth ($\overline{\mathrm{x}}$ = 0.24, 95% CI: 0.182, 0.299; and $\overline{\mathrm{x}}$ = 0.12, 95% CI: 0.057, 0.177; respectively). Therefore, presses and pulses, operating across multiple climate variables and life stages, differentially impacted penguin population growth.

Climatic and oceanographic presses and pulses differentially impacted penguin demographic rates—but what are the expected outcomes for the population given the trajectory of future climate change? Using the empirically parameterized IPM, we simulated penguin population growth over the next 50 y under future climate scenarios, including examining the sole effects of the strongest presses and pulses (Fig. 4). We explored four scenarios: 1) no change in climate variables from present; 2) projected climate change according to the current rate of change in climate variables (±SD); 3) change in the mean strength of the plume (press; ±SD); and 4) change in the frequency of air temperature extremes (pulse; ±SD) (Fig. 5). Under the current projected rate of climate change for all variables, the population growth rate will be higher (but still negative) than if the climate did not change from present ($\overline{\mathrm{x}}$ = 0.95, 95% CI: 0.951, 0.961; and $\overline{\mathrm{x}}$ = 0.93, 95% CI: 0.929, 0.937; respectively), due to the positive effects of the weakening plume on adult and juvenile survival. Confirming this finding, changes in solely the mean strength of the plume (press) had an even greater positive effect on population growth rate ($\overline{\mathrm{x}}$ = 0.97, 95% CI: 0.969, 0.975). In contrast, changes in solely the most impactful pulse, the frequency of air temperature extremes, drove growth rates even lower than other scenarios ($\overline{\mathrm{x}}$ = 0.92, 95% CI: 0.911, 0.928). In all scenarios, the population trajectory eventually approached extirpation (defined as population size near 0), given the population growth rate <1. In sum, opposing effects of presses and pulses of different climate variables lead to the mediation of penguin population decline.

Our work highlights the importance of disentangling how animals are differentially impacted by multiple, co-occurring climate factors across their life cycle. Yet, many studies focus on the effect of single climate factors on one component of demography across an animal’s life cycle (in contrast, see ref. 8, 38, and 39). Further, our work highlights the challenges faced by migratory species that are exposed to multiple presses and pulses throughout their migration cycle and across broad geographies. We showed that migratory Magellanic penguins experienced contrasting effects of climatic and oceanic presses and pulses on demography and population persistence that mediated the rate of population decline. For example, we found that both the presses of the Rio de la Plata plume and the pulses of air temperature extremes affected penguin population persistence in opposing ways. Such results underscore the complexities animal populations face under climate change, as well as both the challenges and opportunities for researchers in accurately predicting population responses to future climate scenarios.

Understanding which climate variables most affect population growth, and where these impacts occur geographically, not only aids prediction of population responses to climate change but can also inform conservation planning. For example, our results revealed that adult survival, followed by juvenile survival, contributes most to population growth (SI Appendix, Fig. S3) and that each of these processes was heavily affected by Rio de la Plata plume conditions which affect prey availability during energy-intensive foraging migrations (Fig. 3). These findings suggest the need to protect migratory foraging areas near the plume, where currently no habitat protection exists, to help protect the penguins’ forage base during migration and buffer against strong plume years (26). Such understanding may also inform prioritization of limited conservation resources. For example, although penguin fecundity is negatively impacted by rainstorms (Fig. 3), these processes do not strongly impact population persistence (Fig. 4) because the population growth rate is relatively insensitive to fecundity (SI Appendix, Fig. S3). In contrast, heat waves had a significant impact on adult survival, which strongly affected population growth, suggesting that conservation resources that mitigate the effects of heat waves would be relatively more cost-effective. Therefore, following our framework, conservation strategies can identify which environmental processes are demographically more important (e.g., heat waves and plume strength) and which may be less important (e.g., rainstorms) to better concentrate conservation efforts.

This research also highlights the relevance of disentangling the importance of climate presses versus pulses on population persistence. Understanding the relative impacts of presses versus pulses may help identify the timescales of conservation strategies (e.g., short term vs. long term) needed to promote persistence. For example, it may be necessary to mitigate short-term exposure of penguins to heat stress during heat waves, but also necessary for long-term management of ocean habitat or fisheries along migratory corridors (40). Furthermore, recognizing the impacts of increasing climate variability, manifested as both presses and pulses, on population persistence indicates the potential value of implementing dynamic management strategies. Dynamic management, whereby management is adjusted dynamically based on changing environmental, biological, or socioeconomic conditions, is increasingly called upon as a means of better coping with environmental variability (41). In the context of conserving migratory predators, proactive seasonal risk assessments based on projected climate conditions could be merged with seasonal management zoning to improve population persistence (42). For example, in seasons or years with predicted poorer ocean productivity and stronger plumes, regulations of fisheries and/or petroleum extraction or transport along penguin migratory corridors could be spatially or temporally adjusted to offset poor foraging conditions and enhance protection of penguins. Thus, identifying how presses and pulses across multiple climate variables differentially impact population growth may help target the types, timescales, and locations of conservation actions that will maximally improve population persistence in a changing climate.

The Punta Tombo colony of Magellanic penguins in Argentina (44°03’ S, 65°13’ W) was established in the 1920s and was once the largest colony of its species, reaching peak size of ~400,000 breeding pairs in the 1960s–1980s (20, 30, 33). Since then, the population has declined to ~150,000 breeding pairs as of 2019. Magellanic penguins in the south Atlantic are migratory from April to August, following the seasonal migration of Argentine anchovy (Engraulis anchoita) to northern Argentina up to southern Brazil (21). Breeding adults arrive at Punta Tombo in September or early October and lay two eggs in mid-October to early November. Chicks hatch in mid-November to early December, and fledge in mid-January to mid-February (22, 27).

To measure annual survival, we banded a total of 45,359 fledglings and 8,600 adults at the Punta Tombo colony from 1982 to 2019. We used stainless-steel bands marked with a unique identification number, which we fit to the left flipper of penguins, or we used small metal tags attached to the foot webbing (44). We note that there is little to no evidence that stainless-steel bands significantly affected the foraging, reproduction, or survival of penguins in our study (44, 45), so bands are unlikely to affect the relationship between climate variables and demographic parameters. We conducted daily searches for banded individuals throughout the breeding season (37). To measure reproduction, we followed individual nests at the same colony from 1983 to 2019 in an area of ~7,200 km² from mid-September to late February, and marked and monitored chicks. We counted chicks as fledged if they were last seen after January 10, weighed ≥ 1,800 g, and had not lost weight during the last few measurements (46). We calculated reproductive success in each year as the number of fledged chicks divided by the number of nests with eggs (maximum eggs per nest = 2).

To measure breeding population density, we established a permanent grid with lines running 1.5 km north to south and 2 km east to west, and counted nests, penguins, eggs, and chicks within a 100-m² circle every 100 m along the grid from 1987 to 2019 (33). We multiplied the estimated breeding pair density in each year by 3.524 km², which is the estimated occupied penguin habitat at Punta Tombo (33). Prior to 2016, many surveys were conducted in early October, when unpaired, reproductively unsuccessful males could have been incorrectly counted as part of a breeding pair. In 1993 (33), 2017, and 2019, we conducted replicate surveys in both early and late October to determine the difference in density estimates during early October when extraneous males could be counted, and late October, when breeding pairs are established in the colony and most pairs have eggs. The average difference between replicate surveys in early vs. late October was 1.36, which we used to adjust breeding pair density estimates prior to 2016. Afterward, population surveys were conducted in late October.

We built a pre-breeding, three-stage, female-only IPM with immigration for the penguin population following Kery & Schaub (47). Males were excluded from this model as the sex ratio at the colony at Punta Tombo is increasingly male-biased and males are not a limiting factor for reproduction or population growth (37, 48). Although Magellanic penguins do not breed until age 4+, and sometimes skip breeding due to poor body condition (33), we assumed a simple stage structure due to evidence that survival among adults of all ages is very similar (37). Additionally, models with additional stage structure (e.g., pre-breeding and non-breeding stages) exhibited similar dynamics to simpler models. Stages included juveniles (N_juv), adults (N_ad), and immigrants (N_imm). We modeled the number of individuals in these stages with binomial and Poisson distributions to account for demographic stochasticity:

where f_t is the fecundity (number of fledged chicks per nest with eggs), φ_juv,t is the apparent juvenile (from fledging to first year) survival probability, φ_ad,t is the apparent adult survival probability of individuals over one year of age, and ω_t is the immigration rate, which is a latent (unobserved) variable expressed as a per capita rate relative to the total population size the year before. f_t is multiplied by 0.5 to reflect an even sex ratio at hatching (22). We assumed immigration and emigration (incorporated within the apparent survival rates) due to evidence of metapopulation dynamics among Magellanic penguin colonies (49, 50). However, we believe that the emigration rate is very small relative to survival. Only 226 banded penguins at our study site have been found alive at other surveyed colonies, equating to a rough estimate of 0.04% emigration.

The observation equation links the observed population count (y) with the true population size (N), assuming Normal observation error:

where τ_y is the observation error, and the distribution is truncated to prevent population values falling below zero. The likelihood for the population count data is denoted as L_CO(y | ϕ_juv, ϕ_ad, ω, N, τ_y, f).

For capture–recapture data, we used the standard Cormark–Jolly–Seber model with a multinomial likelihood that allows for estimating adult and juvenile survival and recapture (p) probabilities. We did not conduct surveys or recapture during 2011, so we manually set the recapture probability, p, for this year to be approximately zero. The likelihood for the capture–recapture data is denoted as L_CR(m | ϕ_juv, ϕ_ad, p), where m is the capture–recapture data summarized in the m-array format. We modeled fecundity data (i.e., the number of fledglings produced) using a Poisson distribution. The likelihood for the fecundity data is denoted as L_PR(J, R | f), where R is the total number of surveyed nests and J is the total number of fledglings. We assumed independence among our population, capture–recapture, and fecundity datasets. Therefore, we constructed the joint likelihood as a product of the likelihoods for the datasets (y, m, J, R):

At Punta Tombo, we recorded daily precipitation and air temperature data. We recorded daily precipitation (±0.1 mm) using a manually emptied plastic rain gauge and daily maximum air temperature using a Sper Scientific min/max thermometer. We calculated the RAIN covariate as the total precipitation between October 15 and December 15, as this was found to have a strong impact on reproductive success (30, 31). We calculated the TEMP covariate as the % of days per breeding season when maximum air temperature was greater than 25 °C, which drove heat deaths (32). Because the locally collected covariates were not measured throughout the year, we also estimated the monthly rate of precipitation and air temperature change at Punta Tombo using TerraClimate, a global ~4 km resolution dataset of climate and climatic water balance (51).

We extracted SSTA for the foraging area during breeding season using the Coral Reef Temperature Anomaly Database (52), which contains global ~4 km resolution SSTA data on a weekly time scale. We extracted these data for an area bounded from 65.6 to 61° W, 45 to 40° S, which takes up >95% of maximum foraging distances for breeding penguins and is important habitat for adults feeding chicks (23, 24, 40). We calculated the SSTA covariate by spatially and temporally averaging the weekly SSTA data from October–February to align with the breeding season.

We also extracted SSTAs during migration for an area from 40 to 23° S and from the coast to the 1,000 m isobath. We excluded pixels representing water deeper than 1,000 m because Magellanic penguins generally migrate over the continental shelf (21, 26). We then conducted a principal component analysis of the SSTA data to quantify the spatial and temporal patterns in the sea surface temperature and extend the index developed by Rebstock et al. (26). The principal pattern of the PCA shows the extent and strength of the Rio de la Plata plume, with positive values indicating a weak, contracted river plume, and negative values indicating a strong, extended river plume. The first three principal components represented 30.4, 14.0, and 6.3% of the variance in the SSTAs. Following Rebstock et al. (26), we only used the first principal component, and temporally averaged the monthly plume index from May–August to align with the migration season.

We fitted our IPM using log-linear and logit-linear regressions to estimate the effect of climatic (rain, air temperature) and oceanographic (SSTA and plume index) variables on demographic parameters (fecundity, juvenile and adult survival). We modeled fecundity as a function of rain, air temperature, SSTA, and plume index with a log-link:

where µ_f is the intercept, α₁ is the slope for the climatic covariate RAIN_t, α₂ is the slope for the climatic covariate TEMP_t, α₃ is the slope for the oceanographic covariate SSTA_t, and α₄ is the slope for the oceanographic covariate PLUME_t prior to the breeding season.

We modeled juvenile and adult survival as a function of climatic and oceanographic variables with a logit-link:

where µ_juv and µ_ad are the intercepts, and PLUME_l_t is the oceanographic covariate PLUME_t, but lagged forward to account for the fact that the Rio de la Plata Plume will affect juvenile and adult survival during the migratory stage of their life cycle. SSTA_t was included to account for indirect effects on fledgling body condition and effects during the beginning stages of migration, where fledglings are near the breeding–foraging grounds of colony. Each covariate was standardized to have zero mean and unit variance to compare the effect sizes (slopes) of α to each other. We computed the mean and 95% CRI for the regression coefficients α.

We modeled the intercepts of the demographic parameters (µ_f, µ_juv, and µ_ad) and slopes of covariates (α) using vaguely informative prior distributions of N(0, 0.34). Observation error (${\tau }_y$) was given a vague prior of U(0, 4*10⁻¹⁰). Bayesian posterior distributions were approximated via Markov chain Monte Carlo (MCMC) algorithms. Following standard practice (47), three independent MCMC chains of 30,000 iterations were used with a burn-in period of 10,000. Final inferences were derived from a sample of 3 × 20,000 iterations that resulted from merging the three chains. Gelman–Rubin convergence diagnostics were below 1.1 for each parameter and chains were well mixed (53). We assessed goodness-of-fit graphically using posterior predictive checks, and calculated Bayesian p values, all of which were near 0.5, indicating that the model fit the data well (47). We performed the analyses using Nimble (54, 55) and program R (56).

We conducted sensitivity analyses to quantify the effects of the presses and pulses of climatic and oceanographic variables on the population persistence of Magellanic penguins, as measured by the geometric mean lambda, which is the product of population growth rates, and more useful for characterizing fluctuating populations than using the arithmetic mean lambda. We ran 1,000 simulations for 100 y each of penguin population trajectories using Nimble and the parameter values we estimated previously. In these simulations, we generated climatic and oceanographic data using normal and gamma distributions which were fit to the prior observed climate data. Normal distributions best fit the air temperature, SSTA, and plume index distributions, whereas gamma distributions best fit the precipitation data as it was positively skewed (Fig. 2C). We did not find any evidence of autocorrelation in the climatic and oceanographic variables based on calculating correlation coefficients at multiple time lags, so we did not add autocorrelation to the randomly drawn values. In each simulation, we randomly drew 100 values from these distributions to represent the yearly climate variables that penguins experienced. We ran sensitivity analyses by varying the mean and standard deviation (for normal distributions) and shape and rate (for gamma distributions) according to +/− the past rate of change in the past climate variables. We did this to generate climate variables within the realistic rate of climate change in the sensitivity analysis. For each simulation in the sensitivity analysis, we also recorded the frequency and magnitude of extreme events in the 100 randomly drawn climate values by calculating the frequency of values greater than the 95% quantile and the maximum value, respectively (8). We did not calculate the frequency of values less than the 5% quantile and the minimum value because all climatic and oceanographic variables were changing towards positive extreme events.

We then calculated partial correlation coefficients (PCC) via linear regression to determine the relationship between the climate parameters in our sensitivity analysis (mean, SD, frequency of extreme events, magnitude of extreme events) and geometric lambda. The sign of the PCC determines the qualitative relationship between the parameters and population growth, and the magnitude indicates the importance of climate variables in contributing to population growth (57). To explore the effects of including extreme events, we calculated the PCC with and without the frequency and magnitude of extreme events variables.

We then projected the population forward under four climate scenarios by running 100 simulations for 50 y. Scenario 1 (no change in climate variables from present) was run using generated climatic and oceanographic data from normal and gamma distributions based on the prior observed climate data. Scenario 2 (projected climate change based on current rate of change) was run using generated data from distributions based on the observed rate of increase of the climate data (±SD) between 1958 and 1988 and 1989 and 2019 for precipitation and air temperature data, and between 1982 and 2000 and 2001 and 2020 for oceanographic data. Scenario 3 (change in the mean strength, or press, of the plume) was run by using generated data for the mean plume from Scenario 2, and then generated data for all other variables from Scenario 1. Scenario 4 (change in the frequency, or pulse, of air temperature extremes) was run by using generated data for frequency of air temperature extremes based on the rate of increase in air temperature extremes, and then generated data for all other variables from Scenario 1. For all scenarios, we calculated the mean and 95% CI across the 100 simulations and calculated geometric lambda.

T.J.C.-W., P.D.B., G.A.R., and B.A. designed research; T.J.C.-W., P.D.B., G.A.R., and B.A. performed research; T.J.C.-W. analyzed data; and T.J.C.-W., P.D.B., G.A.R., and B.A. wrote the paper.

The authors declare no competing interest.