Integrating Causal and Evolutionary Analysis of Life-History Evolution: Arrival Date in a Long-Distant Migrant

In migratory species, the timing of arrival at the breeding grounds is a life-history trait with major fitness consequences. The optimal arrival date varies from year-to-year, and animals use cues to adjust their arrival dates to match this annual variation. However, which cues they use to time their arrival and whether these cues actually predict the annual optimal arrival date is largely unknown. Here, we integrate causal and evolutionary analysis by identifying the environmental variables used by a migratory songbird to time its arrival dates and testing whether these environmental variables also predicted the optimal time to arrive. We used 11 years of male arrival data of a pied flycatcher population. Specifically, we tested whether temperature and normalized difference vegetation index (NDVI) values from their breeding grounds in the Netherlands and from their wintering grounds in Ivory Coast explained the variation in arrival date, and whether these variables correlated with the position of the annual fitness peak at the breeding grounds. We found that temperature and NDVI, both from the wintering and the breeding grounds, explained the annual variation in arrival date, but did not correlate with the optimal arrival date. We explore three alternative explanations for this lack of correlation. Firstly, the date of the fitness peak may have been incorrectly estimated because a potentially important component of fitness (i.e., migration date dependent mortality en route or directly upon arrival) could not be measured. Secondly, we focused on male timing but the fitness landscape is also likely to be shaped by female timing. Finally, the correlation has recently disappeared because climate change disrupted the predictive value of the cues that the birds use to time their migration. In the latter case, birds may adapt by altering their sensitivity to temperature and NDVI.

In migratory species, the timing of arrival at the breeding grounds is a life-history trait with major fitness consequences. The optimal arrival date varies from year-toyear, and animals use cues to adjust their arrival dates to match this annual variation. However, which cues they use to time their arrival and whether these cues actually predict the annual optimal arrival date is largely unknown. Here, we integrate causal and evolutionary analysis by identifying the environmental variables used by a migratory songbird to time its arrival dates and testing whether these environmental variables also predicted the optimal time to arrive. We used 11 years of male arrival data of a pied flycatcher population. Specifically, we tested whether temperature and normalized difference vegetation index (NDVI) values from their breeding grounds in the Netherlands and from their wintering grounds in Ivory Coast explained the variation in arrival date, and whether these variables correlated with the position of the annual fitness peak at the breeding grounds. We found that temperature and NDVI, both from the wintering and the breeding grounds, explained the annual variation in arrival date, but did not correlate with the optimal arrival date. We explore three alternative explanations for this lack of correlation. Firstly, the date of the fitness peak may have been incorrectly estimated because a potentially important component of fitness (i.e., migration date dependent mortality en route or directly upon arrival) could not be measured. Secondly, we focused on male timing but the fitness landscape is also likely to be shaped by female timing. Finally, the correlation has recently disappeared because climate change disrupted the predictive value of the cues that the birds use to time their migration. In the latter case, birds may adapt by altering their sensitivity to temperature and NDVI.

INTRODUCTION
In seasonal habitats, the timing of annual cycle stages such as flowering and bud burst in plants, and reproduction, molt or migration in animals, is matched with a short window in which conditions are optimal for these stages. Organisms often suffer major fitness costs if their annual cycle stages occur outside this "optimal window" and hence, in several species, the timing of these stages is often under strong stabilizing selection (Perrins, 1970;Noordwijk et al., 1995;Brown and Brown, 2000;Bêty et al., 2004;Smith and Moore, 2005). Examples range from birds (Perrins, 1970;Noordwijk et al., 1995;Brown and Brown, 2000;Bêty et al., 2004;Smith and Moore, 2005) to fish (Einum and Fleming, 2000), plants (Kelly and Levin, 1997;Donohue et al., 2005;Koenig et al., 2012), and invertebrates (Philippart et al., 2003). Optimal time windows may be very narrow and vary in date and/or duration from one year to the other. For instance, the timing of the optimal window for breeding may be related to the moment when food abundance is at its peak (Noordwijk et al., 1995;Visser et al., 2006), while the local climate would be an important determinant of optimal arrival date at the breeding grounds in migratory birds (Brown and Brown, 2000;Erni et al., 2005).
Matching the timing of annual cycle stages with the year-toyear variation in the optimal time window is not straightforward. The timing of annual cycle stages often has to be made well in advance of the time when the stages are under selection, since many of these stages require a period of preparation (Noordwijk and Müller, 1994;Visser et al., 2004Visser et al., , 2010. For instance, to take full advantage of an optimum in food availability during the reproductive stage, a bird must have chicks in the nest at this precise moment, but having chicks in the nest depends mainly on when eggs were laid, which happens several weeks earlier in the season (Visser, 2008(Visser, , 2013. Individuals therefore respond to environmental conditions that serve as so-called cues, which affect the timing of these stages. Importantly, these cues are only useful if they are predictive, i.e., when they correlate with or causally affect the fitness landscape for that annual cycle stage (Noordwijk and Müller, 1994;Visser et al., 2004; Figure 1).
Migratory animals are particularly interesting to study the phenology of key annual cycle stages, as their decisions on the timing of some stages need to be made not only at an earlier time than the optimal time window, but also at a different geographical location. Long-distance migrants overwinter far away from their breeding grounds, and hence cannot directly sense the conditions they will face upon arrival and during breeding (Both and Visser, 2001;Visser and Both, 2005). Photoperiodic variation is considered to be the major cue that these animals use to time their migration (Gwinner, 1987(Gwinner, , 1989aDawson, 2002;Coppack and Pulido, 2004;Pulido, 2007). However, the photoperiodic regime is invariable among years and cannot predict annual variation in optimal conditions, making it unreliable to fine-tune the timing of migration to the among-year variation in timing of the optimal window.
The fine-tuning of migration timing is related to environmental conditions (Pulido, 2007). For instance, some birds arrive earlier on their breeding grounds when temperatures are higher before departure from non-breeding areas during spring migration Marra et al., 2005). Moreover, weather conditions can not only affect departure decisions (Gordo, 2007;Bauer et al., 2008;Eikenaar and Schmaljohann, 2014;Deppe et al., 2015;Ouwehand and Both, 2017) but en route decisions as well, leading to minor adjustments in timing of spring arrival (Hüppop and Hüppop, 2003;Hüppop and Winkel, 2006;Bauer et al., 2008;Both, 2010;Haest et al., 2018). It is poorly known if the environmental conditions in the overwintering grounds that affect arrival date also correlate with the conditions at the breeding grounds. However, if they do correlate, then it is possible for environmental variables at the wintering grounds to also predict the annual variation in fitness peaks at the breeding sites. Some studies suggest that large-scale correlations in climatic conditions exist between Europe and Africa Saino and Ambrosini, 2008). This means that individuals may acquire information (e.g., temperature, weather) before and during migration to predict the conditions at the breeding grounds, and thus time their arrival in relation to the optimal time window (Bauer et al., 2008;Both, 2010;Tøttrup et al., 2010).
Identifying the climatic variables that allow birds to predict the fitness peak later in the season is relevant in light of current environmental changes. It is known that many bird species have advanced their migration dates, suggesting that migrants are able to adjust their migration timing to the environmental conditions, but this pattern is not consistent across species (Both and Visser, 2001;Cotton, 2003;Marra et al., 2005;Gienapp et al., 2007;Usui et al., 2016;Mayor et al., 2017). In contrast, the relationship between advancements in timing of migration and shifts in optimal conditions to migrate is poorly known since individual fitness data are rarely available . If the environmental variables the animals use for timing their migration merely correlate with the drivers of selection, rather than causally affecting the fitness landscape, climate change may decrease the reliability of these cues by disrupting these correlations, as there is both spatial and temporal variation in the rate of climate change (Easterling et al., 1997;Vose et al., 2005;Stocker et al., 2013).
Here we analyzed 11 years of data of individual male arrival date and breeding success of a long-distance migratory passerine, the European pied flycatcher (Ficedula hypoleuca) to understand which environmental variables serve as cues for the timing of arrival, and whether these cues are indeed predictors of the annual variation in the optimal arrival date (i.e., the fitness landscape). For that, we explored temperature and the normalized difference vegetation index (NDVI) in the breeding and wintering grounds as potential cues and predictors of optimal arrival date.

Study System and Study Area
Pied flycatchers [Ficedula hypoleuca ( [Pallas], 1764) Muscicapidae] are small long-distance migratory birds that breed in Europe and winter in Africa. Due to their propensity to use nest boxes, it is also possible to precisely monitor their breeding. Long-term data collection was conducted in the southern forests of the Hoge Veluwe National Park (Netherlands; 52 • 02 07 N 5 • 51 32 E). Voucher material of this pied flycatcher population was deposited in the ornithology collection of the Naturalis Biodiversity Center (Leiden, Netherlands) under the inventory numbers RMNH 592347, RMNH 592348, and RMNH 592349.

Male Arrival Dates
Male individual arrival data were collected from 2005 to 2015 (11 years). Male flycatchers settle soon upon arrival in a territory from which they are rarely displaced, and advertise a cavity FIGURE 1 | Link between environmental variables that are correlated with year-to-year variation in arrival date (cues) and their effect on the fitness consequences of arrival date: The year-to-year variation in fitness depends on the interaction between the environmental drivers of selection and individual decisions. Because decisions on arrival date are not made in the same (spatial and temporal) environment of selection, animals use cues to predict conditions of the environment at the time of selection (adapted from Noordwijk and Müller, 1994;Visser et al., 2004Visser et al., , 2010 or nest box to the females by singing continuously at or close to the potential nest site, where they can eventually breed (Potti, 1998;Both et al., 2016). Their advertising behavior made them very conspicuous, which provides very high daily detection probabilities (>80%; Both et al., 2016). Individual arrival date of males was assessed by daily field observations of singing males in the study area from early April until the first week of May. These field observations were made by one to three (typically two) trained observers who walked independently pre-established routes covering the whole study area and visiting all nest boxes. Routes and direction of the routes were alternated daily among observers in order to prevent biases due to among-observer variation in detection probability and accuracy. Although marking all male birds with unique color ring combinations was logistically impossible (due to potential leg injury and to the relatively large population size), we used field descriptions of male plumage characteristics and the presence/absence and position of one metal and one (sometimes two) color rings to establish whether a bird initially singing in a particular nest box was indeed the same that ended up breeding in it (see Potti, 1998;Both et al., 2016 for similar approaches). In the Netherlands, male pied flycatchers display a relatively large individual variation in plumage characteristics, such as the color of dorsal feathers (from female-like brown to almost fully black) and the size and shape of their forehead patch (see Drost, 1936;Lundberg and Alatalo, 1992;Both et al., 2012 and Supplementary Material). This method allows us to assign with high certainty the arrival date to more than half of the singing males in the study area (Visser et al., 2015). Likewise, the reliability of this approach is supported by a strong correlation (r = 0.92) between the arrival dates estimated by visual observations and the arrival dates obtained from geolocator data (Both et al., 2016).
During the chick-rearing phase (see below), males were caught, described again in terms of plumage characteristics and ringed. In most years (but not in 2005 and 2012) we also collected data on "bachelors, " i.e., males that were unpaired when most females had arrived. Those males were captured and identified when they were inspecting and advertising their potential nest site using spring traps installed at the entrance of the nest boxes. Some of those presumed bachelors became breeding birds later in the season. Arrival data of these bachelor males were also used in the analyses except for number of fledglings and recruits.

Individual Breeding Success
Apart from data on timing of arrival, we also collected information on individual breeding success using a standardized nest monitoring protocol (Visser et al., 2015). For all years in which we had arrival data we collected data on whether the male obtained a social female or not (i.e., probability of obtaining a mate, as described above), number of chicks successfully raised (i.e., number of fledglings), and whether individual offspring recruited as a breeding bird (a so-called recruit) in the following year(s) or not. Female pied flycatchers only have a single brood per season, with cases of second broods being rare , but males may attract and pair with a second female after their first female starts to incubate. In those few cases of polygyny in which males were paired with two females with similar laying dates (33 out of 676 observations), both broods were included.

Environmental Variables
Data on Dutch daily mean temperatures were obtained from the Dutch meteorological institute database (KNMI -https:// www.knmi.nl/nederland-nu/klimatologie/, accessed on February 2016). For these, we used data from the Deelen weather station which is directly adjacent to the study area. For African temperatures, we obtained daily air gridded temperatures for Ivory Coast from the NCEP Daily Global Analyses data provided by the NOAA/OAR/ESRL PSD, Boulder, CO, United States, from their Web site at https://www.esrl.noaa.gov/psd/ (Accessed on January 2018). Previous geolocator data pinpoint Ivory Coast as the main wintering location of our study population (Tomotani et al., 2019), which is additionally supported by published data from other geographically close pied flycatcher populations . We then obtained NDVI data from the MODIS/Terra Vegetation Indices 16-Day L3 Global 0.05 Deg CMG V006 (Didan, 2015) and extracted average values for the whole Ivory Coast and for the breeding location at the Hoge Veluwe National Park in the Netherlands. We opted for only using wintering and breeding grounds temperatures because en route temperatures (timing and location) would be highly dependent on individual migration routes and timing of migration, and this information was not available.
Finally, we obtained data on the photoperiodic variation of the Netherlands from the National Oceanic and Atmospheric Administration, United States (NOAA, https://www.noaa.gov/ accessed on February 2016). We considered the civil twilight as the boundary of the effective light phase relevant for the birds (see, Gwinner, 1989a). Because photoperiodic variation on the breeding and wintering grounds has a high correlation, changing in the same direction at the same time, it was not necessary to obtain and model day-length data for Ivory Coast separately (Supplementary Figure 1).

Data Analysis
All analyses were performed in R version 3.4.3 (R Core Team, 2017). To define the minimal adequate models, in all cases we used backward model selection, i.e., sequentially dropping non-significant terms starting with interactions, retaining only significant terms. All statistics are reported at the point that the term was dropped from the model.

Cues Used for Arrival Date
To assess which environmental variables affect arrival date, we followed the method described in Gienapp et al. (2005) and used proportional hazard models (Cox, 1992) with random effects implemented in the R "survival" and "coxme" packages (Therneau, 2015). Proportional hazard models calculate the daily probability of an event to occur, and how this depends on explanatory variables, which can be time-dependent (i.e., variables that change their value during the time an individual is "at risk"). Using this approach is biologically more realistic than using fixed time windows over which these environmental variables are averaged (Gienapp et al., 2005). The proportional hazard approach also allows the inclusion of day-length effects as a proxy for time of year.
First, we defined the most relevant period over which the environmental variables affected the arrival dates. For both Dutch and Ivory Coast temperatures and NDVI values, the value at day t was calculated as the average over six periods of increasing length ending at day t (window duration from 5 to 30 days with a 5-day increment for each period). For example, for the 30th of April, we produced six windows with averaged temperatures and NDVI of the Netherlands and Ivory Coast: from days 25-April to 30-April (5 days), 20-April to 30-April (10 days), 15-April to 30-April (15 days), 10-April to 30-April (20 days), 5-April to 30-April (25 days), and 1-April to 30-April (30 days). In the case of Ivory Coast NDVI and temperatures, we also used lagged shifting windows of the same length (20 days) but ending 20-80 days before day t [based on the time taken during spring migration; ]. Thus, for the same 30th of April example, averages over 20 days were made for four additional periods: 21-January to 10-February (80-day lag), 10-February to 2-March (60-day lag), 2-March to 22-March (40-day lag), and 22-March to 11-April (20-day lag). In total, we compared 24 possible combinations of window durations (5, 10, 15, 20, 25, or 30 days) and window lags (20, 40, 60, 80 days) using the model log-likelihood (see Supplementary Table 1).
After defining the position and duration of the climate window, we tested, also with proportional hazard models, which variables significantly explained the variation in arrival date of the pied flycatchers. We fitted Ivory Coast (with and without lag) and Dutch temperatures, Ivory Coast (with and without lag) and Dutch NDVI and the two-way interactions between day-length and Dutch temperature, day-length and lagged Ivory Coast temperature and day-length and lagged Ivory Coast NDVI to investigate whether the effects of environmental variables differ throughout the season. In this analysis we also included characteristics of the birds that could potentially explain variation in arrival dates: age and the Drost score (variation in male color from 1, completely black males, to 7, brown males; Drost, 1936) which could be correlated with age, metabolic or behavioral differences (Ivankina et al., 2001(Ivankina et al., , 2007Kerimov et al., 2013). Given the difficulties of aging the breeding birds based on plumage characteristics, age was solely assigned based on ringing data records, distinguishing three categories: "2nd-calendar-year birds" (SCY), i.e., birds that hatched the year before (when they were ringed as chicks); "after second-calendar-year-birds" (ASCY), i.e., birds that were at least 2 years old; and "unknown, " i.e., for non-ringed breeders in which age cannot be reliably assigned). The interaction between age and the environmental variables was also included to test whether those variables affected younger and older birds in a different way. Individual was included as random effect to control for multiple observations of the same individual.

Fitness Consequences of Arrival Date
To test whether male fitness depended on arrival date, we regressed fitness against the arrival date in interaction with the environmental variable. To account for the expected nonlinear effect of fitness versus arrival date, we also included arrival date squared. A significant interaction of arrival date and the particular environmental variable indicates that the environmental variable correlates with the date of arrival where fitness is maximized in that year (i.e., the position of the fitness peak).
We used three fitness components to assess the fitness consequences of arrival date: the probability of obtaining a mate, number of fledglings of a male brood(s) and the probability that at least one of these fledglings recruited as a breeding bird in the following years (probability of obtaining a recruit). Tests were done with linear mixed effect models with year and individual as a random effects and "nloptwrap" as optimizer ("lme4" R package, Bates et al., 2015) for each environmental variable and each fitness component separately. We used a logit-link and Binomial error distribution for the probability of obtaining a mate and obtaining a recruit (both modeled as whether or not the male obtained a mate or recruit, respectively) and Poisson distribution for the number of fledglings. Model comparison was performed with parametric bootstrap with 2000 simulations ("PBmodcomp" implementation of the R package "pbkrtest, " Halekoh and Højsgaard, 2014). Sample sizes varied in the analysis of each fitness components, because some components were affected by one of the experimental manipulations done in this population over the past decades (see, for instance, Bauchau and Seinen, 1997;Tomotani et al., 2016).
Finally, we looked at whether the mean arrival date (observed) correlated with the optimal arrival date. In order to extract the latter, we fitted a model of individual fitness (probability of obtaining a mate or number of fledglings) against year, arrival date, arrival date squared, and the interaction between arrival date and year (with individual as random effect) and regressed the estimates for the year × arrival date term (which represents the position of the annual fitness peak) against the mean arrival date of each year.

Cues Used for Arrival Date
The most relevant period over which the environmental variables affected the arrival dates included Ivory Coast temperatures and NDVI with the window moved by a 60-day lag (thus starting, on average, around the 18th of February), as well as a window duration of 25 days for the non-lagged Ivory Coast and Dutch temperatures and NDVI (Supplementary Table 1).
Higher temperatures in both the Netherlands and Ivory Coast were related to earlier arrival; with a significant interaction between temperature and day-length, indicating that the temperature effect was weaker later in the season (coefficients indicate a change in the "hazard" for the event to occur: Netherlands: coefficient for the interaction = −0.01 ± 0.001, coefficient for temperature = 5.88 ± 1.1; Ivory Coast: coefficient for the interaction = −0.14 ± 0.03, coefficient for temperature = 132.09 ± 30.60; Table 1). Finally, the nonlagged Ivory Coast temperature was also related to an earlier arrival (coefficient = 0.55 ± 0.47; Table 1), but without the significant interaction with day-length (Table 1).
Higher NDVI in the Netherlands and higher lagged NDVI in Ivory Coast were related to later arrival (Netherlands: coefficient = −1.46 ± 0.26; Ivory Coast: coefficient = −0.56 ± 0.22; Table 1). On the other hand, higher NDVI in Ivory Coast in the month prior to arrival (not lagged) was related to earlier arrival (coefficient = 1.34 ± 0.17; Table 1).
Male arrival was also explained by age. After second calendar year (ASCY) males tended to arrive earlier than second calendar year birds (SCY; coefficient for SCY = −1.48 ± 0.49; Table 1). We also found a significant interaction between age and lagged Ivory Coast temperatures, with a stronger temperature effect in younger birds (coefficient for the interaction: SCY relative to ASCY = 5.63 ± 2.54; Table 1) as well as between age and lagged Ivory Coast NDVI, again, with a stronger effect in SCY (coefficient for the interaction: SCY relative to ASCY = 0.50 ± 0.70; Table 1). Thus, higher temperatures and NDVI values have a stronger effect on the arrival of younger birds (SCY) than of older birds (ASCY). Drost scores (male color) did not explain variation in arrival dates ( Table 1).

Fitness Consequences of Arrival Date
We observed a significant relationship between arrival date and both the probability of obtaining a mate ( Table 2) and the number of fledglings produced (Table 3), but found no relationship between arrival date and the probability of obtaining a recruit (Table 4). Late arriving males had a much lower probability of obtaining a mate than early arriving males ( Table 2 and Supplementary Figure 2A). The number of fledglings generally declined with arrival date but the decline was steeper in late-arriving males ( Table 3 and Supplementary Figure 2B). Yearly recruitment was low overall with 11-49 recruits per year.
Although we found a relationship between arrival date and fitness, when we tested whether the position of this fitness peak was affected by environmental variables, we did not find an effect: none of the environmental variables correlated with the position of the fitness peak regardless of the fitness measure used (represented by the interaction between arrival date and the environmental variable, Tables 2, 3). In accordance with this result, when we correlated the observed mean arrival date and the calculated optimal arrival date of a given year, we did not observe a relationship. This was the case for both when the probability of obtaining a mate and when the number of fledglings -the two fitness measures that were related to arrival date -were used (for the probability of obtaining a mate: estimate = −0.38 ± 1.84, F 1 , 6 = 0.04, p-value = 0.84; for the number of fledglings: estimate = −0.47 ± 0.46, F 1 , 9 = 1.06, p-value = 0.33; Figure 2). Therefore, there is a mismatch between optimal and mean arrival date and this is partly related to the very low variation in the optimal arrival date between years (Figure 2).

Temperature and NDVI values in the Netherlands and Ivory
Coast explained year-to-year variation in timing of arrival at the breeding grounds. These variables, however, did not correlate with the year-to-year variation in optimal arrival date, i.e., the arrival date at which fitness was maximized in that year. The lack of a relationship between the optimal and mean arrival dates suggests that this is not because we overlooked potential cues. These results are surprising: why would birds use cues for their arrival date that do not lead to an arrival date with a higher fitness return?
If there are no environmental variables that are predictive for the year-to-year variation in the conditions the birds will encounter when arriving at the breeding grounds, the best strategy for arrival date would be to arrive at the same time every year, at a date which on average over years is closest to the optimal arrival date (e.g., Gienapp et al., 2014). For example, the birds could make the arrival decision by relying solely on endogenous rhythms and/or responding directly to the photoperiodic changes in their wintering grounds (Gwinner, 1996), which show no yearto-year variation. But while repeatability of individual arrival dates is moderate in some species (e.g., pied flycatcher: 0.08-0.  Both et al., 2016 for a revision), there is still variation between years (Both et al., 2016). Our data show, however, that environmental variables from both the breeding and wintering grounds do explain year-toyear variation in arrival dates. Other studies have also reported effects of environmental conditions, both en route and at the breeding grounds, on the variation in timing of spring arrival, or, at least, the timing of passage dates, in the pied flycatcher (Ahola et al., 2004;Both et al., 2005;Saino et al., 2007). It is generally accepted that several migratory bird species migrate earlier in warmer springs (Gordo, 2007) and that timing of migration advances with climate change (Walther et al., 2002;Cotton, 2003;Gienapp et al., 2007, but see Both and Visser, 2001;Jenni and Kery, 2003;Tomotani et al., 2018). These studies and our results here suggest that there is plasticity in timing of arrival (but see Pulido, 2007), arguing against the inflexible strategy as the main (or only) mechanism. So, the question remains why the cues that are used by the pied flycatcher in their arrival date do not predict the timing of the fitness peak associated with arrival date. We discuss below three hypotheses: incomplete fitness measure, fitness benefits via females and disruption of cues due to climate change.

Incomplete Fitness Measure
A potential reason why we were unable to detect a correlation between environmental variables that are used as cues and their effect on the position of the fitness peak are our fitness measures. While we are able to precisely determine the breeding success of the birds that we were able to identify upon arrival, we have no data of birds that did not successfully complete their migration or of birds that completed the migration but bred elsewhere. While mortality during migration can be estimated from direct tracking (e.g., Klaassen et al., 2013;Rotics et al., 2016) or extensive re-sighting studies (e.g., Sillett and Holmes, 2002;Leyrer et al., 2013;Lok et al., 2015), data on mortality related to timing of migration, which is what is relevant for our research, is difficult, if not impossible, to obtain, since for obtaining an arrival date the individual must be alive (but see, e.g., Brown and Brown, 2000;Newton, 2007). If there is a departuredate-related mortality, the proper fitness measurement would be a measurement that includes all mortality events between departure from the wintering grounds and breeding (i.e., en route and arrival before breeding), and then the reproductive success, that are related to the departure decision. This could potentially shift the position of the fitness peak that we obtained only using A significant interaction between arrival date and the environmental variable indicates that the variable predicts the position of the fitness peak (i.e., the arrival date at which fitness is maximized). All variables are centered on the grand mean across all years. Statistics are given for each term at the point of exclusion (or inclusion if significant) of the term from the model. In case of significant interactions, statistics of main terms were obtained by removing the term in the presence of the interaction (indicated by an "i" in the "observation" column). P-values in bold represent significance after a Bonferroni correction to account for multiple testing (critical p-value = 0.05/number of tests). Arrival date is in centered April-day × 10 −1 , NDVI = centered NDVI × 10 −7 .ˆ2 = squared.
the reproductive success of the birds that did arrive (Jonzén et al., 2007;Lof et al., 2012).
There are few data on migration-timing-dependent mortality, particularly for small migrants such as pied flycatchers. The few existing examples show that poor weather conditions during migration, or when birds arrive at the breeding grounds, can heavily impact the early arriving birds thereby shaping the indicates that the variable predicts the position of the fitness peak. All variables are centered on the grand mean across all years. Statistics are given for each term at the point of exclusion (or inclusion if significant) of the term from the model. In case of significant interactions, statistics of main terms were obtained by removing the term in the presence of the interaction (indicated by an "i" in the "observation" column). P-values in bold represent significance after a Bonferroni correction to account for multiple testing (critical p-value = 0.05/number of tests). Arrival date is in centered April-day × 10 −1 , NDVI = centered NDVI × 10 −7 .ˆ2 = squared.
optimal timing of arrival (Brown and Brown, 2000;Newton, 2007). This acute effect of weather upon arrival could shift the fitness peak to a later date when spring is cold, if the low temperatures lead to a higher mortality for early arriving birds. Such effect of the temperature at the breeding grounds was previously shown for the pattern of recruitment of pied flycatcher chicks of the same population as the present study: selection for early arrival at the moment of recruitment to the breeding A significant interaction between arrival date and the environmental variable indicates that the variable predicts the position of the fitness peak. All variables are centered on the grand mean across all years. Statistics are given for each term at the point of exclusion (or inclusion if significant) of the term from the model. In case of significant interactions, statistics of main terms were obtained by removing the term in the presence of the interaction (indicated by an "i" in the "observation" column). P-values in bold represent significance after a Bonferroni correction to account for multiple testing (critical p-value = 0.05/number of tests). Arrival date is in centered April-day × 10 −1 , NDVI = centered NDVI × 10 −7 .ˆ2 = squared.
population was stronger in warmer than colder years, likely due to lower survival of early recruiting birds in cold springs (Visser et al., 2015). Thus, part of our missing fitness component could be explained by temperatures in the Netherlands between arrival and breeding for both experienced and first-time breeders. Another cue that explained the variation in arrival dates was the NDVI, which could indicate that the food availability at the wintering grounds and, in turn, the fueling and individual FIGURE 2 | Relationship between optimal and mean arrival dates. Optimal arrival date was extracted from the estimates of the model testing the relationship between probability of obtaining a mate (open dots) or number of fledglings (closed dots) and year, arrival and arrival squared, also including the interaction between year and arrival. Mean arrival date was calculated averaging the arrival dates of all males per year. Both optimal and mean arrival dates are centered on the grand mean. Each point represents 1 year. condition prior to departure play a role in the arrival timing. Timing of migration may indeed be dependent on the individual condition due to carry-over effects from good wintering habitats and this "condition" could potentially modify the relationship between the mortality and the departure date from the wintering grounds (Norris et al., 2004;Saino et al., 2004;Pulido, 2007;Studds and Marra, 2011). For example, higher NDVI could improve the condition of early migrants and lead to lower mortality of early birds shifting the fitness peak to an earlier date (due to a flatter relationship between mortality and departure). Indeed, our results show that higher lagged NDVI values correlate with earlier arrivals (Table 1).
Finally, high temperatures in Ivory Coast were also correlated to a later arrival, suggesting a similar role for the Ivory Coast temperatures as the NDVI. Tøttrup et al. (2012), for example, reported that droughts in the wintering grounds caused delays in arrival date of songbirds; the reason for this was a prolonged stay in stopover sites, suggesting that these birds had a slower fueling rate. Thus, temperatures at the wintering grounds have a similar potential to affect the slope of the mortality and departure relationship thereby shifting the optimal arrival date.
The absence of a correlation between the cues that explain the variation of arrival dates and the position of the fitness peak could, therefore, be caused by these various ways that the unmeasurable date-dependent mortality between departure of the wintering grounds and reproduction at the breeding grounds shift the fitness peak.

Fitness Benefits via Females
This study focused on the causes and fitness consequences of variation in male arrival time, using three fitness components: the probability of obtaining a mate, the number of fledglings and the probability to obtain recruits. The probability of a male to find a mate is related to male arrival as early arriving males will be more likely to attract a female (Canal et al., 2012;Tomotani et al., 2017). On the other hand, the number of fledglings and probability of obtaining recruits are fitness components much more related to females and thus the fitness benefits of early arriving for males must be because for instance early arriving males pair with early laying females (Saino et al., 2004;Cooper et al., 2010;Tomotani et al., 2017) or occupy higher quality territories. Thus, the fitness consequences of male arrival time may (partly) arise through correlation with female timing, especially female breeding time, which is related to reproductive success and to local environmental conditions (Dunn and Winkler, 1999;Townsend et al., 2013;Visser et al., 2015;Tomotani et al., 2018;Shave et al., 2019).
For the purposes of this study, female timing can be seen solely as a part of the causal chain by which male arrival is linked to fitness. It is likely that the male fitness landscape at the breeding grounds is, at least, partly shaped by the timing of breeding of the females that are available at male arrival.

Disruption of Cues Due to Climate Change
Another potential reason for the lack of a correlation between the cues used for arrival date and their effects on the fitness peak date could be that cues have become unreliable due to climate change. The environmental variables that the birds use to time their arrival are certainly not the ones that determine the optimal arrival date, but they only correlate to these causal variables that determine the optimal arrival date (the drivers of selection). Climate change may have altered the correlation between the cues and causal variables, especially because climate change has not affected temperatures in a uniform manner across different regions and moments of the year (Easterling et al., 1997;Vose et al., 2005;Serreze and Francis, 2006). Climate change has also disrupted the relationship between climatic and weather variables at large spatial scales, so called "teleconnections" (Diaz Henry et al., 2002). Consequently, the cues used by the birds to time their migration may actually have become unreliable over time and therefore we now no longer find a correlation between the actual and the optimal arrival dates (Figure 2). Adaptation to climate change, in this case, would require a change in the bird sensitivity (slope of plasticity of arrival date) to the cues so that a correlation is re-established.

Conclusion
The lack of adjustment in timing of breeding and/or migration in some pied flycatcher populations but not in others (Both and Visser, 2001;Ahola et al., 2004;Both et al., 2016;Tomotani et al., 2018) reflect the need of detailed studies on the predictability and use of cues across populations. In the present study we show that the timing of arrival of male pied flycatchers is a plastic trait as it is affected by environmental variables. Our data suggest that, to assess the reliability of cues as predictors for fitness, it is perhaps not sufficient to only use breeding success as fitness measurement, since information on how females shape the fitness landscape may also be important. Alternatively, cues that used to be predictive for the drivers of selection may no longer be so now due to climate change. If this is the case, we would expect selection on the sensitivity to the cues and, if this sensitivity is heritable, pied flycatchers could adapt to this disruption of cues and drivers of selection. Therefore, a better understanding of the mortality related to migration timing is essential to distinguish between these scenarios, allowing a more complete evaluation of the impacts of climate change on migratory animals.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Experimental Committee of the Royal Netherlands Academy of Arts and Sciences (KNAW).

AUTHOR CONTRIBUTIONS
BT, PG, IH, and MV: conceptualization. BT, PG, IH, FP, and MV: methodology. BT, IH, FP, and MT: investigation. BT and PG: formal analyses. BT and MV: data curation, project administration, and funding acquisition. BT, PG, and MV: writing -original draft. BT, PG, IH, MT, FP, and MV: writingreview and editing. PG and MV: supervision. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We are grateful to the board of the National Park "De Hoge Veluwe" for the permission to conduct our research on their property. Chiel Boom and Liam Bailey for the assistance with obtaining and cleaning the climatic data. Finally, we thank numerous researchers, technicians, students, and volunteers who helped to collect and manage the detailed data on arrival and breeding of pied flycatchers.