The cooler the better? Indirect effect of spring-summer temperature on fecundity in a capital breeder

Female reproductive performance is a central component of ungulate population dynamics, and it can be influenced by individual, social and environmental factors. Researchers have often assumed direct effects of different predictors on reproduction, yet more complex relationships should be considered when investigating temporal variations in life history traits within a broader eco-evolution-ary context. In this study, we explored direct effects of individual, social and environmental predictors on female reproductive performance, and investigated potential causal chains among variables. We analysed the variation in fecundity, measured as the probability of being pregnant, in 215 adult female deer Cervus elaphus culled on the Italian Alps, with respect to age, body mass, kidney fat, jaw length, lactation status, population size, temperature and precipitation in spring-summer, temperature and snow depth in winter, and the delayed effect of spring-summer temperature. We used random forest and logistic regression models to select variables whose direct effects best explained variation in fecundity. Path analysis was used to test for alternative hypotheses of direct / indirect effects between pre-selected weather (spring-summer temperature) and individual (age, KFI) predictors. The most important direct predictors of fecundity were age, kidney fat and the interaction between kidney fat and spring-summer temperature. Path analysis supported the hypothesis that higher spring-summer temperature had negative, indirect effects on the probability of being pregnant, mediated by decreasing values of kidney fat index. Our study revealed some complex, cause-effect relationships between weather stochasticity, body condition and reproduction, possibly suggesting a conditional trade-off between opportunity for reproduction and survival, and emphasizing how environmental variations and individual characteristics may interact to shape life history traits in ungulate populations.


INTRODUCTION
Reproduction is a central component of animal population dynamics. Several individual, social and environmental factors such as age, body condition, population density and climatic conditions, are known to affect reproductive performance of females (Gaillard et al. 2000). The importance of each factor may vary among species, owing to different life history strategies (Coulson et al. 2000, Hamel et al. 2010, and within species, depending on different habitat characteristics (Balbontín and Ferrer 2008).
In many species of birds and mammals, female reproduction is strongly related to age in a nonlinear fashion, and different reproductive parameters generally increase from sexual maturity to prime age, and then decline later in life (Nussey et al. 2008). The age effect is intimately associated with the long-term variations in body condition that occur over an individual's lifetime, which influence its reproductive performance (Nussey et al. 2011, Flajšman et al. 2017. Body condition may also undergo great short-term variations, in response to internal and/or external pressures. For example, reproductive events may impose high energetic costs on females and negatively affect individual condition, thus lowering the probability of reproducing in the following breeding season. This pattern has been observed in mammals and in other taxa (Gustafsson and Sutherland 1988, Yurewicz and Wilbur 2004, Hamel et al. 2010, especially in food-limited populations (Clutton-Brock et al. 1989). Density-dependent availability of food supply (Bonenfant et al. 2009) or climate effects may in fact trigger short-term changes in body condition and, in turn variation in reproductive performance. For example, several studies showed evidence for negative responses of reproductive traits to increasing ambient temperature in different taxa, possibly mediated by variations in individual conditions (Grazer and Martin 2012). Individual, social and environmental variables may thus influence female reproductive traits in a complex manner, operating through pathways that include both direct and indirect relationships between different factors, possibly reflecting adaptive responses to optimize the trade-offs between reproduction and survival in different environmental conditions (Sand 1996). A deeper understanding of the adaptive mechanisms underlying variations in female reproductive traits may benefit from explicitly assuming non-independence among predictive variables. Recent studies have supported the occurrence of more complex causal relationships among potential drivers of reproduction: in Norwegian red deer Cervus elaphus, for example, body condition during summer is indirectly affected by climate through plant phenology, as higher spring temperatures accelerate plant development (Mysterud et al. 2008). This, in turn suggests that climatic variables may exert an important indirect effect on deer reproductive performance. Indirect consequences of rising temperatures on life history traits are increasingly suggested to occur also in Alpine ungulates. In highly seasonal, energy-limited mountain temperate environments, fat reserves are mainly deposited before autumn: warmer spring-summer periods may cause a reduction in food acquisition through several mechanisms (Pettorelli et al. 2007, Mason et al. 2014, 2017 and thus negatively impact on body condition and, possibly, reproduction in capital breeders (Rughetti and Festa-Bianchet 2012). Different hypotheses may be put forward to explain this pattern. A first hypothesis suggests that warm temperatures in spring-summer may accelerate plant development and reduce the availability of high-quality food resources over the summer (Pettorelli et al. 2007).
An alternative hypothesis suggests that in mountain dwelling ungulates such as the chamois Rupicapra rupicapra, body condition in autumn is not limited by summer resource availability: rather, the high spring-summer temperatures may reduce the time spent foraging before autumn thus limiting the ability of individuals to acquire resources (Mason et al. 2014). More recently, for the Alpine ibex Capra ibex it has been suggested that warmer temperatures in summer would force animals to thermoregulate by using less productive areas at higher elevations, and thus consume lower quality forage, without compensating their foraging effort (Mason et al. 2017). Despite their potential demographic and evolutionary significance, however, the indirect effects between individual, social and climatic variables on the reproductive performance of female ungulates have received comparatively little attention. Owing to its wide geographical distribution, the direct role of individual, social and environmental factors on female reproductive performance has been largely investigated in red deer. The importance of age on hind reproduction, for example, has been supported in several populations (e.g. Albon et al. 1986, Bertouille andde Crombrugghe 2002), while the role of other individual drivers appears to be less consistent, varying with latitude and habitats. Body mass and fat reserves had positive effects on hind pregnancy probability in Scottish populations (Albon et al. 1986), but no relationship was found in central Europe (Borowik et al. 2016). Body size also showed contrasting effects on pregnancy probability in northern and central Europe (Mitchell and Brown 1974, Albon et al. 1986, Bertouille and de Crombrugghe 2002. The negative effect of lactation on the probability of being pregnant in the following reproductive season is strong in foodlimited environments (Clutton-Brock et al. 1989), but little information is available in rich environments. Similarly, the direct effects of social and environmental variables on hind reproduction show site-dependent variations. Negative density-dependent relationships between pregnancy probability and body mass were found in Scotland (Albon et al. 1983), while no densitydependent effects were found in central Europe (Bonenfant et al. 2002;Borowik et al. 2016). In Scottish deer, summer precipitation, winter temperature and snow negatively influenced female fecundity Clutton-Brock 1988, Langvatn et al. 1996), whereas adult pregnancy ratio in Rocky Mountain elk increased following summers with high precipitation (Proffitt et al. 2014).
More recently, no significant effects of winter and summer temperature on deer fecundity were found in Poland (Borowik et al. 2016). Overall, the direct effects of individual, social and climatic variables on hind reproduction do not show congruent patterns over large geographic scale, and little is known about their potential indirect effects. A deeper understanding of the adaptive responses to environmental conditions should thus account for more complex interactions (cf. Stopher et al. 2014).
In this paper we first investigate the direct effect of individual, social and environmental factors on reproductive performance in an Alpine population of red deer, as limited information is available in this environment. We then explore potential causal pathways affecting female fecundity, accounting for the occurrence of direct and indirect relationships between variables (sensu Shipley 2016). In particular, we test hypotheses about the potential effect of weather conditions on fecundity, already suggested for other Alpine ungulates, which indicate that increasing springsummer temperatures could negatively and indirectly affect reproductive performance through a decline in body condition (Pettorelli et al. 2007, Rughetti and Festa-Bianchet 2012, Mason et al. 2014).

Study area and population
The study site 'Valfurva' lies in the northwestern part of the Stelvio National Park, within the Province of Sondrio, central Italian Alps (10°25′N, 46°27′E). Valfurva is the wintering site of a large population management unit for red deer, defined by tracking GPS-collared individuals (cf. Corlatti et al. 2016), and extends over 4975 ha between 1200 and 2400 m a.s.l. (Fig. 1). About 73% of its surface is dominated by spruce Picea abies, larch Larix decidua and Stone pine Pinus cembra forests, while the remaining 27% consists of open areas with mesic meadows at lower elevations, and Alpine grasslands of Carex spp., Festuca halleri and Sesleria coerulea above the treeline. The climate is alpine continental, with mean temperatures between 15.7°C in July and -2.8°C in January and yearly precipitation of about 765 mm. Between 2011 and 2015, the winter density of the red deer population in Valfurva was about 27.4 ind./km 2 ( 2.5 SD) (Corlatti et al. 2016). The large increase in deer density that occurred in the Park over the last two decades, owing to the good environmental conditions within the Park and possibly the absence of hunting pressure inside the protected area, severely impacted on forest regeneration, on agricultural activities, and on the ecosystem biodiversity. In 2011 the National Park Agency therefore started a culling program aimed at reducing population density: given the initial high density of deer (about 31 ind. /km 2 in winter), the effects of culling were apparent only since 2015 (cf. Corlatti et al. 2016). In 4 years,

Data collection
In the literature on ungulate ecology, different terms such as 'pregnancy probability', 'fecundity' and 'fertility' are often used interchangeably to indicate reproductive performance. For the sake of clarity, hereafter we refer to 'fecundity' as to the probability of being pregnant. To investigate deer fecundity, we focused on adult females (2+ years of age) only, as very few yearlings breed. All females were brought to a control centre within 2 hours from culling. For some individuals, it was not possible to collect all the parameters needed for the analysis (see below), and our sample size reduced to n = 215 adults. Pregnancy status was investigated through the presence / absence of corpora lutea (i.e. endocrine structures that develop from ovarian follicles during the luteal phase of the oestrous cycle) by dissecting ovaries: for nearly all individuals, pregnancy was confirmed by the presence of foetuses, which always occurred in conjunction with corpora lutea, thus supporting the suitability of this last parameter to assess reproductive status. Potential drivers of hind fecundity included individual and external (social and environmental) variables. Individual variables comprised age (estimated by counting the cementum annuli), lactation status (evaluated through the presence of milk, or cutting into the udder), jaw length (in mm, measured by means of an electronic calliper), dressed body mass (in kg) and kidney fat index (KFI, measured following Riney 1955). We considered KFI as the most suitable proxy for individual condition in red deer (Riney 1955), as it was the most direct proxy of body fat reserves available to us. External variables included population size in the previous spring (based on mark-resight estimates: Corlatti et al. 2016) and meteorological variables. On the Alps, deer fertilization occurs at the beginning of October, and fat reserves are mainly deposited in spring-summer: because in this environment warmer temperatures are linked to anticipated vegetation growth (Pettorelli et al. 2007), to predict fecundity in each year t we used the mean daily temperature (in °C) during the entire vegetative season, i.e. from April -when vegetation growth begins -to September -when the seasonal course of shoot biomass drops -at year t. Additionally, we also considered the effects of other climatic variables on fecundity: the cumulative precipitation (in mm) in spring-summer (April-September) at year t; the mean daily temperature (in °C) and the mean snow depth (in cm) during winter (January-March) at year t; the mean daily spring-summer temperature (in °C) at year t-1. Data were retrieved from a meteorological station within the study area (cumulative precipitation in spring-summer at year t-1 was not available for the year 2010, hence this variable was not included in the analysis).

Statistical analysis
To investigate the direct relationships between female fecundity and individual, social or environmental variables, we used both non-parametric and parametric approaches, as different selection criteria can give different results even with the same dataset (Gotelli and Ellison 2013).
Prior to analysis, all continuous explanatory variables were standardized by subtracting each sample's mean and then dividing by the sample's standard deviation, to return comparable coefficients and reduce issue of collinearity in presence of interaction terms. To investigate direct and indirect causal relationships among selected variables (age, KFI, spring-summer temperature), we used path analysis (Wright 1934).  (Kursa and Rudnicki 2010). This nonparametric approach allows to identify ecologically important predictors of fecundity, and offers some advantages over traditional variable selection procedures, as it is more robust to collinearity (Cutler et al. 2007). The interpretation of high-order interactions in random forest -based approaches, however, is not immediate, and the method is not suitable for hypothesis testing or for identifying ecologically important subsets of variables in the way model selection does (Cutler et al. 2007).
Direct relationships: parametric approach. -We proceeded exploring direct effects of additive and interactive combinations of individual and environmental variables on deer fecundity, by applying a parametric, information-theoretic (IT) model selection approach (Burnham and Anderson 2002) on a set of biologically plausible competing models explicitly tied to underlying mechanisms linking climate and individual/social variables to fecundity. Specifically, Binomial linear regression models were fitted with logit link function: Fecundity i represents the pregnancy status (0/1) of individual i at time t. Individual age i at time t was fitted as a quadratic term (to account for non-linear effect on fecundity) in each model, as preliminary analysis showed that it consistently improved models' fit. The fixed variable X 1 i was represented by either individual body mass i , KFI i , jaw length i , lactation status i , or by population ¿ ¿i ¿ , all referred to year t, while the variable X 2 i was represented by mean spring−summer temperature i at year t-1. When assessing the relative importance of variables in the information-theoretic framework, it is important to achieve a balance in the number of models that contain each explanatory variable (Burnham and Anderson 2002). Therefore, to explore the effect of every biologically plausible combination of predictors X 1 i and X 2 i a set of 24 models was generated (Supporting Information Appendix S1: Table S1). To account for temporal variation in the value of body mass and KFI, all the models that included these two variables also included individual shooting date i as a covariate, that is the number of days elapsed from October 30 of each year (i.e. the first day of shooting) to individual culling. A preliminary analysis using AICc showed that fixed-effect models (eqn. 1) consistently outperformed mixed-effect models fitted with an observation-level random intercept (to account for unexplained heterogeneity among subjects: Harrison, 2015). Prior to analysis, a matrix based on Pearson's correlation coefficient (r p ) on standardized continuous variables was built to identify potential issues of collinearity. In case of severe collinearity (r p > 0.7, Dormann et al. 2013), variables were never fitted in the same model, to avoid bias in parameter estimation. Model parameters were estimated using maximum likelihood.
All fitted models (n=24) were subsequently ranked based on their AICc values and retained in the final candidate set if they had ΔAICc ≤ 2 (Burnham and Anderson 2002). The explained variance of models in the candidate set was measured by studying the Nagelkerke's pseudo R 2 , and the goodness-of-fit was assessed by calculating the p-value associated with the Hosmer and Lemeshow test statistic (Hosmer and Lemeshow 2000). Additionally, the predictive accuracy of candidate models was measured using the area under the receiver operating characteristic (ROC) curve (AUC) which, in this study, refers to the ability of a given model to discriminate between pregnant and non-pregnant females. Bootstrapping is arguably the best alternative for obtaining predictive ability measures, as it provides stable estimates with low bias, especially with small sample sizes (Steyerberg et al. 2001). To estimate the AUC for each model in the candidate set, we thus performed an internal validation using 1000 bootstrap samples following Harrel et al. (1996). Causal relationships: path analysis. -Path analysis (Wright 1934) was used to investigate potential causal relationships between selected climate (spring-summer temperature) and individual (age, KFI) variables and their direct /indirect effects on fecundity. Path analysis requires the creation of diagrams that illustrate the hypothesized relationships among the selected variables (Gotelli and Ellison 2013), thus we first built three directed acyclic graphs to represent alternative causal models (Fig. 2) (Kursa and Rudnicki 2010) to investigate variable importance based on random forest classification. The 'glm' function was used to fit logistic regression models; model selection was performed using the 'MuMIn' package (Bartoń 2015), while the bootstrap internal validation was conducted with the 'rms' package (Harrel 2017). Causal relationships using path analysis were tested with the package 'lavaan' (Rosseel 2016).

Direct relationships
The all-relevant feature selection based on the wrapper algorithm around random forest classification showed that age and KFI were the most relevant variables directly related to fecundity in our study population (Table 1, Supporting Information Appendix S2: Fig. S1).
The correlation matrix did not suggest severe issues of collinearity, except between temperature and population size (r p = -0.88) or winter snow (r p = -0.88): these two pairs of variables were thus never included in the same model, to avoid bias in parameter estimation. The model selection procedure retained only one model as candidate to explain variation in fecundity of adult hinds (ΔAICc ≤ 2, Table 2; see also Supporting Information Appendix S1:  3). As the 2AICc cut-off might arguably be considered overly conservative, in the supporting information we show that consistent results are reached when averaging models using less restrictive cut-offs (ΔAICc ≤ 4, see Supporting Information Appendix S1: Table S3).

Causal relationships
Path analysis, performed to discriminate alternative causal relationships among the selected meteorological and individual variables, showed that two of the structures hypothesized in our directed acyclic graphs (Model a and Model c) satisfactorily fitted the correlational structure of the sample data ( 2 value > 0.05, RMSEA value ≤ 0.6, Table 4). The Satorra-Bentler test did not show a significant difference between the two models ( 2 difference = 0.600, p-value = 0.439). Following the principle of parsimony, this result suggests that the additional path assumed in Model c (Fig. 2) may be considered uninformative (as confirmed by the estimates reported in Supporting Information Appendix S2: Fig. S3). Model a was thus selected as the best model, supporting the hypothesis that temperature had an indirect effect on fecundity, mediated by KFI: the standardized path coefficients (Fig. 4) show that increasing spring-summer temperatures negatively impacted on KFI and, in turn, on pregnancy status (indirect effect: estimate = -0.054, p-value = 0.020). The direct non-linear effect of age on fecundity, however, was stronger than the indirect effect of temperature (Fig. 4).

DISCUSSION
The parametric and non-parametric variable selection procedures consistently supported the direct role of age and kidney fat index (KFI) to explain fecundity variation within our study population.
Age had a strong, non-linear effect on fecundity probability, with relatively low values for the extreme age-classes and high values for the intermediate age-classes, whereas increasing values of KFI had positive effects. In both approaches, the direct, independent effects of climatic variables were negligible. Spring-summer temperature was not retained by the non-parametric analysis since our approach only evaluates the importance of individual effects, but its influence became significant in parametric models when fitted as an interaction with KFI. The result of this interaction suggests that a mediation effect might occur between these variables (Baron and Kenny 1986): indeed, air temperature and KFI were assumed to be causally related, as warm spring- summer temperatures might lead to relatively poorer body conditions because of anticipated plant development (Mysterud et al. 2008). The path analysis went beyond the simplistic regression approach, allowing to confirm that increasing spring-summer temperature negatively, indirectly affected hind fecundity through a negative effect on KFI values, thus allowing a deeper understanding of the adaptive response of females to different environmental conditions.
The non-linear relationship between age and adult fecundity has been supported in several deer populations (e.g. Albon et al. 1986, Stewart et al. 2005, Morano et al. 2013 and in other ungulates (e.g. mountain goats Oreamnos americanus, Côté and Festa-Bianchet 2001, Soay sheep Ovis aries, Tavecchia et al. 2005). Younger and older individuals are less likely to become pregnant than prime-aged individuals, suggesting an age-dependent cost of investment in reproduction. While the relatively low pregnancy rates in young age classes may be explained by the necessity to reach a threshold body mass to attain primiparity (Gaillard et al. 2000), there is still limited understanding of the causes of decline in fecundity in old females. In red deer, however, the disposable soma and antagonistic pleiotropy theories of senescence have received some support (Nussey et al. 2006). The  (2013). The limitations in the use of KFI as a ratio index have been discussed by Serrano et al (2008), who suggested the use of residuals of the linear regression log(kidney fat) ~ log(kidney weight) instead, as they have the advantage of being size independent. Other works, however, pointed out that "size independence does not necessarily mean that residual indices predict body fat content better than ratio indices" (Labocha et al. 2014). Given this uncertainty, in this study we reported the results obtained using KFI as a ratio index to ensure comparability with other studies, although preliminary analyses suggested consistent findings when using KFI residuals.
None of the social and environmental variables that we have considered was directly related to temporal changes in fecundity. Spring-summer temperature was found to negatively affect adult deer fecundity only when the interaction with body condition was considered. In absence of randomised or experimentally controlled experiments, however, the use of multiple regression limits the possibility to explore the underlying cause-effect relationships among biological variables. Many studies on deer fecundity used multiple regression to investigate the directs effects of individual, social and environmental drivers on female reproductive traits, and the occurrence of indirect effects between, e.g., climate, density, body condition and fecundity often remained descriptive (but see Stopher et al. 2014). Nonetheless, the occurrence of indirect effects between different factors is assumed to explain variation in fecundity (Bonenfant et al. 2009)  knowledge this has never been quantitatively investigated. The observed variation in fecundity likely reflects an adaptive response of females, to optimize the trade-offs between the opportunity for reproduction and survival under different environmental conditions (cf. Sand 1996). It remains unclear, however, which mechanism might account for the observed relationships: given the predominantly nocturnal behaviour of red deer in our study site, we suggest body condition is most likely limited by summer resource availability (cf. Pettorelli et al. 2007), rather than by the temperature-mediated time constraints on foraging (cf. Mason et al. 2014). Ideally, to disclose which mechanism is at play, individual data on female summer foraging behaviour are needed. This information should be integrated with information on forage quality in the area occupied by females over the same period, for example using values of faecal crude proteins, as the use of NDVI in a population that inhabits both forested and open habitats within and between seasons might be problematic (Borowik et al. 2013), at least without marked individuals.
Path analysis represents an appealing approach to investigate cause-effect mechanisms in biology, yet there are limitations in the use of this methodology. Like other frequentist approaches, path analysis is concerned with finding a model that does not reject the null hypothesis (i.e. the hypothesized correlational structure is consistent with the correlational structure of the sample data). If a model is not rejected, however, we cannot be sure it is the 'true' model, as other models may fit the data equally well (Raykov and Marcoulides 2006). Alternative models may include direct or indirect effects of other variables that we did not take into account, and further research is needed to identify missed factors and disclose their mechanisms. In this respect, the association of model selection and path analysis may be useful, as the first allow to discriminate and identify influential variables, while the second may be used to test hypothesized causal relationships, helping to discriminate between potentially opposing mechanisms that can generate similar patterns.
Finally, it is worth noting that the 'direct' or 'indirect' effect of a variable in a path analysis, should be interpreted as 'relative to the other variables that are explicitly invoked in the causal explanation' (Shipley 2016), not with respect to any other variable that might exist (Shipley 2016 (Carlsson et al. 2018). In fact, our study population showed some evidence of negative consequences of Toxoplasma gondii infection on foetal development, and this effect changed with hind age (Formenti et al. 2015). Whether the inclusion of further parameters will enable a better prediction of deer fecundity, and whether increasing temperatures in the future years could have long-term effects on the life history of our study population, however, still remains to be investigated.
Notwithstanding the caveats in the application of path analysis, our study highlights the importance of considering more complex relationships between individual, social and environmental variables to explain variation in life history traits (cf. Stopher et al. 2014). In this respect, it appears crucial to formulate sound a priori hypotheses on which factors to include in a causal model, and on the direct and indirect relationships between them. Exploratory data analysis,