Matthieu Authier
email: authier@cebc.cnrs.fr

Céline Martin

Aurore Ponchon

Stephanie Steelandt

Ilham Bentaleb

Christophe Guinet

Breaking the sticks: a hierarchical change-point model for estimating ontogenetic shifts with stable isotope data

1. Stable isotopes are increasingly used in ecology to investigate ontogenetic shifts 1 in foraging habitat (via δ 13 C) and in trophic level (via δ 15 N). These shifts are in essence an individual-level phenomenon, requiring repeated measures throughout the life of individuals, that is longitudinal data. Longitudinal data require in turn specifying an appropriate covariance structure. Here we present a hierarchical model to jointly investigate individual ontogenetic shifts in δ 13 C and δ 15 N values.

2. In a Bayesian framework, we used a Cholesky decomposition for estimating a moderately-sized covariance matrix, thereby directly estimating correlations between parameters describing time-series of isotopic measurements. We offer guidelines on how to select the covariance structure.

3. The approach is illustrated with a hierarchical change-point (or broken stick) model applied to a data set collected on Southern Elephant Seals, Mirounga leonina.

Ontogenetic shifts in foraging habitat, following a juvenile and variable stage, were detected and interpreted as fidelity to a foraging strategy; while ontogenetic shifts in trophic level were more likely the result of complete independence from maternal resources followed by a gradual increase in trophic level as seals aged.

4. Specifying both an appropriate covariance and mean structure enabled us to draw strong inferences on the ecology of an elusive marine predator, and has wide applicability for isotopic ecology provided repeated isotopic measurements are available.

Introduction

The use of stable isotopes in ecology is expanding rapidly [START_REF] Kelly | Stable Isotopes of Carbon and Nitrogen in the Study of Avian and Mammalian Trophic Ecology[END_REF][START_REF] Newsome | A Niche for Isotopic Ecology[END_REF][START_REF] West | Stable Isotopes as One of Nature's Ecological Recorders[END_REF][START_REF] Wolf | Ten Years of Experimental Animal Isotopic Ecology[END_REF]. This inexpensive technique has become extremely popular to investigate various phenomena, from migration [START_REF] Hobson | Stable Isotope (δ d and δ 13 c) are Geographic Indicators of Natal Origins of Monarch Butterflies in Eastern North America[END_REF] to diet estimation [START_REF] Semmens | Quantifying Inter-and Intra-Population Niche Variability using Bayesian Hierarchical Stable Isotope Mixing Models[END_REF]. A recent application is the detection of temporal shifts in a species' diet trophic level throughout the life of an individual, that is the detection of ontogenetic shifts [START_REF] Estrada | Use of Isotopic Analysis of Vertebrae in Reconstructing Ontogenetic Feeding Ecology in White Sharks[END_REF][START_REF] Post | Individual Variation in the Timing of Ontogenic Niche Shifts in Largemouth Bass[END_REF]. An ontogenic shift is defined as the patterns in an organism's resource use that develop as it increases in size from birth or hatching to its maximum [START_REF] Werner | The Ontogenetic Niche and Species Interactions in Size Structured Populations[END_REF]. In their review on ontogenetic shifts, [START_REF] Werner | The Ontogenetic Niche and Species Interactions in Size Structured Populations[END_REF] focused on changes in habitat use and trophic level, both of which are apprehended in isotopic ecology via carbon (δ 13 C) and nitrogen (δ 15 N) stable isotopes respectively.

Carbon isotopes are used for identifying carbon sources and fluxes within ecosystems [START_REF] Kelly | Stable Isotopes of Carbon and Nitrogen in the Study of Avian and Mammalian Trophic Ecology[END_REF][START_REF] Peterson | Stable Isotopes in Ecosystem Studies[END_REF][START_REF] West | Stable Isotopes as One of Nature's Ecological Recorders[END_REF]. Natural gradients in carbon isotopes occur between terrestrial and marine food webs [START_REF] Schoeninger | Nitrogen and Carbon Isotopic Composition of Bone Collagen from Marine and Terrestrial Animals[END_REF][START_REF] Hobson | Using Stable Isotopes to Determine Seabird Trophic Relationships[END_REF], between inshore and offshore waters (Rau et al., 1982;[START_REF] Hobson | Using Stable Isotopes to Determine Seabird Trophic Relationships[END_REF], between benthic and pelagic foodwebs [START_REF] France | 13 carbon Enrichment in Benthic Compared to Planktonic Algae: Foodweb Implications[END_REF] or between low and high latitudes water masses (Rau et al., 1982[START_REF] Rau | Latitudinal Variation in Plankton δ 13 c: Implications for co 2 and Productivity in Past Oceans[END_REF]. The nitrogen isotopic ratio is a reflection of the trophic level of organisms [START_REF] Post | Using Stable Isotopes to Estimate Trophic Position: Models, Methods and Assumptions[END_REF][START_REF] Vanderklift | Sources of Variation in Consumer-Diet δ 15 n Enrichment: a Meta-Analysis[END_REF]. Because the lighter isotope is usually more reactive, 14 N is preferentially excreted and the heavier 15 N is preferentially retained, a phenomenon known as fractionation [START_REF] Fry | Stable Isotope Ecology[END_REF]. This differential reactivity results in a predictable enrichment of the ratio of 15 N to 14 N from preys to consumers [START_REF] Kelly | Stable Isotopes of Carbon and Nitrogen in the Study of Avian and Mammalian Trophic Ecology[END_REF].

A large number of studies looking at ontogenic shifts concerns species with "cryptic lifestages", in particular marine organisms such as turtles [START_REF] Reich | The "Lost Years" of Green Turtles: using Stable Isotopes to Study Cryptic Lifestages[END_REF], fish [START_REF] Estrada | Use of Isotopic Analysis of Vertebrae in Reconstructing Ontogenetic Feeding Ecology in White Sharks[END_REF][START_REF] Post | Individual Variation in the Timing of Ontogenic Niche Shifts in Largemouth Bass[END_REF] or marine mammals [START_REF] Drago | Ontogenic Dietary Changes in South American Sea Lions[END_REF][START_REF] Hobson | Stable Isotope Analyses of Tooth Annuli Reveal Temporal Dietary Records: an Example Using Steller Sea Lions[END_REF][START_REF] Mendes | Stable Carbon and Nitrogen Isotope Ratio Profiling of Sperm Whale Teeth Reveals Ontogenetic Movements and Trophic Ecology[END_REF][START_REF] Newsome | Retrospective Characterization of Ontogenic Shifts in Killer Whale Diets via δ 13 c and δ 15 n Analysis of Teeth[END_REF]. In some studies, repeated isotopic measurements were available for the same individual using so-called archive tissues, because they are metabolically inert after synthesis, such as vertebrae [START_REF] Estrada | Use of Isotopic Analysis of Vertebrae in Reconstructing Ontogenetic Feeding Ecology in White Sharks[END_REF], or teeth [START_REF] Hobson | Stable Isotope Analyses of Tooth Annuli Reveal Temporal Dietary Records: an Example Using Steller Sea Lions[END_REF][START_REF] Mendes | Stable Carbon and Nitrogen Isotope Ratio Profiling of Sperm Whale Teeth Reveals Ontogenetic Movements and Trophic Ecology[END_REF][START_REF] Newsome | Retrospective Characterization of Ontogenic Shifts in Killer Whale Diets via δ 13 c and δ 15 n Analysis of Teeth[END_REF]. These studies addressed the estimation of a change-point in the time-series of isotopic measurements, yet they typically pooled data from all individuals to infer a population-level change-point, or ontogenetic shift. For example, [START_REF] Newsome | Retrospective Characterization of Ontogenic Shifts in Killer Whale Diets via δ 13 c and δ 15 n Analysis of Teeth[END_REF] fitted a 4 parameters logistic model to estimate a change in dentin δ 15 N of Californian Killer Whales (Orcinus orca) after weaning. The model is fit at the population level, that is assuming all individuals experienced an ontogenetic shifts at the same age, despite apparent individual heterogeneity in the raw plot (their Figure 2a). Ignoring individual heterogeneity when it is in fact present may hinder our ability to draw accurate inferences [START_REF] Cooch | Occam's Shadow: Levels of Analysis in Evolutionary-Ecology -Where to Next[END_REF][START_REF] Petrovskii | Variation in Individual Walking Behavior Creates the Impression of a Lévy Flight[END_REF]. In addition, the change-point is often treated as known even when it was first estimated from the same data. Unless a profile likelihood approach is used, no confidence interval for the change-point is usually reported, and all subsequent inferences are conditional on the point estimate for the change-point.

Stable isotopes in ecology of wild animals are often hailed as a powerful technique. Yet, inferences are typically drawn from statistical analyses that tend to 1) emphasize testing over estimation and goodness-of-fit [START_REF] Graham | Statistical Significance versus Fit: Estimating the Importance of Individual Factors in Ecological Analysis of Variance[END_REF][START_REF] Martìnez Abraìn | Statistical Significance and Biological Relevance: a Call for a more Cautious Interpretation of Results in Ecology[END_REF]; and 2) focus on the mean response at the expense of variability (but see [START_REF] Hénaux | Tracking Large Carnivore Dispersal using Isotopic Clues in Claws: an Application to Cougars across the Great Plains[END_REF]). In the case of detecting an ontogenic shift, the problem is clearly one of estimation: when does an organism change its habitat use or trophic level? Further questions may arise as to what are the ecological, life-history and ultimately population consequences of such an individual change [START_REF] Werner | The Ontogenetic Niche and Species Interactions in Size Structured Populations[END_REF][START_REF] Graham | A Rapid Ontogenetic Shift in the Diet of Juvenile Yellowfin Tuna from Hawaii[END_REF]. This paper thus deals with the problem of estimating individual ontogenic shifts with longitudinal isotopic data, that is repeated measurements of δ 13 C and δ 15 N on the same organism throughout its life. We present a Bayesian change-point model to jointly estimate individual ontogenic shifts in δ 13 C and δ 15 N.

Our aim is to bring forward to a larger audience the vast literature on change-point models [START_REF] Beckage | Bayesian Change-Point Analyses in Ecology[END_REF][START_REF] Hall | A Change Point Model for Estimating the Onset of Cognitive Decline in Preclinical Alzheimer's Disease[END_REF][START_REF] Muniz-Terrera | Random Change Point Models: Investigating Cognitive Decline in the Presence of Missing Data[END_REF][START_REF] Ghosh | Random Changepoint Modelling of HIV Immunologic Responses[END_REF], and how to fit them using the BUGS language [START_REF] Lunn | WinBUGS -a Bayesian Modelling Framework: Concept, Structure, and Extensibility[END_REF].

Change-point, or broken-stick, models aim at finding an abrupt rupture in a time-series. The time-series is assumed to be the juxtaposition of piece-wise linear homogeneous segments, each segment separated from the next by a change-point. Such models have been used in epidemiology to infer the onset of cognitive decline [START_REF] Hall | A Change Point Model for Estimating the Onset of Cognitive Decline in Preclinical Alzheimer's Disease[END_REF][START_REF] Muniz-Terrera | Random Change Point Models: Investigating Cognitive Decline in the Presence of Missing Data[END_REF], of prostate cancer [START_REF] Bellera | Hierarchical Changepoint Models for Biochemical Markers Illustrated by Tracking Postradiotherapy Prostate-Specific Antigen Series in Men with Prostate Cancer[END_REF] or of HIV immunologic response decline [START_REF] Ghosh | Random Changepoint Modelling of HIV Immunologic Responses[END_REF]. In ecology, [START_REF] Beckage | Bayesian Change-Point Analyses in Ecology[END_REF] used a change-point model to study allometric relationships between tree height and tree diameter or to study seedling recruitment with respect to canopy cover along a transect; while [START_REF] Da-Silva | A Bayesian Change-Point Model for Describing Partial Semelparity of a Neotropical Didelphid Marsupial[END_REF] studied post-reproductive survival in a partially semelparous marsupial. These models are very flexible as they allow specifying different probability distributions to describe different parts of a time series. Change-point models thus seem appropriate to describe ontogenetic shifts (e.g. [START_REF] Post | Individual Variation in the Timing of Ontogenic Niche Shifts in Largemouth Bass[END_REF]), but are not prescriptive. Other models (for example [START_REF] Newsome | Retrospective Characterization of Ontogenic Shifts in Killer Whale Diets via δ 13 c and δ 15 n Analysis of Teeth[END_REF]) may prove useful when investigating ontogenic shifts. Our aims here are to expose the use of powerful statistical tools to help ecologists drawing strong inferences [START_REF] Platt | Strong Inference[END_REF]. We will illustrate our methodology with an example using data on Southern Elephant Seals Mirounga leonina.

Southern Elephant Seal Example

Southern Elephant Seals are marine carnivores with a very elusive lifestyle since they can spend more than 80% of their lifetime at sea [START_REF] Mcintyre | A Lifetime at Depth: Vertical Distribution of Southern Elephant Seals in the Water Column[END_REF]. Where they are foraging remained a mystery until the advent of miniaturized electronic tags [START_REF] Biuw | Variations in Behaviour and Condition of a Southern Ocean Top Predator in Relation to in situ Oceanographic Conditions[END_REF]. Seals from îles Kerguelen (49°30' S, 69°30' E) in the Southern Indian Ocean show a dual foraging strategy: animals forage either in Antarctic waters or in polar frontal waters [START_REF] Bailleul | Looking at the Unseen: Combining Bio-logging and Stable Isotopes to Reveal a Shift in the Ecological Niche of a Deep-Diving Predator[END_REF]. Across the Southern Ocean, δ 13 C decreases with increasing latitude [START_REF] Bentaleb | Carbon Isotopic Fractionation by Phytoplankton in the Southern Indian Ocean: Relationship between δ 13 C of Particulate Organic Carbon and Dissolved Carbon Dioxide[END_REF][START_REF] Trull | Insight into Southern Ocean Carbon Export from the δ 13 c of Particles and Dissolved Inorganic Carbon using the SOIREE Iron Release Experiment[END_REF]. Carbon stable isotopes can thus help identify the foraging areas of marine predators: a relative difference of ≈ 2‰ is expected between the two strategies [START_REF] Cherel | Geographical Variation in the Carbon Stable Isotope Signatures of Marine Predators: a Tool to Investigate Their Foraging Areas in the Southern Ocean[END_REF][START_REF] Jaeger | Seabird Satellite Tracking Validates the Use of Latitudinal Isoscapes to Depict Predators'Foraging Areas in the Southern Ocean[END_REF]. Processes underlying carbon isotopic fractionation in marine foodwebs are briefly reviewed in [START_REF] Mackenzie | Locations of Marine Animals Revealed by Carbon Isotopes[END_REF] and a model for fractionation is described in [START_REF] Rau | A Model of Photosynthetic 13 C Fractionation by Marine Phytoplankton based on Diffusive Molecular CO 2 Uptake[END_REF].

With Southern Elephant Seals, we were interested in answering the following questions:

• Are seals faithfull to a foraging strategy [START_REF] Bradshaw | Loyalty Pays: Potential Life History Consequences of Fidelity to Marine Foraging Regions by Southern Elephant Seals[END_REF])?

• When do they become faithfull?

• Are ontogenic shifts in carbon (foraging habitat) and nitrogen (trophic level) isotopes concomitant?

• Are there notable sex differences?

• Can we detect differences in stable isotope values before and after the 1970s population crash [START_REF] Authier | Population Trends of Female Elephant Seals Breeding on the Courbet Peninsula, îles kerguelen[END_REF]?

Material & Methods

Notations and Assumptions

Throughout we will assume the data are composed of N measurements of δ 13 C and δ 15 N on m different individuals. For the j th individual, there are n j measurement, such that N = m ∑ j=1 n j .

These measurement are collected along some biologically-meaningful ordered scale such as age (or size). This scale is assumed continuous for convenience. We will also posit that a piecewise linear, or broken-stick model, provides an adequate description of the data, although this may be relaxed to consider non-linear functions as well. With the broken-stick model, we will denote by K δ 13 C j (K δ 15 N j ) the age of the j th individual when an ontogenetic shift in foraging habitat (trophic level) occurs.

Model Building

The time-series of isotopic measurements for the j th individuals is then modelled as:

for i ∈ [1 : n j ] δ 13 C i, j = a 1, j + (Age i, j -e a 3, j ) ×      a 2, j + ε i,1 , Age i, j ≤ e a 3, j a 4, j + ε i,2 , Age i, j > e a 3, j (1)
where

                                    
a 1, j = isotopic value at ontogenetic shift a 2, j = slope before the ontogenetic shift

a 3, j = log(K δ 13 C j )
a 4, j = slope after the ontogenetic shift ε i,1 ∼ N(0, σ δ 13 C,1 ) are the residuals before the ontogenetic shift ε i,2 ∼ N(0, σ δ 13 C,2 ) are the residuals after the ontogenetic shift and σ δ 13 C is the residual standard deviation, which is allowed to be different before and after the ontogenetic shift. A logarithmic transformation is used to guarantee positive values for all K δ 13 C j or K δ 15 N j . We implicitly assume that only the consumer, not its prey, can experience an isotopic shift, but the model cannot be used to distinguish between these two alternatives [START_REF] Matthews | A Critical Evaluation of Intrapopulation Variation of δ 13 C and Isotopic Evidence of Individual Specialization[END_REF].

The individual coefficients a k ∈ [1:4] , j are assumed to be exchangeable and drawn from a multivariate normal distribution of vector mean α k ∈ [1:4] and covariance matrix of dimension 4:

           a 1 a 2 a 3 a 4            j ∼ MVN                       α 1 α 2 α 3 α 4            ,            σ 2 1 σ 1,2 σ 1,3 σ 1,4 σ 2,1 σ 2 2 σ 2,3 σ 2,4 σ 3,1 σ 3,2 σ 2 3 σ 3,4 σ 4,1 σ 4,2 σ 4,3 σ 2 4                       (2)
This formulation allows to directly estimate correlations between parameter of interest via the covariance matrix. For example, one could be interested to assess whether an ontogenetic shift occurs later or earlier depending on the steepness of the slope a 2, j . The interpretation of such correlations would depend on the biology of the studied organism.

The same broken-stick model can be applied to δ 15 N: this model then calls for the estimation of two independent covariance matrices each of dimension 4: one for δ 13 C and one for δ 15 N (hereafter referred to as 2x4x4). An obvious question is whether ontogenetic shifts in δ 13 C and δ 15 N are simultaneous or correlated. Answering this question requires the estimation of covariance matrix V of dimension 8, as represented on Figure 1 (this model is referred to as 8x8 hereafter).

Specifying the covariance structure of a model has generally received less attention than specifying its mean response, but the problem is no less relevant [START_REF] Pourahmadi | Covariance Estimation: the GLM and Regularization Perspectives[END_REF].

Estimating a covariance matrix of size greater than 2 is a challenge: in addition to the usual restriction to lie between -1 and 1, correlations are jointly constrained. For example, with a 3 × 3 covariance matrix, ρ 1,2 and ρ 1,3 can take any value between -1 and 1, but ρ 2,3 must then conform to the following constraints for the matrix to be positive-definite and invertible [START_REF] Budden | Generating Valid 4 × 4 Correlation Matrices[END_REF]:

ρ 1,2 ρ 1,3 -(1 -ρ 2 1,2 )(1 -ρ 2 1,3 ) ≤ ρ 2,3 ≤ ρ 1,2 ρ 1,3 + (1 -ρ 2 1,2 )(1 -ρ 2 1,3 )
Estimating a matrix such as represented in Figure 1 presents some additional challenges since some elements are constrained to be 0. We opted for a Cholesky decomposition of V into a diagonal matrix Γ and a lower triangular matrix L with 1s on the diagonal:

V = ΓLL T Γ (3)
There are several Cholesky decompositions, all of which guarantee positive-definiteness [START_REF] Pourahmadi | Choleski Decompositions and Estimation of a Covariance Matrix: Orthogonality of Variance-Correlation Parameters[END_REF], but equation 3 neatly separates standard deviation (Γ) and correlation (LL T ) parameters [START_REF] Barnard | Modeling Covariance Matrices in Terms of Standard Deviations and Correlations, with Application to Shrinkage[END_REF][START_REF] Chen | Random Effects Selection in Linear Mixed Models[END_REF]. It becomes possible to force some correlations to be 0 and impose the desired structure for V.

In a Bayesian framework, priors need to be specified on each of the parameters. We used weakly-informative priors: for parameters on the same scale as the data (α 1 , α 2 and α 4 ) we used normal priors with a large variance. For the parameter governing the distribution of ages at ontogenetic shifts, a logarithmic transformation in equation 1 guarantees positive values for all

K δ 13 C j or K δ 15 N j .
For the parameter α 3 , we used a Student-t prior (with location, scale and degrees of freedom set to 0, 10 and 7 respectively [START_REF] Gelman | A Weakly Informative Default Prior Distribution for Logistic and Other Regression Models[END_REF]). For modelling V, we used the priors similar to those of [START_REF] Chen | Random Effects Selection in Linear Mixed Models[END_REF]: independent Half-Normal priors of mean 0 and standard deviation 1.5 for the elements, γ p ∈ [1:8] , of the diagonal matrix Γ, and independent normal priors of mean 0 and standard deviation 0.5 for the elements, λ p ∈[2:8] ,q <p , of L. A prior covariance matrix of dimension 4 (8) with such a specification is depicted on Figure S1 (Figure S3). This prior gives reasonable values (that is between 0 and 10) for the variances of the a i, j , but can be altered depending on the studied organisms. It is also somewhat conservative as most of the probability mass for variance parameters is put on values less than 5.

This prior thus reflects skepticism for large differences between individuals. Uniform priors were put on the residual standard deviations [START_REF] Gelman | Prior Distributions for Variance Parameters in Hierarchical Models (comment on Article by Browne and Draper)[END_REF].

Model Selection

With hierarchical models, model selection is a challenge and several methods have been suggested, such as DIC [START_REF] Spiegelhalter | Bayesian Measures of Model Complexity and Fit (with Discussion)[END_REF][START_REF] Barnett | Using Information Criteria to Select the Correct Variance-Covariance Structure for Longitudinal Data in Ecology[END_REF]; but there is currently no consensus [START_REF] Jordan | What Are the Open Problems in Bayesian Statistics[END_REF]. We choose to avoid using the DIC because of drawbacks such as lack of invariance to reparametrization [START_REF] Spiegelhalter | Bayesian Measures of Model Complexity and Fit (with Discussion)[END_REF] and the following discussion). In fact, DIC was computed but yielded non-sensical results for the estimated number of parameters when the Cholesky decomposition was used (see Table S2). To select an appropriate model, we focused on Posterior Predictive Checks [START_REF] Gelman | Posterior Predictive Assessment of Model Fitness via Realized Discrepancies[END_REF][START_REF] Berkhof | A Bayesian Approach to the Selection and Testing of Mixture Models[END_REF] wherein each fitted model is used to predict (hypothetical) repetitions of the data set. From this hypothetical dataset, we compared an observed summary statistic (T obs ) to its predicted values (T rep ) and compute a p value :

p value = Pr(T rep > T obs ) (4)
A p value close to 0.5 tells us of a good fit (T rep ≈ T obs ), while an extreme p value (0 or 1) betrays a major model misfit. We chose the range of observed isotopic values as discrepancy statistics to assess model fit. The rational for choosing the range as a test statistic is the following: if a change-point is necessary to describe the time-series of isotopic measurement, the range of predicted value is likely to be underestimated when fitting a model with no change-point. The tip of the broken stick will be missed by a simple linear regression, hence an underestimation of the range. Posterior Predictive Checks can be used to test whether a broken-stick model is justified or to select a covariance structure. For example, we can compare the covariance structure depicted in Figure 1 with a simpler structure where the matrix is block diagonal with no correlation between δ 13 C and δ 15 N (that is, ρ 1,5 = ρ 2,6 = ρ 3,7 = ρ 4,8 = 0 in Figure 1).

Checking Model Fit

Once a model has been selected, it is crucial to check model fit [START_REF] Gelman | Philosophy and the Practice of Bayesian Statistics[END_REF].

Therefore model fit was assessed for each individual using a goodness-of-fit statistic for non-linear models [START_REF] Vonesh | Goodness-of-fit in Generalized Nonlinear Mixed-Effects Models[END_REF][START_REF] Huang | Assessing the Goodness of Fit of Forest Models Es-timated by Nonlinear Mixed-Model methods[END_REF]. This concordance coefficient is denoted r c and varies between -1 and 1, with values ≤ 0 betraying a complete lack of fit [START_REF] Vonesh | Goodness-of-fit in Generalized Nonlinear Mixed-Effects Models[END_REF][START_REF] Huang | Assessing the Goodness of Fit of Forest Models Es-timated by Nonlinear Mixed-Model methods[END_REF]. This concordance coefficient assesses the fit of the model at the individual level [START_REF] Huang | Assessing the Goodness of Fit of Forest Models Es-timated by Nonlinear Mixed-Model methods[END_REF], and is computed as follow, with j denoting an individual:

r c j = 1 - n j ∑ i=1 (µ i, j -δ i, j ) 2 n j ∑ i=1 (δ i, j -δ j ) 2 + n j ∑ i=1 (µ i, j -μ j ) 2 + n j ( δ j -μ j ) 2 (5)
where

                         µ i, j = a 1, j + (Age i, j -K j ) ×      a 2, j , Age i, j ≤ K j a 4, j , Age i, j > K j δ j = E(δ i, j ) = n j ∑ i=1 δ i, j n j μ j = E(µ i, j ) = n j ∑ i=1 µ i, j
n j δ j and μ j are the means of the observed and fitted values respectively, while the numerator in equation 5 is the sum of squared-residuals ε i for the j th individual. In the next section, we will apply the above methodology to a "real-life" case.

Southern Elephant Seal Data

Teeth were collected from elephant seals that died of natural causes on îles Kerguelen. Canines grow continuously throughout the whole life without closing of the pulp cavity, allowing for age determination [START_REF] Laws | A New Method of Age Determination for Mammals[END_REF][START_REF] Laws | Antarctic Seals. Research Methods and Techniques[END_REF]. Canines from 47 males and 20 females were analyzed and sampled for isotopic analysis. 18 teeth were sampled on animals that died before a population crash in the 1970s, while the remaining 49 were sampled in the 2000s, after the population had stabilized [START_REF] Authier | Population Trends of Female Elephant Seals Breeding on the Courbet Peninsula, îles kerguelen[END_REF].

Each tooth was cut longitudinally and observed under diffused light to count growth layers. The alternate pattern of two opaque and two translucent growth layers corresponds to the annual biological cycle of Southern Elephant Seals [START_REF] Laws | A New Method of Age Determination for Mammals[END_REF]. Translucent bands are enriched in vitamin D and synthesized when seals are ashore to breed and to moult, while opaque ones are synthesized when at sea [START_REF] Wilske | Seasonal Variation in Vitamin D Metabolites in Southern Elephant Seals (Mirounga leonina) Females at South Georgia[END_REF]. Within a year, a Southern Elephant Seal comes onshore to breed, returns to the sea, then comes onshore to moult before another trip at sea. Thus each growth layer was assumed to correspond to one forth of a year [START_REF] Martin | Stable Carbon and Nitrogen Isotope Variations in Canine Dentin Growth Layers of Kerguelen Southern Elephant Seals[END_REF]. Each growth layer was sampled for 1 mg of bulk dentin using a Micromill T M sampler (ISEM, Université de Montpellier 2). Organic matter δ 13 C and δ 15 N signatures of the bulk dentine were measured with an elemental analyzer (EA-IRMS, Euro-Vector EA 3000) coupled to a continuous flow mass spectrometer (Optima-Micromass) at the Université de Montpellier 2.

As a recent study raised concerns about non-linear offsets of organic %C, %N and C N after acid treatment [START_REF] Brodie | Evidence for Bias in C and N Concentrations and δ 13 c Composition of Terrestrial and Aquatic Organic Materials Due to Pre-Analysis Acid Preparation Methods[END_REF], we forwent any acid (or demineralization) treatment prior to isotopic measurement. As a result, the measured δ 13 C is a mixture of organic carbon with a small amount of inorganic carbon. To test the impact of the inorganic fraction, Martin et al.

(2011) compared acid-treated and untreated samples but found no differences (±0.02‰). [START_REF] Schulting | Stable Carbon and Nitrogen Isotope Analysis on Human Remains from the Early Mesolithic Site of La Vergne (Charente-Maritime, france)[END_REF] found similar C N ratios between bulk dentin and collagen, with a lower carbon and nitrogen contents in bulk dentin most likely due to the mineral fraction. Here we assumed that the impact of the mineral fraction is negligible. If not, relative trends (see Results) should be unaffected under the assumption of a systematic bias.

Stable isotopic signatures are presented in the usual δ notation (in ‰) relative to Pee Dee Belemnite and atmospheric N 2 for δ 13 C and δ 15 N respectively. Typical precisions for isotopic measurement were 0.20 ‰ for both carbon and nitrogen. We used C N ratio thresholds of bone and tooth collagen (2.9 to 3.6) as criteria for the identification of diagenetic alteration [START_REF] Ambrose | Preparation and Characterization of Bone and Tooth Collagen for Isotopic Analysis[END_REF]; assuming that total dentin, whose organic phase is mainly collagen and water [START_REF] Moyes | Composition of the Mineral Phase of Dentin in Southern Elephant Seal and Antarctic Fur Seal Teeth[END_REF], has the same C N ratio than bone and tooth collagen. 1, 590 samples were analyzed, but 176 were discarded because of anomalous C N ratios, yielding a final sample size of 1, 414 (1, 115 from males and 299 from females) analyses from 67 individuals (47 males and 20 females). The first δ 15 N value of each time-series was also removed as it is clearly a reflection of maternal diet [START_REF] Hobson | Stable Isotope Analyses of Tooth Annuli Reveal Temporal Dietary Records: an Example Using Steller Sea Lions[END_REF][START_REF] Martin | Stable Carbon and Nitrogen Isotope Variations in Canine Dentin Growth Layers of Kerguelen Southern Elephant Seals[END_REF]. Summary statistics of the data are available in Table S1 and depicted in Figure S2. It should be stressed that females are under-represented in this data set, and that samples collected from dead females on beaches were biased toward young females. Thus time-series of isotopic measurement were usually shorter for females (Table S1). We fitted the model defined by equation 1 to these data.

To answer questions about any differences between males and females, or between animals living before and after the population crash, we can easily modify the hierarchical change-point model defined by equation 1 by further specifying that the vector of means (α k ∈ [1:4] ) depends on the sex of seals and whether they lived before or after the population crash:

a k ∈[1:4] , j = α 1,k + α 2,k * Sex j + α 3,k * Crash j + η k, j
where the individual-level residuals η k, j are drawn from a multivariate normal distribution of mean 0 and covariance matrix V (see equation 3).

Software

All models were fitted with winBUGS [START_REF] Spiegelhalter | WinBUGS User Manual version 1.4[END_REF] called from R (R Development Core Team, 2009) with the package R2WinBUGS [START_REF] Sturtz | R2winbugs: a Package for Running WinBUGS from R[END_REF]. We used normal priors for regression parameter on the natural scale and Student priors with 7 degrees of freedom [START_REF] Gelman | A Weakly Informative Default Prior Distribution for Logistic and Other Regression Models[END_REF] for regression parameters on the log scale. Three chains were initialized with overdispersed starting values. After appropriate burn-in (200, 000 iterations) and thinning of the chains (1 value every 200 iterations stored), convergence was assessed using the Gelman-Rubin convergence diagnostic [START_REF] Cowles | Markov Chain Monte Carlo Convergence Diagnostics: a Comparative Review[END_REF] with the coda package [START_REF] Plummer | coda: Output Analysis and Diagnostics for MCMC[END_REF]. Posterior mean (or median when posterior distributions were asymmetric) with 95% Highest Probability Density (HPD) intervals are reported as 2.5% Mean 97.5% following [START_REF] Louis | Effective Communication of Standard Error and Confidence Interval[END_REF]. Inferences are based on a posterior sample of 3, 000 iterations. Annotated BUGS code is available in the Appendix, along with an R script and a simulated data set.

Results

Model Selection and Fit

A hierarchical change-point model provided an adequate fit to the elephant seal isotopic data (Figures 2 &3). Ontogenetic shifts in δ 13 C and δ 15 N values were generally supported, except for short time-series and a few individuals. The broken-stick model provided a better fit than a null model with no change-point. The model with the most complex covariance structure (8x8 model) did not greatly improve predictive ability (Table 1). Moreover, the estimated correlations δ 13 C and δ 15 N were small, with a posterior mean of ≈ 0.1 in absolute magnitude (Figure 1). Results from the hierarchical model with no correlation between δ 13 C and δ 15 N are thus reported, although results from the other hierarchical model were very similar. There was no statistical support for distinguishing between sexes or between individuals sampled before or after the population crash (Supplementary Figures 4 &5): the posterior distribution of regression coefficients for both factors was as diffuse as that of its prior and included 0.

Ontogenetic Shifts

Results for the selected hierarchical change point model are summarized in Tables 2 &3. The residual variances for both isotopes were larger before the ontogenetic shift (Table 2). We found individual heterogeneity in all four parameters a k ∈ [1:4] : all variance components were well estimated (Table 3, Supplementary Figure 3). The estimated age at ontogenetic shift was larger for δ 13 C values (3.2 years) than for δ 15 N values (1.9 years, Table 2). This difference was statistically significant at the 5% level. δ 13 C values at ontogenetic shifts were more variable than δ 15 N values, but the variability in age at ontogenetic shift was similar for the two elements (Table 3). There is a sign reversal in slopes before and after the ontogenetic shift in both carbon and nitrogen isotopes (Table 2): the slope was positive and then negative for δ 13 C and the opposite for δ 15 N. Slopes were more variable before than after the ontogenetic shift for both δ 13 C and δ 15 N values (Table 3). There was respectively a small and no correlation between slopes before and after the change-point in δ 13 C and δ 15 N values (Figure 1).

Discussion

Southern Elephant Seal Foraging Ecology

Using as an example the Southern Elephant Seal, a species with a cryptic life-style, we analyzed stable isotope data with a hierarchical change-point model to draw inferences on its foraging habits and its trophic level. Despite the on-going "biologging" revolution, some questions are still not easily addressed with miniaturized tags [START_REF] Hebblewhite | Distinguishing Technology from Biology: a Critical Review of the Use of GPS Telemetry Data in Ecology[END_REF]. For example, equipping a large enough (in the statistical sense) sample of individuals with expensive data recorders that may be lost is usually not an option. For this reason, carbon and nitrogen stable isotopes are no longer studied in ecology as a complementary "side-kick" to biologging, but in their own right [START_REF] Newsome | A Niche for Isotopic Ecology[END_REF][START_REF] Wolf | Ten Years of Experimental Animal Isotopic Ecology[END_REF]. We were interested in inferring the foraging behaviour of Southern Elephant Seal using repeated measurements of dentin δ 13 C and δ 15 N values over the whole life of individuals. Using a hierarchical change-point model, we estimated ontogenetic change-points in both foraging habitats and in trophic level, and found that there was individual variability in both the trajectory and timing of shifts.

Our modelling approach proved fruitfull to investigate some aspects of the ecology of Southern Elephant Seals. In particular, our selected model answered all five questions we asked. After a juvenile stage characterized by a large residual variance, Southern elephant seals became faithfull to a foraging strategy. Inferences drawn from longitudinal isotopic data are in agreement with those of biologging studies [START_REF] Bradshaw | Loyalty Pays: Potential Life History Consequences of Fidelity to Marine Foraging Regions by Southern Elephant Seals[END_REF], but the former involved a larger sample over a longer time-period than the latter. This commitment to a foraging strategy occurred at an early age, on average at about 3 years, but there was substantial individual heterogeneity (Table 3, Figures 2, S6 &S7). An ontogenetic shift in δ 15 N was also detected, but this shift occurred earlier (around 1.9 year-old on average).

The ontogenetic shifts we identified can be the result of several processes, such as complete independence from maternal resources acquired before weaning [START_REF] Hobson | Stable Isotope Analyses of Tooth Annuli Reveal Temporal Dietary Records: an Example Using Steller Sea Lions[END_REF][START_REF] Polischuk | Use of Stable-Carbon and -Nitrogen Isotopes to Assess Weaning and Fasting in Female Polar Bears and their Cubs[END_REF] or a shift in foraging habitat (interfrontal versus Antarctic waters) and trophic level [START_REF] Bailleul | Looking at the Unseen: Combining Bio-logging and Stable Isotopes to Reveal a Shift in the Ecological Niche of a Deep-Diving Predator[END_REF]. If the estimated shift solely resulted from a decay of maternal resources, we would not expect a difference in residual variances before and after a shift. In the case of Southern Elephant Seals, not only residual variances, but also slope variances were larger before the shift (Tables 2 &3). This pattern may be interpreted as an individual switching from a very variable state to a more stable one, or in other words for carbon isotopes, in seals becoming faithfull to a foraging strategy. The posterior mean for the marginal slope after the ontogenetic shift was negative, which we interpreted as individuals foraging in Antarctic waters. These seals have to haul out on îles Kerguelen for reproduction and moulting, and they are very likely to feed on the way [START_REF] Thums | In-Situ Measures of Foraging Success and Prey Encounter Reveals Marine Habitat-Dependent Search Strategies[END_REF], thus diluting a "pure" Antarctic signature for δ 13 C. Hence a negative slope, as the Antarctic signal becomes preponderant over the years. The estimated individual variability showed that some slopes after the shift were null or slightly positive, which can be a reflection of seals foraging always in the same water mass, for example, in pelagic waters of the Polar Front [START_REF] Bailleul | Looking at the Unseen: Combining Bio-logging and Stable Isotopes to Reveal a Shift in the Ecological Niche of a Deep-Diving Predator[END_REF].

Finally, a few individuals had a large positive slope before the shift and a shift late in life. The large positive slope before the shift may be a reflection of seals foraging on the Kerguelen Plateau [START_REF] Bailleul | Looking at the Unseen: Combining Bio-logging and Stable Isotopes to Reveal a Shift in the Ecological Niche of a Deep-Diving Predator[END_REF], which has an enriched δ 13 C signature compared to pelagic water masses [START_REF] Cherel | Geographical Variation in the Carbon Stable Isotope Signatures of Marine Predators: a Tool to Investigate Their Foraging Areas in the Southern Ocean[END_REF]; before switching to an alternative strategy.

Concerning trophic level inferred from δ 15 N values, the shift occurred on average earlier than for the δ 13 C data (Table 2). Slopes before the shift were negative, yet they reversed sign after.

Their magnitude also halved before and after the shift, with very few individual variability left after the shift (Table 3). This pattern suggested the shift in δ 15 N values to mostly reflect the gradual decay of maternal influence on δ 15 N [START_REF] Hobson | Stable Isotope Analyses of Tooth Annuli Reveal Temporal Dietary Records: an Example Using Steller Sea Lions[END_REF], followed by a gradual elevation in the trophic web as seals grew in size. Growth is indeterminate in these seals: they keep growing until their death although growth is very slow in adults [START_REF] Mclaren | Growth in Pinnipeds[END_REF]. This continuous growth means that older seals can physically catch bigger preys, which may explain why we observed a gradual elevation in trophic levels. Additionally, the large energy stores males must accumulate before the breeding season may also drive a shift toward large and energetically profitable preys. Residual variances were also larger before than after the shift but the decrease was not as dramatic as for δ 13 C values (Table 2). Thus this shift may mostly reflect complete independence from maternal inputs.

This pattern of an elevation in trophic level with age (Figure 2) does not conflict with blood isotopic data for males, but was not expected for females: in a previous study, Bailleul et al.

(2010) collected blood samples on juvenile males and on adult females. This study evidenced an elevation in δ 15 N with increasing snout-to-tail length, a proxy for age, only in juvenile males.

This discrepancy probably results from the imbalance of the female data compared to males: few time-series for females spanned more than 4 years (Table S1, Figures S6 &S7). The limited number time-series spanning more than 4 years means that the male pattern largely dominates the population-level pattern in our hierarchical model. Thus blood isotopic data is more reliable to infer the female pattern [START_REF] Bailleul | Looking at the Unseen: Combining Bio-logging and Stable Isotopes to Reveal a Shift in the Ecological Niche of a Deep-Diving Predator[END_REF], although the dentin isotopic analysis suggested that a few females too underwent an elevation in trophic position as they aged (that is, individuals with increasing slope after the ontogenetic shift; Figures 2, S6 & S7).

Modelling strategy

The explicit modelling of correlations between parameters governing a broken-stick model for both δ 13 C and δ 15 N values allowed us to investigate whether ontogenetic shifts in foraging habitat and trophic level were concomitant. There was a very small positive correlation between the ages at shift. The explicit incorporation of this correlation into the model did not substantially improve its predictive ability for δ 13 C or for δ 15 N values (Table 1). There seemed to be such a large variability in individual trajectories of foraging strategy and trophic level in this population that there is no meaningfull 'average' δ 13 C profile associated with an 'average' δ 15 N profile.

Finally, the hierarchical modelling approach enabled us to assess whether there were differences between sexes and between seals living before and after a population crash. The data at hand suggested none (Figures S4 andS5), but the Bayesian framework is explicit about inferences being drawn conditional on the observed data. Thus, failure to detect any differences in this peculiar data set may stem for the imbalance between males and females (respectively 70% versus 30% of seals), and between animals living before and after the population crash (respectively 28% versus 72% of seals).

We believe that the piecewise linear formulation of our change-point model is biologically sound for this species since the change-points reflect life-history events such as complete independence from maternal resources or commitment to a foraging strategy. This assumed model suggested gradual changes after a shift (non-null slopes), which we deemed to be reasonable with longitudinal isotopic data. The interpretation of isotopic data in ecology crucially depends on the rate of tissue turn-over/synthesis, and the accuracy (not the precision) of isotopic data can be quite crude depending on the sampled tissue. Turn-over rates may be very short for some tissues (for example blood plasma), but one order of magnitude larger for others (for example claws) [START_REF] Carleton | Should We Use One-,or Multi-Compartement Models to Describe 13 c Incorporation Into Animal Tissues?[END_REF]. These rates also scale with body mass [START_REF] Carleton | The Effect of Cold-Induced Increased Metabolic Rate on the Rate of 13 c and 15 n Incorporation in House Sparrows (Passer domesticus)[END_REF], which may allow to use experimentally-estimated rates from one species on similar-sized species. However, this is still somewhat of a blackbox for wild animals [START_REF] Wolf | Ten Years of Experimental Animal Isotopic Ecology[END_REF].

Assumptions are unavoidable, but the Bayesian framework is very flexible, allowing to fit models to peculiar data sets rather than "adjusting the data to fit the model". The broken-stick model we assumed reasonable for Southern Elephant Seal need not be so for other species.

With little modification in the prior specification of the covariance matrix, non-linear functional responses such as a logistic curve, which also has 4 parameters, can be easily fitted. However, a logistic curve carries also assumptions such as symmetry and asymptotic isotopic values at the end of the time scale. Finally, the broken-stick model was useful for estimating individual shifts for Southern Elephant Seals, but it did not accommodate cyclic-patterns discernible during the first years in some individuals (Figure S6). The broken-stick model lumped these cycles into a residual variance which was larger in early life compared to late life.

Conclusions

Carbon and nitrogen stable isotope analyses are a powerful technique to peek into the ecology of cryptic species: even a cursory glance at the plethora of studies using this technique cannot fail to notice how often "stable isotopes revealed" biological surprises. The technique is hailed as powerful, which it is even more so conditional on using statistical analyses specifically designed to investigate a particular question (see for example [START_REF] Hénaux | Tracking Large Carnivore Dispersal using Isotopic Clues in Claws: an Application to Cougars across the Great Plains[END_REF]). Here, we presented a hierarchical model to investigate individual patterns of ontogenetic shifts in foraging habitat and trophic level [START_REF] Werner | The Ontogenetic Niche and Species Interactions in Size Structured Populations[END_REF]. The most important aspect of the model is not the specification of the mean response, which can readily be modified to conform to the biology of the studied species, but of the covariance structure. The methodology we outlined can be useful for researchers interested in drawing inferences at the individual level [START_REF] Cooch | Occam's Shadow: Levels of Analysis in Evolutionary-Ecology -Where to Next[END_REF][START_REF] Semmens | Quantifying Inter-and Intra-Population Niche Variability using Bayesian Hierarchical Stable Isotope Mixing Models[END_REF]. Bayesian methods allow to fit with relative ease complex models, and thereby to accommodate the (usually complex) structure of ecological data [START_REF] Ellison | Bayesian Inference in Ecology[END_REF][START_REF] Clark | Why Environmental Scientists are Becoming Bayesians[END_REF]. This move towards Bayesian methods is not confined to ecology [START_REF] Link | Bayesian Inference with Ecological Applications[END_REF][START_REF] O'hara | Bayesian Approaches in Evolutionary Quantitative Genetics[END_REF] or even the biological sciences [START_REF] Treier | Democracy as a Latent Variable[END_REF][START_REF] Wainer | 14 Conversations About Three Things[END_REF]. Rather, it stems for a growing realization that uncertainties need to be quantified and to flow freely across different levels of an analysis to avoid overconfident claims. As more data become available, more complex models can also be fit to refine our knowledge [START_REF] Gelman | Philosophy and the Practice of Bayesian Statistics[END_REF]. The modelling approach outlined here can be further extended to incorporate, for example, a survival analysis [START_REF] Guo | Separate and Joint Modeling of Longitudinal and Event Time Data Using Standard Computer Packages[END_REF][START_REF] Horrocks | Prediction of Pregnancy: a Joint Model for Longitudinal and Binary Data[END_REF][START_REF] Vonesh | Shared Parameter Models for the Joint Analysis of Longitudinal Data and Event Time[END_REF] of Southern Elephant Seals to assess the life-history consequences of a foraging strategy; thereby harnessing the power of stable isotope analyses.

6 Acknowledgements

We would like to thank all volunteers who helped collecting teeth from dead animals found on îles Kerguelen. We thank the Museum National d'Histoire Naturelle (Paris) for kindly providing teeth collected before the population crash. We thank Hubert Vonhof and Els Ufkes for discussing the results. We are also indebted to Christophe Barbraud, Emmanuelle Cam, Luca Börger, two anonymous reviewers and the associate editor for helpful and constructive comments that greatly improved the manuscript. σ ε,1 and σ ε,2 are respectively the residual standard deviations before and after the shift; α 1 and K δ the isotopic value and age at the shift respectively, and α 2 and α 4 the slopes before and after the shift respectively.

8

Figure 1 :

 1 Figure 1: Covariance matrix for a joint broken-stick model of δ 13 C and δ 15 N values. Light gray

Figure 2 :

 2 Figure 2: Broken-stick model fitted to 4 individual time-series of isotopic measurements. Each row corresponds to a different individual. δ 13 C (δ 15 N) profiles are depicted on the left (right)

Figure 3 :

 3 Figure 3: Assessing the fit of the selected model (2x4x4). Distributions of individual-level concordance coefficients, r c are reported for both δ 13 C (x-axis) and δ 15 N (y-axis) values.

Table 1 :

 1 This study is part of a national research program (no. 109, H. Weimerskirch and the observatory Mammifères Explorateurs du Milieu Océanique, MEMO SOERE CTD 02) supported by the French Polar Institute (Institut Paul Emile Victor, IPEV). The Territoire des Terres Australes et Antarctiques Françaises (TAAF),the TOTAL Foundation and ANR-VMC 07 IPSOS-SEAL program contributed to this study. The ethics committee of the French Polar Institute (IPEV) approved this study. All animals in this study were cared for in accordance with its guidelines. Posterior Predictive Checks. The statistic considered is the range of isotopic values and the reported p values are the probability that the predicted range exceeds the observed one. The percentage of individuals with a 0.1 < p value < 0.9 is reported for both carbon and nitrogen isotopic time-series. Broken-stick models decreased the proportion of individuals with extreme p values : a broken-stick model was appropriate for most individuals. There was however little support for an increase in covariance complexity: overall changes in δ 13 C were not correlated with changes in δ 15 N.

	691	7 Tables
		Model δ 13 C δ 15 N
		8x8	0.85 0.73
		2x4x4 0.84 0.69
		Null	0.69 0.58

Table 2 :

 2 Estimated marginals from a broken-stick model fit to the Southern Elephant Seal data.

[START_REF] Phillips | Estimating the Timing of Diet Shifts Using Stable Isotopes[END_REF][START_REF] Popa-Lisseanu | Bats' Conquest of a Formidable Foraging Niche: the Myriads of Nocturnally Migrating Songbirds[END_REF], and more specifically of changes in