Marl Ène Gamelon
email: marlene.gamelon@ntnu.no

Jar Tufto

Anna L K Nilsson

Kurt Jerstad

Ole W Røstad

Nils C Stenseth

Bernt-Erik Saether

Environmental drivers of varying selective optima in a small passerine: A multivariate, multiepisodic approach

Keywords: laying date, fluctuating environment, selection episode

In changing environments, phenotypic traits are shaped by numerous agents of selection. The optimal phenotypic value maximizing the fitness of an individual thus varies through time and space with various environmental covariates. Selection may differ between different life-cycle stages and act on correlated traits inducing changes in the distribution of several traits simultaneously. Despite increasing interests in environmental sensitivity of phenotypic selection, estimating varying selective optima on various traits throughout the life cycle, while considering (a)biotic factors as potential selective agents has remained challenging. Here,

we provide a statistical model to measure varying selective optima from longitudinal data. We apply our approach to analyze environmental sensitivity of phenotypic selection on egg-laying date and clutch size throughout the life cycle of a white-throated dipper population. We show the presence of a joint optimal phenotype that varies over the 35-year period, being dependent on altitude and temperature. We also find that optimal laying date is density-dependent, with high population density favoring earlier laying dates. By providing a flexible approach, widely applicable to free-ranging populations for which long-term data on individual phenotypes, fitness, and environmental factors are available, our study improves the understanding of phenotypic selection in varying environments.

Environments, through variation in habitats, competition, or predation, are heterogeneous at both temporal and spatial scales. Such changes in biotic and abiotic conditions impose selection on wild populations [START_REF] Bell | Fluctuating selection: the perpetual renewal of adaptation in variable environments[END_REF], leading phenotypic traits to be constantly shaped and reshaped by the environment and numerous agents of natural selection (Endler 1986). Although often ignored, population density is one of these potential selective agents [START_REF] Saether | Demographic routes to variability and regulation in bird populations[END_REF]. For a long time, evolution has been considered too slow and too weak to leave a signature in ecological dynamics [START_REF] Slobodkin | Growth and regulation of animal populations[END_REF]). However, it is now widely accepted that rates of evolution can be rapid and strong [START_REF] Pelletier | The evolutionary demography of ecological change: linking trait variation and population growth[END_REF]Ozgul et al. 2009[START_REF] Ozgul | Coupled dynamics of body mass and population growth in response to environmental change[END_REF][START_REF] Pemberton | Evolution of quantitative traits in the wild: mind the ecology[END_REF][START_REF] Bell | Fluctuating selection: the perpetual renewal of adaptation in variable environments[END_REF][START_REF] Schoener | The newest synthesis: understanding the interplay of evolutionary and ecological dynamics[END_REF]. Interestingly, observations of phenotypic selection in free-ranging populations also indicate that evolution may vary through space [START_REF] Hedrick | Genetic polymorphism in heterogeneous environments[END_REF][START_REF] Endler | Geographic variation, speciation, and clines[END_REF][START_REF] Hereford | A quantitative survey of local adaptation and fitness tradeoffs[END_REF]Siepielski et al. 2013) and time [START_REF] Siepielski | It's about time: the temporal dynamics of phenotypic selection in the wild[END_REF][START_REF] Bell | Fluctuating selection: the perpetual renewal of adaptation in variable environments[END_REF][START_REF] Morrissey | Directional selection in temporally replicated studies is remarkably consistent[END_REF]. A landmark case study of varying selection is the beak size variation in Darwin's finches in response to droughts. While drought events have favored large beaks well adapted to large seeds, high precipitation have selected for smaller beak sizes particularly useful for consuming small and soft seeds [START_REF] Grant | Unpredictable evolution in a 30-year study of Darwin's finches[END_REF]) such that the optimal phenotypic value maximizing fitness is moving as a result of fluctuating environmental conditions [START_REF] Charlesworth | Directional selection and the evolution of sex and recombination[END_REF][START_REF] Visser | Estimating the variation, autocorrelation' and environmental sensitivity of phenotypic selection[END_REF]Chevin et al. 2015). Chevin et al. (2015) provided a method for estimating varying phenotypic selection from measurements of a fitness-related trait across time. It assesses the support for stabilizing selection and for an optimal phenotype possibly influenced by environmental covariates (hereafter called environmental sensitivity of selection sensu Chevin et al. (2010)) and random effects autocorrelated across years. Using great tits (Parus major) as a case study, this work showed autocorrelated variations in the optimal egg-laying date that maximizes the number of offspring surviving to the fledgling stage. In addition, the optimal date was well predicted by spring temperature. This study left some questions unanswered, in particular how to measure varying phenotypic selection on multiple correlated traits and also throughout multiple episodes of selection.

Indeed, natural selection does not operate on a single trait but acts jointly and correlatively on multiple characters, and the environment causes this complex selection to change in a more or less predictable way. This very fundamental and widely acknowledged vision of how adaptation to changing environment proceeds still fails to be detected in a comprehensive way. The classical approach of Lande and Arnold (1983) allows estimating variation in linear and quadratic selection gradients acting on multiple traits including correlational selection (Phillips and Arnold 1989b;[START_REF] Sinervo | Correlational selection and the evolution of genomic architecture[END_REF] over time (Engen et al. 2012). This multivariate selection analysis thus provides important information on the direction, shape, and strength of selection acting on multiple phenotypic traits over time by linking relative fitness to trait values. However, relating such phenotypic selection gradients to environmental factors may lead to incomplete representations of how the environment interacts with the trait-fitness relationships [START_REF] Hunter | Quantification and decomposition of environment-selection relationships[END_REF]. Indeed, this variation not only reflects variation in the fitness function (i.e., the relationship between individual expected fitness and individual phenotype, see [START_REF] Walsh | Individual fitness and the measurement of univariate selection[END_REF]), but is complicated by the response to selection generated by the varying fitness function and phenotypic changes caused by other evolutionary forces. As a simple illustration, if the mean phenotype tracks a varying optimum almost perfectly (e.g., through plasticity), little variation in selection gradients will be detected using the traditional method of Lande and Arnold (1983). Conversely, an evolutionary force such as genetic drift will generate varying gradients, even if the fitness function is constant (Chevin and Haller 2014). In other words, considering phenotypic selection coefficients such as selection gradients or differentials alone does not necessarily allow one to characterize how the fitness function has changed, because changes in the distribution of phenotype can change selection coefficients, independently of changes in the fitness function.

Natural selection does not operate on a single episode but the strength and the direction of selection on a trait may change from one life-cycle stage to another (Engen et al. 2011;Chevin et al. 2017). This was recognized more than thirty years ago by [START_REF] Arnold | On the measurement of natural and sexual selection: theory[END_REF] who highlighted the need to measure selection through separate episodes of selection across the life cycle. However, if selection is estimated separately for each lifecycle segment as in Engen et al. (2012), this leads to a loss of parsimony when different episodes are similarly influenced by the same environmental covariates or random processes.

Here, we analyze fluctuating fitness functions through time and space in a Norwegian white-throated dipper population (Cinclus cinclus) by extending the approach from Chevin et al. (2015).

In particular, we explore the dynamics of selective optima through time and space on two key fitness-related traits in such a small passerine [START_REF] Newton | Population limitation in birds[END_REF], namely egg-laying date and clutch size. For many breeding females, information on egg-laying date and clutch size is available annually, thus providing the required data to develop a multivariate approach. We evaluate the effects of biotic and abiotic factors such as weather conditions, altitude, and densities as well as random unobserved drivers on spatio-temporal variation in selective optima of the two traits. To make efficient use of all the data, we estimate varying phenotypic selection through several episodes of selection in a single joint model, from egg to fledgling stage (episode 1), from fledgling to recruit stage (episode 2), and via mother survival (episode 3), thanks to the availability of individual-based data from birth to death.

Importantly, instead of modeling variation in selection gradients, we model fluctuations in the fitness function directly. Parameter estimates from our model can thus more easily be related to the theory on varying selection, both in time [START_REF] Bull | Evolution of phenotypic variance[END_REF]Lande and Shannon 1996;[START_REF] Lande | Expected relative fitness and the adaptive topography of fluctuating selection[END_REF][START_REF] Visser | Estimating the variation, autocorrelation' and environmental sensitivity of phenotypic selection[END_REF] and space [START_REF] Kirkpatrick | Evolution of a species' range[END_REF]. While much of this theory predominantly deals with simple life histories with discrete, nonoverlapping generations, our approach provides a statistical model for estimating varying selection acting in an age-structured population by applying a stochastic trait-dependent Leslie matrix [START_REF] Caswell | Matrix population models[END_REF], including covariates such as population density. Further theoretical work will be needed to understand the evolutionary response to selection described by our modeling approach, through some integrated measure of overall lifetime fitness. Still, in our dipper case study, even in the absence of such theory, we find that the observed pattern in mean phenotypic trait values, qualitatively

VARYING SELECTION IN THE WILD

behaves as expected in response to our estimated spatial and temporal variation in selective optima.

Methods

STUDY SPECIES AND DATA COLLECTION

The studied population is located in the river system of Lyngdalselva in southern Norway (58 • 08'-58 • 40'N, 6 • 56'-7 • 20'E). The white-throated dipper is a short-lived passerine bird distributed in mountainous regions across the Palearctic. It depends on open water for foraging and running water for nesting. The amount of ice during the winter thus influences the availability of feeding and breeding habitats explaining why survival and fecundity rates are affected by mean winter temperature (temp) (December-February) of the whole region called Sørlandet (http://www.yr.no/sted/Norge/Vest-Agder/Audnedal/Konsmo6 051/klima.vinter.html) [START_REF] Saether | Population dynamical consequences of climate change for a small temperate songbird[END_REF][START_REF] Loison | Disentangling the sources of variation in the survival of the European dipper[END_REF]Nilsson et al. 2011a;[START_REF] Gamelon | Interactions between demography and environmental effects are important determinants of population dynamics[END_REF]. Demographic rates are also strongly density-regulated and annual estimates of the number of breeding females in the population (variable dens) are available from a previous study [START_REF] Gamelon | Interactions between demography and environmental effects are important determinants of population dynamics[END_REF]. For each year t = 1979, 1980, . . . , 2013, all breeding sites were visited during the nest building period to identify breeding pairs and record occupied nests. Assuming that covariates are missing at random, we based the analysis on a subset of i = 1, 2, . . . , 546 breeding events with nonmissing covariates out of a total of 1880 observed breeding events. The altitude (alt i ) (ranging from sea level to an altitude of about 600 m above sea level) and the identity of each breeding territory (k = 1, 2, . . . , 167) were recorded. During visits in the breeding season, ringed mothers were identified and unringed mothers given a ring to allow future identifications ( j = 1, 2, . . . , 375). For each breeding pair, the egg-laying date (date of first egg laying) and the clutch size (y 0i ) was determined (Nilsson et al. 2011b). We express egg-laying date as the number of days elapsed since 1 st of March (z i ). On average, 20 days later, fledglings that had survived were ringed and their number recorded (y 1i ). Finally, the next season, that is, on average 343 days later, a number y 2i of ringed fledglings were recorded as recruited to the breeding population if they were caught breeding. Mothers caught again breeding the following year were recorded as having survived (y 3i = 1); otherwise, they were considered as dead (y 3i = 0) (Fig. 1). Therefore, survival on episodes 2 and 3 ( p 2i and p 3i ) correspond to apparent survival, that is the probability for a female to survive and stay in the population until the next breeding season. Survival on the first episode ( p 1i ) corresponds to true survival, because there is no possible dispersal during the first episode. Note that the annual recapture rate is high during the studied period, ranging between 88 and 92% (estimates obtained in a previous work [START_REF] Gamelon | Interactions between demography and environmental effects are important determinants of population dynamics[END_REF])), meaning that virtually all the females alive were caught breeding. Moreover, age of the mothers (a i ) was determined. The oldest breeding female recorded in our population was 10 years of age. Also, multiple-clutching was sometimes observed in that population. Thus, we also recorded the total number of eggs produced by a female a given year.

STATISTICAL MODEL

Using the above data, our aim is to estimate varying selective optima, extending the method of Chevin et al. (2015) to multiple traits (individual laying dates and clutch sizes) and to multiple episodes of selection (egg-to-fledgling survival s = 1, fledglingto-recruit survival s = 2, and adult female survival s = 3) through associated survival probabilities p si , s = 1, 2, 3 (Fig. 1). A more technical discussion of other differences from the Chevin et al. (2015) method is given in Appendix A. We emphasize that our model of stabilizing selection does not necessarily imply a fitness optimum within the range of observed phenotypic values in any particular generation. Instead, as in theoretical models such as Lande and Shannon (1996); [START_REF] Hansen | Stabilizing selection and the comparative analysis of adaptation[END_REF]; [START_REF] Bürger | Evolution of genetic variability and the advantage of sex and recombination in changing environments[END_REF], overall individual fitness is a strictly decreasing function on both sides of some optimal trait value possibly located outside this range. The particular model we implement, at least when survival is low such that the model becomes approximately Gaussian, also implies that selection acts to reduce the phenotypic variance and that selection gradients (sensu Lande and Arnold 1983) change linearly with the distance from the assumed optimum.

Before describing in detail how we implement stabilizing selection, we describe more generally how covariates and random effects (some of which induce varying selection) need to enter into the model given that the different episodes differ in duration. As recommended by [START_REF] Ergon | The utility of mortality hazard rates in population analyses[END_REF], we model all three survival probabilities p si only indirectly via effects of covariates and random effects on the hazard function. More specifically, we allow a nonconstant hazard (instantaneous mortality rate) with respect to age but assume that covariates and random effects act with a constant multiplicative effect on this rate (via a log link) within each episode of selection. Thus, the hazard for an offspring or an adult female at age a is λ si (a) = λ 0 (a) exp(η si ).

(1)

Here, η si is a (nonlinear) predictor containing fixed and random effects on the hazard associated with the i'th breeding event during selection episode s = 1, 2, 3, and λ 0 (a) is a baseline hazard affecting all individuals, possibly varying with age a within each interval. The survival probabilities associated with each of the three episodes of selection are given by where a s-1 and a s is the age at the beginning and end of selection episode s. Note how λs is the mean of the possibly nonconstant baseline hazard λ 0 (a) during selection episode s.

p si = exp - as as-1 λ 0 (a) exp(η si ) da = exp (-exp(η si ) λs (a s -a s-1 )), (2)
To model selection on the laying date z i and on clutch size y 0i , we in turn assume that the (time-averaged) hazard of individual i during each selection episode has the form exp(η si ) λs = exp η

(β) si + 1 2ω 2 s z i -η (θ) si 2 , ( 3
)
with the parameter ω s determining the rate of proportional increase in the hazard with increasing deviations of the phenotypic laying date z i from the optimal laying date η (θ) si . ω s thus corresponds to the width of the fitness peak (smaller ω s causes stronger stabilizing selection). Here

η (θ) si = θ 0s +
+ σ s u ti + τ s v ki + κ s w ji + x ai ,s (5)
is another linear subpredictor determining the hazard at the optimal laying date containing effects of covariates as well as a number of random effects (details are given in the next 4 paragraphs and Appendix A). The above regression coefficients must not be confused with the selection gradient β as defined by Lande and Arnold (1983). All parameters possibly differ between episodes s = 1, 2, 3, but can also be constrained to the same value for different subsets of episodes. Importantly, this facilitates the formulation of more parsimonious model alternatives in cases where the evidence for any difference between episodes is small. Note that the log of mean baseline hazard has been absorbed in the possibly episode-dependent intercept β 0s = ln λs in (5). Before going through the details of the linear predictors in 4 and 5, note first that the expected number of recruits produced by a given female (the fecundities f a in the first row of a prebreeding census Leslie matrix), assuming that a single clutch is laid, are given by products of clutch size y 0 and the survival probabilities for the two first episodes of selection f a (y 0i , z i ) = y 0i p 1i p 2i . These fecundities are important fitness components and correspond to the number of young produced during the breeding season in year t that have survived until the next year t + 1 (see Fig. 1). Even without a quadratic effect of clutch size in (5), provided that increasing clutch sizes translates to a reduction in overall survival during episodes s = 1 or 2 (β clutchsize,s sufficiently positive), it follows that these fecundities are maximized for some intermediate clutch size, as expected through the tradeoff between offspring number and offspring survival [START_REF] Smith | The optimal balance between size and number of offspring[END_REF]. Note that the model may predict an optimum located outside the range of observed phenotypic values. Thus, our model specifies a joint optimum for the two phenotypic traits laying date and clutch size for which f a (y 0 , z) is maximized. To obtain a more standard measure of the strength of stabilizing selection acting jointly on both traits, we evaluated the matrix of second derivatives of ln f a (y 0 , z) at the joint optimum to obtain the parameters of the Gaussian approximation of the fitness function (appearing in many theoretical models, e.g., Phillips and Arnold (1989a); [START_REF] Chevin | Genetic constraints on adaptation to a changing environment[END_REF]Tufto (2017)). We report the widths ω y0 and ω z of this approximation with respect to each trait (in units of number of eggs and number of days, respectively), analogous to standard deviations of Gaussian distributions.

Temporal covariates and temporal random effects appearing in 4 and 5 translate to variation from year to year in the optimal laying date and clutch size, respectively. To model possibly correlated fluctuations in the joint optimum as in [START_REF] Chevin | Genetic constraints on adaptation to a changing environment[END_REF] as well as autocorrelation across time (as in Lande and Shannon 1996;Lande 2009;[START_REF] Visser | Estimating the variation, autocorrelation' and environmental sensitivity of phenotypic selection[END_REF]Chevin et al. 2017), the random effects representing yearly variation in overall survival u t (5) and variation in the optimal laying date ζ t (4) are assumed to follow a first-order vector autoregressive VAR(1) process

u t ζ t = u t-1 ζ t-1 + w t , ( 6
)
where is a 2 × 2 matrix of autoregressive coefficients and w t is bivariate normal N (0, ) white noise. This only specifies the autocorrelation matrix function (see Wei 2006, ch. 16.1) of the process (u t , ζ t ). But as long as the variance of u t is small, optimal clutch size will be approximately linearly dependent on u t and so the autocorrelation matrix function of the joint optimal clutch size and laying date will be almost identical to that of (u t , ζ t ).

Correlation between u t and ζ t can arise either through , or both having nonzero off-diagonal entries. If and are both diagonal, this simplifies to two independent AR(1) processes and if all entries of are zero, u t and ζ t are simple independent and identically normally distributed (iid) white noise processes. We parameterize this part of the model in terms of , the white noise correlation ρ = 12 / √ 11 22 , the stationary variance σ 2 ζ of ζ t , and with the stationary variance of u t in (5) set equal to one but with separate parameters σ s representing the potentially different effects of u t on the three selection episodes. Note also that additional correlation in variation of the joint optimum as well as autocorrelation across time can be induced through the temporal fixed effect covariates appearing in each linear predictor. Correlated optima discussed above are distinct from correlational selection. Two kinds of correlational selection can be accommodated in our non-Gaussian model for f a (y 0i , z i ), either by adding z i y 0i as a covariate in (5) making the optimal clutch size dependent on laying date (first kind), or by adding clutch size y 0i as an additional covariate in (4) making the optimal laying date dependent on clutch size (second kind).

To model the effect of the age of breeding female, we consider models where the log hazard during the different episodes is either independent of age (a single term β 0,s in ( 5)), linearly dependent on age (an additional term β age,s age i in ( 5)) or nonlinearly dependent on age. Instead of modeling non-linear age-dependence parametrically (using for example quadratic, Gompertz, piecewise linear or two-parameter Weibull models [START_REF] Gaillard | Assessing senescence patterns in populations of large mammals[END_REF]; [START_REF] Marzolin | Frailty in state-space models: application to actuarial senescence in the Dipper[END_REF])), we use a nonparametric approach: we model the age effects through a set of correlated random effects, more precisely as a second-order random walk. This is a commonly used method for smoothing data and modeling response functions [START_REF] Green | Nonparametric regression and generalized linear models: a roughness penalty approach[END_REF][START_REF] Rue | Gaussian Markov random fields: theory and applications[END_REF]. It provides a simple and flexible way to model the hazard variations as a function of age, and thus to explore senescence. Under the second-order random walk model, the joint distribution of the age effects x s,1 , x s,2 , . . . , x s,10 (last term in (5), Fig. 2, upper plot, blue curve) is specified by assuming that the second-order differences 2 x s,a are independently normally distributed with zero mean and variance ν s . The parameter ν s (estimated along with the random age effects themselves) thus controls the magnitude of these second-order differences (analogous to the second-order derivative) and hence the smoothness of the resulting function. Further details on the implementation of the model and the two last random effects terms appearing in 5 are given in Appendix A.

Apart from the nonlinear predictor, our model based on the assumption of multiplicative effects on the hazard corresponds to a generalized linear mixed model with a log-log link (or complementary log-log) function with the log of the length of the selection episode included as an offset variable. Surprisingly, this link and its associated assumption of multiplicative effects on the hazard is rarely used in ecology and evolution. However, this is a standard assumption in many models in survival analysis [START_REF] Ergon | The utility of mortality hazard rates in population analyses[END_REF], for example in Cox proportional hazards model, and seems more reasonable from a biological point of view. We must point out, however, that the magnitude of the variation in optimal clutch size induced by temporal covariates and random effects, depends, to some extent, on our choice of link function. For the alternative logit link, the induced variations implied by the model would have been somewhat smaller and would go to zero in the limit of low survival where the logit tends to a log link. The logit and other link functions such as the probit, however, would not lead to regression coefficients having a common interpretation across episodes of selection of different length [START_REF] Ergon | The utility of mortality hazard rates in population analyses[END_REF]). Hence, those link functions would not facilitate the fitting of sometimes more parsimonious null models for which these regression coefficients are constrained to a common value across several selection episodes.

Results and Discussion

Our statistical model applied to the dipper allows us to select the most parsimonious model (Table 1) among all the tested ones (Table 2). In this section, we explore in detail all the effects retained in the best model and discuss their implications.

AGE EFFECTS ON SURVIVAL

The best model indicates no effect of mother age on offspring survival from egg to fledgling stage (episode 1) and from fledgling to recruit stage (episode 2) (see last row in Table 1, Fig. 2). Considering for instance a linear effect of female age on offspring survival does not improve the model fit (models B to E, Table 2). Indeed, model B that includes a different linear age effect for episodes 1 and 2 provides estimates not significantly different from zero ( βage,1 = -0.009 ± 0.039 for episode 1 and βage,2 = -0.029 ± 0.024 for episode 2). Because the probability of rearing a chick is independent of mother age, our findings indicate no senescence in maternal care. Note also that we do not detect senescence in clutch size (results not shown here). Mean survival on the first episode is close to 83% and drops around 7% on the second episode (Fig. 2). This high mean survival on episode 1 compared to episode 2 (about 10 times higher) simply results from the much longer duration of the second episode. At the adult stage (episode 3), survival is age-dependent, increasing from age 1-4 and decreasing from age 4 onwards. This result indicates actuarial senescence (Fig. 2), caused by a progressive loss of cellular and physiological functions late in life [START_REF] Williams | Pleiotropy, natural selection, and the evolution of senescence[END_REF][START_REF] Hamilton | The moulding of senescence by natural selection[END_REF]. Senescence is pervasive in the wild [START_REF] Nussey | Senescence in natural populations of animals: widespread evidence and its implications for bio-gerontology[END_REF]), and we provide here additional evidence in a short-lived species.

ENVIRONMENTAL EFFECTS ON SURVIVAL

As expected for this species that strongly depends on open water for foraging, warmer winters favor adult survival (see third row in Table 1). Therefore, adult females are more likely to survive, stay, and breed the next breeding season in the population when the winters are mild. Apparent survival from fledgling to recruit stage is not affected by winter temperatures (model Y, estimate = 0.043 (SE: 0.042)). In accordance with previous studies [START_REF] Gamelon | Interactions between demography and environmental effects are important determinants of population dynamics[END_REF]Nilsson et al. 2011a;[START_REF] Saether | Population dynamical consequences of climate change for a small temperate songbird[END_REF][START_REF] Loison | Disentangling the sources of variation in the survival of the European dipper[END_REF], high population density increases competition among individuals and thus mortality (and possibly dispersal rate), especially at the adult stage (see fourth row in Table 1). Including an effect of density on survival on other episodes does not improve the model fit (models AA to AI, Table 2). However, it is noteworthy that a model including an additional effect of density on episode 1 is close to the best model (model AC), but the low effect size indicates no density-dependent mortality on this first episode ( βdens,1 = -0.003 ± 0.003). We thus do not find any evidence for density-dependent mortality at the offspring stage, contrary to some other passerine bird species such as great tit (Reed et al. 2013a;[START_REF] Saether | Demographic routes to variability and regulation in bird populations[END_REF] for which low densities are generally associated with high offspring survival. Notice also that high altitudes negatively affect offspring survival from egg to fledgling stage and also from fledgling to recruit stage (see fifth row in Table 1).

SURVIVAL AND VARYING SELECTIVE OPTIMA FOR CLUTCH SIZE

Mortality on episode 1 decreases in large clutches (see second row in Table 1, first axis on Fig. 3A). However, because of high overall survival during the first short selection episode, this does not translate to strong selection for large clutch sizes. Low mortality on episode 1 in large clutches also means that females with large clutches have offspring with the highest survival in the nest. This suggests individual heterogeneity, with some females performing better than others in terms of reproductive success. It is noteworthy that females that lay a large number of eggs, not only during a given breeding event but generally during the whole breeding season, do not pay direct survival costs. Indeed, adding a term for the effect of the total number of eggs laid during the whole breeding season in (5) does not provide any improvement as indicated by model F (Table 2, βeggs,3 = -0.062 ± 0.061). For episode 2, the best model indicates that large clutch sizes increase mortality (Table 1, Fig. 3B). This might be due to malnutrition and reduced parental care [START_REF] Noordwijk | Genetic and environmental variation in clutch size of the great tit (Parus major)[END_REF]. As a result of lower survival during the second episode, this translates to overall survival over the two first episodes combined decreasing with increasing clutch size. Therefore, laying too many eggs is associated with increasing offspring mortality. But obviously, laying too few eggs is not a successful breeding tactic for a female either [START_REF] Lack | The natural regulation of animal numbers[END_REF][START_REF] Boyce | Optimizing great tit clutch size in a fluctuating environment[END_REF][START_REF] Both | Adaptive density dependence of avian clutch size[END_REF]. Here, we provide evidence for an optimal clutch size maximizing offspring survival and more generally overall fecundity rate estimated to be 6.69 eggs. The width of the fitness function with respect to clutch sizes, ω y0 , is estimated to be between 4.89 and 5.89 eggs in different years, being approximately proportional to the optimum in any given year. The fitness function is also somewhat asymmetric with respect to clutch size (Fig. 3C). This high value for the width of the fitness function indicates weak stabilizing selection for optimal clutch size, as also illustrated in Fig. 3C.

We find that this optimal clutch size varies through time between 6.19 and 7.53 eggs (Fig. 4A). This agrees quite well with observed clutch sizes in the population, mainly ranging between 4 and 6 eggs (Table 3). The mean observed clutch size of 4.52 eggs is somewhat smaller than the estimated overall optimum. However, this observed mean is influenced by the left tail of the frequency distribution, which might reflects total failure of some Clutch size 0 1 2 3 4 5 6 7

β dens,1 = 0, β dens,2 = 0, β dens,3 = 0 AJ 1.88 -1 β alt,1 = β alt,2 = 0 AK 2.9 1 β alt,1 = 0, β alt,3 = 0 AL 1.97 1 β alt,1 = 0, β alt,2 = 0 AM 1.97 1 β alt,3 = 0 AN 3.41 0 β alt,2 = 0 AO 0.94 0 β alt,1 = 0 AP 2.11 0 β alt,1 = 0, β alt,2 = β alt,3 = 0 AQ 3.74 0 β alt,1 = β alt,2 = 0, β alt,3 = 0 AR 1.27 0 β alt,1 = β alt,2 = β alt,3 AS 3.95 2 β alt,1 = 0, β alt,2 = 0, β alt,3 = 0 AT 1.03 1 β clutchsize/date,1 = 0 (correlational selection of first kind) AU 2 1 β clutchsize/date,2 = 0 (correlational selection of first kind) AV 2.13 1 θ 0,1 = 0, θ 0,2 = 0, θ 0,3 = 0 AW 1.99 1 θ t,1 = θ t,2 AX 2.14 2 θ clutchsize,1 = 0, θ clutchsize,2 = 0 (correlational selection of second kind) (Continued)
= 0 BD 1.32 1 θ alt,1 = 0, θ alt,2 = 0 BE 1.99 1 σ 2 = 0, σ 3 = 0 BF 2 1 σ 1 = 0, σ 2 = σ 3 = 0 BG 4 2 σ 1 = 0, σ 2 = 0, σ 3 = 0 BH 2.26 2 τ 1 = 0, τ 2 = 0, τ 3 = 0 BI 1.8 1 τ 1 = 0, τ 2 = 0, τ 3 = 0 BJ -0.18 0 τ 1 = τ 3 = 0, τ 2 = 0 BK 0.28 1 τ 1 = τ 3 = 0, τ 2 = 0 BL 5.2 3 κ 1 = 0, κ 2 = 0, κ 3 = 0 BM 3.2 2 κ 1 = 0, κ 2 = 0, κ 3 = 0 BN 1.4 1 κ 1 = κ 3 = 0, κ 2 = 0 BO 1.28 1 κ 2 = κ 3 = 0, κ 1 = 0 BP 3.96 2 σ ζ >
Frequency 57 1 52 69 186 636 192 4 females. Interestingly, the observed mean clutch size maximizes total fecundity rate (Fig. 3C). According to our model, the temporal variation in optimal clutch size is induced by random variations in survival during episode 2 modeled by the random effect term σ 2 u t in (5) (σ s = 0.144 ± 0.050 for episode s = 2, Table 1) such that larger clutches turn out to be favored in years with high survival. Modeling u t as an autoregressive process does not improve the model (Table 2, model G), that is, we find no evidence for autocorrelation in these variations.

Finally, the models including correlational selection between egg-laying date and clutch size on maximum survival on the first two episodes do not perform better than the best model (models I and J, Table 2). There is thus no evidence for correlational selection between clutch size and egg-laying date on offspring survival.

It should be noted that we have treated clutch size as a trait on which selection operates, rather than as a fitness component. This is reasonable because in many bird populations, recruitment (function of clutch size) has a limited influence on the population growth rate (see [START_REF] Saether | Demographic routes to variability and regulation in bird populations[END_REF] for a comparative analysis). Instead, a large proportion of the temporal variance in population change of temperate passerines is due to variation in survival [START_REF] Gould | Estimation of temporal variability of survival in animal populations[END_REF]Saether et al. 2004), dipper being no exception [START_REF] Loison | Disentangling the sources of variation in the survival of the European dipper[END_REF]. It is therefore relevant not to consider clutch size as a fitness component but rather as a trait under selection.

SURVIVAL AND VARYING SELECTIVE OPTIMA FOR EGG-LAYING DATE

We find that egg-laying date is under stabilizing selection. Indeed, removing stabilizing selection in (3) for episodes 1 and 2 does not provide any improvement in AIC (model K, Table 2). Interestingly, ω s , that determines the strength of stabilizing selection, is similar for episodes 1 and 2 (see 11th row in Table 1). Indeed, estimating two different values ( ω1 = 33.37 ± 9.1 and ω2 = 49.11 ± 12.8) does not improve the model fit (model L, Table 2). Given the much longer duration of episode 2, however, most of the selection happens during this episode. The approximate width of the overall fitness function ω z varies between 27.76 and 32.01 days reflecting strong stabilizing selection on egg-laying date, as illustrated in Fig. 3C. Adding stabilizing selection for episode 3 does not improve the model fit (model M,

Clutch size 0 2 4 6 8 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••••••••••••••••••••••••••••••••• •• • • • • • • • • • • •
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • B
Observed phenotype Overall optimum Optimum for episode 1 Optimum for episode 2 1). In plot B, the red and green curves are the estimated optimal laying dates in terms of survival during episodes 1 and 2, respectively, and the black curve the optimal laying date for both episodes combined. Note that the overall optima (black curve) nearly coincide with the optima of the episode 2. All optima are estimated at the average altitude and for the most frequent clutch size (five eggs). Gray lines represent mean phenotypic values and the size of the gray dots the frequencies of different phenotypes in the total population. The minor tick marks in plot B are located as in Fig. 3.

Table 2), meaning laying eggs early or late in the breeding season has no effect on mother survival. In contrast, in a tropical parrot, the green-rumped parrotlet (Forpus passerinus), there is strong selection on egg-laying date through adult survival, such that females that breed early in the season exhibit lower survival to the next breeding season [START_REF] Tarwater | Opposing selection and environmental variation modify optimal timing of breeding[END_REF]. This discrepancy between our finding and Tarwater and Beissinger (2013) might be explained by the contrasting life-history strategies of these two species. Indeed, parrotlets have a much slower pace of life than the dipper. In slow-living species, adult survival has the highest contribution to population growth rate and is thus expected to be particularly canalized (Gaillard and Yoccoz 2003). This explains why viability selection is important in slow-living species like parrotlets compared to short-lived species such as dipper. As a consequence, our results indicate that the strong stabilizing selection on egg-laying date operates only through offspring survival.

The timing of egg laying and more generally the timing of reproduction is critical in many species [START_REF] Price | Directional selection and the evolution of breeding date in birds[END_REF] because hatching/birth should match with good environmental conditions in terms of weather and/or food resources. Otherwise, offspring survival may be jeopardized. Our findings provide evidence for an optimal laying date in the dipper with an overall mean egg-laying date maximizing survival from egg to fledgling stage estimated to be θ0,1 = 40.9 ± 9.6 days after the 1 st of March (i.e., mid-April, Fig. 3A, see sixth row in Table 1) and with an overall mean egg-laying date maximizing survival from fledgling to recruit stage estimated to be θ0,2 = 42.1 ± 5.5 days after the 1 st of March (Fig. 3B, Table 1). This translates to a mean egg-laying date maximizing total fecundity rate estimated to be April 12 (Fig. 3C). This mean optimal date is slightly earlier than the mean laying date actually observed over the study period (April 22 for n = 741 clutches for which information on egg-laying date was available, second axis, Fig. 3C).

Observed egg-laying dates have advanced at a rate of 0.15 ± 0.07 days/year during the 35-year period (Fig. 4B, gray dots). Interestingly, we find a significant trend toward earlier optimal dates at a rate of θt,s = 1.37 ± 0.59 days/year (Table 1).

Removing the trend worsens model fit considerably (model N, Table 2). This estimate seems somewhat large and would imply an advance in the optimum of 47 days over the course of the study. The lower confidence limit of 7.24 days for this advance (based on approximate normality of θt,s ) seems more reasonable and comparable to the observed phenotypic change of 8.8 days on average (ranging from 4 to 17 days), recorded in other bird populations over a 25-year period [START_REF] Crick | UK birds are laying eggs earlier[END_REF].

Importantly, in addition to the linear trend, the estimated optimal laying date (at the average altitude) varies over years between February 8 and May 24 (Fig. 4B, black line) as a result of the effects of winter temperature and population density. Removing the effect of environmental covariates on the optimal laying date worsens model fit (models S and T, Table 2) thus providing strong evidence for environmental sensitivity of selection. In particular, when the winter following the breeding season is warm, the optimal laying date maximizing survival from fledgling to recruit stage (i.e., episode 2) is delayed by θtemp,2 = 8.9 ± 3.5days/ • C (see eighth row in Table 1). One can hypothesize that when the subsequent winter is mild, offspring produced at late laying dates during the previous breeding season are disproportionately more likely to survive, thus generating selection for later laying dates. While including the same effect of winter temperature on the optimal laying date of episode 1 led to a slight improvement

• • • • • • •• • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • •• • • •• ••• • • • • • • • • • • • • • • • • • •• •• • • • • • • • • • • • • • • • • • •• •• • •
• •• • • • • • • ••••••
• • • • • •• • • • •• • • • •• •• • • ••• • • • • • •• • •• • • •• • • •
• •• • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • ••
• • • • • • • • • • • • • •• • • • • • • • • •• • • • • • • • • • • • • • • •• •• •• • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • •• ••• • • • •• •• • • • • •• • • • • • • • • •• • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • •• • •• • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • •
Mean phenotype value Estimated optimum in AIC, such a model would clearly be biologically unrealistic as survival on episode 1 has to be causally independent of the following winter conditions (Fig. 1). This model alternative was thus excluded from consideration (see Burnham and Anderson 2002, ch. 6.8.7). Under our best model, different optima for episodes 1 and 2 are therefore estimated (Fig. 4B, red and green curves). The overall optimum (black curve) is approximately an average of the optima for each of the two episodes, weighted by the respective strength of stabilizing selection for each episode (this approximation would be exact if the fitness functions for each episode were exactly Gaussian). Indeed, as can be seen from Fig. 4 B, the overall optimum always falls between the green and red curve but much closer to the green curve given the much stronger stabilizing selection during episode two. Interestingly, previous studies on other passerine bird species have shown that selection on egg-laying date often depends on the timing of the peak in caterpillars, the main food resources, itself closely related to spring temperatures (e.g., in great tits, (Chevin et al. 2015;[START_REF] Visser | Shifts in caterpillar biomass phenology due to climate change and its impact on the breeding biology of an insectivorous bird[END_REF]Reed et al. 2013b)). Here, we do not find statistical evidence for the effect of spring conditions (date of ice break-up) on the optimal egg-laying date. Based on the upper confidence limit, the regression coefficient is considerably smaller than its expected value (see Appendix A), suggesting that other environmental factors are more important for the dipper's optimal laying date. This discrepancy between our findings and previous works may be explained by the difference in the biology of these species. While great tit strongly depends on insect avail-ability during spring, food resources are available on a larger time window (over spring and summer) for the dipper.

We also find that optimal egg-laying date is densitydependent, with high densities favoring earlier optimal date (see seventh row in Table 1). Despite a large amount of theoretical and experimental works on the role of density as a selective agent (Charlesworth 1994;[START_REF] Engen | A quantitative genetic model of r-and K-selection in a fluctuating population[END_REF][START_REF] De Lisle | Correlated evolution of allometry and sexual dimorphism across higher taxa[END_REF][START_REF] Shaw | Response to density in a wild population of the perennial herb salvia lyrata: variation among families[END_REF]), empirical evidence in the wild remains scarce (but see [START_REF] Saether | Demographic routes to variability and regulation in bird populations[END_REF] for great tits and [START_REF] Hunter | Quantification and decomposition of environment-selection relationships[END_REF] for Soay sheep (Ovis aries)). In green-rumped parrotlets for instance, late breeding dates are selected for when the density is low, through enhanced adult survival [START_REF] Tarwater | Opposing selection and environmental variation modify optimal timing of breeding[END_REF]. Similarly, a recent study has shown that in North American red squirrels (Tamiasciurus hudsonicus), high population densities might increase the strength of selection for earlier birth dates [START_REF] Fisher | Multilevel and sex-specific selection on competitive traits in North American red squirrels[END_REF], favoring successful recruitment for juveniles. Here, our findings provide additional support for a key role of density as a selective agent in the wild.

Having included the temporal covariates winter temperature and population density, we do not find any latent variations in optimal laying date as indicated by the lack of improvement in AIC when including ζ t (corresponding to random variation in the optimal laying date) as a iid random effect (model P, Table 2). We also considered including ζ t distributed as an AR(1) process (φ 22 = 0, model Q), cross-correlated with u t-1 (corresponding to variation in overall survival) (φ 21 = 0, model H) or crosscorrelated with u t+1 (φ 12 = 0, model R) but neither of these model alternatives led to any improvement in AIC.

Through their estimated joint effect on the optimal laying date, winter temperature and population density induce autocorrelation on the deviations of the optimal laying dates from the estimated linear trend. For episode 1, only influenced by population density and no additional latent random process, the autocorrelation function of the optimum is identical to that of population density, with a significant autocorrelation of 0.58 at lag 1 (Fig. A3). For episode 2, despite also being influenced by population density in addition to winter temperature, the resulting optimum given by θdens,2 dens t + θtemp,2 temp t exhibited no significant autocorrelation, mainly because of the larger effect of winter temperature (exhibiting no autocorrelation, Fig. A2) accounting for 75% of the total variance in the optimum. Since the optimal laying date for episodes 1 and 2 combined almost coincides with the optimum for episode 2 (Fig. 4B, black and green curves), the same applies to the corresponding autocorrelation function.

Finally, in addition to the dependency on winter temperature and population density, we find that optimal and observed egglaying dates depend on altitude (model U, Table 2), occurring later at high altitudes (see tenth row in Table 1, black line on Fig. 5). This result gives clear support for spatial variation in phenotypic selection and provides empirical evidence for adaptation along an environmental gradient. To what extent this can be explained by the altitudinal gradient in temperatures would require finegrained data on local temperature, unfortunately not available. Interestingly, the slope in the mean phenotype based on the egglaying dates in the complete data (gray line in Fig. 5, n = 741) is significantly shallower than the estimated slope in the optimum (black line in Fig. 5, one sided Z -test, P-value = 0.041). Thus, laying tends to occur too late at low altitudes and too early at high altitude relative to the estimated optimal laying date. This suggests maladaptation at the extreme ends of the gradient, possibly caused by gene flow as in [START_REF] Kirkpatrick | Evolution of a species' range[END_REF] model of evolution of species' range or by imperfect plasticity [START_REF] Tufto | The evolution of plasticity and nonplastic spatial and temporal adaptations in the presence of imperfect environmental cues[END_REF]Chevin and Lande 2011;[START_REF] Gienapp | Why climate change will invariably alter selection pressures on phenology[END_REF].

MODEL WITHOUT TEMPORAL COVARIATES

Our statistical approach can be used when environmental covariates involved in the varying optima are missing. As an illustration, we removed population density and winter temperatures as temporal covariates as well as the trend in optimal egg-laying dates and tested different models for the random effects instead (see Table A1). The best model (see Fig. A1 and Table A2) included both latent variations in overall mortality (through the u t term inducing variation in optimal clutch size as before) and in optimal laying date (through a significant ζ t term, model BT vs best model, σζ = 19 ± 12 days), jointly following a vector autoregressive process (see ( 6)). In line with the absence of autocorrelation in the optima induced by the temporal covariates for model A in Table 2 (Fig A3 ), we found no evidence for autocorrelation in optima via u t and ζ t for the model without temporal covariates (models BU and BV in Table A1). Surprisingly, this model includes a negative autoregressive coefficient φ12 = -0.45 ± 0.95, making ζ t-1 (optimal laying date in year t -1) negatively correlated with u t (the hazard in year t). A possible explanation is that an unknown temporal covariate influences the optimal laying date with a delayed effect on survival.

The moving optimal laying dates and clutch sizes estimated with the best models with and without temporal covariates are generally in accordance (Fig. 4 vs Fig. A1). The parameter estimates provided by the model without temporal covariates are close to the ones of our best model for the effects of clutch size and altitude on mortality (comparison between Table 1 andTable A2). However, including temporal covariates improves the precision of most parameter estimates (Table 1 vs Table A2). In particular, the precision of the estimates of optimal egg-laying date strongly differs between the two models, being equal to θ0,1 = 22 ± 29 days after the 1 st of March (i.e., March 30 ± 29 days) when excluding temporal covariates and to 40.9 ± 9.6 days (i.e., April 10 ± 9.6 days) for the first episode and to θ0,2 = 42.1 ± 5.5 days (i.e., April 12 ± 5.5 days) for the second episode when including them. Similarly, ω s for episodes 1 and 2 is estimated to be 45 ± 12 and 81 ± 45 for the best models with and without temporal covariates (translating to widths ω z of the Gaussian approximation of the fitness function ranging from 28 to 32 days and 39 to 62 days, respectively). This difference also explains the shift in the estimated optimal laying dates (Fig. 4B vs Fig. A1B) since the location of the optimum is, to some extent, estimated through extrapolation of the quadratic model of the effect of laying date on the hazard beyond the range of observed phenotypic values (see (3)).

Conclusion

Understanding how life-history traits vary in time and space and determining the selective forces behind this variation is one of the central issues in ecology and evolution. Instead of working within the classical framework of Lande and Arnold (1983), we have extended the statistical glm-like approach used in several works [START_REF] Janzen | Logistic regression for empirical studies of multivariate selection[END_REF]Shaw and Geyer 2010;Chevin et al. 2015), by allowing movements of the underlying fitness function. These movements can be induced through observed environmental covariates or latent processes following different types of plausible autoregressive models. This extension leads to a nonlinear latent variable statistical model, efficiently handled thanks to modern statistical software. A particular advantage of our approach is that all the data for several episodes are utilized in a single joint model. Also, correlational selection on multiple traits and correlated optima are accommodated as possible model alternatives.

The statistical approach we have used leads to models of stabilizing selection that are different from the standard Gaussian fitness function frequently used in theoretical models. Our approach is similar to how a quadratic effect of a trait in the logistic regression approach of [START_REF] Janzen | Logistic regression for empirical studies of multivariate selection[END_REF] translates into a non-Gaussian, plateauing fitness function if survival is high. Our view is that such non-Gaussian functions are more realistic because they reflect the constraint that survival probabilities necessarily have an upper bound of one. The difference from our approach is in the exact link function used (the logit link vs the loglog link corresponding to our proportional hazard assumption). The same type of argument can be made for how non-Gaussian stabilizing selection on clutch size emerges in our model. This being said, if survival is low, for example during selection on laying date in episode s = 2, and not strongly dependent on clutch size, the resulting fitness function is well approximated by a Gaussian function with widths ω z and ω y0 easily derived from the basic parameters of the model (Table 1). It should also be noted that when estimated optima fall outside the range of observed phenotypic values, the existence of an optimum is not an inference drawn from the data alone. Instead, we assume that an optimum exists and this is a reasonable assumption. For instance, based on a priori biological knowledge, we know that the reproductive success of a female that starts breeding too early in the season under harsh winter conditions will be low. The locations of optima are in turn estimated based on the most parsimonious model of the fitness curvature supported by the data.

In our dipper case study, using model selection criteria to choose between a large number of alternative models, we find evidence for varying selective optima on two key life-history traits. Spatio-temporal variation in optimal laying dates is induced by variation in altitude, winter temperature, and population density. Optimal clutch sizes, resulting from an estimated trade-off between offspring number and survival, exhibit similar random variations over time. Selection on these traits mainly operates through low survival from fledgling to recruit stage (episode 2), paralleling the key role of this life stage as a driver of fluctuations in avian population dynamics [START_REF] Saether | Demographic routes to variability and regulation in bird populations[END_REF]. For adult survival, we clearly show that senescence occurs but we find no cost of large clutch sizes or selection on laying dates via adult survival.

Our statistical approach can also be used when environmental covariates involved in the varying optima are missing. Indeed, the moving optimal laying dates and clutch sizes estimated with the models with and without temporal covariates are generally in accordance (Fig. 4 vs Fig. A1). However, including relevant environmental covariates is important to improve the overall precision of parameter estimates and of course, to identify the agents of selection.

Our main objective in the present study has been to estimate the pattern of varying selective optima acting at various life stages. We are not aware of any simple theory for how this translates to variation in selection acting over the whole lifespan of an individual. It is noteworthy that the variation in mean phenotypes through space and time in many respects are qualitatively similar to variation in the estimated optima. However, we cannot conclude that these patterns in mean phenotypes are consistent with the pattern of varying selection estimated by our method. Further work is needed to quantify the expected evolutionary response that is expected from our model. Given the complications of age-structured and density-dependent models (Engen and Saether 2017), this is likely feasible only via simulation-based approaches. given by equations ( 3), (4), and (5) (a function of y 0 , z, v, and w), and a sa s-1 the duration of episode s. In the terminology of generalized linear-mixed models, this gives us the so called marginal as opposed to conditional model [START_REF] Agresti | Categorical data analysis[END_REF], Section 12.2.2) with respect to v and w. Fixed effect covariates were set to their mean values and other temporal random effects were set to zero (in Fig. 3) or to their estimated values (in Fig. 4) (giving us a conditional model with respect to ζ t and u t ). Computing the arithmetic mean fitness in this way, averaging over these random effects rather than considering the fitness function conditional on their values, appears to be the most reasonable approach, at least under the assumption of hard selection (nonlocal density regulation occurring after locally varying selection, see Tufto (2015) Displayed are the meaning of the parameters, their notations, and their estimates (± standard errors) for the three episodes of selection. location. Modeling variation in the elevation of fitness functions across years in this way makes better use of the data as some of the information contained in mean survival in a given year to some extent, depending on the magnitude of the random effect on the fitness maximum, is informative about the location of the optimum.

Another important improvement over the Chevin et al.'s (2015) method is the inclusion of individual level covariates (such as altitude) directly influencing the location of the phenotypic optimum at the individual level. This makes the overall predictor (3) nonlinear also in the fixed effect parameters (the parameter θ alt,s in 4). In the INLA-based approach of Chevin et al. (2015), a term for this could naively be included as an interaction term between trait value and altitude in the linear predictor. But again, this would translate to an unrealistic quadratic relationship between the maximum of the fitness function and the fitness optima experienced by individuals at different altitudes.

EFFECT OF SPRING CONDITIONS ON OPTIMAL EGG-LAYING DATE

As the species mainly feed underwater, we used the timing of ice break-up during spring as a measure of spring conditions. Ice cover break-up was defined as when there no longer was a connective ice layer across the southern end of the lake Lygne located in the middle of the dipper system, where the outlet is 1 st of March was set as day 1 and dates were sequentially numbered. This information was available from 1979 to 2009. From 2010 to 2013, the timing of ice break-up was set to its mean observed between 1979 and 2009, that is 46 days after 1 st of March.

From the best model retained (Table 2), we tested an additional effect of the timing of ice break-up (denoted spring) on optimal egg-laying date, similar for episodes 1 and 2 (i.e., θ spring,1 = θ spring,2 = 0, AIC=0.89). We found that θspring,1 = θspring,2 = -0.229 ± 0.276, in the opposite direction and different from the theoretical value of θ spring,s = 1 expected if the optimal laying date occurs at a fixed number of days after ice break up. While there may still be an effect that is not detected because of low statistical power, the effect would have to be quite small, based on the upper approximate confidence 95% confidence limit of 0.31. We also considered additional model alternatives by excluding the effect of mean winter temperatures (i.e., θ spring,1 = θ spring,2 = 0, θ temp,2 = 0). Once again, this model did not show any improvement ( AIC = 9.16). In addition, we tested an effect of the timing of ice break-up on optimal egg-laying date, different for episodes 1 and 2 (i.e., θ spring,1 = 0, θ spring,2 = 0, AIC = 2.23). We also evaluated the effect of the timing of ice break-up on optimal egg-laying date on episode 1 only (i.e., θ spring,1 = 0, AIC = 0.23) and on episode 2 only (i.e., θ spring,2 = 0, AIC = 1.99). Finally, we tested the effect of the timing of ice break-up on optimal egg-laying date on episode 1 only while excluding the effect of mean winter temperatures (i.e., θ spring,1 = 0, θ temp,2 = 0, AIC = 8.19). None of these models improved the fit.

Figure 1 .

 1 Figure 1. Dipper life cycle. Episode 1 corresponds to the episode from egg to fledgling stage (in red), episode 2 from fledgling to recruit stage (in green), and episode 3 corresponds to adult female stage (in blue).

Figure 2 .

 2 Figure 2. The log of the hazard (per year) (upper plot) and the resulting probability of surviving each episode (lower plot) as a function of mother age. Mean hazard and survival are computed at optimal egg-laying date and at the most frequent clutch size (five eggs), at mean altitude, population density, and winter temperature for each episode.

CFigure 3 .

 3 Figure 3. Estimated (A) survival rate on episode 1 (from egg to fledgling stage, in red), (B) survival rate on episode 2 (from fledgling to recruit stage, in green), and (C) total fecundity rate (i.e., total survival rate from egg to recruit stage × clutch size, in black) as a function of egg-laying dates and clutch sizes. Mean vital rates are computed at average altitude, population density, and winter temperature. The red cross indicates mean clutch size and mean egg-laying date observed during the study period. The minor tick marks are at the 10 th , 20 th , and 30 th of each month.

  0 and φ 11 = 0 and φ 22 = 0 (u t and ζ t ∼ AR(1)) BQ 4 2 σ ζ > 0 and φ 12 = 0 and φ 21 = 0 BR 7.96 4 σ ζ > 0 and φ 11 , φ 12 , φ 21 and φ 22 = 0 Subscripts indicate the selection episode under consideration, that is 1 for the first episode from egg to fledgling stage, 2 for the second episode from fledgling to recruit stage and 3 for the third episode corresponding to the adult stage.

Figure 4 .

 4 Figure 4. Estimated annual optimal clutch size (plot A) and laying dates (plot B) for the best model (Table1). In plot B, the red and green curves are the estimated optimal laying dates in terms of survival during episodes 1 and 2, respectively, and the black curve

Figure 5 .

 5 Figure 5. Estimated optimal (black line) and observed mean (gray line) egg-laying dates as a function of altitude (mean centered), together with observed egg-laying dates (gray dots). Dot size illustrates the number of individuals exhibiting the same laying date at a given altitude. Estimated optimal egg-laying dates are computed at average population density and winter temperature. The minor tick marks are located as in Fig. 3.

Figure A2 .

 A2 Figure A2. Sample autocorrelation-and cross-autocorrelation functions (see e.g., Shumway and Stoffer 2011, Defs. 1.14 and 1.41) for the time series of observed population densities and temperatures. The off-diagonal plots displays estimates of cor r(dens t , temp t-k ).

Figure A3 .

 A3 Figure A3. Estimated autocorrelation functions of optimal laying dates θdens,s dens t + θtemp,s temp t induced by population density and winter temperature for episodes s = 1 and s = 2.

  θ dens,s dens i + θ temp,s temp i + θ alt,s alt i + ζ ti (4)

	is a linear subpredictor determining the optimal laying date during
	episode s containing possible effects of environmental covariates
	such as population density (dens i ), winter temperature (temp i ),
	and altitude (alt i ) as well as a random effect term ζ t as in Chevin
	et al. (2015). Similarly,
	η (β) si = β 0s + β clutchsize,s y 0i
	+ β temp,s temp i + β dens,s dens i + β alt,s alt i + β age,s a i

Table 1 . Parameters estimated with the best model retained (Table 2).

 1

	Selection episode

Table 2 . Model selection-Displayed are the tested models derived from the best model retained. AIC and p are the relative differences to the best model in AIC values and number of parameters p, respectively.

 2

		AIC	p	Description
	A	0	0	Best model, see Table 1
	B	2.44	2	β age,1 = 0, β age,2 = 0 (different linear age effect for episodes 1 and 2)
	C	1 . 9 7	1	β age,1 = 0 (linear age effect for episode 1)
	D	0.66	1	β age,1 = β age,2 = 0 (common linear age effect for episodes 1 and 2)
	E	0 . 4 9	1	β age,2 = 0 (linear age effect for episode 2)
	F	0.94	1	β #eggs, 3
	G	1 . 9 6	1	φ 11 = 0 (u t ∼ AR(1))
	H	2	1	σ ζ > 0 and φ 12 = 0 (u t dependent on ζ t-1 )
	I	2 . 9 7	2	β clutchsize/date,1 = 0, β clutchsize/date,2 = 0 (correlational selection of first kind)
	J	1.56	1	β clutchsize/date,1 = β clutchsize/date,2 = 0 (correlational selection of first kind)
	K	16.73	-6	Remove stab. sel. on z (quadratic term in eq. 3) for episodes 1 and 2
	L	0.5	1	ω 1 = ω 2 (strength of stabilizing selection different for episodes 1 and 2)
	M	2 . 0 2	1	θ 0,3 = 0 (stabilizing selection for episode 3)
	N	5.73	-1	θ t = 0 (remove trend in optimal laying date)
	O	1 . 2 1	1	θ clutchsize,1 = θ clutchsize,2 = 0 (correlational selection of second kind)
	P	2	1	σ ζ > 0 (random effect on optimal laying date)
	Q	4	2	σ ζ > 0 and φ 22 = 0 (ζ t ∼ AR(1))
	R	3.93	2	σ ζ > 0 and φ 21 = 0 (ζ t dependent on u t-1 )
	S	9 . 0 6	-1	θ temp,2 = 0 (remove temperature effect on optimal laying date)
	T	7.31	-1	θ dens,1 = θ dens,2 = 0 (remove density effect on optimal laying date)
	U	9 . 8 8	-1	θ alt,1 = θ alt,2 = 0 (remove altitude effect on laying date)
	V	16.9	0	All age effects modelled as second order random walks
	W	2 . 4 8	-2	β clutchsize,1 = β clutchsize,2 = 0
	X	45.88	-1	β temp,3 = 0
	Y	0 . 7	1	β temp,2 = 0
	Z	24.39	0	β temp,2 = β temp,3
	AA	10.39	-1	β dens,3 = 0
	AB	1.46	1	β dens,2 = 0
	AC	0.98	1	β dens,1 = 0
	AD	11.62	0	β dens,1 = 0, β dens,3 = 0
	AE	14.34	0	β dens,1 = β dens,3 = 0
	AF	7.71	0	β dens,2 = β dens,3
	AG	12.36	0	β dens,2 = 0, β dens,3 = 0
	AH	10.4	0	β dens,1 = β dens,2 = β dens,3
	AI	2.66	2

Table 2 .

 2 Continued.

		AIC	p	Description
	AY	2.03	0	θ dens,1 = 0
	AZ	7.54	0	θ dens,2 = 0
	BA	1.44	1	θ dens,1 = 0, θ dens,2 = 0
	BB	2.96	0	θ alt,1 = 0
	BC	11.38	0	θ alt,2

Table 3 . Observed frequencies of different clutch sizes in the population.

 3

Table A1 . Model selection-Displayed are all the tested models derived from the best model without temporal covariates (i.e., best model without population density and winter temperatures, see Table A2 for a list of parameters included in the best model), the difference AIC (respectively p) in AIC (respectively in num- ber of parameters p) between each model and the best one and their description.

 A1

		AIC	p	Description
	BS	0	0	Best model
	BT	6.74	-1	φ 12 = 0
	BU	1.79	1	φ 11 = 0
	BV	1.76	1	φ 22 = 0
	BW	1.59	1	φ 21 = 0
	BX	7.1	0	φ 11 = 0, φ 12 = 0
	BY	5.35	0	φ 22 = 0, φ 12 = 0

Table A2 . Parameters estimated with an alternative model without any temporal covariates (population density and winter tempera- tures) but with autocorrelated random effects.

 A2

			Selection episode
	Parameters		Egg to fledgling	Fledgling to recruit	Adult female
	Meaning	Symbol	s = 1	s = 2	s = 3	U n i t
	Intercept, mortality	β0,s	1.94 ± 0.40	0.68 ± 0.29	0
	Effect of clutch size on mortality	βclutchsize,s	-0.170 ± 0.080	0.082 ± 0.055	0
	Effect of altitude on mortality	βalt,s	0.00084 ± 0.00053	0.00084 ± 0.00053	0	m -1
	Intercept, optimal egg-laying date	θ0,s	22 ± 29	22 ± 29	0	d
	Effect of altitude on optimal date	θalt,s	0.20 ± 0.12	0.20 ± 0.12	0	dm -1
	Quadratic effect of egg-laying date	ωs	81 ± 45	81 ± 45	0
	SD of random year effect	σs	0	0 .032 ± 0.067	0.032 ± 0.067
	SD of random territory effect	τs	0.181 ± 0.055	0.181 ± 0.055	0.181 ± 0.055
	Beta-binomial dispersion parameter	γs	1.915 ± 0.139	1.072 ± 0.084
	SD of second-order random age effect	νs	0	0	0.34 ± 0.17
	VAR(1) autoregressive parameters	ˆ	0	-0.45 ± 0.95
			0	0
	SD of latent fluctuations in optimum date	σζ	19 ± 12			d

ACKNOWLEDGMENTS

We thank all the persons involved in the field work. We also warmly thank Luis-Miguel Chevin, Ophélie Ronce, Michael Morrissey, and three anonymous referees for helpful comments on previous drafts and Hans J. Skaug for useful discussions. This work was partly supported by the Directorate for Nature Management (Norwegian Environment Agency) and the Research Council of Norway through its Centres of Excellence funding scheme (project number 223257 to C.B.D. and 179569 to C.E.E.S.).

DATA ARCHIVING

All the necessary data and scripts to carry out the analyses presented here are available in the dryad digital repository https://doi.org/ 10.5061/dryad.k871qr2

AUTHOR CONTRIBUTIONS

K.J. and O.W.R. contributed to field work and data collection; A.L.K.N. provided data on egg-laying dates; M.G., J.T., N.C.S., and B.E.S. conceived the study; J.T. developed the statistical model; M.G. and J.T. conducted the analyses; M.G. and J.T. wrote the manuscript and all authors contributed to revisions on later manuscript versions.

Appendix A

STATISTICAL MODEL DETAILS

To accommodate overdispersion in the number of surviving fledglings y 1i and in the number of surviving recruits y 2i , we use a beta-binomial distribution for each episode s = 1, 2, parameterized in terms of the survival probability p si and a dispersion parameter γ s (the factor by which the variance is inflated relative to the simpler binomial model). Conditional on the initial clutch size y 0i , the joint distribution of number of fledglings y 1i and number of recruits y 2i is then

where betabin(y; n, p, γ) is the probability mass function of the beta-binomial distribution. Similarly, the survival of the adult female y 3i associated with breeding event i is modeled as Bernoulli distributed with parameter p 3i . The other random effects v k and w j appearing in (5), based on territory and female identity, are included to model possible positive correlation between number of fledglings and recruits produced by the same territory and female in different years. We parameterized the model such that they are iid standard normal and like u t , their effect on the different episodes s = 1, 2, 3 are potentially different, depending on whether their standard deviations τ s and κ s differ between episodes.

In practice, the expected number of recruits produced as function of clutch size y 0 and laying date z, f a (y 0 , z) (Fig. 3C), was computed by evaluating

using numerical integration (R-package cubature), thus integrating out the above nontemporal random territory and female identity effects v and w. Here f is the standard normal probability density function, η s (y 0 , z, v, w) is the nonlinear predictor Appendix S3). The joint optimum of f a (y 0 , z) (Figs. 4 andA1) and the second derivatives of its log was also computed numerically (R function optim).

To implement the model, software for fitting generalized linear-mixed models cannot be used, because η si in (3) is nonlinear in the parameters and the random effects. Instead, we used Template Model Builder [START_REF] Kristensen | TMB: automatic differentiation and Laplace approximation[END_REF]) (R-package TMB) that provides a general framework for fitting complex, nonlinear, random effects (latent variable) models. Briefly, the user defines the joint likelihood for the data and the random effects as a C++ template function. Based on this template, TMB generates a function computing the Laplace approximation of the marginal likelihood. This is, in turn, maximized numerically to obtain the maximum likelihood estimates of the model parameters. Biologically meaningful models derived from variations of (4) and ( 5) were fitted (Table 2). We selected the best model using the Akaike information criteria AIC [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF] among all tested models and recovered the estimates of all parameters together with their associated standard errors. Explanatory variables alt i , temp i and dens i were mean centered. A2).

An important advantage of the modeling framework provided by TMB is that models that are nonlinear in the parameters and random effects can be fitted with little effort from the user perspective. In contrast, Chevin et al. (2015) used the INLA R-package (Rue et al. 2009) to estimate a log-linear model for the Poisson mean with a constant quadratic term representing stabilizing selection and a random autocorrelated effect on the slope. This translates to autocorrelated fluctuations in the moving optimal phenotype. This approach requires modeling variation in the elevation of the fitness function in different years through a fixed effect on the intercept (the term μ t in their eq. ( 3)). If instead, variation in elevation was modeled through a random effect on the intercept, this would translate to an undesirable and unrealistic quadratic relationship between the expected fitness at the optimum and the location of the optimum. In contrast, the TMB framework we have used here allows more biologically realistic models to be fitted with the inclusion of random effects acting directly on the fitness maximum (or equivalently, on the minimum of the hazard function) and on its Associate Editor: M. Morrissey Handling Editor: Mohamed A. F. Noor