Pascal ; Marrot
email: pascal.marrot@usherbrooke.ca

Anne Charmantier

Jacques Blondel

Dany Garant

Current Spring Warming as a Driver of Selection on Reproductive Timing in a Wild Passerine

Keywords: Natural selection, Cyanistes caeruleus, Spring warming, Laying date, Blue tit, Model averaging Natural selection, Cyanistes caeruleus, Spring warming, Laying date, Blue tit, Model averaging

1. Evolutionary adaptation as a response to climate change is expected for fitness-related traits affected by climate and exhibiting genetic variance. Although the relationship between warmer spring temperature and earlier timing of reproduction is well documented, quantifications and predictions of the impact of global warming on natural selection acting on phenology in wild populations remain rare. If global warming affects fitness in a similar way across individuals within a population, or if fitness consequences are independent of phenotypic variation in key-adaptive traits, then no evolutionary response is expected for these traits.

2. Here we quantified the selection pressures acting on laying date during a 24-year monitoring of blue tits in southern Mediterranean France, a hot spot of climate warming. We explored the temporal fluctuation in annual selection gradients and we determined its temperature-related drivers.

3. We first investigated the month-specific warming since 1970 in our study site and tested its influence on selection pressures using a model averaging approach. Then, we quantified the selection strength associated with temperature anomalies experienced by the blue tit population.

4. We found that natural selection acting on laying date significantly fluctuated both in magnitude and in sign across years. After identifying a significant warming in spring and summer, we showed that warmer daily maximum temperatures in April were significantly associated with stronger selection pressures for reproductive timing. Our results indicated an increase in the strength of selection by 46% for every +1°C anomaly.

Introduction

Global climate change results in an increase in temperature mean and variability and generates new environmental conditions for wild populations [START_REF] Coumou | A decade of weather extremes[END_REF].

Rising temperatures are known to alter the timing of vegetation development [START_REF] Fu | The Impact of Winter and Spring Temperatures on Temperate Tree Budburst Dates: Results from an Experimental Climate Manipulation[END_REF], leading to shifts in phenology and range shifts of many taxa [START_REF] Chen | Rapid Range Shifts of Species Associated with High Levels of Climate Warming[END_REF]) such as plants [START_REF] Svenning | Disequilibrium vegetation dynamics under future climate change[END_REF], birds [START_REF] Crick | The impact of climate change on birds[END_REF]), insects [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF] and mammals [START_REF] Moritz | Impact of a Century of Climate Change on Small-Mammal Communities in Yosemite National Park, USA[END_REF]. These responses to climate warming can be due to non-genetic changes (i.e. phenotypic plasticity) or to evolutionary adaptations [START_REF] Gienapp | Climate change and evolution: disentangling environmental and genetic responses[END_REF].

Whereas evidence of phenotypic plasticity in this context is often reported in the literature [START_REF] Merilä | Climate change, adaptation, and phenotypic plasticity: the problem and the evidence[END_REF], studies showing an evolutionary response to climate change are still rare [START_REF] Gienapp | Climate change and evolution: disentangling environmental and genetic responses[END_REF][START_REF] Merilä | Evolution in response to climate change: In pursuit of the missing evidence[END_REF]). This "missing evidence" [START_REF] Merilä | Evolution in response to climate change: In pursuit of the missing evidence[END_REF]) is worrying because phenotypic plasticity alone is unlikely to sustain population responses to climate change over a long term warming period [START_REF] Hoffmann | Climate change and evolutionary adaptation[END_REF].

Evolutionary adaptation under climate change is expected if climate-induced selective pressures are found for key adaptive traits with significant genetic variance [START_REF] Hoffmann | Climate change and evolutionary adaptation[END_REF]. While most life history traits affected by the recent climate warming were shown to possess some level of genetic variance in the wild [START_REF] Postma | Four decades of estimating heritabilities in wild vertebrate populations: Improved methods, more data, better estimates. Quantitative genetics in the wild[END_REF], strong selective pressures resulting from climate change have always been assumed, but rarely quantified empirically [START_REF] Maccoll | The ecological causes of evolution[END_REF]. Global warming is often expected to generate strong selection pressures because numerous phenotypic mismatches between adaptive optima and population phenotypic mean have been documented following a warming trend [START_REF] Gienapp | Climate change and evolution: disentangling environmental and genetic responses[END_REF]. In particular, theoretical modelling has shown that climate change should increase selection on consumer phenology, even if the environmental changes affecting consumer and resource phenologies are of similar magnitude [START_REF] Gienapp | Why climate change will invariably alter selection pressures on phenology[END_REF]). However, one of the basic assumptions for natural selection to occur is that fitness varies among individuals and that this variance in fitness is related to their phenotypes. Hence, for climate warming to induce a selective pressure on wild populations, it should cause differences in individual fitness depending on specific phenotypes. In other words, if a phenotypic mismatch affects all individuals similarly, climate warming will not result in a new selective pressure, and no evolutionary response is expected. So far, although many shifts in phenology have been shown in response to warming across all living taxa [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF], few studies have attempted to assess changes in selection under climate change [START_REF] Domínguez | Rainfall and flowering synchrony in a tropical shrub: Variable selection on the flowering time ofErythroxylum havanense[END_REF]McAdam & Boutin 2003;[START_REF] Reale | Genetic and plastic responses of a northern mammal to climate change[END_REF][START_REF] Maad | Variable selection in Platanthera bifolia (Orchidaceae): phenotypic selection differed between sex functions in a drought year[END_REF][START_REF] Nussey | Selection on Heritable Phenotypic Plasticity in a Wild Bird Population[END_REF][START_REF] Husby | Contrasting Patterns of Phenotypic Plasticity in Reproductive Traits in Two Great Tit (Parus Major) Populations[END_REF][START_REF] Karell | Climate change drives microevolution in a wild bird[END_REF][START_REF] Anderson | Phenotypic plasticity and adaptive evolution contribute to advancing flowering phenology in response to climate change[END_REF][START_REF] Campbell | Natural selection on floral morphology can be influenced by climate[END_REF][START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF][START_REF] Kingsolver | Climate variability slows evolutionary responses of Colias butterflies to recent climate change[END_REF][START_REF] Visser | Effects of Spring Temperatures on the Strength of Selection on Timing of Reproduction in a Long-Distance Migratory Bird[END_REF]. These studies showed no clear pattern of new selective regime resulting from global warming, with both increases [START_REF] Domínguez | Rainfall and flowering synchrony in a tropical shrub: Variable selection on the flowering time ofErythroxylum havanense[END_REF][START_REF] Maad | Variable selection in Platanthera bifolia (Orchidaceae): phenotypic selection differed between sex functions in a drought year[END_REF][START_REF] Nussey | Selection on Heritable Phenotypic Plasticity in a Wild Bird Population[END_REF][START_REF] Anderson | Phenotypic plasticity and adaptive evolution contribute to advancing flowering phenology in response to climate change[END_REF]) and decreases (McAdam & Boutin 2003;[START_REF] Karell | Climate change drives microevolution in a wild bird[END_REF][START_REF] Campbell | Natural selection on floral morphology can be influenced by climate[END_REF] in the strength of natural selection. Moreover, most studies compared selection pressures between only two years (e.g. [START_REF] Domínguez | Rainfall and flowering synchrony in a tropical shrub: Variable selection on the flowering time ofErythroxylum havanense[END_REF][START_REF] Maad | Variable selection in Platanthera bifolia (Orchidaceae): phenotypic selection differed between sex functions in a drought year[END_REF]) and/or failed to identify a significant climate warming over the study period (e.g. [START_REF] Maad | Variable selection in Platanthera bifolia (Orchidaceae): phenotypic selection differed between sex functions in a drought year[END_REF][START_REF] Anderson | Phenotypic plasticity and adaptive evolution contribute to advancing flowering phenology in response to climate change[END_REF][START_REF] Visser | Effects of Spring Temperatures on the Strength of Selection on Timing of Reproduction in a Long-Distance Migratory Bird[END_REF]. Despite these limitations, previous studies showed an increase in the strength of selection along with temperatures in different birds species such as the Pied Flycatcher, Ficedula hypoleuca [START_REF] Visser | Effects of Spring Temperatures on the Strength of Selection on Timing of Reproduction in a Long-Distance Migratory Bird[END_REF], the Great Tit, Parus major [START_REF] Visser | Warmer springs lead to mistimed reproduction in great tits (Parus major)[END_REF][START_REF] Charmantier | Adaptive Phenotypic Plasticity in Response to Climate Change in a Wild Bird Population[END_REF][START_REF] Husby | Speeding Up Microevolution: The Effects of Increasing Temperature on Selection and Genetic Variance in a Wild Bird Population[END_REF][START_REF] Gienapp | Predicting demographically sustainable rates of adaptation: can great tit breeding time keep pace with climate change[END_REF][START_REF] Reed | Phenological mismatch strongly affects individual fitness but not population demography in a woodland passerine[END_REF][START_REF] Vedder | Quantitative Assessment of the Importance of Phenotypic Plasticity in Adaptation to Climate Change in Wild Bird Populations[END_REF][START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF] and the Blue Tit, Cyanistes caeruleus [START_REF] Phillimore | Passerines may be sufficiently plastic to track temperature-mediated shifts in optimum lay date[END_REF]). Overall, these previous studies explored a potential influence of temperature, partially neglecting the investigation of climate warming (i.e. an exploration of warming patterns over the study period). In order to predict the evolutionary response of populations to climate warming, it is crucial to go beyond an assessment of the impact of temperature on natural selection, and to quantify and predict variation in the strength of selection beyond the current warming period.

Here we investigated temperature-related changes in selection pressures acting on the timing of reproduction in a wild population of blue tits (Cyanistes caeruleus). The Blue tit is an insectivorous passerine living in temperate forests of Europe and in parts of the Middle East, breeding from March to June in southern France [START_REF] Blondel | A thirty-year study of phenotypic and genetic variation of blue tits in mediterranean habitat mosaics[END_REF]. The amount of food availability for offspring (mostly leaf-eating caterpillars) is dependent on the timing of budding of oaks (Quercus sp.), which is highly responsive to ambient temperature [START_REF] Van Noordwijk | Selection for the Timing of Great Tit Breeding in Relation to Caterpillar Growth and Temperature[END_REF]. Thus, like most passerines in temperate forests, blue tits must adjust their laying date to track the amount of food available [START_REF] Visser | Warmer springs lead to mistimed reproduction in great tits (Parus major)[END_REF][START_REF] Thomas | Energetic and Fitness Costs of Mismatching Resource Supply and Demand in Seasonally Breeding Birds[END_REF]). Fitness differences among birds breeding at different dates are strong with earlier breeders achieving higher reproductive success [START_REF] Blondel | A thirty-year study of phenotypic and genetic variation of blue tits in mediterranean habitat mosaics[END_REF][START_REF] Porlier | Variation in phenotypic plasticity and selection patterns in blue tit breeding time: between-and withinpopulation comparisons[END_REF]. Previous studies on a closely related species, the Great tit (Parus major) have shown that the selection differential on breeding dates is correlated with the mean temperature during the period when birds are incubating [START_REF] Van Noordwijk | Selection for the Timing of Great Tit Breeding in Relation to Caterpillar Growth and Temperature[END_REF][START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF]. Warmer springs generated an increased mismatch between great tit breeding phenology and seasonal food peak, leading to intensifying directional selection on laying dates [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][START_REF] Reed | Phenological mismatch strongly affects individual fitness but not population demography in a woodland passerine[END_REF]. Warmer spring temperatures are thus expected to drive patterns of selection for earlier laying date in the blue tit. However, contrary to the previous work exploring the temperature-induced selection on laying date in great tits in North-Western Europe, we explored the impact of warming on blue tits in the Mediterranean area, which is a warming "hot spot" (i.e. a responsive region to climate change, [START_REF] Giorgi | Climate change hot-spots[END_REF]. Indeed, a pronounced decrease in precipitation and increase in warming, especially in summer, is forecasted in this region for the next decades [START_REF] Giorgi | Climate change projections for the Mediterranean region[END_REF]. The relative importance of climate warming on ecological processes may also depend on temperature seasonality. For example, budburst seems to be more sensitive to spring warming than to winter warming [START_REF] Fu | The Impact of Winter and Spring Temperatures on Temperate Tree Budburst Dates: Results from an Experimental Climate Manipulation[END_REF].

Moreover, the influence of climate warming might depend on the type of temperature because patterns of warming are different among mean, minimum and maximum temperatures [START_REF] Rebetez | Monthly air temperature trends in Switzerland 1901-2000 and 1975-2004[END_REF]. In order to investigate in detail these climatic drivers, we explored the increase in monthly minimum, mean and maximum temperatures across the calendar year, and assessed their roles as potential determinants of selection acting on blue tit laying date. In addition to these climatic drivers, we also took into account the selective impact of two non-climatic variables: annual predation rate on nestlings and population density. Predation on nestlings by small carnivores (in particular the weasel Mustela nivalis) has strong but variable impacts on blue tit reproductive success in our study population.

Population density has previously been shown to increase the strength of selection acting on laying date (e.g. in pied flycatchers, Ficedula hypoleuca, [START_REF] Ahola | Selection on laying date is connected to breeding density in the pied flycatcher[END_REF]). These variables were taken into account because they fluctuated across years, potentially shaping the climatic signal on selection acting on laying date. Finally, after assessing the relative impact of climatic and non-climatic drivers of selection on laying date, we quantified the relationship between temperature drivers and selection strength.

Materials and Methods

Study area and field work

Our analyses were based on data from a long-term study of blue tits in the forest of La

Environmental variables

Data on daily minimum (Tn), maximum (Tx) and mean (estimated as (Tn+Tx)/2) temperatures were obtained from a weather station (43°24'N, 03°57'E) located approximately 24km from the breeding site, provided by Meteo France (https://donneespubliques.meteofrance.fr/). Anomalies in temperatures were estimated as the difference between monthly temperatures and temperatures averaged for each month across all years during the 1970-2000 period. Annual predation rate was estimated as the ratio of the number of broods predated by the number of nest boxes occupied by blue tits within each year. Because annual predation rate ranged from 0 to 0.4, and because 80% of annual predation rates were below 0.1, its skewed distribution prevented model convergence. Hence, we considered this variable as a categorical variable taking two different modalities: low and high corresponding to a predation rate below or above 0.1 respectively. Only 4 years corresponded to a high predation rate (2001, 2006, 2009 and 2010, i.e. 16% of the time span of the study). Note that variation in predation rate does not result solely from the natural dynamics of predators, but also from hunting practices and various anti-predator protocols that have changed over the monitoring years in order to limit the impact of predators on the study population. Population density was estimated as the inverse of the mean distance (in km) between each nest box occupied by blue tits and its closest neighbour (including nest boxes occupied by great tits). Instead of using the proportion of occupied nest-boxes, we used this metric to estimate local density around the nest boxes since the number and the density of nest boxes within the population changed across years.

Statistical analyses

All statistical analyses were carried out using the software R (version 3.1.1) (R Core Team 2014).

Climatic analyses

Mean, minimum and maximum monthly temperatures were estimated as the daily temperature averaged during each month. Smoothed monthly temperatures were obtained from a 2-years moving average window (supplementary material, Fig. S1). Temperature changes for each month were estimated as the slope of anomalies in monthly temperatures regressed over time (years) across the period 1970-2015.

Selection analyses

Analyses were conducted on female individual data that were not subject to any experimental manipulation, with a dataset of 1024 breeding observations on 721 female blue tits breeding between 1991 and 2014. Selection acting on laying date was defined as the slope of the regression of relative fitness (i.e. number of fledglings) on laying date [START_REF] Lande | The Measurement of Selection on Correlated Characters[END_REF]. Clutch size was also included as covariate in the fixed part of the model to estimate the selection acting on laying date independently from selection acting on the highly correlated clutch size [START_REF] Perrins | Laying Dates and Clutch Size in the Great Tit[END_REF]. Thus, selection was estimated from the effect of laying date on the proportion of a clutch that fledged (i.e. brood fledging success). Nonlinear selection acting on laying date was not analysed further because quadratic and correlational selection were not significant (results not shown). Relative fitness corresponded to the individual fitness divided by the annual mean fitness of the population, and followed a Gaussian distribution (Fig. S2). All selection analyses were conducted on laying dates standardized annually (zero mean and unit variance within each year) in order to control for an environmental covariance between fitness and the trait across years. As this dataset includes repeated values within each year and for some individuals, year and individual identities were included in the random part of the selection models to control for pseudoreplication. The significance of selection gradients was estimated using F tests using lmerTest [START_REF] Kuznetsova | Package 'lmerTest'. R package version[END_REF] package in R. Temporal variation in selection gradients across years was assessed by including an interaction term between laying date and years (considered as a discrete variable) in the random part of the selection model.

This method allows to account for errors in the estimation of annual selection gradients when exploring their fluctuation [START_REF] Morrissey | Directional Selection in Temporally Replicated Studies Is Remarkably Consistent[END_REF]. Mixed models were fitted using the lmer function from the lme4 package [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF] in R and the significance of random terms was tested by Likelihood Ratio Tests (LRT).

Model averaging approach

The model averaging analysis was conducted following methods provided by [START_REF] Grueber | Multimodel inference in ecology and evolution: challenges and solutions[END_REF] and by [START_REF] Symonds | A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike's information criterion[END_REF]. Based on information theory, the aim of model averaging is to evaluate relative supports for a given hypothesis over others instead of null hypothesis testing [START_REF] Burnham | AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons[END_REF]. In a first step, we determined which variable to include in the models. Focusing on the four months covering the blue tit breeding season (March, April, May and June), we explored the influence of 14 potential drivers (daily minimum, mean and maximum monthly temperatures (for each of the four months), predation rate and population density) on annual selection by including, in the selection models, interaction terms between laying date and the drivers. In a second step, we generated a "model set" by adjusting all possible submodels from the set of predictors of interest. All submodels included laying date, clutch size, and a combination of drivers (interaction terms between environmental variables and laying date). Since 14 environmental variables were investigated, 16384 models could potentially be fitted (2 14 ). Because collinearity among predictors can be a problem in model selection [START_REF] Freckleton | Comparative Methods as a Statistical Fix: The Dangers of Ignoring an Evolutionary Model[END_REF], we constrained the number of possible submodels by avoiding models including temperatures from the same month (e.g. maximum March temperature did not appear in a model including mean or minimum temperatures in March), leading to 1024 models generated. In a third step, we employed the model averaging procedure on our set of models (based on their AIC) in order to obtain robust environmental driver estimates (associated with their standard errors) accounting for model selection uncertainty. This approach is particularly suitable when the Akaike weight of the best model is not considerably higher than the next best model (i.e. it has an Akaike weight of >0.9) [START_REF] Symonds | A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike's information criterion[END_REF], which was the case in our analysis since the Akaike weight of the best model was 0.66. Because employing a model averaging procedure on the full set of models can be spurious (due to models with very low Akaike weight), we selected the 95% confidence set of models [START_REF] Symonds | A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike's information criterion[END_REF]. This procedure led to a set of 17 best models identified (Table S1). Moreover, relative importance of each of the 14 environmental variables was estimated by summing the Akaike weights from each model in which the specific variable appeared. This relative importance can be interpreted as the probability that the variable is a component of the best model [START_REF] Symonds | A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike's information criterion[END_REF]. A variable was considered as driving selection if its relative importance was above 80%. Note that this relative importance should not be interpreted as a p value.

Results

Patterns of temperature increase

The first aim of our analysis was to assess the recent warming in the forest of La Rouvière near Montpellier. Over the last 45 years , the trend in mean annual temperatures in this habitat was +0.44°C / decade (P < 0.001) (Fig. S1). The warming trend was slightly stronger for minimum annual temperature (+0.51°C / decade, P < 0.001) than for maximum annual temperature (+0.33°C / decade, P < 0.001). While warming was negligible (between -0.20°C and +0.20°C per decade) during winter months (Dec.-Feb.), the increase in temperatures in spring (March-May) and summer (June-August) ranged between +0.30°C and +0.80°C per decade (Fig. 1). Warming during autumn (Sept.-Nov.) was also high with trends ranging between +0.20°C and +0.90°C per decade. Minimum monthly temperatures increased all year long except in winter and maximum monthly temperatures increased more strongly in spring-early summer with more than 0.5°C per decade in March, April, May and June (Fig. 1). Our results thus show a strong seasonality in the warming trend experienced by birds at our study site, along with differences between warming in minimum and maximum monthly temperatures.

Selection on laying date and its temporal variation

The fixed effects of the model revealed that the selection gradient on laying date was significantly negative over the 24 years of monitoring (β = -0.049 ± 0.009, F = 26.15, P < 0.001). Clutch size was under a strong and significant positive selection (β = 0.119 ± 0.009, F = 149.34, P < 0.001).

Our analyses of annual variation in selection gradients for laying date revealed highly significant fluctuations (χ² = 21.86, P < 0.001) in sign and magnitude, across the 24-years of monitoring (Fig. 2). For example, while early breeders were strongly favoured in 2008

(estimates obtained from a model of selection fitted in 2008: β = -0.155 ± 0.036, F = 33.84, P < 0.001), selection favouring late breeders was found in 2009 (albeit not significant, β = 0.130 ± 0.068, F = 2.54, P = 0.133) (Fig. 2). Overall, the annual standard deviation in the slope of the selection gradients was relatively strong (σ year (β LD ) = 0.063).

Climate warming and selection

Among the 14 possible drivers tested, the model averaging analysis showed that only daily maximum April temperature and predation rate influenced selection on laying date, with a probability to be a component of the best model (relative importance; e.g. from the fixed effect part of the model) of respectively 89% and 88% (Table 1). Warm maximum April temperatures and low predation rates increased the strength of selection on laying date (since selection acting on laying date is negative overall -earlier breeders are favoured). Note that these two drivers also significantly altered the strength of selection when tested using a null hypothesis statistical framework (F = 12.12, P < 0.001 and F = 16.57, P < 0.001 for the interaction terms between laying date with predation rate and maximum April temperatures respectively).

All other potential drivers had probabilities ≤ 6% of influencing breeding selection patterns.

Given the values of the overall linear selection gradient for laying date (-0.050 ± 0.009) and the interaction term between laying date and maximum April temperature (-0.024 ± 0.013), these results suggest an increase in the strength of selection by 46% (95% confidence intervals (CI): 30%-63%) for a +1°C anomaly (Fig. 3, Table 1). This means that the difference in relative fledging success between breeders laying 38 days apart (5 standard deviation in laying date, i.e. encompassing 97% of the laying dates) are 0.32 and 0.72 for 0°C and +3°C

anomalies respectively (Fig. S3). To illustrate further the selective effect of anomalies in maximum April temperatures, we conducted selection models on three sub-datasets, depending on temperature: 12 years, 3 years and 5 years in which anomalies in April temperature were lower than 1°C, ranged from 1°C to 2°C, and were higher than 2°C respectively. Results from these three models of selection confirmed that higher temperature anomalies translated into stronger selection acting on laying date (β LD = -0.0267 ± 0.0139, β LD = -0.0789 ± 0.0335 and β LD = -0.1307 ± 0.0201 for T°C < 1°C, 1°C < T°C < 2°C and T°C > 2°C respectively; see Fig. 4). While selection acting on laying date was marginally nonsignificant for years experiencing a temperature anomaly lower than 1°C (F = 3.72, P = 0.054), the selection gradient was significant for years experiencing a temperature anomaly higher than 1°C (F = 5.54, P = 0.021 and F = 42.30, P < 0.001, for 1°C < T°C < 2°C and T°C > 2°C respectively). Finally, the selection gradients acting on clutch size were not stronger with greater temperature anomalies (β CS = 0.1321 ± 0.0139, β CS = 0.1310 ± 0.0336 and β CS = 0.0703 ± 0.0196 for T°C < 1°C, 1°C < T°C < 2°C and T°C > 2°C respectively; all P-value were < 0.001).

Discussion

Our analyses confirmed that selection favors earlier breeding in Mediterranean blue tits overall, but also that natural selection showed important fluctuations both in sign and magnitude over the past 25 years. Importantly, we showed that daily maximum April temperature was the main climatic variable correlated with the strength of selection on laying date.

In temperate forests, spring temperatures are tightly linked to the annual laying date adjustment of insectivorous passerines [START_REF] Visser | Temperature has a causal effect on avian timing of reproduction[END_REF]) that are faced with the challenge of matching maximum caterpillar availability to nestling food demand (e.g. [START_REF] Van Noordwijk | Selection for the Timing of Great Tit Breeding in Relation to Caterpillar Growth and Temperature[END_REF][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]). As a result, avian ecologists have emphasized that warmer spring temperatures should contribute to strengthening selection for earlier laying date [START_REF] Crick | The impact of climate change on birds[END_REF][START_REF] Gienapp | Why climate change will invariably alter selection pressures on phenology[END_REF]).

This assumption is potentially misleading because natural selection is caused by a variance in fitness and phenotypic performances of individuals, which have rarely been shown to be driven by climate warming (but see [START_REF] Reed | Phenological mismatch strongly affects individual fitness but not population demography in a woodland passerine[END_REF][START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF]. In other words, a mismatch between trophic levels induced by warming temperatures does not always represent a selection pressure, in particular if all individuals are affected similarly by the warming. Hence, our results represent a rare demonstration of an increase in natural selection strength following a strong spring warming. This warming experienced by blue tits was linked to an increase in the variation in their expected fledging success, depending on their timing of reproduction. This result is concordant qualitatively with those found in two great tit populations monitored for more than 40 years in the Netherlands (e.g. 1), we expect a thermal sensitivity B = -9.6 days/°C in the present Mediterranean blue tit population. This discrepancy between our and previous estimates may have several origins, such as the time period considered for the temperature measure and/or the warming intensity [START_REF] Phillimore | Passerines may be sufficiently plastic to track temperature-mediated shifts in optimum lay date[END_REF]. Indeed, previous work investigated the thermal sensitivity during a period that was preliminary identified to be the best predictor of the laying date using a sliding-window analysis (Bailey & van de Pol 2016a). Moreover, the only estimation of B for blue tits in the literature [START_REF] Phillimore | Passerines may be sufficiently plastic to track temperature-mediated shifts in optimum lay date[END_REF]) revealed a disparity by a factor of three between latitude (B lat ) and longitude (B long ) thermal sensitivity (B lat = -3.84 days/°C and B long = -9.34 days/°C, see Table S2 in [START_REF] Phillimore | Passerines may be sufficiently plastic to track temperature-mediated shifts in optimum lay date[END_REF]. Note also that [START_REF] Phillimore | Passerines may be sufficiently plastic to track temperature-mediated shifts in optimum lay date[END_REF] found no disparity between B lat and B long estimates for great tits, which could suggest a species-dependent thermal sensitivity. Since the period, as well as the type of measure (minimum, maximum or average temperature), varies systematically among studies, it greatly limits the quantitative comparison in selection sensitivity.

Even within a single species, life history responses to temperature variation is known to differ among populations [START_REF] Husby | Contrasting Patterns of Phenotypic Plasticity in Reproductive Traits in Two Great Tit (Parus Major) Populations[END_REF]. Since the Mediterranean area is considered a hotspot of climate disruption [START_REF] Giorgi | Climate change hot-spots[END_REF], similar studies on the force of natural selection across the blue tit distribution range would provide crucial comparative insights. Up to now, the selective influence of climatic warming has rarely been showed in the wild, with some previous studies showing an increase in the strength of natural selection in plants with climatic warming (e.g. [START_REF] Domínguez | Rainfall and flowering synchrony in a tropical shrub: Variable selection on the flowering time ofErythroxylum havanense[END_REF][START_REF] Maad | Variable selection in Platanthera bifolia (Orchidaceae): phenotypic selection differed between sex functions in a drought year[END_REF][START_REF] Anderson | Phenotypic plasticity and adaptive evolution contribute to advancing flowering phenology in response to climate change[END_REF], while others showed a decrease, mainly in animals (McAdam & Boutin 2003;[START_REF] Karell | Climate change drives microevolution in a wild bird[END_REF][START_REF] Campbell | Natural selection on floral morphology can be influenced by climate[END_REF]. Moreover, estimating the drivers of selection remains challenging because of several interactions and covariances between potential drivers and traits shaping the causal paths to selection. Our current approach based on multi-model inference [START_REF] Grueber | Multimodel inference in ecology and evolution: challenges and solutions[END_REF] proved to be powerful to disentangle the different months and modes of temperature that are related to the strength of selection.

Rather than analyzing a general relationship between selective forces and temperature, we explored temporal changes in a variety of temperature modes and seasonality, which could affect spring phenology in several ways. By showing a warming close to zero in winter, our study confirms results obtained in previous studies in Switzerland [START_REF] Rebetez | Monthly air temperature trends in Switzerland 1901-2000 and 1975-2004[END_REF], UK [START_REF] Prior | Analyses of trends in air temperature in the United Kingdom using gridded data series from 1910 to 2011[END_REF] and southern France [START_REF] Chaouche | Analyses of precipitation, temperature and evapotranspiration in a French Mediterranean region in the context of climate change[END_REF]. Trends in temperatures during the blue tit breeding season (around 0.6°C per decade from March to June) were similar to trends found previously for the same region during the same period (+0.7°C per decade) [START_REF] Chaouche | Analyses of precipitation, temperature and evapotranspiration in a French Mediterranean region in the context of climate change[END_REF]. Such temperature warming in spring and summer provides scope for potentially important mismatches across tropic levels [START_REF] Harrington | Climate change and trophic interactions[END_REF].

Although this is not a causal relationship, April maximum temperatures were the main variable impacting the strength of selection on laying date. The causal path behind this relationship is certainly linked to the date of emergence of the main food resource of blue tit nestlings (i.e. the caterpillars). Indeed, because the caterpillars hatching is tracking April temperatures (through the bud burst oaks), blue tits are desynchronised with their main food resource, and selection for early laying is occurring in a warming context. Such mistiming has been shown in great tits in the UK [START_REF] Charmantier | Adaptive Phenotypic Plasticity in Response to Climate Change in a Wild Bird Population[END_REF]) and the Netherlands [START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF], where birds could not track the advancement in caterpillar emergence date. Moreover, spring warming could also decrease the caterpillar biomass by disrupting the synchrony between tree and caterpillar phenology, as shown in European oaks (Quercus robur, [START_REF] Visser | Warmer springs disrupt the synchrony of oak and winter moth phenology[END_REF] and in Sitka spruces (Picea sitchensis, [START_REF] Dewar | Predicted changes in the synchrony of larval emergence and budburst under climatic warming[END_REF].

In this context, it would be relevant to explore the effect of warming seasonality on the caterpillar emergence date and its relationship with selection for early breeding in blue tits.

Unfortunately, only 10 years of data on caterpillar biomass were available in this study site, preventing us from conducting such a long-term analysis.

Maximum April temperatures anomalies increased significantly since 1970 in the forest of La Rouvière (slope = 0.056 ± 0.012, P < 0.001), yet we did not detect any significant pattern in selection gradients across years since 1991 (Fig. 2). This discrepancy is due to the fact that spring temperatures, while getting warmer, did not significantly increase during the 1991- modulating selection pressure remains associated with spring warming [START_REF] Charmantier | Climate change and timing of avian breeding and migration: evolutionary versus plastic changes[END_REF]. Indeed, a desynchronization between the caterpillar peak and the anomalies in spring temperature could modify the evolutionary trajectory of the population, making future predictions challenging.

Several processes other than warming could influence the strength of selection on the timing of reproduction of blue tits and shape its evolutionary trajectory. For instance, the selection gradient acting on laying date was lower under higher annual predation rate in our study site (see Table 1). Although the difference in mean absolute fitness was non-significant between high and low predation years (F = 2.87, P = 0.09), the variance in relative fitness was significantly (F = 1.70, P < 0.001) lower in high predation (=0.064) than in low predation years (=0.109), which could potentially explain the reduction in the magnitude of selection in high predation years. This result, however, is only based on four years of high predation rates and thus needs to be confirmed and refined by a longer term study of predation and selection in order to refine the biologically interpretation of this relationship. Additionally, while climate warming is mainly characterized by an increase in mean temperature, an increase in climate variability among years is also observed and predicted in the future [START_REF] Coumou | A decade of weather extremes[END_REF], especially in the Mediterranean area [START_REF] Giorgi | Climate change projections for the Mediterranean region[END_REF]. Hence, an increase in the variability of selection regimes among years could occur in the future, and could shape the adaptive response to warming. For instance, Kingsolver and Buckley (2003) found that the predicted evolutionary response to mean climate warming is small in an alpine butterfly (Colias meadii) studied since 1980, mainly because of the variability in selection.

While climatic variability is increasing along with warming, the impact of this variability on evolutionary response will be difficult to predict, especially in term of extreme climatic events (Bailey & van de Pol 2016b). For instance, a previous study in the same population showed that when 10% of broods experienced an extremely hot day, selection for earlier breeding increased by 39% [START_REF] Marrot | Multiple extreme climatic events strengthen selection for earlier breeding in a wild passerine[END_REF]). In addition, besides global warming, climate change is also characterized by changes in precipitation regimes. Indeed, heavy rainfalls are expected to be more frequent [START_REF] Giorgi | Climate change projections for the Mediterranean region[END_REF], which could impact the amount of caterpillar available for tits during the breeding season. We did not explore the influence of monthly precipitations on selection because no change in precipitation regime was detected in our study site since 1975 (unpubl. data). Another possible limitation of our study is the fact that selection was investigated using fledging success as a fitness proxy.

Although laying date mainly affects the number of fledglings [START_REF] Perrins | The timing of bird' breeding seasons[END_REF], climate warming could also influence the selection acting on nestling survival until recruitment and/or adult survival [START_REF] Grosbois | Climate impacts on Mediterranean blue tit survival: an investigation across seasons and spatial scales[END_REF]). However, nestling recruitment is linked to winter survival and although the investigation of the effect of winter climate on number of recruits and selection would be interesting, this would require to take into account 36 additional climatic variables in our models (three climatic variables per month), leading to challenging interpretation of results. Moreover, although number of recruits is more closely related to the true fitness than fledging success is, the number of fledglings is more tightly linked with adult reproductive performance than number of recruits (which depends more on the offspring phenotypes).

To conclude, this study reveals that climate warming induces an increase in selection pressure during reproduction for a blue tit population in the Mediterranean region. This finding implies that recent climate change results in strong selection favouring specific phenotypes (e.g. earlier laying dates). This means that over and above the common phenological plasticity observed [START_REF] Charmantier | Climate change and timing of avian breeding and migration: evolutionary versus plastic changes[END_REF], natural selection could allow populations to cope with global warming, if the lag between the increase in selection strength and the evolutionary response remains short [START_REF] Hendry | Whither adaptation?[END_REF]. While increasing natural selection under climate change could be a route for populations to adapt to climate change, it could also be an A:

  laid, March 1 st = 1), clutch size (number of eggs laid) and number of fledglings (number of chicks leaving the nest) are recorded. Parents are captured in nest boxes when chicks are c.9 days or older, and are uniquely marked with metal rings (provided by C.R.B.P.O). Nestlings are ringed when 9-15 days old. Only first clutches that were not manipulated were included in the analyses (second clutches represent less than 1% of total number of clutches). Bird capture, ringing and measurement were performed under individual ringing permits delivered by the Centre de Recherche sur la Biologie des Populations d'Oiseaux, C.R.B.P.O, Paris (France). The monitoring protocol was approved by the Animal Care and Use Committee Languedoc-Roussillon (CEEA-LR-12066).

  2015 period (slope = 0.051 ± 0.035, P = 0.16), arguably due to a smaller sample size for this period (n = 24 versus 45 years). If the warming mechanism involved until now remains stable over time, average maximum April temperature anomalies, compared to the 1970-2000 period, are predicted to reach +3°C in 2050 (predictions extracted from an autoregressive integrated moving average model conducted on the 1970-2015 period, see FigS4). Such spring warming might reinforce the strength of selection acting on laying date by 139% (a multiplication by almost 2.5), meaning that differences in fledging success between blue tits will keep increasing depending on their laying dates. However, these results must be considered with caution because the population should theoretically respond by advancing its breeding timing (a heritable trait in this population, h²= 0.11[START_REF] Delahaie | Conserved G-matrices of morphological and life-history traits among continental and island blue tit populations[END_REF]. Moreover, mean laying date is expected to advance only if the causal factor (the caterpillar peak)

Figure S1 :

 S1 Figure S1: Daily (A) mean, (B) minimum and (C) maximum monthly temperature anomalies over the 1970-2015 period compared to the 1970-2000 period. The red line displays a smooth trend based on a 2-year moving window.

Figure S2 :

 S2 Figure S2: Distribution of relative fitness (i.e. number of fledglings) of female blue tits in La Rouvière, southern France.

Figure S3 :

 S3 FigureS3: Standardized selection gradients on laying date for anomalies in daily maximum April temperatures ranging between 0 and 3°C. Selection gradients were extracted from the best model of the model averaging approach (including an interaction between laying date and predation rate and maximum April temperature). The mean annual standard deviation of laying date was 7.6 days in our population.

Figure

  Figure S4: (A) ARIMA procedure. The graph show temporal autocorrelation in April maximum temperature. Dotted blue lines represent significance thresholds. Temperatures were autocorrelated at time lags of 3, 4 and 5 years. (B) Anomalies in daily maximum April temperatures over the 1970-2050 period compared to the 1970-2000 period. Predictions (blue line) for the 2015-2050 period were extracted from an autoregressive integrated moving average model including a temporal autocorrelation order = 4. The grey area shows the 95% confidence intervals.

  [START_REF] Visser | Warmer springs lead to mistimed reproduction in great tits (Parus major)[END_REF][START_REF] Husby | Speeding Up Microevolution: The Effects of Increasing Temperature on Selection and Genetic Variance in a Wild Bird Population[END_REF][START_REF] Reed | Phenological mismatch strongly affects individual fitness but not population demography in a woodland passerine[END_REF][START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF] and the United Kingdom (e.g.[START_REF] Charmantier | Adaptive Phenotypic Plasticity in Response to Climate Change in a Wild Bird Population[END_REF]. It is informing to compare the strength of the relationship between the directional selection on laying date

	and spring temperature among studies. This relationship, termed the thermal sensitivity of
	selection (B, Chevin, Lande & Mace 2010), has been estimated at = -5.30 days/°C for great
	tits in Wytham Woods (UK, Vedder, Bouwhuis & Sheldon 2013) and at -5.01 days/°C for the
	same species in Hoge Veluwe (Netherlands, Chevin, Visser & Tufto 2015). Considering our
	estimate of the interaction term between maximum April temperatures and laying date (=-
	0.024 ± 0.013, Table

Table 1 :

 1 Summary results of the model averaging approach to assess climatic variables driving selection on laying date.T°C min, T°C mean and T°C max represent the daily minimum, mean and maximum monthly temperatures respectively. Relative importance represents the probability that the variable was included in the best model. In bold, the only variable with a relative importance above 70%. Population density represents the annual mean distance between an occupied nest box and its closest neighbour. Predation rate was estimated as a categorical variable taking two different modalities: low and high corresponding to a predation rate below or above 0.1 respectively.

	Variable in interaction with laying date	Averaged estimates	Relative importance
	T°C min -March	0.0000 (± 0.0006)	1 %
	T°C mean -March	0.0000 (± 0.0004)	< 1 %
	T°C max -March	0.0000 (± 0.0004)	< 1 %
	T°C min -April	0.0000 (± 0.0000)	< 1 %
	T°C mean -April	-0.0016 (± 0.0068)	6 %
	T°C max -April	-0.0244 (± 0.0108)	89 %
	T°C min -May	0.0000 (± 0.0008)	1 %
	T°C mean -May	0.0002 (± 0.0022)	2 %
	T°C max -May	0.0013 (± 0.0058)	6 %
	T°C min -June	0.0000 (± 0.0011)	1 %
	T°C mean -June	0.0001 (± 0.0013)	1 %
	T°C max -June	0.0000 (± 0.0010)	1 %
	Population density	0.0001 (± 0.0015)	1 %
	Predation rate (Low)	-0.0725 (± 0.0351)	88 %

Journal of Animal Ecology: Confidential Review copy Journal of Animal Ecology: Confidential Review copy

Journal of Animal Ecology: Confidential Review copy

Acknowledgements

We thank Samuel Caro, Jarrod Hadfield, Ben Sheldon and an anonymous reviewer for useful comments on previous versions of this manuscript. We are very grateful to Fanie Pelletier, Mark Vellend, Hélène Fréville and Céline Teplitsky for useful discussions. We thank Cyril Bernard and the Forecast team at the CEFE for providing climatic data. We thank everyone who helped maintain the study site and conduct the blue tit monitoring since 1991, in particular Philippe Perret, Marcel Lambrechts, Arnaud Grégoire and Claire Doutrelant. This project was funded by the Agence Nationale de la Recherche (BioAdapt grant ANR-12-ADAP-0006-02-PEPS to A.C.), the European Research Council (Starting grant ERC-2013-StG-337365-SHE to A.C.), the OSU-OREME, and the Natural Sciences and Engineering Research Council of Canada (NSERC discovery grant to D.G.).

Journal of Animal Ecology: Confidential Review copy Journal of Animal Ecology: Confidential Review copy indicator of increased extinction risk since it suggests that part of the population has reduced fitness [START_REF] Bell | Evolutionary rescue and the limits of adaptation[END_REF].

Journal of Animal Ecology: Confidential Review copy Journal of Animal Ecology: Confidential Review copy Directional selection gradients for laying date (solid lines) estimated on three datasets containing only years experiencing three different warming context: (A) 12 years experiencing anomalies in maximum April temperatures lower than 1°C (n = 530 broods); including 4 years (1994, 1998, 2000 and 2013) associated with negative temperatures anomalies, (B) 3 years experiencing anomalies in maximum April temperatures ranging from 1°C to 2°C (n = 104 broods) and (C) 5 years experiencing anomalies in maximum April temperatures higher than 2°C (n = 288 broods). Temperatures anomalies correspond to the temperatures compared to the average temperatures over the 1970-2000 period. Years associated with a high predation rate (2001, 2006, 2009 and 2010) were removed from these datasets to avoid confounding their effect with the climatic signal (see Table 1). Dashed lines represent 95% CIs. The mean annual standard deviation of laying date was 7.6 days in our population.

677x381mm (72 x 72 DPI)

B:

ARIMA procedure:

First, we tested for temporal autocorrelation in temperatures using the acf function in R. Because significant temporal autocorrelation was detected in the temperature series (based on the visualization of the correlogram), an autoregressive model structure was needed.

Commonly used in temporal series analysis, autoregressive models are a special case of statistical models, allowing to take into account temporal autocorrelation at several lags, and integrating a moving average [START_REF] Brockwell | Introduction to Time Series and Forecasting[END_REF]. Because several ARIMA are possible for a time series, we used the auto.arima function from the package forecast [START_REF] Hyndman | Automatic Time Series Forecasting: The forecast Package for R[END_REF] in R, in order to select the best ARIMA describing temperatures changes. The best ARIMA included a degree of differencing of 1 (which allows to take into account the positive trend) and a negative autoregressive order of 4, meaning that temperatures identified as a driver of selection were negatively autocorrelated for up to four years.

Table S1: 95% confidence set of models. These models were extracted from the 95% confidence set of best models (see text). Because laying date, clutch size, as well as individual and years random intercept were constrained to be included in all models, they are not displayed. ∆AIC and w i represent the difference between each model and the best model and the Akaike weights respectively.