Hunting for ecological indicators: are large herbivore skeleton measures from harvest data useful proxies for monitoring?

Hunter-collected data and samples are used as indices of population performance, and monitoring programs often take advantage of such data as ecological indicators. Here, we establish the relationships between measures of skeleton size (lower jawbone length and hind-leg length) and autumn carcass mass of slaughtered individuals of known age and sex of the high Arctic and endemic Svalbard reindeer (Rangifer tarandus platyrhynchus). We assess these relationships using a long-term monitoring dataset derived from hunted or culled reindeer. The two skeleton measures were generally strongly correlated within age class. Both jaw length (R2 = 0.78) and hind-leg length (R2 = 0.74) represented good proxies of carcass mass. These relationships were primarily due to an age effect (i.e. due to growth) as the skeleton measures reached an asymptotic size at 4–6 years of age. Accordingly, strong positive correlations between skeleton measures and carcass mass were mainly evident at the young age classes (range r [0.45–0.84] for calves and yearlings). For the adults, these relationships weakened due to skeletal growth ceasing in mature animals causing increased variance in mass with age—potentially due to the expected substantial impacts of annual environmental fluctuations. As proxies for carcass mass, skeleton measurements should therefore be limited to young individuals. Although body mass is the ‘gold standard’ in monitoring large herbivores, our results indicate that skeleton measures collected by hunters only provide similar valuable information for young age classes, particularly calves and yearlings. In sum, jaw length and hind-leg length function as proxies identical to body mass when documenting the impacts of changing environmental conditions on important state variables for reindeer and other herbivores inhabiting highly variable environments.


Introduction
Human perturbations, such as climate change, habitat fragmentation, pollution, and harvest, threaten species and ecosystems worldwide (Heleno et al. 2020). In the Northern Hemisphere, particularly in the Arctic, the mean annual temperatures are increasing three times faster than in other biomes (IPCC 2021), putting Arctic ecosystems under increasing pressure (Ims and Ehrich 2013;Meltofte 2013). In such a context, monitoring of state variables or indicators that quantify characteristics of ecological processes and their stochastic dynamics are therefore fundamental (Festa-Bianchet et al. 2017;Ims and Yoccoz 2017). However, logistically, economic and practical constraints often limit the type of data available. The remoteness of the Arctic often challenges effective, systematic field data collection, even for harvested populations of management concern (Christensen et al. 2020).
Data collection protocols for monitoring programs of harvested species are commonly designed to inform management about population performance as a basis for decisions related to, for instance, setting hunting quotas and harvest intensities (e.g. Fryxell et al. 2010;Leclerc et al. 2016). Indeed, hunters may also collect valuable information to incorporate in a monitoring context. This includes data from harvested animals (e.g. sex, age, and carcass mass; Olofsson et al. 2008;Gamelon et al. 2012;Rughetti 2016;Teichman et al. 2016), systematic field observations (e.g. individual observations per hunter day; Ericsson and Wallin 1999;Solberg and Saether 1999;Mysterud et al. 2007), or samples of different body parts (e.g. antlers, lower jawbone, hind-legs, uterus, femur bone marrow; e.g. Hewison et al. 1996;Sand and Cederlund 1996;Azorit et al. 2003;Zannèse et al. 2006;Martin et al. 2013;Rivrud et al. 2013;Becciolini et al. 2016;Cook et al. 2021a, b). For large, long-lived terrestrial herbivores, body mass is strongly related to both survival, reproduction, and, ultimately, population dynamics (e.g. Saether 1997;Taillon et al. 2011;Bårdsen et al. 2014;Albon et al. 2017;Bårdsen 2017). Accordingly, the most frequently used proxies or measures of individual states, such as condition or size, are live body mass or dressed/carcass mass (e.g. Saether 1997;Sand 1996;Bårdsen et al. 2014;Albon et al. 2017). The body mass of ungulates is a function of both skeleton size and body composition (Festa-Bianchet 1998), where the latter reflects individual amounts of fat and muscles (Chan-McLeod et al. 1999;Monteith et al. 2013).
Ungulates inhabiting temperate or Arctic areas experience large seasonal variations in environmental conditions and food availability. This seasonality in resource availability generates predictable, annual variation in growth patterns, the timing of reproductive investment, and consumption of internal energy storages. Environmental conditions can influence both body growth (i.e. skeleton measures) (Klein 1964;Nugent and Frampton 1994;Bertouille and Decrombrugghe 1995;Azorit et al. 2003;Toigo et al. 2006;Zannèse et al. 2006;Berger et al. 2018) and body mass (e.g. Sand 1996;Saether 1997;Bårdsen et al. 2014;Albon et al. 2017). In practice, body mass measures are often unavailable unless in resource-intensive protocols, such as capture-mark-recapture programs (Morellet et al. 2007;Albon et al. 2017;Festa-Bianchet et al. 2017). However, for hunted populations, hunters can be involved in collecting different types of samples or data that may be useful indicators of growth both on the level of individuals and populations. For example, jawbones from hunted animals can be used to determine both skeleton size and age (e.g. Bertouille and Decrombrugghe 1995;Nugent and Frampton 1994;Azorit et al. 2003;Suttie and Mitchell 1983;Høye and Forchhammer 2006), while slaughtered mass can be used as proxies for body condition (Olofsson et al. 2008). However, both skeleton size (e.g. Bertouille and Decrombrugghe 1995;Klein 1964;Martin et al. 2013) and body mass (Knott et al. 2005) typically reach an asymptote at a certain age (i.e. maturity or later), after which only body mass can be substantially influenced by annual variation in intrinsic and extrinsic factors (see references above).
In this study, we combine two long-term time series from hunted and scientifically culled individuals of paired data on jaw length, hind-leg length, and carcass mass in wild Svalbard reindeer (Rangifer tarandus platyrhynchus). We assess the relationships between these metrics and evaluate the potential applicability of skeleton measures from material collected by hunters as monitoring indicators. As males are generally larger than females, and because it takes several years for the reindeer to reach adult size, we expect the relationship between skeleton size and carcass mass to vary across age and sex and to reach asymptote at maturity. Specifically, we predicted a stronger relationship for the younger age classes and this relationship to last longer for males as they grow larger than the females. Based on the sex-and age-specific patterns and co-variation among the skeleton measures and carcass mass using data spanning four decades (Hansen et al. 2012), we discuss the use of skeleton measures in future ecological monitoring of large herbivores in highly unpredictable and seasonal environments.

Study species
The endemic Svalbard reindeer are found in small groups or as single animals year-round. Sexual segregation is evident for much of the year, except during the breeding season (Loe et al. 2006). Unlike most other reindeer, the Svalbard reindeer is neither migratory nor nomadic and tends to use small seasonal home ranges (Tyler and Øritsland 1989). Predation on Svalbard reindeer is negligible (but see Derocher et al. 2000;Stempniewicz et al. 2021), harvest offtake is low (Peeters et al. 2021), and there is no competition with other large herbivores or insect harassment (Halvorsen and Bye 1999). Therefore, besides intraspecific competition for food, the population fluctuations are mainly driven by environmental variability (e.g. Solberg et al. 2001;Aanes et al. 2003;Albon et al. 2017;Hansen et al. 2019). The total Svalbard population has doubled since the 1980s, and the current population size is estimated to be approximately 22,000 individuals (Le Moullec et al. 2019).

Hunting system
The Svalbard reindeer have been harvested since the discovery of the archipelago in the fifteenth century. In 1925, hunting was banned because the species was on the brink of extinction (Lønø 1959). Following a recovery, hunting re-opened in 1983 with an annual offtake of 105-238 individuals until present. This represents around 7% of the total population in the areas (units) where hunting is allowed and is assumed to have limited impact on the population dynamics (Peeters et al. 2021). The hunting season lasts from the 15th of August to the 20th of September and is exclusively for Svalbard residents. Hunting licences are distributed across six hunting units ( Fig. 1) and three age and sex categories: calf, yearlings or adult females, and a free choice. Hunters can harvest any animal independent of sex and age in the latter category, whereas calves can be culled on all licenses. Each hunter is obliged to report the sex and age class (according to quota classes), date of the shooting, and the area they made the kill (hunting unit) to hunting authorities. It is also mandatory to submit the lower jawbone (mandible), primarily used for age assessment. Total length of the jawbone is collected as a measure of overall skeleton size. In some years, during the hunting season, the hunters have also provided carcass mass (see Online Resource 1 for sample sizes). As the hunting grounds cover a large area and are far from roads and other transportation means, harvested animals are usually field dressed, skinned, and cut into several parts to facilitate transportation.

Reindeer data
Two different data sources provided individual information on both skeleton size (lower jawbone length and hindleg length: measured with precision to the nearest mm) and carcass mass (defined as total mass after removing the head, skin, viscera, metapodials, and blood; measured with an accuracy of 0.5 kg) from harvested Svalbard reindeer. The first dataset, obtained from the Governor of Svalbard, was from the regular autumn hunt (5 years from August-September 1984-2012: N = 252; range = 13-87 individuals per year). The second dataset was from a scientific culling program (August-October 1994-2007; N = 388; range = 15-57 individuals per year). Measures of hind-leg length were only available for the latter dataset, which mainly consisted of females (70%). See Online Resource 1 for the distribution of samples by age, sex, and month.
We based the ageing of calves, yearlings, and 2-yearolds on the tooth eruption pattern in the lower mandible (de Bie 1977). Older individuals were aged by counting annuli in the cementum of decalcified, stained, sectioned incisor roots (Veiberg et al. 2020). We followed Langvatn (1977) and measured jaw length and hind-leg length, using callipers, as the maximum length in mm of the lower mandible and as the distance from the tuberculosis of the fibular tarsal bone to the distal end of the metatarsus, respectively.

Exploratory analyses
We fitted Generalised Additive Models (GAMs) separately for each sex, using the mgcv library (Wood 2006(Wood , 2021 in R (R Core Team 2021). We used GAMs to explore age effects in our responses: jaw length, hind-leg length and carcass mass. In these analyses, we used cubic regression splines to model potential nonlinear effects, using the identity link and assuming a Gaussian family [R-code; gam (response ~ s (predictor, bs = "cr", k = 4, gamma = 1.4)]. We extracted predictions, including their precision (± 1 SE), from these models using the predict.gam function in the mgcv package and presented this graphically, along with the underlying data. In addition, we report the estimated effective degrees of freedom (edf), their corresponding level of . Mean [± 1 standard deviation (SD)] of A maximum jaw length (mm), B hind-leg length (mm), and C carcass mass (kg). Lower panel: Exploratory plot of Svalbard reindeer carcass mass as a function of skeletal measures and sex. Mean (± 1 SD) for D jaw length (mm) and E hind-leg length (mm). The thick lines in the plots are the predicted curves from the GAM models (solid for males and dotted for females): edf is the effective degrees of freedom, R adj 2 is the adjusted coefficient of determination, and P is the statistical significance for the degree of smoothing in the GAM significance, and the adjusted R 2 for each GAM. Separated by age and sex, we also ran bivariate tests of the degree of linearity in the relationship, using Pearson's product-moment correlations (r), between the three variables.

Confirmatory analyses
Based on our objectives and available sample sizes, we defined three relatively simple a priori models expressing our 'competing' biological hypotheses. In these models, we predicted carcass mass based on jaw and hind-leg length in separate analyses. In addition to these two continuous variables, we added sex and its interactions with jaw length or hind-leg length to the set of candidate models. A visual inspection of models fitted using the lm-function (R's stats package) indicated problems related to heteroscedasticity in our models: increasing residuals for increased predicted values occurred both when jaw and hind-leg length were predictors (Zuur et al. 2010;Cleasby and Nakagawa 2011). Because we expected both the skeleton measures and carcass mass to reach an asymptote with age, we fitted models with different variance structures. Hence, we fitted both Generalised Least Squares (GLS; i.e. a model very similar to standard lm-fitted models) and Linear Mixed Effect (LME) models (random intercepts across 'Hunting units' only: Bates 2000, Zuur et al. 2009) using the nlme package (Pinheiro et al. 2021). The inclusion of the hunting unit as a random effect means that we estimate and account for both between-variability of the response due to this factor (Pinheiro and Bates 2000). Both GLS and LME open the possibility to add variance structures (for details regarding the variance functions applied, see Pinheiro and Bates 2000, 206-225;Zuur et al. 2009, Chapter 4)-hence we used GLS instead of standard linear models. To account for this, we assessed if the added complexity by including commonly applied variance functions (varPower, var-Exp, and varFixed; Pinheiro et al. 2008) was justified (see Online Resource 2 for details on how we defined the variance structures). The second-order Akaike's Information Criterion (AICc: e.g. Burnham and Anderson 2002;Anderson 2008) was used as the model selection criteria, using the aictab function in the AICcmodavg package (Mazerolle 2020). If the resulting Δ i (i.e. the difference in AICc for model i compared to the model with the lowest AICc value) was ≤ 1.5, we selected the simplest model and used that for inference. Model selection occurred in three steps (Zuur et al. 2009). First, based on the most complex fixed effect structure, we fitted models using restricted maximum likelihood (REML; defined as method = 'REML') as the fixed effects were constant. We defined two candidate models: one with (LME) and without (GLS) random effects. Second, we re-fitted the selected model from the first step, again using a REML, and defined a set of candidate models with different variance structures (defined by the weights argument). Third, the selected model from the second step was re-fitted using a maximum likelihood to select among models representing our three biological hypotheses (we fitted them using ML as the fixed effects differed). Finally, we re-fitted the selected model from step three using a REML and used that model for inference. We plotted predicted values from the fitted linear models using the built-in predict functions in R.

Fig. 3
The plot of predicted values from the best-fitted linear mixed effect model describing Svalbard reindeer carcass mass (log-scale) as a function of A jaw length (centred so that the mean value of jaw length is subtracted from all observed jaw lengths to let the intercept represent predicted carcass mass for an average jaw length) and B hind-leg length. See Table 2 for the estimated effect sizes of the selected models and Online Resource 6 for estimated relations without log transformation

Results
Age explained a high percentage of the variation in the skeleton measures and carcass mass for both sexes ( Fig. 2A-C, Table 1). The most rapid development of skeleton measures and carcass mass occurred among calves and yearlings. After this age, males had a higher growth rate, and growth continued until an older age than females, causing increased sexual dimorphism as the animals grew older ( Fig. 2A-C; Online Resource 3). Jaw length reached an asymptote in females at 5 years of age (mean = 215 mm, SE = 1) and in males at 6 years (mean = 247 mm, SE = 2; Fig. 2A; Online Resource 3). Note that the sample sizes also decreased as a function of age, particularly for males (only 10% of males were older than 6 years). The same pattern was similar for hind-leg length, but growth ceased at an earlier age than for the jaw length, reaching a mean of 279 (± 1 SE) mm at 3 years of age for females and 316 (± 6 SE) mm around 4 years of age for males ( Fig. 2B; Online Resource 3). Given these data and models, carcass mass continued to increase across all age classes for males but reached an asymptote of 36.8 (± 1.2 SE) kg at three years of age for females ( Fig. 2C; Online Resource 3). Because of the overall growth patterns with age, carcass mass was closely related to jaw and hind-leg length (Fig. 2D-E) when not accounting for age. Within age and sex classes, we also found consistently high positive correlations between the two skeleton measures and between carcass mass and the skeleton measures for particularly calf and yearlings of both sexes, while the correlations generally dropped at higher ages and were not significant (Table 1 and Online Resource 4).
In the analyses of carcass mass and its relations to the skeleton measures, we selected a LME model over the GLS model. The fixed effects part of the models was kept constant using our a priori defined most complex fixed effect structure (Online Resource 2). We found heteroscedasticity (i.e. the added complexity acquired through adding a variance structure to the model) in the residuals. However, the type of variance structure did differ as we selected the varExp and the varPower in the analyses of jaw length and hind-leg length, respectively (Online Resource 2). Among the set of candidate models for the jaw length as a function of carcass mass, the simplest model, consisting of only this predictor, was selected in the analysis of carcass mass ( Fig. 3A; Table 2). This means that neither the effect of sex nor its interaction with jaw length was strong enough to be justified in a model explaining carcass mass. Table 2 Estimated effects, including 95% confidence intervals (CI), from the selected models where carcass mass of Svalbard reindeer (on log e -scale) was predicted as a function of skeleton measures: a) jaw length and b) hind-leg length. We centred the continuous predictors meaning that the intercept represents the predictions from the model for the average jaw length (a) and hind-leg length (b), respectively. Sex is a factor variable with females (baseline) and males as levels. We used the treatment contrast. This means that the estimated effect of sex (b) represents the difference, or contrast, between the sexes where males were smaller than females (evaluated for the average hind-leg length; see Figs. 2 and 3 for visual representations of the models, re-fitted without centred continuous predictors, and data). We added different variance structures (through the weight argument in the fitted LME models) to account for heteroscedasticity in the residuals by selecting the model with the lowest AICc value In contrast, in the analysis of hind-leg length, we selected a model that included the effect of sex and its interaction with hind-leg length ( Fig. 3B; Table 2). However, we selected the same random structure in both analyses, showing only a small across-area variance (relative to the intercept) in the relationships between carcass mass and the skeleton measures (Table 2). These relatively simple models explained a large proportion of the variance in carcass mass for the models using jaw length (R 2 = 0.78) and hindleg length (R 2 = 0.74) as predictors. Another characteristic of these relationships was that both the selected models included a variance function showing that the variation in carcass mass increased with age. See Online Resource 5 for model outputs, and Online Resource 6 ( Fig. S2) for plots of model predictions on arithmetic scale and sex and agespecific relationships on log-log scale.

Discussion
In this paper, we combined hunted and culled Svalbard reindeer data to establish the relationships between carcass mass and skeleton measures. Unsurprisingly, because of strong patterns of asymptotic growth with age, both skeleton measures can be good proxies for carcass mass, but within age, the relationships quickly weaken with increasing age. Skeleton measurements as proxies for carcass mass should therefore be limited to young individuals.
Age was a good predictor of carcass mass and skeleton measures early in life. After reaching a certain age, this effect reached an asymptote (see Fig. 2). Such a growth pattern corresponds to many large herbivores (e.g. red deer Cervus elaphus [Bertouille andDecrombrugghe 1995, Langvatn et al. 2004]; Sitka black-tailed deer Odocoileus hemionus sitkensis [Klein 1964]; roe-deer Capreolus capreolus [Høye and Forchhammer 2006]). The age at which the skeleton growth reached its asymptote differed between females and males. Among females, hind-leg length and jaw length growth ceased around 3 and 4-5 years of age, respectively. In males, the corresponding age of growth cessation occurred when animals reached 4 and 5-6 years of age. In females, the weakened correlation between the skeleton measurements and the carcass mass (abruptly from 2 years in hind-leg length and more gradually for jaw length) probably reflects how resources are allocated differently following the individual onset of reproduction. Female Svalbard reindeer produce their first calf the summer they turn 2, or more often, 3 years old (Tyler and Øritsland 1987). At older ages, individual variation in females will be influenced by reproductive status (i.e. current and delayed costs of reproduction: Albon et al. 2017;Veiberg et al. 2017), resource variation (Aanes et al. 2000;Solberg et al. 2001), parasitic burden (Albon et al. 2002), or differences in survival rates related to size (i.e. larger individuals survive to older ages (Bårdsen et al. 2011(Bårdsen et al. , 2017. In ungulates, body mass is a good predictor of both survival and reproduction (e.g. Albon et al. 2017;Bårdsen et al. 2011;Veiberg et al. 2017;Taillon et al. 2011) and therefore considered the 'gold standard' when monitoring their individual performance (reviewed in Saether 1997;Gaillard et al. 2000;Morellet et al. 2007). Fluctuations in body mass primarily reflect variation in body fat, the primary energetic store, which is particularly important for Rangifer (Parker et al. 2009;Cook et al. 2021a), and maybe even more critical for this high-Arctic sub-species. Skeletal growth is likely prioritised over acquisition of body fat reserves, so there may not be detectable variation in skeletal growth during the younger years, when body fat and therefore body mass are lower. In such a case, only extreme environmental conditions and nutritional deficiencies may be detected using skeletal measurements. However, for monitoring systems where body mass data are lacking, our results from Svalbard reindeer indicate that skeleton measures, such as jaw length and hind-leg length, may serve as alternative proxies for the growth of younger individuals. At older ages, these relationships vanish, as demonstrated by the increased variance in carcass mass with age, and skeleton measures are not applicable for monitoring. From a statistical point of view, heteroscedasticity in the model residuals, as in our case, represents a violation of assumptions behind standard regression models (Zuur et al. 2009;Pinheiro 2021), which we accounted for in our models. In our analyses, this violation estimates and accounts for a biological process resulting in more residual variance for the large-bodied and adult individuals in our population.
Even when carcass mass (i.e. slaughtered weight) is available from the hunters, skeleton measures may be a preferred proxy for young individuals as they are likely more precise. Slaughter mass is typically influenced by factors related to varying weighing practices between hunters or herders (Olofsson et al. 2008). In addition, skeleton measures are likely less affected by seasonal dynamics (Martin et al. 2013;Becciolini et al. 2016), which could both be an advantage and a disadvantage depending on the objectives of a particular study (e.g. Albon et al. 2017). However, body mass is likely to be a better predictor of individual condition at older ages, as suggested in other studies (e.g. Becciolini et al. 2016;Berger et al. 2018;Langvatn et al. 2004;Suttie and Mitchell 1983;Klein 1964). Moreover, because body mass fluctuates between seasons , the use of the skeletal measures must be strictly related to autumn and not related to other seasons. Further, body mass and likely skeletal growth is influenced by maternal effects (Bernardo 1996;Monteith et al. 2009;Freeman et al. 2013;Michel et al. 2016). Thus, a combination of direct and lagged (potentially for > 1 year) environmental conditions may also influence skeletal growth. For instance, extremely harsh weather conditions that result in declines in body mass could be detected by carcass mass more readily and for more age classes than skeletal measurements. Although less accurate, estimates of carcass mass are more useful than skeletal measurements, at least for older individuals.
Ecosystem monitoring programs should collect a set of indicators or state variables, reflecting different aspects of the targeted species and the relevant ecological processes (Morellet et al. 2007;Ims and Yoccoz 2017). Choosing such indicators is a cost-benefit analysis where different state variables have diverse applicability and accessibility. In our case, complete records of carcass mass and skeleton measures from all harvested Svalbard reindeer would be ideal for monitoring, providing several proxies reflecting different aspects of growth, body size, and condition (e.g. mass by size). In practice, however, carcass mass is both challenging to collect in the field and inaccurate due to varying weighing practices. Therefore, continued collection of skeleton measures should be encouraged for animals in their growth phase, particularly for calves and yearlings. Collecting the samples (hunters) and measurements (lab personnel) are more easily obtained for jaw bones than hind-legs, and their correlations with carcass mass were also generally higher, even with increasing ages. In addition, and importantly, access to the lower jaw provides teeth for age determination. The jaw length is also known to be a part of the skeleton that possesses the highest growth priority (Langvatn et al. 2004). Therefore, we recommend the local management prioritise hunter collection of lower jawbones from at least young individuals of both sexes (i.e. calves and individuals 1-2 years of age). Given the relationships established here, the historical archive of jawbones (from 1983 to present) has opened up an opportunity to track reindeer growth patterns for young individuals (e.g. potential cohort effects; Hamel et al. 2016) and how and why they vary in space and time related to, for example, density-dependent population processes and climatic drivers. See Online Resource 7 for an example of within and between variation for calves and yearlings.
Svalbard is a 'hot spot' for climate change where substantial impacts of annual environmental variation have been documented on reindeer population growth rates (e.g. Hansen et al. 2013Hansen et al. , 2019Albon et al. 2017), age and sexspecific mortality (Peeters et al. 2017), and reproduction (Stien et al. 2012;Albon et al. 2017). Less is known about the mechanistic pathways these processes operate, although individual and annual variation in late winter body mass is clearly of considerable importance ). Environmental variation is thus likely to explain some of the observed variations in age-specific skeleton metrics and growth patterns. Based on four decades of huntercollected skeleton data, this study provides insights into these questions to better understand how body size and growth respond to changing environmental conditions (see Berger et al. 2018 for an example) in this high-Arctic endemic reindeer species.