Impact of selection and domestication on hindlimb bones of modern reindeer populations: Archaeological implications for early reindeer management by Sámi in Fennoscandia

ABSTRACT For centuries, reindeer herding has been an integral part of the subsistence and culture among the Sámi of northern Fennoscandia. Despite the importance of this husbandry in their history, the timing and details of early reindeer domestication are still highly debated. Indeed, identifying domesticated individuals in the archaeological record remains complicated because reindeer are still considered to be in the early phases of the domestication process. In this work, we propose solutions for identifying domestic individuals using 3D geometric morphometrics on isolated elements from the long bones of the hindlimb in modern reindeer populations. These bones are important for understanding both the mobility of reindeer and the effect of load carrying or draught. A good level of distinction between the size and shape variables of these bones was found among subspecies, sex and lifestyles. This demonstrates that the long bones of the hindlimb can provide information on changes in locomotor behaviour induced by the domestication process, such as control and reduction of reindeer mobility by humans. This also demonstrates that analysis in geometric morphometrics is useful for exploring the use of draught reindeer in early Sámi reindeer herding and the implications for understanding reindeer domestication and early reindeer herding strategies.


Introduction
The domestication of the reindeer (Rangifer tarandus Linnaeus, 1758) had a far-reaching impact on the subsistence, lifeways, economy and cosmology of many peoples in northern Eurasia. Nowadays, reindeer herding is widely practiced in most tundra and taiga areas by around 30 indigenous reindeer herder groups, from eastern Siberia and northern Mongolia to northern Fennoscandia (Mirov 1945;Baskin 2000;Reindeer Herding 2021). The vast majority of the reindeer populations in these areas are now considered domesticated (i.e. reindeer used by humans to perform domestic tasks such as transport or pulling) or semi-domesticated (i.e. reindeer bred freely in near-natural conditions under the supervision of humans) (Syroechkovskii 1995;Baskin 2005). Despite the importance of this husbandry in the history of the Arctic people, the exact period and location of the earliest centres of domestication, as well as the nature of the early reindeer management strategies, are yet to be clearly identified from the archaeological records. Although some scholars estimate the earliest onset of reindeer domestication in the Mesolithic, the earliest material evidence of reindeer management, such as sled runners or harness parts found in Siberia, date from 1,500 BC (Murashkin et al. 2016) and ca. 200 BC-160 AD (Losey et al. 2021). This kind of direct archaeological evidence is lacking in Fennoscandia during the early stages of reindeer husbandry, even though historical and ethnographic sources have suggested that draught reindeer were important to early reindeer domestication (Ingold 1986;Bjørklund 2013).
In Fennoscandia, it is generally supposed that reindeer herding by the Sámi people developed gradually from the Late Iron Age (ca. 800-900 AD) onwards (Helskog and Indrelid 2011;Hansen and Olsen 2014), independently of the indigenous cultures in northwestern Siberia (Røed et al. 2008(Røed et al. , 2011. Historically, the Sámi have preserved their nomadic traditions of hunter-fisher-gatherers, while also practicing small-scale reindeer herding (Bjørklund 2013). During these early phases of husbandry, only small domestic herds of three to four individuals were kept close to human settlements and used as decoys for hunting wild reindeer, for milking or performing various tasks, such as pulling sleds and carrying loads in trade or transhumance (Tegengren 1952;Aronsson 1991;Korhonen 2008). Reindeer herding became the basis of social organisation as recently as from 1,400-1,600 AD onwards, driven by various new socio-economic factors (e.g. colonial policies of the emerging nation states of southern Fennoscandia, intensification of the fur trade or expansion of Christianity), which marked a profound change in Sámi societies (Hansen and Olsen 2014;Salmi et al. 2018). These centuries of socio-economic upheaval resulted in considerable chronological and geographical variations in the adoption and intensification of reindeer herding in Fennoscandia (Tegengren 1952;Bergman et al. 2013;Hedman et al. 2015). Thus, this transition cannot be understood as linear; rather it has been dependent on the effectiveness of the economic and political networks in these territories. Thus, the beginnings of the reindeer domestication in Fennoscandia are still widely discussed and have remained difficult to define, as no abrupt change from hunting to pastoralism has been identified in the archaeological record.
Reindeer bones and teeth collected from archaeological sites can provide direct evidence of this important transition in the socioeconomic structure and the history of the Sámi people and provide a unique opportunity to study the reindeer domestication process.
In fact, domestication is often initiated by controlling the behaviour of wild animals -when they are removed from their natural habitat to be placed in an anthropogenic environment -in order to be used as production or working animals (Price 2002;Vigne 2015;Zeder 2015). These animals are then subjected to selective pressures and different environmental stimuli, which induce significant phenotypic and genotypic changes (Vigne et al. 2005). These resulting morphological changes are traditionally used to document domestication, and they are often defined by the retention of paedomorphic traits (e.g. Morey 1992;Evin et al. 2017;Geiger et al. 2017) or by a general reduction in body size (e.g. Albarella 2002;Evin et al. 2013;Neaux et al. 2020). While these changes do not reflect the early stages of domestication but occur later in the process (Arbuckle 2005), a reduction in body size among reindeer has been observed between domestic and wild individuals (Puputti and Niskanen 2009;Pelletier et al. 2020). However, the morphological changes associated with the early stages of domestication in the response of reindeer populations due to artificial selection induced by the Sámi people are largely unknown and therefore remain to be identified.
Identifying domesticated individuals in the archaeological record remains more difficult due to the presence of two interbreeding reindeer subspecies in Fennoscandia with similar morphologies and overlapping body sizes: 1) the mountain reindeer (R.t. tarandus), including wild populations now only present in southern Norway and the Kola Peninsula in Russia, and domestic herds throughout Lapland; and 2) the wild Finnish forest reindeer (R.t. fennicus), occupying the taiga in central and southern Finland, in the provinces of North Karelia, Savonia and Kainuu, as well as in Russian Karelia (Figure 1). However, the geographical distribution of wild populations and the genetic history of domestic herds must have been significantly different in the past compared to today. Previously, wild mountain reindeer were extant throughout the mountainous areas of northern Fennoscandia, and wild forest reindeer were found throughout the taiga zone of northern Finland (Helle 1982). While domestic reindeer have been domesticated from wild mountain reindeer (Røed et al. 2008), Sámi reindeer herding is more likely to have originated in Scandinavian mountain areas rather than in northeastern Fennoscandia (Wallerström 2000;Bergman et al. 2013;Heino et al. 2021). Also, the extant domesticated reindeer populations in northern Fennoscandia were strongly influenced by the introduction of non-native reindeer from Siberia during the 16th and 17th centuries (Røed et al. 2018). These early domestic reindeer can be identified from wild forms of the species utilising destructive methods such as genetic analysis using ancient DNA (e.g. Bjørnstad et al. 2012;Røed et al. 2018;Salmi and Heino 2019;Heino et al. 2021), geochemical analysis using stable isotopes (e.g. Salmi et al. 2015Salmi et al. , 2020aFjellström et al. 2020), or nondestructive methodologies using the morphometric characteristics of the skeleton. Regarding non-destructive methods, previous studies have succeeded in defining the specific characteristics of each subspecies from cranial elements, body proportions, as well as linear measurements on the postcranial skeleton of modern specimens (Nieminen and Helle 1980;Hakala et al. 1996;Puputti and Niskanen 2008;Puputti and Niskanen 2009). However, these methods remain difficult to apply to the fossil record, in which bones are regularly found to be broken. In addition, the assessments of subspecies on bony elements are prone to error due to the presence of a strong sexual dimorphism in both subspecies and/or a significant phenotypic plasticity that can create significant overlaps (Puputti and Niskanen 2009). Thus, traditional methodologies using measurements and/or bone morphology to identify subspecies and domesticated individuals from archaeological deposits have various limitations and do not allow for a sufficiently robust identification of domestic individuals.
In this study, we partially solved these issues using 3D geometric morphometrics on a large set of long bones of the hindlimb (i.e. femur, tibia and metatarsal) from a broad sample of modern Fennoscandian specimens and compared their discriminate potential. Similar work was recently conducted on the long bones of the forelimb and has shown that they could provide information about changes in feeding and locomotor behaviour prompted by captivity and domestication (Pelletier et al. 2020). Unlike the long bones of the forelimb, the bones of the hindlimb are more subject to the constraints associated with body propulsion and are more impacted by external pressures (Schmidt and Fischer 2009;Hanot et al. 2017;Mallet et al. 2019). In this respect, it is historically known that among most Eurasian nomadic peoples, domestic reindeer may be used for pulling and carrying loads, as well as for riding (Mirov 1945;Nieminen and Pietilä 1998;Inamura 2005;Dwyer and Istomin 2008;Korhonen 2008;Anderson et al. 2017). Thus, these bones have a great potential for approaching reindeer domestication issues, particularly through the control and/or reduction of mobility by humans or the use of reindeer as working animals. The purpose of our study was to provide a reliable method of identifying domestic and wild individuals in modern reindeer populations from Finland, taking into account subspecies, sex and lifestyle, for an application to the archaeological contexts of the indigenous Sámi reindeer herders in northern Fennoscandia. A better comprehension of the details of the domestication process and its implications on human-reindeer relationship could be the key to understanding the history of present and past Sámi communities.

Modern reindeer sample
Several intrinsic (e.g. [sub]species, sex, age, genetic factors) and extrinsic (e.g. geographical position, topography, variations in climate and environment) factors can influence the body size and/or morphology of reindeer (e.g. Thomas and Everson 1981;Weinstock 1997Weinstock , 2002Weladji and Holand 2006;Puputti and Niskanen 2009;Pelletier et al. 2020). In order to minimise these biases, particularly the variation due to genetic diversity and environment, we mainly relied on a reindeer sample from central Finland. Thus, the samples studied included complete or partial skeletons of 123 individuals and represented the two subspecies currently present in Fennoscandia: the mountain reindeer (R.t. tarandus, n = 62) and the wild Finnish forest reindeer (R.t. fennicus, n = 51), as well as first-generation hybrids resulting from the crossing of these two subspecies (n = 10). All specimens were adults with fully fused epiphyses and the sex was known (males, ♂ = 70 and females, ♀ = 53). Knowing that mobility reduction and domestication induce differential stress changes that can affect the shape and robusticity of bones (e.g. Pelletier et al. 2020;Harbers et al. 2020aHarbers et al. , 2020b, lifestyle was taken into account depending on whether the individuals were free-ranging (n = 75), captive (n = 28) or were used for racing and pulling (n = 20). Captive and free-ranging specimens are part of the collection of the Biodiversity Unit and working reindeer are donated skeletons archived at the Laboratory of Archaeology, both hosted by the University of Oulu. The details of samples used for each bone are given in Table 1. Hindlimb long bones were digitised from computed tomography (CT) images. CT scans were performed on a clinical CT scanner (Somatom Definition Flash, Siemens Healthcare, Forcheim, Germany) using 120 kVp, 700 eff. mAs, 0.5 s rotation time, 0.6 mm slice thickness and increment, B70f reconstruction kernel, 140 mm reconstruction diameter, 0.35 pitch and 128 × 0.6 collimation. Each bone was scanned individually to avoid beam-hardening artefacts. Only left bones were selected for digitisation but when left bones were not available, right bones were selected instead and mirrored before analysis (i.e. 8/87 for femur; 10/90 for tibia; and 24/116 for metatarsal).

3D geometric morphometrics
Geometric morphometrics (GMM) is a quantitative approach that allows the comparison of bone shapes and the visualisation of significant morphological changes between groups of specimens by means of spatial coordinates of points called landmarks (Adams et al. 2004;Zelditch et al. 2012). In recent years, this methodology has been widely used to explore the morphofunctional changes induced by the domestication process (Evin et al. 2013(Evin et al. , 2015Owen et al. 2014;Drake et al. 2015;Hanot et al. 2017Haruda 2017;Cucchi et al. 2019;Haruda et al. 2019;Pöllath et al. 2019;Harbers et al. 2020aHarbers et al. , 2020bNeaux et al. 2021;Pelletier et al. 2020). Bone shape was therefore quantified by placing a set of landmarks on the 3D models. The analyses were performed on complete bones but the methodology has been adapted by focusing on the proximal and distal parts unaffected by entheseal changes and pathological lesions Salmi and Niinimäki 2016;Salmi et al. 2020b), as well as on the anatomical parts best preserved in the archaeological record in order to be complementary and directly applicable to the fossil material (Owen et al. 2014;Cornette et al. 2015;Pelletier et al. 2020). Due to the difficulty of quantifying the shape of articular surfaces, trochanters or condyles using traditional landmarks and Unlike landmarks, semilandmarks do not have an exact anatomical correspondence on the structure of the epiphyses, and instead were allowed to slide along curves and surfaces in order to minimise the bending energy of the thin plate spline (TPS) interpolation function (Bookstein 1997;Gunz et al. 2005;Gunz and Mitteroecker 2013). After sliding, all specimen coordinates were aligned using Generalised Procrustes Analysis (GPA, Rohlf and Slice 1990;Bookstein 1991Bookstein , 1996. GPA was conducted using the function 'gpagen' of the 'geomorph' package (v.3.1.1, Adams and Otárola-Castillo 2013) in the Rstudio environment (v.1.1.383, R Development Core Team 2011). All configurations were translated and rotated to minimise the overall sum of the squared distances between the corresponding landmarks and semilandmarks. To remove the effects of scale, GPA also computed a unit centroid size as the square root of the summed squared distances from all landmarks and semilandmarks to their centroid (Bookstein 1996;Dryden and Mardia 1998). Size differences were evaluated from logtransformed centroid sizes for the six epiphyses analysed by pooling the specimens according to 1) subspecies, 2) sexes, 3) lifestyles, 4) 'subspecies + sexes' and 5) 'subspecies + sexes + lifestyles', using Kruskal-Wallis tests with an error threshold set at α = 5%. Pairwise comparisons of the populations were performed using multiple Wilcoxon rank tests according to these Figure 2. Location of anatomical landmarks (yellow spheres), curve sliding (blue spheres) and surface sliding (black spheres) semilandmarks placed on the bone epiphyses. A: proximal femur (from top to bottom: cranial, lateral, caudal, medial and proximal views); B: distal femur (from top to bottom: cranial, lateral, caudal, medial and distal views); C: proximal tibia (from top to bottom: cranial, lateral, caudal, medial and proximal views); D: distal tibia (from top to bottom: cranial, lateral, caudal, medial and distal views); E: proximal metatarsal (from top to bottom: cranial, lateral, caudal, medial and proximal views); F: distal metatarsal (from top to bottom: cranial, lateral, caudal, medial and distal views). The definitions of landmarks and semilandmarks are given in Table 2. Table 2. Definition of anatomical landmarks and semilandmarks for each studied bone epiphyses shown in Figure 2. ALM: anatomical landmarks; CSLM: curve semilandmarks; SSLM: surface semilandmarks. Figure 2

Proximal femur (A on
Most proximo-cranial point of the greater trochanter 12:25 Crest between ALM 2 and ALM 3 (n = 14) 2 Most proximo-caudal point of the greater trochanter 26:37 Crest between ALM 3 and ALM 4 (n = 12) 3 Most cranio-lateral point of the convexity of the greater trochanter 38:46 Crest between ALM 4 and ALM 2 (n = 9) 4 Most disto-caudal point of the greater trochanter 47:55 Crest between ALM 3 and ALM 5 (n = 9) 5 Most distal contact point between the intertrochanteric line and the lateral insertion 56:71 Crest between ALM 2 and ALM 6 (n = 16) crest of the gluteus minimus muscle 72:90 Crest between ALM 7 and ALM 8 (n = 19) 6 Most distal point of the trochanteric fossa 91:108 Crest between ALM 8 and ALM 7 (n = 18) 7 Most lateral point of the border of the head 109:116 Crest between ALM 10 and ALM 11 (n = 8) 8 Most medial point of the border of the head 9 Point of maximum of concavity of the fovea capitis No. SSLM 10 Most proximal point of the lesser trochanter 117:127 Greater trochanter patch (n = 11) 11 Most distal point of the lesser trochanter 128:144 Femoral head patch (n = 17) Distal femur (B on Figure 2 Most proximal point of the medial lip of the trochlea 13:32 Crest between ALM 1 and ALM 4 (n = 20) 2 Most proximal point of the trochlear groove 33:50 Line between ALM 2 and ALM 5 (n = 18) 3 Most proximal point of the lateral lip of the trochlea 51:70 Crest between ALM 3 and ALM 6 (n = 20) 4 Most distal point of the medial lip of the trochlea 71:86 Crest between ALM 7 and ALM 10 (n = 16) 5 Distal maximum of curvature of the trochlear groove 87:108 Crest between ALM 10 and ALM 9 (n = 22) 6 Most distal point of the lateral lip of the trochlea 109:126 Crest between ALM 12 and ALM 11 (n = 18) 7 Most medial point of the fossa extensoria 127:142 Crest between ALM 11 and ALM 12 (n = 16) 8 Most cranial point of the fossa extensoria 9 Most lateral point of the fossa extensoria No. SSLM 10 Contact point between the intercondylar line and the lateral condyle 143:156 Trochlea patch (n = 14) 11 Contact point between the intercondylar line and the medial condyle 157:170 Medial condyle patch (n = 14) 12 Most cranial point of the medial condyle 171:187 Lateral condyle patch (n = 17) Proximal tibia (C on Figure 2 Most cranial point of the articular surface of the lateral condyle 12:21 Crest between ALM 1 and ALM 2 (n = 10) 2 Most caudal point of the lateral tubercle of the intercondylar eminence 22:27 Crest between ALM 2 and ALM 3 (n = 6) 3 Most caudal point of the articular surface of the lateral condyle 28:41 Crest between ALM 3 and ALM 4 (n = 14) 4 Most crano-lateral point of the articular surface of the lateral condyle 42:44 Crest between ALM 4 and ALM 1 (n = 3) 5 Most cranial point of the articular surface of the medial condyle 45:56 Crest between ALM 5 and ALM 6 (n = 12) 6 Most caudo-lateral point of the articular surface of the medial condyle 57:66 Crest between ALM 6 and ALM 7 (n = 10) 7 Most proximal point of the medial tubercle of the intercondylar eminence 67:71 Crest between ALM 7 and ALM 5 (n = 5) 8 Most caudal point of the medial condyle 72:82 Crest between ALM 9 and ALM 10 (n = 11) 9 Most proximal point of the lateral part of the tibial tuberosity No. SSLM 10 Most distal point of the lateral part of the tibial tuberosity 83:96 Lateral condyle patch (n = 14) 11 Most cranio-distal point of the lateral condyle 97:105 Medial condyle patch (n = 9) Distal tibia (D on Figure 2 Most distal point of the contact between the medial malleolus and the distal 7:16 Crest between ALM 1 and ALM 2 (n = 10) articular surface 17:24 Crest between ALM 2 and ALM 3 (n = 8) 2 Most distal point of the medial part of the distal articular surface 25:27 Crest between ALM 3 and ALM 4 (n = 3) 3 Most distal point of the contact between the lateral malleolus and the distal 28:36 Crest between ALM 3 and ALM 5 (n = 9) articular surface 37:43 Crest between ALM 6 and ALM 1 (n = 7) 4 Most cranio-distal point of the lateral malleolar sulcus 5 Most cranio-distal point of the lateral part of the distal articular surface No. SSLM 6 Most cranio-distal point of the medial part of the distal articular surface 44:62 Distal articular surface patch (n = 19) Proximal metatarsal (E on Figure 2 Most anterior point of contact between the articular surface for the naviculo-11:18 Crest between ALM 1 and ALM 2 (n = 8) cuboid bone and the articular surface for the lateral cuneiform bone 19:22 Crest between ALM 2 and ALM 3 (n = 4) 2 Most medial point of the articular surface for the lateral cuneiform bone 23:25 Crest between ALM 4 and ALM 5 (n = 3) 3 Most postero-lateral point of the articular surface for the lateral cuneiform bone 26:34 Crest between ALM 5 and ALM 1 n = 9) 4 Most postero-medial point of the articular surface for the naviculo-cuboid bone 35:39 Crest between ALM 1 and ALM 4 (n = 5) 5 Most postero-lateral point of the articular surface for the naviculo-cuboid bone 40:41 Crest between ALM 7 and ALM 8 (n = 2) 6 Top of the most postero-lateral point of the proximal part of the metatarsal bone 42:43 Crest between ALM 8 and ALM 7 (n = 2) 7 Most lateral point of the posterior articular surface for the naviculo-cuboid bone 44:45 Crest between ALM 9 and ALM 10 (n = 2) 8 Most medial point of the posterior articular surface for the naviculo-cuboid bone 46:47 Crest between ALM 10 and ALM 9(n = 2) 9 Most posterior point of the articular surface for the medial cuneiform bone Most proximal point of the anterior outline of the medial articular eminence 23:36 Crest between ALM 2 and ALM 12 (n = 14) 3 Most latero-proximal point of the anterior outline of the medial articular eminence 37:46 Crest between ALM 3 and ALM 11 (n = 10) 4 Most medio-proximal point of the anterior outline of the lateral articular eminence (LAE) 47:56 Crest between ALM 4 and ALM 10 (n = 10) 5 Most proximal point of the anterior outline of the lateral articular eminence 57:70 Crest between ALM 5 and ALM 9 (n = 14) 6 Most latero-proximal point of the anterior outline of the lateral articular eminence 71:78 Crest between ALM 6 and ALM 8 (n = 8) 7 Most latero-proximal point of the posterior outline of the lateral articular eminence (Continued) different categories. To control for the false discovery rate, a multi-comparison correction was applied to the P-values using the 'Benjamini-Hochberg' method (Benjamini and Hochberg 1995). This test was run using the function 'p.adjust' in the 'stats' package in Rstudio. Shape differences between these different groups were estimated using a multivariate analysis of variance (MANOVA) performed on the Procrustes coordinates, with significant interaction (α = 5%) assumed to reflect group differences. Shape variation was visualised using Principal Component Analysis (PCA) based on Procrustes coordinates and via the function 'gm.prcomp' in the 'geomorph' package. To better apprehend variations along the principal axes, we created a 3D digital mesh for each of the elements that were warped towards the Procrustes grand mean using a thin plate spline (TPS) interpolation function (Bookstein 1991). The visualisations of shapes at the extremes of the principal component axes were performed from the surface of the Procrustes mean configuration (Wiley et al. 2005), with magnification by a scale factor of 0.1, using the function 'tps3d' in the 'Morpho' package (details of the process are provided in the documentation of the package v.2.8; Schlager 2020). Finally, allometry was assessed using multivariate regressions of shape variables on the log-transformed centroid sizes using the function 'procD.lm' in the 'geomorph' package (Adams and Otárola-Castillo 2013).

Size variation of skeletal elements
The results of the Kruskal-Wallis tests on log-transformed centroid size data reveal significant differences between all categories, namely, subspecies, sexes, lifestyles, 'subspecies + sexes', as well as 'subspecies + sexes + lifestyles' (P < 0.01; Table 3). All six analysed epiphyses of the hindlimb long bones displayed the same pattern of size differentiation among subspecies, sexes and lifestyles (Supplementary Figs. S1, S2, S3, S4). The pairwise comparisons revealed that forest reindeer (R.t. fennicus) were systematically larger than mountain reindeer (R.t. tarandus) and hybrids (all P < 0.01), but these latter two did not differ significantly from each other (all P > 0.05; Supplementary Fig.  S1). Regarding sexes, males were significantly larger than females (all P < 0.01; Supplementary Fig. S2 Supplementary Fig. S4).
Other particularities could be observed when the specimens were divided into 11 groups, simultaneously comprising the subspecies, sex and lifestyle ( Figure 3, Table 4). Significant differences were always found between the free-ranging male R.t. fennicuswith the largest bone size -and all other groups (all P < 0.01), except between the only captive male R.t. fennicus and the two captive male R.t. tarandus. Contrariwise, the free-ranging female R.t. tarandus had a smaller bone size and were generally  Most proximal point of the posterior outline of the medial articular eminence 87:95 Lateral trochlear patch on MAE (n = 9) 12 Most medio-proximal point of the posterior outline of the medial articular eminence 96:104 Medial trochlear patch on LAE (n = 9) 13 Most proximo-medial point of the ligament insertion fossa 105:112 Lateral trochlear patch on LAE (n = 8) 14 Most proximo-lateral point of the ligament insertion fossa

Shape variation of skeletal elements
With the exception of proximal epiphysis of the femur and distal epiphysis of the tibia, MANOVA analyses revealed at least a significant difference in shape among some groups for all the other elements, although it varied according to the bone or category investigated (Table 5). In contrast to what we have previously noted concerning size differences of the bone elements, shape did not differentiate significantly between the 'subspecies', 'sexes' and 'lifestyles' categories (except for the proximal metatarsal). However, significant differences were noted for the 'subspecies + sex' category for the distal femur and metatarsal. Furthermore, significant differences were found for the proximal parts of the tibia and metatarsal, as well as the distal parts of the femur and metatarsal among the 'subspecies + sexes + lifestyles' category.
• Femur: For the proximal femur epiphysis, the first two axes of the PCA expressed 27.85% of the global variance and did not appear to show preferential variations according to the different groups ( Figure 4A). Despite these significant overlaps, it would appear that females, all categories included (i.e. free-ranging and captive in both subspecies, as well as hybrids) and free-ranging male R.t. tarandus, were more distributed on the positive values of the PC1. This showed a more slender and thin aspect with a rounded femoral head and a shorter and thinner neck. In addition, the greater trochanter convexity was less developed laterodistally and the top was less pronounced. At the opposite, morphological variation along the negative values of the PC1 showed a more massive morphology, where more male specimens were located. In particular, this was expressed by a large femoral head and neck. The greater trochanter convexity was more expanded laterodistally and the top of the greater trochanter was more developed. However, the variations along axis 2 did not allow us to observe preferential group trends. In contrast, the variations for the distal femur epiphysis were more marked between groups and expressed both sexual and subspecific variations, but also lifestyles in the morphospace on the PC1 ( Figure 4B).  more medially orientated and craniocaudally broader, while the trochlear lips were separated by a shallower trochlear groove and the femoral medial patellar margin presents a more elliptical shape. The lateral and medial condyles were more mediolaterally oblique and separated by a wider intercondylar space. Within the subspecies and hybrids, females and captive individuals were distributed more towards the negative values of the PC1 compared to males and free-ranging individuals, respectively. However, like the proximal epiphysis, the variations along the PC2 did not allow for preferential group trends to be observed. • Tibia: The proximal tibia epiphysis generally followed the same variations along PC1 as the distal femur epiphysis for each group ( Figure 5A) The theoretical shape at the PC1 maximum showed a massive morphology with a mediolaterally and craniocaudally broad epiphysis. The central intercondylar area was reduced and the extensor groove was deeper. The lateral and medial articular surface areas were roughly similar with a caudally extended sliding surface for the m. popliteus. The tibial tuberosity was also rounded and elongated. In contrast, negative values of the PC1 displayed a relatively gracile morphology with a relatively large central intercondylar area and a shallow extensor groove. The lateral articular surface was larger than the medial surface, while the tibial tuberosity was slightly mediolaterally broader. The PC2 (11.40% of the global variance) instead showed a sexual distinction, in which males led the variation more towards the positive value of the axis, and females towards the negative values (except for the captive female R. t. fennicus). This resulted in a more slender morphology on the negative values of PC2, and more massive on the positive values. Thus, in females, the proximal epiphysis appeared to be mediolaterally narrower and had a tibial tuberosity that deviated medially. In contrast, variations in the distal tibia epiphysis were expressed differently in the morphospace on both the PC1 and PC2 ( Figure 5B). Captive individuals in both subspecies appeared to be distribute more towards the negative value of the PC1 (15.90% of the global variance), in the same way in males and females, compared to their freeranging counterparts. The negative values of the PC1 expressed more of a mediolaterally broader distal section and a less prominent medial malleolus. The distal articular surface was distally more stretched on the caudal and cranial edges, respectively. On the positive values of the PC1, the morphology was squarer with a prominent medial malleolus and a less distally pronounced distal articular surface on the caudal edge. Despite some overlaps between groups, the PC2 (9.69%) roughly separated R.  Figure 6A). The theoretical shape at the PC1 maximum showed a massive morphology with a craniocaudally broad epiphysis, but the space between the articular surface for the lateral cuneiform bone and the articular surface for the naviculo-cuboid bone was narrower. Thus, the posterior articular surface for the naviculo-cuboid bone was smaller than the articular surface for the medial cuneiform bone. In addition, the postero-lateral part of the proximal epiphysis was markedly more developed proximally. At the opposite, negative values indicated a more slender morphology, less craniocaudally stretched. The space between the articular surface for the lateral cuneiform bone and the articular surface for the naviculo-cuboid bone was broader and the postero-lateral part of the proximal epiphysis was markedly less developed. Also, the posterior articular surface for the naviculo-cuboid bone was more extensive than the articular surface for the medial cuneiform bone. Female and captive individuals tended to cluster more towards the positive values of the PC2 (9.61% of the global variance). This resulted in a flattening of the proximal articular surfaces, which were also more mediolaterally extended.
Regarding the distal epiphysis of the metatarsal, the PCA also showed some minor distinction between groups firstly more distributed R.t. tarandus on the negative values of the PC1, and R.t. fennicus on the positive values of the PC1 (20.65% of the global variance; Figure 6B). The epiphysis was more compressed proximodistally and more stretched craniocaudally in R.t. tarandus, while it was more stretched proximodistally and less craniocaudally in R.t. fennicus. This peculiarity was particularly marked in the captive female R.t. fennicus. Along negative values of the PC2 (8.94% of the global variance), in which there were more captive individuals, the medial and lateral edges of the articular eminences were thinner along the mediolateral axis.

Allometry
For all epiphyses, allometry was significant (all P < 0.01). Although the percentage of shape variance related to size was always relatively low (proximal femur: 2.94%; distal femur: 11.07%; proximal tibia: Figure 6. Scatter plots of the two first axes (PC 1 and PC 2) of principal component analyses performed on the shape data, and visualisation of shape variation via deformation of the mean shape along negative and positive PC1 and PC2 values (magnified by a scale factor of 0.1). A: proximal metatarsal; B: distal metatarsal. The proportion of the total variance expressed by the axes PC1 and PC2, respectively, is indicated in brackets.
6.38%; distal tibia: 3.54%; proximal metatarsal: 2.50%; distal metatarsal: 3.99%), this also indicates that the allometric pattern varies slightly depending on the bone. Multivariate regressions of shape scores against log-transformed centroid size showed that the freeranging male R.t. fennicus had the largest centroid size and was quite distinguishable from the other groups (Figure 7). At the opposite, the captive and free-ranging female R.t. tarandus and captive hybrids had the smallest centroid size. Generally, captive individuals had either lower shape scores or a lower centroid size than free-ranging individuals, as well as in males compared to females, whether in R.t. tarandus, R.t. fennicus or hybrids. The separation between these groups was similar and consistent with previous analyses (i.e. size and shape), although controlling the allometry, morphological variations showed this to be less evident between intermediate groups. The overlaps mainly concerned male R.t. tarandus for all lifestyles combined (i.e. captive, free-ranging and working), female R.t. fennicus (i.e. captive and free-ranging) and male hybrids. However, for similar shape scores, particularly for femur (proximal and distal) and proximal tibia, working male R. t. tarandus showed the highest values in centroid size compared to the free-ranging and captive female R.t. fennicus, as well as the freeranging and captive male R.t. tarandus and hybrids. In contrast for a given size, particularly for the distal tibia and metatarsal (proximal and distal), working male R.t. tarandus generally had lower shape scores than female R.t. fennicus. Finally, clusters of the different groups overlapped less for the femur and the tibia than for the metatarsal.

Importance of subspecies and sex for identifying early domesticated reindeer in Fennoscandia
In this study, we showed that R.t. fennicus have significantly larger bones than R.t. tarandus, which is in accordance with previous studies (Puputti and Niskanen 2009;Pelletier et al. 2020). Our analyses of morphological variation also allowed a relatively good distinction between the two reindeer subspecies, although this distinction appears to be more evident on the femur and tibia than on the metatarsal. For example, the patellar surface of the distal femur presents a more circular section in R.t. fennicus, while it is more elliptical in R.t. tarandus. For the lateral tibial plateau, the morphology is more round in proximal view in R.t. fennicus and more compressed in the medial-lateral direction in R.t. tarandus. This is in accordance with the work of Curran (2012Curran ( , 2015 which reveal adaptations to closed (i.e. taiga-type in the case of R.t. fennicus) or open (i.e. tundra-type in the case of R.t. tarandus) environments. This indicates a relatively strong phylogenetic signal on the shape of the long bones of the hindlimb related to locomotion and adaptation to their environment. However, the few morphological differences between the subspecies for the metatarsal would indicate that this bone probably does not respond similarly to frequency and magnitude of function compared to more proximal limb elements. Indeed, their overall morphology is more likely to reflect a greater involvement in the support of the limbs and therefore of the body mass (Niinimäki et al. in press; Pelletier et al. submitted). It has already been hypothesised that these morphological differences on the long limb bones between both subspecies could be due to behavioural and ecological differences (Pelletier et al. 2020). Indeed, R.t. tarandus is a more gregarious reindeer living in more open tundra areas to the north, while R.t. fennicus has a more complex social organisation in a more closed taiga environment to the south. In this respect, R.t. fennicus is significantly better adapted to taiga conditions with a deep and soft snow coverwhich explains the greater withers height and leg length -while R.t. tarandus lives on hard-packed tundra snow (Nieminen and Helle 1980). Thus, the identification of the subspecies at archaeological sites appears to be an essential prerequisite for understanding the history of the past Sámi communities, since this may reflect different subsistence strategies or cultural interpretations. Indeed, the presence of R.t. fennicus in a deposit could directly induce a subsistence strategy based on hunting, since this subspecies has never visibly been domesticated in these regions, while the presence of R.t. tarandus could involve herding (e.g. Salmi et al. 2018;Salmi and Heino 2019;Heino et al. 2021). In addition, since the two subspecies also have different eco-ethologies, this could affect the hunting strategies of past human societies, including mass hunting in the case of the wild mountain reindeer R.t. tarandus or individual hunting in the case of the wild forest reindeer R.t. fennicus (Rankama and Ukkonen 2001). However, the two reindeer subspecies now share ecological similarities in some Fennoscandian regions, which could make it difficult to accurately assess the environmental effect on bone shape. Although they appear to have had more marked biotopes in the past than today (Helle 1982), this should be evaluated with the utmost caution since it is historically known that the ranges of wild and domestic reindeer have fluctuated greatly under the pressure of various anthropogenic and/or climatic factors (Ingold 1980;Reimers and Colman 2006;Pape and Löffler 2012;Bergman et al. 2013).
In addition to phylogeny, sex also leads a significant part of variations in the shape of the bony elements of the hindlimb. Our results confirmed that male reindeer had significantly larger bone elements than female reindeer, which is explained by a strong sexual dimorphism (Reimers et al. 1983;Weinstock 2000;Puputti and Niskanen 2009;Melnycky et al. 2013;Pelletier et al. 2020). This difference was also significant within subspecies, in which the largest body size was found in male R.t. fennicus, then in male R.t. tarandus, female R.t. fennicus and finally in female R.t. tarandus, which had the smallest bone dimensions. In terms of morphology, this resulted in more massive bones in male individuals and slender bones in females, which is probably due to the different weight-bearing functions of the skeletal elements. In general, identifying sex is relatively easy and is often applied in archaeology to establish sex ratios and identify the hunting methods and/or reindeer exploitation by past Sámi societies (e.g. Puputti and Niskanen 2009;Salmi et al. 2015Salmi et al. , 2021. In fact, this has an even greater implication for finding domestic individuals in the fossil record since domestic males R.t. tarandus are preferentially used for transport and pulling (Korhonen 2008;Salmi and Niinimäki 2016). This is because males are larger and more robust than females, which makes them more capable of carrying out domestic tasks, although female R.t. tarandus could also be present in domestic Sámi herds, especially for milking purposes (Tegengren 1952;Aronsson 1991). However, the composition and structure of the herds has changed over time, particularly the transition to meat production-based reindeer herding. Thus, the highest possible proportion of breeding-age female reindeer became optimal, in contrast to a minimum number of males to serve the females successfully during the rut (Holand 2007).
Another factor that could lead to the misidentification of subspecies and/or sex in archaeological assemblages is the possibility of finding hybrid individuals that are the result of crossbreeding between R.t. fennicus and R.t. tarandus. Since domestic reindeer are mostly left to free range by Sámi herders, this hybridisation between the two subspecies could then occur in areas in which their distribution ranges overlap (Nieminen and Helle 1980;Nieminen and Ojutkangas 1986;Røed et al. 2008Røed et al. , 2011. It is also highly likely that wild individuals have been captured to incorporate them into domestic herds in order to avoid consanguinity (Sommerseth 2011). The main problem of their potential presence in archaeological contexts could be the potentially large overlap between hybrids and their parents in terms of morphometric diversity. Indeed, hybrids can present morphological traits that are more similar to a particular parent, as well as an intermediate morphology and size (Evin et al. 2015;Hanot et al. 2017Hanot et al. , 2019Hanot and Bochaton 2018;Savriama et al. 2018 This resulted in more slender epiphyses, although it appeared to be significantly more evident in the distal femur and proximal tibia. The inclusion of hybrid individuals in our study revealed significant overlaps in both size and bone shape variations, particularly in the case of intermediate specimens such as male hybrids and wild female R.t. fennicus. Thus, their presence on an archaeological site could lead to further confusion in the interpretation of fossil material.

Impact of selection and archaeological perspectives for the identification of early domesticated reindeer
Beyond the impact of phylogeny, ecology and sex on morphometric variation among Finnish modern reindeer populations, there is also a significant impact of lifestyle on bone morphology. Despite the small number of captive male individuals in our sample, individuals bred in captivity generally had smaller bones compared to freeranging reindeer. This reduction in body size of captive individuals relative to their wild counterparts is a typical characteristic of the domestication syndrome and has been long observed in many other domestic species (e.g. Davis 1981;Morey 1992;Dayan 1994;Zohary et al. 1998;Zeder et al. 2006;Rowley-Conwy et al. 2012). In reindeer, it had already been noted that domesticated individuals were smaller than wild ones (Puputti and Niskanen 2009;Pelletier et al. 2020). However, it is important to note that most captive individuals in our sample were born in the wild and do not have a long ancestry in zoo (Pudas, personal communication). This would therefore imply that the effects of a reduction in mobility are observed forthwith and are evident without even a preferential selection. This observation appears to be consistent with recent experimental work carried out on pigs and wild boars (Harbers et al. 2020a;Neaux et al., 2020). Body size reduction in captivity appears to be particularly more evident for wild forest reindeer (R.t.fennicus) than for R.t. tarandus. This could be explained by the fact that there are no longer any completely wild modern R.t. tarandus genetic lineages in Finland following the introgression of domestic reindeer into the wild gene pool in the 19th century (Røed et al. 2011(Røed et al. , 2014 and/or resulting from a strong selection of individuals in breeding for several centuries. Thus, the impact of captivity on these R.t. tarandus individuals was not perceptible in our sample, but could have occurred during the early stages of domestication, as we observed in R.t.fennicus. Captivity also resulted in significant changes in terms of bone morphology, although this was observed more on the femur and tibia than on the metatarsal. Despite some overlaps with their freeranging counterparts, captive female R.t. tarandus differed relatively well from the other categories. Indeed, the bony elements had a smaller size and a generally more slender shape than male R.t. tarandus or R.t. fennicus. In contrast, wild male R.t. fennicus were also relatively easy to identify since their bones were larger and had a more robust shape. The fact that female R.t. tarandus and wild male R.t. fennicus were potentially the easiest to identify from the archaeological record directly suggests opposing socio-economic strategies among past Sámi groups (e.g. Salmi et al. 2018;Salmi and Heino 2019;Heino et al. 2021). As in the early phases of husbandry, these groups kept small domestic herds close to settlements and the identification of small-size female R.t. tarandus could involve herding since females were used for milking (Tegengren 1952). Conversely, the identification of male R.t. fennicus involves wild reindeer hunting and therefore a diametrically opposed human-reindeer relationship.
On the other hand, the identification of male R.t. tarandus proved to be more difficult. However, this identification is of major interest for the application to the fossil record since these reindeer were preferentially used by the Sámi for transport and pulling (Korhonen 2008;Salmi and Niinimäki 2016). Except for the femur and metatarsal distal epiphyses, no significant differences in size were found between free-ranging, captives and working male individuals. Nevertheless, working individuals tended to be larger than free-ranging individuals. This could be explained by the fact that working reindeer are selected for their physical properties and ability to perform domestic work, although other characteristics such as personality and behaviour can be considered to be the most important criteria (Soppela et al. 2020). Thus, size reduction should not be a reliable criterion for identifying male domestic R.t. tarandus. In addition, we did not note any significant differences in shape between free-ranging and working R.t. tarandus. We believe that this could be explained by the fact that, apart from their use for domestic tasks (e.g. pulling sleds, racing), working reindeer are to left free range, just like wild individuals. Activity levels can therefore be similar most of the time, leading to similar size and robustness requirements on hindlimb long bones.
However, there is a noticeable difference between the largest individuals with more robust and generally more enlarged epiphyses, and captive individuals whose reduced activity and smaller body size result in a more slender morphology. For the proximal femur, the epiphysis was more developed mediolaterally in larger male reindeer, particularly the free-ranging R.t. fennicus and most of the working R.t. tarandus. This could indicate a stronger resistance to constraints linked to both body propulsion and weight bearing (Kappelman 1988;Mallet et al. 2019 and references therein). At the same time, the distolateral development of the greater trochanter partially enhances a large lever arm for the muscles attached to the proximal femur, allowing strong hip flexion and abduction. In these larger individuals, the improvement in weight bearing and the amplitude of flexion are favoured by a more deviated position of the knee. This is expressed by a greater lateral torsion of the rotation axis of the trochlea for the distal epiphysis of the femur, as well as a more medially widened tibial plateau. In addition, the enlargement of the tibial tuberosity involves the presence of a stronger and larger patellar ligament, which strengthens the knee joint (Mallet et al. 2019). This general conformation, mainly of the femur and tibia, provides greater stability to the knee joint, in particular among cursorial cervids (Hildebrand 1985). This would be due to repetitive flexion of the knee articulation, more involved during propulsion (Kappelman 1988;Curran 2012). In R.t. fennicus this may be partially explained by body mass and adaptation to the closed environments of the taiga, while in working reindeer it could imply a need for a large range of motion in the hip and knee joint (Salmi et al. 2020b). However, this activity is reduced or even absent in captive individuals. In addition, the distal tibia and distal metatarsal epiphyses appeared to widen mediolaterally among captive individuals. This widening of more caudal limb elements and their distal ends could be a result of the need to strengthen articular areas for prolonged periods of static loading, as has been observed on the forelimb (Pelletier et al. 2020).
Thus, the reduction in body size of wild individuals in captivity could be the first element to consider in identifying domestic individuals in the archaeological record. However, this does not concern all individuals that comprise the herds kept by the herders, as most of them are left to free range. The focus must be on particular individuals, namely, some domestic female R.t. tarandus that could be used for milking and kept close to settlements, or larger male R.t. tarandus that are instead chosen for domestic tasks such as transport or pulling. In addition, the control of mobility combined with their use to perform domestic tasks induces significant changes in the bone morphology of reindeer. As the activity of captive individuals is greatly reduced, general morphology would trend towards greater gracility than their freeranging counterparts. In contrast, the bones of the hindlimb appear to be more adapted to propulsion and weight bearing in working reindeer. However, to our knowledge, there is no evidence to suggest that Fennoscandian reindeer were kept in total captivity from the Iron Age. These particular individuals had to be tamer and kept close to settlements under fairly close supervision by Sámi herders, rather than kept corralled. In any event, these small domestic herds were probably engaging in lower levels of physical activity compared to free-ranging animals. Also, the presence of wild individuals identified from fossil bones, such as wild R.t. fennicus, is not necessarily evidence of the absence of domestication by Sámi. Indeed, Sámi groups had long continued to hunt wild reindeer along with the breeding of domestic individuals (Bjørklund 2013;Hansen and Olsen 2014). For example, in Sámi sacrificial sites, wild and/or domestic reindeer bones were frequently deposited and could reflect cultural changes within these groups (Salmi et al. 2015(Salmi et al. , 2018(Salmi et al. , 2020aSalmi and Heino 2019;Heino et al. 2021). Finally, the domestication process is even more difficult to understand since it has been gradual and does not appear to have been synchronous in the different regions, nor with the same amplitude (Tegengren 1952;Lundmark 2007;Korhonen 2008;Bjørklund 2013). A careful analysis of the size, shape and allometry of reindeer bones, as well as the archaeological context and associated cultural material, are all essential parameters to be taken into consideration in order to better understand early reindeer management by the Sámi in Fennoscandia.

Conclusion
Historically, reindeer have probably been one of the most recently domesticated species by humans. In this respect, many scholars consider them to still be in the early stages of the domestication and could serve as an excellent model species to understand how the initial steps of this process may have occurred. Thus, the identification of relevant biomarkers on the skeleton appears to be an essential prerequisite for documenting the origin of this process. Our study has demonstrated the potential of 3D GMM studies to identify subspecies, sex and lifestyle among modern reindeer populations from Finland. Our results showed significant changes in the size and shape of most of the isolated elements of the hindlimb, allowing a relatively reliable distinction between wild and domesticated individuals. Captivity, which results in less activity, leads to reduced body size and a more thin and slender bone morphology in domesticated individuals than their wild counterparts. In contrast, working reindeer, which are specially selected for their aptitude to perform domestic tasks (e.g. transport and pulling), are larger and more robust; their bones are more adapted to the constraints linked to body propulsion and weight bearing. These plastic changes associated with selection and domestication can be used as a proxy for the early process of reindeer management in the archaeological record and can therefore to shed light on the evolution of socio-economic models of the different Sámi communities of reindeer herders in Fennoscandia. However, caution must be exercised with regard to the correct identification of domestic reindeer due to a domestication process that took place gradually and in a nonsynchronous manner in the different regions, and which did not have the same amplitude. In addition, Sámi populations have simultaneously practiced wild reindeer hunting and small-scale reindeer herding for several centuries. Every parameter such as size, shape and allometry must be thoroughly analysed and coupled with archaeological contexts in order to be able to identify individuals and better understand the morphometric variability of reindeer in Fennoscandia. New studies allowing for a better understanding of the morphometric diversity of reindeer should be carried out in the future (e.g. cross-sections, teeth), complemented by analyses of ancient DNA, stable isotopes, morpho-functional investigations and entheseal changes and pathological lesions. Such studies would allow for the refinement of research on archaeological sites in order to better identify early reindeer domestication and herding practices in time and space.