spring and autumn phenology curves reveals overwinter mule deer ...

Downloaded from rstb.royalsocietypublishing.org on April 14, 2014

Functional analysis of Normalized Difference Vegetation Index curves reveals overwinter mule deer survival is driven by both spring and autumn phenology Mark A. Hurley, Mark Hebblewhite, Jean-Michel Gaillard, Stéphane Dray, Kyle A. Taylor, W. K. Smith, Pete Zager and Christophe Bonenfant Phil. Trans. R. Soc. B 2014 369, 20130196, published 14 April 2014

Supplementary data

"Data Supplement" http://rstb.royalsocietypublishing.org/content/suppl/2014/04/10/rstb.2013.0196.DC1.ht ml

References

This article cites 88 articles, 5 of which can be accessed free

http://rstb.royalsocietypublishing.org/content/369/1643/20130196.full.html#ref-list-1

This article is free to access

Subject collections

Articles on similar topics can be found in the following collections ecology (526 articles) environmental science (269 articles) plant science (77 articles)

Email alerting service

Receive free email alerts when new articles cite this article - sign up in the box at the top right-hand corner of the article or click here

To subscribe to Phil. Trans. R. Soc. B go to: http://rstb.royalsocietypublishing.org/subscriptions


Functional analysis of Normalized Difference Vegetation Index curves reveals overwinter mule deer survival is driven by both spring and autumn phenology rstb.royalsocietypublishing.org

Research Cite this article: Hurley MA, Hebblewhite M, Gaillard J-M, Dray S, Taylor KA, Smith WK, Zager P, Bonenfant C. 2014 Functional analysis of Normalized Difference Vegetation Index curves reveals overwinter mule deer survival is driven by both spring and autumn phenology. Phil. Trans. R. Soc. B 369: 20130196. http://dx.doi.org/10.1098/rstb.2013.0196 One contribution of 9 to a Theme Issue ‘Satellite remote sensing for biodiversity research and conservation applications’. Subject Areas: ecology, environmental science, plant science Keywords: demography, Normalized Difference Vegetation Index, phenology curve, population dynamics, ungulate, winter severity Author for correspondence: Mark A. Hurley e-mail: [email protected]

Mark A. Hurley1,2, Mark Hebblewhite2,3, Jean-Michel Gaillard4, Ste´phane Dray4, Kyle A. Taylor1,5, W. K. Smith6, Pete Zager7 and Christophe Bonenfant4 1

Idaho Department of Fish and Game, Salmon, ID, USA Wildlife Biology Program, Department of Ecosystem and Conservation Sciences, University of Montana, Missoula, MT, USA 3 Department of Biodiversity and Molecular Ecology, Research and Innovation Centre, Fondazione Edmund Mach, San Michele all’Adige, Trentio, Italy 4 UMR CNRS 5558, Laboratoire Biome´trie et Biologie E´volutive, Universite´ Claude Bernard, Lyon 1, 43 boulevard du 11 novembre 1918, 69622 Villeurbanne Cedex, France 5 Department of Botany, University of Wyoming, Laramie, WY, USA 6 Numerical Terradynamics Simulation Group, Department of Ecosystem and Conservation Sciences, University of Montana, Missoula, MT, USA 7 Idaho Department of Fish and Game, Lewiston, ID, USA 2

Large herbivore populations respond strongly to remotely sensed measures of primary productivity. Whereas most studies in seasonal environments have focused on the effects of spring plant phenology on juvenile survival, recent studies demonstrated that autumn nutrition also plays a crucial role. We tested for both direct and indirect (through body mass) effects of spring and autumn phenology on winter survival of 2315 mule deer fawns across a wide range of environmental conditions in Idaho, USA. We first performed a functional analysis that identified spring and autumn as the key periods for structuring the among-population and among-year variation of primary production (approximated from 1 km Advanced Very High Resolution Radiometer Normalized Difference Vegetation Index (NDVI)) along the growing season. A path analysis showed that early winter precipitation and direct and indirect effects of spring and autumn NDVI functional components accounted for 45% of observed variation in overwinter survival. The effect size of autumn phenology on body mass was about twice that of spring phenology, while direct effects of phenology on survival were similar between spring and autumn. We demonstrate that the effects of plant phenology vary across ecosystems, and that in semi-arid systems, autumn may be more important than spring for overwinter survival.

1. Introduction

Electronic supplementary material is available at http://dx.doi.org/10.1098/rstb.2013.0196 or via http://rstb.royalsocietypublishing.org.

A major challenge for the application of remote sensing to monitoring biodiversity responses to environmental change is connecting remote sensing data to large-scale field ecological data on animal and plant populations and communities [1]. Large herbivores, for example ungulates, are an economically and ecologically important group of species [2] with a global distribution and varied life-history responses to climate that are very sensitive to the timing and duration of plant growing seasons [3]. Until recently, monitoring plant phenology and the nutritional influences on ungulate life histories have been impossible at large spatial scales owing to the intense effort necessary to estimate even localized plant phenology. The remote sensing community has largely solved this issue by partnering with ecologists to

& 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.


2

Phil. Trans. R. Soc. B 369: 20130196

across an entire growing season or discriminates between sites. To fill this important gap, the joint use of functional analysis [51] to characterize seasonal variation in NDVI curves and path analyses [52] to assess both direct and indirect effects of plant phenology offers a powerful way to address entangled relationships of plant quality and their effects on population dynamics of ungulates. Pioneering experimental work on elk (Cervus elaphus) [53] has led to a growing recognition that in temperate areas, late summer and autumn nutrition are important drivers of overwinter survival and demography of large herbivores [53,54]. Summer nutrition first affects adult female body condition [54], which predicts pregnancy rates [53–55], overwinter adult survival rates [54,56], litter size [57] as well as birth mass and early juvenile survival [57–59]. The addition of lactation during summer increases nutritional demand and thus is an important component of the annual nutritional cycle [47,60]. Nutrition during winter (energy) minimizes body fat loss [58], but rarely changes the importance of late summer and autumn nutrition for survival of both juveniles and adults [53]. Winter severity then interacts with body condition to shape winter survival of ungulates [54,61] and can, in severe winters, overwhelm the effect of summer/autumn nutrition through increased energy expenditure, driving overwinter survival of juveniles. Like most other large herbivores of temperate and northern areas, mule deer (Odocoileus hemionus) population growth is more sensitive to change in adult female survival than to equivalent change in other demographic parameters. Survival of adult female mule deer, however, tends to vary little [62,63]; see [64] for a general discussion. By contrast, juvenile survival shows the widest temporal variation, often in response to variation in weather [65–67] and population density [68]. This large variation in juvenile survival, especially over winter, often drives population growth of mule deer [58,62,63]. Fawns accumulate less fat than adults during the summer, which increases their mortality because variation in late summer nutrition interacts with winter severity [62,69]. While previous studies have shown that spring plant phenology correlates with early juvenile survival in ungulates, summer survival is not necessarily more important than overwinter survival. Yet, to date, the effect of changes in autumn plant phenology on overwinter juvenile survival remains unexplored. Our first goal was to identify the annual variation of plant primary production and phenology among mule deer population summer range, measured using NDVI curves of the growing season. Second, with annual plant phenology characterized, we assessed both direct and indirect (through fawn body mass) effects of these key periods on overwinter survival of mule deer fawns. We used a uniquely long-term (1998–2011) and large-scale dataset to disentangle plant phenology effects on mule deer survival, encompassing 13 different populations spread over the entire southern half of Idaho, USA, while most previous studies have focused only within one or two populations. These populations represent diversity of elevations, habitat quality and climatological influences. We focused on overwinter fawn survival because previous studies [62,63] have demonstrated that this parameter is the primary driver of population growth. However, the influences of plant phenology during the growing season and of winter severity on winter survival are not independent because they both involve a strong indirect effect of body mass. Mysterud et al. [70] used a path analysis

rstb.royalsocietypublishing.org

provide circumpolar remotely sensed vegetation indices, fuelling the recent explosion of the integration of remote sensing data into wildlife research and conservation [1,4,5]. With satellites such as the Advanced Very High Resolution Radiometer (AVHRR), the Moderate Resolution Imaging Spectroradiometer (MODIS), the Satellite Pour l’Observation de la Terre (SPOT) [6,7], and growing tool sets for ecologists [8], derived metrics are being commonly used to analyse the ecological processes driving wildlife distribution and abundance [5]. Indices such as the Normalized Difference Vegetation Index (NDVI) and the enhanced vegetation index (EVI) strongly correlate with vegetation productivity, track growing season dynamics [9,10] and differences between landcover types at moderate resolutions over broad spatio-temporal scales [6]. Indices extracted from NDVI correlate with forage quality and quantity [5,11,12] and thus have become invaluable for indexing habitat quality for a variety of ungulates [11,13,14]. For example, only this technology can track a landscape scale plant growth stage that ungulates often select to maximize forage quality [15]. Because of this spatial and temporal link to forage quality, NDVI can be predictive of ungulate nutritional status [11], home range size [16], migration and movements [12,14,17]. An increasing number of studies have also linked NDVI to body mass and demography of a wider array of vertebrates. While there have been recent reviews of the link between NDVI and animal ecology [5], few provided examples where autumn phenology was considered. We conducted a brief review of recent studies to expose readers working at the interface of remote sensing and biodiversity conservation to the preeminent focus on spring phenology using a priori defined variables. From the literature review we performed, 16 out of 22 case studies in temperate areas focused on spring, while three used a growing season average, and only three considered both spring and autumn phenology (table 1). Most studies were based on NDVI metrics describing the active vegetation period, such as start, end and duration of growing season (table 1). Moreover, all but one (see table 1, [40]) was based on a priori defined NDVI metrics assumed to provide a reliable description of plant phenology through the growing season. From this empirical evidence so far reported (see table 1 for details), spring phenology appears as an important period in temperate systems. However, recent field studies on ungulates emphasized the critical importance of late summer and autumn nutritional ecology, suggesting vegetation conditions during this period will also influence population performance of large herbivores. Our brief review complements that of Pettorelli et al. [5] and illustrates the importance of considering phenological dynamics over the entire growing season. Despite this focus on spring phenology, the best existing approach is to use a number of standardized growing season parameters derived from NDVI describing the onset, peak and cessation of plant growth. Unfortunately, these useful parameters are often highly correlated. In Wyoming for example, the start of the growing season was delayed and the rate of green-up was slower than average following winters with high snow cover [50], but these ecologically different processes were highly correlated. Thus, an important barrier to understanding the complex influence of growing season dynamics on ungulate survival is how to disentangle correlated plant phenology metrics. Another underappreciated barrier is the challenge of harnessing the time-series nature of NDVI data, which requires specific statistical tools; no previous study has attempted to describe how the NDVI function varies

juvenile and adult survival

females)

body mass (males and

clutch size

white stork Ciconia ciconia

rustica

barn swallow Hirundo

wing length

tail length

red deer Cervus elaphus

Alectoris rufa

red-legged partridge

africana

African elephant Loxodonta

sheep Ovis aries

minor

lesser grey shrike Lanius

Rangifer tarandus

semi-domesticated reindeer

body mass

parasite loading

survival

(lamb/ewe in December)

reproductive performance

population size

juvenile body mass

population density

juvenile autumn mass

birth mass

body fat

carcass mass

caribou Rangifer tarandus

caribou Rangifer tarandus

body mass

protein mass

NDVI in December –February (eastern and southern Africa)

NDVI in October –November (Sahel)

NDVI in December–February (wintering areas)

monthly NDVI

yearly mean NDVI

seasonal maximum NDVI

NDVI in May

NDVI in March –May

NDVI in January –March (wintering areas)

NDVI in May–June (breeding areas)

summed NDVI over the breeding season

average NDVI in June

average NDVI in June

NDVI metrics

winter ranges only)

adult survival: 0

þ

other traits: 0

clutch size: þ

male wing length, male and female tail length,

other metrics/populations: 0

spring NDVI: þ (Spain only)

þ

juvenile survival: þ

NDVI in March – May: 0

NDVI in May: þ

þ

population density: þ (in populations with poor

juvenile mass: 0

þ

other traits: 0

protein mass: þ

outcome

Phil. Trans. R. Soc. B 369: 20130196

species

[24] [25] [26]

Spain Europe (across Italy (breeding area)

areas)

Poland (breeding

eastern Germany

areas)

Africa (wintering

(Continued.)

[27]

[23]

Kenya

population)

[22]

[21]

[20]

northwestern Patagonia

areas)

Kalahari (wintering

Spain (breeding areas)

France

populations)

Norway (across

(Canada)

Que´bec-Labrador

[19]

[18]

Que´bec-Labrador (Canada)

references

location


trait

Table 1. A brief literature survey of the studies that investigated relationships between NDVI metrics and life-history traits linked to performance and population abundance. The literature survey was performed using ISI web of knowledge using the keywords ‘NDVI and survival’, ‘NDVI and body mass’, ‘NDVI and body weight’, ‘NDVI and reproductive success’, ‘NDVI and recruitment’, ‘NDVI and population growth’ and ‘NDVI and population density’. Only studies performed on vertebrate species were retained. For each case study, the table displays the focal trait(s), the focal species, the NDVI metric(s) used, the outcome (‘þ’: positive association between NDVI and performance; ‘2’: negative association between NDVI and performance; ‘0’: no statistically signiﬁcant association between NDVI and performance), the reference and the location of the study.


3

hystricognath rodents

moose Alces alces

wild boar Sus scrofa

kidney mass

body mass

body mass

pregnancy rate

juvenile body mass

reproductive success

juvenile body mass

body condition

roe deer Capreolus capreolus

elk C. elaphus

reindeer Rangifer tarandus

procyonoides

raccoon dog Nyctereutes

capreolus

roe deer Capreolus

red deer C. elaphus

juvenile body mass

urophasianus

Centrocercus

greater sage grouse

collurio

red-backed shrike Lanius

percnopterus

Egyptian vulture Neophora

africana

African elephant Loxodonta

body mass

juvenile survival


survival

juvenile and adult survival

conception rates

barn swallow Hirundo

adult survival

rustica

species

October)

previous growing season

NDVI correlated with bi-weekly forage biomass and quality over the

productivity)

(onset of spring), mS (rate of EVI increase) and mEVI (plant

EVI modelled using a double logistic function. Use of the parameters S

growing season and rate of NDVI increase in spring)

four NDVI metrics (onset of spring, peak NDVI, summed NDVI over the

summed NDVI over the growing season

seven NDVI metrics (PCA)

yearly NDVI (calculated from monthly NDVI)

summed NDVI in August–October

summed NDVI in April –May

NDVI in the 1st of May

NDVI and Max NDVI in May, June, July and August

0 (Trois Fontaines population)

pregnancy

quality: þ juvenile body mass þ female

exposure to higher predicted forage

success: þ

S and mEVI on both mass and reproductive

rate of NDVI increase: 0

peak NDVI and summed NDVI: þ

onset of spring: 2

wild boar: 0

roe deer: 0

þ

2

þ (Chize´ population)

þ

strong covariation among NDVI metrics

þ (trends only)

NDVI in April (Germany) NDVI in May–August

reproductive success: þ (NDVI in September–

survival: þ (NDVI in December to March)

þ

þ

þ

outcome

NDVI in December –March (South Africa)

NDVI in September–October (Sahel)

NDVI in September –June (breeding areas)

yearly NDVI (wintering areas)

seasonal NDVI (wet versus dry seasons)

NDVI in March –May (wintering areas in Africa)

NDVI in December –February

NDVI in September–November

NDVI metrics

Phil. Trans. R. Soc. B 369: 20130196

trait

[28]

[29,30] [31] [32]

[33]

Denmark

Kenya Spain Germany

Idaho

[36]

France South America (across

[38]

[39]

[40]

[14]

Poland

Finland

Norway

Canada

(Continued.)

[37]

Norway

species

[34] [35]

Norway

Utah (USA)

references

location


Table 1. (Continued.)


4

population rate of increase

population density

survival


population abundance

alpine ibex Capra ibex

kangaroos Macropus sp.

azarae

murine rodent Akodon

virginianus

white-tailed deer Odocoileus

Setophage ruticilla

American redstarts

depending on the breed

NDVI for six and 12 months

seasonal NDVI

effects on survival: 0

þ (but not better predictor than rainfall)

þ

rate of NDVI change and maximum change on reproductive success: 2

success: þ

maximum change between May and June

summed NDVI in May–August on reproductive

þ

other metrics: 0

rate of NDVI change: 2

þ (summed NDVI in March the best predictor)

þ

summed NDVI in June –August: 0 or –

NDVI in late May: þ

outcome

rate of NDVI change between May and June

summed NDVI in May–August

NDVI in December–March (wintering areas)

rate of NDVI change

bighorn sheep Ovis

canadensis

summed NDVI over the growing season

summed NDVI in May

maximum NDVI, date of NDVI peak, summed NDVI in March)

ﬁve NDVI metrics (NDVI slope in spring, NDVI maximum slope in spring,

NDVI in December–February (wintering areas in Africa)

summed NDVI in June –August

NDVI in late May

NDVI metrics

americana

mountain goat Oreamnos

juvenile growth

juvenile survival

chamois Rupicapra rupicapra

common swift Apus apus

Delichon urbicum

common house-martin

juvenile body mass

population size

sheep Ovis aries (two

juvenile mass

breeds)

species

Phil. Trans. R. Soc. B 369: 20130196

trait

[43]

[44] [45]

Italy

France Canada

[48] [49]

Australia

[47]

Argentina

(Canada)

Anticosti, Que´bec

(wintering areas)

Carribean – Cuba

areas)

North America (breeding

[46]

[41,42]

Norway

Italy

references

location


Table 1. (Continued.)


5


(a) Study areas The study area spanned approximately 160 000 km2, representing nearly the entire range of climatic conditions and primary productivity of mule deer in Idaho. We focused on 13 populations with winter ranges corresponding to 13 Idaho game management units (GMUs); hereafter, we use GMU synonymous with population (figure 2). There are three main habitat types (called ecotypes hereafter) based on the dominant overstory canopy species on summer range: coniferous forests, shrubsteppe and aspen woodlands. The populations were distributed among the ecotypes (figure 2) with five populations in conifer ecotype (GMUs 32, 33, 36B, 39, 60A), two in shrub-steppe ecotype (GMUs 54, 58) and six in aspen (GMUs 56, 67, 69, 72, 73A, 76). Elevation and topographic gradients within GMUs affect snow depths and temperature in winter, and precipitation and growing season length in the summer, with elevation increasing from the southwest to the northeast. Conifer GMUs ranged in elevation from 1001 to 1928 m, but most were less than 1450 m. Winter precipitation (winter severity) varied widely (from 10 to 371 mm) in coniferous GMUs. Coniferous ecotype summer ranges are dominated by conifer species interspersed with cool season grasslands, sagebrush and understory of forest shrubs. Shrub-steppe GMUs ranged from 1545 to 2105 m, with winter precipitation from 24 to 105 mm. Summer range within shrub-steppe ecotypes was dominated by mesic shrubs (bitterbrush (Purshia tridentata), sagebrush (Artemisia spp.), rabbitbrush (Chrysothamnus spp.), etc.). Aspen ecotype GMUs were located in the east and south with winter use areas ranging from 1582 to 2011 m, with five of the six GMUs above 1700 m with early winter precipitation ranging from 25 to 146 mm. In summer, productive mesic aspen (Populus tremuloides) woodlands were interspersed with mesic shrubs.

(b) Mule deer monitoring We radiocollared mule deer fawns at six months of age in the 13 GMUs (figure 1), resulting in 2315 mule deer fawns from 1998 to 2011. We captured fawns primarily using helicopters to move deer into drive nets [72], but occasionally by helicopter netgun [73] or clover traps [74]. Fawns were physically restrained and blindfolded during processing with an average handling time of less than 6 min. We measured fawn mass to the nearest 0.4 kg with a calibrated spring scale. Collars weighed 320 – 400 g (less than 2% of deer mass) were equipped with mortality sensors and fastened with temporary attachment plates or surgical tubing, allowing the collars to fall off the animals after approximately 8 – 10 months. We monitored between 20 and 34

(c) Defining population ranges of mule deer We used the mule deer winter and summer ranges for each GMU as the main spatial units of analysis, and we extracted NDVI data from summer range and winter weather from winter range for each year from each population. We combined relocation points for all individuals and years in a single study site to estimate a 95% adaptive kernel home range for both summer and winter [76] for mule deer captured within a population. All deer populations were migratory with an average winter range size of 430 km2 and average summer range size of 3360 km2. Migratory periods, 1 April to 1 June and 1 October to 15 November, were excluded from the home range estimates, and remaining animal locations between 1 June and 30 September were used for summer, 1 December to 31 March for winter. Climate and habitat information was then summarized by the aggregate home range of radiocollared deer for winter and summer within each population.

(d) Functional analysis of Normalized Difference Vegetation Index curves We measured growing season phenology for each population-year using 1 km resolution, 7 day composite AVHRR NDVI data obtained from the National Oceanic and Atmospheric Administration (NOAA)-14, -16 and -17 AVHRR, and maintained by the United States Geological Survey (USGS; http://phenology.cr. usgs.gov/index.php) [77]. AVHRR NDVI data extend over the full temporal extent of our mule deer monitoring effort and has been shown to correspond well with MODIS NDVI data [77]. Radiometric sensor anomalies, atmospheric effects and geometric registration accuracies were previously accounted for according to [77]. Further, the data were accompanied by a cloud contamination mask, which was generated using an adaptation of the cloud clearing of AVHRR data (CLAVR) algorithm [76]. We then rescaled the processed data from the USGS 0-200 classification, with 100 corresponding to vegetated/non-vegetated threshold to the standard NDVI scale of 21 to 1. All cloud contaminated pixels were thus removed by applying this previously generated cloud contamination mask, and the resulting data gaps were infilled using a simple temporal interpolation method [10]. Finally, a minimum NDVI threshold value of zero was applied to define periods of little to no photosynthetic activity and filter any pixels containing ice and snow from the analysis. As phenological changes in NDVI only directly represent ungulate forage dynamics in non-forested vegetation types, we extracted NDVI values from only grass and shrub vegetation types (not burned within 5 years), which we characterized using SAGEMAP landcover data (2005 USGS, Forest and Rangeland Ecosystem Science Center, Snake River Field Station, Boise, ID, USA). Masking in this fashion directly parallels nutritional ecology as mule deer are adapted to feeding in open vegetation types and actively select these types during the growing season [78–80]. To encompass the entire growing season for each population-year, but excluding winter anomalies caused by varying snow condition, we restricted NDVI data to

6

Phil. Trans. R. Soc. B 369: 20130196

2. Material and methods

mule deer fawns in each study area for a total of 185 – 253 annually from 1998 to 2011. We monitored fawns with telemetry for mortality from the ground every 2 days between capture and 15 May through 2006, and then once at the first of each month during 2007– 2011. We located missing fawns aerially when not found during ground monitoring. When a mortality signal was detected, we determined cause of death using a standard protocol [75]. In addition, we kept a minimal annual sample of approximately 600 adult females with radiocollars, using the same capture techniques as fawns. We used the composite sample of monthly aerial and mortality locations over the entire study period from these deer to estimate mule deer population ranges.


to separate independent effects of summer versus winter on body mass. We present a novel methodological framework in which we analyse NDVI measurements using functional principal component analysis (FPCA) to discriminate among study areas in Idaho with differing autumn and spring phenology. We then use hierarchical Bayesian path analysis to identify factors of overwinter mule deer survival. Based on previous studies, we expected that plant phenology should be strongly associated with body mass of mule deer at six months of age, and that body mass and winter severity should interact to determine overwinter survival. We expected direct effects of plant phenology on winter survival to be weaker than winter severity because severe conditions may overwhelm nutritional improvements to fawn quality. We also expected early winter severity would affect overwinter fawn survival more than late winter severity [71].


(a)

(b)

harmonic 1

+

A

FPCA component 2 (27 %) spring J A O M J S

N D

+ +

− −

(c) spring

d =1

A M J J A S O ND AM J J A S O ND 0.6

0.6 A M J J A S O ND 0.6

0.2

0.2

0.2

4

3

autumn 1

2 5

A M J J A S O ND 0.6 0.2

A M J J A S O ND 0.6 0.2

Figure 1. Results of FPCA of the typology of NDVI curves in Idaho, USA, from 1998 to 2011, from April (A) to November (N) for each population-year (dot) identifying two key periods, the spring (second FPCA component, the Y-axis) and the autumn components (first FPCA component, X-axis). (a) Variation in NDVI curves among populations and years was best explained by FPCA 1, which explained 48.9% of the variation and characterized primary production from June to October (e.g. summer/autumn). (b) FPCA 2 (Y-axis) characterized primary production in May and June and explained 27% of the seasonal variation. (c) NDVI typology was best characterized by five clusters, shown in different colours, that corresponded to different patterns of spring and autumn primary production, compared to the mean NDVI curve across all of Idaho. For example, typology 5 was characterized by low NDVI intensity in both spring and autumn, typology 3 by high NDVI intensity in both spring and autumn and typology 4 by high NDVI intensity in spring, but low in autumn, etc. 15 March to 15 November. This time period provided a standardized measure of growing season while capturing the variability both within and between populations for comparing curves. We first assessed among population-year variation in NDVI curves to test direct and indirect (i.e. through body mass) effects of changes in plant phenology on overwinter survival of fawns. In most previous studies (see table 1 for a review), ecologists have either used a priori summary statistics of NDVI. Unfortunately, this approach has led to the use of only a few variables

to define the growing season in any ecosystem; thus to more completely assess vegetation phenology, we proposed a new approach to identify the key periods along the NDVI curve. Instead of defining these periods a priori, our approach is based on a multivariate functional analysis of variation in observed NDVI curves. We used FPCA, a type of functional data analysis (FDA) to analyse among-population and among-year variation in NDVI curves. FDA is specifically designed to characterize information

Phil. Trans. R. Soc. B 369: 20130196

−


A

FPCA component 1 (48.9 %) autumn M J J A S O N D A J A O D

7


4

3

2

32

1

36B 58

33 39

54

73A 56

60A 67 69 72

Figure 2. Distribution of the five NDVI typologies shown in figure 1, with corresponding colours (inset) across the 13 mule deer populations (GMUs) in Idaho, USA, from 1998 to 2011. The size of the pie wedge is proportional to the frequency of occurrence of each NDVI typology within that mule deer population. For example, population 56 had all but one population-year occurring in NDVI typology 4 (figure 1) indicating low primary productivity during spring but higher during autumn. in multivariate time series [51]. FPCA techniques are relatively recent [51] and surprisingly rarely used in ecology and remote sensing (but see [81]) even if they offer a very powerful way to analyse temporal ecological data such as NDVI time series. FPCA was applied to NDVI curves to identify spatio-temporal patterns of vegetation changes. While a priori defined metrics estimated from NDVI data have occasionally been analysed using principal components analysis (PCA) [37], standard PCA is not optimal for time-series data. In PCA, weeks would be considered as independent vectors of values, whereas functional PCA (FPCA) explicitly accounts for the chronology of weeks by treating the statistical unit as the individual NDVI curve. This ensures that the patterns identified by FPCA are constrained to be temporal trends within the growing period (i.e. portions of the curve) and not due to few independent NDVI values. FPCA produces eigenvalues (measuring variation explained by each dimension) and principal component scores for sampling units (summarizing similarities among NDVI curves). However, eigenvectors are replaced by eigenfunctions (harmonics) that show the major functional variations associated to each dimension. To facilitate the application of FPCA by ecologists and remote sensing scientists, we have provided in the electronic supplementary materials the data and the full R code (based on the fda package) to reproduce the analysis performed in the paper. As these methods are poorly known in ecology and remote sensing, we have also provided an expanded description of the mathematical theory, but the reader could consult the original books [51,82] for additional information. Sampling units ( population-years) were partitioned using the k-means algorithm applied on the first two principal component scores. We computed the Calinski and Harabasz criterion for partitions between two and 10 groups, and select the optimal number of clusters that maximizes the criterion. We also computed the amount of variation in the first two principal component scores (NDVI curves) that were explained by space (i.e. population) and time (year). This allowed us to understand which source of variation contributed most to differences in growing season dynamics. We then used principal component scores in subsequent analyses as explanatory variables of mule deer fawn mass and survival.

(e) PRISM weather data We characterized winter (1 November to 31 March) weather conditions using 4 km gridded PRISM observations of minimum monthly temperature and total monthly precipitation from 1995 to

(f ) Environmental effects on body mass and overwinter survival of fawns We estimated population- and year-specific estimates of overwinter fawn survival (from 16 December to 1 June) using staggered Kaplan–Meier non-parametric survival models. We then employed path analysis [52] to test the population-level effects of body mass and winter weather, and to tease apart the direct from the indirect effects (through fawn body mass, see figure 3) of key periods of NDVI on overwinter survival. For the path analysis, we transformed our response variable with an empirical logit function [84] because average survival for each population-year is a proportion bounded between 0 and 1 [85]. We used mass of female fawns in December to measure the cohort quality of the birth year [86] and eliminate the effect of sexual size dimorphism [63]. A first, indirect, mechanistic link between environmental conditions early in life and overwinter survival could be that variation plant phenology and nutritional quality affects the body development of fawns, which in turn, drives overwinter survival. An alternative could be that variation in plant phenology is directly related to overwinter survival as a result of the availability and quality of winter forage. Because winter precipitation was recorded in November–December at the same time as the weighing of fawns, we could not test for an indirect effect of winter precipitation through body mass on overwinter survival. Our model included a population effect entered as a random factor on the intercept to account for the repeated measurements of overwinter survival in different years within a population. We used a Bayesian framework to fit the path analyses to our data [87]. We used non-informative normal (mean of 0 and a s.d. of 100) and uniform (range between 0 and 100) priors for the regression coefficients and variance parameters, respectively. Using JAGS [88], we generated 50 000 samples from Monte Carlo Markov chains to build the posterior distributions of estimated parameters after discarding the first 5000 iterations as a burn in. We checked convergence graphically and based on Gelman’s statistics [87]. Estimated parameters were given by computing the mean of the posterior distribution, and the 2.5th and 97.5th percentiles of the distribution provided its 95% credibility interval. We considered a variable as statistically significant if the credibility interval of its posterior distribution excluded 0. We assessed the fit of the model by computing the squared correlation coefficient between observed and predicted values [89]. Finally, to compare the relative effect sizes of the explanatory variables on overwinter survival, we replicated the analyses using standardized coefficients.

3. Results (a) Functional analysis of Normalized Difference Vegetation Index curves FPCA of NDVI data led to the identification of two independent eigenfunctions (hereafter FPCA components), which reflected contrasting patterns of plant phenology in spring and autumn. Both FPCA components corresponded to

Phil. Trans. R. Soc. B 369: 20130196

76

8


5

2011 [83] (available from http://www.prism.oregonstate.edu). Temperature and precipitation data were averaged across the winter range for each population, and then summed (averaged) across months for precipitation (temperature) to produce climate covariates that represented measures of winter severity, respectively. We produced variables for early winter (November–December) and late winter (January–March) for both precipitation and temperature. These variables were highly correlated (r . 0.4); thus we selected the variable with the highest first-order correlation to our response variable, overwinter survival of fawns, as our winter severity index.


survival (6–12 months)

Figure 3. Hierarchical Bayesian path analysis of the effects of spring and autumn growing season functional components (from figure 1) and winter precipitation on mule deer fawn body mass and overwinter survival from 1998 to 2011 in Idaho, USA. This model explained 44.5% of the variation in survival. Beta coefficients and their s.d. are shown, with solid lines indicating the indirect effects of NDVI on survival through their effects on body mass, and dashed lines indicate the direct effects of NDVI on survival. continua of increasing NDVI intensity, in early and late growing seasons, and were used as explanatory variables of overwinter survival of mule deer fawns. The first FPCA component described the late season phenology, after peak value, and accounted for 48.9% of the total variation in NDVI curves. The second FPCA component represented the early season phenology and accounted for approximately half as much variation as the first FPCA component (27%; figure 1). FPCA components can be interpreted as the amount of deviation from the overall average NDVI curve in terms of overall primary productivity at different times within the growing season. For example, high FPCA component 1 scores mean both high primary productivity in open habitats in autumn, but also a longer autumn growing season compared to lower FPCA component 1 scores (figure 1a,c). Similarly, positive values of FPCA component 2 reflect both higher spring primary productivity and early onset of plant growth (e.g. figure 1b,c; type 4 dark green). Combining both continua in a factorial plane allowed us to distinguish five NDVI types of curve in reference to the overall mean trend (figure 1c). For example, NDVI in autumn was close to the average for the NDVI curve type 2 (dark blue, figure 1c), but NDVI in spring was the lowest of all curve types in figure 1c. Conversely, NDVI curve type 3 (light green, figure 1c) has NDVI values above average in both spring and autumn. The NDVI curve type 1 (light blue, figure 1c) has the highest NDVI in autumn, while NDVI curve type 5 (red, figure 1c) had lowest autumn productivity. Generally, a given population displayed one NDVI curve type, with some extreme values belonging to a different type (figure 2, see also the electronic supplemental material, figure S1). Decomposition of the among-population and among-year variance in NDVI curves in fact shows that most (73.8%) of the observed variation in NDVI curves was accounted for by population (i.e. spatial variation), and much less (20.8%) by annual variation within a population, with a high degree of synchrony between populations within a year (only 5.4% of the variation in NDVI curves is unexplained).

The average body mass of female fawns in December was 34.0 kg (s.e. ¼ 2.55). In agreement with our hypothesis, body mass of six-month-old fawns was positively related to NDVI in both spring and autumn (figures 3 and 4). From the estimated standardized regression coefficients, the effect of NDVI in autumn (FPCA component 1) on autumn body mass of fawns (standardized b ¼ 0.694, s.e. ¼ 0.209) was greater than the effect of NDVI in spring (FPCA component 2; standardized b ¼ 0.652, s.e. ¼ 0.206). FPCA component in the autumn explained more variance in body mass than traditional estimates of phenology such as, start, end or peak date of growing season (electronic supplemental material, table S3). The autumn was thus of more importance to the body development of mule deer fawns at the onset of winter than spring (figures 3 and 4). The annual overwinter survival of mule deer fawns averaged 0.55 (s.e. ¼ 0.24, range ¼ 0–0.94) across populations. Our best model accounted for 44.5% of the observed variation in overwinter survival, including the additive effects of autumn body mass of female fawns, early winter precipitation and of spring and autumn NDVI. As expected, when mean body mass reflects the average demographic performance of a given cohort, the annual overwinter survival of fawns was associated positively with the mean cohort body mass in late autumn (figures 3 and 5a). Total precipitation during early winter from November to December (ranging from 11 to 372 mm) was associated with decreased fawn survival (figures 3 and 5b). Once the effect of body mass and winter precipitations were accounted for, spring had negative impacts on the overwinter survival of fawns (figures 3 and 5d), so that survival was lower with higher NDVI during the spring plant growth season. Autumn was not significantly related to overwinter survival beyond the positive effect on body mass. Winter precipitation has the greatest effect size on overwinter survival of fawns (standardized b ¼ 21.138, s.d. ¼ 0.200), followed by spring (standardized b ¼ 20.587, s.d. ¼ 0.217) and autumn (standardized b ¼ 20.369, s.d. ¼ 0.247), while fawn body mass in autumn has the smallest relative effect size (standardized b ¼ 0.350, s.d. ¼ 0.146). The observed relationships between environmental conditions and overwinter survival of fawns differed slightly among populations but differences were not statistically supported (electronic supplementary material, figure S2).

4. Discussion Our results linked variation in observed plant phenology to body mass and survival of juvenile mule deer during winter across populations and years, demonstrating the benefits of connecting remote sensing and biological information to understand consequences of environmental change on biodiversity. We used a new statistical approach to identify plant phenology from NDVI curves encompassing the entire growing season. Previous studies have reported effects of plant phenology on body mass and demographic parameters

9

Phil. Trans. R. Soc. B 369: 20130196

b=

b = 0.350 (0.146)

This suggests that the five NDVI types we identified (figure 1) strongly reflect the distribution of ecotypes and vegetation characteristics among populations (figure 2).

(b) Environmental effects on body mass and overwinter survival of fawns –1. 138 (0.2 00)

body mass

winter precipitation b = –0.369 (0.247)

b = 0.694 (0.209)

b = 0.652 (0.206)

autumn


b = –0.587 (0.217)

spring


(a)

–5

–10 –2

–1 0 PCA axis 1 (late season NDVI)

1

(b)

fawn body mass (in kg)

5

0

–5

–10 –1.5

–1.0 –0.5 0 0.5 PCA axis 2 (early season NDVI)

1.0

Figure 4. Results of hierarchical Bayesian path analysis showing the standardized direct effects of (a) FPCA component 1 from the functional analysis (autumn NDVI) and (b) FPCA component 2 (Spring NDVI) on body mass (kg) of mule deer fawns in Idaho, USA, from 1998 to 2011. in several species of mammals and birds (see table 1 for a review). However, all these studies but one [40] were based on a priori defined metrics mostly focusing on indices of spring phenology; thus spring metrics appear to explain population parameters, but the relative role of late plant growth season has rarely been investigated. Our approach provides a compelling example and motivation for functional analysis of remote-sensing-derived measures of plant growth as a first step to help identify plant phenological periods most affecting population dynamics of animals. Our results emphasized that the relative role of spring versus autumn phenology is unclear for ungulate species adapted to more arid environments. By defining the periods a posteriori, we found that mule deer fawns survived better in populations with higher NDVI during autumn, and thus longer autumn growing seasons. The effect size of autumn NDVI was stronger than the effect size of spring NDVI for predicting six-month-old body mass. Body mass was positively related to overwinter survival, but precipitation during early winter decreased survival with an effect size almost three

Phil. Trans. R. Soc. B 369: 20130196

fawn body mass (in kg)

0

10


5

times as strong as early winter body mass, similar to other studies of winter ungulate survival [63,90,91]. Previous studies on large herbivores reported an effect of the preceding winter conditions when the juvenile was in utero [37,40,70,92] or an effect of spring conditions [37] on body mass. The patterns of variation in NDVI curves translated to spatial variation in plant growth during autumn, and hence mule deer body mass and survival. First, we found almost twice as much variation in the NDVI curves occurred in the autumn (FPCA component 1, figure 1a) compared with spring (FPCA component 2, figure 1a). Thus, plant phenology during the autumn was more variable than spring in our semi-arid system. Second, we found almost three times the variation in NDVI curves was explained by spatial variation among populations in a given year compared with among-year variation. The high proportion of the variance explained among populations indicates that variation among NDVI curves within a population was consistent from year to year and also synchronous between units within a year. These patterns of stronger variation during autumn (versus spring) and among populations (versus among years) contributed to autumn NDVI having double the effect size on body mass, and hence survival. Thus, the most variable period of the growing season (e.g. autumn) had the strongest effect size on mass and survival. These results mirror results from studies of just the spatial variance in survival [93] and suggest that plant phenology may also synchronize population dynamics. With the recent focus on autumn nutrition of elk [53], however, many ungulate managers in North America are focusing increasingly on autumn nutrition. Our results emphasize that, at least for large herbivores, focusing a priori on just one season, spring or autumn, without explicit consideration of the spatio-temporal variation in the entire curve of plant phenology could be misleading. Forage availability for large herbivores varied by vegetation cover type, precipitation and temperature during the growing season [55,94]. Increased rainfall in summer, reflected in increased NDVI in autumn, will promote growth of forbs [94], a highly selected forage for mule deer [94,95], and can promote new growth in autumn germinating annual graminoids (e.g. cheatgrass, Bromus tectorum) and delay senescence, prolonging access to higher quality forage [14]. Increased summer–autumn nutrition improved calf and adult female survival, fecundity rates and age of first reproduction in captive elk [53]. Rainfall during the growing season also increases quality and quantity of winter forage [94], which increases survival of fawns and adult female mule deer [58]. Tollefson et al. [57] showed that summer forage has the greatest impact on mule deer juvenile survival and overall population growth rate in a penned experiment in eastern Washington, USA. In our study area, effects of climate and plant phenology certainly varied across our southeast to northwest gradient (electronic supplementary material), but will require individual-level analyses of single radiocollared mule deer to most clearly separate out local influences on overwinter survival. Therefore, especially in arid or semi-arid systems, we expect that future studies will identify strong signatures of autumn NDVI and climate on demographic parameters of large herbivore populations, similar to our results. One obvious difference between our arid study system and previous studies of NDVI and large herbivores is that NDVI curves were not a classic bell shape. Instead, plants in open-habitats had a left-skewed growth curve, with a rapid green-up in spring, but then a long right tail in the


(a)

(b)

11


2.5 2

juvenile survival

0

2.5

Phil. Trans. R. Soc. B 369: 20130196

juvenile survival

0

2 5.0

4 30

32

34

36

38

0

body mass (kg) (c)

100 200 300 cumulative winter precipitation (mm)

(d) 2.5

2

juvenile survival

juvenile survival

0 0

2

4 2

1 0 PCA axis 1 (late season NDVI)

1

2.5

5.0 1.5

1.0

0.5

0

0.5

1.0

PCA axis 2 (early season NDVI)

Figure 5. Results of hierarchical Bayesian path analysis showing standardized direct effects of (a) body mass (kg), (b) cumulative winter precipitation (in mm) and (c) FPCA component 1 from the functional analysis (autumn NDVI) and (d ) FPCA component 2 (spring NDVI) on the overwinter survival of mule deer fawns in Idaho, USA, from 1998 to 2011. NDVI distribution, and, occasionally, secondary growth peaks in late summer and autumn (e.g. figure 1c). Most other studies that examined NDVI curves found more symmetrical shapes, with a rapid plant green-up and senescence [37,96]. However, Martinez-Jauregui et al. [25] found the classic bell-shaped NDVI curve for Norwegian and Scottish red deer (C. elaphus), but a similarly earlier and flatter NDVI curve in southern Spain. We believe our right-skewed autumn growing season dynamics may be characteristic of arid or semi-arid systems where precipitation and growing seasons cease during summer. Nonetheless, the variability among studies in the shape of the NDVI curves emphasizes the importance of identifying key periods of the growing season a posteriori. One unexpected result from our study was the negative direct effects of spring NDVI on overwinter survival of mule deer fawns, in contrast to the stronger positive effect of both spring and autumn NDVI on body mass, and of body mass on overwinter fawn survival. There could be several competing explanations for this puzzling result. First,

despite the power of path analysis at disentangling complex relationships [52], there could still remain some confounding effects of body mass or winter severity. Although we attempted to control for spatial variability with random effects of study site, there could also be negative covariance between winter severity, which, because spring NDVI is correlated to winter severity of the preceding winter [50], could lead to negative correlation between spring NDVI and subsequent winter severity. The effect of this general relationship may downscale to study site differently if snow depth passes a threshold where few fawns survive regardless of mass, as is the case sporadically in some of our higher elevation study sites [96– 98] that typically display the most productive NDVI curve types. Mysterud & Austrheim [97] provide a very plausible explanation based on how the negative effect of a later spring (axis 2) will increase winter survival through prolonging access to high-quality forage. Alternatively, viability selection operating on mule deer cohorts may explain this pattern [99]. Counterintuitively, if


Mule deer capture and handling methods were approved by IDFG (Animal Care and Use Committee, IDFG Wildlife Health Laboratory) and University of Montana IACUC (protocol no. 02-11MHCFC-031811).

Acknowledgements. J. Unsworth, B. Compton, M. Scott, C. White, J. Shallow and C. McClellan, and M. Elmer provided guidance and logistical support without which this project would not be possible. We thank Idaho Department of Fish and Game Wildlife technicians, biologists and managers for quality data collection and support. We thank S. Running and M. Zhao for valuable discussions about remote sensing, and M. Mitchell, W. Lowe, P. Lukacs, N. Pettorelli, A. Mysterud and one anonymous reviewer for helpful discussion and comments on previous drafts of the paper.

Funding statement. Financial support was provided by Idaho Department of Fish and Game, Federal Aid in Wildlife Restoration grant no. W-160-R-37, NASA grant no. NNX11AO47G, University of Montana, Mule Deer Foundation, Safari Club International, Universite´ Lyon 1, CNRS and Foundation Edmund Mach.

References 1.

2.

3.

4.

Turner W, Spector S, Gardiner N, Fladeland M, Sterling E, Steininger M. 2003 Remote sensing for biodiversity science and conservation. Trends Ecol. Conserv. 18, 306 –314. (doi:10.1016/S0169-5347 (03)00070-3) Gordon IJ, Hester AJ, Festa-Bianchet M. 2004 The management of wild large herbivores to meet economic, conservation and environmental objectives. J. Appl. Ecol. 41, 1021 –1031. (doi:10. 1111/j.0021-8901.2004.00985.x) Senft RL, Coughenour MB, Bailey DW, Rittenhouse LR, Sala OE, Swift DM. 1987 Large herbivore foraging and ecological hierarchies. Bioscience 37, 789. (doi:10.2307/1310545) Pettorelli N, Vik JO, Mysterud A, Gaillard J-M, Tucker CJ, Stenseth NC. 2005 Using the satellite-derived NDVI to assess ecological responses to

5.

6.

7.

environmental change. Trends Ecol. Evol. 20, 503 –510. (doi:10.1016/j.tree.2005.05.011). Pettorelli N, Ryan S, Mueller T, Bunnefeld N, Jedrzejewska BA, Lima M, Kausrud K. 2011 The Normalized Difference Vegetation Index (NDVI): unforeseen successes in animal ecology. Clim. Res. 46, 15 –27. (doi:10.3354/cr00936) Huete A, Didan K, Miura T, Rodriguez EP, Gao X, Ferreira LG. 2002 Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 83, 195–213. (doi:10.1016/S0034-4257(02)00096-2) Running SW, Nemani RR, Heinsch FA, Zhao MS, Reeves M, Hashimoto H. 2004 A continuous satellite-derived measure of global terrestrial primary production. Bioscience 54, 547–560. (doi:10.1641/00063568(2004)054[0547:ACSMOG]2.0.CO;2)

8.

Dodge S et al. 2013 The environmental-data automated track annotation (Env-DATA) system: linking animal tracks with environmental data. Mov. Ecol. 1, 3. (doi:10.1186/2051-3933-1-3) 9. Zhang XY, Friedl MA, Schaaf CB, Strahler AH, Hodges JCF, Gao F, Reed BC, Huete A. 2003 Monitoring vegetation phenology using MODIS. Remote Sens. Environ. 84, 471–475. (doi:10.1016/ S0034-4257(02)00135-9) 10. Zhao MS, Heinsch FA, Nemani RR, Running SW. 2005 Improvements of the MODIS terrestrial gross and net primary production global data set. Remote Sens. Environ. 95, 164– 176. (doi:10.1016/j.rse. 2004.12.011) 11. Hamel S, Garel M, Festa-Bianchet M, Gaillard J-M, Coˆte´ SD. 2009 Spring Normalized Difference Vegetation Index (NDVI) predicts annual variation in

12

Phil. Trans. R. Soc. B 369: 20130196

a novel and powerful approach for studies of the ecological effects of plant phenology, and arose out of the productive collaboration between remote sensing scientists and ecologists. We anticipate the benefits of functional analyses to extend far beyond NDVI, to ecological analyses of variation in the other remotely sensed vegetation indices (e.g. fPAR, EVI), MODIS snow and temperature datasets, and aquatic measures such as sea surface temperature, chlorophyll and other important ecological drivers. In conclusion, in large parts of the world that are semiarid or deserts, our results strongly show that it may not be just spring phenology that matters to ungulate population dynamics. Our new approach using functional analysis of the entire NDVI curve provides a powerful method to identify first key periods within the growing season and then disentangle their respective roles on demographic traits when combined with hierarchical path analysis. Our approach thus allowed us to determine the most likely pathways by which plant growth influenced mule deer overwinter survival of fawns. Finally, and perhaps most importantly, we demonstrated a novel approach to first identify different temporal components of remote sensing datasets that are the key drivers of large-scale population responses, aiding the broad objective of enhancing our ability to monitor responses of biodiversity to environmental change at global scales.


good spring growing conditions enhance summer survival, a large proportion of the cohort will survive until the onset of the winter, including frail [100] individuals that would experience increased mortality during winter [98], and the opposite during harsh springs. As individual early mortality in populations of large herbivores is tightly linked with maternal condition [66], fawns surviving to the winter will be mostly high-quality fawns enjoying high maternal condition. Those fawns would thus be expected to be robust enough to survive winter. Bishop et al. [58] suggested this exact viability selection process for mule deer fawns in Colorado, supporting our interpretation of this counterintuitive spring NDVI effect. Viability selection could also be compounded through the interaction between winter severity and the preponderance of predator-caused mortality in winter [63]. There might also be negative covariance between neonate and overwinter survival [58], driven as we suggest here by different spring and autumn phenology patterns. Regardless, many plausible biological processes exist to explain the effect of early season plant growth on winter survival of fawns. Functional analysis provides a powerful approach to identify the key periods of the growing season from remote sensing data and to assess their differential effects on life-history traits. Our functional analysis applied to year- and population-specific NDVI curves allowed us to identify two distinct components of variation that corresponded closely to contrasting spring and autumn phenology. Of course, many remote sensing studies have used NDVI for decades to examine differences in spring and autumn phenology [6]. Yet, despite the primacy of multivariate approaches in remote sensing, only a few studies have used even standard PCA to examine spatial trends in NDVI [101] or identify NDVI anomalies [102]. Functional analysis allowed us to identify phenological patterns a posteriori and to summarize NDVI curves into only two independent components instead of 5– 12 a priori defined metrics that are strongly correlated (table 1). Moreover, our FPCA axes explained variation similarly or better than pre-defined parameters based on previous studies (e.g. axis 1 versus senescence date; electronic supplementary material, table S3). Functional analysis provides


13.

15.

16.

17.

18.

19.

20.

21.

22.

23.

24.

25.

27.

28.

29.

30.

31.

32.

33.

34.

35.

36.

37. Herfindal I, Sæther B-E, Solberg EJ, Andersen R, Høgda KA. 2006 Population characteristics predict responses in moose body mass to temporal variation in the environment. J. Anim. Ecol. 75, 1110– 1118. (doi:10.1111/j.1365-2656.2006. 01138.x) 38. Mysterud A, Tryjanowski P, Panek M, Pettorelli N, Stenseth NC. 2007 Inter-specific synchrony of two contrasting ungulates: wild boar (Sus scrofa) and roe deer (Capreolus capreolus). Oecologia 151, 232–239. (doi:10.1007/s00442-006-0584-z) 39. Melis C, Herfindal I, Kauhala K, Andersen R, Høgda K-A. 2010 Predicting animal performance through climatic and plant phenology variables: the case of an omnivore hibernating species in Finland. Mammal. Biol. 75, 151–159. (doi:10.1016/j. mambio.2008.12.001) 40. Tveraa T, Stien A, Ba˚rdsen B-J, Fauchald P. 2013 Population densities, vegetation green-up, and plant productivity: impacts on reproductive success and juvenile body mass in reindeer. PLoS ONE 8, e56450. (doi:10.1371/journal.pone.0056450) 41. Nielsen A, Yoccoz NG, Steinheim G, Storvik GO, Rekdal Y, Angeloff M, Pettorelli N, Holand O, Mysterud A. 2012 Are responses of herbivores to environmental variability spatially consistent in alpine ecosystems? Glob. Change Biol. 18, 3050–3062. (doi:10.1111/j. 1365-2486.2012.02733.x) 42. Nielsen A, Steinheim M, Mysterud A. 2013 Do different sheep breeds show equal responses to climate fluctuations?. Basic Appl. Ecol. 14, 137–145. (doi:10.1016/j.baae.2012.12.005) 43. Ambrosini R, Orioli V, Massimino D, Bani L. 2011 Identification of putative wintering areas and ecological determinants of population dynamics of common house-martin (Delichon urbicum) and common swift (Apus apus) breeding in northern Italy. Avian Conserv. Ecol. 6, 3. (doi:10.5751/ACE00439-060103) 44. Garel M, Gaillard J-M, Jullien J-M, Dubray D, Maillard D, Loison A. 2011 Population abundance and early spring conditions determine variation in body mass of juvenile chamois. J. Mammal. 92, 1112– 1117. (doi:10.1644/10-mamm-a-056.1) 45. Pettorelli N, Pelletier F, von Hardenberg A, FestaBianchet M, Coˆte´ SD. 2007 Early onset of vegetation growth vs. rapid green-up: impacts on juvenile mountain ungulates. Ecology 88, 381–390. 46. Wilson S, LaDeau SL, Tøttrup AP, Marra PP. 2011 Range-wide effects of breeding-and nonbreedingseason climate on the abundance of a neotropical migrant songbird. Ecology 92, 1789–1798. (doi:10. 1890/10-1757.1) 47. Simard MA, Coulson T, Gingras A, Coˆte´ SD. 2010 Influence of density and climate on population dynamics of a large herbivore under harsh environmental conditions. J. Wildl. Manage. 74, 1671– 1685. (doi:10.2193/2009-258) 48. Andreo V, Provensal C, Scavuzzo M, Lamfri M, Polop J. 2009 Environmental factors and population fluctuations of Akodon azarae (Muridae: Sigmodontinae) in central Argentina. Austral Ecol. 34, 132–142. (doi:10.1111/j.1442-9993.2008.01889.x)

13

Phil. Trans. R. Soc. B 369: 20130196

14.

26.

Langvatn R, Coulson T. 2009 Are local weather, NDVI and NAO consistent determinants of red deer weight across three contrasting European countries? Glob. Change Biol. 15, 1727–1738. (doi:10.1111/j. 1365-2486.2008.01778.x) Saino N, Sze´p T, Ambrosini R, Romano M, Møller AP. 2004 Ecological conditions during winter affect sexual selection and breeding in a migratory bird. Proc. R. Soc. Lond. B 271, 681–686. (doi:10.1098/ rspb.2003.2656) Schaub M, Kania W, Ko¨ppen U. 2005 Variation of primary production during winter induces synchrony in survival rates in migratory white storks Ciconia ciconia. J. Anim. Ecol. 74, 656 –666. (doi:10.1111/j. 1365-2656.2005.00961.x) Sze´p T, Møller AP, Piper S, Nuttall R, Szabo´ ZD, Pap PL. 2006 Searching for potential wintering and migration areas of a Danish barn swallow population in South Africa by correlating NDVI with survival estimates. J. Ornithol. 147, 245–253. (doi:10.1007/s10336-006-0060-x) Wittemyer G, Barner Rasmussen H, DouglasHamilton I. 2007 Breeding phenology in relation to NDVI variability in free-ranging African elephant. Ecography 30, 42 –50. (doi:10.1111/j.0906-7590. 2007.04900.x) Rasmussen HB, Wittemyer G, Douglas-Hamilton I. 2006 Predicting time-specific changes in demographic processes using remote-sensing data. J. Appl. Ecol. 43, 366–376. (doi:10.1111/j.13652664.2006.01139.x) Grande JM, Serrano D, Tavecchia G, Carrete M, Ceballos O, Diaz-Delgado R, Tella JL, Dona´zar JA. 2009 Survival in a long-lived territorial migrant: effects of life-history traits and ecological conditions in wintering and breeding areas. Oikos 118, 580 –590. (doi:10.1111/j.1600-0706.2008.17218.x) Schaub M, Reichlin TS, Abadi F, Ke´ry M, Jenni L, Arlettaz R. 2012 The demographic drivers of local population dynamics in two rare migratory birds. Oecologia 168, 97 –108. (doi:10.1007/s00442-0112070-5) Guttery MR, Dahlgren DK, Messmer TA, Connelly JW, Reese KP, Terletzky PA, Burkepile N, Koons DN. 2013 Effects of landscape-scale environmental variation on greater sage-grouse chick survival. PLoS ONE 8, e65582. (doi:10.1371/journal.pone.0065582) Pettorelli N, Mysterud A, Yoccoz NG, Langvatn R, Stenseth NC. 2005 Importance of climatological downscaling and plant phenology for red deer in heterogenous environments. Proc. R. Soc. B 272, 2357 –2364. (doi:10.1098/rspb.2005.3218) Pettorelli N, Gaillard J-M, Mysterud A, Duncan P, Stenseth NC, Delorme D, Van Laere G, Toı¨go C, Klein F. 2006 Using a proxy of plant productivity (NDVI) to find key periods for animal performance: the case of roe deer. Oikos 112, 565–572. (doi:10.1111/j. 0030-1299.2006.14447.x) Diaz GB, Ojeda RA, Rezende EL. 2006 Renal morphology, phylogenetic history and desert adaptation of South American hystricognath rodents. Funct. Ecol. 20, 609–620. (doi:10.1111/j. 1365-2435.2006.01144.x)


12.

timing of peak faecal crude protein in mountain ungulates. J. Appl. Ecol. 46, 582 –589. (doi:10. 1111/j.1365-2664.2009.01643.x) Cagnacci F et al. 2011 Partial migration in roe deer: migratory and resident tactics are end points of a behavioural gradient determined by ecological factors. Oikos 120, 1790– 1802. (doi:10.1111/j. 1600-0706.2011.19441.x) Ryan SJ, Cross PC, Winnie J, Hay C, Bowers J, Getz WM. 2012 The utility of Normalized Difference Vegetation Index for predicting African buffalo forage quality. J. Wildl. Manage. 76, 1499– 1508. (doi:10.1002/jwmg.407) Hebblewhite M, Merrill E, McDermid G. 2008 A multi-scale test of the forage maturation hypothesis in a partially migratory ungulate population. Ecol. Monogr. 78, 141–166. (doi:10.1890/06-1708.1) Fryxell JM, Greever J, Sinclair ARE. 1988 Why are migratory ungulates so abundant? Am. Nat. 131, 781–798. Morellet N et al. 2013 Seasonality, weather and climate affect home range size in roe deer across a wide latitudinal gradient within Europe. J. Anim. Ecol. 82, 1326–1339. (doi:10.1111/1365-2656. 12105) Sawyer H, Kauffman MJ. 2011 Stopover ecology of a migratory ungulate. J. Anim. Ecol. 80, 1078 –1087. (doi:10.1111/j.1365-2656.2011.01845.x) Couturier S, Coˆte´ SD, Huot J, Otto RD. 2008 Bodycondition dynamics in a northern ungulate gaining fat in winter. Can. J. Zool. 87, 367–378. (doi:10. 1139/Z09-020) Couturier S, Coˆte´ SD, Otto RD, Weladji RB, Huot J. 2009 Variation in calf body mass in migratory caribou: the role of habitat, climate, and movements. J. Mammal. 90, 442–452. (doi:10. 1644/07-mamm-a-279.1) Tveraa T, Fauchald P, Yoccoz NG, Anker Ims R, Aanes R, Arild Høgda K. 2007 What regulate and limit reindeer populations in Norway? Oikos 116, 706– 715. (doi:10.1111/j.0030-1299.2007.15257.x) Giralt D, Brotons L, Valera F, Krisˇtin A. 2008 The role of natural habitats in agricultural systems for bird conservation: the case of the threatened lesser grey shrike. Biodivers. Conserv. 17, 1997 –2012. (doi:10. 1007/s10531-008-9349-9) Texeira M, Paruelo JM, Jobbagy E. 2008 How do forage availability and climate control sheep reproductive performance? An analysis based on artificial neural networks and remotely sensed data. Ecol. Model. 217, 197– 206. (doi:10.1016/j. ecolmodel.2008.06.027) Wittemyer G. 2011 Effects of economic downturns on mortality of wild African elephants. Conserv. Biol. 25, 1002 –1009. (doi:10.1111/j.1523-1739. 2011.01713.x) Calvete C, Estrada R, Lucientes J, Estrada A, Telletxea I. 2003 Correlates of helminth community in the red-legged partridge (Alectoris rufa L.) in Spain. J. Parasitol. 89, 445– 451. (doi:10.1645/ 0022-3395(2003)089[0445:COHCIT]2.0.CO;2) Martinez-Jauregui M, San Miguel-Ayanz A, Mysterud A, Rodriguez-Vigal C, Clutton-Brock T,


78. Thiel JR. 2012 Forage selection by maternal mule deer: body condition of maternal females, and birth characteristics and survival of neonates. Pocatello, ID: Idaho State University. 79. Hamlin KL, Mackie RJ. 1989 Mule deer in the Missouri River breaks, Montana: a study of population dynamics in a fluctuating environment, p. 401. Bozeman, MT: Montana Fish, Wildlife and Parks. 80. Mackie RJ, Pac DF, Hamlin KL, Ducek GL. 1998 Ecology and management of mule deer and whitetailed deer in Montana, p. 180. Helena, MT: Montana Fish, Wildlife and Parks. 81. Embling CB, Illian J, Armstrong E, van der Kooij J, Sharples J, Camphuysen KCJ, Scott BE. 2012 Investigating fine-scale spatio-temporal predator– prey patterns in dynamic marine ecosystems: a functional data analysis approach. J. Appl. Ecol. 49, 481–492. (doi:10.1111/j.1365-2664.2012.02114.x) 82. Ramsay JO, Hooker G, Graves S. 2009 Functional data analysis with R and MATLAB. New York, NY: Springer. 83. Daly C, Taylor G, Gibson W. 1997 The PRISM approach to mapping precipitation and temperature. In Proc., 10th Conf. on Applied Climatology, Seattle, USA, October 1997, pp. 10 –12. Seattle, WA: American Meteorology Society 84. Warton DI, Hui FKC. 2010 The arcsine is asinine: the analysis of proportions in ecology. Ecology 92, 3–10. (doi:10.1890/10-0340.1) 85. Zar JH. 1995 Biostatistical analysis. London, UK: Prentice Hall. 86. Hamel S, Coˆte´ SD, Gaillard J-M, Festa-Bianchet M. 2009 Individual variation in reproductive costs of reproduction: high-quality females always do better. J. Anim. Ecol. 78, 143–151. (doi:10.1111/j.13652656.2008.01459.x) 87. Gelman A, Hill A. 2007 Data analysis using regression and multilevel/hierarchical models, p. 625. Cambridge, UK: Cambridge University Press. 88. Plummer M. 2003 JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. In Proc. 3rd Int. Workshop on Distributed Statistical Computing (DSC 2003), Vienna, Austria, 20 –22 March 2003. Vienna, Austria: Austrian Association for Statistical Computing and R Foundation for Statistical Computing. 89. Zheng B, Agresti A. 2000 Summarizing the predictive power of a generalized linear model. Stat. Med. 19, 1771–1781. (doi:10.1002/1097-0258 (20000715)19:13,1771::AID-SIM485.3.0.CO;2-P) 90. Bartmann RM. 1984 Estimating mule deer winter mortality in Colorado. J. Wildl. Manage. 48, 262–267. (doi:10.2307/3808485) 91. Bishop CJ, Unsworth JW, Garton EO. 2005 Mule deer survival among adjacent populationsin Southwest Idaho. J. Wildl. Manage. 69, 311–321. (doi:10.2193/ 0022-541x(2005)069,0311:mdsaap.2.0.co;2) 92. Post E, Stenseth NC, Langvatn R, Fromentin J-M. 1997 Global climate change and phenotypic variation among red deer cohorts. Proc. R. Soc. Lond. B 263, 1317– 1324. (doi:10.1098/rspb. 1997.0182)

14

Phil. Trans. R. Soc. B 369: 20130196

63. Hurley MA, Unsworth JW, Zager P, Hebblewhite M, Garton EO, Montgomery DM, Skalski JR, Maycock CL. 2011 Demographic response of mule deer to experimental reduction of coyotes and mountain lions in southeastern Idaho. Wildl. Monogr. 178, 1–33. (doi:10.1002/wmon.4) 64. Gaillard J-M, Yoccoz NG. 2003 Temporal variation in survival of mammals: a case of environmental canalization? Ecology 84, 3294 –3306. (doi:10.1890/ 02-0409) 65. Portier C, Festa-Bianchet M, Gaillard J-M, Jorgenson JT, Yoccoz NG. 1998 Effects of density and weather on survival of bighorn sheep lambs (Ovis canadensis). J. Zool. Lond. 245, 271–278. (doi:10.1111/j.1469-7998.1998.tb00101.x) 66. Gaillard J-M, Festa-Bianchet M, Yoccoz NG, Loison A, Toı¨go C. 2000 Temporal variation in fitness components and population dynamics of large herbivores. Annu. Rev. Ecol. System. 31, 367 –393. (doi:10.1146/annurev.ecolsys.31.1.367) 67. Coulson T, Catchpole EA, Albon SD, Morgan BJT, Pemberton JM, Clutton-Brock TH, Crawley MJ, Grenfell BT. 2001 Age, sex, density, winter weather, and population crashes in soay sheep. Science 292, 1528 –1531. (doi:10.1126/science.292.5521.1528) 68. Bartmann RM, White GC, Carpenter LH. 1992 Compensatory mortality in a Colorado mule deer population. Wildl. Monogr. 121, 3 –39. 69. White GC, Bartmann RM. 1998 Effect of density reduction on overwinter survival of free-ranging mule deer fawns. J. Wildl. Manage. 62, 214–225. (doi:10.2307/3802281) 70. Mysterud A, Yoccoz NG, Langvatn R, Pettorelli N, Stenseth NC. 2008 Hierarchical path analysis of deer responses to direct and indirect effects of climate in northern forest. Phil. Trans. R. Soc. B 363, 2357 –2366. (doi:10.1098/rstb.2007.2206) 71. Hurley MA, Hebblewhite M, Unsworth JW, Zager P, Miyasaki H, Scott MS, White C, Skalski JR. 2010 Survival and population modeling of mule deer. In Federal Aid in Wildlife Restoration Annual Project Report. Boise, ID: Bradly B. Compton. 72. Beasom S, Evans W, Temple L. 1980 The drive net for capturing western big game. J. Wildl. Manage. 44, 478– 480. (doi:10.2307/3807981) 73. Barrett M, Nolan J, Roy L. 1982 Evaluation of hand held net-gun to capture large mammals. Wildl. Soc. Bull. 10, 108–114. 74. Clover MR. 1954 A portable deer trap and catch net. Calif. Fish Game 40, 367 –373. 75. Wade DA, Bowns JE. 1982 Procedures for evaluating predation on livestock and wildlife. College Station, TX: Texas Agricultural Extension Service, Texas Agricultural Experiment Station, Texas A & M University System. 76. Worton BJ. 1989 Kernel methods for estimating the utilization distribution in home-range studies. Ecology 70, 164 –168. (doi:10.2307/1938423) 77. Eldenshink J. 2006 A 16-year time series of 1 km AVHRR satellite data of the conterminous United States and Alaska. Photogramm. Eng. Remote Sens. 72, 1027. (doi:10.14358/PERS.72.9.1027)


49. Pople AR, Grigg GC, Phinn SR, Menke N, McAlpine C, Possingham HP. 2010 Reassessing the spatial and temporal dynamics of kangaroo populations. In Macropods: the biology of kangaroos, wallabies, and rat-kangaroos (eds G Coulson, M Eldridge), pp. 197 –210. Collingwood, Australia: CSIRO. 50. Christianson D, Klaver R, Middleton A, Kauffman M. 2013 Confounded winter and spring phenoclimatology on large herbivore ranges. Landscape Ecol. 28, 427–437. (doi:10.1007/s10980012-9840-2) 51. Ramsay R, Silverman BW. 2005 Functional data analysis. New York, NY: Springer. 52. Shipley B. 2009 Confirmatory path analysis in a generalized multilevel context. Ecology 90, 363–368. (doi:10.1890/08-1034.1) 53. Cook JG, Johnson BK, Cook RC, Riggs RA, Delcurto T, Bryant LD, Irwin LL. 2004 Effects of summer – autumn nutrition and parturition date on reproduction and survival of elk. Wildl. Monogr. 155, 1–61. 54. Monteith KL, Stephenson TR, Bleich VC, Conner MM, Pierce BM, Bowyer RT. 2013 Risk-sensitive allocation in seasonal dynamics of fat and protein reserves in a long-lived mammal. J. Anim. Ecol. 82, 377–388. (doi:10.1111/1365-2656.12016) 55. Stewart KM, Bowyer RT, Dick BL, Johnson BK, Kie JG. 2005 Density-dependent effects on physical condition and reproduction in North American elk: an experimental test. Oecologia 143, 85–93. (doi:10. 1007/s00442-004-1785-y) 56. Bender LC, Lomas LA, Browning J. 2007 Condition, survival, and cause-specific mortality of adult female mule deer in north-central New Mexico. J. Wildl. Manage. 71, 1118– 1124. (doi:10.2193/ 2006-226) 57. Tollefson TN, Shipley LA, Myers WL, Keisler DH, Dasgupta N. 2010 Influence of summer and autumn nutrition on body condition and reproduction in lactating mule deer. J. Wildl. Manage. 74, 974–986. (doi:10.2193/2008-529) 58. Bishop CJ, White GC, Freddy DJ, Watkins BE, Stephenson TR. 2009 Effect of enhanced nutrition on mule deer population rate of change. Wildl. Monogr. 172, 1– 28. (doi:10.2193/2008-107) 59. Lomas LA, Bender LC. 2007 Survival and causespecific mortality of neonatal mule deer fawns, north-central New Mexico. J. Wildl. Manage. 71, 884–894. (doi:10.2193/2006-203) 60. Sadleir RMFS. 1982 Energy consumption and subsequent partitioning in lactating black-tailed deer. Can. J. Zool. 60, 382 –386. (doi:10.1139/ z82-051) 61. Singer FJ, Harting A, Symonds KK, Coughenour MB. 1997 Density dependence, compensation and environmental effects on elk calf mortality in Yellowstone National Park. J. Wildl. Manage. 61, 12 –25. 62. Unsworth JA, PAc DF, White GC, Bartmann RM. 1999 Mule deer survival in Colorado, Idaho, and Montana. J. Wildl. Manage. 63, 315–326. (doi:10. 2307/3802515)


100. Vaupel JW, Manton KG, Stallard E. 1979 The impact of heterogeneity in individual frailty on the dynamics of mortality. Demography 16, 439–454. (doi:10.2307/2061224) 101. Hall-Beyer M. 2003 Comparison of single-year and multiyear NDVI time series principal components in cold temperate biomes. Geosci. Remote Sens., IEEE Trans. 41, 2568–2574. (doi:10.1109/tgrs.2003. 817274) 102. Lasaponara R. 2006 On the use of principal component analysis (PCA) for evaluating interannual vegetation anomalies from SPOT/VEGETATION NDVI temporal series. Ecol. Model. 194, 429– 434. (doi:10.1016/j.ecolmodel.2005.10.035)

15

Phil. Trans. R. Soc. B 369: 20130196

96. Pettorelli N, Pelletier F, von Hardenberg A, FestaBiachet M, Coˆte´ SD. 2007 Early onset of vegetation growth vs. rapid green-up: impacts on juvenile mountain ungulates. Ecology 88, 381 –390. 97. Mysterud A, Austrheim G. In press. Lasting effects of snow accumulation on summer performance of large herbivores in alpine ecosystems may not last. J. Anim. Ecol. (doi:10.1111/1365-2656.12166) 98. Wilson AJ, Nussey DH. 2010 What is individual quality? An evolutionary perspective. Trends Ecol. Evol. 25, 207 –214. (doi:10.1016/j.tree.2009. 10.002) 99. Fisher RA. 1930 The genetical theory of natural selection. Oxford, UK: Oxford University Press.


93. Lukacs PM, White GC, Watkins BE, Kahn RH, Banulis BA, Finley DJ, Holland AA, Martens JA, Vayhinger J. 2009 Separating components of variation in survival of mule deer in Colorado. J. Wildl. Manage. 73, 817–826. (doi:10.2193/2008-480) 94. Marshal JP, Krausman PR, Bleich VC. 2005 Rainfall, temperature, and forage dynamics affect nutritional quality of desert mule deer forage. Rangeland Ecol. Manage. 58, 360 –365. (doi:10.2111/15515028(2005)058[0360:rtafda]2.0.co;2) 95. Hobbs NT, Baker DL, Gill RB. 1983 Comparative nutritional ecology of Montane ungulates during winter. J. Wildl. Manage. 47, 1 –16. (doi:10.2307/ 3808046)