0% found this document useful (0 votes)
3 views52 pages

9584173

Download as pdf or txt
Download as pdf or txt
Download as pdf or txt
You are on page 1/ 52

Effects of climate change on an emperor penguin

population: analysis of coupled demographic and


climate models.
Stéphanie Jenouvrier1,2,4,5 ,
Marika Holland3 ,
Julienne Strœve4 ,
Christophe Barbraud2 ,
Henri Weimerskirch2 ,
Mark Serreze4,5 ,
Hal Caswell1,6
June 21, 2012

Running title: Emperor penguins and climate change

To whom correspondence should be addressed to sjenouvrier@whoi.edu.


phone:+1 508 289 3245

Emails of co-authors: mholland@ucar.edu; stroeve@kryos.colorado.edu;


barbraud@cebc.cnrs.fr; henriw@cebc.cnrs.fr; serreze@kryos.colorado.edu; hcaswell@whoi.edu.

1. Biology Dept., MS-34, Woods Hole Oceanographic Institution, Woods Hole, MA 02543,
USA

2. Centre d’Etudes Biologiques de Chizé, Centre National de la Recherche Scientifique,


F-79360 Villiers en Bois, France

3. Oceanography Section, National Center for Atmospheric Research, Boulder, 80305 CO,
USA

4. National Snow and Ice Data Center, Boulder, 80309 CO, USA

5. Cooperative Institute for Research in Environmental Science, University of Colorado


Boulder, 80309-0449, USA

6. Max Planck Institute for Demographic Research, Rostock, Germany.

Keywords: stochastic matrix population model; stochastic climate forecast; IPCC; un-
certainties; sea ice; seabirds.

1
1 The following table of contents is provided for convenience of reviewers.
2 Contents
3 1 Introduction 4

4 2 Emperor penguin life cycle and sea ice 6

5 3 Influence of sea ice on population growth 13

6 4 Population response to climate change 15

7 5 Discussion 17

8 6 Acknowledgments 26

9 7 Supporting Information legends 35

10 8 Tables 36

11 9 Figures 41

2
12 Abstract

13 Sea ice conditions in the Antarctic affect the life cycle of the emperor penguin

14 (Aptenodytes forsteri ). We present a population projection for the emperor pen-

15 guin population of Terre Adélie, Antarctica, by linking demographic models (stage-

16 structured, seasonal, nonlinear, two-sex matrix population models) to sea ice forecasts

17 from an ensemble of IPCC climate models. Based on maximum likelihood capture-

18 mark-recapture analysis, we find that seasonal sea ice concentration anomalies (SICa )

19 affect adult survival and breeding success. Demographic models show that both deter-

20 ministic and stochastic population growth rates are maximized at intermediate values

21 of annual SICa , because neither the complete absence of sea ice, nor heavy and per-

22 sistent sea ice, would provide satisfactory conditions for the emperor penguin. We

23 show that under some conditions the stochastic growth rate is positively affected by

24 the variance in SICa . We identify an ensemble of 5 general circulation climate models

25 whose output closely matches the historical record of sea ice concentration in Terre

26 Adélie. The output of this ensemble is used to produce stochastic forecasts of SICa ,

27 which in turn drive the population model. Uncertainty is included by incorporating

28 multiple climate models and by a parametric bootstrap procedure that includes pa-

29 rameter uncertainty due to both model selection and estimation error. The median of

30 these simulations predicts a decline of the Terre Adélie emperor penguin population

31 of 81% by the year 2100. We find a 43% chance of an even greater decline, of 90%

32 or more. The uncertainty in population projections reflects large differences among

33 climate models in their forecasts of future sea ice conditions. One such model predicts

34 population increases over much of the century, but overall, the ensemble of models

35 predicts that population declines are far more likely than population increases. We

36 conclude that climate change is a significant risk for the emperor penguin. Our ana-

37 lytical approach, in which demographic models are linked to IPCC climate models, is

38 powerful and generally applicable to other species and systems.

3
39 1 Introduction
40 Given the observed and predicted changes in climate (Solomon et al., 2007), conservation
41 biologists face a major challenge because climate affects all aspects of the life cycle of a species
42 (life history traits, phenology, movement). This results in changes in populations, species
43 distributions, and ecosystems (see reviews by Hughes (2000); McCarty (2001); Parmesan
44 (2006); Parmesan & Yohe (2003); Vitousek (1994); Walther et al. (2002)). The emperor
45 penguin (Aptenodytes forsteri ) is a species that is known to be extremely sensitive to climate
46 change, especially to changes in the sea ice environment (Ainley et al., 2010; Barbraud &
47 Weimerskirch, 2001; Croxall et al., 2002; Forcada & Trathan, 2009; Jenouvrier et al., 2005a,
48 2009b; Trathan et al., 2011; Fretwell & Trathan, 2009). In this paper, we analyze the
49 population responses of the emperor penguin to sea ice conditions and project its fate over
50 the rest of the 21st century using a novel and more comprehensive analysis than previous
51 studies.
52 Climate may affect vital rates (a general term for rates of survival, fertility, development,
53 etc.) in many ways (Ballerini et al., 2009; Stenseth et al., 2002) and possibly in opposite
54 directions (Barbraud & Weimerskirch, 2001). For example, mean winter temperature has a
55 positive effect on adult survival but a negative effect on fecundity of the Eurasian Oyster-
56 catcher (Haematopus ostralegus) (van de Pol et al., 2010). Therefore, to fully understand
57 the population effects of climate change, climate effects on the vital rates should ideally be
58 incorporated into models including the full life cycle (Adahl et al., 2006; van de Pol et al.,
59 2010). Although many studies have related climate to one or a few vital rates (e.g. survival
60 Jenouvrier et al. (2005b, 2003)), few have integrated these effects into the entire life cycle
61 (but see e.g. Jenouvrier et al. (2009b); Wolf et al. (2010); Hunter et al. (2010); Barbraud
62 et al. (2010)).
63 Seasonality can be particularly important, because the effects may occur at different times
64 during the seasonal cycle (Visser et al., 1998). Complex interactions of climate variables have
65 been shown to occur for seasonally breeding species, which may fail to adjust their breeding

4
66 phenology to track the peak of food availability for their young (Both et al., 2006; Moller
67 et al., 2008).
68 The responses of vital rates to climate conditions is only half of the story; to project
69 the population response to climate change, it is necessary to link the demographic models
70 to forecasts of future conditions (e.g. Jenouvrier et al. (2009b)). A primary source for such
71 forecasts is the set of climate model simulations that have contributed to the Intergovern-
72 mental Panel on Climate Change (IPCC) assessment reports (Stock et al., 2011). A growing
73 number of studies have now linked climate-dependent demographic models to these climate
74 models (e.g. seabirds Jenouvrier et al. (2009b); Barbraud et al. (2010); polar bears Hunter
75 et al. (2010)).
76 Our approach is to measure the effect of climate on the vital rates in a complete life
77 cycle model, to incorporate those vital rates into a population model to compute population
78 growth as a function of climate, and then to obtain forecasts of climate trajectories from
79 climate models and use those to drive projections of population growth. In this paper, we
80 extend previous studies (Barbraud and Weimerskirch 2001, Jenouvrier et al. 2005b, 2009b,
81 2010), especially that of Jenouvrier et al. (2009b), in several ways. (1) We obtain rigorous
82 statistical estimates of how sea ice, at different seasons of the year, affects penguin vital rates.
83 (2) We distinguish males and females, recognizing that the sexes differ in their sensitivity
84 to sea ice variations (Jenouvrier et al. 2005b) and that breeding is absolutely dependent on
85 participation by both males and females (Prevost, 1961). (3) We introduce a new method of
86 selecting climate models based on the agreement of their output with both the mean and the
87 variance in observed sea ice. (4) In order to include stochasticity in the climate forecasts,
88 only one or a few realizations of which are available, we developed a new method to estimate
89 stochasticity from time series of global circulation model (GCM) output.
90 Projections from models that are estimated from data are always accompanied by un-
91 certainty. If the results are to be useful for policy makers, it is essential to quantify that
92 uncertainty, and to draw conclusions that remain valid even given the range of uncertainty

5
93 (Hunter et al. 2010). We have used the results of statistical estimation of demographic pa-
94 rameters, and differences among an ensemble of climate models, to quantify this uncertainty.
95 The organization of the paper is as follows. In Section 2 we use a long-term dataset to
96 estimate the effects of sea ice on the vital rates. In Section 3 we evaluate the effect of sea ice
97 on deterministic and stochastic population growth rates for the emperor penguin. In Section
98 4, we compute stochastic sea ice forecasts from a selected set of climate models, and use
99 those forecasts to project population response to future sea ice change. We conclude with
100 discussion in Section 5.

101 2 Emperor penguin life cycle and sea ice


102 Various components of the sea ice environment affect different parts of the emperor
103 penguin life cycle during different seasons (see review by Croxall et al. (2002); Forcada &
104 Trathan (2009); Ainley et al. (2010)). In this section we outline the life cycle of the emperor
105 penguin and define the sea ice variables used in our analysis. Then we discuss the mechanisms
106 by which sea ice affects the vital rates and present the results of estimating these effects using
107 capture-mark-recapture (CMR) analysis.

108 Study population and data


109 Our analysis is based on a long-term data set on breeding emperor penguins at Dumont
110 D’Urville, Terre Adélie, in Antarctica (66◦ 40 S 140◦ 01 E). The colony has been monitored
111 every year, during the breeding season (March–December), from 1962 onwards. For details
112 of the history and data, see Prevost (1961) and Jenouvrier et al. (2005a). From 1962 on-
113 wards, breeding adults, number of eggs, frozen chicks, and surviving chicks at the end of the
114 breeding season have been counted, allowing the estimation of breeding success (Barbraud
115 & Weimerskirch, 2001). From 1968 to 1988, penguins were individually marked using flipper
116 bands. Banding stopped in 1988, but banded birds have been recorded since then. We limit
117 our analysis to the period before 2000 because too few marked birds were seen after that to
118 permit estimation.

6
119 The life cycle and the demographic model
120 Emperor penguins breed on motionless sea ice (i.e. fast ice) during the Antarctic winter.
121 They arrive at the colony sometime in late March to early April while sea ice is thickening,
122 and leave the colony in late December before the ice breaks up. The colony site is usually far
123 from the ocean, and during the breeding season emperor penguins travel long distances to
124 feed in ice-free areas, such as polynyas, within the sea ice cover. They feed on fish (mainly
125 Pleuragramma antarcticum), crustaceans (mainly Euphausia superba and amphipods), and
126 squid (mainly Psychroteuthis glacialis) (Cherel, 2008; Cherel & Kooyman, 1998; Kirkwood
127 & Robertson, 1997). The female lays a single egg in May, which is then incubated by the
128 male for two months while the female leaves the colony to feed. Females return when chicks
129 hatch around July, and both parents take turns feeding the chick until fledging in December.
130 Our demographic model is described in detail in Jenouvrier et al. (2010). It is a stage-
131 classified matrix population model with a projection interval of one year, but the annual
132 projection is based on four seasonal steps. The model has 5 stages: male and female pre-
133 breeders (birds that have yet to breed for the first time), breeding pairs, and male and
134 female non-breeders (birds that have bred before but do not do so in the current year). The
135 formation of pairs is a nonlinear function of the operational sex ratio.
136 Our model does not include density dependence because we expect small density ef-
137 fects in this population relative to effects of environmental variations, and especially sea ice
138 (Supplementary Appendix S1, Herrando-Pérez et al. (2012)).
139 The annual population projection depends on the vital rates: the probability that an
140 individual of a given stage returns to the breeding site, the probability of mating as a
141 function of the availability of potential mates, breeding success (probability that a breeding
142 pair raises offspring given the female lays an egg), the primary sex ratio (fixed at 0.5), the
143 survival of offspring during the first year at sea, and the annual survival of pre-breeders,
144 non-breeders and male and female breeders.
145 We divide the year into four seasons: (1) the non-breeding season from January to March,

7
146 (2) the arrival, copulation and laying period (April–May), hereafter called the laying period,
147 (3) the incubation period (June–July) and, (4) the rearing period (August–December).

148 Sea ice variables


149 Many components of the sea ice environment affect penguins in various way. To avoid
150 examining the effect of all possible factors on vital rates and to untangle the networks
151 of causation among them, we examine the covariation among several factors selected on
152 the basis on the emperor penguin responses to climate (i.e fast ice area and polynya area
153 indices, sea ice concentration (SIC) and sea ice extent (SIE), see Ainley et al. 2010 for a
154 comprehensive review). All these variables are strongly correlated (Supplementary Appendix
155 S2) and we focus our analysis on SIC, including the seasonality in SIC because it drives the
156 emperor penguin life cycle.
157 Sea ice concentration is the fraction of area covered by sea ice. Observed values of SIC
158 from 1979 to 2007 were obtained from passive microwave satellite imagery provided by the
159 National Snow and Ice data Center, using the NASA Team sea ice algorithm (see Cavalieri
160 et al. (1996) and http://nsidc.org/data). Forecasts of SIC from climate models were
161 extracted from 20 models available as part as the WCRP CMIP3 multi-model dataset from
162 1900 to 2100 (see Meehl et al. (2007) and http://esg.llnl.gov/portal).
163 In order to link population models to the output of GCMs, which use relatively coarse
164 spatial grids (100-200 km resolution), we use observed values of SIC over similarly large
165 spatial scales. We averaged SIC, both observed and simulated, over a large sector centered
166 on the colony. This sector included a 20 degree span in latitude, between longitudes 130◦ E
167 and 150◦ E during the breeding season, and between longitudes 120◦ E and 160◦ E during
168 the nonbreeding season. This includes the maximum foraging distances from the colony,
169 of about 100 km during the breeding season and at least 650 km during the non-breeding
170 season (Zimmer et al. 2008).
171 As a variable to describe the sea ice conditions, we use the proportional anomalies in SIC,
172 relative to the mean from 1979 to 2007. We denote this variable as SICa , and calculate it for

8
173 each of the four seasons (Fig. 1). We estimated the vital rates as functions of the observed
174 seasonal SICa (see following section), and used forecasts of seasonal SICa to project future
175 population trajectories (Section 4). However, to provide a comprehensive understanding
176 of the effect of seasonal SICa on vital rates and population growth rate, we present our
177 results as functions of ”annual sea ice” and ”seasonal difference in sea ice”, two variables that
178 accounted for most of the variability in the four seasonal SICa variables (Fig. 2, S2). Annual
179 sea ice is a weighted mean of seasonal SICa over the year. The seasonal difference in SICa is
180 the difference between ice concentration in the non-breeding season and its weighted average
181 over the combined incubation and rearing seasons. The contribution of SICa variations
182 during the laying season to the seasonal SICa difference variations is small. Thus positive
183 values of the seasonal SICa difference correspond to years with positive SIC anomalies in
184 the non-breeding season and negative SIC anomalies in the incubating and rearing seasons.

185 Relationships between sea ice and vital rates


186 Sea ice concentration affects emperor penguin vital rates through various mechanisms,
187 which are not mutually exclusive; see review by Ainley et al. (2010). First, SIC may directly
188 affect foraging; in years with dense sea ice cover, foraging trips to the nearest open water area
189 are on average longer, energetic costs for breeding adults are higher, and the provisioning of
190 chicks is lower (Massom et al., 2009; Zimmer et al., 2008). We would thus expect negative
191 effects of high SICa on breeding success and on adult survival of both sexes during the
192 rearing period.
193 Sea ice concentration is also critical to Antarctic ecosystem function (Thomas & Dieck-
194 mann, 2003). It may indirectly affect the emperor penguin through its effects on other species
195 of the Antarctic food web, either prey or predators (Ainley et al., 2010, 2007; Barbraud &
196 Cotte, 2008; Barbraud & Weimerskirch, 2001; Forcada & Trathan, 2009). However, the
197 mechanisms involved remain unclear and their relative importance is still debated (Ainley
198 et al., 2007; Barbraud & Cotte, 2008).
199 The effects of SIC on primary productivity and krill may cascade up to fish and upper

9
200 predators; e.g. Nicol et al. (2000); Forcada & Trathan (2009); Fraser & Hoffmann (2003).
201 Primary productivity and krill density are known to be related to sea ice extent (SIE) and
202 concentration (Atkinson et al., 2004; Brierley et al., 2002; Loeb et al., 1997). Around Antarc-
203 tica, krill density during summer is positively related to chlorophyll concentrations (Atkinson
204 et al., 2008). In the South West Atlantic, krill densities and recruitment during summer are
205 positively related to SIE of the previous winter (Loeb et al., 1997; Nicol et al., 2000; Atkinson
206 et al., 2004). Other studies have shown non-linear relations between sea ice cover and krill
207 populations (Quetin et al., 2007; Wiedenmann et al., 2009). High krill recruitment occurs
208 over a range of optimum sea ice conditions, suggesting complex mechanisms linking sea ice
209 and krill abundance.
210 Sea ice may also influence the emperor penguin by top-down processes (Ainley et al.,
211 2010, 2007). Reduced sea ice cover may allow greater access to the foraging areas of the
212 emperor penguin by potential predators such as killer whales (Orcinus orca (Ainley et al.,
213 2010, 2007; Pitman & Durban, 2010), although no emperor penguin remains have been
214 detected in the diet of killer whales (Barbraud & Cotte, 2008). Energetically compromised
215 penguins (especially males after their four months fast) may be particularly vulnerable to
216 predation (Ainley personal communication).
217 Parameter estimation
218 The capture-recapture dataset permits the estimation of survival and return probabilities
219 (see Jenouvrier et al. (2005a, 2010) for more details), which may differ between males and
220 females because of their different breeding investment (Barbraud & Weimerskirch, 2001; Je-
221 nouvrier et al., 2005a). We estimated probabilities of survival, and of return to the breeding
222 colony, using a multistate capture-mark-recapture model (CMR) based on sex and repro-
223 ductive status, with some unobservable stages. Adult survival is allowed to differ between
224 sexes, and return probability to differ between breeders and non-breeders. The statistical
225 capture-recapture model is described in detail in Appendix C of Jenouvrier et al. (2010); for
226 a review of these methods see Lebreton et al. (2009).

10
227 We estimated linear and quadratic effects of SICa using models relating the vital rates
228 to covariates with a general logit link function (Lebreton et al., 1992). It is impossible to
229 test the impact of SICa on survival during the first year at sea, or on the survival or return
230 probabilities of pre-breeders, because data on pre-breeders and SICa overlap for only a few
231 years (modern SIC satellite data are only available from 1979 onwards).
232 Estimation and CMR model selection was performed using the program M-SURGE (Cho-
233 quet et al., 2004). We used Akaike’s information criterion (AIC, Akaike (1974)) to compare
234 models, the model with the lowest AIC being the model most supported by the data. We
235 based our inferences on the most plausible set of models using model averaging using AIC
236 weights wi (Burnham & Anderson, 2002). Thus a parameter θ is calculated as

M
!
−1
X
θ = logit wi βi (1)
i=0

237 where M is the number of models and βi is the estimated parameter on the logit scale for
238 model i. The Akaike weights were calculated as

exp (−∆i /2)


wi = M
(2)
X
exp (−∆r /2)
r=0

239 where ∆i = AICi − min AIC, where min AIC is the smallest value of AIC in the model set.
240 Concerns have been raised that the use of flipper bands negatively affects the survival
241 and breeding success of other species of penguins (Saraux et al., 2011). The banding study
242 at Terre Adélie was discontinued in 1988 partly out of such concerns. However, the breeding
243 success of banded and unbanded emperor penguins at Terre Adélie is not significantly differ-
244 ent (Barbraud, unpublished results). Saraux et al. (2011) found that the effects of bands on
245 survival of king penguins disappeared 4.5 years after banding. To reduce possible impacts
246 of the flipper bands, but still retain sufficient sample size, we eliminated capture histories
247 for the first two years after banding. Circumstantial evidence that the possible effects of

11
248 bands are minimal in this population comes from estimates of population growth rate. If the
249 banded birds on which CMR estimates are based had lower survival probabilities, we would
250 expect estimated population growth rates to be consistently lower than those observed. This
251 is not the case: the population growth rate obtained from CMR estimates agrees well with
252 the observed population growth rate (Jenouvrier et al., 2010).
253 To quantify uncertainties resulting from model selection and estimation error, we used
254 the parametric bootstrap procedure introduced by Regehr et al. (2010) and Hunter et al.
255 (2010). To generate a bootstrap sample of a model output in this procedure, a CMR model
256 is first selected with probability proportional to its AIC weight. Then a vector of parameter
257 values is drawn from a multivariate normal distribution with a mean equal to the estimated
258 parameter vector and a covariance matrix obtained from the Hessian matrix of the CMR
259 model, or from the logistic regression model in the case of breeding success. The resulting
260 parameter vector is used to create the population model, and the output calculated from
261 that model. This process is repeated to generate a bootstrap sample, from which confidence
262 intervals can be calculated using the percentile method.
263 Results
264 Supplementary Appendix 3 details some of the results outlined here, including the un-
265 derlying parameter estimates and their uncertainties.
266 Breeding success is a decreasing function of SICa during the rearing period (Figure 3).
267 The CMR model selection procedure reveals effects of SICa on adult survival of both sexes
268 during all four seasons (Table 1). The models with ∆AIC ≤ 4 include effects in all seasons,
269 those with ∆AIC ≤ 3 include effects in the laying, incubation, and rearing seasons, and
270 those with ∆AIC ≤ 2 include effects during the incubation and rearing seasons. There was
271 no support for effects of SICa on the probability of return to the breeding colony (S3).
272 Figure 4 shows adult survival as a function of annual SICa and seasonal differences in
273 SICa . Survival probability is a concave nonlinear function of annual SICa (Figure 4.a,b).
274 The maximum annual adult survival is higher for females than for males (0.96 and 0.93

12
275 respectively). The effect of seasonal differences in SICa (Figure 4.c,d) is small compared to
276 the effect of annual SICa, but it has an important effect on the difference between male and
277 female survival. The effect of seasonal differences are positive when annual SICa < 0, and
278 negative when annual SICa > 0.

279 3 Influence of sea ice on population growth


280 To assess the effect of sea ice conditions on population growth, we use a demographic
281 model in which the parameters are functions of SICa . Because sea ice is strongly seasonal and
282 breeding biology is tied to the seasons, we use a seasonal periodic matrix model (Caswell 2001,
283 Chap. 13) to capture these effects. The model, which is described in detail in Jenouvrier et al.
284 (2010), includes a sequence of seasonal behaviors (arrival to the colony, mating, breeding)
285 and accounts for differences in adult survival between males and females. The model is
286 nonlinear because mating probability depends on the availability of males and females for
287 mating. This frequency dependence is captured by expressing reproduction as a function of
288 the proportional structure of the population.
289 In the model, a matrix Mi projects the population from season i to season i + 1. Since
290 we identify four seasons, M4 projects the population from season 4 to season 1 in the next
291 year. The annual projection matrix is given by the periodic product of the Mi :

A = M4 M3 M2 M1 . (3)

292 where M1 includes the birth process, M2 includes annual mortality process, M3 includes
293 migration to the breeding site, and M4 includes the mating process (Jenouvrier et al. 2010).
294 To establish clear notation, let x(t) represent the SICa values (in all four seasons) in year
295 t, and let p(t) be the proportional population structure in year t. The vector of demographic
296 parameters is a function of sea ice and population composition; we denote it by θ(x, p). Then

13
297 the population in year t grows according to

h i
n(t + 1) = A θ(x(t), p(t)) n(t) (4)

298 Deterministic population growth


299 In a fixed, specified sea ice environment, the population will eventually converge to an
300 equilibrium proportional structure p̂ and (even though it is nonlinear) grow exponentially
301 at a rate λ̂ given by the dominant eigenvalue of the projection matrix evaluated at the
302 equilibrium structure, A[θ(x, p̂)]. To compute p̂, we calculate SICa , during each of the four
303 seasons, from the annual SICa and the seasonal differences in SICa . We then project the
304 population from an arbitrary initial vector until it converges to the equilibrium structure, and
305 use that structure to compute λ̂. For a fixed value of annual SICa , the population growth
306 rate increases (annual SICa < −0.8) or decreases (annual SICa > −0.3) with increasingly
307 positive seasonal differences. Figure 5 shows λ̂ as a function of annual SICa and of seasonal
308 differences in SICa . The population growth rate is maximized at intermediate values of
309 SICa close to 0, and declines at higher or lower values.
310 The range of positive growth is wide (white contours on Figure 5) and λ̂ declines from
311 its maximum more rapidly for negative than for positive annual SICa values.

312 Stochastic population growth


313 To examine the effect of sea ice variability, we calculate the stochastic population growth
314 rate as a function of the means and variances of annual SICa and of seasonal differences in
315 SICa . The growth rates are calculated from stochastic simulations. At each time step t,
316 values for annual SICa and seasonal differences in SICa are drawn from normal distributions
317 with specified means and variances, and are used to parameterize the projection matrix A
318 in eq.(3). We also include stochastic variation, unrelated to sea ice conditions, in breeding
319 success and the probability of return to the colony. For breeding success, we add a normally
320 distributed error term with variance given by the residual variation around the logistic re-

14
321 gression of breeding success on SICa (Figure 3). For probability of return, which is not a
322 function of SICa , we sample repeatedly from the set of measured return probabilities from
323 1970 to 2000.
324 The simulation begins with an arbitrary population vector n0 , and projects the population
325 according to
h i
n(t + 1) = A θ (x(t), p(t)) n(t). (5)

326 The stochastic growth rate is given by

1
log λs = lim log kA [θ(T − 1)] · · · A [θ(0)] p(0)k (6)
T →∞ T

327 We evaluate log λs numerically using T = 50, 000.


328 Figure 6 shows the results and Appendix S4 provides more details on the distribution of
329 sea ice, vital rates and population growth rate for various examples. The stochastic growth
330 rate is maximized at intermediate values of mean annual SICa , and declines as anomalies
331 become very positive or very negative. At very high or very low values of the mean annual
332 SICa the growth rate is improved by increasing variance, but at intermediate values, the
333 effect of variance is negative. The effect of the mean and variance in seasonal differences in
334 SICa are smaller, and negative. Changes in the mean and/or the variance of annual SICa
335 have the potential to greatly reduce the stochastic growth rate.

336 4 Population response to climate change


337 To project the population response to climate change, we use our demographic model
338 to determine the response of the population to future sea ice conditions as forecasted by a
339 select set of GCMs. We obtain forecasts of SICa from a set of GCMs, compute stochastic
340 SICa forecasts from these, and use the results to generate population trajectories from 2000
341 to 2100.

15
342 Stochastic sea ice forecasts
343 GCMs differ in their ability to reproduce sea ice conditions in Antarctica. Thus, from an
344 initial set of 20 climate models (Table 2), we selected those for which the statistical properties
345 of the distribution of SICa output agree well with the observations from 1979 to 2007, in
346 terms of both the median and the standard deviation of the SICa distribution (see S5 for
347 details). From the original set of 20 climate models, five were selected: cccma-cgcm3-1,
348 cccma-cgcm3-1-t63, mpi-echam5, ukmo-hadcm3, and ukmo-hadgem1 (Table 2, see Figure in
349 Appendix S5).
350 These climate models were forced with a middle range emissions scenario (SRES A1B,
351 IPCC, 2000). This scenario assumes a future socio-economic development depending on
352 fossil and non-fossil energy sources in balanced proportions. Under this scenario, the CO2
353 level doubles by 2100, from 360 ppm to 720 ppm.
354 To generate stochastic SICa forecasts, we first obtain output for SICa in each of the
355 four seasons. From this output, we compute smoothed means x̄(t) and smoothed covariance
356 matrices C̄(t), using a Gaussian kernel smoother with smoothing parameter h = 2 (Appendix
357 S5). We then generate stochastic SICa vectors by drawing x(t) as an iid sample from a
358 normal distribution with mean x̄(t) and covariance matrix C̄(t).
359 All five GCMs agree that the smoothed mean SICa will decline by 2100, but the rate
360 of decline varies between climate models and seasons (Figure 7). For example, SICa during
361 the laying season is forecast to decline only by 9% relative to present for the model ukmo-
362 hadcm3, while SICa during the non-breeding season will decline by 71% according to the
363 model cccma-cgm3-t63. There is no clear pattern of change in the smoothed variance (Figure
364 8), but a high variability over time and between models seems likely. Some models predict a
365 decline in the smoothed variance by the end of the century (e.g. cccma-cgm3-t63 during the
366 non-breeding season), while other an increase (e.g. ukmo-hadcm3 during the laying season).

16
367 Stochastic population projections
368 We used each stochastic SICa forecast to generate a sequence of demographic rates from
369 2010 to 2100 (Appendix S6). These rates were used to project the population using as an
370 initial population vector the average equilibrium population structure from all the GCMs
371 in 2010. To evaluate uncertainties in climate, we used 200 stochastic forecasts from each of
372 the five GCMs. To evaluate uncertainties in demography, we use the parametric bootstrap
373 approach to generate a sample of 200 simulations for each sea ice forecast. Thus, we project
374 40,000 population trajectories for each GCM, for a total of 200,000 population trajectories.
375 The population projections exhibit considerable variability (Figure 9). Some projections
376 produce dramatic declines in the number of breeding pairs (e.g., projections from cccma-
377 cgcm3-1), while a few produce large increases (e.g., ukmo-hadcm3). The median of the 40,000
378 trajectories from each GCM has a unique pattern (Figure 9 and 10). Some increase gradually
379 (e.g., mpi-echam 5), while some decline gradually (e.g., ukmo-hadgem1). Some remain stable
380 for a while and then decline (e.g., cccma-cgcm3-1-t63). For each GCM, however, there exists
381 a year beyond which the median projection declines; this tipping year may be late (e.g., 2089
382 for ukmo-hadcm3) or early (e.g., 2038 cccma-cgcm3-1-t63). By the end of the century, the
383 medians of all GCMs except ukmo-hadcm3 project that the number of breeding pairs will
384 decline compared to the minimum number over the last six decades.
385 The median over the entire set of simulations declines to 575 breeding pairs by 2100. Over
386 this set of simulations, the probability of a decline by 90% or more by 2100 is 0.43 (Table
387 3). By 2100, the probability of a decline below the maximum number of breeding pairs ever
388 observed in Terre Adélie since 1962 is 93% (Table 4). Therefore, only 7% of population
389 trajectories included in this set are projected to increase by 2100.

390 5 Discussion
391 If the climate during the rest of this century follows the patterns forecast by the GCMs
392 examined here, the emperor penguin population at Terre Adélie will respond by declining

17
393 towards extinction. Our median projection shows a decline in the number of breeding pairs
394 by 81% over this period, and a good chance (43%) of a more severe decline of 90% or more.
395 By the end of our projection, the population is continuing to decline, regardless of which
396 GCM is used to forecast future sea ice conditions. The range of uncertainty associated with
397 this result might change the details, but not the overall biological conclusions.
398 We arrived at this conclusion by tracing the effects of sea ice concentration from the level
399 of the individual to the level of the population. First we measured the response of the vital
400 rates to sea ice conditions, then incorporated those responses into a demographic model to
401 calculate the population growth, and finally we coupled the demographic model to forecasts
402 of sea ice conditions produced by an ensemble of GCMs.

403 Effects of sea ice on the vital rates


404 The vital rates consist of male and female adult survival, the probability of returning
405 to the colony to breed, and breeding success. Previous studies have shown that breeding
406 success and adult survival have the biggest impacts on population growth rate, and that
407 return probability has only a small effect (Jenouvrier et al., 2005a, 2009a).
408 We found that breeding success declines with increasing values of SICa during the rearing
409 season. Years with high concentration sea ice may require longer foraging trips, reducing the
410 provisioning of chicks and thus breeding success (Massom et al., 2009; Wienecke & Robertson,
411 1997). Not all the variability in breeding success is explained by SICa during the rearing
412 season, probably because it is affected by many factors; e.g., prolonged blizzards and colder
413 temperatures may increase chick mortality (Jouventin, 1975), and premature ice break-out
414 may cause massive fledging failures (Budd, 1962). We include this unexplained variability
415 in our sea ice-dependent demographic model.
416 Sea ice affects adult survival in several ways, more complicated than the linear relation-
417 ships assumed in previous studies (Barbraud & Weimerskirch, 2001; Jenouvrier et al., 2005a).
418 The effect of annual SICa on adult survival is greater than the effect of seasonal difference
419 in SICa . Adult survival is maximized at positive annual SIC anomalies ' 2 and declines

18
420 otherwise. The maximum annual survival is higher for females than males because males
421 are likely more constrained energetically due to their long fasting period (Jenouvrier et al.,
422 2005a). Males incubate the egg, and thus fast for four consecutive months, whereas females
423 return to open water to forage after a two months fasting period. The response of survival
424 to sea ice can be explained by several non-exclusive mechanisms affecting sea ice habitat and
425 constraining the energy expenditure of emperor penguins. During years of concentrated sea
426 ice, foraging trips may be longer, resulting in higher energetic costs (Wienecke & Robertson,
427 1997; Massom et al., 2009). On the other hand, during years of concentrated sea ice, food
428 resources may be higher (Barbraud & Weimerskirch, 2001) and/ or predation lower (Ainley
429 et al., 2007). Seasonal difference in SICa affects the survival difference between females and
430 males, suggesting that the sexes respond differently to SICa during the different seasons.
431 Although the mechanisms remain unclear, this difference is probably linked to contrasted
432 energetic costs during breeding.
433 Male-female survival differences affect the population structure and growth, directly
434 through mortality or indirectly by limiting the availability of mates, and thus reproduction,
435 because penguins are monogamous (Jenouvrier et al., 2010).

436 Effect of sea ice on population growth rate


437 Population growth rate provides a measure of the quality of the environment, in terms of
438 the fitness of a population occupying that environment (Caswell, 2001). The largest effect
439 of sea ice on the population growth of the emperor penguin is due to the annual SICa ; the
440 effects of seasonal differences are smaller, but still appreciable.
441 The deterministic growth rate λ is maximized at intermediate values of annual SICa . The
442 maximum occurs at a value of SICa ≈ 0 (depending slightly on the value of the seasonal
443 difference in SICa ). The optimum is relatively broad, as shown by the wide range of annual
444 SICa values enclosed by the contours for log λ = 0 in Figure 5. This implies that changes
445 in sea ice conditions have little effect on λ until annual SICa becomes quite positive or
446 negative. This intermediate optimum is expected because neither complete absence of sea

19
447 ice, nor heavy and persistent sea ice (i.e. no access to resources through polynya) would
448 provide satisfactory conditions for the emperor penguin (Ainley et al., 2010).
449 A complicated interaction exists between annual SICa and seasonal differences in SICa .
450 When annual SICa is lower than −0.8 , the effect of seasonal differences is positive (i.e., the
451 emperor penguin performs better when sea ice concentration is higher in breeding season
452 and lower in non-breeding season). But, when annual SICa is higher than −0.3, the effect of
453 seasonal differences is negative; the emperor penguin does better when the seasonal pattern
454 is the opposite.
455 The emperor penguin population is more sensitive to negative than to positive annual
456 sea ice anomalies; i.e., λ decreases from its optimal value faster in the negative than the
457 positive direction (Figure 5). This is well illustrated by the dramatic 50% population decline
458 in the late 1970s in Terre Adélie, which coincided with several years of the lowest sea ice
459 extent ever recorded during the last 40 years (Jenouvrier et al., 2005c). Other colonies have
460 disappeared in regions with high temperature and low sea ice duration (Fretwell & Trathan,
461 2009; Trathan et al., 2011). The Dion Islets colony along the west coast of the Antarctic
462 Peninsula (67 52S, 68 43W) declined from 250 breeding pairs in the 1970s to 20 pairs in
463 1999, and was extinct by 2009 (Ainley et al., 2010; Trathan et al., 2011). This extinction
464 coincided with a decline in sea ice duration, resulting from a warming of the west coast of
465 the Antarctic Peninsula at an unprecedented rate (Vaughan et al., 2001).
466 The stochastic population growth rate log λs shows a similar response to mean annual
467 SICa , with a broad maximum when the mean annual SICa ≈ 0. The effect of mean seasonal
468 difference in SICa is smaller and, as in the deterministic case, is negative when evaluated at
469 a mean annual value of SICa = 0 (Figure 6).
470 The effects of the variance in annual SICa depend on the value of the mean annual SICa ,
471 as shown in Figure 6a. Within the range of approximately −3 ≤ SICa ≤ 4, the effect of
472 variance is negative, and log λs is maximized when the variance is zero. But for mean annual
473 SICa outside this range, the effect of variance on log λs is positive.

20
474 It is well known that, all else being equal, temporal variance in the vital rates reduces
475 the stochastic growth rate; covariances among vital rates and temporal autocorrelation can
476 reverse this conclusion (Caswell, 2001; Lewontin & Cohen, 1969; Tuljapurkar, 1990). How-
477 ever, this conclusion does not translate directly to the effect of variance in environmental
478 factors. Because the vital rates are a nonlinear function of sea ice, and the two-sex model is
479 itself nonlinear, the effect of environmental (i.e., SICa ) variance on population growth may
480 be either positive or negative (e.g., Koons et al. (2009)).

481 Effect of future sea ice change


482 Our best projection of the future growth of the Terre Adélie emperor penguin population,
483 under the impact of climate change, is a decline by the year 2100 from approximately 3000
484 breeding pairs to 575 breeding pairs, a decline of 81%, or an average rate of change of −0.018
485 per year. This projection is the median of a large set of simulations that incorporate as many
486 sources of uncertainty as possible.
487 The median is a smoothed pattern based on 200,000 population trajectories. Each of
488 those trajectories fluctuates, with increases, decreases, and periods of relative stability, until
489 it reaches a tipping year, after which it declines. The tipping year varies among trajectories,
490 depending on the forecasted sea ice, the responses of the vital rates, and the impact of
491 sex-specific adult survival on the demography. The decline in the median of the 200000
492 population trajectories accelerates after 2040 because more of the population trajectories
493 are likely to have reached their tipping year as time goes on.
494 These population projections required three steps.

495 1. We had to extract the biologically relevant GCM outputs, in our case the seasonal
496 SICa values, on appropriate spatial scales. The choice of spatial scales is an important
497 issue in studies of climate change using GCMs, which project sea ice variables, such as
498 concentration, thickness (SIT) and extent, over a greater spatial scale than the scale
499 of emperor penguin habitat requirements. For emperor penguins, the size of polynyas
500 (Ainley et al., 2005), sea ice thickness, area of fast ice (Massom et al., 2009), the timing

21
501 of ice breakup and formation are meaningful variables with respect to the life cycle
502 (see review in Ainley et al. (2010)) but are measurable only at small spatial scales. By
503 focusing on SIC we hope to reduce the number of correlated covariates (Grosbois et al.,
504 2008). At large spatial scales, SIC is strongly correlated to SIT and open water area,
505 and thus SIC is a good measure of the sea ice environment experienced by penguins.
506 Further climate research will be required to downscale the sea ice projections of GCMs
507 in Antarctica both spatially (e.g. RCMs) and temporally (e.g. daily data to calculate
508 the timing of ice breakup) as climate model output is typically saved and available
509 with a monthly resolution.

510 2. We had to model temporal variance in SICa from the GCM output. Ideally, this would
511 be obtained from multiple stochastic realizations of each GCM, but such output does
512 not exist. Thus, we obtained smoothed temporal means, variances, and covariances
513 from the output, and used those to parameterize variability in SICa at each time step.
514 The forecasts of smoothed temporal means and variances of SICa differ strongly among
515 climate models. Differences are present in numerous aspects of the climate models,
516 including their physical parameterizations and spatial resolution. As such, it is difficult
517 to attribute differences in the projected sea ice change to a single factor. However, as
518 noted by Lefebvre & Goosse (2008), numerous mechanisms do appear to play a role.
519 These include the influence of local simulated climate conditions at the end of the 20th
520 century, aspects of the atmospheric Southern Annular Mode response (e.g. Fyfe &
521 O.A. (2006)), and changes in the simulated response of the Southern Ocean circulation.
522 While our climate model selection process does reduce the uncertainty associated with
523 the first of these, differences in simulated feedbacks and climate response will still play
524 an important role. This points to a further need for climate model development and
525 enhancement.

526 3. It was critical to incorporate uncertainty in our projections; a projection of an 81%


527 decline without an associated range of uncertainty would be useless. Our uncertainty

22
528 analysis included components from the demographic parameters and their response to
529 sea ice conditions, including both sampling error and the uncertainty due to model
530 selection. We also included uncertainty in future sea ice conditions by using an ensem-
531 ble of GCMs, chosen on the basis of their agreement with past sea ice observations.
532 Figure 9 shows the variation among the models in their forecasts of SICa , and Figure
533 10 provides a powerful graphical summary of the resulting uncertainty in population
534 forecasts. Models cccma-cgcm3-1 and cccma-cgcm3-1-t63 predict the most rapid de-
535 clines, with the population reduced to 20 or 11 breeding pairs, respectively, by 2100. In
536 contrast, model ukmo-hadcm-3 makes more optimistic sea ice forecasts, and projects
537 an increase in the population by 2100.

538 Stochasticity and uncertainty naturally lead to significant variation among projections
539 of future emperor penguin population growth. Among the 200,000 population trajectories,
540 there are some examples in which the population does not decline, or even increaes, but the
541 central tendency is an unambiguous and serious population decline of 81%. More important,
542 the probability distributions of projected population size in 2100 show that declines are far
543 more likely than increases or stability. In addition, the median population trends predicted
544 by all five GCMs during the last decade (2090-2100), even the outlier ukmo-hadcm3, are
545 negative. Thus the difference between climate models, i.e. uncertainty in the sea ice forecasts,
546 affects the timing of the population median decline, but not whether that median decline
547 occurs or not (Figure 10).
548 The conclusion that the emperor penguin population will decline dramatically by the end
549 of this century raises the issue of possible adaptation. Emperor penguins might adapt to
550 the new sea ice conditions or, more likely, disperse to locations where sea ice conditions are
551 more favorable (see review Ainley et al. (2010); Forcada & Trathan (2009) and discussion in
552 Trathan et al. (2011); Jenouvrier et al. (2009b)). Future studies should quantify potential
553 refuges for the species (but see Ainley et al. (2010)) and consider potential evolutionary
554 responses. Thus, we encourage ecologists to collect information on phenotypic traits and

23
555 their heritability.
556 The median projection of 575 breeding pairs in 2100 is close to the projection of ≈ 400
557 breeding pairs obtained previously using a simpler modelling approach (Jenouvrier et al.
558 2009). Consistencies among population projections from different demographic and climate
559 models increase confidence in the assessment of the impact of climate change.

560 Methodological notes


561 Our approach to predicting and understanding the effects of climate change requires the
562 measurement of effects of climate on individuals, the integration of those effects into a demo-
563 graphic model, and the connection of the demographic model to climate forecasts. Aspects of
564 this approach have been applied to some other species, including seabirds (Barbraud et al.,
565 2010; Wolf et al., 2010), polar bears (Hunter et al., 2010), the oystercatcher (van de Pol
566 et al., 2010), and our previous analysis of the emperor penguin (Jenouvrier et al., 2009b).
567 In this study, we have gone beyond these studies in several ways.
568 Because of the length and quality of the emperor penguin data set, in this analysis we
569 were able to identify detailed effects of sea ice on a seasonal basis through the breeding cycle,
570 and to include, for the first time, the nonlinear effects of sex-specific mortality patterns using
571 a two-sex demographic model. In the case of the emperor penguin, the extreme conditions
572 of its breeding cycle make the presence of both parents essential for successful reproduction.
573 Hence, sex-specific climate effects are particularly important (cf. Jenouvrier et al. (2010)).
574 In general, however, we expect that even in less dramatic conditions, the different roles of
575 the sexes may often cause climate change impacts to differ between males and females.
576 Because of the ability to extract relevant sea ice output from the GCMs, we were able
577 to include both stochasticity and uncertainty, and to draw conclusions about the projected
578 fate of the Terre Adélie emperor penguin population even in the face of that uncertainty.
579 Because climate models predict that both the mean and variability of climate will change
580 (Solomon et al., 2007), it is important to include these stochastic effects.
581 Two factors not included in our models deserve mention: density dependence and demo-

24
582 graphic stochasticity. It is almost reflexive to ask about density dependence in an analysis
583 such as ours. This is a result of exposure to simple scalar population models. In a stage-
584 structured population, especially in a highly seasonal environment, and in which different
585 stages occupy vastly different environments during the course of the year, population growth
586 is never likely to be a function of something as simple as “density.” In the case of the em-
587 peror penguin, it is a priori unlikely that density has strong effects, and in Appendix 1, 3
588 we describe our analysis of the effects of the number of breeding pairs on the vital rates; we
589 found no support for these effects.
590 Demographic stochasticity refers to fluctuations arising from the random outcome of
591 probabilities of survival and reproduction, applied to individuals in a population. It can be
592 analyzed using multi-type branching process models (Chapter 15 of Caswell 2011). Demogr-
593 paphic stochasticity is unimportant in large populations, but reduces the stochastic growth
594 rate and increases extinction probability in small populations. As the emperor penguin popu-
595 lation declines, demographic stochasticity will, at some point, become important. In general,
596 however, this requires population sizes of only a few tens to a hundred individuals. To the
597 extent that demographic stochasticity becomes important at the end of our simulations, our
598 projections overestimate the persistence of the population.
599 We encourage further collaborations between ecologists and climatologists. The devel-
600 opment of data archive resources such as those provided by the National Snow and Ice
601 Data Center (NSIDC) and the Program for Climate Model Diagnosis and Intercomparison
602 (PCMDI) (Meehl et al., 2007) have allowed unprecedented access to observed and modeled
603 climate data. We believe that the participation of climatologists is critical for selecting the
604 most appropriate set of climate models, emissions scenarios, climate variables, and the use
605 of climate ensembles versus single climate outputs.
606 An important future step will be to incorporate evolutionary responses. A theoretical
607 framework to link life history traits and population dynamics (eco-evolutionary dynamics)
608 is emerging (Hoffmann & Sgrò, 2011; Pelletier et al., 2009; Kinnison & Hairston, 2007; Reed

25
609 et al., 2010), and studies of the potential evolutionary responses to climate change and their
610 population consequences have recently been initiated (Jenouvrier & Visser (2011); Reed et al.
611 (2011)). We believe that the way climatologists have approached the problem, using coupled
612 climate system models in which climate systems components (e.g. ocean, atmosphere, sea
613 ice, land surface, ice sheets, biogeochemistry and more) are gradually connected, can provide
614 a valuable example to ecologists toward an integrative climate- eco-evolutionary framework.

615 6 Acknowledgments
616 The penguin data come from a long-term study supported by Expéditions Polaires
617 Françaises, by Institut Paul Emile Victor (Programme IPEV 109), and by Terres Australes
618 et Antarctiques Françaises. We acknowledge the efforts of all the wintering fieldworkers
619 involved in the long-term monitoring programs in Terre Adélie since 1963, and thank Do-
620 minique Besson and Karine Delord for the management of the database. We acknowledge the
621 modeling groups, the Program for Climate Model Diagnosis and Intercomparison (PCMDI)
622 and the WCRP’s Working Group on Coupled Modelling (WGCM) for their roles in making
623 available the WCRP CMIP3 multi-model dataset. Support of this dataset is provided by
624 the Office of Science, U.S. Department of Energy. HW, CB and SJ acknowledge support
625 from the REMIGE (Behavioural and demographic REsponses of Indian Ocean Marine top
626 predators to Global Environmental changes) program funded by ANR (Agence Nationale de
627 la Recherche) Biodiversité (ANR Biodiv 011). This work was carried out during the tenure
628 of Marie-Curie and CIRES Visiting fellowships by S.J. MH acknowledges support through
629 the National Science Foundation. HC acknowledges support from NSF Grant DEB-0816514,
630 from the WHOI Arctic Research Initiative, and from the Alexander von Humboldt Founda-
631 tion. The authors thank the two anonymous reviewers and David Koons for their valuable
632 comments and suggestions on previous manuscript version.

26
633 References
634 Adahl E, Lundberg P, Jonzen N (2006) From climate change to population change: the need
635 to consider annual life cycles. Global Change Biology, 12, 1627–1633.

636 Ainley D, Ballard G, Ackley S, et al. (2007) Paradgim lost, or is top-down forcing no longer
637 significant in the antarctic ecosystem? Antarctic Science.

638 Ainley D, Clarke E, Arrigo K, Fraser B, Kato A, Barton K, Wilson P (2005) Decadal-scale
639 changes in the climate and biota of the pacific sector of the southern ocean, 1950s to the
640 1990s. Antarctic Science, 17, 171–182.

641 Ainley D, Russell J, Jenouvrier S, Woehler E, Lyver P, Fraser W, Kooyman G (2010) Antarc-
642 tic penguin response to habitat change as earths troposphere reaches 2 degrees above
643 preindustrial levels. Ecological Monographs, 80, 49–66.

644 Akaike H (1974) A new look at the statistical model identification. IEEE Transactions on
645 Automatic Control, 19 (6):, 716723.

646 Atkinson A, Siegel A, Pakhomov E, Rothery P (2004) Long-term decline in krill stock and
647 increase in salps within the southern ocean. Nature, 432, 100–103.

648 Atkinson A, Siegel V, Pakhomov E, et al. (2008) Oceanic circumpolar habitats of antarctic
649 krill. Marine Ecology Progress Series, 362, 1–23.

650 Ballerini T, Tavecchia G, Olmastroni S, Pezzo F, Focardi S (2009) Nonlinear effects of winter
651 sea ice on the survival probabilities of adélie penguins. Oecologia, 161, 253–265.

652 Barbraud C, Cotte C (2008) Short note: Paradigms need hypothesis testing: no evidence
653 for top-down forcing on adélie and emperor penguin populations. Antarctic Science, 20,
654 2.

27
655 Barbraud C, Rivalan P, Inchausti P, Nevoux M, Rolland V, Weimerskirch H (2010) Con-
656 trasted demographic responses facing future climate change in southern ocean seabirds.
657 Journal of Animal Ecology, 80, 89–100.

658 Barbraud C, Weimerskirch H (2001) Emperor penguins and climate change. Nature, 411,
659 183–186.

660 Both C, Bouwhuis S, Lessells CM, Visser ME (2006) Climate change and population declines
661 in a long-distance migratory bird. Nature, 441, 81–83.

662 Brierley A, Fernandes P, Brandon M, et al. (2002) Antarctic krill under sea ice : Elevated
663 abundance in a narrow band just south of ice edge. Science, 295, 1890–1892.

664 Budd GM (1962) opulation studies in rookeries of the emperor penguin aptenodytes forsteri.
665 Proceedings of the Zoological Society of London, 139, 365.

666 Burnham K, Anderson D (2002) Model Selection and Multimodel Inference: A Practical
667 Information-theoretic Approach. Springer.

668 Caswell H (2001) Matrix population models, vol. Second. Sinauer, Sunderland, Mas-
669 sachusetts.

670 Cavalieri D, Parkinson C, Gloersen P, Zwally HJ (1996) Sea ice concentrations from nimbus-
671 7 smmr and dmsp ssm/i passive microwave data. Boulder, Colorado USA: National Snow
672 and Ice Data Center. Digital media.

673 Cherel Y (2008) Isotopic niches of emperor and adélie penguins in adélie land, antarctica.
674 Marine Biology, 154, 813–821.

675 Cherel Y, Kooyman G (1998) Food of emperor penguin (Aptenodytes forsteri) in the western
676 ross sea, antarctica.hu. Marine Biology, 130, 335–344.

28
677 Choquet R, Reboulet A, Pradel R, Giménez O (2004) M-surge: new software specifically
678 designed for multistate capture-recapture models. Animal Biodiversity and Conservation,
679 27, 207–215.

680 Croxall J, Trathan P, Murphy E (2002) Environmental change and antarctic seabirds popu-
681 lations. Science, 297, 1510–1514.

682 Forcada J, Trathan PN (2009) Penguin responses to climate change in the southern ocean.
683 Global Change Biology, 15, 1618–1630.

684 Fraser W, Hoffmann E (2003) A predator’s perspective on causal links between climate
685 change, physical forcing, and ecosystem response. Marine Ecology Progress Series, 265,
686 1–15.

687 Fretwell PT, Trathan PN (2009) Penguins from space: faecal stains reveal the location of
688 emperor penguin colonies. Global Ecology and Biogeography, 18, 543–552.

689 Fyfe J, OA S (2006) Simulated changes in the extratropical southern hemisphere winds and
690 currents. Geophysical Research Letters, 33, L06701.

691 Grosbois V, Gimenez O, Gaillard J, et al. (2008) Assessing the impact of climate variation
692 on survival in vertebrate populations. Biological Reviews, 83, 357–399.

693 Herrando-Pérez S, Delean S, Brook B, Bradshaw CJA (2012) Decoupling of component and
694 ensemble density feedbacks in birds and mammals. Ecology.

695 Hoffmann AA, Sgrò CM (2011) Climate change and evolutionary adaptation. Nature, 470,
696 479–485.

697 Hughes L (2000) Biological consequences of global warming: is the signal already apparent?
698 Trends in Ecology and Evolution, 15, 56–61.

29
699 Hunter C, Caswell H, Runge M, Regehr E, Amstrup S, Stirling I (2010) Climate change
700 threatens polar bear populations: a stochastic demographic analysis. Ecology, 91, 2883–
701 2897.

702 Jenouvrier S, Barbraud C, Weimerskirch H (2003) Effects of climate variability on the tem-
703 poral population dynamics of southern fulmars. Journal of Animal Ecology, 72, 576–587.

704 Jenouvrier S, Barbraud C, Weimerskirch H (2005a) Long-term contrasted responses to cli-


705 mate of two antarctic seabirds species. Ecology, 86, 2889–2903.

706 Jenouvrier S, Barbraud C, Weimerskirch H (2005b) Sea ice affects the population dynamics
707 of adélie penguins in terre adélie. Polar Biology, 29, 413–423.

708 Jenouvrier S, Barbraud C, Weimerskirch Hand Caswell H (2009a) Limitation of population


709 recovery: a stochastic approach to the case of the emperor penguin. Oikos, 118, 1292–
710 1298.

711 Jenouvrier S, Caswell H, Barbraud C, Holland M, Stroeve J, Weimerskirch H (2009b) De-


712 mographic models and ipcc climate projections predict the decline of an emperor penguin
713 population. Proceedings of the National Academy of Sciences, 106, 1844–1847.

714 Jenouvrier S, Caswell H, Barbraud C, Weimerskirch H (2010) Mating behavior, population


715 growth, and the operational sex ratio: A periodic two-sex model approach. The American
716 naturalist, 175, 739–752.

717 Jenouvrier S, Visser M (2011) Climate change, phenological shifts, eco-evolutionary responses
718 and population viability: toward a unifying predictive approach. Journal of Biometeorol-
719 ogy, 458, 1–15.

720 Jenouvrier S, Weimerskirch H, Barbraud C, Park YH, Cazelles B (2005c) Evidence of a shift
721 in cyclicity of antarctic seabirds dynamics link to climate. Proceedings of the Royal Society
722 of London.B., 272, 887–895.

30
723 Jouventin P (1975) The Biology of Penguins., chap. Mortality parameters in emperor pen-
724 guins Aptenodytes forsteri., p. 435446. Macmillan Press, London.

725 Kinnison M, Hairston N (2007) Eco-evolutionary conservation biology: contemporary evo-


726 lution and the dynamics of persistence. Functional Ecology, 21, 444–454.

727 Kirkwood R, Robertson G (1997) The foraging ecology of female emperor penguins in winter.
728 Ecological Monograph, 67, 155–176.

729 Koons DN, Pavard S, Baudisch A, Metcalf CJE (2009) Is life-history buffering or lability
730 adaptive in stochastic environments? Oikos, 118, 972–980.

731 Lebreton J, Burnham K, Clobert J, Anderson D (1992) Modeling survival and testing bio-
732 logical hypotheses using marked animals: a unified approach with case studies. Ecological
733 Monographs, 62, 67–118.

734 Lebreton JD, Nichols JD, Barker RJ, Pradel R, Spendelow JA (2009) Modeling individ-
735 ual animal histories with multistate capture-recapture models. Advances in Ecological
736 Research, 41, 87–173.

737 Lefebvre W, Goosse H (2008) Analysis of the projected regional sea-ice changes in the south-
738 ern ocean during the twenty-first century. Climate dynamics, 30, 59–76.

739 Lewontin RC, Cohen D (1969) On population growth in a randomly varying environment.
740 Proceedings of the National Academy of Sciences, pp. 1056–1060.

741 Loeb V, Siegel V, Holm-Hansen O, Hewitt R, Fraser W, Trivelpiece S (1997) Effects of sea-ice
742 extend and krill or salp dominance on the antarctic food web. Nature, 387, 897–900.

743 Massom R, Hill K, Barbraud C, Adams N, Ancel A, Emmerson L, Pook M (2009) Fast ice
744 distribution in adélie land, east antarctica: interannual variability and implications for
745 emperor penguins aptenodytes forsteri. Marine Ecology Progress Series, 374, 243–257.

31
746 McCarty J (2001) Ecological consequences of recent climate change. Conservation Biology,
747 15, 320–331.

748 Meehl GA, Covey C, Delworth T, et al. (2007) The wcrp cmip3 multi-model dataset: A
749 new era in climate change research. Bulletin of the American Meteorological Society, 88,
750 1383–1394.

751 Moller AP, Rubolini D, Lehikoinen E (2008) Populations of migratory bird species that
752 did not show a phenological response to climate change are declining. Proceedings of the
753 National Academy of Sciences, 105, 16195–16200.

754 Nicol S, Pauly T, Bindoff N, et al. (2000) Ocean circulation off east antarctica affects ecosys-
755 tem structure and sea-ice extent. Nature, 406, 504–507.

756 Parmesan C (2006) Ecological and evolutionary responses to recent climate change. Annual
757 Review of Ecology, Evolution, and Systematics, 37, 637–669.

758 Parmesan C, Yohe G (2003) A globally coherent fingerprint of climate change impacts across
759 natural systems. Nature, 421, 37–42.

760 Pelletier F, Garant D, Hendry A (2009) Eco-evolutionary dynamics. Philosophical Transac-


761 tions B, 364, 1483–1489.

762 Pitman RL, Durban JW (2010) Killer whale predation on penguins in antarctica. Polar
763 Biology, 33, 1589–1594.

764 Prevost J (1961) Ecologie du manchot empereur. E.P.F., 222.

765 Quetin L, Ross R, Fritsen C, Vernet M (2007) Ecological responses of antarctic krill to
766 environmental variability: can we predict the future? Antarctic Science, 19, 253–266.

767 Reed T, Schindler D, Waples RS (2010) Interacting effects of phenotypic plasticity and
768 evolution on population persistence in a changing climate. Conservation Biology, 25,
769 56–63.

32
770 Reed TE, Schindler DE, Hague MJ, Patterson DA, Meir E, Waples RS, Hinch SG (2011)
771 Time to evolve? potential evolutionary responses of fraser river sockeye salmon to climate
772 change and effects on persistence. Plos One, 6, e20380.

773 Regehr EV, Hunter CM, Caswell H, Amstrup SC, Stirling I (2010) Survival and breeding of
774 polar bears in the southern beaufort sea in relation to sea ice. Journal of Animal Ecology,
775 79, 117–127.

776 Saraux C, Bohec CL, Durant JM, et al. (2011) Reliability of flipper-banded penguins as
777 indicators of climate change. Nature, 469, 203–206.

778 Solomon S, D Qin MM, Chen Z, Marquis M, Averyt K, Tignor M, Miller H (2007) Contribu-
779 tion of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel
780 on Climate Change, 2007. Cambridge University Press, Cambridge, United Kingdom and
781 New York, NY, USA.

782 Stenseth N, Mysterud A, Ottersen G, Hurell J, Chan K, Lima M (2002) Ecological effects
783 of climate fluctuations. Science, 297, 1292–1295.

784 Stock CA, Alexander MA, Bond NA, et al. (2011) On the use of ipcc-class models to assess
785 the impact of climate on living marine resources. Progress in Oceanography, 88, 1–27.

786 Thomas D, Dieckmann G (2003) Sea ice: an introduction to its physics, chemistry, biology
787 and geology. Oxford: Blackwell Science, 240-266 pp.

788 Trathan PN, Fretwell PT, Stonehouse B (2011) First recorded loss of an emperor penguin
789 colony in the recent period of antarctic regional warming: Implications for other colonies.
790 Plos One, 6, e14738.

791 Tuljapurkar S (1990) Delayed reproduction and fitness in variable environments. Proceedings
792 of the National Academy of Sciences, 87, 1139.

33
793 van de Pol M, Vindenes Y, Saether BE, Engen S, Ens BJ, Oosterbeek K, Tinbergen JM
794 (2010) Effects of climate change and variability on population dynamics in a long-lived
795 shorebird. Ecology, 91, 1192–1204.

796 Vaughan D, Marshall G, Connolley W, King J, Mulvaney R (2001) Devil in the detail.
797 Science, 293, 1777–1779.

798 Visser M, van Noordwijk AJand Tinbergen J, Lessells C (1998) Warmer springs lead to
799 mistimed reproduction in great tits (Parus major ). Proceedings of the Royal Society B:
800 Biological Sciences, 265, 18671870.

801 Vitousek P (1994) Beyond global warming: Ecology and global change. Ecology, 75, 1861–
802 1876.

803 Walther G, Post E, Convey P, et al. (2002) Ecological responses to recent climate change.
804 Nature, 416, 389–395.

805 Wiedenmann J, Cresswell K, Mangel M (2009) Connecting recruitment of antarctic krill and
806 sea ice. Limnology and Oceanography, 54, 799–811.

807 Wienecke BC, Robertson G (1997) Foraging space of emperor penguins Aptenodytes forsteri
808 in antarctic shelf water in winter. Marine Ecology Progress Series, 159.

809 Wolf SG, Snyder MA, Sydeman WJ, Doak DF, Croll DA (2010) Predicting population
810 consequences of ocean climate change for an ecosystem sentinel, the seabird cassin’s auklet.
811 Global Change Biology, 16, 1923–1935.

812 Zimmer I, Wilson R, Gilbert C, Beaulieu M, Ancel a, Plotz J (2008) Foraging movements of
813 emperor penguins at pointe géologie, antarctica. Polar Biology, 31, 229–243.

34
814 7 Supporting Information legends
815 Appendix 1 discusses the effect of density dependence.
816 Appendix 2 describes the selection of sea ice variables and details the principal component
817 analysis.
818 Appendix 3 details vital rates estimation (model selection, estimates and 95% confidence
819 intervals) for: (2.1) breeding success, (2.2) adult survival, and (2.3) probabilities of return
820 to the colony.
821 Appendix 4 details the effect of sea ice variability on demography. Figures show the dis-
822 tributions of annual SICa , breeding success, and male and female adult survival for females
823 and males, along with the resulting distribution of the deterministic growth rate λ. This
824 deterministic rate can be thought of as approximating the growth of the population between
825 time t and t + 1, although this is not always true (see Appendix).
826 Appendix 5 shows the sea ice projected by each climate model for each season of the
827 penguin’s life cycle and the climate selection procedure (section 4.1). It also details our
828 novel approach to obtain stochastic sea ice forecasts from single climate output (section
829 4.2).
830 Appendix 6 describes the projections of vital rates in the future and shows that the range
831 of variation in the forecast vital rates from 2010 to 2100 is plausible.

35
832 8 Tables

36
Table 1: Results of model selection for annual adult survival, including linear or quadratic
effects of SICa at each of four seasons: non-breeding, laying, incubating and rearing. The
first and second columns give the effects included in the model for females and males re-
spectively. The number of parameters in the model is K, ∆AIC is the difference in Akaikes
information criterion (AIC) between each model and the model with the smallest AIC (i.e.
best supported by the data). AIC weights represent the relative likelihood of a model and
are used to create the averaged model; only models for which the cumulative sum of AIC
weights is 0.98 are included.

Effect for females Effect for males K ∆AIC AIC weight


incubating incubating 135 0.00 0.14
rearing2 rearing2 137 0.78 0.09
rearing2 incubating 136 1.17 0.08
incubating incubating2 136 1.25 0.07
rearing2 time-invariant 135 1.59 0.06
rearing2 incubating2 137 2.11 0.05
rearing time-invariant 134 2.22 0.04
rearing incubating 135 2.46 0.04
rearing laying 135 2.93 0.03
rearing2 rearing 136 3.53 0.02
time-invariant time-invariant 133 3.63 0.02
rearing rearing2 136 3.71 0.02
incubating time-invariant 134 3.73 0.02
rearing non breeding 135 3.94 0.02
rearing2 laying 136 4.00 0.02
rearing rearing 135 4.12 0.02
incubating rearing2 136 4.13 0.02
incubating rearing 135 4.14 0.02
rearing incubating2 136 4.15 0.02
rearing incubating2 136 4.15 0.02
rearing2 non breeding 136 2.97 0.03
time-invariant incubating 134 4.46 0.01
laying incubating 135 4.52 0.01
incubating2 incubating2 137 4.62 0.01
incubating laying 135 4.81 0.01
time-invariant non breeding 134 4.99 0.01
laying time-invariant 134 5.25 0.01
time-invariant rearing 134 5.25 0.01
incubating non breeding 135 5.49 0.01
incubating rearing2 136 5.50 0.01
laying incubating2 136 5.75 0.01
time-invariant incubating2 135 5.88 0.01
non breeding time-invariant 134 6.02 0.01
incubating2 rearing2 137 6.10 0.01
time-invariant laying 134 6.31 0.01
laying rearing2 136 6.34 0.01
37
Table 2: Selection of General Circulation Models for climate. Each model was evaluated by
comparing the statistical properties of its sea ice output to observed data from 1979–2007.
Agreement is indicated by an x; the GCMs selected are shown in bold.

Model Non-breeding Laying Incubating Rearing


bccr-bcm2-0
cccma-cgcm3-1 x x x x
cccma-cgcm3-1-t63 x x x x
cnrm-cm3 x x
csiro-mk3-0 x x
gfdl-cm2-0 x x
gfdl-cm2-1 x
giss-aom x x x
giss-model-e-r
iap-fgoals1-0-g
inmcm3-0
ipsl-cm4 x x
miroc3-2-hires x x x
miroc3-2-medres x x x
miub-echo-g
mpi-echam5 x x x x
mri-cgcm2-3-2a x
ncar-ccsm3-0 x x x
ukmo-hadcm3 x x x x
ukmo-hadgem1 x x x x

38
Table 3: Probability that the emperor penguin population in Terre Adélie will decline by more
than 90% from 2010 to 2040, 2060, 2080, and 2100, when sea ice follows the forecasts of each of
the five climate models selected.

models 2040 2060 2080 2100


cccma-cgcm3-1 0.0168 0.2366 0.7625 0.9903
cccma-cgcm3-1-t63 0 0.0205 0.6783 0.9997
ukmo-hadcm3 0 0 0 0.0001
ukmo-hadgem1 0 0.0001 0.0088 0.1276
mpi-echam5 0 0 0 0.0181
entire set 0 0.0514 0.2899 0.4272

39
Table 4: Probability that the emperor penguin population in Terre Adélie will decline below a
specific number of breeding pairs (threshold) by 2100, when sea ice follows the forecasts of each of
the five climate models selected. Thresholds are based on the minimum and maximum number of
observed breeding pairs (N obs) during specific time periods (3 first columns), or specific numbers
(last two columns).

threshold min(N obs1979−2010 ) max(N obs1979−2010 ) max(N obs1962−2010 )


2303 3482 6236 8000 10000
cccma-cgcm3-1 1.00 1.00 1.00 1.00 1.00
cccma-cgcm3-1-t63 1.00 1.00 1.00 1.00 1.00
ukmo-hadcm3 0.21 0.43 0.77 0.87 0.93
ukmo-hadgem1 0.94 0.98 1.00 1.00 1.00
mpi-echam5 0.75 0.90 0.99 1.00 1.00
entire set 0.75 0.84 0.93 0.96 0.98

40
833 9 Figures
834 Figure 1: Observed proportional anomalies in sea ice concentration (SICa ) relative to
835 the mean from 1979 to 2007, for each of four seasons of the penguin life cycle. The grey line
836 shows SICa = 0 and represents the mean SIC from 1979 to 2007.
837 Figure 2: (a) Seasonal SIC anomalies (color lines, see legend) as a function of the annual
838 SICa , with the seasonal differences in SICa set to zero. (b) Seasonal SIC anomalies as a
839 function of the seasonal differences in SICa , with the annual SICa set to zero.
840 Figure 3: Breeding success as a function of proportional anomalies in SIC during the
841 rearing season. The line is a logistic regression fit to the data points shown.
842 Figure 4: Annual adult survival as a function of annual SICa and seasonal differences
843 in SICa . (a c) and (b d) show the survival of females and males respectively. Upper panels
844 (a b) show survival as a function of annual SICa for negative (=-1, grey line), zero (thick
845 black line), and positive (=+1, dotted line) values of seasonal differences in SICa . Lower
846 panels (c d) show survival as a function of seasonal differences in SICa for negative (=-4,
847 grey line), zero (thick black line), and positive (=+4, dotted line) values of annual SICa .
848 Figure 5: The deterministic population growth rate (log λ) as a function of annual SICa
849 and seasonal differences in SICa . The white contours indicate log λ = 0. The color bar
850 shows the values of log λ.
851 Figure 6: The stochastic population growth rate log λs as a function of the mean and
852 variance of sea ice variables. (a) log λs as a function of the mean and variance of annual
853 SICa with seasonal differences in SICa set to zero. The black line stands for a zero variance,
854 the grey line for the observed variance, and the dotted line for a variance equal to twice the
855 observed. (b) log λs as a function of the mean and variance of seasonal differences in SICa ,
856 with annual SICa equal to zero.
857 Figure 7: (a) Smoothed mean of SICa during the four seasons, from each of the five
858 selected GCMs (color lines, see legend). The black line shows the observed SICa from 1979
859 to 2007 and the red line shows the zero value, i.e. represents averaged SIC from 1979 to

41
860 2007.
861 Figure 8: Smoothed standard deviation of SICa during the four seasons, from each of
862 the five selected GCMs (color lines, see legend on figure 7). The red line shows the ob-
863 served standard deviation. The grey box represents 0.5 and 1.5 times the observed standard
864 deviation; these values were used in selecting climate models (Appendix 4).
865 Figure 9: Projections of the emperor penguin population based on SICa forecasts from
866 an ensemble of five GCMs. The black line gives the observed number of breeding pairs
867 from 1979 to 2010. For each GCM, three random population trajectories are shown (thin
868 colored lines), along with the median (thick colored line) and the 95% envelope (grey area),
869 from 40,000 stochastic simulations. The median and 95% envelope are also shown from the
870 combined 200,000 simulations for the set of 5 GCMs.
871 Figure 10: Summary of the projections of the emperor penguin population based on SICa
872 forecasts from an ensemble of five GCMs. The thick colored lines (see legend on figure 9) give
873 the median and the grey area is the 95% envelope from the combined 200,000 simulations for
874 the set of 5 GCMs. Also shown are the probability density functions for simulated population
875 size in 2100, for each GCM.

42
Figure 1: Observed proportional anomalies in sea ice concentration (SIC) relative to the
mean from 1979 to 2007, for each of four seasons of the penguin life cycle. The grey line
shows SICa = 0 and represents the mean SIC from 1979 to 2007.

43
Figure 2: (a) Seasonal SIC anomalies as a function of the annual SICa , with the seasonal
differences in SICa set to zero. (b) Seasonal SIC anomalies as a function of the seasonal
differences in SICa , with the annual SICa set to zero.

44
Figure 3: Breeding success as a function of proportional anomalies in SIC during the rearing
season. The line is a logistic regression fit to the data points shown.

45
Figure 4: Annual adult survival as a function of annual SICa and seasonal differences in
SICa . (a c) and (b d) show the survival of females and males respectively. Upper panels
(a b) show survival as a function of annual SICa for negative (=-1, grey line), zero (thick
black line), and positive (=+1, dotted line) values of seasonal differences in SICa . Lower
panels (c d) show survival as a function of seasonal differences in SICa for negative (=-4,
grey line), zero (thick black line), and positive (=+4, dotted line) values of annual SICa .

46
Figure 5: The deterministic population growth rate (log λ) as a function of annual SICa and
seasonal differences in SICa . The white contours indicate log λ = 0. The color bar shows
the values of log λ

47
Figure 6: The stochastic population growth rate log λs as a function of the mean and variance
of sea ice variables. (a) log λs as a function of the mean and variance of annual SICa with
seasonal differences in SICa set to zero. The black line stands for a zero variance, the grey
line for the observed variance, and the dotted line for a variance equal to twice the observed.
(b) log λs as a function of the mean and variance of seasonal differences in SICa , with annual
SICa equal to zero.

48
Figure 7: Smoothed mean of SICa during the four seasons, from each of the five selected
GCMs (color lines, see legend ). The black line shows the observed SICa from 1979 to 2007
and the red line shows the zero value, i.e. represents averaged SIC from 1979 to 2007.

49
Figure 8: Smoothed standard deviation of SICa during the four seasons, from each of the five
selected GCMs (color lines, see legend on figure 7). The red line shows the observed standard
deviation. The grey box represents 0.5 and 1.5 times the observed standard deviation; these
values were used in selecting climate models (Appendix 4).

50
Figure 9: Projections of the emperor penguin population based on SICa forecasts from an
ensemble of five GCMs. The black line gives the observed number of breeding pairs from 1979
to 2010. For each GCM, three random population trajectories are shown (thin colored lines),
along with the median (thick colored line) and the 95% envelope (grey area), from 40,000
stochastic simulations. The median and 95% envelope are also shown from the combined
200,000 simulations for the set of 5 GCMs.

51
Figure 10: Summary of the projections of the emperor penguin population based on SICa
forecasts from an ensemble of five GCMs. The thick colored lines (see legend on figure 9) give
the median and the grey area is the 95% envelope from the combined 200,000 simulations for
the set of 5 GCMs. Also shown are the probability density functions for simulated population
size in 2100, for each GCM.

52

You might also like