Effects of weather and hunting on wild reindeer population dynamics in Hardangervidda National Park

Funding information Horizon 2020 research and innovation programme, Grant/Award Number: 641762 Abstract Wild reindeer have a range that extends across the circumpolar region. In the last few decades, however, populations of wild reindeer have been on the decline. The reasons for these declines are poorly understood, but are suggested to be linked to both local and global climatic factors, disease, and human interference. Hardangervidda plateau in Norway is home to the largest wild reindeer population in Europe, and is at the southern end of its European range. This population is therefore of particular importance, particularly in the light of climate change. We investigated howweather and hunting have affected the wild reindeer population in Hardangervidda over the last two decades. Our findings suggest that the wild reindeer population in Hardangervidda is most affected by winter temperature and hunting, where colder temperatures and lower harvest rates typically result in higher growth rates. We did not find significant evidence for linear density dependence. Our results show trends across Hardangervidda, and give an indication of how region-wide weather and hunting pressure can affect the wild reindeer population. As new data emerge, future investigations should look into the existence and nature of density dependence and the influence of other weather and human disturbance related factors.

Hardangervidda plateau in Norway is home to the largest wild reindeer population in Europe, and is at the southern end of its European range. This population is therefore of particular importance, particularly in the light of climate change.
We investigated how weather and hunting have affected the wild reindeer population in Hardangervidda over the last two decades. Our findings suggest that the wild reindeer population in Hardangervidda is most affected by winter temperature and hunting, where colder temperatures and lower harvest rates typically result in higher growth rates. We did not find significant evidence for linear density dependence. Our results show trends across Hardangervidda, and give an indication of how region-wide weather and hunting pressure can affect the wild reindeer population. As new data emerge, future investigations should look into the existence and nature of density dependence and the influence of other weather and human disturbance related factors.

K E Y W O R D S
caribou, GPS collar, probabilistic forecast, Rangifer tarandus, Ricker model

| INTRODUCTION
Wild reindeer, or caribou (Rangifer tarandus), are keystone species of the circumpolar region, because they influence ecosystem processes such as nutrient cycling and primary production (e.g., Olofsson, Stark, & Oksanen, 2004). While they have a large range that extends across northern Europe, Siberia and North America, the species was recently categorized as vulnerable by the International Union for Conservation of Nature (IUCN), due to a 40% decline over the last two to three decades (Gunn, 2016). The exact extent and reasons for the decline are not well understood and the causes vary across the distribution of the species, but are suggested to be related to landscape change and fragmentation, unregulated hunting and a warming climate (Gunn, 2016). Norway has also seen a contraction in the range size of reindeer, and human infrastructure, diseases, hunting and climate change, are possible factors underlying this pattern. The recent discovery of Chronic Wasting Disease (CWD) in a Norwegian wild reindeer population in 2016 is also a cause for great concern for the stock of European wild reindeer (Hansen et al., 2017).
Reindeer populations are known to be affected by a variety of climatic factors. For example, high population growth rates have been linked to dry winters, although climate effects are probably more important at high population densities (Aanes, Saether, & Øritsland, 2000). A large scale climatic pattern, the North Atlantic Oscillation, has also been shown to affect body mass (Post & Stenseth, 1999) and reindeer population fluctuations (Forchhammer, Post, Stenseth, & Boertmann, 2002;Uboni et al., 2016). However, reindeer are clearly also affected by local weather patterns. As reindeer are dependent on a relatively short summer season to gain body weight to last them through the winter, any factors that may reduce their ability to feed during this time will inevitably have a negative effect on their survival (Reimers, 1997). Insect harassment is one factor that is particularly detrimental to optimal feeding during the summer, and it has been shown in a Norwegian population of reindeer that hotter summers with higher levels of insect harassment are correlated with lower body mass (Colman, 2000). In addition to insect harassment, the timing of greening is also very important, particularly in determining calf survival. Wild reindeer are shown to have a calving time that is highly synchronized with the onset of greening (Post, Bøving, Pedersen, & MacArthur, 2003;Skogland, 1989). However, climate change is also expected to shift the timing and the duration of greening, particularly at high latitudes (Molau, 1997;Walther et al., 2002). Migratory herbivores are believed to be able to compensate for this temporal advancement by tracking phenological variation across landscapes, however, climate warming can reduce spatial variability in plant phenology, consequently causing a decline in calf production (Post, Pedersen, Wilmers, & Forchhammer, 2008). Indeed, trophic mismatch has already been documented in reindeer populations, where the timing of calving has not kept up with the advancement of greening as spring temperatures have risen over time .
One climatic factor that may control reindeer populations is snow cover and quality. Snow cover strongly affects lichen biomass and lichen heath development (Odland, Sandvik, Bjerketvedt, & Myrvold, 2014;Skogland, 1978), which reindeer greatly depend upon as a food source in the relatively long winter period (Heggberget, Gaare, & Ball, 2002). Lichen cover is also known to be reduced by high precipitation in Hardangervidda (Odland et al., 2014). In addition to affecting the availability of forage, snow cover also plays a role in the accessibility of lichen heaths. When snow thaws and re-freezes causing layers of ice in the snow, digging for food becomes particularly problematic. Thus, snow depth is not the single most important factor in determining accessibility to winter forage, but rather the combination of both the depth and the consistency of the snow. As a result, the locations of winter grazing pastures are often much more widely distributed than those of the summer pastures, because access to winter grazing pastures is dependent on local snow cover conditions (Iuell & Strand, 2007). It has also been suggested that the suitability of winter pastures determines the effect that hunting has on population regulation, where reindeer with access to good winter conditions are regulated by hunting, and those with access to poor conditions are regulated by bottom-up processes (Tveraa et al., 2007). Thus, although hunting and other human impacts most definitely contribute to reindeer population fluctuations, it is clear that anthropogenic and climatic factors frequently interact.
The reindeer population on the Hardangervidda plateau is the largest wild reindeer population in Europe, and is therefore important for its ecological value, but also for its economical and recreational value for hunters and landowners (Bjerketvedt, Reimers, Parker, & Borgstrøm, 2014). Thus, a loss of this population would have a negative impact not only for the ecosystem, but also for the people that depend on these animals for their livelihood. The reindeer population on the Hardangervidda plateau is also at the southern end of its European range. Studying the changes in environmental attributes of this region, particularly those which are expected to change the most under climate change, are therefore of particular interest when considering the future of this population. Although previous studies have addressed the effect of weather on wild reindeer populations, it is important to consider local populations, as there can be large variations in response to local weather patterns. Given the importance of this population in a European context, it is therefore of particular importance to gain detailed knowledge at a local level. The aim of this paper is therefore to investigate the effect of different weather variables, and how these affect reindeer population fluctuations on the Hardangervidda plateau.

| Study site
The wild reindeer population studied in this paper inhabits the Hardangervidda plateau in southern central Norway (60 6 0 N, 7 28 0 E). Due to its area of over 8,000 km 2 , there is a large variation in climate, particularly along an east-west gradient, where the west tends to have an oceanic climate and the east is more continental. In general, the climate is characterized by alpine conditions with an average annual temperature of 1 C, and precipitation ranging between 700 and more than 2,000 mm per annum (www.met.no). The average elevation of the plateau is around 1,200 m.a.s.l., but with peaks ranging between 900 and 1,800 m.a.s.l. The west generally has more rugged terrain, while the eastern and central areas are more flat. The plateau is characterized by numerous rivers and lakes, with an open landscape largely above the tree line, and dominated by middle alpine vegetation (Tveitnes, 1980).

| Population data
The reindeer population in Hardangervidda, southern central Norway, has been monitored for the last six decades, with regular population counts taken to estimate the population. However, despite the fact that reindeer form very large herds, these counts have varied in their efficacy and sampling effort, and thus, population estimation remains difficult (Strand, Gaare, Solberg, & Wilmann, 2004). This is largely due to the large range that the animals use on the Hardangervidda plateau. In recent years, however, population counts have been vastly improved with the use of helicopters and the implementation of GPS collars in 2002. Population fluctuations can be traced back to hunting over the last decades and these most likely control the Hardangervidda population, but it is also clear that this is only the case if there is sufficient forage (Strand et al., 2004). Hardangervidda has a clear climatic east-west gradient, where precipitation is considerably higher in the west. As a result of this climate gradient, winter grazing grounds tend to be in the eastern parts of the national park, although annual variation in snow quantity and quality determines how much of the most central parts of the region can be used by the animals.
The Wild Reindeer Centre provided population count and hunting data from 1994 to 2015. Two counts are performed each year. During both counts, reindeer are counted in classes; calves (0-0.5 years old), cows and bulls. One count occurs in February-March. This population count is partial (i.e., only a portion of the population is actually counted) and aims to estimate the population structure of the winter population. During this count, 2-5 small airplanes are flown systematically over Hardangervidda with 3 km spacing between paths. GPS-tagged individuals are not used during the count itself, but are used after the count to check if groups of animals with tagged individuals have been found, that is, to check for sightability. The success of these counts is mostly dependent on snow, cloud and light conditions and therefore counts were only performed in good weather conditions.
Another count is performed in July to estimate the calf/cow ratio, where only female groups are actively searched for, and bulls are probably underestimated.
Apart from the very first years, monitoring was performed from a helicopter, and all groups that were spotted were photographed and counted. Flight times of 4-6 hours are standard, although this can vary from year to year. The Hardangervidda population of reindeer is consistently concentrated in the southwest of the park in early spring and summer, so animals are searched for systematically in areas where most cow/calf groups are known to be found. Beyond this, flights are also made into other regions of the national park, based on information from reported sightings, previous experience from other counts, and GPS-tagged individuals. In general, the number of GPS-tagged individuals in a given year can affect the count if these individuals are used to find groups of animals. However, the number of tagged individuals with an updated position directly before the count has not varied much between years. The success of these counts can be dependent to some extent on light conditions, group size, and weather conditions. However, the success of the counts is considered to be comparable from year to year (S. E. Lund, personal communication, January 31, 2019).
In Hardangervidda, in each year, hunting licenses are issued, and when an animal has been killed, this is reported back to the park authority. This reporting is believed to yield a reliable estimate of the number killed through hunting in practice. The hunting season lasts from the 20th of August to the 30th of September, which is after the annual population counts take place. Intuitively, the number of animals killed in the autumn of year i would be expected to heavily affect the population change between year i and year i + 1 (in fact, the park authority uses hunting quotas as a method of controlling the population). Yearly hunting quotas are also available for the whole study period, with details on sex and age classes. The hunting regime of this population results in approximately equal numbers of males and females being killed. From 1995 to 2015 an average of 41% of the individuals killed were female, 39% were male, and the rest were calves (11% male and 9% female). In terms of age classes, the majority of individuals killed between 1995 and 2015 were above 2.5 years of age (62%, where 35% were female and 27% were male). The 1.5-year age class makes up 18% of the killed individuals (6% female, 12% male). Although the number of large mammal predators in southern Norway is growing (lynx, wolverine, wolves and brown bear), the number of individuals killed by predation is negligible as predators have been eradicated and are strictly controlled on the Hardangervidda plateau.
The winter count aims to estimate population structure rather than size, however, a good fraction of the actual population is counted, and we assume that the population structure of counted animals is representative of the whole population. Since July counts are reliable only for females and calves, we added the rate of bulls as calculated during the successive population structure counts to estimate yearly population abundance. We cannot assume that male and female culling is balanced; therefore, we also took into account the removal of bulls during the hunting season. We estimated the number of bulls in the summer as follows: Bulls_ jul = bulls_winter=cows_winter ð Þ * cows_ jul ð + cows_cullÞ -bulls_cull: where "*_ jul" are bulls and cows counted in July, "*_winter" are those counted in the winter, and "*_cull" is the number of killed animals. In 1998 and 2014 population structure counts were not performed due to bad weather conditions. The bull to cow ratio in winter is marginally serially correlated (r = .57), therefore, we linearly interpolated the ratio for the missing years, when population structure counts were not performed, thus, obtaining a complete time-series of estimated population abundance ( Figure 1).

| Weather data
The Norwegian Meteorological Institute provided gridded datasets (1 × 1 km resolution) of mean daily temperature and total accumulated precipitation (version 1.1). This dataset is produced by interpolating weather data from climate stations, and using a 100 m resolution Digital Elevation Model (for more information see: The Norwegian Meteorological Institute, 2010; Tveito, Bjørdal, Skjelvåg, & Aune, 2005;Tveito & Førland, 1999). All weather data were clipped to the extent of the Hardangervidda reindeer area (e.g., Falldorf, Strand, Panzacchi, & Tømmervik, 2014) and summarized per year and season, depending on the variable in question.
Snow cover is a potential driver of changes in the reindeer population as it affects the ability of the animals to be able to reach lichen on the ground. However, snow cover is not necessarily a problem in itself, as reindeer are able to dig through at least 150 cm of snow to get to their food, if the snow has not formed a crust or become particularly packed (O. Søbekkseter, reindeer herder, personal communication, March 2017). Thus, melting and refreezing (icing) is an important factor to consider. Given the data available, the mechanism in which the snow melts and refreezes is difficult to quantify. However, when the mean temperature stays far below zero, icing is less likely to occur. In cold temperatures, precipitation is also much less likely to fall as rain, which can freeze to cause the same icing effect as melting snow. We have therefore chosen the mean temperature in the coldest months (January and February), the mean temperature in the late winter/early spring (March and April), as well as the number of days in which the temperature exceeds 0 C in February and March as possible variables affecting the population.
We have chosen the average summer temperature in the warmest months (July and August) as a potential predictor variable for relative population change because insect harassment, which can negatively affect bodyweight and mortality, tends to be greater in warmer temperatures (Colman, 2000). Furthermore, warm temperatures in themselves can be detrimental for reindeer as they are adapted to cold temperatures. Good quality and quantity of summer forage allows the animals to be suitably healthy going into the winter season, and thus may be expected to decrease winter mortality. The number of growing degreedays (the sum of the daily difference between the average F I G U R E 1 Estimated reindeer counts in Hardangervidda between 1994 and 2015. The winter population count (blue line) shows the partial count that occurs in the winter, aimed at estimating the population structure. The summer count (red line) refers to the count done in July, which estimates the calf/cow ratio, where female groups are actively searched for and bulls are probably underestimated. These counts are used to obtain the total population estimate (black line) [Color figure can be viewed at wileyonlinelibrary.com] temperature and 5 C) during the summer may therefore be expected to have some effect on the population. However, the number of growing degree days and the mean summer temperature are strongly correlated (r = .73, p = .003), and thus, if high summer temperatures genuinely have a negative effect on the population through insect harassment, this effect may be partially canceled out if the number of growing degree days has a positive effect through improved quantity and quality of summer forage. We have therefore chosen to include growing degree-days (from June to September) as a predictor variable as well.

| Population models and forecasting
In this article, the effects of density dependence, hunting and various weather conditions on the reindeer population of Hardangervidda are investigated. First, the effect of each potential predictor variable on the population was assessed and a formal statistical test was applied (see section on Testing for Density Dependence and Other Drivers of Population Change). Second, a probabilistic forecasting framework was defined in which the following year's population is predicted. While the former approach allows for the testing of single predictor variables, the latter allows for the testing of combinations of variables and interactions simultaneously.
where n i is population size in the ith year, a, b and σ are parameters to be selected and ε i is a random draw from a standard Gaussian distribution. If the parameter b makes a significant contribution to forecast performance, the population is said to be density dependent. The Stochastic Ricker Model can be rewritten in the form where log(n i + 1 /n i ) is referred to as the relative population change. The Stochastic Ricker Model is thus a generalized linear model with the relative population change as the dependent variable and the current population as a predictor variable.
The Stochastic Ricker Model can be extended to include variables beyond just the current population size. The Modified Stochastic Ricker (MSR) Model (Jacobson, Provenzale, von Hardenberg, Bassano, & Festa-Bianchet, 2004) is defined by where v i,j is the value of the jth variable in year i and a, b 1 ,…,b J and σ are parameters to be selected. The MSR model allows any suitable variable to be a linear predictor variable of the relative population change.

| Testing for density dependence and other drivers of population change
A number of different variables were considered as potential drivers of changes in the population. Commonly, large population sizes can cause increased mortality in animal populations because of increased competition for scarce resources. It is thus generally considered that a habitat imposes a limit on the size of the population. To assess whether there is evidence of this effect in the reindeer population, the current population size was plotted against the relative population change to look for obvious (perhaps nonlinear) patterns. In addition to density dependence, it is also desirable to investigate how other variables relate to the population growth rate. We therefore considered the relationship between the relative population change and both the harvest rate (number of animals reported killed in year i as a proportion of the recorded population count in that year) and the five weather variables chosen as potential predictors. These are: mean temperature in January and February, mean temperature in March and April, mean temperature in July and August, the number of days in which the temperature exceeded zero degrees Celsius in February and March, and the number of growing degree-days. It should be noted that we found no temporal trend either in the population estimates (Mann-Kendall test, τ = −0.032, p = .871) or in the relative population change (Mann-Kendall test, τ = 0.074, p = .673) and therefore the risk of spurious correlations caused by trends in the variables (Royama, 1981) can be considered negligible. Cook's distance (Cook, 1977) was calculated in each case to test for influential observations. An observation was considered "influential" if its Cook's Distance exceeded one. An influential observation is found corresponding to the year 2010 and therefore the analysis was performed with and without this year included.
The linear relationship between each of the above variables and the relative population change was also assessed by applying Pollard's randomization test (Pollard, Lakhani, & Rothery, 1987). This test (see Supporting Information 1 for details), was first designed to detect density dependence, that is, to test whether the size of the current population linearly affects the relative population change, but is easily extended to test the effects of other variables as well. Potential lagged effects for weather variables were also tested. It is worth noting that the test only considers linear relationships and, by design, will not detect some nonlinear relationships. In all cases, significance was determined using the standard significance levels in statistics. A significant result at the 10, 5 and 1% levels are defined to be "weakly significant," "significant" and "strongly significant," respectively. Due to concerns about finding significant results by chance through testing too many variables, adjusted p-values were calculated using the Westfall-Young permutation procedure (Westfall, Young, & Wright, 1993), which adjusts the p-values, while accounting for correlation between the tested variables.

| Probabilistic forecasting
A probabilistic forecasting framework was defined to predict population counts up to 1 year ahead given the previous count and other predictor variables. The population models are stochastic and therefore naturally produce a probabilistic forecast distribution of the relative population change. These, in turn, were used to form forecasts of the actual population sizes using Monte-Carlo simulation with a large number of points (10,000). The simulations were then converted back to probabilistic forecast distributions using kernel density estimation with a kernel width chosen using Silverman's rule of thumb.
A number of combinations of predictor variables were chosen as candidate models, with the limitation of using at most three variables in each model, to avoid overfitting. In the case of weather variables, only those found to be significant at the 10% level (i.e., weakly significant) according to the Westfall-Young permutation procedure were considered as potential predictor variables. As for first-order interactions, we chose to include those between the current population size (density dependence) or the proportion of animals killed and the weather variable and that between the population size and the proportion of animals killed. A two-sided t-test for Pearson's correlation coefficient between the proportion of animals killed and the population size was conducted and not found to be significant (r = .20, p = .39). Models including both the proportion of animals killed and density dependence were therefore considered.
To measure the performance of each set of variables as a predictor of the population, each combination was used as an input to the Modified Stochastic Ricker model and the skill of the forecasts was assessed probabilistically. This, however, only gives a measure of the relative skill, that is, is only useful as a measure of the relative performance of each combination of variables as an input to the MSR Model. As a result, each set of forecasts says little about the absolute value of the forecasts, that is, whether the models have any value at all. In this article, we define a number of benchmark models with which to compare the value of the model-based forecasts. The benchmark models considered were: • Null Model-Stochastic Ricker/Gompertz model with no explanatory variables. • Climatology-A probability distribution of pastobserved populations. • Persistence-Probabilistic forecasts centred around the count from the previous year.
The benchmark model with the highest support according to the model selection statistics was chosen and defined as model 0. More details of each benchmark model are given below.

| Climatology
The climatology is defined as a distribution derived purely from past-observed population estimates. If no useful information (including the last population size) exists regarding the next estimated population size, the climatology can be considered the best possible probabilistic forecast, as long as it is a robust estimate of the underlying distribution of possible population sizes. As a result, if a model-based probabilistic forecast does not outperform the climatology on average, it is of little value since a better score could be achieved by issuing the climatology as the forecast.
The estimated climatology is shown in Supporting Information 2 (see Supporting Information 1 for details of how this was calculated). Although the estimated climatology is bimodal, it is likely that such a pattern occurs because of the small correlated sample rather than an underlying physical phenomenon.

| Persistence forecasts
An alternative benchmark model can be formed on the basis that the population of reindeer will be largely similar to that of the previous year. Although the population dynamics are usually expected to be more complex than this, the persistence provides an alternative benchmark for any model-based forecasts. Since forecasts are probabilistic in nature, some methodology was therefore required to convert single population estimates into predictive probability distributions. Forecasts of the population in the ith year were defined to take the form where σ represents the SD of the forecasts. The value of σ was chosen using leave-one-out cross-validation with maximum likelihood. Since a population cannot be negative, the forecasts were renormalized such that the integral over (0, ∞) is equal to one.

| Model selection
To find the most informative model in terms of predicting the reindeer population, some formal model selection technique was required. It was decided to use the corrected version of Akaike's Information Criterion (AICc). It is also useful to estimate the probability that the best model selection statistic over all candidate models occurred by chance, that is, could have occurred even if each of the explanatory variables is uncorrelated with the dependent variable. If this probability is low, confidence in the results can be increased in that at least one of the models is likely to have some predictive value. A sanitycheck test was therefore applied to estimate this probability. The test is described in Supporting Information 1 and will be demonstrated further in an upcoming paper.
All analyses were performed using Matlab R2015 v8.5.0. Conditional plots (explanatory variables plotted against the response variable, keeping all other variables constant) of variables included in the best performing models were made with the visreg package F I G U R E 2 Plot of the mean temperature in degrees Celsius over January and February in year i + 1 against the relative population change  (Breheny & Burchett, 2017) of the R ver. 3.5.0 statistical software (R Core Team, 2018).

| Testing for density dependence
We did not detect density dependence in the Hardangervidda reindeer population. Each population for a given year was plotted against the relative population change (where known) over the following year along with the linear regression line (Supporting Information 3). There was no obvious nonlinear relationship, and thus only linear density dependence was tested. The linear regression line has a negative slope but was not found to be significant at any standard level according to the two-sided randomization test (p = .113).

| Effects of hunting
There is little evidence that a larger proportion of animals killed reduces the population. The slope of the linear regression line between the harvest rate and the relative population change was found to be negative (Supporting Information 4), which one would expect if hunting activity tends to reduce the following year's population. The two-sided randomization test yielded a p-value of .109, which is not significant at any standard level. To gain stronger insights into the true effect of hunting on the population, a longer time series would be required.

| Weather effects
Warm winter temperatures had a negative effect on the population. As described in the methods section, we tested for influential observations while performing the linear regression for each variable against the relative population change. An influential observation was found in the case of the mean January and February temperature ( Figure 2). The point denoted with the star has a Cook's Distance greater than one and is therefore categorized as influential. The solid line in the figure shows the linear regression line when that year is included in the analysis and the dashed line shows the line when that point is removed. Parameter values are shown where applicable and are in bold when that parameter is significant at the 5% level. "JF temp," "FM days >0" and GDD correspond to the mean January/February temperature, the number of days in which the temperature went above zero degrees Celsius in February and March, and the growing degree days, respectively. DD and Killed indicate whether density dependence and the proportion of animals killed are included as variables and DD Int. and Killed Int. indicate whether the interaction between the climate variable and density dependence or the proportion of animals killed is included. Each variable is normalized to have mean zero and unit standard deviation and this is reflected in the parameter values.
The year 2010 saw exceptionally cold winter weather, and given the large difference between the January and February temperatures in that year and those in the next coldest year, it is not clear whether this deviation from the regression line resulted because the linear relationship does not extend to such low temperatures, or whether it was caused by some confounding factor in that year (for example, in 2016, over 300 reindeer were struck and killed by a single bolt of lightning). Without an obvious explanation, it was difficult to conclude that the regression line should extend to include this point. It is important to emphasize, however, that this point is not necessarily an outlier but that the question of how the population reacts to unusually cold winters is an open question. Applying the randomization test to all years yielded a p-value of .0142 and .0005 when the year 2010 was left in and out of the analysis, respectively. This provides strong evidence that warm winter temperatures have a negative effect on the population.
The p-values of the linear relationship between each variable and the relative population change are shown in Table 1, both with and without the year 2010 included. The adjusted p-value, to account for multiple testing, is shown in brackets in each case. Here, once multiple testing is taken into account, the mean temperature in January and February is strongly significant (once the year 2010 is omitted) and the number of days above 0 C in February and March is weakly significant. Growing degree-days, March/April and July/August temperature were all not found to be significant at any standard level F I G U R E 3 Conditional plots of the variables in a top performing model (M5) -(a) mean winter temperature in January and February (JF Temp, degrees Celsius), (b) harvest rate, and (c) their interaction plotted against the population growth rate (HR = harvest rate). The gray area represents the 95% confidence interval. The top performing model (M3) is very similar to M5, excluding the interaction, and is therefore not depicted (also see Supporting Information 5). The results of testing lagged weather variables as predictors of the relative population change found no significant results after accounting for multiple tests (see Supporting Information 6).

| Model selection
The AICc for each of the chosen candidate models is shown in Supporting Information 7. Out of the benchmark models, the Stochastic Ricker Persistence performed the best and was thus used as the zero value. Any negative AICc thus outperforms the Stochastic Ricker Persistence forecasts on average. The fact that a persistence approach outperformed the climatology is not surprising since there is clearly an association between consecutive population counts and thus it is clear that knowing the previous population is useful in predicting the following year's population count. The AICc selects the number of days in which the temperature went above zero degrees Celsius in February and March, the proportion of animals killed, and the interaction between the two as the best performing combination of variables. Applying the sanity check test defined in Supporting Information 1 using 4,096 randomized data sets, the probability of observing the best AICc value or lower by chance was estimated to be 0.83%.
As previously noted, the winter of 2009/2010 was far colder than the rest of the years in the data set (see Figure 2) and the observation of the mean temperature in January/February was categorized as an influential observation. Since no other years were close to being as cold, it is difficult to say what the effect of such temperatures might be on the population. Therefore, we performed the analysis without data from that year included, leaving the effect of very cold winter temperatures as an open question. The  Table 2. Here, the best two models according to the AICc include both the January/February temperature and the proportion of animals killed with the latter including the interaction between the two (Figure 3). The third best performing model includes the number of days in which the temperature went above zero degrees Celsius in February and March, the proportion of animals killed and the interaction between the two (Figure 4). The probability of finding an AICc value as low or lower than the best model by chance is estimated to be 0.34% and thus it is unlikely that these results occurred by chance. This helps build confidence in the predictive value of the bestfitting model.

| Analysis of population trends
Over the first few years of the data set, there was a notable decline in the reindeer population, showing a fall from around 20,000 animals in 1997 to just 6,000 in 2003. The population then slowly recovered back up to around 12,000 in 2015. To attempt to understand these trends, we looked more carefully at the behavior of our predictor variables over the length of the data set. The model selection results showed that high mean January/February temperatures and a high proportion of animals killed tends to be negatively correlated with the relative population change. We therefore considered the effect of these two variables ( Figure 5).
Since the proportion of animals killed and the mean January/February temperature are both negatively correlated with the relative population change, we expect years in which one or both of the variables are below the median to coincide with negative observations of the relative population change (i.e., a decline in the population which accounts for half of all years). This is indeed the case. For eight of the years in which the population declined, the mean January/February temperature fell above its median observed value and, in six, the proportion of animals killed was above its median. In 5 years in which the population declined, both the mean January/ February temperature and the proportion of animals killed were above their median observed values.

| DISCUSSION
Our findings suggest that the population of wild reindeer in Hardangervidda was most affected by winter temperature and hunting. Colder temperatures in January and F I G U R E 5 (a) Reindeer population counts in each year of the data set. A black asterisk highlights that there was a drop in the population between the year on the x-axis and the following year. (b) and (c) show the proportion of animals killed and the mean January/ February temperature relevant to the population in the year shown on the x-axis and the following year, respectively. A black asterisk highlights that the relevant variable was above the median in that year (i.e., was in the top 10 highest observations over the data set) February, as well as fewer days with temperatures above zero degrees in February and March tended to result in higher growth rates. The same was true for lower harvest rates. Although we found growing degree-days and higher spring and summer temperatures to have a negative linear relationship with the population growth rate, these patterns were not found to be significant at any standard level. We did not find significant evidence for linear density dependence in the reindeer population.
Towards the beginning of the data set, the proportion of animals killed, and to a lesser extent, the mean January/ February temperature, tended to be relatively high, likely explaining the sharp decline in the population over this period. After 2003, in which the population started to increase again, the proportion of animals killed fell sharply (probably as a direct reaction to decisions aimed at stabilizing the population decline seen in previous years), while mean January/February temperatures also tended to be relatively low over this time. The combination of the two factors probably contributed to the recovery of the population over that period. In general, our results suggest that hunting is one of the main drivers of population dynamics, and is probably the reason that we did not detect density dependence in the population. Thus, in order to prevent severe population declines, hunting quotas should be determined in accordance with winter weather.
Cold winter temperatures may be beneficial to reindeer populations due to a reduction in the incidence of icing, which makes the snow pack impenetrable, and thereby reduces food availability. In contrast, other authors have found that reindeer population growth increases as a result of more periods of warm weather (Tyler, Forchhammer, & Øritsland, 2008). In the year in which the average winter temperature was the lowest (2010, see Figure 2), the population was also negatively affected. Although this is only a single data point, and therefore should be treated with caution, this result suggests that the Hardangervidda reindeer population may exhibit a similar response given more years with colder than average winter temperatures. However, the wild reindeer population in Hardangervidda is at the southern end of its range, and, as such, is likely to be more limited by higher temperatures in general. The timing of warm periods in the winter is also likely to be an important factor. For example, warm periods in the autumn, or generally during times where the snow depth is not as great, are likely to be less deleterious than warm periods in the early spring where icing events make it impossible for the animals to reach suitable forage. The positive effect of cold temperatures for the Hardangervidda reindeer population has implications under future climate change, where temperatures are predicted to increase, and a shift in the timing of seasons is anticipated. Weather variables are likely to become more important for the population if movement is restricted, because the animals become more dependent on a variety of winter grazing areas that can support them through winters with many warm spells. Reindeer management should therefore consider the accessibility of habitat in the regions bordering Hardangervidda, particularly north of the national park. This is a challenge in the light of findings of Chronic Wasting Disease in the reindeer population in that region, where a mixing of populations is not desirable (Hansen et al., 2017).
Although we did not find significant evidence of a linear relationship between summer temperature or growing degree-days and population growth, both variables show a negative trend. It is possible that the negative effects (e.g., insect harassment) cancel out the positive effects of warmer summer temperatures (e.g., more productivity, longer growing season), or that these effects vary spatially. For example, animals that find an easy escape from insect harassment in more topographically complex areas could benefit from better summer grazing conditions that come with a higher number of growing degree-days, while animals in other areas suffer more. Additionally, although warming generally results in more plant productivity, it has also been shown to reduce the variability in plant phenology, which in turn, has been shown to result in a decline in reindeer calf production . If indeed it is true that summer temperatures have conflicting positive and negative effects on the population, the overall relationship may be highly nonlinear and extremely difficult to capture given the small sample size in the study. There are also other factors that can lead to a less clear result, such as variability in wind speed and cloudiness, which have been suggested to affect insect harassment along with temperature (Weladji, Holand, & Almøy, 2003). There is limited evidence from our results that growing degree-days are an important factor in controlling the reindeer population, and more data would be required to investigate this variable. However, growing degree-days and the onset of the summer grazing season are undoubtedly important for the survival of reindeer calves. For example, studies of other reindeer populations have shown that there is a drastic reduction in calf production when growing seasons begin earlier and the animals do not adjust their calving range . However, the calving time in Hardangervidda is at the end of May, which is slightly later than in other Norwegian reindeer areas (Reimers, 2002;Reimers, Klein, & Sørumgård, 1983), so this particular population may not be as vulnerable.
Our results show trends across Hardangervidda, and give an indication of how region-wide weather patterns and hunting pressure can affect the wild reindeer population. The dataset used in this study is relatively small and thus care needs to be taken into the interpretation of the results, as it may be that the sample size is not large enough to yield statistical evidence of such effects. However, as more data is collected and the reindeer counts become more complete, predictions for population growth will also become more reliable. Further investigation into the existence and nature of density dependence and the influence of other climatic and human disturbance related factors would be useful, particularly given the emergence of new data points in the future.