The onset in spring and the end in autumn of the thermal and vegetative growing season affect calving time and reproductive success in reindeer.

A developing trophic mismatch between the peak of energy demands by reproducing animals and the peak of forage availability has caused many species' reproductive success to decrease. The match-mismatch hypothesis (MMH) is an appealing concept that can be used to assess such fitness consequences. However, concerns have been raised on applying the MMH on capital breeders such as reindeer because the reliance on maternal capita rather than dietary income may mitigate negative effects of changing phenologies. Using a long-term dataset of reindeer calving dates recorded since 1970 in a semidomesticated reindeer population in Finnish Lapland and proxies of plant phenology; we tested the main hypothesis that the time lag between calving date and the plant phenology in autumn when females store nutrient reserves to finance reproduction would lead to consequences on reproductive success, as the time lag with spring conditions would. As predicted, the reproductive success of females of the Kutuharju reindeer population was affected by both the onset of spring green-up and vegetative senescence in autumn as calves were born heavier and with a higher first-summer survival when the onset of the vegetation growth was earlier and the end of the thermal growing season the previous year was earlier as well. Our results demonstrated that longer plant growing seasons might be detrimental to reindeer's reproductive success if a later end is accompanied by a reduced abundance of mushrooms.

The ongoing large changes in climate have already caused unprecedented consequences on the phenology of many plant and animal species. In the Arctic tundra ecosystem, a warming 2–3 times higher than the increase of the global mean surface temperature (IPCC 2013) has led to longer vegetative growing seasons (Oberbauer et al. 2013); and a 50% increase in the above ground vegetation biomass (Epstein et al. 2000). In response to these plant phenological changes, many animal species have shifted their timing of reproduction (bird: Visser et al. 1998; Frederiksen et al. 2004; Bourret et al. 2015; amphibian: Blaustein et al. 2001; fish: Asch 2015; mammal: Réale et al. 2003; Moyes et al. 2011; Lane et al. 2012; marine species: review by Poloczanska et al. 2013). However, there is increasing evidence of a mismatch, due to these phenological changes, between the peak of resource demands by reproducing animals and the peak of forage availability (Post and Forchhammer 2008). Although in ruminants, and especially ungulates, the peak of available high-quality forage is of better importance than simply the peak of forage biomass as they are more selective on forage quality than forage quantity (Tveraa et al. 2013). Also, the question on whether potential shifts in timing of reproduction are sufficient to track a rapidly changing environment is still a debated question with contrasted answers. For instance, climatic changes were shown to affect the reproduction of ungulate species either negatively or positively. Negatively because deep snow cover and “locked pastures” under an impenetrable layer of ice restrict access to the field layer and to forage during winter for herbivores like Rangifer species (Aanes et al. 2000; Tyler 2010); thus generating reduced reproductive rates down to 49% (Helle and Kojola 2008), major die-offs and population declines (Tveraa et al. 2007). Positively because the timing of flowering was advanced in parallel with earlier snow melting. Also, deep snow may lead to an extended period of access to newly emergent high-quality forage (Post and Stenseth 1999), which was shown to be favorable to both red deer Cervus elaphus and Soay sheep Ovis aries L. (Mysterud et al. 2003), to caribou (Forchhammer et al. 2002) and more recently to reindeer (Tveraa et al. 2013). Those contrasted findings highlight the need for a better understanding of how quickly species are able to adapt to their changing environment and if long-term fitness consequences can be detected.

The match–mismatch hypothesis (MMH; Cushing 1990) was developed to address these issues. The MMH is used to estimate the gap between the phenology of a species at the higher level (e.g., grazers) and that of species at the immediate lower level (e.g., vegetation). If there is a time lag between the food requirements and the food availability of the species at the higher level, that is, a mismatch, then the survival and reproduction rate of this species are expected to be low (Durant et al. 2007). For example, the gosling body mass and structural size at fledging of a snow goose Chen caerulescens atlantica population was reduced when the trophic mismatch between hatch dates of young and date of peak in plant nitrogen availability was high (>9 days; Doiron et al. 2015). The MMH is an appealing concept regarding the ongoing climate change as we could better describe the changes of the time lag between the peak of resource demands by reproducing animals and the peak of forage availability (Post and Forchhammer 2008). But to date, the MMH has been rarely applied on ungulates and with contradictory results. For instance, the temporal mismatch found between the vegetation peak and the average birth date of the Soay sheep in Hirta, Outer Hebrides in Scotland did not affect the lambs’ survival (Durant et al. 2005), whereas in the Low Arctic Greenland, the increased time lag between the onset of the plant growing season and the caribou’s timing of calving has contributed to a reduced production and survival of caribou calves (Post and Forchhammer 2008).

Criticisms have been raised on applying the MMH on capital breeders such as caribou and reindeer (Taillon et al. 2013) because the reliance on maternal body reserves rather than dietary income to ensure fetal development could act as a buffer against negative effects of between-year variability in plant phenology (Veiberg et al. 2016; Gustine et al. 2017), curtailing potential effects of spring phenological mismatches. As such, Veiberg et al. (2016) recently found that the influence of spring phenology was less important than that of maternal winter body condition on reproductive success of the Svalbard reindeer Rangifer tarandus platyrhynchus. As previously documented, the peak of calving time for many Rangifer populations usually occurs well (up to 3 weeks) before the spring green-up (Tveraa et al. 2013; Veiberg et al. 2016; Gustine et al. 2017) in order for the timing of high-energetic demands (i.e., lactation) to coincide with the timing of peak nutrient availability/quality in forage species (Klein 1990; Post et al. 2003). An increase in the preexisting time lag between parturition time and the emergence of vegetation would then diminish both mother and calf’s ability to exploit high-quality forage during a period of high nutritional requirements (i.e., lactation, replenishing winter fat reserves, and calf physical growth; Post and Forchhammer 2008), reinforcing the importance of maternal capital for reproduction (Barboza and Parker 2008). Gustine et al. (2017) therefore suggested that it might be more relevant to consider potential lagged effects in the periods when females build up reserves to ensure reproduction in the next breeding season (e.g., in the late summer and autumn; Thompson and Barboza 2014) rather than at parturition. Indeed, the major determinant for calving date was observed to be the conception date (Holand et al. 2002; Clements et al. 2011), therefore, inadequate nutrition in autumn would be detrimental to females’ body condition and cause a subsequent delay in calving dates (Cook et al. 2004).

To circumvent these issues and following the suggestion made by Gustine et al. (2017), we therefore examined the potential climate-induced lagged effects at the individual-level between calving date and 1) start of vegetation emergence in spring and 2) end of the vegetative growing season in autumn, which is new. We also determined if the potentially lagged effects of plant phenology led to fitness consequences, as recently found in other reindeer populations (Tveraa et al. 2013; Veiberg et al. 2016). Using a long-term dataset of reindeer calving dates recorded since 1970 in a population in Finnish Lapland and proxies of plant phenology, we specifically tested the following 3 hypotheses: 1) The reported rate of advancement in calving dates (Paoli et al. 2018) would match the rate of change in spring green-up in the study area, with no changes in time lag over time as already reported in another caribou population (Gustine et al. 2017); 2) The calving time would be more affected by the plant phenology in autumn than in spring, mediated by effects on females’ nutrition and body condition at the time of conception (Cook et al. 2004; Veiberg et al. 2016); 3) The time lag between the end of the vegetative season in autumn and the forthcoming calving time may be expected to diminish over time, because of a reported delay in autumn leaf coloring by up to 2.4 days°C−1 in Europe from 1971 to 2000 (Menzel et al. 2006).

Materials and Methods

Study area and reindeer population

Calving dates were obtained from an experimental herd of semidomesticated reindeer of about 100 animals in every year at the Kutuharju field reindeer research station in Kaamanen (kept by the Reindeer Herders’ Association), northern Finland (69°N, 27°E). The study area is composed mainly by open birch Betula spp. and pine Pinus sylvestris forests, bogs and lakes and the landscape varies between 185 and 370 m above the sea level. The long-term book-keeping of the herd demography and the unique identification of animals fitted with collars and ear tags since their birth allowed reindeer to be individually recognizable and of known age. Most of the year (in summer and during the rut) the herd is free ranging in 2 large fenced enclosures, the north-west section (Lauluvaara ∼13.8 km2) and the south-east section (Sinioivi ∼15 km2). When the mating season is over around late October, the animals are gathered and taken to a winter grazing area (15 km2) where they can graze freely on natural pastures. In late winter, the animals receive supplemental feed (pellets and hay) in addition to natural forage as part of a common management practice that has started since the 1980s in the northern part of Finland (Helle and Jaakkola 2008). The significant between-years variability in females’ body weight (1-way analysis of variance: F(44,2,094) = 21.3, P < 0.001) makes us believe that the regular supplemental feeding alone was not able to sustain individuals’ body weight at a stable level and therefore, females were still dependent on natural pastures. No detailed information was available on the duration or the amount of supplemental feeding previously given every year to the animals. By the end of winter, females are transferred into a calving enclosure (∼0.5 km2) where newborn calves are captured, weighed, sexed, and marked with ear tags. The enclosure was surveyed daily, so that calving date is known for all individuals and has been recorded since 1970. Handling of live mammals and data collection were done in agreement with the Animal Ethics and Care certificate provided by Concordia University (Protocol number 30000303) and by the Finnish National Advisory Board on Research Ethics.

Calving dates and population variables

The calving dates from females that have been artificially fed for experimental purpose, anytime during the year preceding the calving season (i.e., more than just the regular extra feeding given in late April) were excluded from the dataset (see Paoli et al. 2018). Given the reported trend of heavier (and older) females giving birth earlier (Cameron et al. 1993; Flydal and Reimers 2002), an artificial feeding by improving the females’ physical condition could lead to earlier calving dates, independently from plant phenology. This artificial feeding could therefore buffer the effects of plant phenology on calving date. All calendar dates were converted into Julian days since 1 January for analysis purposes. In total, calving dates spanning over 45 years were available from 1970 to 2015. A lowercase “t” indicated that the variable was estimated from the same calendar year as calving, whereas a lowercase “t−1” indicated that the variable was estimated from the previous calendar year as calving. The documented effects of females’ physical condition, population sex ratio, and population density on calving date (Ropstad 2000; Forchhammer et al. 2001; Flydal and Reimers 2002; Holand et al. 2002; Mysterud et al. 2009) were also controlled for.

The proportion of males during the previous mating season (“PM−1”) was estimated as the number of males divided by the number of females over 1 year of age present in a specific enclosure. Between 1996 and 2013 (except 1998), the proportion of males was estimated per enclosure as the herd was subjected to a number of experiments including manipulation of the proportion of males, leading to the simultaneous use of the 2 large enclosures, Sinioivi and Lauluvaara. Every year, all animals are gathered in corrals just before the rut period (in September) and then released in either Sinioivi or Lauluvaara for the entire rutting season. As the identity of the animals involved in any experiment was known, as well as their presence in each enclosure (recorded during the release), every calving date was related to the corresponding proportion of males in that enclosure the previous mating season. In addition to the proportion of males, we also estimated the population density (“DENS−1”) per enclosure-year as the total number of individuals present in a specific enclosure during a given year.

All factors linked to physical condition in reindeer interact with each other so that older individuals tend to be heavier (Ropstad 2000). Therefore, an individual body condition index (BCI) was used to allow 1) the effects of individual body weight controlling for age to be taken into account and 2) multicollinearity between these 2 highly correlated variables to be avoided. This BCI was estimated by a measure of individual body weight after the effect of age is controlled—the age-specific residual body mass (see Festa-Bianchet et al. 1997; Weladji et al. 2003), calculated by subtracting from each individual’s body weight the average body weight of all individuals of the same age, across all years. To also account for the reported senescence in female reindeer from this population (Weladji et al. 2010), we extracted the residuals from the linear and quadratic forms of the relationship between females’ body weight and females’ age. The females’ BCI was averaged across 3 key periods in their reproductive cycle: in autumn (BCI Autumn−1, September–October–November) as a female’s conception date (and therefore calving date) is influenced by her physical condition at the time of mating (Clutton-Brock et al. 1982; Adams and Dale 1998; Flydal and Reimers 2002); in winter (BCI Winter, from December [the previous calendar year] to March [the same calendar year as calving date]) as winter climatic conditions during pregnancy can cause females’ nutritional stress and are critical for the fetal growth and development (Barboza and Parker 2008); and in early spring (BCI Spring, April–May) as climatic conditions in April–May were already found to affect calving date (Paoli et al. 2018) and that most of the energy requirements for gestating females occur during the third trimester of pregnancy (Parker et al. 2009).

Plant phenology data

The starting and ending day of the thermal growth season in each year were estimated using daily average temperatures from 1970 to 2015 obtained from 3 weather stations (Utsjoki, Ivalo airport and Nellim) in northern Finland (68°N, 27°E) downloaded with the “FMIDownloader” of the Finnish Meteorological Institute (https://tuomassalmi.com/FMI-weather-downloader/). To obtain the daily average temperatures at our study site with as much reliability as possible, the daily recorded temperature values from the 3 nearest weather stations were averaged with their weighted mean by the distance from the weather station to our study site. The distance between our study site and each of the weather stations was precisely assessed using their respective GPS coordinates and the Great Circle longitude–latitude calculations tool (http://www.cpearson.com/excel/LatLong.aspx). Then, we defined the onset of the thermal growing season in each year as the last day of the first 5-day period, after the last winter/spring frost, during when the daily mean temperatures all remained above 5°C (Jones et al. 2002; Walther and Linderholm 2006). The last day of frost was defined as the last day in late winter/spring when the minimum temperature was negative (<0°C, Menzel et al. 2003). The end of the thermal growing season was computed as the first day of the first 10-day period when the mean daily temperatures fell <5°C (Carter 1998; Walther and Linderholm 2006). The length of the thermal growing season (LTGS) was assessed by the day of end (“ThermalEnd”) minus the day of start (“ThermalStart”) of the thermal growth season.

Vegetation phenology in the study area was assessed using the phenology of deciduous birch Betula pubescens as this species is commonly present in the lowland regions of the study area and is the dominant tree species in most parts of the study area (Karlsen et al. 2008). The deciduous growth of birch also allows to well-defined phenophases, phenomena which are not observed easily in conifers, making birch species an ideal phenological indicator (Karlsen et al. 2008). The vegetation phenology data were obtained from the Muddusjärvi station situated at ∼19.3 km from the study area and monitored by the Natural Resources Institute of Finland (Luke), Finland. Data are part of the Pan European Phenology Project PEP725 and were downloaded through their website (http://www.pep725.eu/index.php) and were available from 1997 to 2015. In the present study, 2 phenophases were used as indices to reflect, respectively, the onset in spring and end in autumn of the vegetative growing season: “onset for sprout of birch leaves” (“VegOnset”) and “50% autumnal coloring of birch leaves” (“VegEnd”). Onset of leafing is the date when the first leaves sprout and a “mouse ear” is visible (BBCH10 code according to Meier 1997, p. 622). The second phenophase “50% autumnal coloring of birch leaves” is the date when half of the leaves on each trees have turned yellow (BBCH94 code). The chosen birch phenophases represent well the general greening and coloring of the region’s vegetation (Karlsen et al. 2008). The length of the vegetative growing season (“LVGS”) was assessed by the day of VegEnd minus the day of VegOnset of birch. We conducted analyses on 2 study periods depending on the availabilities of the plant phenology variables: the whole study period (1970–2015) and the last part of the study period (1997–2015). The correlations between the plant phenology variables were assessed using Pearson’s correlation test.

Statistical analyses

Temporal trends

To evaluate the long-term trends of the population variables (population density and proportion of males) and plant phenology variables in the Kutuharju area, we fitted Linear Models (LMs) with the year of study entered as a covariate (continuous variable) in the models. Variation in calving date, our response variable, was analyzed using a Linear Mixed-effects Model (LMM), by running the lmer-function in the R package lme4 (Bates et al. 2015, ). Year only was entered as a continuous fixed-effect parameter in the model, and year as a multilevel random effect to control for repeated measures (Kruuk et al. 1999; Milner et al. 1999). Unstandardized value of the temporal trend was reported and the parameter estimate was derived using the restricted maximum likelihood estimates (REML) as recommended for mixed-effect models (Bolker et al. 2009). The same method was applied for the 3 metrics of females’ BCI as calving date and females’ BCI had several values each year. For all models, the temporal trends were considered statistically significant if 95% confidence intervals (CIs) of the parameter estimates excluded 0.

Plant phenology variables

We then evaluated whether calving date was affected by the start, end, and duration of both the thermal and vegetative growing seasons by running LMMs with the population density (DENS−1) and proportion of males (PM−1) entered as fixed-effects in a “basic model” to account for their known effects on calving date, with calving date as the response variable. In total, 9 predictor variables were tested in the models: ThermalStart−1, ThermalEnd−1, LTGS−1, VegOnset−1, VegEnd−1, LVGS−1, BCI Autumn−1, BCI Winter, and BCI Spring. Given that the plant phenology variables had only 1 value per year and that we applied the “detrending” method (see below), only the females’ identity was entered as a multilevel random effect in the models to account for within-individual dependency (Pinheiro and Bates 2000). If statistically significant temporal trends were found for any of the explanatory variables, we applied the “detrending” method to avoid the potential confounding effect of the year on the phenology-calving date relationship, which can happen simply because both variables change across years (Iler et al. 2017). To do so, we extracted the residuals from independent linear regressions with year as the predictor from the population, phenological, or calving date variables (Iler et al. 2017). All the predictor variables considered were also centered and standardized (X– = 0, standard deviation [SD] = 1) to be on a comparable scale. We also tested the interactions between the metrics of females’ BCI (BCI Autumn−1, BCI Winter, and BCI Spring) and plant phenology variables (ThermalStart−1, ThermalEnd−1, VegOnset−1, VegEnd−1). Since our predictor variables were beforehand centered and standardized, we could directly interpret their main effects even when involved in interactions and thus avoided the potential misinterpretation of main effects between models with and without the interaction term (Gelman 2008; Schielzeth 2010; Grueber et al. 2011).

A model selection was then performed to find combinations from all the explanatory variables used providing the most probable models to explain variation in calving date and this was based both on the Akaike Information Criterion (AIC) and Akaike weights (AIC weights) to compare the relative performance of the models tested (Anderson et al. 2001; Burnham and Anderson 2002). The delta AICc () was calculated to provide a measure of each model relative to the best model (with the lowest AIC value). All models within a ΔAIC of 2 units were retained as competing models because a substantial evidence was given to the model if (Burnham and Anderson 2002). To account for model selection uncertainty and if >1 model were retained as best models in explaining the data then the estimates of the coefficients of parameters in all models with ΔAIC < 2 were averaged, following the model averaging approach (Schielzeth 2010; Grueber et al. 2011; Mazerolle 2017). The model.avg function in the R package AICcmodavg (Mazerolle 2017, www.r-project.org) was used to calculate the model-averaged parameter estimates, as well as their 95% CIs based on our entire list of candidate models. The model selection of mixed models was performed using the maximum likelihood methods, whereas the parameter estimates were obtained from the selected model using the REML as recommended for mixed-effects models (Bolker et al. 2009; Zuur et al. 2010). The effect size of the predictor variables included in the competing models was considered important if the 95% CIs excluded 0 and only the important variables were further discussed. The analyses were repeated over the 2 study periods, with a total of 30 models performed for the whole study period and 84 for the last part of the study period.

Time lags

To test our hypotheses of changes in time lags between calving date and spring or autumn phenology, we first estimated such time lags in number of days as follow: , , , . Since the onset of thermal or vegetative growing season occurs after calving in the same calendar year, we simply calculated the difference between ThermalStart or VegOnset and calvdate for TLThermalStart and TLVegOnset (and given that the plant phenology variables are expressed in Julian days). Given that ThermalEnd−1 and VegEnd−1 were estimated from the previous year to calving, we first calculated the number of days remaining from ThermalEnd−1 or VegEnd−1, respectively, until the end of the year (; 31 December being the last day of the year, it is equivalent to the 365th Julian day), and added this number to the number of days until calvdate (which is equivalent to the calving date as expressed in Julian days since 1 January; ), to estimate, respectively, TLThermalEnd and TLVegEnd. The temporal trends in the time lags were then evaluated by running LMMs with year entered as a continuous fixed-effect parameter in the models, and year as a random effect on the intercept as several values were obtained in each year.

Fitness consequences

We also evaluated whether the plant phenology variables or the measures of time lags could have fitness consequences on calves. To do so, two females’ reproductive success attributes were used, such as the birth weight of her calves (“BirthWeight”) and the first-summer survival of her calves (“Surv”). Given that early-born calves are also heavier (Eloranta and Nieminen 1986), we have corrected the calves’ birth weight by their respective calving date, by extracting the residuals of the LM of calves’ birth weight fitted against calves’ calving date. The causes of death of a calf excluded from the analyses were those with “no information” or “slaughtered”. The values for the survival of a calf ranged from 0: dead during calving season to 1: survived to autumn (at the time of the round-up in September). Then, we performed LMMs for BirthWeight, with females’ identity only entered as a random factor on the intercept and plant phenology variables added 1 by 1 in the model and their influence assessed by looking at their estimate and the associated 95% CI. The same protocol was used with Surv but by fitting generalized linear mixed-effects models (GLMMs), with a binomial error structure and a logit link. Analyses were performed in R 3.6.0 (R Development Core Team 2019).

Results

Temporal trends and correlation tests

During the whole study period (1970–2015), the end of the thermal growing season (ThermalEnd) remained unchanged, whereas the start of the thermal growth (ThermalStart) occurred earlier over time (Figure 1A), resulting in a significant lengthening of the thermal growing season in the study area (Table 1). During the last part of the study period (1997–2015), neither the onset of the vegetative growing season (VegOnset), nor the end in autumn (VegEnd) changed over time and as a consequence, the LVGS did not vary over these years (Table 1). During the same study period (1997–2015), however, the start of the thermal growing season has advanced by 15 days, whereas the end date did not vary (Table 1).

Figure 1.

Lagged effects of plant phenology on calving date of the Kutuharju reindeer population in Kaamanen, northern Finland. For the whole study period (1970–2015), earlier calving dates were observed when (A) the start (ThermalStart−1) and (B) end (ThermalEnd−1) of the previous thermal growing season were earlier, and when (C) the females’ physical condition in winter (December–March) was higher (BCIWinter). The females’ physical condition in winter was deteriorated when ThermalEnd−1 was later (D). For the last part of the study period (1997–2015), earlier calving dates were observed when (E) the onset (VegOnset−1) and (F) end (VegEnd−1) of the previous vegetative growing season were earlier, and when (G) the females’ physical condition in spring (April–May) was higher (BCISpring). A decrease in the females’ physical condition in spring was observed when VegOnset−1 occurred later (H). All the dates are expressed in Julian day starting from 1 January (JD). The 95% CI band around the fitted line is also presented. Each point represents the averaged value of the response variable per value of the predictor variable for graphical clarity.

Table 1.

Temporal trends of plant phenology variables, mean calving date (calvdate) and population variables from the reindeer herd of the Kutuharju field reindeer research station in Kaamanen, northern Finland

Variables	Estimate	SE	95% CI	Total change over the study period
1970–2015
Thermal growing season
ThermalStartt	−0.24	0.08	−0.41, −0.08	–11 days
ThermalEndt	0.09	0.13	−0.17, 0.35	+4.1 days
LTGSt	0.34	0.14	0.04, 0.63	+15 days
Calving season
calvdatet	−0.14	0.04	−0.23, −0.06	–6.5 days
Population variables
BCIAutumnt−1	0.29	0.04	0.21, 0.38	From −8.4 to 4.8
BCIWintert	0.32	0.04	0.24, 0.39	From −8.8 to 5.6
BCISpringt	0.31	0.07	0.16, 0.45	From −9.1 to 4.8
PMt−1	0.0002	0.001	−0.001, 0.002	From 16% to 17%
DENSt−1	−0.04	0.19	−0.42, 0.34	–1.8 individuals
1997–2015
Thermal growing season
ThermalStartt	−0.84	0.28	−1.44, −0.25	–15 days
ThermalEndt	0.04	0.43	−0.86, 0.95	+0.8 days
LTGSt	0.89	0.47	−0.10, 1.87	+16 days
Vegetative growing season
VegOnsett	−0.89	0.44	−1.83, 0.06	–16 days
VegEndt	−0.41	0.34	−1.13, 0.30	–7.5 days
LVGSt	0.007	0.76	−1.65, 1.66	+0.12 days
Calving season
calvdatet	−0.51	0.13	−0.77, −0.25	–9.2 days
Population variables
BCIAutumnt−1	0.36	0.08	0.21, 0.52	From −1.8 to 4.8
BCIWintert	0.38	0.06	0.26, 0.50	From −1.3 to 5.5
BCISpringt	0.40	0.13	0.13, 0.66	From −2.4 to 4.8
PMt−1	0.006	0.003	−0.0007, 0.01	From 9.5% to 21%
DENSt−1	−1.39	0.30	−2.02, −0.76	–25 individuals

The parameter estimates (with SE) from LMs with year as a fixed covariate indicate the direction of the change over time of the variables over the whole study period (1970–2015) or just the last part of the study period (1997–2015). ThermalStart and ThermalEnd were the start and the end, respectively, of the thermal growing season. VegOnset and VegEnd represented, respectively, the onset and the end of the vegetative season, estimated by birch phenology. BCIAutumn was the females’ BCI averaged over the autumn (September–October–November), BCIWinter the females’ BCI averaged in winter (from December to March) and BCISpring stood for the females’ BCI in spring (April–May, see text for details). The estimates in bold type indicated a statistically significant temporal trend (whose 95% CI excluded 0). DENS, population density; LVGS, length of the vegetative season; PM, proportion of males in the herd; SE, standard error.

The start of the thermal growing season was highly correlated to the start of the vegetative growing season estimated with birch onset of leafing (Table 2), with VegOnset occurring on average 2.5 days earlier than ThermalStart. The correlation between the end of the thermal growing season and end of the vegetative growing season in autumn was lower but still significant (Table 2), with VegEnd being followed by ThermalEnd by on average 21 days. The correlation between LTGS and LVGS was non-statistically significant (Table 2). The lack of statistical power by having too few years for the vegetative growing season might explain why no statistically significant temporal trends were detected for either VegOnset, VegEnd, or LVGS (Table 2). Following the reported temporal trends in Table 1, the calving date, BCI of females in autumn, winter, and spring, ThermalStart, and LTGS were detrended and then standardized for further analyses, whereas the other variables of the whole study period (proportion of males, population density, and ThermalEnd) were only standardized. The calving date, females’ BCI in autumn, winter, and spring, population density, and ThermalStart were detrended for the last study period and then standardized for further analyses, whereas ThermalEnd, LTGS, VegOnset, VegEnd, LVGS, and proportion of males were only standardized. Calving dates have occurred 6.5 days earlier (Figure 1B) during the whole study period and 9.2 days earlier during the last part of the study period.

Table 2.

Correlation coefficients between different plant phenology variables used in the study to understand the effects of plant phenology on calving date of a semidomesticated reindeer population in Kaamanen, northern Finland from 1970 to 2015

Variables	r	df	95% CI
ThermalStartt – ThermalEndt	0.06	44	−0.24, 0.34
ThermalStartt – VegOnsett	0.68	13	0.26, 0.88
ThermalStartt – VegEndt	0.35	16	−0.14, 0.70
ThermalEndt – VegOnsett	0.08	13	−0.45, 0.57
ThermalEndt – VegEndt	0.55	15	0.10, 0.82
VegOnsett – VegEndt	−0.001	12	−0.53, 0.53
LTGSt – LVGSt	0.38	12	−0.19, 0.76

Pearson’s correlation tests were used to determine the correlation coefficients between the plant phenology variables. The estimates in bold type indicated a statistically significant correlation between the 2 variables (whose 95% CI excluded 0). ThermalStart, ThermalEnd, and LTGS estimated, respectively, the start, end, and LTGS in the study area. VegOnset, VegEnd, and LVGS depicted, respectively, the onset, end, and LVGS in the study area, estimated by birch phenology (see text for details).

Plant phenology variables and time lags

The mean calving date was on 19 May, with the earliest reported peak calving date having occurred on 10 May 2010 and the latest on 30 May 1973. The average calving occurred 16 days prior to ThermalStart (mean = 16.2 days, 95% CI [4.58, 28.5]) and 2 weeks prior to VegOnset (mean = 14.4 days, 95% CI [2.59, 24.1]). The ThermalEnd occurred almost 4 months and half after mean calving date the same calendar year (mean = 134 days, 95% CI [118, 152]), and the end of the vegetative season in autumn almost 4 months after mean calving date (mean = 115 days, 95% CI [106, 124]).

For the whole study period (1970–2015), a total of 1,644 observations were analyzed and two competing models were selected (Table 3). The model-averaged estimates of the 2 best models revealed that calving date occurred earlier when the proportion of males in the herd (PM−1) was higher and the population density (DENS−1) was lower (Table 4), as expected from the literature. Calving dates were also occurring earlier when the start of the thermal growing season the previous calendar year (ThermalStart−1), as well as the end of the previous thermal growing season (ThermalEnd−1) were earlier (Figure 1A, B) and when the BCI of females in winter (BCIWinter) was higher (Figure 1C, Table 4). The interaction term between ThermalEnd−1 and BCIWinter was also statistically significant, indicating that a later ThermalEnd−1 by compromising the females’ BCI in winter (Figure 1D) contributed to delay the calving dates the following calving season (Table 4).

Table 3.

Competing LMMs of the effects of plant phenology on calving date of a semidomesticated reindeer population in Kaamanen, northern Finland

Rank	Models	AIC	df	AIC weights	ΔAIC
	1970–2015
1	zDENSt−1 + zPMt−1 + dBCIWintert × zThermalEndt−1 + dThermalStartt−1	4373.78	9	0.71	0.00
2	zDENSt−1 + zPMt−1 + dBCIWintert × zThermalEndt−1 + dBCIWintert × dThermalStartt−1	4375.56	10	0.29	1.78
	1997–2015
1	dDENSt−1 + zPMt−1 + dBCISpringt × zVegOnsett−1 + zVegEndt−1	1658.53	9	0.62	0.00
2	dDENSt−1 + zPMt−1 + dBCISpringt × zVegOnsett−1 + dBCISpringt × zVegEndt−1	1659.50	10	0.38	0.98

The analyses were repeated over the 2 study periods: the whole study period (1970–2015) and the last part of the study period (1997–2015). All models included the female identity as a random factor on the intercept, as well as the population density (DENS) and the proportion of males in the herd (PM) as fixed effects. A lowercase “d” in front of the variable indicated that the variable was used “detrended,” whereas a lowercase “z” indicated that the variable was used standardized. ThermalStart and ThermalEnd represented, respectively, the start and end of the thermal growing season. The onset and end of the vegetative growing season were depicted by VegOnset and VegEnd, respectively. BCIWinter was the BCI of females averaged over the winter (from December to March), whereas BCISpring was the BCI of females, averaged in spring (April–May). A cross “×” indicated an interaction term between the 2 variables. The models presented in the table are the competing models retained in explaining calving date, that is, with ΔAIC < 2 (see text for details).

Table 4.

Model-averaged estimates of fixed effects from the LMMs of calving date of a reindeer population in relation to plant phenology in Kaamanen, northern Finland.

Variable	Estimate	Unconditional SE	Nbr models	95% CI
1970–2015
zDENSt−1	0.23	0.03	2	0.17, 0.29
zPMt−1	−0.09	0.02	2	−0.14, −0.04
dBCIWintert	−0.20	0.03	2	−0.25, −0.14
dThermalStartt−1	0.10	0.02	2	0.05, 0.15
zThermalEndt−1	0.12	0.02	2	0.07, 0.16
dBCIWintert × zThermalEndt−1	0.06	0.02	2	0.02, 0.10
dBCIWintert × dThermalStartt−1	−0.01	0.02	1	−0.06, 0.03
1997–2015
dDENSt−1	−0.08	0.07	2	−0.21, 0.06
zPMt−1	−0.06	0.04	2	−0.14, 0.03
dBCISpringt	−0.33	0.06	2	−0.45, −0.22
zVegOnsett−1	0.10	0.05	2	0.0006, 0.20
zVegEndt−1	0.47	0.04	2	0.38, 0.55
dBCISpringt × zVegOnsett−1	−0.18	0.06	2	−0.31, −0.06
dBCISpringt × zVegEndt−1	−0.05	0.04	1	−0.13, 0.04

The estimates in bold type were deemed important (whose 95% CI excluded 0) in explaining calving date. “Nbr models” is the number of models (out of the 2 best models in Table 3) including that variable. The analyses were repeated over the 2 study periods: the whole study period (1970–2015) and the last part of the study period (1997–2015). All models included the female identity as a random factor on the intercept. A lowercase “d” in front of the variable indicated that the variable was used “detrended,” whereas a lowercase “z” indicated that the variable was used standardized. The predictor variables were the following: DENS, PM, start and end of the thermal growing season (ThermalStart and ThermalEnd, respectively), onset and end of the vegetative growing season (VegOnset and VegEnd, respectively), BCI of females averaged over the winter (from December to March, BCIWinter) and BCI of females, averaged in spring (April–May, BCISpring). A cross “×” indicated an interaction term between the 2 variables. DENS, population density; PM, proportion of males in the herd; SE, standard error.

For the last part of the study period (1997–2015), a total of 606 observations were analyzed and 2 competing models were also selected (Table 3). The model averaging showed that an earlier onset (VegOnset−1) and end (VegEnd−1) of the vegetative growing season the previous calendar year resulted in earlier calving dates the next calving season (Figure 1E, F, Table 4). Also, a higher females’ BCI in spring (BCISpring) just before the calving season induced earlier calving dates (Figure 1G) and the statistically significant interaction term between BCISpring and VegOnset−1 indicated that a later VegOnset−1 the previous spring season by compromising the females’ BCISpring the next year (Figure 1H) lead to late calving dates (Table 4).

Among the 4 measures of time lags, none presented a statistically significant temporal trend, meaning that the time lags between calving date and start of the thermal (ThermalStart) and vegetative (VegOnset) growing season remained unchanged over the years, as well as between calving date and end of the thermal (ThermalEnd−1) and vegetative (VegEnd−1) growing season, being for one or the other of the 2 study periods.

Fitness consequences

The corrected calves’ birth weight (BirthWeight) was increasing over time (b = 0.02, 95% CI [0.01, 0.03]), as well as the calves’ first-summer survival (Surv= 0.06, 95% CI [0.05, 0.07]). BirthWeight was therefore used detrended in the following analyses, whereas Surv was unaltered to preserve its binomial structure. For the whole study period (1970–2015), the BirthWeight was higher when the start of the previous thermal growing season (ThermalStart−1) was earlier, when the length of the previous thermal growing season (LTGS−1) was shorter and when the time lag between calving date and onset of the thermal growing season (TLThermalStart) the previous calendar year was longer (Table 5). A longer TLThermalStart can be due to an earlier calving date or a later ThermalStart−1. Given that ThermalStart−1 was occurring earlier over time and that BirthWeight was positively affected by ThermalStart−1, it would indicate that an earlier calving date more likely caused a higher birth weight of calves. The calves’ first-summer survival (Surv) was higher when ThermalStart was earlier, when ThermalEnd−1 was earlier and finally when LTGS−1 was shorter (Table 5). None of the measures of time lags (TLThermalStart or TLThermalEnd) affected Surv.

Table 5.

Fitness consequences of the lagged effects of plant phenology on the birth weight (BirthWeight) and first-summer survival (Surv) of calves of the Kutuharju reindeer herd in northern Finland.

	BirthWeightt			Survt
	b	95% CI		b	95% CI
1970–2015			1970–2015
dThermalStartt−1	−0.07	−0.11, −0.03	dThermalStartt	−0.20	−0.32, −0.07
zTLThermalStart	0.09	0.04, 0.13	zThermalEndt−1	−0.13	−0.25, −0.01
dLTGSt−1	−0.04	−0.08, −0.002	dLTGSt−1	−0.39	−0.52, −0.26
1997–2015			1997–2015
zVegOnsett−1	−0.09	−0.15, −0.02	zVegOnsett	−0.40	−0.71, −0.11
zThermalEndt−1	−0.18	−0.25, −0.12
zLTGSt−1	−0.17	−0.23, −0.10	zTLThermalStart	0.34	0.08, 0.60
zTLThermalEnd	0.20	0.14, 0.27	zTLThermalEnd	−0.31	−0.60, −0.02

The analyses were repeated over the 2 study periods: the whole study period (1970–2015) and the last part of the study period (1997–2015). All models included the female identity as a random factor on the intercept. A lowercase “d” in front of the variable indicated that the variable was used “detrended,” whereas a lowercase “z” indicated that the variable was used standardized. The predictor variables were the following: start and end of the thermal growing season (ThermalStart and ThermalEnd, respectively), onset of the vegetative growing season (VegOnset), LTGS, time lag between calving date and ThermalStart (TLThermalStart) and time lag between calving date and ThermalEnd the previous calendar year (TLThermalEnd).

For the last part of the study period (1997–2015), heavier new-born calves (corrected by the calving date) were observed when ThermalEnd−1 was earlier, LTGS−1 was shorter, the onset of the vegetative growing season (VegOnset−1) was earlier, and the TLThermalEnd was longer (Table 5). A longer TLThermalEnd indicates either a later calving date or an earlier ThermalEnd−1. Given that BirthWeight is negatively correlated to ThermalEnd−1 and that early-born calves are also heavier (b = −0.02, 95% CI [−0.03, −0.02]), it would most likely indicate that an earlier ThermalEnd−1 lead to heavier calves. A higher calves’ first-summer survival was observed when VegOnset was earlier, the time lag between calving date and ThermalStart (TLThermalStart) was longer and the time lag between calving date and ThermalEnd−1 (TLThermalEnd) was shorter (Table 5). Given that Surv was negatively correlated with VegOnset and that early-born calves have a higher first-summer survival (b = −0.06, 95% CI [−0.07, −0.04]), the positive correlation between TLThermalEnd and Surv would most likely be due to an earlier calving date. A shorter TLThermalEnd indicates either an earlier calving date or a later ThermalEnd−1. An earlier calving date would more likely contribute to a higher Surv as Surv was negatively affected by ThermalEnd−1 in the whole study period. As expected from the literature, heavier new-born calves had a higher first-summer survival (b = 0.94, 95% CI [0.79, 1.09]).

Discussion

As predicted, the time lag between calving date of the Kutuharju reindeer population and onset of the thermal or vegetative growing season did not change over time. As already reported in other ungulate populations, the onset in spring of both the thermal and vegetative growing season affected calving date (Figure 1A, E), with a lagged effect of 1 year (Table 4). Calving date was also affected by the end, in autumn, of the previous thermal and vegetative growing season (Figure 1B, F). Although the end of the thermal and vegetative growing season did not significantly change over time (Table 1), a delay in ThermalEnd caused females’ physical condition in winter to be deteriorated (Figure 1D) and then caused a delay in calving date the following calving season (Table 4). While an earlier ThermalStart and VegOnset was beneficial to females’ reproductive success with an increase in both calves’ birth weight and calves’ first-summer survival, a delay in ThermalEnd was detrimental (Table 5). This is consistent with the suggestion that potential lagged effects of plant phenology arising in late summer and autumn, when females store nutrient reserves to finance reproduction (Barboza and Parker 2008) can also lead to consequences on reproductive success than at parturition due to the capital breeder strategy of Rangifer and its reliance on maternal capital (Taillon et al. 2013; Gustine et al. 2017; Williams et al. 2017).

Keeping up with the onset of spring phenology

Contrary to temporal trends reported in other northern ecosystems (Oberbauer et al. 2013; Park et al. 2016) and despite warmer temperatures in April and May in our study area over the last 45 years (see Paoli et al. 2018), the VegOnset has not advanced (Table 1). An earlier ThermalStart was, however, reported for the 2 study periods, with a statistically significant advancement of 11 days from 1970 to 2015 (−0.24 days.year−1) and of 15 days from 1997 to 2015 (−0.84 days.year−1, Table 1). Despite this temporal trend and as predicted in our first hypothesis, there was no evidence of a developing time lag between the parturition date and the onset of spring vegetative phenology in this semidomesticated reindeer population, as previously reported in Alaskan caribou (Gustine et al. 2017). Although the rate of advancement for calving date was much lower for those 2 study periods, −0.14 days.year−1 for the whole study period and −0.51 days.year−1 for the last part of the study period (Table 1), the results indicated that the pace of advancement of reindeer calving date of the Kutuharju herd was so far keeping up with the pace of advancement in onset of the spring green-up in that area. A high correlation was found between the mean temperature in April–May in our study area and both the onset of the thermal and vegetative growing season (respectively, r = −0.50, 95% CI [−0.69, −0.24] and r = −0.86, 95% CI [−0.95, −0.57]). Moreover, earlier calving dates in the same reindeer population have recently been found to be associated with warmer temperatures in April–May (see Paoli et al. 2018). The mean temperature in April–May being a good predictor of onset of the spring vegetative season, earlier calving dates would have thus concurred with earlier spring green-up. The physiological mechanisms of how such fine adjustment in reindeer calving date is however uncertain. The females’ physical condition in spring was improved when the onset of the vegetative growing season 1 year ago was earlier (Figure 1H). A lagged effect of plant phenology in spring the previous year acting on BCISpring the next year might thus be involved. Perhaps a compensatory mechanism, as highlighted in numerous species (Rachlow and Bowyer 1991; Cameron et al. 1993; Flydal and Reimers 2002), has favored shortened gestation lengths when females’ physical condition in late pregnancy was improved following warmer temperatures in April–May (b = 0.68, 95% CI [0.40, 0.97]). On top of that, an earlier VegOnset−1 by allowing the females to calf earlier the previous year might have benefit them as being able to nurse on a high-quality forage for a longer period of time (Festa-Bianchet 1988; Rachlow and Bowyer 1994). As a cascading effect, their body condition the next spring was improved and they could again calf earlier (Figure 1G). The previous year’s reproductive status has already been shown to affect variation in birth date (Guinness et al. 1978) but this study demonstrates that, as suspected, it was indeed acting through female body condition. Furthermore, a bottom-up trophic effect caused by spring plant phenology the previous year might explain why the effect on females’ physical condition in April–May could be observed 1 year after (as shown for the AO index on reindeer population growth in Aanes et al. 2002).

Calving date constrained by the plant phenology in autumn

Confirming our second hypothesis, the calving date was also affected by ThermalEnd−1 and VegEnd−1 (Figure 1B, F). The timing of birth in a Pyrenean chamois Rupicapra pyrenaica pyrenaica population was shown to be shaped by onset of autumn more than onset of spring and concurs with our findings (Kourkgy et al. 2016). Although this semidomesticated reindeer population receives supplemental feeding in late winter to buffer the harshness of winter conditions, it suggests that both wild and domesticated animal populations could also be influenced by the end of the plant growing season in autumn. More studies on autumn phenology are however needed to state this with certainty. The timing of parturition in ungulate species is related to either autumn or over-winter physical condition of females (Cameron et al. 1993; Cook et al. 2004) and it is assumed that variable calving dates might result from nutritionally related differences in gestation length (Bergerud 1975; Berger 1992). In our study, the females’ physical condition in winter (from December to March), when females rely on maternal capita to support gestation costs, was improved when the end of the previous vegetative growing season in autumn was later (Figure 1D). Caribou and reindeer are capital breeders (Taillon et al. 2013), with females building nutritional reserves in late summer/autumn to finance reproduction (Barboza and Parker 2008). As such, the end of the vegetative season in the autumn by influencing females’ physical condition in winter and therefore maternal body stores for fetal development (Barboza and Parker 2008; Williams et al. 2017) would also be of great importance in influencing future calving dates.

Contrary to our third hypothesis, however, the time lag between calving date and either ThermalEnd−1 or VegEnd−1 did not significantly change over time. Too few years available to measure VegEnd−1 might explain why, despite a high rate of change (−0.41 days.year−1), the temporal trend was not statistically significant. Furthermore, the correlation of 55% between ThermalEnd−1 and VegEnd−1 confirmed that the autumn phenology is much harder to estimate from climatic variables as temporal change in the autumn seems less apparent and with more heterogeneous patterns (Walther et al. 2002; Walther and Linderholm 2006). The thermal growing season indices can only estimate the real growing season to some extent (Walther and Linderholm 2006), given that the birch phenological data was limited to a shorter period. The use of various indices to estimate the thermal growing season can also lead to remarkable differences in the long-term (Walther and Linderholm 2006). That VegOnset−1 and VegEnd−1 better explained the variation in calving date during the second part of the study period (1997–2015) than ThermalStart−1 and ThermalEnd−1 (Table 4) also points out that actual observations on birch phenology are more relevant at depicting biological events than indices of the thermal growing season as some discrepancies might be observed between the thermal and vegetative growing season (Walther and Linderholm 2006).

Heavier calves and calves with a better first-summer survival were observed when the start of the thermal growing season was earlier for the whole study period (1970–2015) and when the onset of the vegetation growth was earlier for the last part of the study period (1997–2015). A 1-year lagged effect of ThermalStart−1 and VegOnset−1 on BirthWeight was observed (Table 5), certainly because as stated earlier, an earlier VegOnset−1 or ThermalStart−1 would allow the females to calf earlier the previous year, to free themselves from the “burden” of lactation earlier in the summer so that they are able to replenish their body fat reserves faster in summer (Skogland 1983). As a calf’s birth weight is highly dependent on her mother’s food intake and body weight (Bergerud 1975; Cameron et al. 1993; Côté and Festa-Bianchet 2001; Solberg et al. 2007), an earlier VegOnset−1 or ThermalStart−1 would have contributed to improve the females’ body weight at the time of calving (Pettorelli et al. 2005) and, in turn, to produce heavier calves. A higher calves’ survival was also reported when the onset of spring green-up was earlier the same year (ThermalStart or VegOnset, Table 5). Early-born calves are also heavier and then enjoy a higher survival because they benefit not only a longer growing season but also from the high-quality resources available for growth (Bunnell 1982; Festa-Bianchet 1988; Feder et al. 2008). In short summer growth environments, the transfer of energy and protein through maternal milk in sufficient quantity and in a short amount of time is critical for calf’s growth and survival (Sadleir 1969, p. 321). Ruminants, such as reindeer, are believed to be more selective on forage quality than forage quantity (Tveraa et al. 2013) because females calving earlier produce milk from forage of early phenology, with lower fiber content and higher digestibility (Rachlow and Bowyer 1991). A milk with a higher protein content from a high-quality forage is thus crucial for females’ reproductive success (Festa-Bianchet 1988; Rachlow and Bowyer 1994) by increasing the quality and quantity of resources available for calves’ growth. As previously shown in reindeer (Tveraa et al. 2013), the females’ reproductive success in our study was thus positively affected by the onset of spring green-up, certainly through a positive effect of forage quality on the females’ physical condition during late gestation and, in turn, on calving date and milk production (Rachlow and Bowyer 1991; Cameron et al. 1993).

Despite longer thermal growing seasons benefiting the calves’ birth weight and survival (Table 5), a “too late” end of the thermal growing season seemed to negatively affect those females’ reproductive success components (Table 5). Indeed, a decrease in calves’ birth weight was reported for the last part of the study period (1997–2015) and in calves’ first-summer survival for the whole study period (1970–2015) when ThermalEnd−1 was occurring later. This finding seems counter-intuitive given that a later ThermalEnd−1 would favor a longer growing period for the calves before facing their first winter. Since the main determinant of calving date in ungulates is the conception date (Holand et al. 2002; Clements et al. 2011), an inadequate nutrition in autumn would be detrimental to females’ body condition and cause a subsequent delay in calving dates (Cook et al. 2004). Late-born calves were then reported to be lighter and with a lower first-summer survival. Mushrooms are an important and preferable source of food for reindeer during the late summer (end of July, beginning of August throughout September) up to early winter (October–December) (Kitti et al. 2006; Inga 2013), with their abundance having a considerable effect on reindeer body condition (JK: personal communication). The ideal growing conditions for mushrooms demand humidity and relatively cool temperatures. A warming trend in September–October (b = 0.04, 95% CI [0.02, 0.06]) coinciding with a later ThermalEnd−1 (r = 0.61, 95% CI [0.39, 0.77]) might therefore be detrimental to mushrooms development due to heat stress. As a consequence, the mean autumn fruiting date for both mycorrhizal and saprotrophic fungi in Norway is now later (Kauserud et al. 2008, 2012), with an accelerated delay of fruiting in the last 20 years that has coincided with global warming (Kauserud et al. 2008). In those exceptionally warm years, when the mushroom abundance is reduced and delayed due to deteriorated conditions for their growth, females might not be able to gain as much fat as in regular years and a decreased physical condition in autumn was then observed following a later ThermalEnd−1 (b = −0.06, 95% CI [−0.09, −0.03]). The level of fatness and protein reserves attained by autumn control conception rate (Cameron et al. 1993; Reimers 1997; Barboza and Parker 2008) and thus affect reproductive effort (Weladji et al. 2002). The coincidence of warmer temperatures in September–October with a late ThermalEnd−1 would have thus caused an impaired physical condition of females near the breeding time leading to a delay in calving date (Figure 1B) and a decrease in calves’ BirthWeight and Surv the following calving season (Table 5). In summary and contrary to our expectations, a lengthening of the thermal growing season (Table 1), consistent with long-term trends reported in northern ecosystems (Oberbauer et al. 2013; Park et al. 2016) is not necessarily beneficial for reindeer if a later end is accompanied by a reduced abundance of mushrooms. Our study can therefore be added to the growing body of literature on the MMH and how species are able to adapt to climate change. The consequences of developing lagged effects of plant phenology on species’ reproductive success are, however, still much contrasted among ungulates, an area that deserve further investigation.