Methods

Osteometric Study of Metapodial Bones and Phalanges as Indicators of the Behavioural Ecology of Modern Reindeer <i>Rangifer tarandus</i> and Implications for Reconstruction of Paleo Mobility

Authors: {'first_name': 'Ana', 'last_name': 'Gal\xc3\xa1n L\xc3\xb3pez'},{'first_name': 'Sandrine', 'last_name': 'Costamagno'},{'first_name': 'Ariane', 'last_name': 'Burke'}

Abstract

Paleolithic reindeer (Rangifer tarandus) played an important role for human populations in western and central Europe during much of the Paleolithic period. In southwestern France and in particular during the Magdalenian, reindeer frequently figures among the privileged prey of hunter-gatherer groups. However, and despite numerous attempts to reconstruct the migratory behavior of Paleolithic reindeer, there is no agreement on the degree of mobility of this prey. Modern ethological data indicate that reindeer herds adopt different mobility strategies depending on the type of habitat and the topography of the environment. Thus, our project (Emorph) aims to explore morphometric criteria (through metapodial bones and phalanges) in combination with cutting-edge methodologies like Machine Learning to identify the extent of reindeer migrations. Based initially on the study of modern caribou populations with distinct migratory behaviors, the results obtained could be applied to several Magdalenian assemblages from southwestern France in the future, with the goal of reconstructing the mobility of these tardiglacial reindeer.

Keywords: ecomorphologyosteometryreindeerMagdalenianmigrationMachine Learning 
DOI: http://doi.org/10.5334/oq.106
 Accepted on 18 Apr 2022            Submitted on 15 Jul 2021

Introduction

Palaeolithic reindeer movements are key to understanding human hunting strategies and the mobility of human groups in northern Eurasia during the Late Pleistocene. In Western Europe, for example, reindeer was such a keystone species that the Magdalenian period was initially named “l’Âge du Renne” (Lartet & Christy 1864, 1875; Weinstock, 2000). Whether the dominance of reindeer reflects selective hunting strategies or simply their local abundance has been the subject of hot debate, however (Grayson & Delpech, 2002; Mellars, 2004; Costamagno and Mateos, 2007; Costamagno et al. 2016; Fontana & Chauvière, 2018), and reconstructing the degree of mobility of Pleistocene reindeer herds is still matter of contention.

There have been numerous attempts to reconstruct the migratory behaviour of reindeer during the Upper Palaeolithic in Western Europe (Bouchud et al. 1953; Bouchud, 1954; Lacorre, 1956; Bouchud, 1966; Bahn, 1977; Spiess, 1979; White, 1980; Delpech, 1983, 1987; Gordon, 1988a, 1988b, 1990; Deplano, 1994; Weinstock, 2000a, 2000b; Fontana, 2000; Enloe, 2001; Weladji et al. 2002; Sorensen et al. 2007; Banks et al. 2008; Morin, 2008; Kuntz & Costamagno, 2011; Kuntz, 2011; Costamagno et al. 2016). Some authors (Lacorre, 1956; Bahn, 1977; Gordon, 1988a, 1988b, 1990) have proposed that reindeer migrated long distances across Europe. In southwestern France, for example, relatively long-distance migrations have been proposed along a North-South axis, between wintering grounds in the Gironde, Perigord and Quercy to summer ranges in the Pyrenees (see Map 1 in Supplementary Information). A second hypothesis, based on osteometric and dental data as well as seasonality studies, suggests relatively small-scale migrations on an East-West axis, between the Massif Central and the Gironde (see Map 1 in Supplementary Information) (Delpech, 1983, 1987; Kuntz & Costamagno, 2011). Thirdly, based on dental data and reindeer antlers, other authors have proposed that reindeer herds were present all year round in the Perigord (Bouchud et al. 1953; Bouchud, 1954; Deplano, 1994; Fontana, 2000, 2012, 2017; Fontana & Chauvière, 2018).

Seasonal data accumulated since the 1990s seem to indicate the absence of large-scale migrations on a North-South axis (Kuntz, 2011). Furthermore, we cannot discard the possibility that reindeer underwent evolution during the Late Glacial under selective pressure of significant climatic changes, which could have affected their migratory behaviour (Delpech, 1983, 1987; Costamagno et al. 2016).

In the absence of consensus about the expected migratory behaviour of Late Pleistocene reindeer in southwestern Europe, we propose the implementation of an actualistic approach from an ecomorphological perspective. Extant reindeer migrate seasonally to avoid excessive predation on the calving grounds, pests and the depletion of local food resources (Bergerud, 2000). However, reindeer herds adopt different mobility strategies according to climate conditions, habitat type and topography. For example, mid- to long-distance migrations take place in open habitats such as tundra and steppe (Miller, 1990; Schaefer et al, 2000; Goddard, 2009; Seip, 2010) whereas animals living in closed habitats, in or near woodlands tend to migrate shorter distances (Thomas & Grey, 2002; COSEWIC, 2011). Since Pleistocene environments do not necessarily have modern analogues, = ethological data can be used to propose hypothetical migration patterns, but another source of information is required to test the hypotheses. A link between migration distance (degree of mobility) and the functional morphology of limb bones of a variety of taxa has been proposed elsewhere in the context of ecomorphological research, which relies on osteometric analyses.

The use of linear measurements in osteometric analyses has been successfully applied to taxonomic classification (Prummel & Frisch, 1986; Puputti & Niskanen, 2008; Castaños et al. 2012; Gruvier et al. 2015;) and to determining the age and sex of individuals to establish herd structure (Davis, 1987, 1985, 2000; Berteaux & Guintard, 1995; Weinstock, 2000; Greenfield, 2002; Munro et al. 2011; Arceredillo et al. 2011; Zeder, 2001; Davis et al. 2012; Galán López & Domínguez-Rodrigo, 2014). Osteometric analyses can be used to reconstruct the herd structure of a local assemblage and, in gregarious, migratory species, this information can be used to reconstruct seasonal herd mobility on a regional scale which in turn helps tracking human mobility (Delpech 1983; Weinstock; 2005; Costamagno & Kuntz, 2011). This information cannot tell us how far individual herds migrated, however.

Ecomorphological studies use osteometric analysis to reconstruct the habitat type within which a taxon lived (Spencer 1997; Kappelman, 1988; Bishop, 1994; Plummer & Bishop, 1994; Kappelman et al. 1997; DeGusta and Vrba, 2002, 2003, 2005; Kovarovic, 2004; Kovarovic & Andrews, 2007; Plummer et al. 2008; Van Asperen, 2010). African bovids, for example, have been studied by DeGusta & Vrba, (2002, 2005a, 2005b), who use the astragalus and proximal, intermediate and distal phalanges, and Kappelman (1987, 1991), who uses the femur to reconstruct bovid habitats. Plummer and Bishop (1994) use metapodials to group bovids according to three broad habitat types (open, intermediate, and closed). Plummer et. al (2008) demonstrated the use of astragali for ecomorphological analysis for 36 extant African antelope species. Van Asperen (2010) studied metapodia and first phalanges for an ecomorphomogical study of caballoid horses. Kovarovic & Andrews (2007) collected and successfully tested long bones, carpals, tarsals and phalanges from mainly extant bovids and also cervids and tragulids from known habitat types.

As far as reindeer are concerned, several studies on the Fennoscandian subspecies (Rangifer tarandus tarandus and Rangifer tarandus fennicus) have been carried out in relation to its morphology (e.g. cross-sectional morphology, entheseal changes) and identification in the archaeological record (particularly to distinguish the different stages of its domestication). They have provided an important and valuable understanding of this subject, which is key in reconstructing the mobility and lifestyles of northern peoples (Niinimäki, & Salmi, 2016; Hull, 2019; Hull et al. 2020; Salmi et al. 2020; Pelletier et. al. 2020; Nomokonova et. Al. 2021; Hull et.al, 2021; Pelletier et.al, 2021a; Pelletier et. al. 2021b).

The aim of the present study is to test if it is possible to determine the degree of mobility of reindeer using linear measurements from metapodial bones and phalanges and if so, create a referential framework that can be applied to the study of the archaeological record (archaeological sites from South-Western France) to test hypotheses about reindeer migration patterns during the Magdalenian (18,000–14,000 cal BP). To this end we will study the North American reindeer, or caribou. The present study is part of a larger project (Emorph) in which we address several methodologies (linear measurements-machine learning and Geometric Morphometrics [GMM]) to investigate reindeer ecomorphology. Furthermore, although we know that GMM has proved a very efficient method to study morphology, we are aware that access to 3D scanners is not always possible, and in many cases, archaeological faunal collections cannot be taken from museums, so the objective was also to explore a methodology accessible to all.

Extant North American Reindeer

Rangifer tarandus (known as reindeer in Eurasia and called caribou in North America) is the most abundant large land-based mammal in northern North America (Hummel & Ray, 2010), where it occupies a wide range of biomes from boreal or coniferous forest to arctic tundra and polar deserts (Bergerud, 2000). Reindeer have a limb structure well suited to migration in complex environments (Zhang et al, 2017, 2020) and they can adapt to several terrains, e.g. tundra, ice, snow, mature forest (Wareing et al. 2011). Their concave hooves are adapted to the harsh and often treacherous northern environment. Their hoof is composed of two large crescent-shaped toes and two dewclaws allowing reindeer to walk on all their toes when deep snow covers the ground. In the fall, their feet become harder, and their edges become sharper so that they can easily break through ice layers to search for food, and act as paddles for swimming, or as ice picks when walking along steep, rocky, and icy mountain slopes (Hummel & Ray, 2010).

Four extant subspecies of North American reindeer, or caribou, are recognized: Rangifer tarandus granti (Grant’s caribou, mainly found in Alaska), Rangifer tarandus groenlandicus (Barren-Ground caribou, Nunavut and the NWT), Rangifer tarandus pearyi (Peary Caribou, Nunavut, NWT), and Rangifer tarandus caribou (which includes Eastern caribou/Migratory woodland caribou, mountain caribou and boreal caribou/forest-dwelling caribou). North American caribou can also be classified into four main “ecotypes”. An ecotype is a population or group of populations adapted to a particular set of environmental conditions (COSEWIC, 2011). These are: Migratory Tundra, Boreal Forest, Mountain, and Peary caribou ecotypes (Festa-Bianchet et al. 2011). The picture is somewhat complicated by the fact that herds may inhabit different ecotypes within a single subspecies.

The setting within which a herd lives is correlated with the distance of their seasonal migrations. Rangifer tarandus caribou, which lives in different environments, performs medium, long, and short distance movements. In Québec, caribou living south of 55ºN are forest-dwellers and move relatively short distances seasonally (Miller, 2003; Nuttal, 2004), often between 80–100 kilometres (Thomas and Grey, 2002; Seip & McLellan, 2010; Schaefer, 2010;) and 50–150 kilometres (depending on the herd) (Mallory and Hillis, 1998). The George River Caribou Herd (GRCH) and Leaf River Caribou Herd (LRCH) on the other hand, are mixed forest-tundra dwelling populations of R.t.caribou, covering migration distances ranging between 1120 up to 1770 kilometres on average per year (Leblond et al. 2016) moving from open-tundra to boreal forest habitats (COSEWIC, 2017). In British Columbia, R. t. caribou herds are distributed in different ecotypes, including: Mountain, Northern, and Boreal habitats, while in Alberta, they live in Mountain and Boreal habitats. Some R. t. caribou herds perform altitudinal migrations (vertically) while others perform longitudinal (horizontally; planarly) migrations (Cavedon et. al. 2019; Edmonds, 1988; Hegel & Russel, 2013; Stuart-Smith et. al. 1997; Weaver, 2016). In some cases, they can reach overall migration distances between 140–240 kilometres during the spring migration, depending upon the location of their calving area (Weaver, 2006), while others travel shorter distances (Seip, 2010); e.g., under 100 kilometres (Cumming & Beange, 1987; Saher, 2005) shifting their seasonal ranges vertically in response to snowfall conditions, forage availability and to avoid predation (Heard & Vagt, 1998; Taillön, 2013). In addition, some herds move twice a year, while others can make four migrations each year.

Rangifer tarandus pearyi inhabit the Canadian Arctic Archipelago (except Baffin Island). During the Wisconsin Glaciations, they probably survived in at least one glacial refugium in the High Arctic, as not all the islands were glaciated (Seip & McLellan, 2010). Their habitat is treeless arctic tundra (Geist, 1998; Bergerud, 2000; Gunn, 2010; Taillön, 2013). Peary caribou herds are made up of tens to hundreds of individuals, very different from other tundra populations. They perform seasonal migrations between islands, crossing the sea ice and swimming (Miller & Gunn, 1978; Miller, 1990, 2002; COSEWIC, 2011; Jenkins et al. 2011) reaching migration distances of up to 500 kilometres (Miller, 1990; Geist, 1998). They may move even further from islands to the mainland in severe winters (Miller, 1990; COSEWIC, 2002) and make erratic, large-scale movements among islands (COSEWIC, 2015). In their annual cycle, they either spend most of their time on one island or travel between two or more islands. For example, it has been estimated that seasonal migrations within the Prince of Wales-Somerset-Boothia complex (Canadian Arctic Archipelago) can cover between 300–500 kilometres or more (COSEWIC, 2004).

Rangifer tarandus granti are mainly found in Alaska and adjacent territories in Canada. Very similar to Rangifer tarandus groenlandicus (Barren-Ground), they are often considered a single subspecies –R.t. groendlandicus. The Barren-Ground caribou is the caribou “par excellence in popular imagination” (Gunn, 2010). These tundra caribou include the largest caribou herds in North America, numbering hundreds of thousands of individuals, undertaking long distance migrations travelling back and forth between northern tundra in the spring, where they remain until the fall when they migrate to winter ranges in the boreal forest (Seip & McLellan, 2010; Taillön, 2013; Bergerud, 2000, 2013). They undertake long distances during their annual migration, ranging between 800 and 5055 kilometres depending on the herd and the year (Fancy et al. 1989; Ferguson & Messier, 2000; Saher, 2005; COSEWIC, 2016; Joly, 2019). They live in both prostrate dwarf shrub tundra and upright shrub tundra. Further, tundra is characterized by the presence of permafrost, which means that most of the area is wet in summer, resulting in an abundant lichen cover (Gunn, 2010); as a result, lichens are the main source of fodder (Taillön, 2013).

Biologists in Canada generally distinguish two mobility patterns for caribou: migratory and sedentary (Festa & Bianchet et al. 2011). Migratory reindeer travel more than 200 kilometres while those considered to be sedentary travel distances of less than 200 kilometres (Wittmer et al. 2006). Therefore, it is important to clarify that the term “sedentary” referring to reindeer, as used in the present work, which follows Wittmer (Wittmer et al. 2006), is removed from what is meant in an archaeological context. Here, sedentary caribou refers to herds that travel less than 200 kilometres (Wittmer et al. 2006). As a mobile animal, current caribou can live in large groups, migrate extensively, and migrate between seasonal calves and breeding grounds. However, some caribou groups are smaller, and the migration between calves and wintering grounds is shorter. (Hummel & Ray, 2010).

Material and Methods

Sample

For this research we compiled a sample of 54 metacarpals, 58 metatarsals, 135 proximal phalanges, 88 intermediate phalanges and 78 distal phalanges belonging to the four subspecies of caribou described above (Table 1). Specimens with bone pathologies were excluded from the study. Because of the specialization and elongation of metapodial bones in artiodactyls, they have been used in a number of morphometric studies assessing body size, locomotor behaviour, and habitat preference (Scott, 1985; Purdue, 1987; Plummer & Bishop, 1994; Scott, 2004; Reese, 2015).

Table 1

Number of metacarpals (Mtc), metatarsals (Mtt), first phalanges (Ph1), second phalanges (Ph2) and third phalanges (Ph3) for every population included in the present study.


CARIBOU’S NAME SUBSPECIES HABITAT MOBILITY PATTERN MTC MTT PH1 PH2 PH3

Peary Caribou Rangifer tarandus pearyi Arctic Tundra (Open) Migratory 13 16 14 0 0

Barren-Ground Rangifer tarandus groenlandicus Tundra
(Open)
Migratory 11 6 10 4 13

Grant Caribou Rangifer tarandus granti Tundra
(Open)
Migratory 3 2 15 8 13

Woodland Caribou Rangifer tarandus caribou Boreal Forest
(Closed)
Sedentary 6 7 21 14 17

Eastern or Migratory Woodland Caribou Rangifer tarandus caribou Tundra
(Open)
Migratory 10 11 27 22 17

Mountain Caribou Rangifer tarandus caribou Mountain
(Closed)
Migratory 5 8 28 15 11

Mountain Caribou Rangifer tarandus caribou Mountain
(Closed)
Sedentary 6 8 17 17 7

Samples were obtained from biologists in the Ministry of Forests, Lands, Natural Resource Operations and Rural Development, British Columbia, who provided us with mountain caribou samples and related information about herds (type of herd, sex if known, and location). Peary caribou, Barren-Ground and Grant caribou specimens were collected from the Canadian Museum of Nature (Ottawa), which owns one of the largest collections of caribou in the country. Eastern migratory and short distance woodland caribou were obtained with the help of biologists from the Forest and Wildlife branch of the Québec Government, the Prehistory and Bioarchaeology lab from Laval University and the archaeozoology laboratories of the University of Montréal.

Due to the fact that sedentary populations are classified as “Endangered”, the acquisition of samples from these herds was more difficult. Since they are barely extant in museums, it was only possible to acquire them through biologists when a caribou died in a traffic accident or was found dead in the forest, resulting in uneven sample sizes.

Linear measurements of podial elements were taken using a digital calliper following Von den Driesch (1976) and Klein et al. (2010) for metacarpals and metatarsals and DeGusta & Vrba (2005) for phalanges; additional measures (HPAS, WPAS, WM) were also taken.

In phalanges, to avoidthe effects of entheseal changes (Hull et. al. 2020), these were taken at the extreme of the bones, away from entheses.

Description as follows

Metacarpals and metatarsals (Figures 1 and 2):

Measurements taken on metacarpals and described in the present work. See description Methods section of text
Figure 1 

Measurements taken on metacarpals and described in the present work. See description Methods section of text.

GL: Greatest bone length (Von den Driesch, 1976).

Bp: Greatest breadth of the proximal end (Von den Driesch, 1976).

Dp: Greatest depth of the proximal end (Von den Driesch, 1976).

SD: Smallest breadth of the diaphysis (Von den Driesch, 1976).

Bd: Greatest breadth of the distal end (Von den Driesch, 1976).

Dd: Greatest depth of the distal end (Von den Driesch, 1976).

TMLMIN: Minimum mediolateral diameter of the medial trochlea (Klein et al. 2010).

Measurements taken on metatarsals and described in the present work. See description Methods section of text
Figure 2 

Measurements taken on metatarsals and described in the present work. See description Methods section of text.

GL: Greatest bone length (Von den Driesch, 1976).

Bp: Greatest breadth of the proximal end (Von den Driesch, 1976).

Dp: Greatest depth of the proximal end (Von den Driesch, 1976).

SD: Smallest breadth of the diaphysis (Von den Driesch, 1976).

Bd: Greatest breadth of the distal end (Von den Driesch, 1976).

Dd: Greatest depth of the distal end (Von den Driesch, 1976).

TMLMIN: Minimum mediolateral diameter of the medial trochlea (Klein et al. 2010).

Measurements taken on first phalanges and described in the present work. See description Methods section of text
Figure 3 

Measurements taken on first phalanges and described in the present work. See description Methods section of text.

GL: Greatest bone length (Von den Driesch, 1976).

Bp: Greatest breadth of the proximal end (Von den Driesch, 1976).

WI: Intermediate Width (DeGusta & Vrba, 2005).

Bd: Greatest breadth of the distal end (Von den Driesch, 1976).

LM: Midline Length. The minimum proximal distal dimension along the dorsal line (DeGusta & Vrba, 2005).

HI: Intermediate Height took on the ventral dimension of the shaft at midshaft (DeGusta & Vrba, 2005).

HD: Distal Height-the midline dorsal ventral dimension just proximal to the distal articular surface (DeGusta & Vrba, 2005).

HP: Proximal Height took on the ventral dimension of the proximal articular end, perpendicular to its major proximal distal axis (DeGusta & Vrba, 2005).

Dp: Greatest depth of the proximal end (Von den Driesch, 1976).

HPAS: Height Proximal Articular Surface. Height of the lateral articular facet of the proximal end (the authors).

WPAS: Width Proximal Articular Surface. Width of the lateral articular facet of the proximal end (the authors).

Measurements taken on second phalanges and described in the present work. See description Methods section of text
Figure 4 

Measurements taken on second phalanges and described in the present work. See description Methods section of text.

GL: Greatest bone length (Von den Driesch, 1976).

LM: Midline Length. The minimum proximal distal dimension along the dorsal line.

Bd: Greatest breadth of the distal end (Von den Driesch, 1976).

LS: Superior Length- The proximal distal dimension of the dorsal surface, measured from the most proximal midline point of the dorsal surface of the proximal end to the most distal midline point on the distal articular surface (DeGusta & Vrba, 2005).

LI: Inferior Length- The proximal distal dimension of the ventral surface, measured from the most proximal midline point of the ventral surface of the proximal end to the most distal midline point on the distal articular surface (DeGusta & Vrba, 2005).

HD: Distal Height -The dorsal ventral dimension of the distal end, measured just proximal to the distal articular surface (DeGusta & Vrba, 2005).

HL: Lateral Height- The dorsal ventral dimension of the lateral portion of the proximal articular facet, measured from the most ventral point to the most dorsal point (DeGusta & Vrba, 2005).

HM: Medial Height -The dorsal ventral dimension of the medial portion of the proximal articular facet, measured from the most ventral point to the most dorsal point (DeGusta & Vrba, 2005).

Dp: Greatest depth of the proximal end (Von den Driesch, 1976).

Measurements taken on third phalanges and described in the present work. See description Methods section of text
Figure 5 

Measurements taken on third phalanges and described in the present work. See description Methods section of text.

Ld: Length of the dorsal surface (Von den Driesch, 1976).

LI: Inferior Length (DeGusta & Vrba, 2005).

HT: Total Height-The maximum ventral dorsal dimension of the proximal end (DeGusta & Vrba, 2005).

WB: Basal Width-The maximum medio-lateral dimension taken at the ventral base of the proximal articular facet (DeGusta & Vrba, 2005).

WM: Medium Width, taken in the middle of the proximal articular surface (the authors).

Statistical and Machine Learning analyses

RStudio 1.3.1093(www.r-project.org) software was used to perform the statistical analyses. First, an exploratory analysis (EDA) was carried out (using “tidyverse”, “ggplot2”, “GGally” packages in R). All the variables were scaled before any analysis to ensure that they were not dependent on body size. In the exploratory analysis, multivariate normality was tested using the Mardia test (“QuantPsyc” R package), and the impact of three predictor variables, habitat (mountain, boreal forest, tundra), mobility (migratory and sedentary) and sub-species (Rangifer tarandus granti, groenlandicus, pearyi, and caribou) was tested. After exploratory analysis and based on the natural distribution of the data, the habitat category was simplified to include two classes: “open” for tundra and “closed” for mountain and boreal forest in the predictive analyses.

Then, a MANOVA test (“stats” R library) for metapodials and a permutation MANOVA (“vegan” R library) in the case of phalanges, were carried out considering all of the measurements taken for each type of bone in order to determine the weight of each predictor variable (mobility, habitat, subspecies) on the data distribution. A second MANOVA was carried out on the distal measurements of metapodial bones to test the differences between the variables, as this is the most commonly encountered part in archaeological sites. To complete the exploratory analysis, Principal Component Analyses (PCA) (“factoextra” R library) were performed to observe data distribution. Individuals were grouped according to subspecies, mobility pattern, and habitat (Table 2). Having verified that mobility was a good predictor, predictive analyses were carried out with this variable.

Table 2

Codes used in this study.


CODE MEANING

Rtcaribou_MS Rangifer tarandus caribou Mountain Sedentary

Rtcaribou_BS Rangifer tarandus caribou BorealForest Sedentary

Rtpeary_TM Rangifer tarandus peary Tundra Migratory

Rtgroenlandicus_TM Rangifer tarandus groenlandicus Tundra Migratory

Rtgranti_TM Rangifer tarandus granti Tundra Migratory

Rtcaribou_MM Rangifer tarandus caribou Mountain Migratory

Rtcaribou_TM Rangifer tarandus caribou Tundra Migratory

Regarding first, second, and third phalanges, differences between anterior and posterior elements were tested using a one-way MANOVA and the phalanges were classified according to the methodology that it is described as follows.

Lastly, we applied Machine Learning (ML) techniques, testing five models to identify the best possible classification algorithms for our purposes. Machine learning techniques have started to be applied in the field of archaeozoology and it has been demonstrated they can represent powerful tools for classification purposes (Lefebvre et al. 2016; Cifuentes-Alcobendas et al, 2019; González-Molina et. al 2020; Domínguez-Rodrigo et. al. 2020).

A total of 70% of the sample was used for the training models, except for metacarpals, metatarsals and anterior proximal phalanges, for which60% was used due to lower sample sizes. The remaining 30% (40% for metacarpals, metatarsals, and anterior proximal phalanges) of the sample was retained as the test sample. This avoided dealing with the bias/variance trade-off (González-Molina et. al 2020; Domínguez-Rodrigo et. al. 2020). Variables were scaled prior to any analysis to avoid any bias. Although some of the variables in phalanges showed non-normal distribution, the machine learning algorithms used in the present study do not require any data processing to deal with normality, skewness, or collinearity (Domínguez-Rodrigo & Baquedano, 2018).

One of the advantages of ML techniques is that they enable the user to estimate the performance of the model. Here, we selected 10-fold cross validation, which is standard for estimating model performance. 10-fold CV randomly divides the data into 10 completely separate random and similar sized partitions. After the process of training and evaluating the model has occurred 10 times with different combinations, the average performance across partitions is reported (Domínguez-Rodrigo & Baquedano, 2018).

When all the model runs are completed, the best ones are selected combining the least amount of error and the greatest accuracy. Cost values are tested vis-à-vis accuracy with the caret function “tuneLength” set to 10. The tuning parameter selected for measuring model performance was the “Kappa” parameter. The “Kappa statistic” measures accuracy by calculating the probability of a correct prediction occurring by chance alone. Kappa values range from a maximum value of 1 (perfect agreement between the model’s predictions and the true values) to –1 (rare occurrence or imperfect agreement). Thus, Kappa values of 0.3–0.6 show reasonable agreement. Values higher than these indicate a high agreement between the expected accuracy and the documented one (Domínguez-Rodrigo & Baquedano, 2018; González-Molina et. al 2020).

Beyond the classification accuracy (the number of correct assessments/number of total assessments) and the Kappa statistic, there are other measures for evaluating the performance of the algorithms such as: “Sensitivity”, the proportion of correctly predicted positives to all true positive events. High sensitivity would suggest a low type II error rate or high statistical power and “Specificity”, or the proportion of correctly predicted negatives to all true negative events. High specificity would suggest a low type I error rate (Domínguez-Rodrigo & Baquedano, 2018).

Several Machine Learning models were tested in order to find the best algorithms (the one with the best accuracy classification values). To account for typical bone fragmentation patterns, ML analyses were carried out on all the variables and then proximal and distal variables separately. Mobility pattern (migratory and sedentary) and habitat (open and closed) – when relevant – were tested. The five algorithms that produced the best classification results are described below:

Support Vector Machines (SVM). SVMs are powerful and highly flexible supervised learning models for classification that classifies unlabelled data, achieving a high accuracy in this task. This method seeks to create the input space to a higher dimension via a kernel function (radial and polynomial in this study) and in that transformed feature space, find a hyperplane that will result with the maximal margin of separation between the two classes (González-Molina et. al 2020). Here, the “e1071” and “caret” R libraries were used.

Neural Network (NN). A neural network is a series of algorithms that work similarly to the neural networks of the human brain. It creates nodes, which hierarchically build a network of synthesized information through nonlinear regression methods. Weighted combinations of the inputs are created and put through some function (in this case, the sigmoid function) to produce the next layer of inputs. This next layer goes through the same process to produce either another layer or to predict the output, which is the final layer (Lantz, 2013). Here, “neuralnet” and “caret” R libraries were used to carry out it.

Random Forest (RF). RF is an ensemble-based method used for classification and regression (Breiman 2001). It combines the base principles of bagging with random feature selection to add diversity to the decision tree models. After an ensemble of trees is generated, the model uses a vote to combine the trees’ predictions. This method avoids overfitting through the use of OOB (“out of the bag”) technique, estimating how many iterations are needed to minimize the OOB error (Lantz, 2013; Kuhn & Johnson, 2013; Domínguez-Rodrigo & Baquedano, 2018; González-Molina et. al 2020). Here, “randomForest” and “caret” libraries were used and the final value used for the selected model was mtry = 5.

Extreme Gradient Boosting Tree (XgbTree). This algorithm is an implementation of gradient boosting machines (Friedman 2001) and is an ensemble method where errors are minimized by a gradient descent algorithm to produce a final model. Here it is used with a tree as the base learner. For that, it uses decision trees composed of the series of binary questions and the final predictions happen at the leaf. The trees are constructed iteratively until a stopping criterion is met. In this study, the “gbm” and “xgboost” libraries were used.

Stacking. Stacking, or stacked generalisation, is an ensemble method that combines the predictions of several base models in order to construct a new, optimal predictive model and improve the overall performance (Wolpert, 1992). The idea of ensemble learning is to build a prediction model by combining the strengths of a collection of simpler base models (Trevor et al. 2008, Clark, 2013) and increase the predictive power of a classifier (Nti et al. 2020). Different combinations of five models (Linear Discriminant Analysis (lda), Extreme Gradient Boosted Tree (xgbTree), Random Forests (rf), Support Vector Machines (svmRadial) and Gradient Boosting (glmboost)) were used here (not necessarily all at the same time) to find the optimal one. One of the advantages of stacking is the avoidance of overfitting when sample size is relatively small and of course, it improves the overall predictive performance of each of the models included in the stack (Nti et al. 2020). Here, “caret” and “caretEnsemble” R libraries were used.

Blending. This ensemble approach is very similar to stacking technique. The only difference is that, while stacking uses a test dataset for prediction, blending uses a holdout dataset from the training set to make predictions. That is predictions take place on only the validation dataset from the training set (Nti et al. 2020). The outcome of the predicted dataset and validation dataset is used for building the final model for predictions on the test dataset. Same models as explained above were also implemented in this case. Here, “caret” and “caretEnsemble” R libraries were also used.

Results

Exploratory Data Analysis (EDA) is included in the supplementary data to this article. The Mardia test showed that the variables from metacarpals and metatarsals were normally distributed (p > 0.05), while variables for first, second, and third phalanges had non-normal distribution (p < 0.05).

Two-way MANOVA (Table 3) for metacarpals using all the variables showed significant values for mobility pattern and subspecies categories. When only distal variables are considered, both mobility pattern (p = 0.0002) and habitat (p = 0.01) are significant, although they are more important in the case of mobility predictor. For metatarsals (Table 3), the three predictors are significant, including habitat (p < 0.001) and pattern (p < 0.001) but only on distal measurements. The relative weight of habitat is very low, however, while mobility (Pattern) carries the most weight (Table 3). So, in this case, habitat predictor was also included in later machine learning analysis.

Table 3

Two-way Manova on metatarsal and metacarpals. Significant level <0.05.


MANOVA

DF TEST STATISTIC (WILKS’LAMBDA) F P

Metacarpus

Subspecies 3 0.07825 7.6557 1.310e–13

    Pattern 1 0.37433 9.3123 9.941e–07

    Habitat 2 0.72408 0.9760 0.4852

Metatarsus

Subspecies 3 0.0589 9.9272 <2.2e–16

    Pattern 1 0.48365 6.5583 2.749e–05

    Habitat 2 0.51972 2.3780 0.007576

Regarding phalanges, due to the fact that they were non-normally distributed, a one-way permutation Manova (PERMANOVA) was carried out. First of all, separation between anterior and posterior phalanges was tested (Table 4). That was significant for first (p = 0.003) and second (p = 0.001) phalanges, and non-significant for third phalanges (p = 0.834). Thus, Support Vector Machines (SVM) algorithm successfully identified anterior and posterior (Table 5) first phalanges using all measurements in 97,67% of cases (using Radial Kernel and after tuning the model. See supplementary material). When proximal variables were included, Neural Network classified correctly 93% of them. However, distal variables poorly performed with only 65% correctly classified using SVM algorithm. Anterior and posterior second phalanges were successfully separated in 94% of the cases using all the variables in the Blending test. Proximal or distal measurements tested separately did not produce satisfactory results.

Table 4

One way PERMANOVA to test statistically significant differences for anterior and posterior phalanges. Significant level <0.05.


PERMANOVA

DF MEAN SQUARE F R2 P

Anterior posterior 1st phalanx 1 4.497 10.81 0.070 0.003

Anterior posterior 2nd phalanx 1 6.638 18.07 0.135 0.001

Anterior posterior 3rd phalanx 1 0.090 0.28 0.006 0.834

Table 5

Performance of machine learning algorithms depending on the skeletal part to separate anterior and posterior first and second phalanges.


BONE MEASUREMENTS ALGORITHM ACCURACY KAPPA 95%CI SENSITIVITY SPECIFICITY BALANCED ACCURACY

First Phalanx All Support Vector Machines 97.7% 0.95 0.87–0.99 1 0.95 0.97

Proximal Support Vector Machines 90.7% 0.81 0.77–0.97 0.89 0.91 0.90

Proximal Neural Network 93% 0.85 0.80–0.98 0.89 0.95 0.92

Distal Support Vector Machines 65% 0.24 0.49–0.78 0.31 0.91 0.61

Second Phalanx All Blending 94% 0.87 0.80–0.99 0.85 1 0.92

Proximal Extreme Gradient Boosting Tree 65% 0.26 0.47–0.80 0.50 0.76 0.63

A two-way permutation Manova was then carried out on anterior and posterior first and second phalanges respectively (Table 6) and third phalanges altogether reporting significant values for subspecies predictor in the six cases, mobility in almost all cases (except posterior first phalanx and anterior second phalanx, and almost significant for posterior second phalanx). Habitat was not significant in second (anterior and posterior) phalanx.

Table 6

Two-way PERMANOVA on phalanges. Significant level <0.05.


PERMANOVA

DF MEAN SQUARE F R2 P

Anterior First phalanx

Subspecies 3 6.07 23.43 0.51 0.001

    Pattern 1 2.26 9.48 0.06 0.002

    Habitat 2 0.97 4.06 0.05 0.009

Posterior First phalanx

Subspecies 3 7.71 31.53 0.51 0.001

    Pattern 1 0.56 2.32 0.01 0.118

    Habitat 2 1.88 7.71 0.08 0.002

Anterior Second phalanx

Subspecies 2 4.71 15.31 0.42 0.001

    Pattern 1 0.71 2.31 0.03 0.121

    Habitat 2 0.21 0.70 0.01 0.529

Posterior Second phalanx

Subspecies 1 3.11 8.33 0.13 0.003

    Pattern 1 1.10 2.94 0.04 0.067

    Habitat 1 0.35 0.94 0.01 0.377

Third phalanx

Subspecies 2 2.26 14.02 0.26 0.001

    Pattern 1 0.56 3.48 0.03 0.024

    Habitat 2 1.13 7.05 0.13 0.001

PCAs were carried out on the metacarpals and metatarsals using all the variables in order to explore data and morphological characteristics. In metacarpals PCA, PC1 and PC2 accounted 92.4% of the variance. The variables that most affect the axes are: Bp (greatest breadth of the proximal end), SD (smallest breadth of the shaft), and Bd (greatest breadth of the distal end). They show wider distal ends in the case of the sedentary group. If we observe the mobility pattern, the confidence ellipses (95% confidence around the mean point of each category) individuals considered as sedentary, cluster in the same area apart from the rest of the individuals (Figure 6).

Principal Component Analysis on metacarpals. The seven groups take into account subspecies, habitat and mobility
Figure 6 

Principal Component Analysis on metacarpals. The seven groups take into account subspecies, habitat and mobility.

The displayed sample documents that the separation provided by the variables used is not dependent on sex or body size. Although adult males are bigger that females (Taillon, 2013), here it can be seen how they cluster together regardless of sex.

In PCA on metatarsals (Figure 7), PC1 and PC2 retained 91.9% of the total variance. SD (smallest breadth of the diaphysis) is almost uncorrelated with the other variables. Distal variables are affecting more the sedentary group, while proximal measurements do it on migratory groups. As observed in metacarpals, sedentary individuals show wider ends than migratory ones and Dd (G = greatest depth of the distal end.) and Bd (greatest breadth of the distal end) have a high contribution on the axes (see Supplementary information). However, in this case the dichotomy between open andclosed habitats is even clearer than in the metacarpal. (Figure 8).

Principal Component analysis on metatarsals. The seven groups take into account subspecies, habitat and mobility
Figure 7 

Principal Component analysis on metatarsals. The seven groups take into account subspecies, habitat and mobility.

Principal Component Analysis on metatarsals. According to habitat (dichotomy open-close)
Figure 8 

Principal Component Analysis on metatarsals. According to habitat (dichotomy open-close).

Regarding the classification analyses, almost all the machine learning algorithms showed an accuracy over 80% set for both mobility and habitat. In relation to mobility pattern (Table 7), metacarpals show accuracies higher than 90% to classify migratory and sedentary patterns, XgbTree performed the highest score (95,2%) on proximal variables, while stacking show a 90,4% with all the variables and 90% for distal measurements using both neural network and blending algorithms.

Table 7

Details on the performance of Machine Learning algorithms according to mobiliy pattern (migratory/sedentary) as predictor variable.


BONE MEASUREMENTS ALGORITHM ACCURACY KAPPA 95%CI SENSITIVITY SPECIFICITY BALANCED ACCURACY

Metacarpus Complete Stacking 90,4% 0.76 0.69–0.98 0.87 1 0.93

Metacarpus Proximal XgbTree 95,2% 0.85 0.76–0.99 1 0.80 0.90

Metacarpus Distal Neural Network and Blending 90% 0.73 0.68–0.98 0.87 1 0.93

Metatarsus All Random Forest 86,9% 0.59 0.66–0.97 1 0.50 0.75

Metatarsus Proximal Stacking 82% 0.42 0.61–0.95 1 0.33 0.66

Metatarsus Distal Support Vector Machines 91,3% 0.77 0.71–0.98 0.94 0.83 0.88

Anterior First Phalanx All Stacking 90% 0.80 0.70–0.98 0.92 0.87 0.9

Posterior First Phalanx All Stacking 91,3% 0.77 0.71–0.98 0.94 0.83 0.88

Anterior Second Phalanx All Neural Network 88,2% 0.76 0.63–0.98 1 0.81 0.90

Proximal Blending 88% 0.74 0.63–0.98 0.9 0.83 0.87

Distal Stacking 76% 0.48 0.50–0.93 0.81 0.66 0.74

Posterior Second Phalanx All Neural Network 93,3% 0.86 0.68–0.99 1 0.88 0.94

Third Phalanx All Extreme Gradient Boosting Tree 92% 0.78 0.73–0.99 1 0.71 0.85

Proximal Stacking 80% 0.52 0.59–0.93 0.83 0.71 0.77

A somewhat lower performance is observed in metatarsals (Table 7) with all (86.9%) and proximal measurements (82%) applying RF and Stacking methods. However, distal variables reach a 91.3% through SVM.

First phalanges (both anterior and posterior) had a successful performance only when all the variables were considered in the analysis (Table 7). Thus, the anterior first phalanx classified correctly a 90% of the testing sample with Stacking methods. Posterior first phalanges had a similar accuracy value (91.3%) for Stacking analyses.

Anterior and posterior second phalanges (Table 7) achieved high accuracy values (88% and 93.3% respectively) for NN, although accuracy decreased in proximal (Blending test-88%) and distal (Stacking-76%) measurements.

Finally, third phalanges (Tabl. 7) reached successfully a 92% of accuracy in the test set through XgbTree algorithm. This accuracy decreased when only measurements from joint surfaces were used (Stacking-80%).

Considering habitat (open/closed) when it was relevant (according to significant values from MANOVA and PERMANOVA), high accuracy values can be observed (Table 8). SVM achieved a 93.8% of accuracy on distal variables in metacarpals, 91.6% in anterior first phalanx with all the variables and 95.8% in third phalanx with all the measurements as well. Stacking method reached a 95.6% in posterior first phalanx with all measurements and finally, Blending tests performed lower scores on distal (82%) and proximal (86%) measurements for posterior first phalanx and 76% on all the variables for anterior second phalanx.

Table 8

Details on the performance of Machine Learning algorithms according to habitat (open/close) as predictor variable.


BONE MEASUREMENTS ALGORITHM ACCURACY KAPPA 95%CI SENSITIVITY SPECIFICITY BALANCED ACCURACY

Metatarsus Distal Support Vector Machines 93,75% 0.87 0–69–0.99 1 0.9 0.95

Anterior First Phalanx All Support Vector Machines 91,6% 0.83 0.73–0.98 1 0.83 0.91

Posterior First Phalanx All Stacking 95,6% 0.91 0.78–0.99 0.91 1 0.95

Distal Blending 82% 0.64 0.61–0.95 0.91 0.72 0.82

Proximal Blending 86% 0.74 0.66–0.97 0.75 1 0.87

Anterior Second Phalanx All Blending 76% 0.54 0.50–0.93 1 0.55 0.77

Third Phalanx All Support Vector Machines 95.8% 0.91 0.78–0.99 1 0.92 0.96

Archaeological Example

This methodology can be applied on Upper Palaeolithic archaeological sites from South-western France with the aim of tracing mobility patterns and habitat type from fossil reindeer remains in the future. Here, we present an explicit application of ML algorithms with a small sample of metapodial bones randomly selected from our archaeological corpus.

Tables 9 and 10 show the results of applying the algorithms with accuracies greater than 90% to a testing data set composed of 29 reindeer metapodial bones (16 metacarpals and 13 metatarsals) from the archaeological site of Enlène (Ariège, Pyrénées, France, see Map.1-Supplementary Information) from Middle Palaeolithic (Beds 1 and 3 13760±70 – 13620±70) (Bégouën et al. 2019). These remains come from “La Salle du Fond.” Although we know that both adult and juvenile carcasses were hunted by the Magdalenians in Enlène (Bégouën et al. 2019), only the adult specimens were included in the study. XgbTree for proximal section (Bp and Dp measurements) in metacarpals (Accuracy 95.2%) and Neural Network (Accuracy 90%) for metacarpal distal measurements (Bd,Dd, TMLMIN) to evaluate mobility pattern and SVM for metatarsal distal measurements (Bd,Dd, TMLMIN) to test mobility pattern (Accuracy 91.3%) and habitat type (Accuracy 93.8%).

Table 9

Classification results for proximal metacarpal according to XgbTree algorithm.


BONE TYPE BONE REFERENCE ALGORITHM MIGRATORY SEDENTARY CLASSIFICATION

Mtc proximal R11e–874 XgbTree 0.9418254 0.05817461 Migratory

Mtc proximal S6d–716 XgbTree 0.9418254 0.05817461 Migratory

Mtc proximal Q12f–972 XgbTree 0.9418254 0.05817461 Migratory

Mtc proximal N8d–1237 XgbTree 0.9418254 0.05817461 Migratory

Table 10

Classification results for distal metacarpal (mobility pattern) and distal metatarsal (Mobility pattern and Habitat) according to SVM and NN algorithms. See Results section to key accuracies and Kappa values.


BONE TYPE BONE REFERENCE ALGORITHM MOBILITY PATTERN CLASSIFICATION HABITAT CLASSIFICATION

Mtc distal T6c–200 NN Migratory

Mtc distal O9d–658 NN Sedentary

Mtc distal N10h NN Migratory

Mtc distal Q12d–591 NN Migratory

Mtc distal O8e inf–1837 NN Migratory

Mtc distal Q12f–971 NN Migratory

Mtc distal R11c–441 NN Migratory

Mtc distal O8c–717 NN Sedentary

Mtc distal S5d–590 NN Migratory

Mtc distal M8e–484 NN Migratory

Mtc distal M8e–474 NN Migratory

Mtc distal Q9d–763 NN Migratory

Mtt distal P9c–116 SVM Migratory Open

Mtt distal R10c–494 SVM Migratory Open

Mtt distal T6c–163 SVM Migratory Open

Mtt distal R11e–1745 SVM Migratory Open

Mtt distal O8e inf–1745 SVM Migratory Open

Mtt distal O8e inf–1824 SVM Migratory Open

Mtt distal O9c–253 SVM Migratory Open

Mtt distal S10f inf–95 SVM Migratory Open

Mtt distal T11f–7 SVM Migratory Open

Mtt distal O8d–1467 SVM Migratory Open

Mtt distal N8c–387 SVM Migratory Open

Mtt distal N7c–176 SVM Migratory Open

Mtt distal M8c –57 SVM Migratory Open

In this way, regarding mobility pattern, 27 out of 29 metapodial bone remains (both proximal and distal) were classified as migratory and only 2 distal metacarpals as sedentary using Extreme gradient boosting Tree, Neural Network and Support Vector Machines with respect to habitat type. Support Vector Machines classified all distal metatarsals as “Open”.

Discussion

As indicated above, it is important to bear in mind that the definition of “sedentary” provided here is far from its meaning in an archaeological context. For some authors (Bouchud et al. 1953, Bouchud 1966, Deplano 1994, Fontana 2000, 2012, 2017), fossil reindeer antlers from La Madeleine (Bonnisent, 1993) and bone remains from layer IX of the Flageolet II (Deplano, 1994) reveal an annual presence of reindeer in the Périgord area; its sedentary lifestyle is linked to environmental characteristics (Fontana, 2000).

Research on entheseal changes in phalanges of Fennoscandian reindeer (Rangifer tarandus tarandus and Rangifer tarandus fenicus) provides a valuable insight into changes related to reindeer activity (Niinimäki and Salmi, 2016), behaviour and foraging strategies (Hull et al. 2020) that occur in these populations, with important implications from the zooarchaeological and palaeoecological point of view. For example, these entheseal changes between R.t.tarandus and R.t.fenicus are more pronounced in first and third phalanx. Consequently, measurements in our study were taken outside of the entheses zones to prevent any influence on them.

Locomotor adaptations are linked to mobility and habitat preferences; thus, an animal’s postcranial morphology allows us to infer these attributes through actualistic approaches. Several authors using different taxa, such as bovids (Plummer and Bishop, 1994; DeGusta and Vrba, 2003; 2005) or suids (Bishop, 1994; 1999), have demonstrated that it is possible to study a taxon’s ecomorphology in the past, even when these are different from their modern descendants. Analyzing the ecological preferences of extant fauna, the present work contributes to testing hypotheses about the relationship of taxa to habitats and mobility in the past. Hence, and according to the results of the present study, linear measurements in combination with Machine Learning algorithms allow us to infer mobility pattern (migratory/sedentary) and habitat type (open/closed). To go further, and infer more specific habitats or degrees of mobility, a larger sample is needed.

The combination of linear measurements and discriminant function analysis (DFA) has been traditionally used to classify and infer habitat type from unknown archaeological specimens (Bishop, 1994; DeGusta and Vrba, 2003, 2005a, 2005b; Kappelman, 1988; Kovarovic and Andrews, 2007, Klein et al. 2010). However, a high rate of reclassification has not always been achieved. Here, we go a step further looking for a high level of accuracy to infer both mobility pattern and habitat type.

As with other research (Pelletier et al. 2020) on reindeer morphological changes, it has been proven that forelimb long bones and phalanges provide information on changes in locomotor behavior linked to mobility. Through 3D GMM, Pelletier et al. (2020) identified incipient domestication markers in reindeer. They observed how individuals bred in captivity had smaller bone elements and a thinner and more slender morphology than free-ranging individuals, linked to control and reduction of mobility and food by humans. DeGusta and Vrba (2005), in relation to habitat type, observed in bovids that Forest taxa had longer phalanges relative to their width, while Open-cover taxa presented shorter ones in relation to width. We have observed that longer metapodial bones (GL) are linked to sedentary individuals and closed habitats (see Supplementary Information). Similarly, longer phalanges correspond to closed habitats and shorter phalanges to open habitats.

Ecomorphological studies on suids (Bishop, 1994) from Pliocene and Pleistocene hominid sites in East Africa revealed that they preferred open, intermediate and closed habitats, inferring that these three habitats were available at these sites. As a result, algorithms and linear measurements will enable us to analyze Magdalenian faunal assemblages and infer the presence of Rangifer tarandus herds on open/closed or migratory/sedentary mobility to evaluate the implications for Upper Palaeolithic human groups on habitat use and hunting strategies.

As was pointed out by Bishop (1994), fossil bones are often incomplete, and it is important to analyze epiphyses since they preserve some locomotor anatomy. In the same way and in accordance with Kovarovic and Andrews (2007), there are several advantages in surveying a number of elements in this kind of osteometric analysis. On the one hand there is the possibility to increase the sample size when the method is applied to the archaeological record. On the other, there is the chance to compare one element with others to ensure that they provide aligned results and offer a bigger picture to show trends more effectively. That is why we emphasize the importance of obtaining high accuracy values using only a few measurements in combinations with algorithms (for example, to evaluate mobility on proximal metacarpals, we obtained a successful 95.2% of accuracy using only two linear measurements), and the use of several bone sections (if possible) to cross different results.

In accordance with Bishop (1999) and Kappelman (1988), the relationship of locomotor anatomy with substrate type can be explored in modern animals. Thus, the method presented above provides a useful tool for inferring fossil reindeer habitat type (open/closed) and two major mobility patterns (mobile reindeer, which cover distances greater than 200 kilometres, and less migratory reindeer that move less than 200 kilometres). As a hypothesis, the use of extant Rangifer tarandus to infer a more precise mobility degree (e.g. altitudinal movements) on fossil reindeer (in this case from Magdalenian period) and habitat type may be possible by using 3D Geometric Morphometrics techniques. This is because linear measurements can be insufficient to capture the geometry of the original object and some aspects of the shape can be lost (Adams et al. 2004; Zelditch et al. 2004). In this regard, a comparison between both methodologies could turn out to be interesting.

Principal Component Analysis on metacarpal bones shows a clear separation in Rangifer tarandus according to mobility pattern while on metatarsals this separation is more evident with regard to an open/closed habitat dichotomy. Moreover, sedentary individuals have longer shafts and a wider distal condyle than migratory individuals (see Supplementary Information), whereas these are lighter in keeping with other studies (Courturier et al. 2010).

For Machine Learning analyses, metacarpals produced higher accuracy than metatarsals for mobility pattern classification, which is probably linked to the fact that forelimbs carry 60% of static body weight and absorb ground impact as the body is thrown forward by the hindlimbs (Fletcher, 2019). However, we obtained better results for metatarsal bones when habitat is tested as a predictor. Thus, metacarpal proximal variables and XgbTree achieved the best accuracy to predict mobility pattern.

Regarding phalanges, it is possible to successfully classify migratory/sedentary patterns using posterior second phalanges and NN algorithm with more than 93% accuracy when all the variables are taken. However, the second phalanx turns out to be a weak classifier if we use habitat as a predictor and mobility on distal measurements, having the lowest classification rates (76%). Both anterior and posterior first phalanges allow us to classify both habitat type and mobility pattern with more than 90% accuracy when all the measurements are taken. Finally, third phalanges using all the variables and XgbTree algorithm show quite good results (92% accuracy) in order to infer mobility pattern and specially, open/closed habitat type (95.8%) through Support Vector Machines. These last results were expected, since third phalanx is the place of contact between the animal and the substrate (Curran, 2012).

Furthermore, high “Kappa” metric values support the reliability of the different models presented here, due to the strong level of agreement produced by the algorithms. High specificity and sensitivity values reinforce that, because of the very low possibilities to get Error type I and II. Although our reference sample is unbalanced (more migratory than sedentary individuals), this may not be a problem according to Balanced accuracy value. However, as mentioned above, a larger sample will be needed for further analyses.

For the last few decades, many studies have focused on the role on reindeer in Upper Palaeolithic subsistence systems, particularly during Magdalenian, applying very different techniques and methods. Long migrations (North-South axis) and herd following theory are proposed based on ethnographical models, dental cemento-chronology and antler analyses (Lacorre, 1956; Bahn, 1977; Gordon 1988a, 1988b, 1990). On the opposite side, we find those who support reindeer having been completely sedentary in the Périgord region through the study of antlers and bone assemblages (Bouchoud et al. 1953; Bouchoud, 1954; Deplano, 1994; Fontana, 2000; 2017). Finally, relatively short migrations along an East to West axis based on osteometric and dental data as well as seasonality studies (Delpech, 1983, 1987; Kuntz and Costamagno, 2011) are suggested. Seasonal data that have multiplied since the 1990s seem to indicate the absence of large-scale migrations on a North-South axis (Kuntz, 2011).

Our results show that subspecies structure the data significantly, but habitat and mobility are also significant predictors of morphological difference. That was expected since the very different habitats of modern Canadian caribou creates genetic isolation and phenotypic differences arise. Thus, it has been proven and supported by statistical evidence that both habitat and mobility affect bone morphology and this does not occur for behavioural reasons. Moreover, as discussed above, habitat types (especially closed/open vegetation types) are closely linked to migration distances covered and tundra caribou are invariably mobile whilst boreal forest and mountain caribou can use a variety of different habitat types while travelling only short distances. In other words, two of our predictors interact. Nevertheless, we are able to distinguish between long and short migration distances and, in some instances, closed and open habitats and this ecomorphological information will be useful for archaeological analyses we intend to carry out in future.

This methodology developed here can be applied on Upper Palaeolithic archaeological sites from South-western France with the aim to trace mobility patterns and habitat type from fossil reindeer remains to explore how they affected human hunting strategies and socio-economic decisions, looking for a better understand of their behaviour and identifying the precise role of reindeer in Magdalenian economy.

Conclusion

Here we have shown that osteometric data in combination with Machine Learning algorithms can successfully classify modern caribou assemblage (and create a referential framework) according to mobility pattern (migratory/sedentary) and habitat (open/closed) with an accuracy more than 95%, irrespective of subspecies with the use of a few linear measurements. This classification system, carried out with several algorithms, avoids any bias that can be produced in the analysis and therefore offers a way to test and classify reindeer mobility in the past. Thus, proximal and distal measurements on metacarpals, distal metatarsals, and all the measurements on anterior first phalanx, posterior second phalanx, and third phalanx, provide a successful means to predict mobility pattern. Furthermore, third phalanx, posterior first phalanx, and the metatarsal distal end are very good predictors regarding habitat type. Moreover, sedentary individuals have longer metacarpals and a wider distal condyle than migrating individuals, while migrating individuals are slenderer, a result in agreement with other studies (Courturier et al. 2010).

It has been proven that the often fragmentary condition of most bones in archaeological assemblages, should not be an obstacle to reconstruction of the mobility pattern, as it is possible to classify using separate proximal and distal measurements. Hence, these results are very encouraging for future archaeological applications.

The results of our study illustrate how linear measurements of metapodial bones and phalanges can be used to infer reindeer mobility and habitat type as locomotor adaptations are in turn linked to habitat preferences. Therefore, these useful methodologies are good starting point in the creation of a frame of reference to successfully determine the type of reindeer mobility in the Magdalenian that will allow us a better understanding of human hunting strategies and identify the precise role of reindeer in their socio-economic decisions.

Additional Files

The additional files for this article can be found as follows:

Supplementary material

Suplementary Information. DOI: https://doi.org/10.5334/oq.106.s1

Dataset

Raw measurements of metapodial bones, first, second and third phalanges as well as the R Script used in the present study. DOI: https://doi.org/10.5334/oq.106.s2

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 Research and Innovation Program (Marie Sklodowska Curie Actions, Grant Agreement n ° 794925) and Marie Curie Alumni Association through a Micro-Media Grant.

We would like to thank Claire Saint-Germain (Ostéoteque de l’Université de Montréal, Quebec, Canada), Maeve Winchester (Ministry of Forests, Lands, Natural Resources Operations and Rural Development of British Columbia, Canada), Kamal Khidas (Canadian Museum of Nature, Ottawa, Canada), Guillaume Gingras (Forest and Wildlife Department of the Government of Quebec, Canada) and James Woollett and Steeve Côté (Laval University, Quebec, Canada) who provided us with samples of caribou to complete this study. We also thank the reviewers of this manuscript. ABGL would like to thank A. Martín Esteban, S. Galán, N. Lupiañez Corpas and M. Domínguez-Rodrigo for their valuable comments and support.

Competing Interests

The authors have no competing interests to declare.

References

  1. Arceredillo, D, Gómez-Olivencia, A and García-Pérez, A. 2011. Three statistical methods for sex determination in extant and fossil carprines: assessment of the Rupicapra long bones. Journal of Archaeological Science, 38: 2450–2460. DOI: https://doi.org/10.1016/j.jas.2011.05.013 

  2. Bahn, PG. 1977. Seasonal migration in South-West France during the Late Glacial Period. Journal of Archaeological Science, 4: 245–257. DOI: https://doi.org/10.1016/0305-4403(77)90092-9 

  3. Banks, WE, D’Errico, F, Peterson, AT, Kageyama, M and Colombeau, G. 2008. Reconstructing ecological niches and geographic distributions of Rangifer tarandus and Cervus elaphus during the LGM. Quaternary Science Reviews, 27: 2568–2575. DOI: https://doi.org/10.1016/j.quascirev.2008.09.013 

  4. Bergerud, AT. 2000. Caribou. In: Demaris, S and Krausman, PR (eds.), Ecology and management of large mammals in North America, 658–693. New Jersey, Prentice-Hall. 

  5. Berteaux, D and Guintard, C. 1995. Osteometric study of the metapodials of Amsterdam Island feral cattle. Acta Theriologica, 40(1): 97–110. DOI: https://doi.org/10.4098/AT.arch.95-10 

  6. Bishop, L. 1994. Pigs and the ancestors: hominids, suids and environments during the Plio-Pleistocene of East Africa. Ph.D. Dissertation, Yale University. 

  7. Bishop, LC. 1999. Suid paleoecology and habitat preferences at African Pliocene and Pleistocene hominid localities. In: Schrenk, TG (ed.), African Biogeography, Climate Change, & Human Evolution. New York: Oxford University Press. 

  8. Bonnisent, D. 1993. Choix et exploitation des bois de Rennes sur le site de la Madeleine (Dordogne), Memoire de D.E.A., Universit. de Bordeaux I, 48 p. 

  9. Bouchud, J. 1954. Dents de Rennes, Bois de Rennes et migrations. Bulletin de la Societé Prehistorique Francaise, 51: 340–345. DOI: https://doi.org/10.3406/bspf.1954.3108 

  10. Bouchud, J. 1966. Essai sur le Renne et climatologie du Paleolithique moyen et superieur. Imprimerie Magne, Perigueux, 297. 

  11. Bouchud, J, Cheynier, A and Guillien, Y. 1953. Dents de Renne et migrations. Bulletin de la Societé Prehistorique Francaise, 50: 127–132. DOI: https://doi.org/10.3406/bspf.1953.3014 

  12. Breiman, L. 2001. Random forests. Machine Learning 45: 5–32. DOI: https://doi.org/10.1023/A:1010933404324 

  13. Castaños, JP, Murelaga, X, Alonso-Olazabal, A, Ortega, LA and Zuluaga, MC. 2012. Osteometric Analysis of the Scapula and Humerus of Rangifer tarandus and Cervus elaphus: a Contribution to the Discrimination of Late Pleistocene Cervids. Acta Palaeontologica Polonica, 59(4): 779–786. DOI: https://doi.org/10.4202/app.2012.0027 

  14. Cavedon, M, Gubili, C, Heppenheimer, E, Von Holdt, B, Mariani, S, Hebblewhite, M, Hegel, T, Hervieux, D, Serrouya, R, Steenwegs, R, Weckworth, B and Mussiani, M. 2019. Genomics, environment and balancing selection in behaviourally bimodal populations: The caribou case. Molecular Ecology, 28: 1946–1963. DOI: https://doi.org/10.1111/mec.15039 

  15. Cifuentes-Alcobendas, G and Dominguez-Rodrigo, M. 2019. Deep learning and taphonomy: high accuracy in the classification of cut marks made on fleshed and defleshed bones using convolutional neural networks. Scientific Reports, 9: 18933. DOI: https://doi.org/10.1038/s41598-019-55439-6 

  16. Clark, M. 2013. An Introduction to Machine Learning with applications in R. Indiana: Centre for Social Research, University of Notre Dame. 

  17. COSEWIC. 2002. Assessment and Update Status Report on the Woodland Caribou Rangifer tarandus Caribou in Canada. pp. xi + 98. Ottawa: Committee on the Status of Endangered Wildlife in Canada. 

  18. COSEWIC. 2004. Assessment and update status report on the Peary caribou Rangifer tarandus pearyi and the barren-ground caribou Rangifer tarandus groenlandicus (Dolphin and Union population) in Canada. Ottawa: Committee on the Status of Endangered Wildlife in Canada. x + 91 pp. 

  19. COSEWIC. 2011. Designatable Units for Caribou (Rangifer tarandus) in Canada. Ottawa: Committee on the Status of Endangered Wildlife in Canada. 

  20. COSEWIC. 2015. Assesment and Status Report on the Peary Caribou, Rangifer tarandus pearyi in Canada. Ottawa: Committee on the Status of Endangered Wildlife in Canada. 

  21. COSEWIC. 2016. Assessment and status report on the Caribou Rangifer tarandus, Barren-ground population, in Canada. Ottawa: Committee on the Status of Endangered Wildlife in Canada. 

  22. COSEWIC. 2017. Assessment and status report on the Caribou Rangifer tarandus, Eastern Migratory population and Torngat Mountains population, in Canada. Ottawa: Committee on the Status of Endangered Wildlife in Canada. 

  23. Costamagno, S, Barshay-Szmidt, C, Kuntz, D, Laroulandie, V, Pétillon, JM, Boudadi Maligne, M, Langlais, M, Maylle, JB and Chevallier, A. 2016. Reexamining the timing of reindeer disappearance in southwestern France in the larger context of late glacial faunal turnover. Quaternary International, 414: 34–61. DOI: https://doi.org/10.1016/j.quaint.2015.11.103 

  24. Costamagno, S and Mateos, A. 2007. Milieu animal de part et d’autre de la chaene pyreneenne: implication sur les modes de subsistance au Magdalenien. In: Cazals, N, Gonzalez Urquijo, J et X (eds.), Terradas, Frontieres naturelles et frontieres culturelles dans les Pyrenees prehistoriques, 53–73. Santander. 

  25. Courturier S, Otto, R, Côté, SD, Luther, G and Mahoney, SP. 2010. Body size variations in caribou cotypes and relationships with demography. Journal of Wildlife Management, 74: 395–404. DOI: https://doi.org/10.2193/2008-384 

  26. Cumming, HG and Beange, DB. 1987. Dispersion and movements of woodland caribou near Lake Nipigon, Ontario. Journal of Wildlife Management, 51: 69–79. DOI: https://doi.org/10.2307/3801634 

  27. Curran, SC. 2012. Expanding ecomorphological methods: geometric morphometric analysis of Cervidae post-crania. Journal of Archaeological Science, 39: 1172–1182. DOI: https://doi.org/10.1016/j.jas.2011.12.028 

  28. Davis, SJM. 1985. The large mammal bones. In: Mazar, A (ed.), Excavations at Tell Qassile 2. Qedem 20, Jerusalem, Hebrew University. 

  29. Davis, SJM. 1987. The archaeology of animals. London: Batsford. 

  30. Davis, SJM. 2000. The effect of castration and age on the development of the Shetland sheep skeleton and a metric comparison between bones of males, females and castrâtes. Journal of Archaeological Science, 27: 373–390. DOI: https://doi.org/10.1006/jasc.1999.0452 

  31. Davis, SJM, Svensson, EM, Albarella, U, Detry, C, Götherström, A, Pires, AE and Ginja, D. 2012. Molecular and osteometric sexing cattle metacarpals: a case of study from 15th century A.D. Beja, Portugal. Journal of Archaeological Science, 39: 1445–1454. DOI: https://doi.org/10.1016/j.jas.2011.12.003 

  32. Degusta, D and Vrba, E. 2003. A method for inferring paleohabitats from the functional morphology of bovid astragali. Journal of Archaeological Science, 30: 1009–1022. DOI: https://doi.org/10.1016/S0305-4403(02)00286-8 

  33. Degusta, D and Vrba, E. 2005a. Methods for inferring paleohabitats from discrete traits of the bovid prostcranial skeleton. Journal of Archaeological Science, 32: 1115–1123. DOI: https://doi.org/10.1016/j.jas.2005.02.011 

  34. Degusta, D and Vrba, E. 2005b. Methods for inferring paleohabitats from the functional morphology of bovid phalanges. Journal of Archaeological Science, 32: 1099–1113. DOI: https://doi.org/10.1016/j.jas.2005.02.010 

  35. Delpech, F. 1983. Les faunes du Paléolithique supérieur dans le Sud-Ouest de la France. Paris: Cahiers du Quaternaire 6. 

  36. Delpech, F. 1987. L’environnement animal des magdaléniens, In: Otte, M. (Éd.), Le Magdalénien en Europe, Liège, ERAUL, Actes du Colloque de Mayence, XI° congrès U.I.S.P.P., 5–30. 

  37. Deplano, S. 1994. Etude de la faune des grands mammifères de la couche IX de l’Abri du Flageolet3 Dordogne. Approche taphonomique et palethnographique. Memoire de Maitrise en Prehistoire, Universitéde Paris. 

  38. Domínguez-Rodrigo, M, Cifuentes-Alcobendas, G, Jiménez-García, B, Abellán, N, Pizarro-Monzo, M, Organista, E and Baquedano, E. 2020. Artificial intelligence provides greater accuracy of modern and ancient bone surface modifications. Scientific Reports, 10: 18862. DOI: https://doi.org/10.1038/s41598-020-75994-7 

  39. Edmonds, E. 1988. Population status, distribution, and movements of woodland caribou in west central Alberta. Canadian Journal of Zoology, 66: 817–826. DOI: https://doi.org/10.1139/z88-121 

  40. Enloe, JG. 2001. Seasonality, Strategy and Site Function: Reindeer Hunting at Verberie. Journal of Human Evolution, 40(3): A7–A8. 

  41. Fancy, SG, Pank, LF, Whittenk, KR and Regelin, WL. 1989. Seasonal movements of caribou in arctic Alaska as determined by satellite. Canadian Journal of Zoology, 67(3): 644–650. DOI: https://doi.org/10.1139/z89-093 

  42. Ferguson, MAD and Messier, F. 2000. Traditional winter range: Population Dynamics and Physical Condition. The Journal of Wildlife Management, 64: 168–178. DOI: https://doi.org/10.2307/3802987 

  43. Festa-Bianchet, M, Ray, JC, Côté, SD and Gunn, A. 2011. Conservation of caribou (Rangifer tarandus) in Canada: an uncertain future. Canadian Journal of Zoology, 89: 419–434. DOI: https://doi.org/10.1139/z11-025 

  44. Fletcher, TF. 2019. Designed for running. College of veterinary medicine. University of Minessota. 

  45. Fontana, L. 2000. La chasse au Renne au Paléolitique Supérieur: nouvelles voies de recherche. Páleo, 12: 141–164. DOI: https://doi.org/10.3406/pal.2000.1600 

  46. Fontana, L. 2012. Archéologie et anthropologie des relations homme-animal dans les sociétés anciennes de chasseurs-collecteurs. Les enjeux d’une recherche intégrée en archéozoologie, In: de Beaune, SA & Fracfort, HP (dir.), L’archéologie à découvert. Paris: CNRS Editions, 72–78. DOI: https://doi.org/10.4000/books.editionscnrs.11242 

  47. Fontana, L. 2017. The four seasons of reindeer: non-migrating reindeer in the Dordogne region (France) between 30 and 18k? Data from the Middle and Upper Magdalenian at La Madeleine and methods of seasonality determination. Journal of Archaeological Science: Reports, 12: 346–362. DOI: https://doi.org/10.1016/j.jasrep.2017.02.012 

  48. Fontana, L and Chaurière, FX. 2018. Economie et nomadisme au Pléniglaciaire supérieur et au Tardiglaciaire en Europe de l’ouest: le Systeme Renne. In: Djindjian, F (ed.), La préhistoire de la France, Harmann Editeurs, Histoire et archéologie, 141–148. DOI: https://doi.org/10.3917/herm.djind.2018.01.0141 

  49. Friedman, JH. 2001. Greedy function approximation: a gradient boosting machine. Annals of statistics, 1189–1232. DOI: https://doi.org/10.1214/aos/1013203451 

  50. Galán López, AB and Domínguez-Rodrigo, M. 2014. A biometric analysis of the pelvic acetabulum as an indicator of sex in bovids. Comptes Rendus Palevol, 13: 561–567. DOI: https://doi.org/10.1016/j.crpv.2014.04.003 

  51. Geist, V. 1998. Deer of the World: Their Evolution, Behavior and Ecology. Mechanisburg: Stackpole Press. 

  52. Goddard, AD. 2009. Boreal caribou in Northeastern British Columbia. Peace Region Technical Report, 1–27. 

  53. González-Molina, I, Jiménez-García, B, Maillo-Fernández, J-M, Baquedano, E and Domínguez-Rodrigo, M. 2020. Distinguishing Discoid and Centripetal Levallois methods through machine learning. PLoS ONE, 15(12): e0244288. DOI: https://doi.org/10.1371/journal.pone.0244288 

  54. Gordon, BC. 1988a. Des Hommes et des Rennes dans la Préhistoire Française Magdalénienne: Résultats. Archaeozoologia, 2: 227–242. 

  55. Gordon, B. 1988b. Of Men and Reindeer Herds in French Magdalenian Prehistory, 390. Oxford: BAR International Series. DOI: https://doi.org/10.30861/9780860545040 

  56. Gordon, B. 1990. More on herd-following hypothesis. Current Anthropology, 31: 399–400. DOI: https://doi.org/10.1086/203858 

  57. Grayson, F and Delpech, F. 2002. Specialized early Upper Palaeolithic hunters in southwestern France? Journal of Archaeological Science, 29: 1439–1449. DOI: https://doi.org/10.1006/jasc.2002.0806 

  58. Greenfield, HJ. 2002. Sexing fragmentary Ungulate Acetabulae. In: Ruscillo, D (Ed.), Recent Advances in Ageing and Sexing Animal Bones, 68–86. Proceeding of the 9th ICAZ Conference. Durham. Oxford, UK. Oxbow Press. DOI: https://doi.org/10.2307/j.ctvh1ds02.9 

  59. Gunn, A. 2010. Migratory Tundra Caribou. In: Hummel, M and Ray, JC (eds.), Caribou and the North: a shared future, 240–267. Toronto: Dundurn Press. 

  60. Heard, DC and Vagt, KL. 1998. Caribou in British Columbia: a 1996 status report. Rangifer Special Issue, 10: 117–123. DOI: https://doi.org/10.7557/2.18.5.1548 

  61. Hegel, TM and Russel, K. 2013. Status of northern mountain Caribou (Rangifer tarandus caribou) in Yukon, Canada. Rangifer, 33: 59–70. DOI: https://doi.org/10.7557/2.33.2.2528 

  62. Hull, EH. 2019. Metric and non-metric guides for the determination between fore-and hindlimb phalanges of Rangifer tarandus. Rangifer, 39(1): 11–26. DOI: https://doi.org/10.7557/2.39.1.4630 

  63. Hull, E, Niinimäki, S and Salmi, AK. 2020. Differences in entheseal changes in the phalanges between ecotypes of Fennoscandian reindeer. International Journal of Osteoarchaeology, 30(5): 666–678. DOI: https://doi.org/10.1002/oa.2897 

  64. Hull, E, Semeniuk, M, Poulakka, H-L, Kynkäänniemi, S-M and Niinimäki, S. 2021. Tendons and ligaments of the Rangifer tarandus metapodial and hoof. Polar Biology, 44: 1803–1816. DOI: https://doi.org/10.1007/s00300-021-02919-z 

  65. Hummel, M and Ray, JC. 2010. Caribou and the North: a shared future. Toronto: Dundurn Press. 

  66. Jenkins, D, Campbell, M, Hope, G, Goorts, J and Mcloughlin, P. 2011. Recent trends in abundance of Peary Caribou (Rangifer tarandus pearyi) and Muskoxen (Ovibos moschatus) in the Canadian Arctic Archipelago, Nunavut. Department of Environment, Government of Nunavut, Wildlife Report No. 1, Pond Inlet, Nunavut. 

  67. Joly, K, Gurarie, E, Sorum, MS, Kaczensky, P, Cameron, MD, Jakes, AF, Borg, BL, Nandintsetseg, D, Hopcraft, JG, Buuveibaatar, B, Jones, PF, Mueller, T, Walzer, C, Olson, KA, Payne, JC, Yadamsuren, A and Hebblewhite, M. 2019. Longest terrestrial migrations and movements around the world. Scientific Reports, 9: 15333. DOI: https://doi.org/10.1038/s41598-019-51884-5 

  68. Kappelman, J. 1988. Morphology and locomotor adaptations of the bovid femur in relation to habitat. Journal of Morphology, 198: 119–130. DOI: https://doi.org/10.1002/jmor.1051980111 

  69. Kappelman, J, Plummer, T, Bishop, L, Duncan, A and Appleton, S. 1997. Bovids as indicators of Plio-Pleistocene paleoenvironments in East Africa. Journal of Human Evolution, 32: 229–256. DOI: https://doi.org/10.1006/jhev.1996.0105 

  70. Klein, RG, Francicus, RG and Steele, TE. 2010. Morphometric identification of bovid metapodials to genus and implications for taxon-free habitat reconstruction. Journal of Archaeological Science, 37: 389–401. DOI: https://doi.org/10.1016/j.jas.2009.10.001 

  71. Kovarovic, K. 2004. Bovids as palaeoenvironmental indicators. An ecomorphological analysis of bovid postcranial remains from Laetoli, Tanzania. PhD diss., University of London. 

  72. Kovarovic, K and Andrews, P. 2007. Bovid postcranial ecomorphological survey of the Laetoli paleoenvironments. Journal of Human Evolution, 52: 663–680. DOI: https://doi.org/10.1016/j.jhevol.2007.01.001 

  73. Kuhn, M and Johnson, K. 2013. Applied Predictive Modelling. New York: Springer. DOI: https://doi.org/10.1007/978-1-4614-6849-3 

  74. Kuntz, D. 2011. Ostéométrie et migration(s) du renne (Rangifer tarandus) dans le Sud-Ouest de la France au cours du dernier Pléniglaciaire et du Tardiglaciaire (21 500–13 000 cal. BP). PhD: Universite de Toulouse le Mirail – Toulouse II. 

  75. Kuntz, D and Costamagno, S. 2011. Relationships between reindeer and man in southwestern France during the Magdalenian. Quaternary Research, 238: 12–24. DOI: https://doi.org/10.1016/j.quaint.2010.10.023 

  76. Lacorre, F. 1956. Les Migrations des Rennes dans la Province Préhistorique des Eyzies. Bulletin de la Societe Prehistorique Francaise, 53: 302–310. DOI: https://doi.org/10.3406/bspf.1956.3336 

  77. Lantz, B. 2013. Machine Learning with R. Birmingham: Packt Publishing. 

  78. Lartet, E and Christy, H. 1864. Cavernes du Périgord. Objets gravés et sculptés des temps préhistoriques dans l’Europe occidentale. Revue Archéologique, Paris. 

  79. Lartet, E and Christy, H. 1875. Reliquae Aquitanicae: Being Contributions to the Archaeology and Paleontology of Perigord and Adjoining Provinces of Southern France. London: Williams and Norgate. DOI: https://doi.org/10.5962/bhl.title.25391 

  80. Leblond, M, St Laurent, MH and Côté, SD. 2016. Caribou, water and ice-fine-scale movements of a migratory arctic ungulate in the context of climate change. Movement Ecology, 4: 14. DOI: https://doi.org/10.1186/s40462-016-0079-4 

  81. Lefebvre, A, Rochefort, GY, Santos, F, Le Denmat, D, Salmon, B and Pétillon, J-M. 2016. A Non-Destructive Method for Distinguishing Reindeer Antler (Rangifer tarandus) from Red Deer Antler (Cervus elaphus) Using X-Ray Micro-Tomography Coupled with SVM Classifiers. PLoS ONE, 11: E0149658. DOI: https://doi.org/10.1371/journal.pone.0149658 

  82. Mallory, FF and Hillis, TL. 1998. Demographic characteristics of circumpolar caribou populations: ecotypes, ecological constraints/releases, and population dynamics. Rangifer, 9–60. DOI: https://doi.org/10.7557/2.18.5.1541 

  83. Mellars, PA. 2004. Reindeer specialization in the early Upper Palaeolithic: the evidence from south west France. Journal of Archaeological Science, 31: 613–617. DOI: https://doi.org/10.1016/j.jas.2003.10.010 

  84. Miller, FL. 1990. Inter-island movements of Peary caribou: a review and appraisement of their ecological importance. In: Harington, CR (ed.), Canada’s missing dimension: science and history in the Canadian Arctic Islands, 2: 608–632. Ottawa, Ontario: Canadian Museum Nature. 

  85. Miller, FL. 2002. Multi-island seasonal home range use by two Peary Caribou, Canadian High Arctic Islands, Nunavut, 1993–94. Arctic 55: 133–142. DOI: https://doi.org/10.14430/arctic697 

  86. Miller, FL. 2003. Caribou. In: Wild Mammals of North America: Biology, Management and Conservation. Feldhamer, GA & Thompson, BC & Chapman, JA (eds.). Baltimore: Johns Hopkins University Press. 

  87. Miller, FL and Gunn, A. 1978. Inter-island movements of Peary Caribou south of viscount Melville sound, Northwest Territories. Canadian Wildlife Service, 92: 327–333. 

  88. Morin, E. 2008. Evidence for declines in human population densities during the early Upper Paleolithic in Western Europe. Proceedings of the National Academy of Sciences, 105(1): 48–53. DOI: https://doi.org/10.1073/pnas.0709372104 

  89. Munro, ND, Baroz, G and Mill, AC. 2011. An exploration of character traits and linear measurements for sexing mountain gazelle (Gazella gazella) skeletons. Journal of Archaeological Science, 38: 1253–1265. DOI: https://doi.org/10.1016/j.jas.2011.01.001 

  90. Niinimäki, S and Salmi, AK. 2016. Entheseal changes in free-ranging versus zoo reindeer—Observing activity status of reindeer. International Journal of Osteoarchaeology, 26(2): 314–323. DOI: https://doi.org/10.1002/oa.2423 

  91. Nomokonova, T, Losey, RJ, Kosintsey, PA and Plekhanov, AV. 2021. Reindeer Demographics at Iarte VI, Iamal Peninsula, Arctic Siberia. Environmental Archaeology, 1–10. DOI: https://doi.org/10.1080/14614103.2021.1995260 

  92. Nti, IK, Adekoya, AF and Asubam, B. 2020. A comprehensive evaluation of ensemble learning for stock-market prediction. Journal of Big Data, 7: 20. DOI: https://doi.org/10.1186/s40537-020-00299-5 

  93. Nuttal, M. 2004. Encyclopedia of the Arctic. New York: Routledge. 

  94. Pelletier, M, Kotiaho, A, Niinimäki, S and Salmi, A-K. 2020. Identifying early stages of reindeer domestication in the archaeological record: a 3D morphological investigation on forelimb bones of modern populations from Fennoscandia. Archaeological and Anthropological Sciences, 12: 169. DOI: https://doi.org/10.1007/s12520-020-01123-0 

  95. Pelletier, M, Kotiaho, A, Niinimäki, S and Salmi, A-K. 2021b. Impact of selection and domestication on hindlimb bones of modern reindeer populations: Archaeological implications for early reindeer management by Sámi in Fennoscandia. Historical Biology. DOI: https://doi.org/10.1080/08912963.2021.1947268 

  96. Pelletier, M, Niinimäki, S and Salmi, A-K. 2021a. Influence of captivity and selection on limb long bone cross-sectional morphology of reindeer. Journal of Morphology, 282: 1533–1556. DOI: https://doi.org/10.1002/jmor.21403 

  97. Plummer, TW and Bishop, LC. 1994. Hominid paleoecology at Olduvai Gorge, Tanzania, as indicated by antelope remains. Journal of Human Evolution, 27: 47–75. DOI: https://doi.org/10.1006/jhev.1994.1035 

  98. Plummer, TW, Bishop, LC and Hertel, F. 2008. Habitat preference of extant African bovids based on astragalus morphology: operationalizing ecomorphology for palaeoenvironmental reconstruction. Journal of Archaeological Science, 35: 3016–3027. DOI: https://doi.org/10.1016/j.jas.2008.06.015 

  99. Prummel, W and Frisch, HJ. 1986. A guide for the distinction of species, sex, and body side in bones of sheep and goat. Journal of Archaeological Science, 13: 567–577. DOI: https://doi.org/10.1016/0305-4403(86)90041-5 

  100. Puputti, AK and Niskanen, M. 2008. The estimation of body weight of the reindeer (Rangifer tarandus L.) from skeletal measurements: preliminary analyses and application to archaeological material from 17th- and 18th-century northern Finland. Environmental Archaeology, 13: 153–164. DOI: https://doi.org/10.1179/174963108X343272 

  101. Purdue, JR. 1987. Estimation of body weight of white-tailed deer (Odocoileus virginianus) from central Illinois. Journal of Ethnobiology, 7: 1–12. 

  102. Reese, NE. 2015. The Ecomorphology of White-tailed Deer Lower Limb Bones Through the Holocene in Central North America. Master of Science Thesis. Southern Illinois University Edwardswille. 

  103. RSTUDIO TEAM. 2020. RStudio: Integrated Development for R. Boston, MA: RStudio, PBC. http://www.rstudio.com/. 

  104. Saher, DJ and Schmiegelow, FK. 2005. Movement pathways and habitat selection by woodland caribou during spring migration. Rangifer, 25: 143–154. DOI: https://doi.org/10.7557/2.25.4.1779 

  105. Salmi, AK, Niinimäki, S and Pudas, T. 2020. Identification of working reindeer using palaeopathology and entheseal changes. International Journal of Paleopathology, 30: 57–67. DOI: https://doi.org/10.1016/j.ijpp.2020.02.001 

  106. Schaefer, JA, Bergman, CM and Luttich, SN. 2000. Site fidelity of female caribou at multiple spatial scales. Landscape Ecology, 15: 731–739. DOI: https://doi.org/10.1023/A:1008160408257 

  107. Scott, KM. 1985. Allometric trends and locomotor adaptations in the Bovidae. Bulletin of the AMNH; v. 179, article 2. 

  108. Scott, RS. 2004. The comparative paleoecology of late Miocene Eurasian hominoids. Ph.D. Dissertation, The Unversity of Texas, Austin, Texas, USA, 457p. 

  109. Seip, D and Mclellan, B. 2010. Mountain Caribou. In: Hummel, M and Ray, JC (eds.), Caribou and the North: a shared future, 289–306. Toronto: Dundurn Press. 

  110. Sorensen, K, Mühldorff, R and Petersen, EB. 2007. The Scandinavian reindeer (Rangifer tarandus L.) after the last glacial maximun: time, seasonality and human exploitation. Journal of Archaeological Science, 34: 914–923. DOI: https://doi.org/10.1016/j.jas.2006.09.014 

  111. Spencer, L. 1997. Dietary adaptations of Plio-Pleistocene Bovidae: implications for hominid habitat use. Journal of Human Evolution, 32: 201–228. DOI: https://doi.org/10.1006/jhev.1996.0102 

  112. Spiess, A. 1979. Reindeer and caribou hunters – an archaeological study. New York: Academic Press. 

  113. Stuart-Smith, AK, Bradshaw, CJA, Boutin, S, Hebert, DM and Rippin, AB. 1997. —Woodland caribou relative to landscape patterns in northeastern Alberta. Journal of Wildlife Management, 61: 622–633. DOI: https://doi.org/10.2307/3802170 

  114. Taillön, J. 2013. Caribou. In: Prescott, J, Ferron, J and Taillön, J (eds.), Sur la piste de nos cervidés: orignal, cerf de Virginie, caribou. La Macaza, Québec: Orinha. 

  115. Thomas, DC and Gray, DR. 2002. Assessment and update status report on the Woodland Caribou Rangifer tarandus caribou in Canada. Ottawa: Committee on the Status of Endangered Wildlife in Canada. 

  116. Trevor, H, Tibshirani, R and Friedman, J. 2008. The Elements of Statistical Learning. New York: Springer. 

  117. van Asperen, EN. 2010. Ecomorphological adaptations to climate and substrate in late Middle Pleistocene caballoid horses. Palaeogeography, palaeoclimatology, palaeocology, 297: 584–596. DOI: https://doi.org/10.1016/j.palaeo.2010.09.007 

  118. von den Driesch, A. 1976. A guide to the measurement of animal bonesfrom archaeological sites. Bulletin 1. Cambridge: Peabody Museum, Harvard University. 

  119. Wareing, K, Tickle, PG, Stokkan, KA, Codd, JR and Sellers, WI. 2011. The musculoskeletal anatomy of the reindeer (Rangifer tarandus): fore and hind limbs. Polar Biology, 1–8. DOI: https://doi.org/10.1007/s00300-011-1017-y 

  120. Weaver, JL. 2006. Big Animals and Small Parks: Implications of Wildlife Distribution and Movements for Expansion of Nahanni National Park Reserve. Wildlife Conservation Society Canada Conservation Report No1.. Toronto, Ontario. 

  121. Weinstock, J. 2000a. Late Pleistocene reindeer populations in Middle and Western Europe, an osteometrical study of Rangifer tarandus. BioArchaeologica, 3. 

  122. Weinstock, J. 2000b. Demography through osteometry: sex ratios of reindeer and hunting strategies in the Late Glacial site of Stellmoor, northern Germany. Archaeozoologia, 11: 187–198 

  123. Weladji, RB, Klein, RG, Holand, O and Mysterud, A. 2002. Comparative response of Rangifer tarandus and other northern ungulates to climatic variability. Rangifer, 22: 33–50. DOI: https://doi.org/10.7557/2.22.1.686 

  124. White, R. 1980. The Upper Palaeolithic occupation of the Perigord: a topographic approach to subsistence and settlement. PhD University of Alberta. Alberta. 

  125. Wittmer, HU, Mclellan, BN and Hovey, FW. 2006. Factors influencing variation in site fidelity of woodland caribou (Rangifer tarandus caribou) in southeastern British Columbia. Canadian Journal of Zoology, 84: 537–545. DOI: https://doi.org/10.1139/z06-026 

  126. Wolpert, DH. 1992. Stacked generalization. Neural networks, 5(2): 241–259. DOI: https://doi.org/10.1016/S0893-6080(05)80023-1 

  127. Zeder, MA. 2001. A metrical analysis of a collection of modern goats (Capra hircus aegargus and Capra hircus hircus) from Iran and Iraq: implications for the study of Caprine domestication. Journal of Archaeological Science, 28: 61–79. DOI: https://doi.org/10.1006/jasc.1999.0555 

  128. Zhang, LG, Han, D, Pang, H, Yu, G, Cao, Q, Wang, C, Kong, L, Chengjin, W, Dong, W, Li, T and Li, J. 2020. Fore-limb joints contribute to locomotor performance in reindeer (Rangifer tarandus) by maintaining stability and storing energy. PeerJ, 8: 10278. DOI: https://doi.org/10.7717/peerj.10278 

  129. Zhang, LG, Qiao, Y, Ji, QL, Ma, SS and Li, JQ. 2017. Macro-microscopic research in reideer (Rangifer tarandus) hoof suitable for efficient locomotion on complex grounds. Nephron Clinical Practice, 61(2): 223–229. DOI: https://doi.org/10.1515/jvetres-2017-0029