Direct and indirect effects of environmental drivers on reindeer reproduction

: The impact of climate change on the dynamics of populations has been well docu-mented and is widespread. However, weather variability influences populations both directly and indirectly, and is mediated by species interactions. This complexity may impede proper climate impact assessments. Hence, predicting the consequences of climate change may require including processes that occur both with time lags and across trophic levels. Based on our current understanding of the mechanisms linking local climate and trophic interactions in tundra ecosystems, we used a state-space formulation of a mediation model that allowed for assessing the relative contribution of direct and indirect environmental (weather and trophic) effects on reindeer Rangi fer tarandus reproductive success. Our study showed that the mediator effect of body condition caused delayed but predictable effects of weather, plant productivity, and reindeer densities on reproductive success. Furthermore, these predictors also affected reproductive success directly and with the same sign, suggesting that direct and indirect effects pulled in the same direction with respect to their combined total effect on reproductive success. Hence, poor weather conditions not only affect calf production negatively the same year, but also increase the likelihood of poor reproductive success the subsequent year. The results support the expectation that calf slaughter mass (as a proxy for herd body condition) is an important indicator of the state of reindeer herds with respect to their production potential and resilience to weather events and climate change. Finally, the model framework employed in the present study can be further developed as a potential vehicle for near-term forecasting, and thereby constitutes a useful tool for adaptive management.


INTRODUCTION
The impact of climate change on the dynamics of populations is widespread across ecosystems (Post et al. 2009, Hansen et al. 2013, Ims et al. 2019, Malhi et al. 2020. Climate change, as well as weather variability on shorter temporal scales, can affect populations both directly through physiological impacts and indirectly via trophic interactions in the food webs to which they belong (Wootton 1994, Abrams & Matsuda 1996. Indirect effects of climate change on population dynamics are also often mediated through direct effects on different life stages and life-history parameters, potentially promoting time lags in the response on population dynamics. For ungulates, body mass throughout the animals' lifetime is considered a key life-history parameter that ultimately affects survival and fecundity (Peters 1983 ABSTRACT: The impact of climate change on the dynamics of populations has been well documented and is widespread. However, weather variability influences populations both directly and indirectly, and is mediated by species interactions. This complexity may impede proper climate impact assessments. Hence, predicting the consequences of climate change may require including processes that occur both with time lags and across trophic levels. Based on our current understanding of the mechanisms linking local climate and trophic interactions in tundra ecosystems, we used a state-space formulation of a mediation model that allowed for assessing the relative contribution of direct and indirect environmental (weather and trophic) effects on reindeer Rangi fer tarandus reproductive success. Our study showed that the mediator effect of body condition caused delayed but predictable effects of weather, plant productivity, and reindeer densities on reproductive success. Furthermore, these predictors also affected reproductive success directly and with the same sign, suggesting that direct and indirect effects pulled in the same direction with respect to their combined total effect on reproductive success. Hence, poor weather conditions not only affect calf production negatively the same year, but also increase the likelihood of poor reproductive success the subsequent year. The results support the expectation that calf slaughter mass (as a proxy for herd body condition) is an important indicator of the state of reindeer herds with respect to their production potential and resilience to weather events and climate change. Finally, the model framework employed in the present study can be further developed as a potential vehicle for near-term forecasting, and thereby constitutes a useful tool for adaptive management.
KEY WORDS: Mediation · Arctic · Heterogeneous effects · Ungulates · Body condition · Reproductive success · SEM · Structural Equation Model 1992, Saether 1997, Gaillard et al. 2000. Many studies have revealed how these life-history traits depend on both current and previous environments (Lindström 1999, Metcalfe & Monaghan 2001, Harrison et al. 2011). Such indirect effects may have complex impacts on population dynamics (Mousseau & Fox 1998, Lummaa & Clutton-Brock 2002, Therefore, it is increasingly accepted that predicting the population dynamic consequences of environmental change in complex natural systems warrants an understanding of processes that occur both with time lags and across trophic levels (Bjørnstad & Grenfell 2001, Evans et al. 2013, Urban et al. 2016, Gellner et al. 2020. This is especially pertinent for many harvested species that reside in the middle of food webs, affected by trophic levels both above and below, as well as by processes within the same trophic level (Mellard et al. 2022 in this Special). However, so far, few studies have evaluated such food web effects (Ogilvie et al. 2017, Antiqueira et al. 2018. Large herbivores in tundra ecosystems are affected by several environmental drivers (Aanes et al. 2000, Tveraa et al. 2013, Hansen et al. 2019) including predation, food availability, and harvest rates (Helle & Kojola 2008, Hebblewhite & Merrill 2011, Tveraa et al. 2013, Serrouya et al. 2019. Moreover, the impact of these environmental drivers may vary spatially (Tveraa et al. 2007, Gaillard et al. 2013, Serrouya et al. 2019, and while some environmental drivers are known to affect demographic rates directly (Tveraa et al. 2013), others introduce time lags in the response. For instance, environmental effects on the body condition of female reindeer Rangifer tarandus in one year have carry-over effects onto next year's reproduction (Fauchald et al. 2004).
Based on our current understanding of the mechanisms linking local climate, plant productivity, and large herbivores and predators in the tundra (Tveraa et al. 2013(Tveraa et al. , 2014, we examined how measures of predator abundance, food availability, and weather directly and indirectly affect calf production in reindeer. Reindeer are the most widespread and ecologically/ economically important large herbivore in sub-Arctic and Arctic ecosystems. Here, we used the slaughter weights of reindeer calves in the autumn as an indicator of the body condition of reindeer in the herd in the autumn, and hence, a variable expected to be associated with subsequent reproductive success in this species (Fauchald et al. 2004, Tveraa et al. 2016. Moreover, autumn calf body mass has previously been found to be strongly influenced by environmental drivers (Tveraa et al. 2013). We made use of a rich data set containing 20 yr and 43 herding dis-tricts, covering ~60 000 km 2 , with data on abundance and calf production of semi-domestic reindeer in Northern Norway. We assessed the spatial variability in, and indirect effects, of environmental drivers through their impact on our measure of body condition in the herd the previous autumn. Our ap proach allowed us to assess the relative contribution of direct and indirect (i.e. time-lagged) effects on the total effect of a suite of environmental drivers. Fig. 1 provides a graphical depiction of the timeline of the hypothesized direct and indirect drivers of reindeer reproductive success that we aimed to estimate.

Semi-domesticated reindeer
Body mass is an important determinant of reproductive success in large ungulates (Saether 1997), and in reindeer in particular (Eloranta & Nieminen 1986, Lenvik & Aune 1988, Fauchald et al. 2004, Bårdsen et al. 2008. Previous studies have shown that body mass and reproductive success are affected by density (Bårdsen & Tveraa 2012), onset of spring greening and peak plant productivity (Bårdsen & Tveraa 2012), and snow depth in winter (Helle & Kojola 2008, Hendrichsen & Tyler 2014. Additionally, the slaughter mass of calves in the previous autumn is a proxy of animal body condition in the herd the previous autumn (see Section S1.2 and Fig. S2 in Supplement 1 at www.int-res.com/articles/ suppl/c086p179_supp1.pdf for validation of this measure), and affects reproductive success as measured by the number of calves counted in the summer relative to the number of female reindeer (Tveraa et al. 2014). Hereafter, we refer to calf slaughter mass in the autumn as 'body condition'; i.e. a variable representative for the body condition of reproducing females in a herd prior to reproduction.
From the national reindeer database (see Supplement 1, Section S1.1), for each district (Fig. 2) and year, we extracted the number of calves marked by the owners in summer or early autumn, the number of adult male and female reindeer, as well as calves recorded in the herd the previous winter, and the average slaughter weights of calves. Moreover, to assess the effect of density dependence on reproductive success and body condition, we calculated the total number of reindeer in the herd at the end of the previous winter as the sum of adult females and males as well as calves (i.e. ~1 yr old) in late winter. Because the area of the pastures varies across dis-tricts, we divided the total number of adults by area of the summer pastures to determine number of adults per km 2 .
Herds with summer pastures on islands ( Fig. 2) have previously been found to have heavier calves on average (Holand et al. 2010) than those with summer pastures on mainland Norway. Moreover, there is an expected spatial difference in climate according to a coast−inland gradient. To allow for such variation in the statistical analyses, we characterized the reindeer herding districts with respect to the oceaniccontinentality of their summer pastures using the classes 'island districts' (n = 10), 'island−continental districts' (n = 4) and 'continental districts' (n = 29) (i.e. 3 herding zones). Districts classified as 'island− continental districts' are borderline cases, as their summer pastures are situated at protruding peninsulas on the mainland (Fig. 2).

Snow depth
To represent conditions during late winter, a critical period during the winter for the productivity of the herds (Hendrichsen & Tyler 2014), we acquired data on daily snow depth (mm) from the Norwegian meteorological institute (MET Norway). The snow depth product is based on calculations from a collection of observational gridded data sets for temperature and precipitation that covers the Norwegian mainland (Lussana et al. 2016). The gridded data sets are based on the observations from the MET Norway's Climate Database, and the observations are interpolated on a regular grid (1 × 1 km; see Lussana et al. 2016). From this gridded data set, we calculated the average daily maximum snow depth (mm) in April in the winter pastures of the reindeer management regions (n = 6) per year as a measure of winter conditions during late winter. Hence, as we do not know the detailed space use of the districts in winter, except that they use winter pastures defined by spatial extent of their region, all districts within a region were given the same average value for snow depth.

Plant productivity and start of spring
To represent changes in peak plant productivity and phenology, we acquired MODIS vegetation indices from the US Geological Survey (Didan 2015). It has been shown that vegetation index-derived phenology agrees to a large extent with end-ofsnowmelt for the start of the growing season and 181 Fig. 1. Conceptual model showing the timeline of density dependence (DD) and environmental drivers and how they affect reproductive success of reindeer (i.e. calves per female), both directly and indirectly. The signs (+/-) denote the expected a priori effects of DD and environmental drivers on body condition and reproductive success, respectively. Note that indirect effects of DD and environmental drivers are from the year before (t − 1) and mediated through the effect of body condition prior to reproduction (t). Direct effects of environmental drivers happen the same year as the reindeer reproduce start-of-snowing for the end of the growing season (Jin et al. 2017). We chose the enhanced vegetation index (EVI; 16 d L3 Global 250 m), which is derived from atmospherically corrected reflectance in the red, near-infrared, and blue spectrum (Huete et al. 2002). We fitted a double logistic function to the seasonal EVI data from each pixel (Tveraa et al. 2013). Onset of spring was estimated using a parameter for the day in spring when average EVI was halfway between the seasonal minimum and maximum EVI, and peak plant productivity was estimated using a parameter for the average maximum EVI in summer. The pixel level parameter estimates were aggregated at the reindeer district level using the average values for all pixels within the summer pastures of the reindeer districts in Finnmark.

Lynx and wolverine family groups
In Norway, populations of lynx are monitored by recording the number of family groups (i.e. lynx females with kittens). Most of the recordings are conducted on fresh snow during winter. Recordings end in late February, because of lynx hunting and since lynx from different family groups start to congregate during the breeding season in late March, and therefore risk being wrongly assigned to family groups. All observations of tracks and individuals are verified by the Norwegian Nature Inspectorate (SNO) as belonging to a family group. Wolverines are monitored by recording the number of litters born every year. During late winter and spring, SNO rangers visit known dens with previously recorded pups as well as extensively search for new dens in areas with no previously known dens. We acquired data on lynx and wolverines from ROVBASE, a national database under the aegis of the Norwegian Environmental Agency, containing monitoring data on carnivores in Norway (https://www.rovbase.no/). Lynx and wolverines occupy large home ranges in boreal and Arctic ecosystems (Herfindal et al. 2005, Persson et al. 2010. In order to reduce the effect of small-scale variation in recordings of lynx and wolverines on their abundance, we used the sum of the number of family groups/reproductions recorded in the 6 management regions for both of these predators (Tveraa et al. 2014). We compiled the data for reindeer, lynx, and wolverines, along with environmental drivers, for the period 2000−2018. . While the traditional approach to mediation analysis is often based on the effect of explanatory variables in aggregate regression models that assume homogeneity of responses, such models may not be sufficient to assess the presence of mediation for different sites or populations. Hence, multisite studies open new opportunities in such analyses of food webs, as impact parameters associated with environmental drivers can be allowed to vary across sites with respect to both direct and indirect effects (Raudenbush & Schwartz 2020). To assess the effect of different covariates on reindeer reproductive success, we used a statespace formulation of a mediator model that allowed for the mediating effects of body condition on reproductive success as well as both measurement and process error (see Section S1.3 in Supplement 1 for JAGS model code). Moreover, it allowed for the modeling of heterogeneous effects of predictors. The process model consists of 2 sub-models; the first, where the mediator, body condition, is modeled as a function of different temporal covariates at t − 1, and the second, where reproductive success is a function of temporal covariates at t, including body condition from the first sub-model (Fig. 1). In the model, all covariates are scaled (over all districts and time points), with mean = 0 and variance = 1, to ease convergence and interpretation of relative effect sizes, and for the calculation of indirect and total effects (Enders & Tofighi 2007, van de Pol et al. 2011. We modeled all the covariates as average temporal effects (X t ), i.e. the average value of a covariate across sites per year. This was done to minimize the spatial noise in some of the predictors (e.g. snow amount, reindeer density), with the goal of improving estimation of their effect. Hence, predictor effects, including heterogeneous effects, are based on the range of the temporal predictor variables. Note that a model with the full predictor effect, using data for both districts and years, provided model coefficients consistent with the model using only average temporal effects.

Statistical analyses
In the model likelihood, the number of calves was assumed to follow a lognormal distribution, i.e. the logarithm of the number of calves has a Gaussian distribution with mean, log(ρ s,t ), and variance, σ 2 m.err : Further, we modeled ρ s,t as the product of the number of recorded females and a latent reproductive success rate, ƒ, in each district and year: In the body condition sub-model, we assumed scaled body condition, bc s,t , to be normally distributed with mean, μbc s,t , and variance, σ 2 bc.p.err : bc s,t~n orm(μbc s,t , σ 2 bc.p.err ) where σ 2 bc.p.err is the variance of the process error for bc. Then, μbc s,t was modeled as a function of temporal covariates: where b.bc type is the herding zone fixed effect denoting the location of districts in 3 different climatic areas (continental, island and island−continental zones), rD bc is a district random effect [~Norm(0, σ 2 rD.bc )] and β.bc x is a vector of coefficients associated with temporal covariates X.bc t . We included the following covariates in the body condition sub-model: average maximum snow in April the year before (t − 1), onset of summer the year before (t − 1), peak plant productivity the year before (t − 1) and the number of adult reindeer per km 2 the year before (t − 1) as a proxy for density dependence (DD) (Fig. 1). Preliminary analyses and plotting of the data showed that the relationship between the number of adults (i.e. DD) and body condition varied strongly across districts. Hence, we modeled the parameter for DD with a random slope (i.e. β.bc.DD s × NoAdult st-1 , where β.bc.DD s ~Norm(mu β , σ 2 β )) and used hyperpriors (i.e. a probability distribution on a hyperparameter, where hyperparameter is defined as a parameter of a prior distribution) for the coefficients, to borrow strength across districts for the district-level estimates.
In the reproductive success sub-model, we modeled latent reproductive success, ƒ s,t , by constraining parameters using regression on the logit-link scale and we assumed the logit of reproductive success to be normally distributed with mean, μƒ s,t , and variance, σ 2 ƒ.ρ.err : logit(ƒ s,t )~norm(μƒ s,t , σ 2 ƒ.ρ.err ) where σ 2 ƒ.ρ.err is the variance of the process error for ƒ. We then modeled μƒ s,t as a function of temporal covariates, including body condition, μbc s=1:Ndist,t , i.e. the response from the body condition sub-model: where b.ƒ type is the herding zone fixed effect (the 3 different climatic areas; continental, island and island−continental), rD ƒ is a district random effect [~Norm(0, σ 2 rD.f )] and β.ƒ x is a vector of coefficients associated with covariates Xƒ t . We included the following covariates in the reproductive success submodel: average maximum snow in April the same year (t), onset of summer the same year (t), peak plant productivity the same year (t), number of adult reindeer per km 2 (i.e. DD) the same year (t), as well as the number of family groups or reproductions of lynx and wolverines, respectively ( Fig. 1). Preliminary analyses and plotting of the data showed that the relationship between body condition and reproductive success (i.e. calves/ females) varied strongly across districts. Hence, we modeled the parameter for body condition with a random slope (i.e. β.ƒ.bc s × μbc s=1:Ndist,t ) and used hyperpriors for the coefficients, to borrow strength across districts for the district-level estimates.
Our analysis was performed using JAGS (Plummer 2003), which uses Markov chain Monte Carlo simulations to estimate posterior probability distributions. All effect sizes from the analysis are given by the mean of the posterior distribution and the 95% credible interval (CI), if not otherwise stated. Direct, indirect, and total effects of the different covariates were calculated from the posterior parameter distributions and are given as the mean and the 2.5, 50, and 97.5 th percentiles (i.e. the mean, median, and the 95% CI).

Calculation of indirect and total effects
Given the heterogeneous effects for DD in the body condition (bc) sub-model and for body condition in the reproductive success sub-model, the product of average effects can be used for estimating average effects if the coefficients are uncorrelated, i.e. if the covariance between βDD t−1 [s] and βbc t [s] is zero, then by definition of the covariance: As the correlation between βDD t−1 and βbc t was very low (ρ = −0.14), we focused on the average, acrossdistrict effects for the presentation of indirect and total effects (but see Table S3 in Supplement 2 at www.intres.com/articles/suppl/c086p179_supp2.xlsx for results concerning district-specific direct, in direct, and total effects). We calculated indirect and total effects using the posterior distributions. For the non-heterogeneous predictors, indirect effects were calculated as: where βX denotes one of the predictors summer onset, plant productivity, or winter snow amount and where (βbc t [iter,]) denotes the mean across districts for each iteration (iter) in the posterior chains.
The heterogeneous indirect effect of DD t−1 was calculated as:
Most of the environmental predictors significantly influenced body condition. Increased amount of snow in April had a small negative effect, later onset of spring had a strong negative effect, and there was a tendency for a positive effect of increased peak plant productivity on body condition (Fig. 4). Moreover, density dependence had on average a strong negative effect on body condition (Fig. 4, Table S1 in Supplement 2).
Most of the environmental predictors also influenced reindeer reproductive success (Fig. 4). Increased amount of snow in April had a negative effect, later onset of spring had a negative effect, increased peak plant productivity tended to have a positive effect, and there was evidence for a strong effect of density dependence on reindeer reproductive success. Body condition had on average a strong positive effect on subsequent reproductive success (Fig. 4, Table S1). Increased number of wolverines had a negative effect, and an increase in lynx abundance had a weak negative effect.
The strong positive effect of the mediator, body condition, on subsequent reproductive success implies that most environmental covariates had both an indirect and a direct effect on reproductive success (Fig. 4). For each of the environmental variables, the estimates of direct and indirect effects all had the same sign (Fig. 4), suggesting that they pulled in the same direction with respect to their effects on reproductive success.
The combined indirect and direct effects resulted in strong negative effects of snow amount in April, spring onset, and density dependence on reindeer reproductive success, while the total effect of peak plant productivity was slightly weaker (Fig. 4, see Tables S1 & S2 in Supplement 2 for  more details and Table S3 for effects on the level of districts). There was, however, variation among the different drivers concerning the relative contribution of direct and indirect effects to the total effect on reproductive success. The effect of winter snow amount was mostly through its direct effect (~87%), while the relative contribution of direct (~54%) and indirect (~46%) effects of spring onset was similar. The effect of peak plant production on reproductive success was slightly dominated by its direct effect (~64%), while the effect of density dependence was slightly dominated by its indirect effect (~60%).
There was spatial variation among districts in the effect of density dependence on body condition (−1.45 ± 1.08, [−5.18; 0.76]), though most districts showed negative effects (Fig. 5a). There was also evidence for spatial variation in the effect of body condition on reproductive success (0.54 ± 0.18, [0.11; 0.80]), however, these effects were consistently positive (Fig. 5b).  Table S1 in Supplement 2 for the districts assigned to each of the 6 regions

DISCUSSION
If we are to understand and predict how current and future climate changes affect species in natural ecosystems, we need to understand how local climate (i.e. weather) influences populations both directly and indirectly (Mellard et al. 2022, Box 3). This information pertains to both species interactions and withinspecies processes that influence important future life-history events. Our mediation analysis showed strong and consistent direct effects of density de pendence, body condition, the amount of late-winter snow, and the onset of spring. Despite strong control of the major predators of reindeer in these areas, we found a negative effect of the number of wolverines, but not lynx. Of the indirect effects, the strongest impact was negative density dependence, followed by effects of the onset of spring and the amount of snow in late winter. Importantly, the indirect effects contributed substantially to the total effect of the environmental drivers on reproductive success.

Local climate, predators, and reproductive success
In recent decades, global climate change has caused a shift towards earlier spring and an increase in plant productivity in tundra, changes expected to increase further in the coming decades (Elmendorf et al. 2012, Ims et al. 2013. Our results suggest that reindeer in this region could experience improved reproductive success due to higher quantity and quality of forage associated with increasingly earlier spring and higher plant productivity. However, while plant productivity has increased and spring arrives earlier on average (Karlsen et al. 2009, Park et al. 2016), the average reproductive success of reindeer has decreased over the last 2 decades in this part of Fennoscandia. This trend suggests that strong negative density dependence combined with a negative effect of deteriorating late-winter conditions (i.e. late-winter snow conditions) likely both directly and indirectly counteract the positive effect of increased food availability (cf. Tews et al. 2007, Tveraa et al. 2013. Moreover, although plant biomass tends to increase in a warming climate, the abundance of palatable plant functional groups, such as herbs and deciduous shrubs, may have become depleted over time because of a consistently high grazing/ browsing pressure (Bråthen et al. 2007(Bråthen et al. , 2017. The degree to which predation by large carnivores affects free-ranging reindeer has been widely debated over the last decades throughout Fennoscandia. This topic has been an especially prominent part of the conflict concerning claims for compensation due to losses in the reindeer industry (Tveraa et al. 2014). Wolverines (Gustine et al. 2006, Mattisson et al. 2016) and lynx (Mattisson et al. 2011) are known to be predators of reindeer/ caribou. Hobbs et al. (2012) found that both wolverines and lynx reduced harvest in the Swedish reindeer industry; similar results were also found by Tveraa et al. (2014). In both the Swedish and Norwegian study, effect sizes were higher for lynx than wolverines. In contrast, in our study focusing on the tundra biome, we found the estimated effect size of wolverines to be higher than the effect size of lynx. This is not surprising, as lynx Fig. 5. Distribution of the heterogeneous effect estimates (i.e. across districts) of (a) effects of density dependence on reindeer body condition and (b) effects of body condition (i.e. measured as calf slaughter weights) on reproductive success, respectively mainly inhabit forested biomes, while the range of the wolverine includes both boreal forest and Arctic tundra.

Heterogeneous direct and indirect effects
The traditional approach to mediation analysis is based on the effect of explanatory variables in aggregate regression models that assume homogeneity of responses. Here, we used a Bayesian model that incorporated heterogeneity in the mediating behavior as well as in direct effects across districts. We found a strong average negative indirect effect of density dependence and a positive direct effect of body condition on reproductive success. However, there was substantial spatial variation in the strength of these effects (Fig. 5). This variation was especially prominent for the indirect effect of density dependence, with many districts showing weak or no effect of density dependence (Table S1). This result seems surprising given the high densities of reindeer in most parts of Finnmark. However, measurement error in the counts of reindeer could result in attenuation of the estimated effects of both density dependence and body condition (Lebreton & Gimenez 2013). Given that measurement error is unlikely to be a homogeneous process (i.e. varying among districts and across time) and involves both bias and precision to un known degrees, we could not model it in our analyses. Hence, the contribution of measurement error, as well as spatial variation in population fluctuations, need to be further explored.
Our approach combined a hierarchical model (districts within regions) with a mediation model to identify heterogeneous direct and indirect effects. Such an approach has become commonplace in social sciences and epidemiology (e.g. Bloom et al. 2017, Raudenbush & Schwartz 2020, where largescale randomized experiments and observational studies have helped to identify heterogeneous treatment effects. The hierarchical structure of the model, with a common distribution for districts within a region, presents the advantage of borrowing strength when estimating individual district effects, but also assumes that districts are exchangeable within a region, i.e. that districts have similar statistical properties be sides the heterogeneous effects. While this assumption is often made, its validity is seldom assessed. Epidemiologists have shown that this can lead to biases in estimated effects if the assumption is not fulfilled (Greenland & Robins 1986, Robins & Greenland 1992, and should therefore be considered in future studies using these methods in ecology.

Implications for reindeer herding
Stochastic effects of weather are usually considered to operate on ungulate populations without delays (e.g. Aanes et al. 2000). Our study showed that mediator effects, operating through body condition, cause delayed and predictable effects of weather on reproductive success and, in turn, the number of calves that will be available for harvesting. Poor weather conditions not only affect calf production negatively the same year, but also increase the likelihood of poor reproductive success in the subsequent year. Hence, body condition, as measured by calf slaughter mass, appears to be an important indicator of the state of reindeer herds-both with respect to their production potential and their resilience to weather events and climate change. We also note that the effects of reindeer density on body condition varied markedly among districts. Identifying the mechanisms behind this variation will likely be important as it may indicate whether changes in reindeer densities will affect herd productivity.

Perspectives on improved scope for prediction
The model framework employed in the present study can be further developed for use in near-term forecasting. It may thereby become a useful tool for adaptive management of reindeer in this region (Pouyat et al. 2010, Hobbs et al. 2015, Hobday et al. 2016, Dietze et al. 2018, Nichols et al. 2019. Indeed, predicting responses to environmental disturbances and change is a long-standing goal in ecology (Mouquet et al. 2015, Petchey et al. 2015 and is particularly relevant for adaptive management of harvested populations (Lande et al. 2003). It has been argued that a deeper understanding of the multitude of ways environmental drivers affect food webs will im prove our ability to predict the future states of populations (Wootton & Emmerson 2005, Montoya et al. 2009). Some studies also suggest that simple models might be adequate for making predictions regarding population states under equilibrium conditions, while predictions under transient conditions demand more complex models (Stephens et al. 2002). However, the benefit of including indirect effects or specific mechanistic understanding is not obvious and largely untested (Farrer et al. 2014, Bircher et al. 2015, Gerber & Kendall 2018, Yates et al. 2018, Henden et al. 2020. Hence, we believe predictions of the future state of reindeer reproductive success in Finnmark will benefit from a multi-model approach, where predictive ability is assessed from several models that vary in the spatial and temporal complexity of drivers (cf. Henden et al. 2020) and complexity of the model structure.