Cross-sectional properties of reindeer long bones and metapodials allow identification of activity patterns

Habitual loading patterns of domesticated animals may differ due to human influence from their wild counterparts. In the early stages of human-reindeer interaction, cargo and draft use was likely important, as well as corralling tame reindeer. This may result to changes in loading as increased (working) or decreased (captive) loading, as well as foraging patterns (digging for lichen from under the snow versus fed working and/or captive reindeer). Our aim is to study whether differences in activity modify variation in bone cross-sectional properties and external dimensions. Our material consists of donated skeletons of modern reindeer: 20 working reindeer (19 racing and one draft), 24 zoo reindeer, and sample of 78 free-ranging/wild reindeer as a reference group. We used general linear modelling to first establish the total variation in cross-sectional properties among wild and free-ranging reindeer, and then to infer how differences in loading modify observed variation among zoo and working reindeer. According to our results, direction of greater bone quantity as well as external dimensions in of radioulna of female reindeer differs from female reference group, likely relating to foraging behavior. External dimensions of humerus differ in working and zoo male reindeer compared to male reference group. Increased robusticity of long bones, especially of tibia among working male reindeer, may indicate increased loading, and increased cortical area of long bones may indicate sedentary lifestyle among female reindeer. The results of this study can be used to understand early stages of reindeer domestication by observing reindeer activity patterns from archaeological material.


Introduction
Habitual loading patterns representing physical activity of a human or an animal person are considered to result in modifications in bone mass and its distribution (Shackelford et al. 2013;Cameron and Pfeiffer 2014;Ruff and Larsen 2014;Ruff et al. 2015). Due to sustained interaction with humans, physical activity patterns in the skeletons of domesticated animals can be different from the patterns observed in the skeletons of their wild counterparts (Flensborg and Kaufmann 2012). The different patterns typically emerge as a response to different loading environment, such as working, living conditions, being fed, or captives under human influence compared to activities in the wild. Reindeer management in historic as well as modern times differs from other domesticated animals, where reindeer roam free most of the year except biannual roundups. Due to elusive nature of identifying the early stages of human-reindeer interaction and reindeer domestication (Salmi and Heino 2018), distinguishing the effects of modern human influence on reindeer activity might allow investigating the earlier contacts between human and reindeer.
Reindeer have had a major role in the subsistence, worldview, and religion of many peoples in the Eurasian Arctic reaching to modern times (Gordon 1990;Helskog 2011). It has been argued that cargo and draft use of reindeer were important in the early stages of reindeer domestication (e.g., Ingold 1986;Bjørklund 2013), which probably began in Scandinavia during the Late Iron Age, ca. 800-900 AD onwards (Bergman et al. 2013;Bjørklund 2013). Tame reindeer were also kept at the homestead (Bately 2007: 56), likely corralled or tethered. Thus, human-reindeer interaction causes changes in habitual loading as increased (working) or decreased (captive) loading, as well as foraging patterns (digging for lichen from under the snow versus fed working and/or captive reindeer). These changes can potentially be observed in bone cross-sectional properties, namely bone strength and direction of greater quantity of bone.
Bone functional adaptation-or mechanostat theoryexplains changes in bone robusticity with activity, as bone is adapted to "normal" level of mechanical loading with removal or accrual of bone mass to appropriate to the level of mechanical loading. Mechanical loading occurs in terms of gravitational and muscle loading (Lanyon and Rubin 1984;Frost 1987;Lanyon 1987Lanyon , 1996Heinonen et al. 2002;Judex and Carlson 2009;Robling 2009), where the former seems to be more important than the latter (Niinimäki et al. 2019). In case of working reindeer, both loading elements are present as increased muscle loading from pulling/pushing a load and increased gravitational loading from walking/running motion while pulling/pushing a load. In addition, frequency and magnitude of loading contribute to bone robusticity (Lanyon and Rubin 1984;Lanyon 1987Lanyon , 1996Umemura et al. 2002;Ruff et al. 2006;Nikander et al. 2006;Niinimäki et al. 2017Niinimäki et al. , 2019, where magnitude of loading is increased in working reindeer due to additional weight compared to non-working reindeer. Activity-related modifications in bone mass and its distribution are based on beam theory (Lieberman et al. 2004;Ruff et al. 2006;Stock and Shaw 2007) and mechanostat theory (Frost 1987;Turner 1998), the latter later termed as functional muscle-bone unit (Schoenau and Frost 2002). According to beam theory, the bone is strongest the further the bone tissue is distributed from the cross-section centroid. Thus, bone adaptation to habitual loading can be observed in the distribution of bone in its cross-section and as modifications of external dimensions of bone shaft towards the principal or usual direction of stress (Forwood and Burr 1993;Ruff et al. 2006;Shaw and Stock 2009). There are several studies where mobility has been studied successfully from bone cross-sections (e.g., Cameron and Pfeiffer 2014;Ruff and Larsen 2014;Ruff et al. 2015). In case of reindeer, the time period during which the free-ranging reindeer must dig for lichen is approximately 4 months of the year (i.e., months when the ground is snowcovered), digging for up to 8 h per day (Itkonen 1948: 82;Helle 1982: 47-59;Nieminen and Pietilä 1999: 20-21;Korhonen 2008: 40). Thus, it can be expected that forelimb bone elements may exhibit differences in distribution of bone within a cross-section between fed reindeer and freely grazing reindeer. Working as cargo and draft animal may also be reflected in distribution of bone mass in reindeer bone due to greater requirement for forward propelling motion.
The aim of this paper is to distinguish potential mediators for contemporary variation in reindeer bone mass and its distribution using bone cross-sectional data of working (racing and draft) and zoo reindeer compared to a sample of freeranging/wild reindeer specimens. We first establish the baseline variability in cross-sectional properties using the material on wild and free-ranging reindeer, after which we infer how differences in loading could modify that variation in the groups of working and zoo reindeer. We expect to find that the wild and free-ranging reindeer would show high heterogeneity in bone elements whereas the different loading types would result in directional modification of bones and thus more homogeneous variation in the cross-sectional properties. Based on beam theory, we expect that the distribution of bone in forelimb bone elements between fed and freely grazing reindeer is different. Working in reindeer may potentially also result to different distribution of bone quantity compared to free-ranging/wild reindeer due to requirement of forward propelling motion while pulling a load. We further expect that, based on bone functional adaptation, working activities should result to increased loading and thus to greater bone robusticity, whereas captivity should result to decreased loading and thus to decreased bone robusticity, in comparison to the loading environment of free-ranging/wild reindeer. The results of this study are essential in identifying reindeer activity patterns from archaeological material in view of its domestication by humans. As our material comprises of modern reindeer (working, zoo, and free-ranging/wild), we also evaluate the similarities in the activity patterns between contemporary working reindeer to historic working reindeer based on ethnographic data.

Material
Our sample comprises a total of 122 donated skeletons that we categorized in three activity types: working, zoo, or free-ranging/wild. All the skeletons are not complete (i.e., not all bone elements are present for most of the individuals) due to practical decisions taken during the skeleton collection work.
The contemporary free-ranging and wild reindeer (males N = 44, females N = 34) were used as the reference sample for activity changes. The skeletons were collected and curated by the Biodiversity Unit at the University of Oulu. This material consists of two subspecies of Rangifer tarandus in Finland: semidomesticated reindeer Rangifer tarandus tarandus (referred to hereafter as free-ranging) and wild forest reindeer Rangifer tarandus fennicus, as well as some known interbreeds (hybrids) in zoo sample (but not all zoo individuals were hybrids). These two subspecies have overlapping distribution in Finland and can freely interbreed. The free-ranging reindeer, as their name suggests, roam freely most of the year, except round-ups in the summer for earmarking young calves and in the autumn for selecting individuals for slaughter. Freeranging reindeer are mostly feeding on natural pastures, but some free-ranging reindeer may have received additional fodder brought in specific feeding places, which is a common practice among reindeer herders in Finland. Wild forest reindeer are not managed. Wild forest reindeer has slenderer anatomy with longer legs compared to shorter and stouter semidomesticated subspecies (Nieminen and Helle 1980). The potential effects of size difference between the two subspecies were considered in the statistical analyses by including a size variable. The working reindeer (racing reindeer, and one reindeer that worked as a draft reindeer in tourism: all males, N = 20) of Northern Ostrobothnia region was covering the southern parts of reindeer herding area of Finland. These skeletons were collected by S.N. and A.-K.S. and are currently archived at the Laboratory of Archaeology at the University of Oulu. The following information for the working reindeer in this sample was retrieved from the reindeer herders by S.N., and the described practices are comparable to those elsewhere in Finnish reindeer herding district (Soppela et al. in press). Modern working reindeer work approximately 4 months of the year, i.e., during months with snow-covered ground. This is roughly the same period of the year free-ranging/wild reindeer dig for lichen. Racing reindeer compete for speed, where training includes stamina and speed training. In the competition event, the reindeer is harnessed to pull its driver (weighing 60-65 kg) on skis. Racing reindeer training begins by building stamina with exercises totaling to 3 h per week, starting in January. During competition season, there is no other exercise, but competition event is followed by short recovery exercise, such as walking, rounding up to an hour per week. Thus, working (training or competition) hours for racing reindeer vary from 1 to 3 h per week. This information pertains to the working reindeer in this sample; training habits vary between reindeer owners (Soppela et al. in press). Draft reindeer working in tourism begin building stamina in November, quickly succeeded by working as soon as there is snow. Working hours are a couple hours per day, most days a week (a day or two per week is reserved for rest, or as much as needed for the welfare of a reindeer), pulling a weight between 150 and 250 kg (sledge and maximum of two adults) for rounds of 300 m or 1500 m. Average years of working career in our sample were between 3 and 4 years. The working reindeer in our data were free-ranging (i.e., digging for lichen) before they were selected to training, and after that they were fed, and corralled when not working during winter. During summer, usually all working reindeer graze freely. As this one draft reindeer in our sample was not an outlier compared to racing reindeer sample, these two activities were pooled into "working reindeer." It should also be noted that there are castrated and non-castrated males in both free-ranging and working reindeer samples. However, there were too few individuals with known castrations status to be included in the detailed analyses or to study the potential effects of castration.
The zoo reindeer (mainly from Oulu University zoo, but two individuals from Ranua zoo; males N = 8, females N = 16) were collected and curated by the Biodiversity Unit at the University of Oulu. Zoo reindeer spent their lifetime in captivity. Zoo reindeer management include corralling and feeding. The reindeer in Oulu University zoo were corralled in 570 m 2 enclosure with flat terrain. There were known hybrids among the zoo individuals. For statistical analysis, due to small number of individuals per subspecies group, hybrids were pooled with Rangifer tarandus tarandus, and this was justified by visual examination of variability of bone cortical areas relative to bone length, where the relationship in hybrids was similar to Rangifer tarandus tarandus. These subspecies interbreed in the wild at the southern parts of reindeer herding area, and hybrids are inevitably included in the archaeological samples as well.

Methods
Definitive work by Christopher Ruff (e.g., Ruff and Hayes 1983;Ruff 1992Ruff , 2003Ruff et al. 1993Ruff et al. , 1994Ruff et al. , 2015 has provided the baseline on how to standardize a method for studying cross-sectional properties of long bones. However, as different species have evolved to specific environments and requirements of movement, these orientations should be observed species-specifically. Firstly, it is important to define an appropriate orientation of bone for obtaining cross-sectional properties as this defines the planes for calculated second moments of areas on x and y axes. This is especially relevant for observing the direction of greater distribution of bone mass as inferences of activity are based on these indicators. In the Appendix 1, we describe comprehensively how we oriented the reindeer bones for obtaining the cross-sectional properties. Second, the selection of cross-sectional location is usually defined to a percentage of a bone length, most commonly to interarticular (IA) length (Ruff and Hayes 1983;Sládek et al. 2010;Davies and Stock 2014;Mongle et al. 2015). IA length is considered to better represent functional muscle levers and thus to be more relevant for studying activity than maximum bone length (Ruff and Hayes 1983;Sládek et al. 2010). As IA have not been defined for reindeer long bones, we provide the definition of IA for all long bones and metapodials (Appendix 1). Only radioulna interarticular length is similar to previously described measurement: physiological length as described by von den Driesch (1976).
Bones were scanned using a peripheral quantity computed tomography (pQCT) scanner (Stratec XCT Research SA+), and it provides information on several bone properties of the scanned cross-section. All bones elements were scanned at 50% of IA length; examples of bone specific cross-sections at 50% are presented in Fig. 1. Bone geometric properties include cortical area (CA) which is the cross-sectional area of cortical bone at the cross-section and moments of areas (I x and I y ). A measure of the bending rigidity of the bone (J) is computed by summing the area moments of inertia by summing moments of areas I x and I y . Direction of greater quantity of bone was calculated as a ratio between moments of areas (I x /I y ). We also measured external dimensions at these 50% cross-sectional locations using small sliding calipers on anteroposterior and mediolateral planes and calculated a shape ratio using these external measurements (anteroposterior/mediolateral) As this methodology is eventually intended to be used in archaeological material to assess activity status, external dimensions were included in order to study their applicability in cases when obtaining cross-sectional properties is not possible.
Body mass contributes to intrinsic loading on bones (Seeman et al. 1996;Schiessl et al. 1998;Ruff 2000;Schoenau and Frost 2002;Brianza et al. 2007), where bone tissue is distributed proportionally further from the crosssectional centroid with increased body mass (Brianza et al. 2007). Thus, it is crucial to control for body size effects on bone robusticity before inferring the effects of activity. Due to absence of body mass information for several individuals in this sample, we used bone specific IA bone length as a size proxy in the analyses (Ruff et al. 1993;Feik et al. 1996).
We performed the statistical analyses in IBM SPSS Statistics for Windows, version 25.0 (IBM Corp. 2017). First, we used the reference group of wild and free-ranging reindeer to describe the baseline variation and to test the influence of natural sources (size, sex, subspecies) of variability in each response variable (CA, J, direction of greater bone quantity and external dimensions of each bone element). This was performed by fitting general linear model on each response variable and assessing sex and subspecies as fixed factors and body size (IA length) as a linear covariate.
Based on the analyses for background variation in bone elements (Table 1) and their visual examination, the response variables showed weak to strong bimodality of distribution in response variables due to different sexes. As also the sampling of activity groups was biased regarding sex (i.e., working reindeer consisted only males), we handled sexes separately in the following analyses on activity effects. The effect of subspecies on the observed variation was negligible most cases (similar range and distribution of response values for both subspecies) and thus the subspecies were pooled in further analyses for those variables indicated in Table 1.
To test whether the response variables differ due to differences in reindeer activity patterns (within sex), we performed general linear model (GLM) for each response variable. In the models, the activity type (working, zoo, or reference group of free-ranging/wild individuals) was set as a fixed factor in the models. Body size proxy (relevant bone length) was included as a linear covariate in the models to simultaneously account for strongly size-dependent variation in the response variables. Finally, for those response variables where activity was found a significant covariate, the equality of variance (heterogeneity or homogeneity) between activity groups were tested with Kruskal-Wallis test.

Sources of background variation in bone response variables
The bone elements showed relatively high heterogeneity among the wild reindeer. Baseline variation in bone robusticity variables (CA and J) in the reference group of wild Fig. 1 Examples of reindeer long bone and metapodial crosssections at 50% of interarticular length and free-ranging reindeer indicated significant variation due to body size (represented as IA bone length) where increasing body size was associated with greater values in bone robusticity variables (Table 1). Examination of the baseline variation in the orientations of greater bone quantity and external dimension was mainly unaffected by body size, as would be expected due these variables being indices and thus readily controlled for size. Body size was significant only for the external dimensions of tibia (Table 1).
However, there were significant main effect of sex indicating that the males had on average larger values in bone robusticity variables relative to female reindeers, and this observation was not related to greater body size in males (Table 1). Further visual examination of the range and variance of different bone elements showed that the sex-specific average values as well as their variance resulted in a bimodal distribution of measured variables. While sex was not significant covariate for external humeral dimension (Table 1), distribution of the values of this variable indicated towards bimodal distribution based on sex, and therefore final model was performed for sexes separately.
The variation (mean, variance) due to subspecies was negligible in most cases. The main effect of subspecies was significant only for a few response variables (Table 1). Apart from metacarpal CA, subspecies status did not have significant effect on bone robusticity variables (Table 1). When examining distribution of metacarpal CA values, bimodal distribution was observed according to sex but not subspecies; therefore, metacarpal CA was performed the same statistical treatment as for other bone robusticity (CA and J) variables. For bone quantity orientation variables, subspecies was significant in humerus and metacarpal orientation of greater bone quantity, and metacarpal external dimensions (Table 1). Thus, subsequent analyses of activity effects for these response variables were performed within Rangifer tarandus taradus individuals only (Table 1). As subspecies was not significant to most variables (Table 1), hybrids included in the zoo sample was considered not to bias results on variation based on activity.

Activity effects on bone response variables
To evaluate whether the bone response variables would reveal human influence in the activities, we present parameter estimates resulting from GLM analyses to estimate the average change in all response variables when the zoo reindeer or the working reindeer are contrasted to the reference group of freeranging/wild reindeer (Table 2). Furthermore, we will illustrate the differences in mean bone properties among all three activity types as predicted by the final models (in terms of Table 1 General linear model statistics assessing the effects of sex and subspecies as fixed factors and body size (IA length) as a linear covariate on response variables within wild and free-ranging reindeer, i.e., reference group of reindeer outside human interaction. Final column sums response variables indicating significant variation due to activity from  (Tables 1 and 2). The results from GLM showed that activity contributed significantly to the variation observed in some of the bone response variables (Table 2). Human influence resulted in homogeneity of variance and bias towards x-axis/ mediolateral orientation in the bone orientation variables (when significant) compared to the references group (Table 2, Fig. 2). For bone robusticity variables, human influence resulted in greater values in CA, but greater J values in working (male) reindeer and smaller J female zoo reindeer, when significant (Table 2, Figs. 3 and 4). Regarding heterogeneity versus homogeneity of variance for bone robusticity variables, human influence resulted in greater heterogeneity of variance (when significant) apart from tibia for working (male) reindeer where variance exhibited greater homogeneity compared to the reference (male) group (Table 2).
Activity mainly did not have significant effects on variation in bone orientation (Table 2). Variation due to working and zoo environment in male reindeer in external orientation of bone in humerus was smaller (confirmed with Kruskal-Wallis test) and biased towards x-axis compared to male reference group (Table 2b, Fig. 2). In females, variation due to zoo environment in the direction of greater bone quantity and external orientation of bone in radioulna was smaller (confirmed with Kruskal-Wallis test) and biased towards x-axis compared to female reference group (Table 2a and b, Fig. 2).
Activity contributed also to the variation in the bone robusticity variables. In males, variation due to working in J and CA for humerus, femur, and tibia was biased towards greater values compared to male reference group (Table 2c and d). Variance attributable to working reindeer was greater for humerus and femur bone robusticity variables, and this observation was confirmed with Kruskal-Wallis for all but humerus J (Figs. 3 and 4). In tibia, variance due to working was smaller compared to male reference group (confirmed with Kruskal-Wallis test; Figs. 3 and 4). In addition, variation in J and CA in radioulna was biased towards greater values in working males compared to male reference group, but this this observation reached only near statistical significance (Table 2c and d). In females, variation due to zoo environment in CA for humerus, radioulna, femur, and tibia was biased towards greater values; in addition, variance in CA was greater among zoo females (confirmed statistical significance with Kruskal-Wallis test for all but radioulna) compared to female reference group (Table 2c and d, Fig. 3). In addition, variation in J in radioulna and near-significantly also variation in J in humerus was biased towards smaller values in zoo environment compared to female reference group (Table 2c and d, Figs. 3 and 4).

Discussion
When baseline variation was observed within wild and freeranging subspecies of Rangifer tarandus bone variables, subspecies status did not bias bone robusticity variables as could be expected as the differences in body morphology are very slight, and no differences in intensity of activity would be expected. Therefore, it is evident that bone robusticity can be utilized to identify sustained human interaction in different reindeer subspecies in archaeological material. Body size (as interarticular bone length) and sex did contribute to the variation observed in the bone robusticity variables. Effects of sex sand body size were more varied for bone orientation variables, for some even subspecies status.
Direction of greater bone quantity and/or orientation of external shape of humerus and radioulna allowed observations of human influence (feeding, working, and captive) relative to reference group (free-ranging/wild) of reindeer. Greater relative bone robusticity of long bones, especially tibia, allowed identification of working reindeer from free-ranging/wild reindeer among males, and greater CA of long bones allowed identification of female zoo reindeer from female free-ranging/wild reindeer; in addition, zoo environment resulted to smaller bending and torsion rigidity in zoo female radioulna compared to female reindeer reference group.
Based on beam theory, we expected to find changes in the orientation of greater quantity of bone and external orientation of bone in the forelimb bone elements induced by differences in reindeer foraging behavior. As expected, the orientation of greater quantity of bone of radioulna in fed (i.e., zoo) female reindeer differed from freely grazing (i.e., free-ranging/wild reference group) female reindeer (Fig. 2). The orientation of greater quantity in radioulna was more mediolateral in fed female reindeer (zoo) compared to freely grazing reindeer. Differences in orientation of bone within a cross-section are considered a result of activity patterns, usually commenced in childhood, and continued through adult life (Ruff et al. 2006). In case of reindeer, foraging as digging for lichen commence in their first year of life, and in case of free-ranging/wild reindeer, continues throughout life. Free-ranging/wild reindeer dig for lichen, spending many hours per day digging during snow-covered ground of the year in this activity, whereas working and zoo reindeer are fed during winter in addition to which they are kept in captivity. In case of working reindeer, in ethnographic and modern-day reindeer data, training for pulling a sledge/toboggan begins usually when a reindeer is 3-4 years old (Itkonen 1948;Niinimäki, pers. comm.;Soppela et al. in press), i.e., at final stages of skeletal maturity (Takken Beijersbergen and Hufthammer 2012). It should be noted that the working reindeer in our data were free-ranging before they started training. As the working reindeer differed significantly from free-ranging/wild reindeer, this suggests that the observed differences likely emerged later in their Table 2 General linear model statistics for the effects of activity on a) direction of greatest bone quantity (I x /I y ); and b) bone external dimensions (anteroposterior/mediolateral); and with and a proxy for body size (bone element specific interarticular length) for c) cortical area (CA); and d) bending and torsion rigidity (J), sex-specifically. The fixed treatment effects are presented parameter estimates (β) with standard error, t-tests used for statistical significance and 95 % confidence interval limits  working life. Differences between fed and freely grazing reindeer were evident also as increased robusticity at the entheses for elbow flexor muscles among freely grazing reindeer, proposed to result from greater demands for elbow flexor muscles in the latter (Niinimäki and Salmi 2016). Based on beam theory, we expected to observe changes in the direction of greater quantity of bone in working (male) reindeer compared to male reindeer reference group (i.e., free-ranging/wild) as a result of requirement for forward propelling motion while pulling a load in the former. However, we observed changes only in the external orientation of bone of humerus as biased towards more mediolateral orientation in working reindeer compared to male reindeer reference group (Fig. 2). It is possible, at least in part, that bias towards mediolateral orientation observed for humerus in working reindeer compared to male reindeer reference group (Fig.  2a, b) may be enhanced due to greater demands for forward propelling motion. Mediolateral broadening as a result of pulling has been identified of the metapodials in draft cattle, especially broadening of the medial half (Bartosiewicz and Gál 2013: 144-145), but no activity-associated changes could be observed for the metapodials-or any other bone element-in this reindeer sample. This cautions against direct comparisons between species in bone shape orientation with activity, as it is likely that species-specific joint configurations have their baseline effect on bone shape.  Fig. 2 The estimated marginal means ± SE for model predicted differences in orientation of bone quantity between the activity types, sex-specifically and when activity was significant (Tables 1 and 2) for a) humerus external orientation of males; b) radioulna cross-sectional orientation of females; and c) radioulna external orientation of females Based on bone functional adaptation theory, we expected to find a bias towards an increase in bone robusticity variables of the working reindeer compared to the free-ranging/wild reindeer to match the increased frequency and magnitude of loading of the former. As expected, we observed a bias towards increased bone robusticity values (CA and J) in working reindeer compared to male reindeer reference group in long bones (Figs. 3 and 4). However, this was not observed for metapodials (Figs. 3 and 4). It is possible that metapodials do not respond similarly to frequency and magnitude of loading compared to more proximal limb elements. Furthermore, the bias observed towards increased bone robusticity in tibia was associated with smaller variation in working reindeer compared to male reindeer reference group, suggesting that observations in greater bone strength in tibia may allow reconstructing the activity patterns in archaeological material. Interestingly, contrary to assumed decreased loading in zoo reindeer, we found that female zoo reindeer exhibited a bias Fig. 3 The estimated marginal means ± SE for model predicted differences in cortical area (CA) between the activity types, sex-specifically and when activity was significant (Tables 1 and 2) Fig. 4 The estimated marginal means ± SE for model predicted differences in bending and torsion rigidity (J) between the activity types, , sex-specifically and when activity was significant (Tables 1 and 2) for humerus, radioulna, femur, and tibia for males in top row and for humerus and radioulna for females in bottom row towards greater CA values in long bones (Fig. 3) compared to female reindeer reference group. A possible explanation for this might be that prolonged standing/sedentary lifestyle of zoo reindeer may result to increased axial loading and/or muscle loading, the latter especially evident in shoulder/hip bracing muscles. This was further supported by the bias towards smaller values in J among zoo females compared to female reference group, which is most likely a result of zoo reindeer being fed. Increased axial loading experienced by zoo reindeer is supported by the geometric shape of zoo reindeer shoulder and elbow joints. Joint morphology has modified to a more sedentary lifestyle, resulting in greater joint stability for standing (Pelletier et al. 2020). Greater muscle loading experienced by zoo reindeer in shoulder bracing muscles is supported by greater entheseal ruggedness observed in Subscapularis muscle attachment in zoo reindeer relative to free-ranging/wild reindeer (Niinimäki and Salmi 2016).
Potentially useful bone traits for identifying fed and freely grazing female reindeer are both external shape and direction of greater quantity of bone in radioulna (Fig. 2). Potentially useful indicators of working (male) status (relative to freeranging/wild male reindeer) are external dimensions and in humerus and bone robusticity of tibia (Fig. 2). Potential indicator of sedentary behavior in female reindeer may be increased CA values in long bones and decreased J values in radioulna.
Comparability of modern reindeer and reindeer from archaeological context Itkonen (1948) describes training of transport reindeer, which is similar to the training in our modern working reindeer (Niinimäki, pers. comm.), as well as in general for modern working reindeer in Finland (Soppela et al. in press). Training starts by tethering the reindeer (to get used to the rope as well as humans), followed by walking next to a human on a leash (Itkonen 1948: 419;Soppela et al. in press;Niinimäki, pers. comm.). Then, the reindeer is accustomed to the harness, and then to pulling a load. The weight of the load is gradually increased, and the trainer starts to teach steering and turning (Itkonen 1948: 421-422;Soppela et al. in press;Niinimäki, pers. comm.). Of course, while training in ethnographic accounts and modern-day working reindeer are similar, it may be different in archaeological material.
There are potential differences in frequency and magnitude of loading between modern-day working reindeer and ethnographic data. In traditional Sámi reindeer pastoralism, reindeer were used to pull sledges/toboggans (with people or goods) and carrying loads (Itkonen 1948: 388-391;Näkkäläjärvi and Pennanen 2000;Bjørklund 2013). Cargo reindeer carried loads of ca. 25-35 kg (Itkonen 1948: 388-391;Näkkäläjärvi and Pennanen 2000). This is somewhat lighter load than in our modern working reindeer (a skier weighing 60-65 kg; sledge and passengers weighing up to 150-250 kg). Furthermore, it is possible that the intensity of cargo use was lower compared to 1-3 h per week of training/competition for the racing reindeer in our material. However, as our material did not permit identifying a threshold of loading for detectable changes in the skeleton of working reindeer, we can say that, on average, 3h per week, approximately 4 months a year, and working career of 3 to 4 years work should be evident as increased bone robusticity of long bones, as well as shape changes in humerus and tibia.
According to ethnographic data, cargo reindeer were usually castrated males, but sometimes females were used (Itkonen 1948: 418-419). In modern reindeer racing, only male reindeers (castrated and uncastrated) are used, while draft reindeer in tourism are mainly castrated males (Soppela et al. in press). Thus, a potential source of variation in bone robusticity indicators (via muscle mass) may result from castration status. Testosterone in males result to larger muscles (Schoenau et al. 2000;Ruff 2003), but, in case of reindeer, the lack of testosterone/smaller amount of testosterone due to castration has an opposite effect on muscle mass. As castrated male reindeers do not compete for females during the rut, they pertain the muscle and fat mass over winter. Indeed, some racing reindeer in our sample have been castrated as an attempt to increase bulk (Niinimäki, pers. comm.) and it is a usual practice in Finnish working reindeer management (Soppela et al. in press). As reasons for castration may vary (e.g., to increase bulk, to promote spike growth to sparsely spiked antlers), the age at castration varies greatly in our material, which causes further problems in testing for castration effects. Fortunately, there were castrated males in both working and free-ranging reindeer in our material to balance the activity groups. However, the generalization of our results is limited by the heterogeneous material regarding castration status, but also sex and subspecies. Due to their confounded effect with working status, our material did not permit studying these effects with plausible inference, and our results apply only to the activity status. However, the results suggest that the variation due to body size in our models was able to account for the random variability due to sex or subspecies.

Conclusions
In our current material, the external orientation of humerus allows identification of working and zoo status among male reindeer, and bone shapes of radioulna allow identification of foraging behavior among female reindeer. There was a statistically significant bias towards increased values in bone robusticity variables of long bones in working (male) reindeer compared to male free-ranging/wild reindeer; however, only bone robusticity in tibia was combined with smaller variation among working (male) reindeer compared to male reindeer reference group. In addition, long bone CA values of female zoo reindeer were biased towards greater values compared to female free-ranging/wild reindeer, probably indicative of increased axial and/or muscle loading.
Of bone robusticity variables observed for long bones and metapodials, bone robusticity of tibia seemed to be the most promising skeletal features to identify working status among male reindeer. External dimensions of humerus may be used for identifying foraging behavior and activity status for male reindeer in archaeological material.
Appendix 1 Description of reindeer humerus, radioulna, metacarpal, femur, tibia, and metatarsal morphology, definition of interarticular measurements, and guidelines for bone orientation for scanning

Description of the bone
Reindeer humeral head is relatively flat, due to limited range of motion of the humerus caudally in the shoulder joint. There are extremely prominent ridges in the bone shaft due to muscle attachments, but distal and proximal bone ridges level at mid-shaft provide relatively unaffected observation point for humeral cross-sections.

Interarticular measurement
Interarticular length (Fig. 1) is measured between the most proximal surface of the humeral head and the most distal projection of the medial trochlea.

Positioning humerus for pQCT scanning
Humerus is laid anterior side up, medial side to left (as it would appear in the scout view of the scanner). The bone is oriented to distal articular surface mid-plane parallel to supporting surface (Fig. 2); add support to lateral side as needed to achieve this. This will provide the plane relative to which orientation of the shaft shape is observed. Proximal end of the humerus is elevated as needed to achieve direct perpendicular cross-section of the anterioposterior mid-shaft relative to bone long axis (Fig. 3). Make sure bone shaft long axis is direct perpendicular in anteroposterior, mediolateral, and superoinferior plane at the point of imaged cross-section.

Description of the bone
Greater trochanter is always proximal to most proximal point of femoral head. There are differences in femoral head torsion in the transverse (anteroposterior) plane relative to the bone long axis. Due to variation in AP orientation of the femoral head-neck axis, coronal plane (parallel to supporting surface) for scanning was defined according to distal articular area.

Interarticular measurement
Interarticular length is measured between the most distal projection of medial epicondyle and the most proximal surface of the femoral head (Fig. 4).

Orientation for pQCT scanning
Femur is laid anterior side up, medial side to left (as it would appear in the scout view of the scanner). The bone is oriented to distal articular surface as condyles rest on the supporting surface (Fig. 5). This will provide the plane relative to which orientation of the shaft shape is observed. Proximal end of the femur is elevated as needed to achieve direct perpendicular cross-section of the anterioposterior mid-shaft relative to bone long axis (Fig. 6). Make sure bone shaft long axis is direct perpendicular in anteroposterior, mediolateral, and superoinferior plane at the point of imaged cross-section.

Description of the bone
Ulnar head projects distinctively towards posterior aspect of the bone from proximal articular surface. In addition, ulnar head is angled laterally in the coronal plane relative to bone long axis, and there is interindividual variation in this angle.

Interarticular measurement
Interarticular length (Fig. 7) is measured between the middle of the medial ridge in proximal radial articular surface and the ridge between the articular surfaces for carpal bones in the distal radioulna.

Orientation for pQCT scanning
Radioulna is laid anterior side up, medial side to left (as it would appear in the scout view of the scanner). The bone is oriented to proximal articular surface mid-plane (Fig. 8), parallel to supporting surface. This will provide the plane relative to which orientation of the shaft shape is observed. Distal end of the radioulna is elevated as needed to achieve direct perpendicular cross-section of the anterioposterior mid-shaft relative to bone long axis (Fig. 9). Secure bone firmly in place using molding clay, as projection of the ulnar head need to be below bone long axis and it requires strong hold. Additional support may be needed to the medial side of the distal end of the bone (Fig. 12) to maintain proximal articular mid-plane levelled. Make sure bone shaft long axis is direct perpendicular in anteroposterior, mediolateral, and superoinferior plane at the point of imaged cross-section.

Description of the bone
Reindeer tibia has prominent tibial tuberosity at the proximal shaft, but the bone shaft distally to tibial tuberosity is flat with minimal to no curvature. Rudimentary fibula is located at the posterolateral aspect of the proximal tibia. Bone proximal end angles laterally from bone long axis in the coronal plane, and this angle varies between individuals.

Interarticular measurement
Interarticular length (Fig. 10) is measured between middle of medial proximal articular surface and midpoint of distal articular area, located at mid-ridge in the articular area for talus.

Orientation for pQCT scanning
Tibia is laid anterior side up, medial side to left (as it would appear in the scout view of the scanner). The bone is oriented to proximal articular surface mid-plane parallel to supporting surface (Fig. 11). In many cases, this is achieved without modifications by placing posterior projection of proximal tibial condyles to supporting surface, but there is potential source of variation due to the expression of rudimentary fibula. This will provide the plane relative to which orientation of the shaft shape is observed. Distal end of tibia is elevated, if necessary, to achieve direct perpendicular cross-section of the anteroposterior mid-shaft along the bone long axis (Fig. 12). Additional support may be needed to the medial side of the distal end of the bone to maintain proximal mid-plane levelled. Make sure bone shaft long axis is direct perpendicular in anteroposterior, mediolateral, and superoinferior plane at the point of imaged cross-section.

Metacarpal and metatarsal
Description of the bones Reindeer metacarpal and metatarsal bone distal end consists of two trochlea with trochlear ridge in between medial and lateral trochlea. Anterior aspect of the metacarpal is flat, but in metatarsal, there are two prominences anteromedially and anterolaterally. In the posterior aspect of both bones, there are prominences on the posteromedially and posterolaterally, where the posteromedial prominence is more protruding.

Interarticular measurement
Interarticular length for metacarpal (Fig. 13) and metatarsal (Fig. 14) is measured from the middle of medial proximal articular area (ankle bone and wrist bones) for metacarpal and metatarsal to the most distal point of the medial trochlea on its medial aspect, not from the trochlear ridge for metacarpal and for metatarsal.

Orientation for pQCT scanning
The bone is laid anterior side up, medial side to left (as it would appear in the scout view of the scanner; note the exception for metatarsal). Metacarpal distal articular (Fig. 15) and metatarsal distal articular (Fig. 16) mid-planes are oriented parallel to supporting surface. This will provide the plane relative to which orientation of the shaft shapes are observed. Add support to proximal and/or distal end to achieve direct perpendicular cross-section of the anterioposterior mid-shaft relative to long axis of metacarpal (Fig. 17) and metatarsal (Fig. 18). Make sure bone shaft long axis is direct perpendicular in anteroposterior, mediolateral, and superoinferior plane at the point of imaged cross-section.