# Reconstructing past migratory behaviour of reindeer
(Rangifer tarandus): Insights from geometric
morphometric analysis of proximal phalanx morphology
from extant caribou populations

Ana Belén Galán López, Maxime Pelletier, Emmanuel Discamps

► To cite this version:

Ana Belén Galán López, Maxime Pelletier, Emmanuel Discamps. Reconstructing past migratory behaviour of reindeer (Rangifer tarandus): Insights from geometric morphometric analysis of proximal phalanx morphology from extant caribou populations. PLoS ONE, 2023, 18 (8), e0285487.

HAL Id: halshs-04183926

https://shs.hal.science/halshs-04183926

Submitted on 21 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Reconstructing past migratory behaviour of reindeer (*Rangifer tarandus*): Insights from geometric morphometric analysis of proximal phalanx morphology from extant caribou populations

Ana Belén Galán López¹*ªa, Maxime Pelletier²ªb, Emmanuel Discamps¹

1 TRACES UMR 5608, CNRS-Université de Toulouse-Jean Jaurès, Toulouse, France, 2 Archaeology, History, Culture and Communication Studies, University of Oulu, Oulu, Finland

ªa Current address: TRACES UMR 5608, CNRS, Université de Toulouse-Jean Jaurès, Maison de la Recherche, Toulouse, France

ᵇb Current address: Archaeology, History, Culture and Communication Studies, Faculty of Humanities, University of Oulu, Oulu, Finland

* ana.galan-lopez@univ-tlse2.fr

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

## OPEN ACCESS

**Citation:** Galán López AB, Pelletier M, Discamps E (2023) Reconstructing past migratory behaviour of reindeer (*Rangifer tarandus*): Insights from geometric morphometric analysis of proximal phalanx morphology from extant caribou populations. PLoS ONE 18(8): e0285487. https://doi.org/10.1371/journal.pone.0285487

**Editor:** Radu Iovita, New York University, UNITED STATES

**Received:** December 15, 2022

**Accepted:** April 25, 2023

**Published:** August 9, 2023

**Copyright:** © 2023 Galán López et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All the 3D models used in the present study as well as the information related to them are available on NAKALA (https://doi.org/10.34847/nkl.fa00p31c).

**Funding:** 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 DeerPal (French National Research Agency, ANR-18-CE03-0007). The funders had no role in study design,

## Abstract

Reindeer mobility patterns vary widely in modern ecosystems, notably between more open or more wooded environments. This renders the reconstruction of past reindeer mobility patterns challenging, while being at the same time key if archaeologists want to better understand the role that reindeer herds played in the subsistence and territorial organisation of Prehistoric hunter-gatherer societies. Studying the morphology associated with different habitats and mobility patterns can be a useful method for understanding past reindeer behaviour. To access paleoecological information, the relationship between locomotor anatomy and substrate type can be explored in modern animals and transposed to the past, as previous studies have proven that an animal's environment affects bone morphology. In this study, 3D Geometric Morphometrics are used to explore the impact of extant reindeer habitat type and mobility pattern on phalanx morphology. Results obtained reflect on the potential archaeological application of such an approach for paleoecological reconstructions. Size and shape of phalanx vary significantly, yet complex to interpret in light of interplaying factors such as subspecies, sexual dimorphism and the influence of migration costs, snow cover and substrate type. If direct application to the archaeological record remains preliminary, this first study permits to highlight promising avenues for future research.

## Introduction

During the Late Pleistocene, the reindeer (*Rangifer tarandus*) was one of the key prey species for many Palaeolithic hunter-gatherer groups [e.g. 1–3]. In Western Europe, several archaeological sequences documented hunting episodes focused on reindeer herds and whose studies
data collection and analysis, decision to publish, or
preparation of the manuscript.

**Competing interests:** The authors have declared
that no competing interests exist.

have highlighted the central importance of this species in the subsistence of past human groups
since the Mousterian to the Magdalenian [e.g. 4–22]. Several authors have suggested that
describing the migration patterns of reindeer is crucial for a better understanding of hunter-
gatherer groups depending on them, from prehistory to present day in the Arctic [Murray
[Unpublished], 23]. Thus, many studies have tried to reconstruct Palaeolithic reindeer migra-
tion patterns [14, 20, 24–41]. Research hypotheses varied from following-herd models, seden-
tary behaviour to migrations carried out on a smaller scale. All, however, pointed out the
difficulty of identifying reindeer migration behaviour with the scarcity of investigation tools at
hand. This is even more difficult to grasp considering the large ecological plasticity of modern
reindeer populations [42, 43].

Extant *Rangifer* (known as caribou in North America, and called reindeer in Eurasia) migrate seasonally to avoid excessive predation on their calving grounds, pests as well as the depletion of local food supplies [44]. Nevertheless, reindeer herds follow different mobility strategies depending on climate conditions, topography and habitat type. In open habitats such as tundra and steppe, mid and long-distance migrations occur [45–48], while animals living in close habitats (e.g. woodlands) tend to move shorter distances [49, 43]. Due to the fact that Pleistocene environments do not necessarily have modern analogues and that reindeer has long been extirpated from Western Europe [7], ethological data can hardly be used to predict past migration patterns, and further palaeoecological data is thus needed. Furthermore, variability in reindeer mobility throughout the Palaeolithic would be expected if one considers the extensive climatic and environmental changes documented for terrestrial ecosystems during the Upper Pleistocene [e.g. 50–52]. For example, with the climate warming at the end of the Late Glacial that lead to reindeer extirpation from southwestern Europe [3, 53], reindeer groups became more fragmented geographically, and different behavioural responses might be expected (e.g. more or less sedentary groups, migration along the Pyrenees, planarly vs. altitudinal migration). Yet, the migration behaviours of reindeer herds during this period are largely unknown, impeding the interpretation of Late Glacial socioeconomic decisions and mobility strategies.

In this study, we seek to develop a new palaeoecological tool for the reconstruction of rein-
deer mobility, taking their bone remains into account in regards to the archaeological record.
Since locomotor adaptations are closely related to habitat preferences and mobility patterns
[39–41, 54–58], behaviour can be inferred from the study of animal morphology. Ecomorpholo-
gical research has proposed a relationship between the degree of mobility (or habitat prefer-
ences) and the functional morphology of limb bones for a wide range of taxa [54, 56, 58–67],
and thus made it possible to successfully reconstruct their habitat [e.g. 54, 56, 58, 60, 61, 63–66].

In this regard, caribou's concave hooves are an interesting adaptation to the harsh and treacherous conditions of northern environments. Their feet are the interface between the environmental substrate and their bodies, and as migratory animals, they rely on their hooves to carry them long distances [68]. In deep snow, caribou walk on their toes, which act like snowshoes. In the fall, their feet begin to harden, and they develop particularly sharp edges that can easily break through ice layers to find food. Likewise, their feet can be used as paddles for swimming or as ice picks when navigating steep, rocky, and icy mountain slopes. The large and flexible feet of caribou also make them extremely useful on uneven, spongy ground when it is not snowing [69]. Furthermore, phalanges are rather dense and strong and are often preserved in archaeological sites.

Due to its use both in livestock practices and domestication, reindeer has been a subsistence
resource for many past and contemporary populations. Thus, several studies on the Fennos-
canian subspecies (*Rangifer tarandus tarandus* and *R. t. fennicus*) have been carried out in
relation to its postcranial morphology (e.g. cross-sectional morphology, enthesal changes)
and identification in the archaeological record (particularly to distinguish the different stages
of its domestication). They have provided an important understanding of this subject, which is
key in reconstructing the mobility and lifestyles of northern peoples [68, 70–77]. In parallel,
the increased development of geometric morphometrics (GMM) provides a powerful set of
techniques for distinguishing subtle morphological differences among members of the same
species. It was successfully applied to the study of ecomorphological patterns in several
instances [67, 75–80].

Building from these previous works, we here extend the application of 3D GMM to phalan-
ges and to a larger sample of wild sedentary and (latitudinal or altitudinal) migratory popula-
tions, in order to explore the potential link between proximal phalanges morphology and
reindeer mobility and habitat. The choice of the proximal phalanx is due to the fact that it is a
often element in the archaeological record of the potential sites to be studied.

Thus, the main purpose of this work is to determine if there are size and shape differences
in proximal phalanx morphology between subspecies (*Rangifer tarandus caribou*, *Rangifer tar-
*andus granti*, *Rangifer tarandus groenlandicus*, *Rangifer tarandus pearyi*) that have different
mobility patterns (altitudinal, longitudinal, sedentary) and live in different habitats (mountain,
boreal forest and tundra). Moreover, since *Rangifer tarandus caribou* subspecies lives in moun-
tainous areas, boreal forests, tundra, and undertakes altitudinal and longitudinal migrations or
has a sedentary behaviour [43, 81], this subspecies (see Sampling section for further details)
will also be studied separately.

Although morphological studies of reindeer often consider sex a key variable [75, 76, 82] or
even other ungulate osteometric studies [83, 84], the presence of sexual dimorphism will be
just visually explored in size (distribution of the log-centroid sizes) in the subspecies *Rangifer
tarandus caribou* for which there is a limited but available sample of males and females.

Finally, as previous studies have shown that allometry is significant [77], this study will
explore the effect of size on the shape of proximal phalanges of both forelimbs and hindlimbs.

**Material and methods**

Sampling extant caribou populations

*Rangifer tarandus* has a wide circumpolar distribution. For this study, we acquired proximal phalanges of individuals belonging to the four subspecies of North American caribou (see Fig 1) (mainly from Canada). These subspecies occupy different habitat settings (boreal, montane and arctic environments), sometimes overlapping, and adopt different mobility strategies. Across this large area, caribou exhibit a wide variation in ecology, genetics, behaviour and morphology. Biologists generally distinguish two main mobility patterns for caribou: migratory and sedentary [42]. However, in mountainous settings some caribou populations perform these migratory movements altitudinally [85], calling it altitudinal or elevational migration.

The subspecies selected for this study are the following:

1) *R. t. caribou* (which includes Eastern caribou/Migratory woodland caribou, mountain caribou and boreal caribou/woodland caribou) inhabits different environments and performs short, medium and long-distance movements. In Quebec, Woodland caribou living south of 55° N are forest-dwellers and seasonally move relatively short distances compared to their tundra counterparts [86, 87], between 50 and 150 kilometres depending on the herd [49, 48, 88, 89]. Eastern or Migratory Woodland caribou group two big herds, George River and Leaf River, that, despite belonging to Woodland caribou, behave like Barren-Ground or tundra caribou, performing migration on long distances ranging between 1120–1770 kilometres on average per year [90], from open tundra settings to boreal forest habitats [91]. Mountain caribou inhabits mountain and boreal settings. Some herds perform altitudinal migrations, with overall
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Location of the *Rangifer tarandus* subspecies in North America. Origin of the specimens used in the present study (notice that this distribution area partially coincides with the current *Rangifer tarandus* distribution in North America, which is wider). Map base by Qyd, Wikimedia Commons, CC Public domain, https://commons.wikimedia.org/wiki/File:Canada-provinces_layout.png.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g001

distances between 140–240 kilometres only in spring migration, depending upon the location of their calving areas, while others remain relatively sedentary, travelling shorter distances under 100 kilometres [81, 92].

2) The Peary caribou (*R. t. pearyi*) inhabits the treeless arctic tundra [44, 93–95] in the Arc-tic Archipelago of Canada (except Baffin Island). These insular caribou notably differ from other tundra populations, mostly by their size (smallest of all *Rangifer tarandus* subspecies). Each season, they perform migrations between islands, by crossing the sea ice and/or swimming [43, 45, 96–98] covering migration distances of up to 500 kilometres [45, 93]. As an example, it has been estimated that seasonal migrations within the Prince of Wales-Somerset-Boothia complex can cover between 300 to 500 kilometres, if not more [99].

3 & 4) The Grant's caribou (*R. t. granti*) and the barren-ground caribou (*R. t. groenlandicus*), referred as Barren-Ground caribou, are alternatively considered to be one or two
subspecies [69]. The first one is mainly found in Alaska and its adjacent Canadian territories, the second one in Canada only (Nunavut, Northwest Territories). These tundra caribou include some of the largest herds in North America, numbering hundreds of thousands of individuals that make very long-distance migrations, ranging between 800 and 5055 kilometres (annual migration) depending on the herds and the year [81, 100–103]. Their habitats are mainly prostrate dwarf shrub tundra and upright shrub tundra [94].

For this research we compiled a sample of 34 fore- and 44 hindlimb proximal (first) phalan- ges belonging to the four subspecies of caribou described above (Tables 1 and 2). All of them were adults and specimens with bone pathologies were excluded from the study. 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, Bar- ren-Ground and Grant's 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.

To facilitate data analysis, each individual was given a unique number (from DPpha 1 to DPpha X, with DP standing for “DeerPal”, the acronym of the project). A detailed table listing

<table><caption>Table 1. Sample used in the present study.</caption><thead><tr><th></th><th>Altitudinal</th><th>Planarly</th><th>Sedentary</th><th>Total</th></tr></thead><tbody><tr><td colspan="5"><strong>(a)</strong></td></tr><tr><td colspan="5"><strong>Forelimb Phalanges</strong></td></tr><tr><td><strong>Rangifer tarandus caribou</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Mountain</td><td>6</td><td></td><td>5</td><td>11</td></tr><tr><td>Tundra</td><td></td><td>6</td><td></td><td>6</td></tr><tr><td>Boreal Forest</td><td></td><td></td><td>7</td><td>7</td></tr><tr><td><strong>Rangifer tarandus ganti</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Tundra</td><td></td><td>4</td><td></td><td>4</td></tr><tr><td><strong>Rangifer tarandus groenlandicus</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Tundra</td><td></td><td>2</td><td></td><td>2</td></tr><tr><td><strong>Rangifer tarandus pearyi</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Tundra</td><td></td><td>4</td><td></td><td>4</td></tr><tr><td><strong>Total</strong></td><td>6</td><td>16</td><td>12</td><td>34</td></tr><tr><td colspan="5"><strong>(b)</strong></td></tr><tr><td colspan="5"><strong>Hindlimb Phalanges</strong></td></tr><tr><td><strong>Rangifer tarandus caribou</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Mountain</td><td>12</td><td></td><td>4</td><td>16</td></tr><tr><td>Tundra</td><td></td><td>10</td><td></td><td>10</td></tr><tr><td>Boreal Forest</td><td></td><td></td><td>10</td><td>10</td></tr><tr><td><strong>Rangifer tarandus ganti</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Tundra</td><td></td><td>2</td><td></td><td>2</td></tr><tr><td><strong>Rangifer tarandus groenlandicus</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Tundra</td><td></td><td>2</td><td></td><td>2</td></tr><tr><td><strong>Rangifer tarandus pearyi</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Tundra</td><td></td><td>4</td><td></td><td>4</td></tr><tr><td><strong>Total</strong></td><td>12</td><td>18</td><td>14</td><td>44</td></tr></tbody></table>

List of the fore-(a) and hindlimb (b) phalanges according to subspecies (*R.t. caribou*, *R.t. granti*, *R.t. groenlandicus* and *R.t. pearyi*), habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal).

https://doi.org/10.1371/journal.pone.0285487.t001
<table><caption>Table 2. Detailed sample.</caption><thead><tr><th>Forelimb Phalanx</th><th>Male</th><th>Female</th><th>Unknown</th><th>Total</th></tr></thead><tbody><tr><td><em>Rangifer tarandus caribou</em></td><td>4</td><td>3</td><td>17</td><td><strong>24</strong></td></tr><tr><td><em>Rangifer tarandus granti</em></td><td>4</td><td>0</td><td>0</td><td><strong>4</strong></td></tr><tr><td><em>Rangifer tarandus groenlandicus</em></td><td>0</td><td>0</td><td>2</td><td><strong>2</strong></td></tr><tr><td><em>Rangifer tarandus pearyi</em></td><td>3</td><td>0</td><td>1</td><td><strong>4</strong></td></tr><tr><td><strong>Total</strong></td><td><strong>11</strong></td><td><strong>3</strong></td><td><strong>20</strong></td><td><strong>34</strong></td></tr></tbody><thead><tr><th>Hindlimb Phalanx</th><th>Male</th><th>Female</th><th>Unknown</th><th>Total</th></tr></thead><tbody><tr><td><em>Rangifer tarandus caribou</em></td><td>3</td><td>7</td><td>26</td><td><strong>36</strong></td></tr><tr><td><em>Rangifer tarandus granti</em></td><td>2</td><td>0</td><td>0</td><td><strong>2</strong></td></tr><tr><td><em>Rangifer tarandus groenlandicus</em></td><td>2</td><td>0</td><td>0</td><td><strong>2</strong></td></tr><tr><td><em>Rangifer tarandus pearyi</em></td><td>4</td><td>0</td><td>0</td><td><strong>4</strong></td></tr><tr><td><strong>Total</strong></td><td><strong>11</strong></td><td><strong>7</strong></td><td><strong>26</strong></td><td><strong>44</strong></td></tr></tbody></table>

Detailed of the fore- and hindlimb phalanges according to subspecies (*R.t. caribou*, *R.t. granti*, *R.t. groenlandicus* and *R.t. pearyi*) and sex (male, female, unknown).

https://doi.org/10.1371/journal.pone.0285487.t002

the habitat and mobility pattern of each phalanx, the number of our unique ID, as well as population information is provided in the S1 Table.

The fact that sedentary populations are classified as “endangered” by Canadian Government makes it more difficult to obtain samples from these herds (since there are very few specimens in museums, the only way to obtain them was through samples sent by biologists when a caribou was killed in an accident or found dead) resulting in uneven samples. Unfortunately, sex was unknown for almost all individuals (n = 14 forelimb phalanges; n = 18 hindlimb phalanges).

Groups were considered according to Subspecies (*Rangifer tarandus caribou*, *Rangifer tarandus granti*, *Rangifer tarandus groenlandicus*, *Rangifer tarandus pearyi*), Mobility (Altitudinal for populations migrating vertically, Planarly for populations migrating mostly on a latitudinal or longitudinal plane, and Sedentary) and Habitat (Mountain, Boreal Forest and Tundra) (see Table 1). In order to analyse *R.t. caribou* subspecies, specimens were grouped according to habitat and mobility pattern (see Table 3).

## Methods

### Data acquisition and 3D geometric morphometrics

Phalanges were individually scanned using a NextEngine 3D laser Scanner (at the Anthropology Department of the Université de Montréal, Canada), in order to obtain a highly accurate

<table><caption>Table 3. Acronyms used in the present study.</caption><thead><tr><th>Name</th><th>Acronym</th></tr></thead><tbody><tr><td>Rtcaribou_BorealForest_Sedentary</td><td>Rtcaribou_BForest_Sed</td></tr><tr><td>Rtcaribou_Mountain_Altitudinal</td><td>Rtcaribou_Mount_Alt</td></tr><tr><td>Rtcaribou_Mountain_Sedentary</td><td>Rtcaribou_Mount_Sed</td></tr><tr><td>Rtcaribou_Tundra_Planarly</td><td>Rtcaribou_Tundra_Plan</td></tr><tr><td>Rtgranti_Tundra_Planarly</td><td>Rtgranti_Tundra_Plan</td></tr><tr><td>Rtgroenlandicus_Tundra_Planarly</td><td>Rtgroenlandicus_Tundra_Plan</td></tr><tr><td>Rtpeary_Tundra_Planarly</td><td>Rtpeary_Tundra_Plan</td></tr></tbody></table>

List of acronyms used in the present study according to subspecies (*R.t. caribou*, *R.t. granti*, *R.t. groenlandicus* and *R.t. pearyi*), habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal).

https://doi.org/10.1371/journal.pone.0285487.t003
3D surface model of each bone. All the 3D models are openly available on NAKALA (https://doi.org/10.34847/nkl.fa00p31c).

Geometric morphometrics (GMM) is a quantitative approach which allows the comparison of bone shapes and the visualization of significant morphological differences between groups of specimens while retaining the element of shape information related to size [104, 105].

Because of the uncertainty in terms of the anatomical position of the curve and surface semilandmarks, sliding step was carried out in curves along their tangent vectors and surface points within their tangent planes in order to minimize the bending energy of the landmark configurations and the template during the Thin-Plate Spline (TPSE) interpolation [106, 107]. After that, the resulting landmark configurations were rotated, translated and scaled using the Generalized Procrustes Analysis (GPA) [108, 109]. This analysis eliminates the effect of position and orientation of the configurations, and limits the effect of size [110–112], yielding Procrustes shape coordinates that were subsequently analysed.

Considering that proximal phalanges are often found fragmented in Palaeolithic sites, we focused our work on proximal epiphyses, in order to allow future applications to a larger part of the archaeological record. The round shape of articular surfaces of these phalanges make it difficult to define true anatomical landmarks, and thus semilandmarks were preferred to define the complex morphology of the joint [106]. We placed landmarks and semilandmarks on areas not affected by enthesal changes and paleopathologies [70, 72]. Two landmark configuration templates were made with a total of 4 landmarks (LM), 43 curve semilandmarks (CSLM) and 25 surface semilandmarks (SSLM) for forelimb phalanges and 4 LM, 45 CSLM and 25 SSLM for hindlimb phalanges, that were digitized to get the form of the proximal epiphyses (Fig 2, Tables 4 and 5). Digitalization was carried out using ViewBox 4.1.0.12 (dHAL software, Kifissia, Greece).

## Statistical analyses: Size and shape differences

First, the four subspecies were analysed altogether. Then, since *R. t. caribou* subspecies constituted a greater percentage (70.5%-forelimb phalanges and 81.8%-hindlimb phalanges) of the sample, the same analyses were only conducted on individuals from this subspecies.

Thus, based on log-transformed centroid sizes for forelimb and hindlimb phalanges, the size differences were assessed using Kruskal-Wallis tests (threshold set at p< 0.05) for the specimens pooled by “Subspecies”, “Habitat”, “Mobility”, and the interactions “Subspecies

<figure><img src="image_5.png" /><figcaption>**Fig 2. Landmarks and semilandmarks location.** Location of anatomical landmarks (red), curve semilandmarks (blue) and surface semilandmarks (green) placed on proximal epiphyses for fore- (left) and hindlimb (right) first phalanges. The definitions of landmarks are given in Tables 4 and 5.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g002
<table><caption>Table 4. List of anatomical landmarks and semilandmarks (forelimb phalanx).</caption><thead><tr><th>Landmark</th><th>Description</th></tr></thead><tbody><tr><td>1</td><td>Abaxial facet of the proximal phalanx, on the edge of the intersection of the curve, just above the small medial articular facet. Insertion tuber edge.</td></tr><tr><td>2</td><td>Most antero-medial point of the lateral glenoid cavity</td></tr><tr><td>3</td><td>Most postero-medial point of the lateral glenoid cavity</td></tr><tr><td>4</td><td>Most postero-distal point of the sagittal groove</td></tr><tr><th colspan="2">Curves (CSLM)</th></tr><tr><td>5:16</td><td>Crest between LM 1 and LM 2.</td></tr><tr><td>17:30</td><td>Crest between LM 2 and LM 3.</td></tr><tr><td>31:38</td><td>Crest between LM 3 and LM 2.</td></tr><tr><td>39:47</td><td>Crest between LM 39 and CSLM 14.</td></tr><tr><th colspan="2">Surfaces (SSLM)</th></tr><tr><td>48:57</td><td>Lateral articular surface patch.</td></tr><tr><td>58:72</td><td>Medial articular surface patch.</td></tr></tbody></table>

Definition of anatomical landmarks and semilandmarks for forelimb proximal phalanx shown in Fig 2.

https://doi.org/10.1371/journal.pone.0285487.t004

+Habitat", "Subspecies+Mobility" and "Mobility+Habitat". Then, multiple Wilcoxon-rank tests were used to compare these groups according to these categories. P values were corrected using the "Benjamini-Hochberg" method to decrease the false discovery rate [113].

A Principal Component Analysis (PCA) based on Procrustes Coordinates was conducted to observe shape variation in the morphospace. Using a Thin-Plate Spline (TPS) interpolation function [108] we generated a 3D digital mesh for each articular surface respectively that was wrapped toward the Procrustes grand mean to understand variations along the principal axes. From the surface of the Procrustes mean configurations [114], shapes at the extremes of the axes of the principal components were visualized and magnified by a scale factor of 0.1 (in order to improve the visualization of the morphological differences), using the function "tps3d" in the "Morpho" package [115].

<table><caption>Table 5. List of anatomical landmarks and semilandmarks (hindlimb phalanx).</caption><thead><tr><th>Landmark</th><th>Description</th></tr></thead><tbody><tr><td>1</td><td>Abaxial facet of the proximal phalanx, on the edge of the intersection of the curve, just above the small medial articular facet. Insertion tuber edge.</td></tr><tr><td>2</td><td>Most antero-medial point of the lateral glenoid cavity</td></tr><tr><td>3</td><td>Most postero-medial point of the lateral glenoid cavity</td></tr><tr><td>4</td><td>Most postero-distal point of the sagittal groove</td></tr><tr><th colspan="2">Curves (CSLM)</th></tr><tr><td>5:18</td><td>Crest between LM 1 and LM 2.</td></tr><tr><td>19:32</td><td>Crest between LM 2 and LM 3.</td></tr><tr><td>33:40</td><td>Crest between LM 3 and LM 2.</td></tr><tr><td>41:49</td><td>Crest between LM 4 and CSLM 16.</td></tr><tr><th colspan="2">Surfaces (SSLM)</th></tr><tr><td>50:58</td><td>Lateral articular surface patch.</td></tr><tr><td>59:74</td><td>Medial articular surface patch.</td></tr></tbody></table>

Definition of anatomical landmarks and semilandmarks for hindlimb proximal phalanx shown in Fig 2.

https://doi.org/10.1371/journal.pone.0285487.t005
A multivariate analysis of variance (MANOVA) on PC scores was used to estimate shape differences between these groups. Additionally, in order to test the assignment accuracy for each category, a Canonical Variate Analysis (CVA) was computed on PC scores (data dimensionality reduction) by calculating the cross-validated correct classification and the Kappa metric (see S1 File). We avoided overfitting [116] by never exceeding 95% of the variance cumulative and by never using more PCs than the number of specimens from the smallest group or category.

As sex variable was known in some specimens from *R. t. caribou* subspecies, sexual size dimorphism was only explored in this subspecies through the Kruskal-Wallis test and if significant, then multiple Wilcoxon-rank test. We also examined shape sexual dimorphism using Principal Component Analysis, Thin-Plate-Spline, and a MANOVA test, as described previously.

Multivariate regressions of shape variables on log-transformed centroid sizes were used to analyse allometry.

All these analyses were performed on RStudio v.4.1.1 [117] using “geomorph” [118] and “Morpho” packages [115].

# Results

## Size variation considering the four subspecies altogether

**Forelimb proximal phalanx.** Kruskal-Wallis tests show significant size differences in almost all of the variables: subspecies (p<0.002), habitat (p<0.05), mobility (p<0.05), and most of the interactions between these variables (Table 6). Mobility+Habitat interaction did not yield significant differences. When subspecies category is included in the interaction for both mobility and habitat, there are always significant differences (p<0.05). However, it must be noted that, for all subspecies expect *R. t. caribou*, only one type of habitat (tundra) and mobility (planarly) are represented in our sample, a composition that affects the statistical

<table><caption>Table 6. Kruskal-Wallis results.</caption><thead><tr><th></th><th>Df</th><th>Chi-square</th><th>P value</th></tr></thead><tbody><tr><td colspan="4"><strong>Forelimb Proximal Phalanx</strong></td></tr><tr><td>Subspecies</td><td>3</td><td>15.132</td><td><strong>0.0017</strong></td></tr><tr><td>Habitat</td><td>2</td><td>6.1572</td><td><strong>0.046</strong></td></tr><tr><td>Mobility</td><td>2</td><td>7.4951</td><td><strong>0.0235</strong></td></tr><tr><td>Subspecies*Habitat</td><td>5</td><td>15.361</td><td><strong>0.0089</strong></td></tr><tr><td>Subspecies*Mobility</td><td>5</td><td>16.699</td><td><strong>0.0051</strong></td></tr><tr><td>Mobility*Habitat</td><td>3</td><td>7.623</td><td>0.054</td></tr><tr><td colspan="4"><strong>Hindlimb Proximal Phalanx</strong></td></tr><tr><td>Subspecies</td><td>3</td><td>14.235</td><td><strong>0.0026</strong></td></tr><tr><td>Habitat</td><td>2</td><td>7.0519</td><td><strong>0.0294</strong></td></tr><tr><td>Mobility</td><td>2</td><td>8.9936</td><td><strong>0.0111</strong></td></tr><tr><td>Subspecies*Habitat</td><td>5</td><td>14.992</td><td><strong>0.0104</strong></td></tr><tr><td>Subspecies*Mobility</td><td>5</td><td>16.933</td><td><strong>0.0046</strong></td></tr><tr><td>Mobility*Habitat</td><td>3</td><td>9.67</td><td><strong>0.0215</strong></td></tr></tbody></table>

Results of the Kruskal-Wallis tests on log-transformed centroid size for fore and hindlimb phalanges according to subspecies (*R. t. caribou*, *R. t. granti*, *R. t. groenlandicus* and *R. t. pearyi*), habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal) and their interactions. A significant contribution was considered for *P* value <0.05 **(in bold)**.

https://doi.org/10.1371/journal.pone.0285487.t006
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Boxplots showing the size variation. Log-transformed centroid size boxplot distribution on forelimb phalanges according to subspecies (*R. t. caribou*, *R. t. granti*, *R. t. groenlandicus* and *R. t. pearyi*), habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g003

results obtained by Kruskall-Wallis tests. Fig 3 highlights the fact that size differences are mainly due to the subspecies variable (i.e. *R. t. groenlandicus* and *R. t. pearyi* are smaller).

Pairwise comparisons did not show size differences when they were displayed by subspecies +mobility+habitat. However, although it is not statistically significant, larger sizes correspond to male specimens which performs altitudinal migrations (*R. t. caribou*, mountain setting) (Table 7) as it can be observed in boxplots (Fig 3).

**Hindlimb proximal phalanx.** Kruskal-Wallis' results (see Table 6) on log-transformed centroid sizes yielded significant differences (p<0.05) among all the categories: Subspecies, Habitat, Mobility, Subspecies+Mobility, Subspecies+Habitat and Mobility+Habitat. Yet, again, this statistical test is affected by our sample composition. Pairwise comparisons (Table 8) by subspecies+mobility+habitat showed size differences (p<0.05) between *R. t. pearyi* and three groups of *R. t. caribou* (BForest_Sedentary; Mountain_Altitudinal and Tundra_Planarly), the former being smaller (Fig 4).

## Size variation within *Rangifer tarandus caribou*

*R. t. caribou* is the only subspecies in our sample that include specimens of varying habitats and mobility. Kruskal-Wallis yielded no statistical differences (Table 9) on size for *Rangifer tarandus caribou* subspecies. Therefore, no pairwise comparisons were carried out. A very small size variation is perceptible in the boxplots (see S2 File) when both habitat and mobility are considered, even if specimens that carried out altitudinal migration have slightly bigger phalanges. These small differences could be explained, at least in part, by differences in sex ratios between subpopulations (see Table 2), although this would be difficult to explore without a larger sample size.

## Shape variation results considering the four subspecies altogether

**Forelimb proximal phalanx.** The first 2 PCs explain 46.42% of the morphological variance (Fig 5). Shape variation along PC1 accounted for 28.65% of the variance while PC2, 17.77%. Although there are overlaps according to both Habitat and Mobility type (Sedentary,
<table><caption>Table 7. Pairwise comparisons of log-transformed centroid-sizes.</caption><thead><tr><th>Subspecies</th><th>Ricaribou_BorealForest_Sedentary</th><th>Ricaribou_Mountain_Sedentary</th><th>Ricaribou_Mountain_Altitudinal</th><th>Ricaribou_Tundra_Sedentary</th><th>Ricaribou_Tundra_Planarly</th><th>Riganti_Tundra_Planarly</th><th>Riganti_Tundra_Planarly</th><th>Rigroenlandicus_Tundra_Planarly</th></tr></thead><tbody><tr><td>Ricaribou_Mountain_Altitudinal</td><td>0.530</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ricaribou_Mountain_Sedentary</td><td>0.943</td><td>0.211</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ricaribou_Tundra_Planarly</td><td>0.943</td><td>0.660</td><td>0.801</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Riganti_Tundra_Planarly</td><td>0.801</td><td>0.161</td><td>0.831</td><td>0.943</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Rigroenlandicus_Tundra_Planarly</td><td>0.161</td><td>0.161</td><td>0.171</td><td>0.161</td><td>0.191</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ripeary_Tundra_Planarly</td><td>0.099</td><td>0.099</td><td>0.105</td><td>0.099</td><td>0.123</td><td>-</td><td>0.399</td><td>-</td></tr></tbody></table>

Pairwise comparisons of log-transformed centroid sized for forelimb phalanx between the different groups according to subspecies (R.t. caribou, R.t. granti, R.t. groenlandicus and R.t. pearyi), habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal) using Wilcoxon rank sum test after “Benjamini-Hochberg” correction. A significant contribution was considered for *P* value <0.05 (in bold).

https://doi.org/10.1371/journal.pone.0285487.t007
<table><caption>Table 8. Pairwise comparisons of log-transformed centroid-sizes.</caption><thead><tr><th>Subspecies</th><th>Rticaribou_BForest_Sedentary</th><th>Rticaribou_Mountain_Sedentary</th><th>Rticaribou_Mountain_Altitudinal</th><th>Rticaribou_Mountain_Sedentary</th><th>Rticaribou_Tundra_Sedentary</th><th>Rticaribou_Tundra_Planarly</th><th>Rigroen_Tundra_Planarly</th><th>Rigroen_Tundra_Planarly</th></tr></thead><tbody><tr><td>Rticaribou_Mountain_Altitudinal</td><td>0.460</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Rticaribou_Mountain_Sedentary</td><td>0.460</td><td>0.180</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Rticaribou_Tundra_Planarly</td><td>1</td><td>0.415</td><td>0.826</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Rigroen_Tundra_Planarly</td><td>0.415</td><td>0.125</td><td>0.602</td><td>0.602</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Rigroen_Tundra_Planarly</td><td>0.415</td><td>0.125</td><td>0.602</td><td>0.826</td><td>1</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Rticaribou_Tundra_Planarly</td><td>0.041</td><td>0.041</td><td>0.125</td><td>0.041</td><td>0.245</td><td>0.245</td><td>0.245</td><td>0.245</td></tr></tbody></table>

Pairwise comparisons of log-transformed centroid sized for hindlimb phalanx between the different groups according to subspecies (R.t. caribou, R.t. granti, R.t. groenlandicus and R.t. pearyi), habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal) using Wilcoxon rank sum test after "Benjamini-Hochberg" correction. A significant contribution was considered for *P* value <0.05 (in bold).

https://doi.org/10.1371/journal.pone.0285487.t008
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Boxplots showing the size variation. Log-transformed centroid size boxplot distribution on hindlimb phalanges according to subspecies (*R.t. caribou*, *R.t. granti*, *R.t. groenlandicus* and *R.t. pearyi*), habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g004

Planarly and Altitudinal), specimens whose migration is altitudinal showed their variation mainly along negative values of PC1. No pattern is evident for sedentary individuals, with Rtcaribou_BForest_Sed falling mainly towards positive PC1 values and negative PC2, while Rtcaribou_Mount_Sed are on negative PC1 and positive PC2 values. Planarly migraters are spread across the morphospace: Rtcaribou_Tundra_Plan are distributed along PC1 positive and PC2 negative values, Rtpeary_Tundra_Plan at the center of the morphospace, whereas Rtgranti_Tundra_Plan and Rtgroenlandicus_Tundra_Plan are mainly towards negative values of PC1.

Towards the positive values of PC1, proximal morphology is more massive towards cranial side and mediolaterally narrower (Fig 6). Axial ridge resulted to be straighter and abaxial artic- ular surface broader. However, towards negative values, proximal joint appeared to be more prominent mediolaterally than craniocaudally. Middle groove was deeper towards the negative values.

<table><caption>**Table 9.** Kruskal-Wallis results on log-transformed centroid size.</caption><thead><tr><th></th><th>Df</th><th>Chi-square</th><th>P value</th></tr></thead><tbody><tr><td colspan="4"><strong>Forelimb Proximal Phalanx</strong></td></tr><tr><td>Habitat</td><td>2</td><td>0.19013</td><td>0.9093</td></tr><tr><td>Mobility</td><td>2</td><td>2.04</td><td>0.3606</td></tr><tr><td>Mobility *Habitat</td><td>3</td><td>2.2869</td><td>0.515</td></tr><tr><td colspan="4"><strong>Hindlimb Proximal Phalanx</strong></td></tr><tr><td>Habitat</td><td>2</td><td>0.62883</td><td>0.7302</td></tr><tr><td>Mobility</td><td>2</td><td>2.9569</td><td>0.228</td></tr><tr><td>Mobility *Habitat</td><td>3</td><td>3.8991</td><td>0.2726</td></tr></tbody></table>

Results of the Kruskal-Wallis tests on log-transformed centroid size for fore and hindlimb phalanges for *Rangifer tarandus caribou* subspecies according to habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal) and their interaction. A significant contribution was considered for *P* value <0.05 (in bold).

https://doi.org/10.1371/journal.pone.0285487.t009
<figure><img src="image_8.png" /><figcaption>**Fig 5. Principal component analysis.** Scatter plots of the two first axes (PC1 and PC2) of principal component analyses performed on the shape data (forelimb phalanx) according to subspecies (*R.t. caribou*, *R.t. granti*, *R.t. groenlandicus* and *R.t. pearyi*), habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal). The result of the same PCA is expressed in two scatterplots (for a better visualization). On the left, the highlighted individuals correspond to *R.t. caribou* subspecies, and on the right, those of the remaining subspecies.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g005

Thus, the abaxial view is more massive towards the positive values of PC1 (mainly individu-
als from tundra habitat and planarly migration) and narrower towards the negative values of
the same axis (mostly Rtcaribou_Mountain_Altitudinal specimens). In the caudal view, the
morphology tends to be more inwardly curved towards positive values of PC1 and straighter
towards negative ones. In the same way, in the axial view, a greater protuberance is observed
in the lower part of the articular surfaces towards the positive values of the PC1 axis and a flat-
ter configuration towards the negative ones.

Regarding PC2, the main differences are observed in the proximal view of the articular sur-
face. In the axial view, a narrower configuration is observed towards the positive values of this

<figure><img src="image_9.png" /><figcaption>**Fig 6. Visualization of shape variation.** Shape variation visualized by deforming the mean shape along negative and positive PC1 and PC2 values (magnified by a scale factor of 0.1). Forelimb phalanges.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g006
<table><caption>Table 10. MANOVA tests' results.</caption><thead><tr><th></th><th>Df</th><th>Pillai</th><th>approx F</th><th>num Df</th><th>den Df</th><th>Pr(&gt;F)</th></tr></thead><tbody><tr><td colspan="7"><strong>Forelimb Proximal Phalanx</strong></td></tr><tr><td>Subspecies</td><td>3</td><td>1.3329</td><td>1.3993</td><td>36</td><td>63</td><td>0.1204</td></tr><tr><td>Habitat</td><td>2</td><td>1.1725</td><td>2.4794</td><td>24</td><td>42</td><td><strong>0.0048</strong></td></tr><tr><td>Mobility</td><td>2</td><td>0.9322</td><td>1.5279</td><td>24</td><td>42</td><td>0.1124</td></tr><tr><td>ssp+habitat+mobility</td><td>6</td><td>3.1578</td><td>1.9443</td><td>72</td><td>126</td><td><strong>0.0005</strong></td></tr><tr><td>Subspecies*Habitat</td><td>5</td><td>2.7339</td><td>2.1113</td><td>60</td><td>105</td><td><strong>0.0004</strong></td></tr><tr><td>Subspecies*Mobility</td><td>5</td><td>2.4937</td><td>1.7412</td><td>60</td><td>105</td><td><strong>0.0065</strong></td></tr><tr><td>Mobility*Habitat</td><td>3</td><td>1.5963</td><td>1.9901</td><td>36</td><td>63</td><td><strong>0.0083</strong></td></tr><tr><td colspan="7"><strong>Hindlimb Proximal Phalanx</strong></td></tr><tr><td>Subspecies</td><td>3</td><td>1.5507</td><td>2.2165</td><td>42</td><td>87</td><td><strong>0.0009</strong></td></tr><tr><td>Habitat</td><td>2</td><td>1.2061</td><td>3.1468</td><td>28</td><td>58</td><td><strong>0.0001</strong></td></tr><tr><td>Mobility</td><td>2</td><td>1.1497</td><td>2.8006</td><td>28</td><td>58</td><td><strong>0.0004</strong></td></tr><tr><td>ssp+habitat+mobility</td><td>6</td><td>3.0831</td><td>2.1894</td><td>84</td><td>174</td><td><strong>7.503e-06</strong></td></tr><tr><td>Subspecies*Habitat</td><td>5</td><td>2.7631</td><td>2.5588</td><td>70</td><td>145</td><td><strong>9.987e-07</strong></td></tr><tr><td>Subspecies*Mobility</td><td>5</td><td>2.7067</td><td>2.4449</td><td>70</td><td>145</td><td><strong>3.094e-06</strong></td></tr><tr><td>Habitat* Mobility</td><td>3</td><td>1.526</td><td>2.1446</td><td>42</td><td>87</td><td><strong>0.0014</strong></td></tr></tbody></table>

MANOVA tests' results for forelimb and hindlimb proximal phalanges according to subspecies (R.t. caribou, R.t. granti, R.t. groenlandicus and R.t. pearyi), habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal) and their interactions. A significant contribution was considered for *P* value <0.05 (in **bold**).

https://doi.org/10.1371/journal.pone.0285487.t010

axis. The theoretical shape at the PC2 minimum showed a proximal joint that is larger cranio-caudally and narrower on an axial-abaxial axe and the opposite it is shown towards its maximum shape and occurs more often on mountain specimens.

MANOVA results (Table 10) produced significant differences in shape according to Habitat (p = 0.0048); Ssp+Habitat+Mobility (p = 0.0005); Subspecies*Habitat (p = 0.0004); Subspecies*Mobility (p = 0.0065) and Mobility*Habitat (p = 0.0083). Contrary, there is no shape difference according to Subspecies and Mobility. When analyses are displayed by pooling according to Ssp+Habitat+Mobility significant differences were found between: on one hand, Rtcaribou_BForest_Sed, and, on the other hand, Rtcaribou_Tundra_Plan (p = 0.016) and Rtpeary_Tundra_Plan (p = 0.025) (see Table 11).

CVA results for classification (Table 12) obtained 82.35% overall accuracy according to habitat and a lower percentage in accuracy when it was performed according to mobility pattern (73.52%). These percentages are a bit higher when calculated according to 2 mobility groups (migratory/sedentary), 82.3%. Thus, classification accuracy in this case improves when mobility patterns are grouped into two categories rather than multiple categories. In all cases, Kappa statistic (see S1 File) was higher than 0.6, which indicates a substantial agreement beyond chance.

## Hindlimb phalanges

First two PCs accounts for 48.61% of the morphological variance (Fig 7). Rtcaribou_BForest_Sed are distributed towards positive and negative values of PC1, while Rtcaribou_-Mount_Sed fell towards positive values of PC1 and Rtcaribou_Mount_Alt towards negative values of PC2. R.t. peary as well as R.t. granti and groenlandicus, fall mainly towards positive values of PC2, whereas Rtcaribou_Tundra_Plan were falling towards both positive and negative values of PC1.
<table><caption>Table 11. MANOVA pairwise comparisons results.</caption><thead><tr><th></th><th>Pillai</th><th>approx F</th><th>num DF</th><th>den DF</th><th>Pr(&gt;F)</th></tr></thead><tbody><tr><td>Rtcaribou_BorealForest_Sedentary—Rtcaribou_Mountain_Altitudinal</td><td>0.68108</td><td>2.8475</td><td>12</td><td>16</td><td>0.55760</td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtcaribou_Mountain_Sedentary</td><td>0.75484</td><td>4.1054</td><td>12</td><td>16</td><td>0.10423</td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtcaribou_Tundra_Planarly</td><td>0.81225</td><td>5.7684</td><td>12</td><td>16</td><td><strong>0.01683</strong></td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtgranti_Tundra_Planarly</td><td>0.61398</td><td>2.1207</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtgroenlandicus_Tundra_Planarly</td><td>0.60321</td><td>2.0270</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtpeary_Tundra_Planarly</td><td>0.80121</td><td>5.3740</td><td>12</td><td>16</td><td><strong>0.02508</strong></td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtcaribou_Mountain_Sedentary</td><td>0.48980</td><td>1.2800</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtcaribou_Tundra_Planarly</td><td>0.74740</td><td>3.9450</td><td>12</td><td>16</td><td>0.12702</td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtgranti_Tundra_Planarly</td><td>0.40312</td><td>0.9005</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtgroenlandicus_Tundra_Planarly</td><td>0.40754</td><td>0.9172</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtpeary_Tundra_Planarly</td><td>0.63921</td><td>2.3622</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtcaribou_Mountain_Sedentary—Rtcaribou_Tundra_Planarly</td><td>0.65022</td><td>2.4786</td><td>12</td><td>16</td><td>0.96660</td></tr><tr><td>Rtcaribou_Mountain_Sedentary—Rtgranti_Tundra_Planarly</td><td>0.54285</td><td>1.5833</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtcaribou_Mountain_Sedentary—Rtgroenlandicus_Tundra_Planarly</td><td>0.41165</td><td>0.9329</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtcaribou_Mountain_Sedentary—Rtpeary_Tundra_Planarly</td><td>0.51080</td><td>1.3922</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtcaribou_Tundra_Planarly—Rtgranti_Tundra_Planarly</td><td>0.75086</td><td>4.0185</td><td>12</td><td>16</td><td>0.11596</td></tr><tr><td>Rtcaribou_Tundra_Planarly—Rtgroenlandicus_Tundra_Planarly</td><td>0.52656</td><td>1.4829</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtcaribou_Tundra_Planarly—Rtpeary_Tundra_Planarly</td><td>0.57170</td><td>1.7798</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtgranti_Tundra_Planarly—Rtgroenlandicus_Tundra_Planarly</td><td>0.58240</td><td>1.8595</td><td>12</td><td>16</td><td>1.00000</td></tr><tr><td>Rtgranti_Tundra_Planarly—Rtpeary_Tundra_Planarly</td><td>0.68400</td><td>2.8861</td><td>12</td><td>16</td><td>0.52725</td></tr><tr><td>Rtgroenlandicus_Tundra_Planarly—Rtpeary_Tundra_Planarly</td><td>0.47736</td><td>1.2178</td><td>12</td><td>16</td><td>1.00000</td></tr></tbody></table>

Results of the MANOVA pairwise comparisons for forelimb proximal phalanx according to the different categories. A significant contribution was considered for *P* value <0.05 (in bold).

https://doi.org/10.1371/journal.pone.0285487.t011

The PC1 minimum configuration (Fig 8) seemed to be more stretched craniodorsally and narrower mediolaterally. Contrary, towards the positive values of PC1 or maximum configuration expressed a theoretical morphology that was broader mediolaterally and more elongated dorsally in the lateral articular surface.

<table><caption>Table 12. Canonical Variate Analysis (CVA) results.</caption><thead><tr><th></th><th>PCs</th><th>CCV (%)</th><th>Kappa</th></tr></thead><tbody><tr><td colspan="4"><strong>Forelimb Proximal Phalanx</strong></td></tr><tr><td>Habitat</td><td>6</td><td>82.35</td><td>0.72</td></tr><tr><td>Mobility</td><td>7</td><td>73.52</td><td>0.59</td></tr><tr><td colspan="4"><strong>Hindlimb Proximal Phalanx</strong></td></tr><tr><td>Habitat</td><td>11</td><td>84.09</td><td>0.75</td></tr><tr><td>Mobility</td><td>12</td><td>86.36</td><td>0.79</td></tr></tbody></table>

Canonical Variate Analysis (CVA) indicating the number of PC scores (PCs) used in the analysis, the percentage of correct cross-validated classification according to habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal) for fore and hindlimb phalanges.

https://doi.org/10.1371/journal.pone.0285487.t012
<figure><img src="image_10.png" /><figcaption>**Fig 7. Principal component analysis.** Scatter plots of the two first axes (PC1 and PC2) of principal component analyses performed on the shape data (hindlimb phalanx) according to subspecies (*R.t. caribou*, *R.t. granti*, *R.t. groenlandicus* and *R.t. pearyi*), habitat (mountain, tundra and boreal forest) and mobility type (sedentary, planarly and altitudinal). The result of the same PCA is expressed in two scatterplots (for a better visualization). On the left, the highlighted individuals correspond to *R.t. caribou* subspecies, and on the right, those of the remaining subspecies.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g007

Towards positive values of PC1 (Fig 8), the morphology tends to be more curved towards the abaxial side. On this same side, the configuration is more massive towards positive than negative values, however in the caudal and axial views the morphology is wider and more pronounced towards negative than positive values of the PC1 axis.

In PC2 (Fig 8), the morphological configuration in this axis stands out for being more pronounced towards its positive values. The abaxial view shows a greater protuberance in the lower part of the articular surfaces towards the positive values of the PC2 axis and a flatter configuration towards the negative ones. However, in axial view, it tends to bulge more towards negative values, while being flatter towards positive values of the axis.

<figure><img src="image_11.png" /><figcaption>**Fig 8. Visualization of shape variation.** Shape variation visualized by deforming the mean shape along negative and positive PC1 and PC2 values (magnified by a scale factor of 0.1). Hindlimb phalanges.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g008
<table><caption>Table 13. MANOVA pairwise comparisons results.</caption><thead><tr><th></th><th>Pillai</th><th>approx F</th><th>num DF</th><th>den DF</th><th>Pr(&gt;F)</th></tr></thead><tbody><tr><td>Rtcaribou_BorealForest_Sedentary—Rtcaribou_Mountain_Altitudinal</td><td>0.65869</td><td>3.3084</td><td>14</td><td>24</td><td>0.102651</td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtcaribou_Mountain_Sedentary</td><td>0.51214</td><td>1.7996</td><td>14</td><td>24</td><td>1.000000</td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtcaribou_Tundra_Planarly</td><td>0.64209</td><td>3.0754</td><td>14</td><td>24</td><td>0.158975</td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtgranti_Tundra_Planarly</td><td>0.72532</td><td>4.5268</td><td>14</td><td>24</td><td><strong>0.012574</strong></td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtgroenlandicus_Tundra_Planarly</td><td>0.59728</td><td>2.5425</td><td>14</td><td>24</td><td>0.451379</td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtpeary_Tundra_Planarly</td><td>0.75467</td><td>5.2734</td><td>14</td><td>24</td><td><strong>0.004002</strong></td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtcaribou_Mountain_Sedentary</td><td>0.36959</td><td>1.0050</td><td>14</td><td>24</td><td>1.000000</td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtcaribou_Tundra_Planarly</td><td>0.69112</td><td>3.8358</td><td>14</td><td>24</td><td><strong>0.039837</strong></td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtgranti_Tundra_Planarly</td><td>0.69538</td><td>3.9134</td><td>14</td><td>24</td><td><strong>0.034828</strong></td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtgroenlandicus_Tundra_Planarly</td><td>0.48744</td><td>1.6303</td><td>14</td><td>24</td><td>1.000000</td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtpeary_Tundra_Planarly</td><td>0.62676</td><td>2.8786</td><td>14</td><td>24</td><td>0.232115</td></tr><tr><td>Rtcaribou_Mountain_Sedentary—Rtcaribou_Tundra_Planarly</td><td>0.47110</td><td>1.5269</td><td>14</td><td>24</td><td>1.000000</td></tr><tr><td>Rtcaribou_Mountain_Sedentary—Rtgranti_Tundra_Planarly</td><td>0.59613</td><td>2.5304</td><td>14</td><td>24</td><td>0.462520</td></tr><tr><td>Rtcaribou_Mountain_Sedentary—Rtgroenlandicus_Tundra_Planarly</td><td>0.45291</td><td>1.4192</td><td>14</td><td>24</td><td>1.000000</td></tr><tr><td>Rtcaribou_Mountain_Sedentary—Rtpeary_Tundra_Planarly</td><td>0.45560</td><td>1.4347</td><td>14</td><td>24</td><td>1.000000</td></tr><tr><td>Rtcaribou_Tundra_Planarly—Rtgranti_Tundra_Planarly</td><td>0.74778</td><td>5.0826</td><td>14</td><td>24</td><td><strong>0.005313</strong></td></tr><tr><td>Rtcaribou_Tundra_Planarly—Rtgroenlandicus_Tundra_Planarly</td><td>0.53954</td><td>2.0087</td><td>14</td><td>24</td><td>1.000000</td></tr><tr><td>Rtcaribou_Tundra_Planarly—Rtpeary_Tundra_Planarly</td><td>0.69918</td><td>3.9845</td><td>14</td><td>24</td><td><strong>0.030831</strong></td></tr><tr><td>Rtgranti_Tundra_Planarly—Rtgroenlandicus_Tundra_Planarly</td><td>0.67240</td><td>3.5185</td><td>14</td><td>24</td><td>0.069901</td></tr><tr><td>Rtgranti_Tundra_Planarly—Rtpeary_Tundra_Planarly</td><td>0.44625</td><td>1.3815</td><td>14</td><td>24</td><td>1.000000</td></tr><tr><td>Rtgroenlandicus_Tundra_Planarly—Rtpeary_Tundra_Planarly</td><td>0.54770</td><td>2.0759</td><td>14</td><td>24</td><td>1.000000</td></tr></tbody></table>

Results of the MANOVA pairwise comparisons for hindlimb proximal phalanx according to the different categories. A significant contribution was considered for *P* value <0.05 (in **bold**).

https://doi.org/10.1371/journal.pone.0285487.t013

Regarding PC2 (Fig 8), we find the opposite as in PC1. The theoretical morphology at the minimum appeared to be more massive mediolaterally in the proximal articular surface and with a deeper groove. However, it was narrower abaxially at the PC2 maximum.

MANOVA analyses yielded significant differences (Table 10) in shape across all categories: Subspecies (p = 0.0009); Habitat (p = 0.0001); Mobility (0.0004); Ssp+Habitat+Mobility (p = 7.503e-06); Subspecies*Habitat (p = 9.987e-07); Subspecies*Mobility (p = 3.094e-06); Mobility*Habitat (p = 0.0014). Multiple comparisons (Table 13) according to Ssp+Habitat+ Mobility yielded significant results between Rtcaribou_BForest_Sed on one hand, Rtgranti and Rtpeary_Tundra_Planarly on the other (p = 0.012 and 0.004 respectively); Rtcaribou_-Mount_Alt and Rtcaribou and Rtgranti_Tundra_Plan (p = 0.039 and 0.034) and finally between Rtcaribou_Tundra_Plan and Rtgranti and Rtpeary_Tundra_Plan (p = 0.005 and 0.03).

Regarding CVA cross-validated results (Table 12), accuracy values obtained were over 80%. According to habitat (mountain, tundra, boreal forest), the overall accuracy was 84.09%, and 86.36% for mobility pattern (planarly, altitudinal and sedentary). In all cases Kappa metric was higher than 0.79 (see S1 File), what means a good agreement between predicted and true labels.
## Shape variation amongst *Rangifer tarandus caribou*

In forelimb phalanges, the first two PCs account for 49.76% of the total variance (Fig 9A). Both Altitudinal and Sedendary Rtcaribou_Mountain are distributed in the morphospace towards positive and negative values of PC2, while Rtcaribou_BForest are distributed mainly towards positive values of PC1 and Rtcaribou_Tundra_Planarly along positive values of PC2. Positive PC1 values reflect the characteristics of Rtcaribou_BForest_Sedendary individuals, which means a more prominent axial articular surface cranially compared to the negative PC1 values.

Individuals with positive PC2 values (sedendary specimens) display a joint morphology (Fig 10) that is more stretched abaxially and elongated in the craniocaudal axis; however, individuals with negative PC2 values have a joint that is wider in the abaxial axis (Rtcaribou_Tundra_Planarly and Rtcaribou_Mountain_Altitudinal).

In the case of hindlimb phalanges, the analysis of the PCA reveals a complex configuration. PC1 and PC2 account for 50% of the variance (Fig 9B). Rtcaribou_tundra and Mountain are mainly distributed along PC1 with positive and negative values. In contrast, specimens of Rtcaribou_BForest_Sedendary are mainly distributed along positive values of PC1 and PC2, overlapping with some Rtcaribou_Mountain_Altitudinal individuals.

Thus, specimens toward positive values of PC1 are shown as having an axial articular surface that is more elongated than those towards negative values (Fig 11). Those specimens towards PC2 positive values show a stretched morphology on the abaxial axis and a more elongated axial surface towards the caudal side, whereas those that have negative values have the opposite configuration (wider abaxially and narrower craniocaudally).

In both fore- and hindlimb phalanges (Table 14), there are significant differences (p<0.05) according to habitat, mobility, and the interactions between the two. Nonetheless, pairwise comparisons revealed differences between all habitat categories and between altitudinal and planarly migration types in forelimb phalanx (Tables 15 and 16).

All habitat categories (Table 17) and mobility (Table 18) only for Altitudinal-Planarly showed significant differences for hindlimb phalanges. Performance according to Habitat+-Mobility (Table 19) yielded differences between Rtcaribou_BForest_Sedendary and Rtcaribou_Mountain_Altitudinal (p = 0.01), BForest_Sedendary and Tundra_Planarly (p = 0.02) and Mountain_Altitudinal and Tundra_Planarly (p = 0.03).

## Sexual dimorphism on size (*R.t. caribou*)

Despite the fact that sex is frequently considered an important variable in reindeer morphological studies [75, 77, 83, 88] size and shape sexual dimorphism can hardly be addressed in our

<figure><img src="image_12.png" /><figcaption>**Fig 9.** Principal component analysis. Scatter plots of the two first axes (PC1 and PC2) of principal component analyses performed on the shape data for *Rangifer tarandus caribou* according to habitat (mountain, tundra and boreal forest) and mobility type (sedendary, planarly and altitudinal); (a) Forelimb; (b) hindlimb.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g009
<figure><img src="image_13.png" /><figcaption>**Fig 10. Visualization of shape variation.** Shape variation visualized by deforming the mean shape along negative and positive PC1 and PC2 values (magnified by a scale factor of 0.1). Forelimb phalanges from *R.t. caribou*.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g010

sample, due to the low number of specimens of known sex (Table 2). Only the *R. t. caribou* sample have specimens of both sexes. The lack of enough specimens prevents further analyses due to low statistical power to detect meaningful differences between groups. In contrast, box-plots were used to visualize log-centroid size distribution. Thus, when the variation in log-transformed centroid size is displayed in boxplots (see S2 File) according to male/female it can be observed that males tend to be larger in size than females, although there is an important overlap. Considering these results and our limited sample in terms of specimens of known sex, we discarded the sex variable for further analyses.

## Allometry

For forelimb phalanx, the percentage of variance related to size was low, only a 5.90% and the allometry was significant (p = 0.04). However, for hindlimb phalanx, the percentage of variance related to shape was even lower, at 2.34%, and allometry was not significant (p = 0.38).

Multivariate regression (Fig 12) of shape scores against log-transformed centroid size showed a clear separation of *R. t. caribou* with the rest of the subspecies (except *R. t. granti*) in forelimb phalanx. For hindlimb phalanx (Fig 12), *R. t. pearyi* is completely separated from the rest of the subspecies, which overlap among them.

## Discussion

In the present study, we have examined the morphology of caribou proximal phalanx and its relationship with mobility and habitat type in order to potentially apply into archaeological
<figure><img src="image_14.png" /><figcaption>**Fig 11. Visualization of shape variation.** Shape variation visualized by deforming the mean shape along negative and positive PC1 and PC2 values (magnified by a scale factor of 0.1). Hindlimb phalanges from *R.t. caribou*.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g011

record for palaeoecological reconstructions. Despite the fact that a wider and more complete sample for further analyses is needed, this work represents a first step in this regard.

As we have already mentioned along our study, the main limitations of our sample are the size and the lack of information on the sex of the specimens. 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 samples sent by biologists acquired when a caribou died in a traffic accident or was

<table><caption>**Table 14. MANOVA tests’ results.**</caption><thead><tr><th></th><th>Df</th><th>Pillai</th><th>approx F</th><th>Num Df</th><th>den Df</th><th>Pr(&gt;F)</th></tr></thead><tbody><tr><td colspan="7"><strong>Forelimb Proximal Phalanx</strong></td></tr><tr><td>Habitat</td><td>2</td><td>1.4987</td><td>4.6508</td><td>18</td><td>28</td><td><strong>0.000148</strong></td></tr><tr><td>Mobility</td><td>2</td><td>1.2746</td><td>2.7334</td><td>18</td><td>28</td><td><strong>0.008286</strong></td></tr><tr><td>Habitat* Mobility</td><td>3</td><td>1.8801</td><td>2.6115</td><td>27</td><td>42</td><td><strong>0.002546</strong></td></tr><tr><td colspan="7"><strong>Hindlimb Proximal Phalanx</strong></td></tr><tr><td>Habitat</td><td>2</td><td>1.3199</td><td>3.7197</td><td>24</td><td>46</td><td><strong>6.441e-05</strong></td></tr><tr><td>Mobility</td><td>2</td><td>1.1897</td><td>2.814</td><td>24</td><td>46</td><td><strong>0.001268</strong></td></tr><tr><td>Habitat* Mobility</td><td>3</td><td>1.6198</td><td>2.2494</td><td>36</td><td>69</td><td><strong>0.001934</strong></td></tr></tbody></table>

MANOVA tests’ results for forelimb and hindlimb proximal phalanges for *Rangifer tarandus caribou* according to habitat (mountain, tundra and boreal forest) and mobility pattern (sedentary, planarly and altitudinal) and their interaction. A significant contribution was considered for *P* value <0.05 (in bold).

https://doi.org/10.1371/journal.pone.0285487.t014
<table><caption>Table 15. MANOVA pairwise comparisons results (habitat-forelimb phalanx).</caption><thead><tr><th>Habitat</th><th>Pillai</th><th>approx F</th><th>num DF</th><th>den DF</th><th>Pr(&gt;F)</th></tr></thead><tbody><tr><td>BorealForest—Mountain</td><td>0.72827</td><td>3.8712</td><td>9</td><td>13</td><td><strong>0.041324</strong></td></tr><tr><td>BorealForest—Tundra</td><td>0.73151</td><td>3.9354</td><td>9</td><td>13</td><td><strong>0.038709</strong></td></tr><tr><td>Mountain—Tundra</td><td>0.80401</td><td>5.9257</td><td>9</td><td>13</td><td><strong>0.006563</strong></td></tr></tbody></table>

Results of the MANOVA pairwise comparisons for *Rangifer tarandus caribou* (forelimb phalanx) according to habitat type (mountain, tundra, boreal forest). A significant contribution was considered for *P* value <0.05 **(in bold)**.

https://doi.org/10.1371/journal.pone.0285487.t015

found dead in the forest), resulting in uneven sample sizes. Regarding sex, most of the specimens from museum collections were collected at the end of the 19th or beginning of the 20th century and they did not specify the sex, so we can only know it from those that were collected in situ by biologists nowadays. Because only few individuals within the sample population are of known sex, we cannot directly evaluate the relationship between tested values and sex regarding to habitat or mobility.

Specimens with bone pathologies [71, 72] were excluded from the study to avoid this bias in our sample. Despite the limitations imposed by our sample size, there is a consonance between the statistical analyses and the results in the configuration of the morphospace. Due to this reason, we considered including the subspecies category as a preliminary step while waiting for a larger sample size.

Studying joint surfaces allows for the use of complete or fragmented bones, increasing the sample size that can be used in archaeological records [67]. Furthermore, when compared to previous studies with linear measurements on reindeer phalanges, a 3D GMM study has the advantage of allowing a deeper study in the morphology of these phalanges, and is able to accurately capture shape variation, even with a small sample size (e.g. altitudinal migration was difficult to explore using linear measurements, cf. [41]). For validation of the methods presented here, we recommend similar studies on caribou populations in other regions so that we can assess the preliminary variations we observed.

## Interpreting size variation: A complex interplay of sexual dimorphism and migration costs?

Size variations are mainly imputable to subspecies differences. *R. t. caribou* are statistically significantly larger than *R. t. pearyi*, which was expected, since the latter is amongst the smallest *R. t.* subspecies [119]. Peary Caribou, inhabit the Canadian Arctic Archipelago (island syndrome) and northern of the Boothia Peninsula [94, 119]. They are adapted to the Arctic environment having a compact body size for conserving heat (Allen's rule) [120] and to the

<table><caption>Table 16. MANOVA pairwise comparisons results (mobility-forelimb phalanx).</caption><thead><tr><th>Mobility</th><th>Pillai</th><th>approx F</th><th>num DF</th><th>den DF</th><th>Pr(&gt;F)</th></tr></thead><tbody><tr><td>Altitudinal—Planarly</td><td>0.73395</td><td>3.9849</td><td>9</td><td>13</td><td><strong>0.03682</strong></td></tr><tr><td>Altitudinal—Sedentary</td><td>0.49471</td><td>1.4142</td><td>9</td><td>13</td><td>0.82844</td></tr><tr><td>Planarly—Sedentary</td><td>0.76317</td><td>4.6545</td><td>9</td><td>13</td><td><strong>0.01933</strong></td></tr></tbody></table>

Results of the MANOVA pairwise comparisons for *Rangifer tarandus caribou* (forelimb phalanx) according to mobility pattern (sedentary, planarly, altitudinal). A significant contribution was considered for *P* value <0.05 **(in bold)**.

https://doi.org/10.1371/journal.pone.0285487.t016
<table><caption>Table 17. MANOVA pairwise comparisons results (habitat-hindlimb phalanx).</caption><thead><tr><th>Habitat</th><th>Pillai</th><th>approx F</th><th>num DF</th><th>den DF</th><th>Pr(&gt;F)</th></tr></thead><tbody><tr><td>BorealForest—Mountain</td><td>0.69207</td><td>4.1204</td><td>12</td><td>22</td><td><strong>0.005958</strong></td></tr><tr><td>BorealForest—Tundra</td><td>0.67136</td><td>3.7452</td><td>12</td><td>22</td><td><strong>0.010688</strong></td></tr><tr><td>Mountain—Tundra</td><td>0.62338</td><td>3.0345</td><td>12</td><td>22</td><td><strong>0.034812</strong></td></tr></tbody></table>

Results of the MANOVA pairwise comparisons for *Rangifer tarandus caribou* (hindlimb phalanx) according to habitat type (mountain, tundra, boreal forest). A significant contribution was considered for *P* value <0.05 **(in bold)**.

https://doi.org/10.1371/journal.pone.0285487.t017

limited plant growth with a highly compressed growing season and long periods of snow-cov- ered frozen standing vegetation [119]. All these facts explain their lower centroid sizes.

At an intra-subspecies scale, size differences are not significant between populations of *R.t. caribou*. Forelimb phalanges of specimens that perform altitudinal migration tend to be slightly larger than those belonging to planarly or sedentary migratory behaviour, but it is yet unknown if this is linked to real population differences or unbalanced sex ratios in our sample.

If confirmed, these size differences might be related to the migration costs. Migration of long distances might increase access to high quality food or reduce the risk of predation but it also involves an energetic cost that may be reflected in body conditions of animals [121, 122]. The proportional cost of horizontal or planarly travel that tundra caribou performs increases with decreasing body size and long movements, particularly in winter, have a negative impact on caribou body size [123, 124]. However, at present, there is no evidence of a statistically sig- nificant influence of habitat or mobility on phalanx size.

Size differences tend to be more pronounced on forelimb than hindlimb members. This could be explained by the fact that forelimb supports more stress as they carry greater share of the body weight [75–77, 125–128]. The use of forelimb more actively is also linked to the search of food, breaking through ice layers in many cases [69, 70, 75]. It is, however, important to keep in mind that these results do not satisfy the thresholds of statistical significance, and therefore should be interpreted carefully.

## Shape variation: The impact of subspecies, substrate type and snow cover?

Despite the fact that significant shape differences exist between populations, it is difficult to identify a specific morphology that corresponds to a given habitat or mobility type.

If we consider CVA, in terms of classification instead of associated morphological descrip- tion, CVA is specifically designed to identify linear combinations of variables that best dis- criminate between groups. It maximizes the separation between groups while minimizing the variation within groups, and thus can be more effective in identifying the underlying differ- ences between groups. For forelimb phalanges, the results suggest that habitat may be a

<table><caption>Table 18. MANOVA pairwise comparisons results (mobility-hindlimb phalanx).</caption><thead><tr><th>Mobility</th><th>Pillai</th><th>approx F</th><th>num DF</th><th>den DF</th><th>Pr(&gt;F)</th></tr></thead><tbody><tr><td>Altitudinal—Planarly</td><td>0.67123</td><td>3.7430</td><td>12</td><td>22</td><td><strong>0.01072</strong></td></tr><tr><td>Altitudinal—Sedentary</td><td>0.57331</td><td>2.4633</td><td>12</td><td>22</td><td>0.09655</td></tr><tr><td>Planarly—Sedentary</td><td>0.57472</td><td>2.4776</td><td>12</td><td>22</td><td>0.09405</td></tr></tbody></table>

Results of the MANOVA pairwise comparisons for *Rangifer tarandus caribou* (hindlimb phalanx) according to mobility pattern (sedentary, planarly, altitudinal). A significant contribution was considered for *P* value <0.05 **(in bold)**.

https://doi.org/10.1371/journal.pone.0285487.t018
<table><caption>Table 19. MANOVA pairwise comparisons results (considering all the categories).</caption><thead><tr><th></th><th>Pillai</th><th>approx F</th><th>num DF</th><th>den DF</th><th>Pr(&gt;F)</th></tr></thead><tbody><tr><td>Rtcaribou_BorealForest_Sedentary—Rtcaribou_Mountain_Altitudinal</td><td>0.68956</td><td>3.8871</td><td>12</td><td>21</td><td><strong>0.01932</strong></td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtcaribou_Mountain_Sedentary</td><td>0.52888</td><td>1.9645</td><td>12</td><td>21</td><td>0.50795</td></tr><tr><td>Rtcaribou_BorealForest_Sedentary—Rtcaribou_Tundra_Planarly</td><td>0.67718</td><td>3.6709</td><td>12</td><td>21</td><td><strong>0.02695</strong></td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtcaribou_Mountain_Sedentary</td><td>0.34718</td><td>0.9307</td><td>12</td><td>21</td><td>1.00000</td></tr><tr><td>Rtcaribou_Mountain_Altitudinal—Rtcaribou_Tundra_Planarly</td><td>0.67132</td><td>3.5743</td><td>12</td><td>21</td><td><strong>0.03136</strong></td></tr><tr><td>Rtcaribou_Mountain_Sedentary—Rtcaribou_Tundra_Planarly</td><td>0.37552</td><td>1.0523</td><td>12</td><td>21</td><td>1.00000</td></tr></tbody></table>

Results of the MANOVA pairwise comparisons for hindlimb proximal phalanx according to the different categories. A significant contribution was considered for *P* value <0.05 (in **bold**).

https://doi.org/10.1371/journal.pone.0285487.t019

stronger predictor of classification accuracy than mobility pattern. However, on hindlimb phalanges, these ones indicate that both habitat type and mobility pattern are important predictors of the observed patterns.

Thus, these outcomes suggest that the morphology of caribou forelimb phalanges is more heavily influenced by the type of habitat in which they inhabit than their mobility pattern. This could mean that the forelimbs are adapted for specific types of terrain (e.g. rocky and mountainous terrain). In contrast, the results for the hindlimb phalanges indicate that both habitat type and mobility pattern are significant predictors of the observed patterns, suggesting that they are better suited for more dynamic activities (e.g. running). Thus, the different predictors of classification accuracy between the two limb sets likely reflect the different functional roles of each set of phalanges. The results are consistent with the MANOVA results, which also indicate that habitat variable plays an important role in the configuration of the shape of the phalanx morphology.

Results on forelimb phalanx revealed, along PC1, on one hand, a joint with a more massive morphology towards cranial side and medio-laterally narrower (mainly present in Rtcaribou_BForest_Sed and partially on Rtcaribou_Tundra_Plan), and on the other hand, proximal joints that appeared to be more pronounced mediolaterally rather than craniocaudally (mostly

<figure><img src="image_15.png" /><figcaption>**Fig 12.** Multivariate regression plots. Multivariate regression plots performed on shape data (regression scores) and log-transformed centroid size (Log-Centroid Size) according to subspecies (*R.t. caribou*, *R.t. granti*, *R.t. groenlandicus* and *R.t. pearyi*); Forelimb (left) and hindlimb (right) phalanges.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0285487.g012
observed in mountain *R.t. caribou* and *R.t. granti*, *R.t. groenlandicus* and *R.t. pearyi*). There is however no clear distinction according to habitat or mobility. In spite of the difficulty of establishing the functional implications of these morphological differences, it could suggest adaptations to different types of movement or stresses placed on the forelimb phalanges in different environments. For example, the more massive joint towards the cranial side and medio-laterally narrower morphology observed in individuals inhabiting the boreal forest may provide greater stability when moving through snow, whereas the more pronounced mediolateral morphology observed in mountain individuals may provide greater flexibility when traversing uneven terrain.

Towards the maximum configuration of PC2 we observe a morphology that is wider in the axial-abaxial axe and narrower craniocaudally (especially for sedentary *R.t. caribou* from the mountains, and also for planarly migrating *R.t. granti* from the tundra) and the opposite configuration along negative values of PC2, with specimens that are both sedentary or migratory, inhabiting mountain, boreal forest or tundra. At this point, PC2 is playing a smaller role in determining the shape of the joint.. In other words, there is no clear distinction between the morphology of populations based on their migration behaviour or habitat.

Forelimb phalanx from sedentary *R.t. caribou* inhabiting the mountains overlaps to a lesser extent in the morphospace with the other populations (especially when the *R.t. caribou* subspecies is analyzed separately) and reflects a more characteristic morphology, contrasting with the boreal forest population, which is also sedentary. The mountain caribou's forelimb phalanx is likely to be more robust and adapted to their mountainous environment, whereas the boreal forest population's forelimb phalanx is likely to be more adapted to the conditions found in the boreal forest (e.g. a more varied terrain, with greater amounts of wetland; a more diverse array of vegetation, including shrubs, grasses, and forbs; and a greater variety of food source). In both cases, it is worth clarifying that when we refer to sedentary populations we do so to talk about those that travel distances less than 200 km annually [see 41], so in this case substrate and snow coverage potentially have an important impact on the morphological configuration of the joint of the proximal forelimb phalanx. Our sample of boreal forest caribou comes mainly from Québec (Canada), where the forests are covered in snow more than half of the year [88]. These types of boreal forest populations also seek places with snow that can be used as refuges from predators. Sedentary montane populations, however, remain in alpine areas to feed on lichens in windswept areas where snow is absent or shallow [48]. These differences in substrate type and snow cover may also explain the slight overlap in morphological configuration between sedentary and altitudinally migrating mountain individuals: even though both inhabit mountainous ecosystems, they move in different ways (altitudinal migration vs. sedentary). While for the sedentary mountain caribou we have already discussed the characteristics of their habitat, for those that perform altitudinal migration, at the end of autumn, when the snow begins, they migrate from upper to lower elevations where the snow is shallow enough, avoiding open habitats where the snow is harder, again beginning its travel to higher elevations around February towards subalpine forests [48, 129]. This suggests that the functional morphology of this particular subspecies (*Rtcaribou_Mountain_Sedentary*) is well adapted to its environment, allowing to thrive in its mountainous habitat (characterized by a diverse range of elevations, and the presence of high-altitude alpine tundra, ranging from rolling hills and valleys, to steep rocky slopes, and high peaks).

Morphological differences are even harder to classify on hindlimb phalanges. All populations are distributed along PC1, and PC2 poorly distinguishes populations, indicating that the morphology of the hindlimb phalanx of each subspecies is largely dependent on the habitat and mobility pattern of the subspecies. As a result, the morphological variation of hindlimb phalanx is not so evident to interpret.
MANOVA tests showed statistically significant differences in relation to Habitat for both fore- and hindlimb phalanges, and in relation to Subspecies and Mobility for hindlimb phalanx. Shape differences were more pronounced on hindlimb than forelimb phalanges, and more marked when pooled by Habitat in both cases. This may suggest that the habitat in which *Rangifer tarandus* lives, has a significant impact on the joint shape of forelimb first phalanges, as has been proven in different Bovidae species [60, 67, 130]. To understand habitat's influence, it is important to note the variety of environments and substrates that caribou inhabit, ranging from closed boreal forest to treeless arctic tundra passing through the mountain with a mosaic of landscapes ranging from open alpine areas to subalpine areas and low-elevation forests [69]. There is no doubt that the snow adaptation impacts the shape of the proximal joint, not only during migration, but also when the front hooves are used to dig through the snow to find lichen during snowy seasons [68].

**Shape variability at the intra-subspecies scale: The key influence of habitat?**

Based on pairwise comparisons, in hindlimb phalanx, main differences on shape are found between the different subspecies, a fact that motivated the realization of the same analyses within the *Rangifer tarandus caribou* subspecies. Such intra-subspecies analysis confirmed that phalanx shape is statistically different according to Habitat and Mobility, even if only *R.t. caribou* specimens are considered.

When *R. t. caribou* was examined separately, MANOVA showed differences in habitat and mobility, as well as their interaction, in fore- and hindlimb phalanges. In pairwise comparisons, all habitat categories showed significant differences, regardless of whether they were fore- or hindleg and demonstrating (as determined by the p values) how crucial the impact of habitat on shape is. Regarding mobility, there were differences for forelimb phalanges between planarly on one hand, altitudinal and sedentary on the other, and only between altitudinal and planarly in the hindlimb ones. Considering the interaction, pairwise comparisons did not yield significantly different results in foreleg phalanx, but in hindleg, between *Rtcaribou_BForest_Sed* and *Mountain_Altitudinal* and *Tundra_Planarly* (respectively). Consequently, habitat has again the greatest impact on the shape, although to a lesser extent the type of mobility also has an impact. The altitudinal migration in the mountains occurs vertically, through rugged and mountainous substrates in response to snowfall conditions, forage availability, and to avoid predators. Migration in tundra occurs horizontally, with a characteristic vegetation of dwarf and upright shrubs [94, 129, 131]. These differences are not only in the substrate and the type of movement, but also in the distance traveled by the herds. *Rtcaribou_Tundra_Planarly* travels distances between 1,220 and 1,770 km annually, while *Rtcaribou_Mountain_Altitudinal* can travel between 140-240 km only in the spring migration [81, 90, 92, 131]. Furthermore, the marked difference in the hindlimb phalanx may also be explained by the implications of propulsion (along with what was previously discussed), resulting in a higher level of stress among individuals who move vertically instead of horizontally due to the substrate conditions, in addition to the fact that in quadrupeds, as it has been seen in other species, e.g. rhino, horse [126, 127, 132] the hind limbs bear relatively less weight than the fore limbs and play an additional propulsive role during locomotion.

The results show that subspecies structure the data significantly, especially for size, but habitat and mobility also have significant effects on phalanx morphology. Shape analyses on *Rangifer tarandus caribou* subspecies yielded shape differences that are especially significant in both fore and hindlimb phalanges according to habitat or mobility categories. Consequently, if the subspecies variable is excluded, we can examine how habitat has a greater impact on the morphological configuration. This can be explained by the wide variety and different habitats
used by *Rangifer tarandus caribou*: those that inhabit mountain environments move among windswept alpine ridges, subalpine forests, and low-elevation forests [69]. In the Boreal Forest, in which the snow is a key feature, they seek habitats with thin and soft snow cover with a more sedentary lifestyle [88]. Those inhabiting the tundra will find many types of flowering plants, mosses, and lichens, as well as shrubs, grasses, and sedges which form swards and tussocks scattered across the tundra [94]. Thus, e.g. Rtcaribou_BForest_Sed tends to have a forelimb phalanx joint morphology with a more cranially prominent axial articular surface, compared with sedentary mountain individuals. This morphology could be linked to the adaptations that enable it to better handle the specific demands of the boreal forest environment (e.g. deep snow and the ability to dig through the snow to reach vegetation; rugged terrain and extreme cold). For example, the cranially prominent axial articular surface may allow for a greater range of motion in the joint, or it may provide additional surface area for muscle attachment, allowing for more precise and powerful movements. However, further research would be necessary to fully understand the functional implications of this morphological difference.

## Preliminary distribution on sexual dimorphism

Even though the sample size by sex is limited for both forelimb and hindlimb phalanges, there are no differences between male and female specimens from the same subspecies (*R. t. caribou*) when log-centroid sizes were displayed on boxplots to explore their distribution. In this regard, Hull [71] noted in her study to distinguish forelimb and hindlimb phalanges, that, even though differences between sexes were appreciated, they also included a large overlap, so she concluded they cannot be included in the study without additional context. The boxplots show that male phalanges are slightly larger than female ones, but they overlap. These non-significant results do not imply that these differences between the two sexes do not exist in our sample, as has been demonstrated for other anatomical parts [82], it can simply reflect that phalanges are not an indicator of this characteristic and/or a larger sample and further studies are needed.

## Allometry

Allometry is also more significant on forelimb than on hindlimb (absent in our sample) phalanges, similar to what has already been observed in other works. Previous studies on reindeer [75, 77] show that allometry was generally stronger on the forelimbs than the hindlimbs as well. It has been observed in other species as well, for example the rhinoceros [127, 128] and some bovids [130]. The fore-phalanx allometry is not very significant, but it may be compensated by the rest of the forelimbs. Thus, a predictable morphological and biomechanical response of the bone is necessary to resist the stresses caused by body mass, which is heavier (as noted above) on forelimbs due to antlers (which are heavier in males), resulting in additional forces and constraints on the foreleg bones.

## Conclusion

Reconstructing the interactions between Palaeolithic human groups and reindeer is often complicated by a general lack of data on the ethology and ecology of these past animal communities. At the end of the Late Glacial, reindeer gradually disappeared from southwestern Europe with the climate warming. As they became more fragmented geographically, different responses should be expected, such as variation in their migratory patterns. As the migration behaviours of reindeer herds during this period are largely unknown, it is difficult to interpret Late Glacial socio-economic organization and mobility strategies.
Due to the fact that locomotor morphological adaptations are closely related to habitat preferences and mobility [67, 78, 133, 134], it is possible to infer behavioural responses using an animal's cranial and post-cranial morphology, in the case of our study, using first phalanx. Besides, phalanges are one of the most useful bones for palaeocological purposes because these elements are commonly preserved in fossil assemblages. Together with 3D geometric morphometrics, we were able to examine their morphology from an anatomical region of interest and visualize morphological variation. The present study has contributed to our understanding of the caribou proximal phalanx's functional morphology and its relationship to mobility and habitat type (adapted to provide support on uneven terrain and harsh habitats), but it still cannot be easily applied to archaeological bone remains. A complex interplay of factors has been shown, with the impact on phalanx morphology of subspecies status, sexual dimorphism, substrate type, snow cover, migration costs, habitat, etc. This diversity of interrelated factors makes interpretation difficult with the limited sample of modern reindeer phalanges currently available. Continued efforts in the 3D scanning of phalanges should lead to further progress.

Although future studies with larger sample size are needed, this contribution represents a necessary first step to explore reindeer phalanx morphology and size differences for potential archaeological applications.

**Supporting information**

**S1 Table. Specimens information.** Detailed information related to Rangifer tarandus sample.
(DOCX)

**S1 File. Kappa parameter concept.** Detailed description of this concept.
(DOCX)

**S2 File. Boxplots on log-transformed centroid size.** Additional boxplots exploring size variation.
(DOCX)

**Acknowledgments**

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 would also like to thank N. Torres-Tamayo for her helpful comments and the anonymous reviewers whose suggestions have undoubtedly improved this manuscript. ABGL would like to thank A. Martín Esteban and N. Lupiañez Corpas for their valuable comments and support.

**Author Contributions**

**Conceptualization:** Ana Belén Galán López, Maxime Pelletier, Emmanuel Discamps.

**Data curation:** Ana Belén Galán López.

**Formal analysis:** Ana Belén Galán López.

**Funding acquisition:** Ana Belén Galán López, Emmanuel Discamps.

**Investigation:** Ana Belén Galán López.
**Methodology:** Ana Belén Galán López, Maxime Pelletier, Emmanuel Discamps.

**Project administration:** Emmanuel Discamps.

**Resources:** Ana Belén Galán López.

**Supervision:** Emmanuel Discamps.

**Validation:** Ana Belén Galán López.

**Writing – original draft:** Ana Belén Galán López.

**Writing – review & editing:** Ana Belén Galán López, Maxime Pelletier, Emmanuel Discamps.

## References

1. Lartet E. Nouvelles recherches sur la coexistence de l'homme fossile et des grands mammifères fosiles réputés caractéristiques de la dernière période géologique. Annales des Sciences naturelles. II. Zoologie. 4e série. 1861; XV: 177-253.

2. Discamps E, Jaubert J, Bachellerie F. Human choices and environmental constraints: deciphering the variability of large game procurement from Mousterian to Aurignacian times (MIS 5-3) in southwestern France. Quat Sci Rev. 2011; 30: 2755-2775.

3. Costamagno S, Barshay-Szmidt C, Kuntz D, Laroulandie V, Pétillon JM, Boudadi-Maligne M, et al. Reexamining the timing of reindeer disappearance in southwestern France in the larger context of late glacial faunal turnover. Quat Int. 2016; 414: 34-61.

4. Gaudzinski S, Roebroeks W. Adults only. Reindeer hunting at the Middle Palaeolithic site Salzgitter Lebenstedt, Northern Germany. J Hum Evol. 2000; 38: 497-521. https://doi.org/10.1006/jhev.1999.0359 PMID: 10715194

5. Gaudzinski S, Street M. Reconsidering hunting specialisation in the German Magdalenian faunal record. In: Costamagno S, Laroulandie V, editors., Mode de vie au Magdalenien: les apports de l'archeozoologie. Actes du XIVe Congres UISPP, 2-8 septembre 2001, Liège, BAR International Series 1144. Archaeopress, Oxford. 2003. pp. 11-21.

6. Costamagno S. Strategies de chasse et fonction des sites au Magdalenien Dans le sud de la France. PhD Thesis. Université Bordeaux I. 1999.

7. Costamagno S, Meignen L, Maureille B, Vandermeersch B. Les Pradelles (Marillac-le-Franc, France): a Mousterian reindeer hunting camp? J Anthropol Archaeol. 2006; 25: pp.466-484.

8. Castel JC, Kuntz D, Chauvière FX, Gerbe M, Juillar F. L'exploitation du monte animal au Paléolithique supérieur en Quercy. Actes de la session C67, XV UISPP International Congress, Lisbon, sept. 2006, Paleo. 2013; 4: 395-418.

9. Niven L, Steele TE, Rendu W, Mallye JB, McPherron S, et al. Neandertal mobility and large-game hunting: The exploitation of reindeer during the Quina Mousterian at Chez-Pinaud Jonzac (Charente-Maritime, France). J Hum Evol. 2012; 63:624-635.

10. Sekhr A. Étude archéozoologique des niveaux aurignaciens (couches 14 à 6) et de la base des niveaux gravettiens (niveaux X à T4) de l'abri Pataud (Les Eyzies, Dordogne). Thèse de doctorat, Muséum National d'Histoire Naturelle de Paris. 1998.

11. Letourneux C. Etude taphonomique et archeozoologique des niveaux aurignaciens de la grotte des Hyères (Brassempouy, Landes): qui a mangé quoi? Archeologie des Pyrenees Occidentales. 2005; 2: 85-102.

12. Soulier MC. Entre alimentaire et technique: l'exploitation animale aux debuts du Paleolithic supérieur. Stratégies de subsistance et chaînes opératoires du traitement du gibier à Isturitz, La Quina aval, Roc-de-Combe et Les Abeilles. Ph.D. thesis. Université de Toulouse 2-Le Mirail, France. 2013.

13. Delpech F. L'environnement animal des Européens au Paléolithique supérieur. In Desbrosse R, Thévenin A, editors. Préhistoire de l'Europe. Des origines à l'Age du Bronze Actes des Congrès nationaux des Sociétés historiques et scientifiques, Lille 2000. CTHS, Paris. 2003. 125: pp. 271-289.

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

15. Lacarrière J. Les ressources cynégétiques au Gravettien en France. Acquisition et modalités d'exploitation des animaux durant la phase d'instabilité climatique précédant le dernier maximum glaciaire, Thèse, Université Toulouse Jean Jaurès. 2015.

16. Enloe JG, David F. Le remontage des os par individus: le partage du renne chez les magdaléniens de Pincevent. Bull. de la Soc. Prehist. Francaise. 1989; 86: 275-281.
17. Enloe JG. Subsistence organization in the Early Upper Paleolithic: reindeer hunters of the Abri du Fla-geolet, couche V. In: Knecht H., Pike-Tay A, White R, editors. Before Lascaux: the Complex Record of the Early Upper Paleolithic. CRC Press. 1993. pp. 101–115.

18. Enloe JG. Seasonality and age structure in remains of Rangifer tarandus: Magdalenian hunting strat-egy at Verberie. Anthropozoologic. 1997; 25: 95–102.

19. Enloe JG. Acquisition and processing of reindeer in the Paris Basin. In: Costamagno S, Laroulandie V, editors. Mode de vie au Magdalénien: les apports de l'archeozoologie. Actes du XIVe Congrès UISPP, 2-8 septembre 2001, Liège, BAR International Series 1144. Archaeopress, Oxford. 2003; pp. 23–31.

20. Fontana L. La chasse au Renne au Paléolitique Supérieur: nouvelles voies de recherche. Páleo. 2000; 12: 141–164.

21. Bignon O. Chasser les chevaux à la fin du Paléolithique dans le Bassin parisien: Stratégie cynégétique et mode de vie au Magdalénien et à l'Azilien ancien. Vol. 1747. BAR International Series 1747. Oxford: Archaeopress, 2008.

22. Kuntz D. Données nouvelles sur les stratégies d'acquisition et de traitement des carcasses au cours du Magdalénien supérieur dans la vallée de l'Aveyron: l'exemple de La Magdeleine La Plaine (Penne, Tarn), Préhistoire du Sudouest. 2006; 13: 151–166.

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

24. Lacorre F. Les Migrations des Rennes dans la Province Préhistorique des Eyzies. Bull. de la Soc. Prehist. Francaise 1956; 53:302-310.

25. Bahn PG. Seasonal migration in South-West France during the Late Glacial Period. J Archaeol Sci. 1977; 4: 245–257.

26. Gordon B. Des Hommes et des Rennes dans la Préhistoire Française Magdalénienne: Résultats. Archaeozoologica. 1988a; 2: 227-242.

27. Gordon B. Of Men and Reindeer Herds in French Magdalenian Prehistory, vol. 390. Oxford: BAR International Series. 1988b.

28. Gordon B. More on herd-following hypothesis. Curr Anthropol. 1990; 31: 399–400. http://www.jstor.org/stable/2743265

29. Delpech F. L'environnement animal des magdaléniens, In: Otte M, editor. Le Magdalénien en Europe, Liège, ERAUL, Actes du Colloque de Mayence, XI congrès U.I.S.P.P. 1987. pp.5-30.

30. Kuntz D, Costamagno S. Relationships between reindeer and man in southwestern France during the Magdalenian Quat Int. 2011; 238: 12–24.

31. Bouchud J, Cheynier A, Guillien Y. Dents de Renne et migrations. Bulletin de la Société Prehistorique Francaise. 1953.

32. Bouchud J. Dents de Rennes, Bois de Rennes et migrations. Bulletin de la Société Prehistorique Francaise. 1954.

33. Spiess A. Reindeer and caribou hunters—an archaeological study. New York: Academic Press. 1979.

34. Deplano S. Etude de la faune des grands mammifères de la couche IX de l'Abri du Flageolet3 Dor-dogne. Approche taphonomique et palethnographique, Memoire de Maitrise en Prehistoire, Universi-téde Paris. 1994.

35. Fontana L. 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, editors. L'archéologie à découvert. Paris: CNRS Editions. 2012. pp.72-78.

36. Fontana L. 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. J Archaeol Sci Rep. 2017; 12: 346–362.

37. Fontana L, Chaurière FX. Economie et nomadisme au Pléniglaciaire supérieur et au Tardiglaciaire en Europe de l'ouest: le Systeme Renne, In: Djindjian F, editor. La préhistoire de la France, Harmann Edi-teurs, Histoire et archéologie. 2018. pp 141–148.

38. Kuntz D. 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 Tou-louse le Mirail—Toulouse II. 2011. Available from: https://theses.hal.science/tel-00725310

39. Galán López AB, Burke A, Costamagno S. The ecomorphology of Caribou (Rangifer tarandus): a geo-metric morphometric study. Open Res Eur. 2021; 1:99 https://doi.org/10.12688/openreseurope.13782.1
40. Galán López AB, Costamagno S., Burke A. Evaluer le degré de mobilité des rennes paléolithiques à partir de l'étude ostéométrique de métacarpiens: intérêt de l'approche actualiste et des algorithmes d'apprentissage automatique. In: Nicoud E, Balasse M, Desclaux E, Théry-Parisot I, editors. Biodiversités, environnements et sociétés depuis la Préhistoire: nouveaux marqueurs et approches intégrées. Actes 41es Rencontres Internationales d'archéologie et d'histoire de Nice Côte d'Azur. Octobre 2020. CNRS. 2021.

41. Galán López AB, Costamagno S, Burke A. Osteometric Study of Metapodial Bones and Phalanges as Indicators of the Behavioural Ecology of Modern Reindeer (Rangifer tarandus) and Implications for Reconstruction of Paleo Mobility. Open Quat, 2022; 8(1): p.8. https://doi.org/http%3A//doi.org/10.5334/oq.106

42. Festa-Bianchet M, Ray JC., Côté SD, Gunn A. Conservation of caribou (Rangifer tarandus) in Canada: an uncertain future. Can J Zool. 2011; 89: 419-434.

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

44. Bergerud AT. Caribou. In: Demaris S, Krausman PR, editors. Ecology and management of large mammals in North America. New Jersey: Prentice-Hall; 2000. pp. 658-693.

45. Miller FL. Inter-island movements of Peary caribou: a review and appraisal of their ecological importance. In: Harington CR, editor. Canada's missing dimension: science and history in the Canadian Arctic Islands. Canadian Museum Nature, Ottawa, Ontario.1990; Vol. 2. pp. 608-632.

46. Schaefer JA, Bergman CM, Luttich SN. Site fidelity of female caribou at multiple spatial scales. Landscape Ecology. 2000; 15: 731-739.

47. Goddard AD. Boreal caribou in Northeastern British Columbia. Peace Region Technical Report. 2009.

48. Seip D, McLellan B. Mountain Caribou. In: Hummel M, Ray JC, editors. Caribou and the North: a shared future. Toronto: Dundurn Press. 2010. pp 289-306.

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

50. Genty D, Blamart D, Ouahdi R, Gilmour M, et al. Precise dating of Dansgaard-Oeschger climate oscillations in western Europe from stalagmite data. Nature. 2003; 421, 833-837, https://doi.org/10.1038/nature01391 PMID: 12594510

51. Sanchez Goñi MF, Landais A, Fletcher WJ, Naughton F, Desprat S, Duprat JM. Contrasting impacts of Dansgaard-Oeschger events over a western European latitudinal transect modulated by orbital parameters. Quat Sci Rev. 2008; 27(11-12): 1136-1151 https://doi.org/10.1016/j.quascirev.2008.03.003

52. Discamps E, Royer A. Reconstructing palaeoenvironmental conditions faced by Mousterian hunters during MIS 5 to 3 in southwestern France: A multi-scale approach using data from large and small mammal communities. Quat Int. 2017; 433: 64-87. https://doi.org/10.1016/j.quaint.2016.02.067

53. Sommer RS, Kalbe J, Ekström J, Benecke N, Liljegren R. Range dynamics of the reindeer in Europe during the last 25,000 years. J Biogeogr. 2014; 41:298-306. https://doi.org/10.1111/jbi.12193

54. Plummer TW, Bishop LC. Hominid paleoecology at Olduvai Gorge, Tanzania, as indicated by antelope remains. J Hum Evol. 1994; 27: 47-75.

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

56. DeGusta D, Vrba E. A method for inferring paleohabitats from the functional morphology of bovid astragali. J Archaeol Sci.2003; 30: 1009-1022.

57. Bignon O, Baylac M, Vigne JD, Eisenmann V. Geometric morphometrics and the population diversity of Late Glacial horses in Western Europe (Equus caballus arcelini): phylogeographic and anthropological implications. J Archaeol Sci. 2005; 32: 375-391.

58. Van Asperen EN. Ecomorphological adaptations to climate and substrate in late Middle Pleistocene caballoid horses. Palaeogeogr Palaeoclimatol Palaeoecol. 2010; 297: 584-596. https://doi.org/10.1016/j.palaeo.2010.09.007

59. Spencer L. Dietary adaptations of Plio-Pleistocene Bovidae: implications for hominid habitat use. J Hum Evol. 1997; 32: 201-228. https://doi.org/10.1006/jhev.1996.0102 PMID: 9061557

60. Kappelman J. Morphology and locomotor adaptations of the bovid femur in relation to habitat. J Morphol. 1988; 198(1):119-30. https://doi.org/10.1002/jmor.1051980111 PMID: 3199446

61. Kappelman J, Plummer T, Bishop L, Duncan A, Appleton S. Bovids as indicators of Plio-Pleistocene paleoenvironments in east Africa. J Hum Evol. 1997; 32:229-56. https://doi.org/10.1006/jhev.1996.0105 PMID: 9061558
62. Bishop L. Pigs and the ancestors: hominids, suids and environments during the Plio-Pleistocene of East Africa. Ph.D. Dissertation, Yale University.1994.

63. DeGusta D, Vrba E. Methods for inferring paleohabitats from discrete traits of the bovid prostcranial skeleton. J Archaeol Sci. 2005; 32: 1115–1123.

64. DeGusta D, Vrba E. Methods for inferring paleohabitats from the functional morphology of bovid phalanges. J Archaeol Sci. 2005; 32: 1099–1113.

65. Kovarovic K, Andrews P. Bovid postcranial ecomorphological survey of the Laetoli paleoenvironment. J Hum Evol. 2007; 52(6):663–80. https://doi.org/10.1016/j.jhevol.2007.01.001 PMID: 17353031

66. Plummer TW, Bishop LC, Hertel F. Habitat preference of extant African bovids based on astragalus morphology: operationalizing ecomorphology for palaeoenvironmental reconstruction. J Archaeol Sci.2008; 35: 3016–3027.

67. Curran SC. Expanding ecomorphological methods: geometric morphometric analysis of Cervidae post-crania. J Archaeol Sci. 2012; 39: 1172–1182.

68. Hull E, Semeniuk M, Puolakka HL, et al. Tendons and ligaments of the Rangifer tarandus metapodial and hoof. Polar Biol. 2021; 44: 1803–1816. https://doi.org/10.1007/s00300-021-02919-z PMID: 34720376

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

70. Niinimäki S, Salmi AK. Entheseal changes in free-ranging versus zoo reindeer—Observing activity status of reindeer. Int J Osteoarchaeol. 2016; 26(2): 314–323.

71. Hull EH. Metric and non-metric guides for the determination between fore-and hindlimb phalanges of Rangifer tarandus. Rangifer. 2019; 39(1): 11–26.

72. Hull E, Niinimäki S, Salmi AK. Differences in entheseal changes in the phalanges between ecotypes of Fennoscandian reindeer. Int J Osteoarchaeol. 2020; 30: 666–678. https://doi.org/10.1002/oa.2897

73. Salmi AK, Niinimäki S, Pudas T. Identification of working reindeer using palaeopathology and enthe-seal changes. Int J Paleopathol. 2020; 30: 57–67. https://doi.org/10.1016/j.ijpp.2020.02.001 PMID: 32464526

74. Nomokonova T, Losey RJ, Kosintsev PA, Plekhanov AV. Reindeer Demographics at larte VI, lamal Peninsula, Arctic Siberia, Environmental Archaeology. 2021; https://doi.org/10.1080/14614103.2021.1995260

75. Pelletier M, Kotiaho A, Niinimäki S, et al. Identifying early stages of reindeer domestication in the archaeological record: a 3D morphological investigation on forelimb bones of modern populations from Fennoscandia. Archaeol Anthropol Sci. 2020; 12, 169. https://doi.org/10.1007/s12520-020-01123-0 PMID: 32704330

76. Pelletier M, Niinimäki S, Salmi AK. Influence of captivity and selection on limb long bone cross-sec-tional morphology of reindeer. J Morphol. 2021; 282(10):1533–1556. https://doi.org/10.1002/jmor.21403 PMID: 34323317

77. Pelletier M, Kotiaho A, Niinimäki S, Salmi AK. 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. 2022; 34:5, 802–820, https://doi.org/10.1080/08912963.2021.1947268

78. Curran SC. Hominid paleoecology and cervid ecomorphology. Ph.D. University of Minnesota. 2009.

79. Viscosi V, Cardini A. Leaf morphology, taxonomy and geometric morphometrics: A simplified protocol for beginners. PLoS One; 2011: e25630. https://doi.org/10.1371/journal.pone.0025630 PMID: 21991324

80. Hanot P, Guintard C, Lepetz S, Cornette R. Identifying domestic horses, donkeys and hybrids from archaeological deposits: a 3D morphological investigation on skeletons. J Archaeol Sci. 2017; 78:88–98. https://doi.org/10.1016/j.jas.2016.12.002

81. Saher D.J & Schmiegelow FK. Movement pathways and habitat selection by woodland caribou during spring migration. Rangifer. 2005; 25: 143–154.

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

83. Greenfield HJ. Sexing fragmentary Ungulate Acetabulae. In: Ruscillo D, editor. Recent Advances in Ageing and Sexing Animal Bones. Proceeding of the 9th ICAZ Conference. Durham. Oxbow Press, Oxford, UK.2002. pp. 68–86.

84. Galán López AB, Domínguez-Rodrigo M. A biometric analysis of the pelvic acetabulum as an indicator of sex in bovids. C R Palevol. 2014; 13: 561–567. https://doi.org/10.1016/j.crpv.2014.04.003
85. Cavedon M, Gubili C, Heppenheimer E, et al. Genomics, environment and balancing selection in behaviourally bimodal populations: The caribou case. Mol Ecol. 2019; 28: 1946–1963. https://doi.org/10.1111/mec.15039 PMID: 30714247

86. Miller FL. Caribou. In: Wild Mammals of North America: Biology, Management and Conservation. George A. Feldhamer, Bruce C. Thompson, Joseph A. Chapman, editors. Baltimore: Johns Hopkins University Press. 2003.

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

88. Schaefer JA. Boreal Forest Caribou. In: Hummel M, Ray JC, editors. Caribou and the North: a shared future. Toronto: Dundurn Press. 2010. Pp.289–306.

89. Mallory F, Hillis T. Demographic characteristics of circumpolar caribou populations: ecotypes, ecological constraints, releases, and population dynamics. Rangifer. 1998; 18(5): 49–60. https://doi.org/10.7557/2.18.5.1541

90. Leblond M, St-Laurent MH, Côté SD. Caribou, water, and ice—fine-scale movements of a migratory arctic ungulate in the context of climate change. Mov Ecol. 2016; 4: 14. https://doi.org/10.1186/s40462-016-0079-4 PMID: 27099756

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

92. Cumming HG, Beange DB. Dispersion and Movements of Woodland Caribou Near Lake Nipigon, Ontario. J. Wildl. Manage. 1987; 51(1): 69–79. https://doi.org/10.2307/3801634

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

94. Gunn A. Migratory Tundra Caribou. In: Hummel M, Ray JC, editors. Caribou and the North: a shared future. Toronto: Dundurn Press. 2010. pp. 240–267.

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

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

97. Miller FL. Multi-island seasonal home range use by two Peary Caribou, Canadian High Arctic Islands, Nunavut, 1993–94. Arctic. 2002; 55: 133–142.

98. Jenkins D, Campbell M, Hope G, Goorts J, Mcloughlin P. 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. 2011.

99. COSEWIC. 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. Committee on the Status of Endangered Wildlife in Canada. Ottawa. 2004.

100. Fancy SGet.al. Satellite telemetry: a new tool for wildlife research and management. United States Fish and Wildlife Service. 1988: 172: 1–54.

101. Ferguson MAD, Messier F. Traditional winter range: Population Dynamics and Physical Condition. J Wildl Manage. 2000; 64: 168–178.

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

103. Joly K, Gurarie E, Sorum MS, et al. Longest terrestrial migrations and movements around the world. Sci Rep- 2019; 9: 15333. https://doi.org/10.1038/s41598-019-51884-5 PMID: 31654045

104. Adams DC, Rohlf FJ, Slice DE. Geometric morphometrics: ten years of progress following the ‘revolution’. Ital J Zool. 2004; 71(1): 5–16. https://doi.org/10.1080/11250000409356545

105. Zelditch ML, Swiderski DL, Sheets HD, Fink WL. Geometric Morphometrics for Biologists: a Primer. Elsevier Academic Press, Boston. 2004.

106. Bookstein FL. Landmark methods for forms without landmarks: morphometrics of group differences in outline shape. Medical Image Analysis. 1997; 13:225–243. https://doi.org/10.1016/s1361-8415(97)85012-8 PMID: 9873908

107. Gunz P, Mitteroecker P. Semilandmarks: a method for quantifying curves and surfaces. Hystrix. 2013; 24(1):103–109. https://doi.org/10.4404/hystrix-24.1–6292

108. Bookstein FL. Morphometric tools for landmark data: geometry and biology. Cambridge: Cambridge University Press. 1991.

109. Rohlf FJ, Slice DE. Extensions of the Procrustes method for the optimal superimposition of landmarks. Syst Biol. 1990; 39(1):40–59. https://doi.org/10.2307/2992207
110. Bookstein FL. Biometrics, biomathematics and the morphometric synthesis. Bulletin of Mathematical Biology. 1996; 58(2):313-365. https://doi.org/10.1007/BF02458311 PMID: 8713662

111. Dryden IL, Mardia KV. Statistical shape analysis. New-York: Wiley. 1998.

112. Dryden I.L. Shape and object data analysis. Biom. J. 2014; 56: 758-760. https://doi.org/10.1002/bimj.201300220 PMID: 24740922

113. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Stat Method. 1995; 57 (1):289-300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x

114. Wiley DF, Amenta N, Alcantara DA, Ghosh D, Kil YJ, Delson E, et al. Evolutionary morphing. IEEE Visual. 2005; 5: 431-438. https://doi.org/10.1109/VISUAL.2005.1532826

115. Schlager S. Morpho: calculations and vizualisations related to Geometric Morphometrics. Accessed 2020 June 11. https://cran.r-project.org/web/packages/Morpho/Morpho.pdf

116. Kovarovic K, Aiello LC, Cardini A, Lockwood CADiscriminant Function Analysis in archaeology: are classification rates too good to be true? J Archeol Sci. 2011; 38: 3006-3018.

117. Rstudio Team. RStudio: Integrated Development for R. RStudio, PBC, Boston, MA http://www.rstudio.com/. 2020.

118. Adams DC, Otárola-Castillo E. Geomorph: an R package for the collection and analysis of geometric morphometric shape data. Meth Ecol Evol. 2013; 4(4):393-399. https://doi.org/10.1111/2041-210X.12035

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

120. Allen J. A. The influence of physical conditions in the genesis of species. Radic. Rev. 1877; 1: 108-140.

121. Hebblewhite M, Merrill EH. Multiscale wolf predation risk for elk: does migration reduce risk? Oecologia. 2007; 152(2):377-87. https://doi.org/10.1007/s00442-007-0661-y PMID: 17287955

122. Alerstam T, Hedenström A, Åkesson S. Long-distance migration: evolution and determinants. Oikos, 2003; 103: 247-260. https://doi.org/10.1034/j.1600-0706.2003.12559.x

123. Fancy SG, White RG. Energy Expenditures by Caribou while Cratering in Snow, J Wildl Manage. 1985; 49: 987-993.

124. Courturier S, Otto R., Côté SD, Luther G, Mahoney SP. Size variations in caribou ecotypes and relationships with demography. J Wildl Manage. 2010; 74: 395-404.

125. Polly PD. Limbs in mammalian evolution. Chapter 15. In: Hall BK, editor. Fins into Limbs: Evolution, Development, and Transformation. Chicago: University of Chicago Press Books. 2007, pp. 245-268.

126. Barone R. Anatomie comparée des mammifères domestiques. Tome 2: Ostéologie. Paris: Vigot Frères. 2010.

127. Mallet C, Cornette R, Billet G, Houssaye A. Interspecific variation in the limb long bones among modern rhinoceroses—extent and drivers. PeerJ. 2019; 7: e7647. https://doi.org/10.7717/peerj.7647 PMID: 31579585

128. Mallet C, Billet G, Houssaye A, Cornette R. A first glimpse at the influence of body mass in the morphological integration of the limb long bones: an investigation in modern rhinoceroses. J. Anat. 2020; 237: 704-726. https://doi.org/10.1111/joa.13232 PMID: 32519813

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

130. Etienne C, Filippo A, Cornette R, et al. Effect of mass and habitat on the shape of limb long bones: A morpho-functional investigation on Bovidae (Mammalia: Cetartiodactyla) J. Anat. 2021; 238: 886-904. https://doi.org/10.1111/joa.13359

131. Weaver JL. 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. 2006.

132. Hildebrand M. Analysis of Vertebrate Structure. John Wiley & Sons. New York. 1974.

133. Kölher M. Skeleton and Habitat of recent and fossil Rumiants. Verlag Friedrich Pfeil. München. 1993.

134. Kovarovic K. Bovids as palaeoenvironmental indicators. An ecomorphological analysis of bovid post-cranial remains from Laetoli, Tanzania. PhD: University of London. 2004. Available from: https://discovery.ucl.ac.uk/id/eprint/1364862/1/415213.pdf