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.

KEYWORDS:

• 3d geometric morphometrics
• captivity
• domestication signal
• draught animals
• Rangifer tarandus
• Sámi archaeology

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, 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 socio-economic 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. 2015, 2020a; Fjellström et al. 2020), or non-destructive 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.

Figure 1. Current geographical distribution of the two reindeer subspecies, including wild and domestic populations in Fennoscandia (from Suominen and Olofsson 2000).

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.

Material and methods

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 1997, 2002; Weladji 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. 2020a, 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).

Table 1. Detail of specimens studied from the Zoological Museum of Oulu according to the subspecies (R.t. tarandus, R.t. fennicus and hybrid), sex (male = ♂ and female = ♀) and lifestyle (free-ranging, captive and working)

Proximal femur	Free-ranging	Captive	Working	TOTAL
R.t. tarandus	8 ♂/5 ♀	2 ♂/9 ♀	9 ♂/0 ♀	19 ♂/14 ♀
R.t. fennicus	19 ♂/20 ♀	1 ♂/4 ♀	0 ♂/0 ♀	20 ♂/24 ♀
Hybrids	0 ♂/0 ♀	6 ♂/4 ♀	0 ♂/0 ♀	6 ♂/4 ♀
TOTAL	27 ♂/25 ♀	9 ♂/17 ♀	9 ♂/0 ♀	45 ♂/42 ♀
Distal femur
R.t. tarandus	8 ♂/5 ♀	2 ♂/9 ♀	9 ♂/0 ♀	19 ♂/19 ♀
R.t. fennicus	19 ♂/20 ♀	1 ♂/4 ♀	0 ♂/0 ♀	20 ♂/26 ♀
Hybrids	0 ♂/0 ♀	6 ♂/4 ♀	0 ♂/0 ♀	6 ♂/4 ♀
TOTAL	27 ♂/25 ♀	9 ♂/17 ♀	9 ♂/0 ♀	45 ♂/42 ♀
Proximal tibia
R.t. tarandus	7 ♂/6 ♀	2 ♂/8 ♀	9 ♂/0 ♀	18 ♂/14 ♀
R.t. fennicus	21 ♂/20 ♀	0 ♂/5 ♀	0 ♂/0 ♀	21 ♂/25 ♀
Hybrids	0 ♂/0 ♀	6 ♂/4 ♀	0 ♂/0 ♀	6 ♂/4 ♀
TOTAL	28 ♂/26 ♀	8 ♂/17 ♀	9 ♂/0 ♀	45 ♂/43 ♀
Distal tibia
R.t. tarandus	7 ♂/6 ♀	2 ♂/8 ♀	9 ♂/0 ♀	18 ♂/14 ♀
R.t. fennicus	21 ♂/21 ♀	1 ♂/5 ♀	0 ♂/0 ♀	22 ♂/26 ♀
Hybrids	0 ♂/0 ♀	6 ♂/4 ♀	0 ♂/0 ♀	6 ♂/4 ♀
TOTAL	28 ♂/27 ♀	9 ♂/17 ♀	9 ♂/0 ♀	46 ♂/44 ♀
Proximal metatarsal
R.t. tarandus	18 ♂/11 ♀	2 ♂/8 ♀	15 ♂/0 ♀	35 ♂/19 ♀
R.t. fennicus	15 ♂/19 ♀	1 ♂/6 ♀	0 ♂/0 ♀	16 ♂/25 ♀
Hybrids	0 ♂/0 ♀	6 ♂/4 ♀	0 ♂/0 ♀	6 ♂/4 ♀
TOTAL	33 ♂/30 ♀	9 ♂/18 ♀	15 ♂/0 ♀	57 ♂/48 ♀
Distal metatarsal
R.t. tarandus	18 ♂/13 ♀	2 ♂/8 ♀	15 ♂/0 ♀	35 ♂/21 ♀
R.t. fennicus	22 ♂/21 ♀	1 ♂/6 ♀	0 ♂/0 ♀	23 ♂/27 ♀
Hybrids	0 ♂/0 ♀	6 ♂/4 ♀	0 ♂/0 ♀	6 ♂/4 ♀
TOTAL	40 ♂/34 ♀	9 ♂/18 ♀	15 ♂/0 ♀	64 ♂/52 ♀

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 morpho-functional changes induced by the domestication process (Evin et al. 2013, 2015, 2017; Owen et al. 2014; Drake et al. 2015; Hanot et al. 2017, 2018; Haruda 2017; Cucchi et al. 2019; Haruda et al. 2019; Pöllath et al. 2019; Harbers et al. 2020a, 2020b; Neaux 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 (Niinimäki and Salmi 2016; 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 the lack of homologous anatomical structures, semilandmarks have been included on curves and surfaces to help capture the three-dimensional structure of the epiphyses (Bookstein 1997). Thereby, 3D templates with 64 anatomical landmarks (ALM), 444 curve semilandmarks (CSLM) and 162 surface semilandmarks (SSLM) have been digitised and warped to capture the form of the six bony epiphyses (Figure 2, Table 2): proximal femur (11 ALM, 8 curves, 105 CSLM, 28 SSLM), distal femur (12 ALM, 7 curves, 130 CSLM, 45 SSLM), proximal tibia (11 ALM, 8 curves, 71 CSLM, 23 SSLM), distal tibia (6 ALM, 5 curves, 37 CSLM, 19 SSLM), proximal metatarsal (10 ALM, 9 curves, 37 CSLM, 13 SSLM) and distal metatarsal (14 ALM, 6 curves, 64 CSLM, 34 SSLM). Digitisation and warping were carried out using ViewBox v.4.0.1.7 software (dHAL software, Kifissia, Greece).

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

	Proximal femur (A on Figure 2)
No.	ALM	No.		CSLM
1	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)
No.	ALM	No.		CSLM
1	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)
No.	ALM	No.	CSLM
1	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)
No.	ALM	No.	CSLM
1	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)
No.	ALM	No.	CSLM
1	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	No.	SSLM
10	Most anterior point of the articular surface for the medial cuneiform bone	48:53	Lateral proximal articular surface patch (n = 6)
		54:60	Medial proximal articular surface patch (n = 7)
	Distal metatarsal (F on Figure 2)
No.	ALM	No.	CSLM
1	Most medio-proximal point of the anterior outline of the medial articular eminence (MAE)	15:22	Crest between ALM 1 and ALM 13 (n = 8)
2	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
8	Most proximal point of the posterior outline of the lateral articular eminence
9	Most medio-proximal point of the posterior outline of the lateral articular eminence	No.	SSLM
10	Most latero-proximal point of the posterior outline of the medial articular eminence	79:86	Medial trochlear patch on MAE (n = 8)
11	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

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.

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 1991, 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 log-transformed 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 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).

Results

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). Male R.t. fennicus, female R.t. fennicus, male R.t. tarandus and female R.t. tarandus all showed significant differences between them for all bony elements (all P < 0.01), except between female R.t. fennicus and male R.t. tarandus for the distal parts of the femurs and tibias, as well as the proximal and distal parts of the metatarsals (P = 0.144, P = 0.526, P = 0.968 and P = 0.181, respectively; Supplementary Fig. S3). Thereby, male R.t. fennicus had the largest bony elements, female R.t. tarandus had the smallest ones, and finally male R.t. tarandus were slightly larger than female R.t. fennicus. In hybrids, males did not differ significantly from female R.t. fennicus and male R.t. tarandus, and females did not differ from female R.t. tarandus. Regarding lifestyles, free-ranging and working individuals were not statistically different (all P > 0.05). However, captive individuals were significantly smaller than free-ranging and working individuals (all P < 0.01; Supplementary Fig. S4).

Table 3. Results of the Kruskal-Wallis tests on the size for each element analysed according to the different categories. A significant contribution was considered for P-value < 0.05 (in bold)

	Df	Chi-squared	P-value
Proximal femur
subspecies (ssp.)	2	13.879	9.690e-04
sex	1	47.572	5.303e-12
lifestyle	2	20.542	3.462e-05
ssp.*sex	5	70.021	1.015e-13
ssp.*sex*lifestyle	10	71.564	2.210e-11
Distal femur
subspecies (ssp.)	2	22.771	1.136e-05
sex	1	32.485	1.201e-08
lifestyle	2	20.264	3.979e-05
ssp.*sex	5	63.658	2.127e-12
ssp.*sex*lifestyle	10	66.953	1.712e-10
Proximal tibia
subspecies (ssp.)	2	16.063	3.251e-04
sex	1	38.137	6.594e-10
lifestyle	2	24.273	5.361e-06
ssp.*sex	5	62.942	2.992e-12
ssp.*sex*lifestyle	9	65.454	1.177e-10
Distal tibia
subspecies (ssp.)	2	23.952	6.293e-06
sex	1	28.380	9.967e-08
lifestyle	2	15.373	4.591e-04
ssp.*sex	5	59.691	1.408e-11
ssp.*sex*lifestyle	10	62.733	1.098e-09
Proximal metatarsal
subspecies (ssp.)	2	22.803	1.118e-05
Sex	1	24.469	7.552e-07
Lifestyle	2	14.706	6.407e-04
ssp.*sex	5	63.051	2.841e-12
ssp.*sex*lifestyle	10	70.261	3.948e-11
Distal metatarsal
subspecies (ssp.)	2	38.821	3.717e-09
sex	1	23.758	1.093e-06
lifestyle	2	15.465	4.384e-04
ssp.*sex	5	74.684	1.083e-14
ssp.*sex*lifestyle	10	81.043	3.134e-13

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. fennicus – with 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 significantly different from all other populations (all P < 0.01), except captive individuals in both subspecies and hybrids (all P > 0.05). In R.t. fennicus, captive female individuals were significantly smaller than their wild counterparts (i.e. free-ranging), while this size difference was not significant between captive and free-ranging females in R.t. tarandus. However, the small number of captive male individuals for both subspecies did not allow us to identify any size differences with other groups, particularly between captives and free-ranging male individuals. Furthermore, there were important overlaps in the size range for many bony elements, and very few significant differences were found, especially between free-ranging female R.t. fennicus, captive male R.t. fennicus, free-ranging male R.t. tarandus, captive male R.t. tarandus, working male R.t. tarandus and captive male hybrids. However, free-ranging female R.t. fennicus appear to have smaller bony elements than male R.t. tarandus, while working reindeer appear to be male R.t. tarandus with the largest bones. Finally, captive female and male hybrids appear to have an intermediate size compared to captive females and males of both subspecies.

Table 4. Multi-test comparisons (P-values) of log-transformed centroid sizes for each element analysed between the different reindeer groups (pairwise Wilcoxon rank tests after the ‘Benjamini–Hochberg’ correction). A significant contribution was considered for P-value < 0.05 (in bold)

	R.t. fennicus ♀ free-ranging	R.t. fennicus ♀ captive	R.t. fennicus ♂ free-ranging	R.t. fennicus ♂ captive	R.t. tarandus ♀ free-ranging	R.t. tarandus ♀ captive	R.t. tarandus ♂ free-ranging	R.t. tarandus ♂ captive	R.t. tarandus ♂ working	Hybrid ♀ captive
Proximal femur
R.t. fennicus ♀ captive	0.172	-	-	-	-	-	-	-	-	-
R.t. fennicus ♂ free-ranging	1.596e-09	0.001	-	-	-	-	-	-	-	-
R.t. fennicus ♂ captive	0.181	0.500	0.500	-	-	-	-	-	-	-
R.t. tarandus ♀ free-ranging	0.001	0.306	3.697e-04	0.447	-	-	-	-	-	-
R.t. tarandus ♀ captive	7.322e-06	0.263	7.322e-06	0.306	0.710	-	-	-	-	-
R.t. tarandus ♂ free-ranging	0.005	0.040	3.963e-05	0.330	0.005	5.028e-04	-	-	-	-
R.t. tarandus ♂ captive	1.000	0.652	0.173	0.764	0.306	0.160	1.000	-	-	-
R.t. tarandus ♂ working	0.002	0.015	3.963e-05	0.306	0.004	3.697e-04	0.341	1.000	-	-
Hybrid ♀ captive	0.005	0.770	0.001	0.500	0.664	0.782	0.011	0.376	0.008	-
Hybrid ♂ captive	0.852	0.295	5.028e-04	0.393	0.041	0.003	0.260	1.000	0.173	0.087
Distal femur
R.t. fennicus ♀ captive	0.092	-	-	-	-	-	-	-	-	-
R.t. fennicus ♂ free-ranging	5.951e-07	0.001	-	-	-	-	-	-	-	-
R.t. fennicus ♂ captive	0.187	0.524	0.190	-	-	-	-	-	-	-
R.t. tarandus ♀ free-ranging	0.011	0.306	3.235e-04	0.458	-	-	-	-	-	-
R.t. tarandus ♀ captive	5.309e-06	0.011	5.309e-06	0.306	0.634	-	-	-	-	-
R.t. tarandus ♂ free-ranging	0.703	0.264	1.239e-05	0.330	0.047	0.002	-	-	-	-
R.t. tarandus ♂ captive	1.000	0.898	0.047	0.764	0.306	0.080	0.959	-	-	-
R.t. tarandus ♂ working	0.042	0.047	1.433e-04	0.306	0.008	3.232e-04	0.032	0.962	-	-
Hybrid ♀ captive	0.006	0.203	0.001	0.524	0.528	0.975	0.025	0.386	0.018	-
Hybrid ♂ captive	0.685	0.595	2.071e-04	0.403	0.069	0.016	0.937	1.000	0.203	0.081
Proximal tibia
R.t. fennicus ♀ captive	0.084	-	-	-	-	-	-	-	-	-
R.t. fennicus ♂ free-ranging	5.342e-07	5.289e-04	-	-	-	-	-	-	-	-
R.t. fennicus ♂ captive	-	-	-	-	-	-	-	-	-	-
R.t. tarandus ♀ free-ranging	0.008	0.123	2.359e-04	-	-	-	-	-	-	-
R.t. tarandus ♀ captive	2.359e-04	0.008	1.448e-05	-	0.503	-	-	-	-	-
R.t. tarandus ♂ free-ranging	0.234	0.057	5.289e-04	-	0.010	0.001	-	-	-	-
R.t. tarandus ♂ captive	0.870	0.870	0.035	-	0.201	0.129	0.577	-	-	-
R.t. tarandus ♂ working	0.008	0.009	0.001	-	0.003	5.289e-04	0.110	0.289	-	-
Hybrid ♀ captive	0.008	0.106	9.410e-04	-	0.669	0.866	0.013	0.343	0.008	-
Hybrid ♂ captive	0.619	0.508	5.861e-04	-	0.051	0.005	0.369	0.870	0.107	0.069
Distal tibia
R.t. fennicus ♀ captive	0.107	-	-	-	-	-	-	-	-	-
R.t. fennicus ♂ free-ranging	2.367e-05	0.003	-	-	-	-	-	-	-	-
R.t. fennicus ♂ captive	0.185	0.436	0.556	-	-	-	-	-	-	-
R.t. tarandus ♀ free-ranging	2.973e-04	0.025	1.858e-04	0.393	-	-	-	-	-	-
R.t. tarandus ♀ captive	1.282e-05	0.006	1.282e-05	0.370	0.587	-	-	-	-	-
R.t. tarandus ♂ free-ranging	0.544	0.621	8.761e-04	0.393	0.015	0.005	-	-	-	-
R.t. tarandus ♂ captive	0.957	0.889	0.151	0.719	0.151	0.111	0.905	-	-	-
R.t. tarandus ♂ working	0.408	0.021	0.002	0.355	0.002	7.541e-04	0.151	0.653	-	-
Hybrid ♀ captive	0.136	0.393	0.002	0.512	0.393	0.368	0.300	0.393	0.054	-
Hybrid ♂ captive	0.630	0.393	0.002	0.393	0.007	0.003	0.691	0.889	0.274	0.225
Proximal metatarsal
R.t. fennicus ♀ captive	0.007	-	-	-	-	-	-	-	-	-
R.t. fennicus ♂ free-ranging	4.570e-04	0.002	-	-	-	-	-	-	-	-
R.t. fennicus ♂ captive	0.306	0.393	1.000	-	-	-	-	-	-	-
R.t. tarandus ♀ free-ranging	1.208e-05	0.080	1.423e-05	0.269	-	-	-	-	-	-
R.t. tarandus ♀ captive	1.423e-05	0.070	3.739e-05	0.330	1.000	-	-	-	-	-
R.t. tarandus ♂ free-ranging	0.077	0.211	5.234e-05	0.193	9.576e-05	6.161e-05	-	-	-	-
R.t. tarandus ♂ captive	0.729	0.737	0.116	0.748	0.193	0.091	1.000	-	-	-
R.t. tarandus ♂ working	1.000	0.005	1.460e-04	0.215	1.208e-05	3.739e-05	0.062	0.662	-	-
Hybrid ♀ captive	0.003	0.269	0.003	0.537	1.000	1.000	0.027	0.386	0.003	-
Hybrid ♂ captive	0.572	0.220	0.013	0.393	0.003	0.003	0.580	0.737	0.696	0.081
Distal metatarsal
R.t. fennicus ♀ captive	0.036	-	-	-	-	-	-	-	-	-
R.t. fennicus ♂ free-ranging	2.613e-05	7.301e-04	-	-	-	-	-	-	-	-
R.t. fennicus ♂ captive	0.294	0.374	1.000	-	-	-	-	-	-	-
R.t. tarandus ♀ free-ranging	1.778e-06	0.050	2.980e-07	0.238	-	-	-	-	-	-
R.t. tarandus ♀ captive	5.126e-06	0.008	4.698e-06	0.339	0.752	-	-	-	-	-
R.t. tarandus ♂ free-ranging	0.007	0.486	1.012e-06	0.200	3.331e-04	1.006e-05	-	-	-	-
R.t. tarandus ♂ captive	0.478	1.000	0.061	0.733	0.203	0.091	0.802	-	-	-
R.t. tarandus ♂ working	0.435	0.027	2.613e-05	0.215	5.876e-06	2.613e-05	0.020	0.607	-	-
Hybrid ♀ captive	0.004	0.367	0.001	0.478	0.317	0.200	0.057	0.367	0.004	-
Hybrid ♂ captive	0.367	0.478	0.002	0.374	0.009	0.002	0.703	0.889	0.367	0.131

Figure 3. Boxplots of the variation in log-transformed centroid size according to the subspecies (R.t. tarandus, R.t. fennicus and hybrids), sex (male = ♂ and female = ♀) and lifestyle (free-ranging, captive and working). A: proximal femur; B: distal femur; C: proximal tibia; D: distal tibia; E: proximal metatarsal; F: distal metatarsal. The numbers in square brackets represent the number of bones analysed by elements in each group.

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). The PC1 (18.68% of the global variance) isolates free-ranging R.t. fennicus individuals fairly well on the positive part from R.t. tarandus and hybrids. In addition, the separation of individual captive and free-ranging reindeer within their respective subspecies (i.e. R.t. fennicus, R.t. tarandus and hybrids) was also relatively well visible. Thus, this resulted in a more massive morphology at the maximum of PC1, as in the free-ranging male R.t. fennicus. The distal epiphysis was broad with well-developed condyles. The medial condyle surface was comparable in size to the lateral condyle surface, both being separated by a narrow intercondylar space. The trochlea was wide with large medial and lateral lips separated by a deep trochlear groove. The femoral medial patellar margin showed a more circular section. The negative values of the PC1 exhibited a more slender morphology, as in the captive female R.t. tarandus and hybrids. The distal epiphysis was 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 PC1 (15.18% of the global variance) distinguished roughly free-ranging R.t. fennicus (males and females) on the positive part with free-ranging and captive female R.t. tarandus and female hybrids on the negative part. The captive female R.t. fennicus and working male R.t. tarandus also fell towards the negative values of the PC1. 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 free-ranging 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.t. fennicus on the negative part and R.t. tarandus on the positive part. Mountain reindeer R.t. tarandus appeared to have a distal epiphysis with a more rectangular general morphology and with larger but shallower grooves than R.t. fennicus, with a rather square distal articular surface. The lateral malleolar sulcus was also larger and more developed in R.t. fennicus (negative values of the PC2).

• Metatarsal: Despite the overlaps, the proximal part of the metatarsal appeared to show some distinctions on the PC1 (15.71% of the global variance). Indeed, R.t. fennicus tended to spread more towards positive values, while R.t. tarandus were more localised towards negative values (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.

Table 5. Results of the MANOVA tests on the shape of each element analysed according to the different categories. A significant contribution was considered for P-value < 0.05 (in bold)

	Df	Pillai	approx F	num Df	den Df	Pr
Proximal femur
subspecies (ssp.)	2	1.986	3.469	168	4	0.115
sex	1	0.997	3.547	85	1	0.403
lifestyle	2	1.982	2.586	168	4	0.182
ssp.*sex	5	4.770	1.282	405	25	0.231
ssp.*sex*lifestyle	10	9.063	1.273	760	100	0.065
Distal femur
subspecies (ssp.)	2	1.962	1.216	168	4	0.487
sex	1	0.998	5.468	85	1	0.330
lifestyle	2	1.975	1.879	168	4	0.288
ssp.*sex	5	4.853	2.046	405	25	0.017
ssp.*sex*lifestyle	10	9.339	1.859	760	100	8.222e-05
Proximal tibia
subspecies (ssp.)	2	1.988	3.899	170	4	0.095
sex	1	0.970	0.380	86	1	0.892
lifestyle	2	1.989	4.153	170	4	0.085
ssp.*sex	5	4.712	0.997	410	25	0.538
ssp.*sex*lifestyle	9	8.299	1.367	702	81	0.039
Distal tibia
subspecies (ssp.)	2	1.973	1.663	174	4	0.338
sex	1	0.995	2.163	88	1	0.502
lifestyle	2	1.977	1.954	174	4	0.273
ssp.*sex	5	4.808	1.493	420	25	0.113
ssp.*sex*lifestyle	10	8.974	1.108	790	100	0.263
Proximal metatarsal
subspecies (ssp.)	2	1.991	4.540	204	4	0.073
sex	1	0.993	1.361	103	1	0.607
lifestyle	2	1.996	9.155	204	4	0.021
ssp.*sex	5	4.810	1.278	495	25	0.234
ssp.*sex*lifestyle	10	9.321	1.461	940	100	0.009
Distal metatarsal
subspecies (ssp.)	2	1.993	4.741	226	4	0.068
sex	1	0.986	0.627	114	1	0.791
lifestyle	2	1.986	2.499	226	4	0.192
ssp.*sex	5	4.922	2.861	550	25	0.001
ssp.*sex*lifestyle	10	9.438	1.599	1,050	100	0.002

Figure 4. 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 femur; B: distal femur. The proportion of the total variance expressed by the axes PC1 and PC2, respectively, is indicated in brackets.

Figure 5. 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 tibia; B: distal tibia. The proportion of the total variance expressed by the axes PC1 and PC2, respectively, is indicated in brackets.

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.

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: 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 free-ranging 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 free-ranging 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.

Figure 7. Multivariate regression plots performed on shape data (regression scores) and log-transformed centroid size (logCS) according to the subspecies (R.t. tarandus, R.t. fennicus and hybrids), sex (male = ♂ and female = ♀) and lifestyle (free-ranging, captive and working). A: proximal femur; B: distal femur; C: proximal tibia; D: distal tibia; E: proximal metatarsal; F: distal metatarsal.

Discussion

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 (2012, 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 cover – which 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. 2015, 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. 2008, 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. 2017, 2019; Hanot and Bochaton 2018; Savriama et al. 2018). In our study, hybrid individuals are morphometrically intermediate to their parents. Specifically, female hybrids have an intermediate size compared to female R.t. fennicus and R.t. tarandus, in the same way that male hybrids have an intermediate size compared to male R.t. fennicus and R.t. tarandus, and which are also closer to female R.t. fennicus. Regarding morphology, the hybrids (both sexes included) were closer to R.t. tarandus than to R.t. fennicus. 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 free-ranging 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, 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 free-ranging 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 free-ranging 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, 2018, 2020a; Salmi 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 non-synchronous 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.

Author’s contributions

MP and AKS conceived and designed the study. AK conducted the CT scans in collaboration with MP and SN. MP collected, analysed and interpreted the morphometric data and wrote the original draft. AKS, AK and SN reviewed and edited the paper. All authors read and approved the final manuscript.

Acknowledgments

We are grateful to Tuula Pudas and Jaakko Tauriainen who prepared the working reindeer skeletons and who granted us access to the Zoological Museum collections of the University of Oulu. Thanks to Jarmo Kynkäänniemi and Mika Lappalainen who provided racing and working reindeers for our study. We would also like to thank the Department of Radiology, Oulu University Hospital (OYS) for allowing us access to the computed tomography scanner for the development of 3D models. Finally, we are very grateful to the reviewers for their constructive comments that greatly improved the quality of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 756431) and the Academy of Finland (project numbers 275635 and 308322);