Advancing plant phenology causes an increasing trophic mismatch in an income breeder across a wide elevational range

. To avoid trophic mismatches, large herbivores should reschedule the breeding timing to keep in synchrony with the recent advance of plant phenology caused by climate change. To do so, individuals must track resources along two main non-exclusive axes: time (rescheduling the breeding timing) and space (i.e., shifting their range). While the vast majority of studies have focused on the study of speci ﬁ c populations inhabiting small and homogeneous environments, little is known about the response of species across large regions with highly heterogeneous topographies. Here, we studied roe deer ( Capreolus capreolus ), an income breeder that has demonstrated a lack of phenotypic plasticity on birth timing in a densely forested lowland area, and evaluated its ability to track resources along a wide elevational gradient (288 – 2366 m a.s.l.). We used data from 8986 parturition events collected in Switzerland during the period 1971 – 2015. We hypothesize that the mismatch between roe deer parturition dates and peak resource availability will increase over the study period. In addition, we expect this mismatch to be larger in low areas compared to higher elevation areas. Contrary to previous studies, we did ﬁ nd a consistent, but small, trend toward more advance parturition dates through time across all elevations. However, this rate of advance was less than that of the plant phenology indicators, resulting in an increasing mismatch at all elevations. Parturition dates changed on average at a rate between 7.5 times slower than the growing season start and 5 times slower than the ﬂ owering start. Thus, at the lowest elevations, parturition timing has already fallen outside the optimal time window for high-quality forage, while at high elevations, parturitions are still synchronized with the peak availability of quality forage. Our study emphasizes the need for large-scale study designs along climatic gradients to help understand the complex relationships that shape adaptation across environments. This would allow for better understanding of how populations will respond to the new environmental conditions imposed by climate change.


INTRODUCTION
Rising temperatures as a consequence of climate change have led to an overall spring phenology advancement for many organisms across the Northern Hemisphere (Menzel 2003, Cleland et al. 2007, Roberts et al. 2015. However, the magnitude of this response might differ among interacting species (Parmesan andYohe 2003, Root et al. 2003). Therefore, phenological mismatches may occur if a different magnitude of response disrupts previously synchronous trophic interactions (Visser et al. 1998, Koh et al. 2004. A decoupling between the date of key lifehistory events and resource availability may lead to strong negative consequences on species demography and population viability (Brook et al. 2009, Bronson 2009, Miller-Rushing et al. 2010, Usui et al. 2017.
For plant consumers, such as large herbivores, breeding timing is the period of peak energy demand for females which should match the optimal time in terms of resource availability (Post et al. 2003(Post et al. , H€ am€ al€ ainen et al. 2017). Otherwise, a mismatch can reduce annual reproductive success (Clutton-Brock et al. 1987, Plard et al. 2014, Paoli et al. 2018) and future reproductive performance of both females and newborns (Clutton-Brock et al. 1983, Festa-Bianchet et al. 2000. For most large herbivores, the leading edge of the growing season has been defined as the optimal timing to give birth (Merkle et al. 2016). The forage maturation hypothesis states that the early vegetative stages are nutritionally superior to later stages of plant development such as flowering, despite being lower in biomass production Langvatn 1992, Hebblewhite, Merrill, andMcDermid 2008). During flowering and seed production, plants produce fewer leaves and more stems, increasing the fiber content and reducing the protein content, and thus reducing their quality (Albon and Langvatn 1992, Duru 1997, Bumb et al. 2016. Therefore, early development stage with low vegetative biomass is nutritionally superior to mature, highbiomass vegetation (Fryxell 1991).
In order to track the rapid changes in plant phenology due to climate change, long-lived mammals must have adaptive mechanisms such as phenotypic or behavioral plasticity (within individual lifetimes; Przybylo et al. 2000, Lane et al. 2012, Hoy et al. 2018, as microevolutionary responses in animals with long generation times are suggested to occur at much slower rates than current climate change. Recent studies investigating temporal shifts in breeding timing of large herbivores as response to climate have found two reaction types: (1) one not being able to keep track of the changing environment because its reproductive cycle mainly depends on photoperiodicity, and (2) another that shows temporal flexibility because it relies in other cues such as temperature and food availability to keep synchrony with the environment. Most capital breeders belong to the second type as they have shown some degree of phenotypic plasticity, for example, by advancing conception dates (e.g., Pyrenean chamois Rupicapra pyrenaica pyrenaica: Kourkgy et al. 2016;red deer Cervus elaphus: Pel aez et al. 2017) or shortening gestation length (e.g., red deer: Moyes et al. 2011; reindeer Rangifer tarandus tarandus: Mysterud et al. 2009;semidomesticated reindeer: Paoli et al. 2018). In fact, recent studies have demonstrated that capital can be used, not only to lessen the need for matching optimal timing in parturition, but also to facilitate temporal synchrony between resource need and supply (Williams et al. 2017). On the contrary, income breeders exhibited limited potential to track early-spring plant growth because their reproduction relies strongly on photoperiodic cycles (e.g., West Greenland caribou Rangifer tarandus groenlandicus: Forchhammer 2008, Kerby andPost 2013b; roe deer Capreolus capreolus: Plard et al. 2014). Income breeding under stable conditions allows animals to synchronize their energy demands with supply (J€ onsson 1997) following a bet-hedging strategy in which parturition timing coincides with a long-term average of favorable conditions (Post and Forchhammer 2008). Thus, capital breeding seems to be beneficial under changing environmental conditions, whereas income breeding is superior under predictable environment conditions. However, in the context of climate change, the income breeding strategy may turn out to be maladaptive (Kerby and Post 2013b).
Roe deer stand out as a species particularly unable to modify plastically parturition dates, as it has been demonstrated in a densely forested lowland area of France (Plard et al. 2014). However, little is known about the response of roe deer in other habitats such as topographically heterogeneous landscapes composed by a mosaic of forests, croplands, and pastures . Topographically heterogeneous landscapes integrate a larger range of climatic conditions per area (B€ untgen et al. 2017(B€ untgen et al. , Rehnus et al. 2018) than flat landscapes (Bellard et al. 2012), and thus, females have access to a wide range of plant phenology phases at small spatial scales. Hence, heterogeneous landscapes may help individuals to minimize the mismatch between parturition timing and vegetation flush (Fryxell and Sinclair 1988) by following a multi-range tactic as a response to resource heterogeneity and unpredictability (Couriot et al. 2018). Finally, while populations in topographically homogeneous areas might be at risk of falling outside the optimal climate niche, populations living in heterogeneous areas might still find suitable environmental conditions on nearby higher elevations.
Here, we investigated the response of roe deer parturition timing to track changes in plant phenology during 45 yr  along a distinct elevational gradient (288-2366 m a.s.l.) in Switzerland. We hypothesize that (H1) the mismatch between roe deer parturition dates and peak resource availability will increase over the study period because of the lack of phenotypic plasticity of roe deer to track temporal changes in plant phenology and (H2) the mismatch between roe deer parturition dates and peak resource availability will be larger in low areas compared to higher elevations.

Study area
The study area comprises three biogeographic regions that cover the whole elevational and climatic range of roe deer in Switzerland (Christen et al. 2018). These adjacent and connected regions are, from north to south and from low elevation to high, the Swiss Plateau and the Pre-Alps (SP/PA; mean elevation 639 m a.s.l.; ranging from 240 to 2361 m a.s.l.), the Northern Alps (NA; mean elevation 1389 m; ranging from 370 to 4045 m a.s.l.), and the Eastern Central Swiss Alps (ECA; mean elevation of 2111 m a.s.l.; ranging from 560 to 3981 m a.s.l.; Fig. 1; Gonseth et al. 2001). Details of the study area have been published elsewhere (Pel aez et al. 2020 ).

Roe deer fawn data collection
Data originated from the project "Rehkitzmarkierung Schweiz" (i.e., roe deer marking in Switzerland), which has been running since 1971. As of 2015, 15380 roe deer fawns had been marked by hunters, gamekeepers, and researchers in the study area with varying marking efforts in number of marked fawns over 45 yr. Although the application of an ear mark is simple, only people experienced in handling wildlife were allowed to participate in this project. Fawns were found by repeatedly observing the mother during late pregnancy until the day of parturition. For each fawn marked, the following information was recorded: the date of marking, the estimated age in days, the geographic coordinates, and the herbaceous layer height at the marking site (<20, 20-50, >50 cm). Fawn age was assessed by examining the umbilical cord and observing fawn behavior during marking (Jullien et al. 1992). Parturition date was then calculated by subtracting the estimated age of the fawn (in days) from the date of capture. Marking site elevation, in meters, was derived from a digital elevation model (DEM) with a 25 m resolution (Zimmermann and Kienast 1999) using ArcGIS 10.6 (ESRI 2018). Finally, elevation was classified into 250-meter intervals.

Environmental variables
Although roe deer has been described as a forest-dwelling and browser species, open herbaceous patches such as meadows constitute an important diet source for roe deer during the breeding season and are actively selected by females (Tixier and Duncan 1996). Indeed, their recent colonization of modern agricultural landscapes throughout Europe has demonstrated a better performance in these habitats where they have access to high-quality food such as crops and wild forbs (Tixier and Duncan 1996, Hewison et al. 2009. Therefore, we characterized the optimal timing for roe deer births, as the period between the start of the growing season (GSS) and the start of the herbaceous flowering (FS). We selected two different variables to define this period: (1) the date after six consecutive days with mean daily temperatures above 5°C as a proxy for GSS (Zimmermann and Kienast 1999) and (2) the date of the first cutting of hay as a proxy for FS. Reasons to select the date of first cut of hay as proxy for FS were as follows: (1) Hay is cut when the ratio of forage quantity/quality is at its maximum, just before the start of flowering, (2) hay cut information matched the full temporal and elevational range of our roe deer data, and (3) it represents the phenological state of the plant community, not just of a single species. Similarly, previous studies on roe deer used the date of flowering as a proxy of spring plant phenology (Plard et al. 2014).
GSS was calculated from the temperature measurement network of MeteoSwiss (Federal Office of Meteorology and Climatology 2018). Dates of the first cutting of hay were obtained from the national phenology network of MeteoSwiss (www.meteoswiss.admin.ch/home/measureme nt-and-forecasting-systems/land-based-stations/ swiss-phenology-network.html). We selected the most representative stations for our study area (those inside the minimum convex polygon including 99% of roe fawns locations). In addition, we only used stations with at least 3 yr of data collected per decade throughout the study period. In total, we used data from 59 phenological stations and 98 meteorological stations located at elevations ranging from 330 to 1900 m a.s.l. Station elevations were also classified into elevational intervals of 250 m. The minimum number of years collected by any station included in our analysis was 24 yr.

Statistical analysis
Spatial scale of analysis.-We performed analysis at three different scales: (1) inter-regional scale, (2) regional scale, and (3) local scale. For the inter-regional-scale analysis, we pulled all data from the study area. For the regional analysis, we divided the three datasets (parturition date, GSS, and FS) by biogeographic region (SP/PA, NA, and ECA). Finally, at the local scale, we performed a spatial analysis to find small areas (maximum of 7.5 km radius) with the highest concentration of parturition recorded (at least 200 parturitions) and a maximum of 250 m a.s.l. difference between the highest and the lowest parturition recorded (see Data S1 for more details). We assumed that, within these small areas, there would be little genetic differentiation among individuals, as the 85 percentile of adult individuals in Switzerland dispersed below this distance (Euclidean distance, unpublished data). Finally, we selected the two closest meteorological or phenological stations for each small area. Parturition date and phenological information (GSS, FS) were expressed in Julian days (JD).

Temporal changes in plant phenology
To assess temporal changes in plant phenology at inter-regional, regional, and local scales, we performed linear models with GSS and FS as response variables. Predictors were year, elevational interval, and the interaction between the two. In addition, to provide an overall estimate of the rate of advance of GSS and FS during the study period, we calculated the average rate of advance across elevational intervals for each spatial scale and then calculated the average rate among the three spatial scales.

Temporal changes in parturition date
We used parturition dates instead of date of birth of individual fawns to avoid pseudo-replication problems (sensu Hurlbert 1984) caused by the presence of twins or triplets (i.e., two or three fawns from the same mother in the same year). To ensure that parturition dates were calculated accurately, we only used information from fawns aged ten days or younger because aging precision decreases with increasing fawn age (Jullien et al. 1992; mean fawn age 4.69 AE 2.79 d). We only included marking sites with a data resolution ≤1 ha. Ultimately, we used information from 8986 parturition events recorded during the period 1971-2015. Finally, for the three scales of analysis, we calculated the median parturition date by year and elevational interval.
To assess temporal changes in parturition dates, we fitted linear regressions with a Gaussian error distribution using median parturition date as the response variable. As predictors, we included year, elevational interval, and their interaction. In addition, we included the number of parturitions used to calculate the median as a weighting variable to adjust the relative influence of each year based on the number of parturition dates recorded. Finally, we calculated the average rate of advance in parturition date across elevational intervals for each spatial scale and then calculated the average rate among the three spatial scales.

Mismatch quantification
Results of the models were used to compare differences in the temporal trend between peak parturition and peak resource availability (GSS and FS) as an indicator of mismatch.
As additional indicator, we analyzed the temporal variation of the herbaceous layer height at the moment of fawn marking. For this analysis, we only used information on individual fawns found in meadows, as the development of the ground vegetation layer is related to tree canopy closure and light conditions in forests (N = 6219 fawns; 71% of the total fawns found with known vegetation height). Therefore, we created a binomial response variable by giving value one to the marking sites with low herbaceous layer height (<20 cm) and value zero to the marking sites with medium (20-50 cm) or high height (>50 cm). We then fitted a generalized linear model with a binomial link function including elevational interval, year, and their interaction as predictors. Because this corresponded to a smaller subset of our data, we only analyzed it at the inter-regional and regional scales.
All analyses were performed using R 3.2.1 statistical software (R Development Core Team 2018). We used the function lm to fit the models and the plot function to check model assumptions through diagnostic plots (Pinheiro and Bates 2000). Both functions were part of the stats package.

Temporal changes in plant phenology
Results showed that GSS and FS advanced significantly through time independently from the scale studied ( Fig. 2; Appendix S1: Fig. S1). Furthermore, the rates of advance of GSS and FS were very similar across different elevation intervals and scales, advancing on average at a rate of 0.45 and 0.32 d per year, respectively (Table 1). Thus, the date of the GSS and FS advanced on average 20 and 14 d, respectively, over the 45 yr of study.

Temporal changes in parturition date
A slight trend toward more advanced parturition dates through time was found at different scales and elevational intervals (Table 1; Fig. 2; Appendix S1: Figs. S1, S2). On average, the advance in parturition date was 0.06 d per year (ranging between 0.05 and 0.06 depending on the scale studied; Table 1). However, note that significance levels varied among elevation intervals and scales, but the direction of this trend ❖ www.esajournals.org signaled almost in all cases toward an advance in parturition date (Table 1; Appendix S1: Fig. S2). Across the three different scales studied, we calculated a total of 25 independent slopes for the relationship between median parturition date and time; for 23 of them, this relationship was negative ( Table 1). The results obtained by using mean parturition dates instead of the median were very similar (Appendix S1: Table S1). Therefore, parturition date is now about 3 d earlier than 45 yr ago.

Mismatch quantification
The magnitude of advance in parturition dates through time was much smaller than that of plant phenology (GSS and FS) at all elevations and scales. Parturition dates changed at a rate 7.5 times slower than GSS and five times slower than FS. Thus, due to the slower rate of advance of parturition date in relation to plant phenology, we observed that, below 1000 m a.s.l., mean parturition timing has already fallen outside the optimal time window for high-quality forage ( Fig. 2; Appendix S1: Fig. S1). In contrast, parturitions occurring at elevations above 1000 m a. s. l. still occurred within the optimal period for high-quality forage ( Fig. 2; Appendix S1: Fig. S1).
The faster rate of advance of plant phenology compared to parturition dates was also supported by the measurements of vegetation height at marking sites. Throughout the study period, the percentage of fawns found in low-height meadows (<20 cm) decreased from 61% to 20% at 1750 m a.s.l., from 33% to 10% at 1250 m a.s.l., and from 13% to 5% at 750 m a.s.l (inter-regional scale; Fig. 3). We found no significant changes at the lowest elevation interval because the percentage of fawns born at sites with low vegetation was already <7% at the beginning of the study. Overall, the same trend was obtained at the regional scale (Appendix S1: Fig. S3).

DISCUSSION
Here, we found a consistent advance toward earlier parturition dates across all scales and most elevational ranges that could be associated with a response to climate change. However, this advance was much slower than that of plant phenology, creating a trophic mismatch between parturition dates and peak resource availability over the 45 yr of study . Furthermore, at low elevations, the phenological decoupling has already caused mean parturition dates to fall outside the temporal window of optimal forage quality, while at high elevations, parturitions are still synchronized with peak resource availability.
Clear evidence was found for a strong climate change signature in plant phenology, resulting in an earlier spring plant growth in the study area Table 1. Parameter estimates and associated standard error (SE) from the models performed to assess temporal variation (days per year) of median parturition date, the start of the growing season (GSS), and the start of the herbaceous flowering (FS) for the selected elevational intervals at the inter-regional and local scales for the period 1971-2015. *** P < 0.001; **P < 0.01; *P < 0.05; †P < 0.1; ns, P > 0.1. ‡ The ellipsis indicates not enough data to perform lm.
(À0.45 and À0.32 d per year for start of the growing season and the start of the herbaceous flowering, respectively). This result is in line with previous studies performed across Switzerland (Klein et al. 2016) and Europe (e.g., Menzel et al. 2006). In addition, vegetation height measurements at the marking sites revealed that the percentage of fawns found in low-height meadows progressively decreased throughout the study period. This indicates that roe deer fawns were born each year at a later stage of herbaceous development and therefore further away from the leading edge of the growing season and closer to the flowering start. Despite the increasing mismatch observed between plant phenology and parturition date, we did find a consistent advance toward earlier parturition dates across all scales and most elevational ranges (although with varying levels of significance; see Yoccoz 1991). Nevertheless, plant phenology advance was between five and seven times larger than the rate observed for roe deer parturition dates, which was on average À0.06 d per year. Similarly, in Trois Fontaines, a densely forested lowland area of France, the different rate of advance between plant phenology (À0.60 d per year) and roe deer parturition dates (mean birth date À 0.05 AE 0.06 d per year; P = 0.421 and median À 0.087 AE 0.070, P = 0.232) also caused a phenological mismatch (Plard et al. 2014). Interestingly, the minor advance in roe deer parturition timing observed in our study area was of the same magnitude and direction of that obtained by Plard et al. (2014), although for this study the trend was not significant. Our results suggest that, although very slow and at population level, roe deer parturition timing has changed in the direction of advancing peak resource availability. This is an important biological finding for the species which, at least in heterogeneous environments, is showing some degree of response to climate change. This highlights the importance of performing large-scale studies across wide environmental gradients to detect subtle phenological trends on key life-history events that may otherwise remain undetected.
Unfortunately, in comparison with other Northern Hemisphere large herbivores studied, roe deer is the species that show the smallest response to climate change (À0.06 d per year). West Greenland caribou, another income breeder, advanced À0.11 (AE 0.03 SE) days per year over 33 yr (Kerby and Post 2013a). Capital breeders like red deer have advanced parturition by À0.42 (AE 0.08 SE) days per year over 26 yr in Scotland (Moyes et al. 2011). In Finland, semidomesticated reindeer, which are also described as capital breeders, advanced parturitions À0.15 (AE 0.04 SE) days per year over 45 yr (Paoli et al. 2018). For Pyrenean chamois in France, a 10-d advance in the onset of autumn or spring plant phenology led to an advance of four and one day in birth dates, respectively (Kourkgy et al. 2016). For all the mentioned studies, the observed changes in parturition date have been attributed mainly to phenotypic plasticity (Moyes et al. 2011, Kerby and Post 2013a, Kourkgy et al. 2016, Paoli et al. 2018. Only recently, a study demonstrated that adaptive evolution likely explained some of the advance in red deer parturition dates during the last four decades on the Isle of Rum, Scotland (Bonnet et al. 2019), suggesting that up to a total of À0.05 d of delay per year (95% CI À0.10 to 0.02) could be explained by evolution.
Unfortunately, we lack information to infer the mechanism responsible for the small advance in parturition date observed for roe deer. An overall advanced parturition date has been reported for older females (Plard et al. 2014) or high-quality females (note that quality here is defined as a set of time-invariant attributes positively linked to individual fitness; ). However, we have no information on population demography or on female age, body weight, and longevity in our study area. Therefore, we could not assess whether (1) female quality at population level varied through time as a response of changes in population density or land use (i.e., increased resource availability) or (2) an increase in the mean age of female at population level was responsible for the observed trend. However, the large sample size used in this study should render the effect of changes in age structures or population demography. Thus, due to the consistency of the results across scales and elevations, we do not expect that our results are biased by these variables. Alternatively, microevolutionary processes may be the cause of the slight advance observed in parturition date. However, a previous study found no statistical support for heritability of birth date for female roe deer based on 28 daughter-mother pairs, even though birth date was highly variable within a population and a trait under strong directional selection (Plard et al. 2014). Therefore, the study predicted that a microevolutionary response of parturition date to climate change would be limited for roe deer. Nevertheless, this question is still under discussion as a much larger sample size is required to confirm these results (Plard et al. 2014). Finally, an alternative explanation would be that the observed advance corresponded to migration or dispersion movements toward higher elevations. This would mean that, instead of rescheduling breeding timing in situ to counteract the phenological mismatch via phenotypic plasticity or microevolution, roe deer tracked appropriate conditions in space by moving toward higher elevations (i.e., behavioral plasticity). Such an adaptive behavior is likely to occur in topographically heterogeneous ecosystems due to their smallscale spatial heterogeneity (B€ untgen et al. 2017). Indeed, in mountainous areas, roe deer population exhibited partial migration, with a percentage of the population being migratory and the rest being stationary (Mysterud 1999, Cagnacci et al. 2011. For this species, migration and dispersion movements may constitute a good strategy to match parturition dates and peak resource availability during the breeding season. However, we were unable to discriminate between both options, microevolution or behavioral plasticity. Therefore, more research is needed to understand the mechanisms that led to the small response observed along the elevational gradient. Only then, better population-specific predictions on the rate of response of roe deer birth timing to climate change would be possible.

ACKNOWLEDGMENTS
We thank all gamekeepers, hunters, and researchers who participated in the marking of fawns over the study period. We also acknowledge the Federal Office for the Environment for the permission to use data from the project "Rehkitzmarkierung Schweiz." Marta Pel aez acknowledges financial support from the Spanish Ministry of Education (FPU13/00567 and EST16/ 00095). We acknowledge J.-M. Gaillard for his advice on an earlier version of the manuscript and E. Gleeson for language editing. Maik Rehnus and Marta Pel aez contributed equally to this work LITERATURE CITED