# A long-term perspective on Neanderthal environment and subsistence: insights from the dental micro-texture analysis of hunted ungulates at Combe-Grenal

(Dordogne, France)

Emilie Berlioz, Eugénie Capdepon, Emmanuel Discamps

► To cite this version:

Emilie Berlioz, Eugénie Capdepon, Emmanuel Discamps. A long-term perspective on Neanderthal environment and subsistence: insights from the dental micro-texture analysis of hunted ungulates at Combe-Grenal (Dordogne, France). PLoS ONE, 2023, 18 (1), e0278395. https://doi.org/10.1371/journal.pone.0278395

HAL Id: hal-03908371

https://hal.science/hal-03908371

Submitted on 14 Jun 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>
# A long-term perspective on Neanderthal environment and subsistence: Insights from the dental microwear texture analysis of hunted ungulates at Combe-Grenal (Dordogne, France)

Emilie Berlioz¹,²*, Eugénie Capdepon¹, Emmanuel Discamps¹

1 UMR5608 TRACES, Team SMP3C, Toulouse, France, 2 Grupo I+D+i EVOADAPTA, Universidad de Cantabria, Santander, Spain

* emilie.berlioz@unican.es

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

## Abstract

Large bovids and cervids constituted major components of the European Middle Palaeolithic faunas and hence a key resource for Neanderthal populations. In paleoenvironmental reconstructions, red deer (*Cervus elaphus*) occurrence is classically considered as a tree-cover indicator while Bovinae (*Bison priscus* and *Bos primigenius*) and reindeer (*Rangifer tarandus*) occurrences are typically associated with open landscapes. However, insights into the ecology of extant ungulate populations show a more complex reality. Exploring the diet of past ungulates allows to better comprehend the hunting strategies of Palaeolithic populations and to reconstruct the modifications through time of past landscapes. By reflecting what animals have eaten during the last days or weeks of their life, dental microwear textures of herbivores link a population and its environment. Here we analyzed, via Dental Microwear Texture Analysis (DMTA), the diet of 50 *Bos/Bison*, 202 *R. tarandus* and 116 *C. elaphus* preyed upon by the Neanderthals that occupied Combe-Grenal rock-shelter, one of the most important Mousterian archaeo-sequences in southwestern France considering its long stratigraphy, abundance of faunal remains and the variations perceptible in Palaeolithic material culture. Grazers and mixed-feeders are the most represented dietary categories among Combe-Grenal's guild of herbivores, highlighting the availability, along the sequence, of open landscapes. The absence of clear changes in the use of plant resources by hunted ungulates through time, even though palaeoenvironmental changes were well-documented by previous studies along the sequence, is interpreted as resulting from the hunting of non-randomly selected prey by Neanderthals, preferentially in open environments. Thus, these results provide further insight into the hunting strategies of Neanderthals and modify our perception of potential links between subsistence and material culture. Combe-Grenal hunters "stayed in the open" through millennia, and were not forced to switch to hunting tactics and material technology adapted to close encounters in forested environments.

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

### OPEN ACCESS

**Citation:** Berlioz E, Capdepon E, Discamps E (2023) A long-term perspective on Neanderthal environment and subsistence: Insights from the dental microwear texture analysis of hunted ungulates at Combe-Grenal (Dordogne, France). PLoS ONE 18(1): e0278395. https://doi.org/10.1371/journal.pone.0278395

**Editor:** Julien Louys, Griffith University, AUSTRALIA

**Received:** July 26, 2022

**Accepted:** November 15, 2022

**Published:** January 18, 2023

**Copyright:** © 2023 Berlioz et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the article and its Supporting Information files (Appendices 1, 2, and 3), as well as in the following NAKALA depository: https://doi.org/10.34847/nkl.f4b0q9rm.

**Funding:** ED: ANR-18-CE03-0007 (funder: Agence Nationale de la Recherche) https://anr.fr/Project-ANR-18-CE03-0007 NO, The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. EB:
L'Oréal-UNESCO FWIS Award 2019 (no grant number) (funder: L'Oréal-UNESCO) https://www. forwomeninscience.com/NO, The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

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

**Introduction**

When considering prehistoric societies whose subsistence relied heavily on the hunting of large mammals, knowledge of the animal communities they preyed upon is key to better con- textualize human–environment interactions [1]. Reindeer (*Rangifer tarandus*), red deer (*Cer- vus elaphus*), steppe bison (*Bison priscus*) and aurochs (*Bos primigenius*) constituted key- resources for hunter-gatherer groups in Eurasia, and were often at the heart of Palaeolithic subsistence economies. Understanding the ethology and ecology of these preyed species (social, spatial, seasonal behavior, feeding ecology) and their long-term evolution is essential to better comprehend human behaviors during the Palaeolithic as these factors ultimately condi- tion human acquisition, subsistence and mobility strategies. It also enlarges our understanding of how the global climatic oscillations have affected local paleoenvironments and human groups.

The impact of habitat anthropization in present days on most modern terrestrial ecosystems, resulting in the use of peri-optimal niches and refugee areas as a consequence of habitat loss and human activities notably for several large mammals [2–4], strongly skew our understanding of the ecological optimum and plasticity of these extant species. It is, however, this ecological knowledge acquired on extant communities that is so often used as habitat proxies to interpret past environmental conditions. This situation leads to undervaluing the potential ecological plasticity of ungulates, leading in turn to flawed paleoecological and paleoenvironmental inferences based on fossil ungulate occurrences [5] and limited understanding of the interactions of these key prey species with prehistoric societies. For example, European Bison (*Bison bonasus*) populations changed their ecology from a grazing diet in open landscapes toward a mixed-feeding behavior in more forested refugial habitats during the Holocene (11.7 ka cal. BP—present) as a response to human landscape transformations [6–9]. Such a change constitutes a perfect illustration of the complementarity between paleoecological and neoecological points of view when it comes to deciphering the ecological optimum and plasticity of a species [10, 11]. In that respect, studies relying only on the occurrence or frequency of macro-faunal species in fossil assemblages most likely miss a high proportion of paleoecological and paleoenvironmental information. To overcome these limitations, taxonomy-free analytical techniques taking the ecological plasticity into account are required. Among these, 3D Dental Microwear Texture Analyses (DMTA; [12–14]) constitute an efficient way to assess food type consumed and thus explore past ecologies, vegetal structure of local habitats and to infer paleoenvironments. This approach is based on the study of dental textures resulting from the ingested food items during the last days or weeks of the life of an animal [15–17].

To enhance zooarchaeological and paleoenvironnemental interpretations, we selected a key sequence from the European Middle Palaeolithic (Combe-Grenal, Dordogne, France) that recorded Neanderthal prey choices in the long term (over several dozens of millennia), and that has been at the center of research on past environments and subsistence strategies for decades [e. g. 18–26].

In Dordogne as in many other regions of Western Europe, this period is characterized by rapid and acute climatic oscillations [27, 28] that affected plant communities [29, 30] and, in turn, shaped micro-mammal and large mammal faunal communities [1, 31]. In the Combe-Grenal sequence (ca. 150,000 to 45,000 BP) as in other archeological sites, Neanderthals adapted their subsistence strategies to changes in animal populations. Reindeer (*Rangifer tarandus*), red deer (*Cervus elaphus*) and large bovids (steppe bison *Bison priscus* and/or aurochs *Bos primigenius*) constituted key resources for hunter-gatherer groups, providing both food (such as meat, marrow or grease) and raw materials (such as skin, bone, horn or antler). In this contribution, we present a long-term cross-species DMTA study of these three taxa
(reindeer, red deer and Bovinae) that are present throughout the Combe-Grenal sequence. Exploring the evolution of the feeding behavior of these herbivores is a first-choice approach to enlarge our understanding of the impact of climate change on human-ungulate-environment relationships during the Middle Palaeolithic.

With this study, we aim at fulfilling three objectives: (i) explore the paleoecology (notably the feeding preferences) of bovids and cervids hunted by Neanderthals at Combe-Grenal, (ii) contribute to a better understanding of Combe-Grenal's paleoenvironment and its evolution through time, and (iii) provide new insights on the subsistence strategies of Neanderthal populations, from the perspective of ungulate paleoecology. We test whether or not the evolution of the feeding ecology of preyed ungulate species through time reflect palaeoenvironmental changes throughout the sequence. This would allow us to better understand the impact of well-known global climate changes on the local paleoenvironments of Neanderthal hunting grounds.

# Material & methods

## Archeological site: Combe-Grenal

Combe-Grenal (Domme, southwestern France; Fig 1) is a rock-shelter situated in close proximity to the Dordogne River valley in the Périgord region. The 13m-long sequence was excavated during the 20th century by, successively, D. Peyrony, E. Peyrony and F. Bordes, before being the focus of a new field project since 2014 led by J.-Ph. Faivre and E. Discamps [22, 24, 25]. F. Bordes' excavations from 1953 to 1965 over 200 m² were the most extensive, and contributed in placing Combe-Grenal as one of the most important Middle Palaeolithic sequences of Western Europe. F. Bordes defined 65 layers in three superimposed terraces ("paliers") spanning MIS 6 to MIS 3 [22, 24, 25]. For this study, we focus on the uppermost terrace (layers 1 to 36) that yielded the most information on Middle Palaeolithic cultural and paleoenvironmental variability [26].

Bordes' collections are characterized by abundant faunal remains [18, 21], (12,000 specimens, belonging to 27 different species and representing more than 550 individuals) witnessing several clear faunal turnovers [1, 22]. At Combe-Grenal, 29 Neanderthal remains, some bearing cutmarks, were also discovered [24, 32], alongside more than 144,000 lithic artefacts relating to all of the main technological flake production systems recognized in the Mousterian of southwestern Europe (notably Levallois, Quina and Discoid; cf. [26] for a recent review). These rich industries contributed to the definition of the different "Mousterian facies" as defined by F. Bordes. The succession of technological systems defines a clear stratigraphic succession and put Combe-Grenal at the heart of the debate on Mousterian variability in Western Europe [26].

Recent fieldwork has provided key comparative data highlighting important recovery bias in previous excavations; proportions of the largest herbivores are highly over-estimated in at least some of Bordes' layers, as a consequence of preferential selection of teeth and larger bone fragments during excavation [22]. Bordes' collections are however still invaluable by their size and stratigraphic extent. Provided that the number of exploitable specimens is sufficient, it is worth mentioning that, unlike other approaches that depend on the relative abundance of remains, DMTA is unaffected by such recovery bias, and thus allows us to reliably establish a correlation between the feeding behavior of a population and the type of environment it occupies. For this study, we sampled Bordes' layers 4 to 36 (Fig 2).

That part of the sequence has not yet been adequately dated by radiometric methods. According to the most recent comparative chronological data [1, 26, 31] and pending the acquisition of new absolute dates, it has been tentatively correlated with MIS 5 or the beginning of MIS4 (layers 36 to 27, 115-65 ka cal. BP?), the end of MIS4 and Heinrich Stadial 6
<figure><img src="image_4.png" /><figcaption>**Fig 1.** **Combe-Grenal rock-shelter: Geographic, climatic and environmental context.** a) Localisation: red dot: Combe-Grenal; grey dots: key regional environmental proxies. Sources of the maps: https://srtm.csi.cgiar.org/; IGN. b) Environmental and climatic context of Combe-Grenal, based on independent environmental proxies [27–30]: i) δ¹³C isotopic records from Villars Cave stalagmites; ii) Pollen percentages of the boreal and altantic forests iii) Sea Surface Temperatures records from the marine core MD04-2845. Heinrich events 4–6 are indicated by light grey intervals. It is worth noting that, as the part of Combe-Grenal sequence studied here has not yet been adequately dated by radiometric methods, correlations between climatic events, environmental changes and the site stratigraphy remain uncertain.</figcaption></figure>

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

(layers 26 to 17, 65–60 ka cal. BP?) and the first half of MIS3 (layers 16 to 4, 60–43 ka cal. BP?). As illustrated in Fig 1, MIS 4 and 3 are characterized by rapid cycles of forest cover expansion and contraction in response to the millenial-scale climatic changes of the Dansgaard-Oeschger and Heinrich events, with a progressive trend towards more open environments [27–30].

Fig 2 summarizes contextual data available for layers 36 to 4. Material culture shows the succession of lithic industries dominated by Levallois, Quina and Discoid production systems [see 26 for a synthesis], as well as a change in mineral pigment use, with layers dominated by manganese oxides or iron oxides [33] (Fig 2A). Paleontological analysis of these layers by J.-L. Guadelli [20, 21] (layers 4–35) and G. Laquay [18], (layer 36) provided detailed information on the herbivores hunted by Neanderthals. Generally speaking, Combe-Grenal upper terrace sees a succession of faunas dominated by red deer (layers 36 to 34), Bovinae (layer 33), followed by a transition (layers 32 to 29) towards reindeer-dominated faunas (layers 28 to 17; Fig 2B). In layers 16 to 11, horse and large bovids are the most abundant, but new fieldwork highlighted the over-representation of large ungulates in these layers as a consequence of recovery bias [22]. Large bovids cannot always be identified to species, but a few specimens of both *Bison priscus* (n = 77) and *Bos primigenius* (n = 24) have been formerly identified in almost all layers [18, 21], (Fig 2B). When both species are present in the same layer, *Bison priscus* is more abundant in 10 out of 12 cases, but sample sizes are low (N≤17).
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Contextual data and distribution of DMTA samples for layers 36 to 4 of Combe-Grenal. a) Summary of material culture evidence for lithic industries [26] and mineral pigments [33], b) ungulate proportions (only layers with a total ungulate NISP [number of identified specimens] of more than 50 are figured), presence/absence of Bison priscus and Bos primigenius, and total ungulate NISP in italics [18, 21], c) number of DMTA samples included in the present study, both by "blocks" as defined in this study (cf. text) and for each of Bordes' layers.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0278395.g002

## Archeological material

For this study, we analyzed the teeth of 116 red deer (*Cervus elaphus*), 202 reindeer (*Rangifer tarandus*) and 50 Bovinae (Fig 2C). Both the steppe bison *Bison priscus* and the aurochs *Bos primigenius* are present at Combe-Grenal. The morphological proximity between the dental material of these two bovids makes specific determination difficult. Therefore, the majority of the bovid dental remains fall into the category "Bos/Bison" [21].
Specimens belonging to the four species were selected from the Middle Palaeolithic levels 4 to 36 of Combe-Grenal Bordes' collections, on the basis of the quality of the preservation of their dental facets. These four ungulate species, present in most of Combe-Grenal sequence, are characterized by a great abundance associated with a large paleobiogeographical distribution, the combination of both characteristics suggesting the use of diversified habitats and an important ecological plasticity. Bordes' collections are stored at the Musée National de la Préhistoire (MNP) of Les Eyzies, France. No permits were required for the described study, which complied with all relevant regulations.

## Extant species as comparative dataset

Seven extant ungulate populations with well-known diets were also included in the analyses to serve as a reference framework to interpret the diet of fossils (Table 1). The DMTA of all reference populations have been the subject of previous studies [5, 7, 34–36].

## DMTA implementation

Lower second molars were preferentially selected for the analysis. When this tooth was not available or too altered, the 1st or 3rd lower molar was selected. Details on the fossil sample are provided in S1 Appendix. Each tooth was carefully cleaned with a cotton swab soaked in ethanol. Once totally dry, the disto-lingual facet of the protoconid of lower molars was preferentially molded with a high-resolution polyvinylsiloxane elastomer (Regular Body President, ref. 6015–ISO 4823, medium consistency, polyvinylsiloxane addition type; Coltene Whaledent; see details on the procedures in the video here: http://anr-trident.prd.fr/v/; https://doi.org/10.5281/zenodo.7305567). The mold was then scanned directly with a Leica DCM8 white light confocal surface profilometer with a 100× objective (numerical aperture = 0.90, working distance = 0.9 mm) at the PAE2O “ArchéoScience” platform (UMR5608 TRACES, Toulouse, France). For each specimen we acquired a 251 × 333 µm surface (2584 × 1945 points, z-step: 0.2 µm) at the center of the dental facet. With LeicaMap 8.0.9, we then sampled a 200 × 200 µm sub-surface, leveled and mirrored in Z. The few missing points (<3%) were replaced by a smooth shape, using an algorithm calculated from the neighboring points (NMP ratio after cleaning: 0%). We applied the macro developed by Merceron et al. [37] to remove abnormal peaks. Remaining bigger artifacts and exogenous particles were removed manually and replaced by a smooth shape before proceeding to a final leveling of the surface. Photosimulations and false-color elevation maps of the processed data are available in S2 Appendix. All the data, both raw scans acquired with the Leica DCM8 (.plux files) and pre-treated data (.mnt files), are available on the NAKALA repository at the following URL: https://doi.org/10.34847/nkl.f4b0q9rm. Texture parameters were generated using a Scale Sensitive Fractal

<table><caption>**Table 1. Reference populations used in the present study, classified by their dietary category.** N = sample size. In the figures, these populations are referred with a population number, given in the last column.</caption><thead><tr><th>Dietary category</th><th>Species</th><th>Locality</th><th>N</th><th>References</th><th>n°</th></tr></thead><tbody><tr><td rowspan="3">Grazer</td><td><i>Bos taurus</i></td><td>Camargue Natural Regional Park, France</td><td>44</td><td>[35]</td><td>2</td></tr><tr><td><i>Cervus elaphus</i></td><td>Bauges Natural Regional Park, Bauges, France</td><td>23</td><td>[5, 36]</td><td>3</td></tr><tr><td><i>Rangifer tarandus</i></td><td>Knutshø, Forollhogna National Park, Norway</td><td>52</td><td>[5, 34]</td><td>7</td></tr><tr><td rowspan="3">Mixed-Feeder</td><td><i>Bison bonasus</i></td><td>Białowieża primeval Forest, Poland</td><td>19</td><td>[7]</td><td>1</td></tr><tr><td><i>Cervus elaphus</i></td><td>Domaniale forest of Chateauroux, France</td><td>29</td><td>[5]</td><td>5</td></tr><tr><td><i>Rangifer tarandus</i></td><td>Hardangervidda National Park, Norway</td><td>68</td><td>[5, 34]</td><td>6</td></tr><tr><td>Browser</td><td><i>Cervus elaphus</i></td><td>Białowieża primeval Forest, Poland</td><td>23</td><td>[5, 7]</td><td>4</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0278395.t001
<table><caption>**Table 2.** Values for DMTA-SSFA parameters for the three ungulate taxa from Combe-Grenal and for extant reference populations. N = sample size; s.e.m. = standard error of mean. **Periods:** Ante-Quina: levels 27 to 36, Quina: levels 17 to 26, Post-Quina: levels 4 to 16. **DMTA-SSFA parameters:** Asfc: complexity; epLsar: anisotropy; Smc: Scale of maximum complexity; HAsfc 9, 36, 81: heterogeneities of the complexity (3x3, 6x6, 9x9 cells).</caption><thead><tr><th rowspan="2">Fossil populations</th><th rowspan="2">Period</th><th rowspan="2">N</th><th colspan="2">Asfc</th><th colspan="2">epLsar (10<sup>-3</sup>)</th><th colspan="2">Smc</th><th colspan="2">HAsfc 9</th><th colspan="2">HAsfc 36</th><th colspan="2">HAsfc 81</th></tr><tr><th>mean</th><th>s.e.m.</th><th>mean</th><th>s.e.m.</th><th>mean</th><th>s.e.m.</th><th>mean</th><th>s.e.m.</th><th>mean</th><th>s.e.m.</th><th>mean</th><th>s.e.m.</th></tr></thead><tbody><tr><td rowspan="3">Bos/Bison (N = 50)</td><td>Ante-Quina</td><td>11</td><td>2.99</td><td>0.61</td><td>4.35</td><td>0.66</td><td>1.40</td><td>0.27</td><td>0.280</td><td>0.041</td><td>0.450</td><td>0.037</td><td>0.559</td><td>0.050</td></tr><tr><td>Quina</td><td>5</td><td>2.42</td><td>1.08</td><td>4.96</td><td>1.39</td><td>11.75</td><td>9.50</td><td>0.418</td><td>0.135</td><td>0.604</td><td>0.138</td><td>0.851</td><td>0.202</td></tr><tr><td>Post-Quina</td><td>34</td><td>2.11</td><td>0.14</td><td>5.19</td><td>0.43</td><td>1.70</td><td>0.18</td><td>0.304</td><td>0.027</td><td>0.464</td><td>0.025</td><td>0.602</td><td>0.032</td></tr><tr><td rowspan="3">Cervus elaphus (N = 116)</td><td>Ante-Quina</td><td>42</td><td>1.66</td><td>0.12</td><td>5.29</td><td>0.35</td><td>12.22</td><td>4.50</td><td>0.365</td><td>0.023</td><td>0.563</td><td>0.042</td><td>0.683</td><td>0.026</td></tr><tr><td>Quina</td><td>47</td><td>1.42</td><td>0.12</td><td>5.22</td><td>0.32</td><td>27.71</td><td>7.48</td><td>0.448</td><td>0.047</td><td>0.631</td><td>0.062</td><td>0.810</td><td>0.072</td></tr><tr><td>Post-Quina</td><td>27</td><td>1.87</td><td>0.19</td><td>5.77</td><td>0.45</td><td>20.73</td><td>11.37</td><td>0.353</td><td>0.023</td><td>0.505</td><td>0.024</td><td>0.669</td><td>0.031</td></tr><tr><td rowspan="3">Rangifer tarandus (N = 202)</td><td>Ante-Quina</td><td>57</td><td>1.83</td><td>0.11</td><td>5.03</td><td>0.23</td><td>4.54</td><td>1.47</td><td>0.362</td><td>0.023</td><td>0.553</td><td>0.037</td><td>0.685</td><td>0.045</td></tr><tr><td>Quina</td><td>94</td><td>2.09</td><td>0.12</td><td>4.21</td><td>0.19</td><td>6.31</td><td>2.61</td><td>0.373</td><td>0.018</td><td>0.538</td><td>0.024</td><td>0.684</td><td>0.026</td></tr><tr><td>Post-Quina</td><td>51</td><td>2.30</td><td>0.15</td><td>4.52</td><td>0.31</td><td>15.29</td><td>5.04</td><td>0.312</td><td>0.017</td><td>0.501</td><td>0.022</td><td>0.641</td><td>0.029</td></tr><tr><td colspan="2">Extant populations</td><td>N</td><td>mean</td><td>s.e.m.</td><td>mean</td><td>s.e.m.</td><td>mean</td><td>s.e.m.</td><td>mean</td><td>s.e.m.</td><td>mean</td><td>s.e.m.</td><td>mean</td><td>s.e.m.</td></tr><tr><td colspan="2">1: Bison bonasus (Białowieża Forest, Poland)</td><td>19</td><td>2.03</td><td>0.25</td><td>3.89</td><td>0.48</td><td>4.25</td><td>1.31</td><td>0.381</td><td>0.052</td><td>0.603</td><td>0.086</td><td>0.781</td><td>0.104</td></tr><tr><td colspan="2">2: Bos taurus (Camargue, France)</td><td>44</td><td>1.54</td><td>0.14</td><td>5.29</td><td>0.31</td><td>1.25</td><td>0.10</td><td>0.421</td><td>0.046</td><td>0.559</td><td>0.041</td><td>0.722</td><td>0.062</td></tr><tr><td colspan="2">3: Cervus elaphus (Bauges, France)</td><td>23</td><td>1.27</td><td>0.12</td><td>6.11</td><td>0.52</td><td>12.92</td><td>5.32</td><td>0.497</td><td>0.106</td><td>0.641</td><td>0.098</td><td>0.833</td><td>0.120</td></tr><tr><td colspan="2">4: Cervus elaphus (Białowieża Forest, Poland)</td><td>23</td><td>2.82</td><td>0.36</td><td>3.11</td><td>0.26</td><td>5.82</td><td>1.82</td><td>0.443</td><td>0.084</td><td>0.674</td><td>0.062</td><td>0.925</td><td>0.081</td></tr><tr><td colspan="2">5: Cervus elaphus (Chateauroux, France)</td><td>29</td><td>2.13</td><td>0.15</td><td>4.29</td><td>0.37</td><td>8.29</td><td>2.59</td><td>0.452</td><td>0.045</td><td>0.651</td><td>0.061</td><td>0.873</td><td>0.082</td></tr><tr><td colspan="2">6: Rangifer tarandus (Hardangervidda, Norway)</td><td>68</td><td>2.30</td><td>1.08</td><td>3.15</td><td>1.54</td><td>9.00</td><td>21.23</td><td>0.460</td><td>0.294</td><td>0.769</td><td>0.342</td><td>0.822</td><td>0.398</td></tr><tr><td colspan="2">7: Rangifer tarandus (Knutsho, Norway)</td><td>52</td><td>1.76</td><td>0.64</td><td>4.21</td><td>2.03</td><td>17.93</td><td>35.38</td><td>0.410</td><td>0.206</td><td>0.698</td><td>0.300</td><td>0.714</td><td>0.311</td></tr></tbody></table>

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

Analysis (SSFA) with the same software. In this study we considered area-scale fractal complexity (Asfc; Tables 2 and 3, S1 Appendix), exact proportion of the length-scale anisotropy of the relief (epLsar; Tables 2 and 3, S1 Appendix), Scale of maximum complexity (Smc; Tables 2 and 3, S1 Appendix) and heterogeneity of the area-scale fractal complexity at 9, 36 and 81 cells (HAsfc 9, HAsfc 36, HAsfc 81; Tables 2 and 3, S1 Appendix). These parameters are further described in Scott [13] and Scott et al. [12].

## Statistics

We tested both inter- and intra-specific dietary variations through time. Low sample sizes for the three studied ungulates in most archeological layers of Combe-Grenal prevented us from statistically testing for texture variations between them (Fig 2C). In order to counterbalance this issue, the choice has been made to group together Bordes' layers into archeological “Blocks” that are overall similar in terms of faunal associations and lithic industry prior to statistical analyses [1, 22, 26], (Fig 2). Therefore, block A corresponds to levels 4 to 10 of Bordes (Discoid/Levallois, mixed faunas), Block B to levels 11 to 15 (Discoid, dominated by large ungulates?), Block C to level 16 (no DMTA sample), Block D to levels 17 to 19 (Quina, mostly reindeer), Block E to level 20 (Quina rich in denticulates, mostly reindeer), Block F to levels 21 to 26 (Quina, mostly reindeer), Block G to 27 to 32 (Levallois, mixed faunas), Block H to level 33 (Levallois, large bovids?), Block I to levels 34 to 36 (Levallois, red deer). We performed an alternative analysis exploring statistical differences between three larger time periods defined in previous works [31, 38]: “Ante-Quina” [levels 27–36], Quina [levels 17–26] and “Post-Quina” [levels 4–16].

Statistical analyses were done following Smith & Warren [39] and Zuur et al. [40]. For each ungulate species, we tested the variation of its dental textures through time (Block or Period).
<table><caption>**Table 3.** Values for DMTA-SSFA parameters in each block for the three ungulate taxa from Combe-Grenal. N = sample size; s.e.m. = standard error of mean. Block A: levels 4 to 10; Block B: levels 11 to 15; Block C: level 16 (no DMTA sample), Block D: levels 17 to 19; Block E: level 20; Block F: levels 21 to 26; Block G: levels 27 to 32; Block H: level 33; Block I: levels 34 to 36. **DMTA-SSFA parameters:** Asfc: complexity; epLsar: anisotropy; Smc: Scale of maximum complexity; HAsfc 9, 36, 81: heterogeneities of the complexity (3x3, 6x6, 9x9 cells).</caption><thead><tr><th rowspan="2">Fossils</th><th rowspan="2">Block</th><th rowspan="2">N</th><th colspan="2">Asfc</th><th colspan="2">epLsar (10-3)</th><th colspan="2">Smc</th><th colspan="2">HAsfc 9</th><th colspan="2">HAsfc 36</th><th colspan="2">HAsfc 81</th></tr><tr><th>mean</th><th>s.e.m.</th><th>mean</th><th>s.e.m.</th><th>mean</th><th>s.e.m.</th><th>mean</th><th>s.e.m.</th><th>mean</th><th>s.e.m.</th><th>mean</th><th>s.e.m.</th></tr></thead><tbody><tr><td rowspan="10">Bos/Bison</td><td>ALL</td><td>50</td><td>2.34</td><td>0.19</td><td>4.98</td><td>0.35</td><td>2.64</td><td>0.97</td><td>0.310</td><td>0.024</td><td>0.475</td><td>0.023</td><td>0.617</td><td>0.032</td></tr><tr><td>A</td><td>9</td><td>2.19</td><td>0.34</td><td>4.46</td><td>0.93</td><td>1.67</td><td>0.16</td><td>0.380</td><td>0.066</td><td>0.516</td><td>0.061</td><td>0.656</td><td>0.077</td></tr><tr><td>B</td><td>25</td><td>2.08</td><td>0.14</td><td>5.45</td><td>0.48</td><td>1.71</td><td>0.24</td><td>0.277</td><td>0.026</td><td>0.445</td><td>0.027</td><td>0.582</td><td>0.033</td></tr><tr><td>C</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>E</td><td>1</td><td>0.92</td><td></td><td>8.43</td><td></td><td>3.35</td><td></td><td>0.162</td><td></td><td>0.384</td><td></td><td>0.508</td><td></td></tr><tr><td>F</td><td>4</td><td>2.80</td><td>1.31</td><td>4.09</td><td>1.41</td><td>13.85</td><td>11.96</td><td>0.482</td><td>0.153</td><td>0.659</td><td>0.163</td><td>0.936</td><td>0.235</td></tr><tr><td>G</td><td>10</td><td>3.08</td><td>0.67</td><td>4.17</td><td>0.69</td><td>1.40</td><td>0.30</td><td>0.268</td><td>0.043</td><td>0.449</td><td>0.041</td><td>0.559</td><td>0.055</td></tr><tr><td>H</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>I</td><td>1</td><td>2.1183</td><td></td><td>6.2216</td><td></td><td>1.361</td><td></td><td>0.3965</td><td></td><td>0.4688</td><td></td><td>0.5627</td><td></td></tr><tr><td rowspan="10">C. elaphus</td><td>ALL</td><td>116</td><td>1.61</td><td>0.08</td><td>5.37</td><td>0.21</td><td>20.48</td><td>4.35</td><td>0.396</td><td>0.022</td><td>0.577</td><td>0.030</td><td>0.731</td><td>0.032</td></tr><tr><td>A</td><td>15</td><td>1.72</td><td>0.25</td><td>6.29</td><td>0.59</td><td>36.18</td><td>19.85</td><td>0.347</td><td>0.028</td><td>0.486</td><td>0.027</td><td>0.680</td><td>0.040</td></tr><tr><td>B</td><td>12</td><td>2.06</td><td>0.28</td><td>5.13</td><td>0.69</td><td>1.42</td><td>0.17</td><td>0.362</td><td>0.039</td><td>0.528</td><td>0.044</td><td>0.655</td><td>0.051</td></tr><tr><td>C</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>D</td><td>5</td><td>2.29</td><td>0.49</td><td>4.12</td><td>0.82</td><td>1.10</td><td>0.25</td><td>0.438</td><td>0.074</td><td>0.548</td><td>0.058</td><td>0.662</td><td>0.063</td></tr><tr><td>E</td><td>5</td><td>1.84</td><td>0.45</td><td>4.64</td><td>1.26</td><td>11.00</td><td>9.68</td><td>0.294</td><td>0.093</td><td>0.449</td><td>0.071</td><td>0.559</td><td>0.079</td></tr><tr><td>F</td><td>37</td><td>1.25</td><td>0.11</td><td>5.44</td><td>0.36</td><td>33.56</td><td>9.20</td><td>0.470</td><td>0.057</td><td>0.667</td><td>0.077</td><td>0.864</td><td>0.089</td></tr><tr><td>G</td><td>19</td><td>1.73</td><td>0.19</td><td>4.65</td><td>0.57</td><td>16.88</td><td>8.72</td><td>0.393</td><td>0.046</td><td>0.622</td><td>0.088</td><td>0.711</td><td>0.049</td></tr><tr><td>H</td><td>2</td><td>1.99</td><td>0.74</td><td>8.45</td><td>0.46</td><td>1.18</td><td>0.18</td><td>0.439</td><td>0.021</td><td>0.656</td><td>0.026</td><td>0.707</td><td>0.002</td></tr><tr><td>I</td><td>21</td><td>1.56</td><td>0.15</td><td>5.57</td><td>0.39</td><td>9.06</td><td>4.39</td><td>0.333</td><td>0.018</td><td>0.501</td><td>0.021</td><td>0.655</td><td>0.028</td></tr><tr><td rowspan="10">R. tarandus</td><td>ALL</td><td>202</td><td>2.07</td><td>0.07</td><td>4.52</td><td>0.14</td><td>8.08</td><td>1.82</td><td>0.354</td><td>0.012</td><td>0.533</td><td>0.016</td><td>0.673</td><td>0.019</td></tr><tr><td>A</td><td>15</td><td>2.45</td><td>0.38</td><td>5.31</td><td>0.64</td><td>32.61</td><td>12.09</td><td>0.316</td><td>0.040</td><td>0.522</td><td>0.049</td><td>0.679</td><td>0.059</td></tr><tr><td>B</td><td>36</td><td>2.24</td><td>0.15</td><td>4.19</td><td>0.33</td><td>8.07</td><td>4.68</td><td>0.310</td><td>0.018</td><td>0.492</td><td>0.024</td><td>0.625</td><td>0.032</td></tr><tr><td>C</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>D</td><td>9</td><td>2.73</td><td>0.61</td><td>2.77</td><td>0.38</td><td>1.92</td><td>0.53</td><td>0.271</td><td>0.038</td><td>0.448</td><td>0.040</td><td>0.582</td><td>0.044</td></tr><tr><td>E</td><td>9</td><td>2.43</td><td>0.29</td><td>3.82</td><td>0.50</td><td>3.02</td><td>1.52</td><td>0.435</td><td>0.053</td><td>0.585</td><td>0.041</td><td>0.736</td><td>0.051</td></tr><tr><td>F</td><td>76</td><td>1.97</td><td>0.12</td><td>4.43</td><td>0.21</td><td>7.22</td><td>3.22</td><td>0.377</td><td>0.021</td><td>0.543</td><td>0.029</td><td>0.690</td><td>0.031</td></tr><tr><td>G</td><td>43</td><td>1.83</td><td>0.13</td><td>4.87</td><td>0.26</td><td>5.36</td><td>1.94</td><td>0.377</td><td>0.029</td><td>0.560</td><td>0.046</td><td>0.692</td><td>0.055</td></tr><tr><td>H</td><td>5</td><td>1.56</td><td>0.29</td><td>5.88</td><td>0.65</td><td>1.43</td><td>0.33</td><td>0.341</td><td>0.069</td><td>0.594</td><td>0.109</td><td>0.761</td><td>0.137</td></tr><tr><td>I</td><td>9</td><td>2.00</td><td>0.30</td><td>5.35</td><td>0.60</td><td>2.37</td><td>0.56</td><td>0.303</td><td>0.041</td><td>0.497</td><td>0.060</td><td>0.611</td><td>0.081</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0278395.t003

Among each time frame considered (Block or Period), we also tested for inter-specific differences in dental textures. Data exploration and statistics were conducted with R v. 4.0.3 and the packages “outliers”[41], “ggplot2”[42], “lawstat”[43], “onewaytests”[44], “MASS”[45], “Hmisc”[46], “corrplot”[47], “regplot”[48], “AICcmodavg”[49], “emmeans”[50], “stats”[51], “lattice”[52], “FSA”[53]. The database and R markdown script are respectively available in S1 and S3 Appendices. After a visual data exploration (notably frequency plots and box-plots), [40] and Brown and Forsythe tests for homogeneity of variance [39] to test the applicability of parametric tests, we performed BoxCox data-transformations whenever needed (among reindeer, for inter-block differences in Smc; among red deer, for inter-block differences in HAsfc 81; for block H, inter-specific differences in epLsar; S3 Appendix). Generalized Linear Models (GLM) were then performed with the objective to select, based on the lower AIC, the model that best explains the distribution of data (“Block”, “Period” or “random”; “species” or “random”; AICs are available in S3 Appendix). Whenever best models showed
statistical differences, post hoc Tukey HSD pairwise tests allowed us to identify inter-group
significant differences. In the situation where the conditions for applying parametric tests were
not satisfied (for Smc when testing for inter-specific differences during the Quina and inter-
specific differences among the block F; and inter-period and inter-block differences in epLsar
for reindeer, S3 Appendix), a Kruskal-Wallis test followed by a Dunn's post-hoc test were used
to explore the inter-group differences. We did not discuss significant differences between
groups represented by less than 5 individuals.

**Results**

**1-Main represented dietary categories**

The three populations of extant ruminants assigned as grazers (lower part of Table 2, Fig 3; populations 2, 3, 7) are characterized by low complexities (Asfcs) and heterogeneities of the complexity (HAsfcs 9) coupled with high anisotropies (epLsar). Low to intermediate anisotropies (epLsar) coupled with middle to high complexities (Asfcs) and intermediate to high heterogeneities of the complexity (HAsfcs 9) are typical of the browsing dietary category, represented here by the red deer population from Białowieża primeval forest (Table 2; population 4). The mixed-feeder are represented here by three populations including the red deer sample from Chateauroux (population 5), the population of Bison bonasus from Białowieża (population 1) and the Rangifer tarandus population from Hardangervidda (population 6). They are characterized by intermediate values of their texture parameters (Fig 3; lower part of Table 2).

Results for fossils are discussed based on the herbivore dietary ecospace defined here on the
basis of the seven extant European ungulate populations with well-known diet (Tables 1 and 2,
Fig 3; see also [5]). Bos/Bison and Rangifer tarandus are characterized by medium to high
mean Asfcs and epLsar values. Heterogeneity (HAsfcs 9) is generally low. For Cervus elaphus, the
Asfcs varies from low to intermediate values, coupled with medium to high epLsar and low
HAsfcs 9 in most cases.

**2-Inter-specific differences between contemporaneous ungulate populations**

During both Ante-Quina (*Asfcs*, *Smc*), Quina (*Asfcs*, *epLsar*, *Smc*, *HAsfcs 81*) and Post-Quina (*epLsar*, *Smc*), the models based on the “species” variable explain better the distribution than the “random” models (lower AIC; S3 Appendix). Post-hoc tests allowed us to identify significant differences (p-value <0.05) between bovids and cervids during the Ante-Quina (bovids have higher *Asfcs* during this period; Table 2; Fig 4, S1 and S3 Appendices) and between red deer and reindeer during the Quina (*Asfcs* is lower while *epLsar* and *Smc* are higher for red deer during this period; Table 2; Fig 4, S1 and S3 Appendices).

In Blocks B (*epLsar*), D (*HAsfcs* 9), E (*HAsfcs* 36 and *HAsfcs* 81), F (*Asfcs*, *epLsar*, *Smc*, *HAsfcs* 81) and G (*Asfcs* and *Smc*) the variable “species” is more efficient than random to explain the distribution (smaller AIC, S3 Appendix). As identified with post-hoc tests, interspecific differences are located between red deer and reindeer in Block D (*HAsfcs* 9 is higher for red deer; Table 3) and Block F (*epLsar* and *Smc* are higher and *Asfcs* lower for red deer; N<5 for Bovidae); Table 3; Fig 5, S1 and S3 Appendices). In Block G, *Cervus elaphus* differs in *Asfcs* from the two other taxa (Table 3; Fig 5, S1 and S3 Appendices).

**3-Intra-specific dietary variations through time**

For almost all texture parameters and for the three taxa, as tested with the GLM models, data
better fits a random distribution than a distribution explained by Blocks (S3 Appendix). Based
<figure><img src="image_6.png" /><figcaption>Fig 3. Distribution of the three ungulate taxa from Combe-Grenal (all layers grouped) based on the complexity (Asfc), anisotropy (epLsar) and heterogeneities of the complexity (HAsfc 9 and HAsfc 81) of the dental textures.</figcaption></figure>

Each taxon and reference population is represented by its mean and standard error of mean. Green circle: Cervus elaphus; Blue square: Rangifer tarandus; Red triangle: Bison priscus/Bos primigenius. Reference populations with **well-known diet**: 1: Bison bonasus, Białowieża, Poland; 2: Bos taurus, Camargue, France; 3: Cervus elaphus, Bauges, France; 4: Cervus elaphus, Białowieża, Poland; 5: Cervus elaphus, Chateauroux, France; 6: Rangifer tarandus, Hardangervidda, Norway; 7: Rangifer tarandus, Knutshø, Norway.

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

on GLM (S3 Appendix), the periods better explain the distribution for *Asfc* (*Cervus elaphus*, *Rangifer tarandus*), *epLsar* (*Rangifer tarandus*), *Smc* (*Bos/Bison*), *HAsfc 9* (*Rangifer tarandus*), *HAsfc 36* and *HAsfc 81* (*Bos/Bison*). However, only few of these differences are supported by post-hoc tests. Bovids present higher HAsfc 81 values during the Quina (Table 2, Fig 4, S1 and S3 Appendices). The Kruskal-Wallis analysis followed by a Dunn post-hoc test testing for inter-Period differences for the *epLsar* of *Rangifer tarandus* support a significant difference
Reference populations

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

Ante-Quina (blocks I to G)

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

Quina (blocks F to D)

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

Post-Quina (blocks C to A)

<figure><img src="image_10.png" /></figure>
**Fig 4.** Distribution (mean and standard error of mean) for the three taxa for each period considered in this study, depending on complexity (Asfc) and anisotropy (epLsar (x10⁻³); left) and on heterogeneities of the complexity (HAsfc 9 and HAsfc 81; right). Reference populations are illustrated on the top of the figure and used as reference for the Ante-Quina, Quina and Post-Quina plots. **They are identified as follows:** 1: Bison bonasus, Bialowieża, Poland; 2: Bos taurus, Camargue, France; 3: Cervus elaphus, Bauges, France; 4: Cervus elaphus, Bialowieża, Poland; 5: Cervus elaphus, Chateauroux, France; 6: Rangifer tarandus, Hardangervidda, Norway; 7: Rangifer tarandus, Knutshø, Norway. Bos/Bison is represented by red triangles, Cervus elaphus by green circles and Rangifer tarandus by blue squares.
**Periods:** Ante-Quina: levels 27 to 36, Quina: levels 17 to 26, Post-Quina: levels 4 to 16.

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

between Ante-Quina and Quina (Table 2; Fig 4, S1 and S3 Appendices). Although not significant, we also observed a tendency toward an increase of *Asfc* and a decrease of *epLsar* through time for both cervids from Block H to Block D.

Discussion

Ungulates from Combe-Grenal: Paleoecology

When compared with the values for the seven reference populations categorizing browsing to grazing dietary categories (Table 1; Fig 3), the dental textures of the 116 *Cervus elaphus* from Combe-Grenal reflect a grazing diet, a feeding behavior that is well-known for extant relatives of this species living in open habitats [54, 55]. The extant red deer is indeed a polymorphic species. Its body mass variations range from one to five from one population to another, a particularity that has to be linked with the diversity of the habitats that are occupied by the species [56]. Through Eurasia, the red deer is indeed currently present from Southern Spain to Northern Norway and from the Atlantic coast of Western Europe to the Caspian sea in a large variety of habitats [57]. As shown by Gebert and Verheyden-Tixier [54] the main driver of its interpopulation dietary variations is the characteristics of its habitat, which makes of the feeding ecology of this cervid a very appropriate proxy for the paleoenvironmental reconstructions in archeological sites [5, 58].

The wild reindeer is nowadays probably the most dietary plastic of European deer, its ecology and body size being dependent on the characteristics of its habitat [56, 59, 60]. Its current distribution is circumpolar [56, 61]. It can tolerate extreme cold and its morphology facilitates its movements in snow [62, 63]. This deer can be found in both mountainous environments and oceanic coasts, in the tundra as in the taiga. The reindeer consumes what is available in its environment according to the seasonal plant phenology [56, 59, 64, 65]. As it is able to consume a resource until it is depleted, it also plays an important role in shaping the vegetation and landscapes [66–68]. When these resources are available, reindeer favors food with a high nutritional value. However, this species is very opportunistic, giving to its feeding behavior a potential interest for paleoenvironmental reconstructions. Knowing this, the mixed-feeding diet of the 202 *Rangifer tarandus* from Combe-Grenal witnesses the availability of a large amount of herbaceous monocotyledons, associated with a significant part of ligneous material, in the vicinity of Combe-Grenal.

The steppe bison *Bison priscus* is a large bovid which taxonomic relationships remain under debate [69–71 and references therein]. This species was one major element of the “mammoth steppe” mammalian assemblage [72, 73] and one of the key prey for European human populations in numerous fossil localities throughout Eurasia [74]. Its vast distribution, from the early middle Pleistocene to the early Holocene (that is to say between 0.781 and 0,0117 Ma), covered Eurasia and North America [72, 75]. During the last glaciation at the Pleistocene-Holocene transition, the species abundance drastically decreased to the point of its extinction. Climate and anthropogenic impacts are among the factors considered responsible for its extinction [75, 76]. The diversity of habitats (from cold and opened steppes to more forested habitats), of possible feeding behaviors [including both herbaceous monocotyledons and lichens: 77, 78]
<figure><img src="image_11.png" /><figcaption>A. Complexity (Asfc)</figcaption></figure>

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

Reindeer

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

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

Red deer

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

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

Large bovids

4.0

Browsers

2.0

Grazers

0.0

<caption>B. Anisotropy (epLsar)</caption>

I

H

G

F

E

D

C

B

A

10.0

Grazers

6.0

4.0

Browsers

0.0

<caption>C. Heterogeneity 9 cells (HAsfc 9)</caption>

I

H

G

F

E

D

C

B

A

0.6

0.4

0.2

0.0

1

2

3

4

5

6

7

Reference populations

I

H

G

F

E

D

C

B

A

Post-Quina

Ante-Quina

Quina
**Fig 5.** Barplots representing the mean and standard error of mean for Blocks A to I of Combe-Grenal and for extant European ungulate reference populations with well-known diet. **Reference populations:** 1: Bison bonasus, Białowieża, Poland; 2: Bos taurus, Camargue, France; 3: Cervus elaphus, Bauges, France; 4: Cervus elaphus, Białowieża, Poland; 5: Cervus elaphus, Chateauroux, France; 6: Rangifer tarandus, Hardangervidda, Norway; 7: Rangifer tarandus, Knutshø, Norway. **Temporal Blocks:** Block A: levels 4 to 10, Block B: levels 11 to 15, Block C: level 16, Block D: levels 17 to 19, Block E: level 20, Block F: levels 21 to 26, Block G: 27 to 32, Block H: level 33, Block I: levels 34 to 36. **Periods:** Ante-Quina: levels 27 to 36, Quina: levels 17 to 26, Post-Quina: levels 4 to 16. Significant differences (see also S3 Appendix) are illustrated by a red "*", for groups with N≥5 individuals.

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

and of mobility behaviors of *Bison priscus*, capable of both migration and sedentarity [74], reflects the important ecological plasticity of the species. The aurochs (*Bos primigenius*) is generally considered as the wild ancestor of *Bos taurus*, the domestic cattle [79 and references therein]. The first occurrences of aurochs date back to the Pleistocene [between 2–1.5 Ma; 81]. Its geographical range included northern Africa and almost the whole Eurasia [80]. Overhunting and suitable habitat reduction have undoubtedly played a major role in their drop until their extinction, that most likely occurred in Poland during the 17th century [80]. While several studies highlight the preference of *Bos primigenius* for forested habitats under mild and humid climatic conditions [80–82 and references therein], the feeding ecology of this bovid suggests a degree of dietary plasticity that allows it to adapt its diet over time in response to environmental changes [83]. The 50 Bovinae from Combe-Grenal are mixed feeders. Here again, such results support the occurrence of an abundant herbaceous monocotyledon layer, coupled with ligneous material.

Inside each of the temporal Blocks, significant inter-specific dietary differences (considered for N≥5 individuals only; Fig 5) are scarce, except during the Ante-Quina (Block G) between *Cervus elaphus* and the other taxa and the Quina (blocks D and F) when *Cervus elaphus* and *Rangifer tarandus* diets statistically differ. Before the Quina, both cervids are more engaged in grazing while bovids present dental textures reflecting the ingestion of a large amount of ligneous material, similar to the reference population of *Cervus elaphus* from Białowieża primeval forest (Fig 5). Block F reflects the difference between the grazing diet of *Cervus elaphus* and the mixed-feeding diet of *Rangifer tarandus*. The heterogeneity of the complexity (*HA*sfc 9) is generally interpreted as reflecting the diversity of the ingested food items [12, 13]. In Block D, *HA*sfc 9 is lower for *Rangifer tarandus*.

Similar dental textures in most cases for the four contemporaneous and sympatric ungulate populations indicate no obvious niche-partitioning between these large herbivores.

## Comparison with previous paleoenvironmental interpretations

Considering the well-known, previously discussed, ecological plasticity of the four ungulates investigated, able to adapt their feeding ecology to changes in resource availability in their habitat, we could expect DMTA analysis of these taxa to track environmental changes known for the region from the MIS5 to the MIS3 (Fig 1). However, in contradiction to this hypothetical proposition, the dental textures do not reflect clear variations in the feeding ecology of the four ungulates from Combe-Grenal through time. This archeological site has however frequently been cited as an example of a sequence that document marked environmental changes during the Middle Palaeolithic. Previously acquired paleoenvironmental data can be summarized as follows:

i. Clear changes in the climate and vegetation have been documented in Western France from MIS 5 to MIS 3 [27–30], (Fig 1), but the environmental conditions in the immediate vicinity of Combe-Grenal are largely unknown.

ii. Changes in the proportions of hunted ungulates at Combe-Grenal [20, 21] have however been shown to track regional vegetation changes known for the period [1, 31].
iii. Data on rodent species are scarce, but highlight a slight diachronic pattern, notably the presence of the woody *Microtus* (*Terricola*) *subterraneus* during the “Ante-Quina” (including in layer 31) but not in the following Quina and “Post-Quina” periods [1, 84].

iv. Isotopic analyses of *Equus* teeth from layers 30 to 4, that are contemporaneous to the ungulates analysed in the present study, were interpreted by the authors as evidencing marked changes in temperature and aridity throughout the sequence [23]. In particular, they highlighted values in dentine $\delta^{13}$C during our period of interest that are in favour of an open environment, and ample variations in both dentine $\delta^{15}$N and enamel phosphate $\delta^{18}$O through time [23]. For several layers of interest to our study [4, 9, 16, 20, 21, 32, 33, 35] isotope results were however unavailable in the study of Richards et al. [23].

We observe for the two deer a global tendency toward a decrease in the ingestion of abrasive herbaceous monocotyledons from blocks H to D (layers 33 to 17). This decrease is supported by both the progressive increase of complexity (Asfc; Fig 5A; Table 3) and a decrease of anisotropy (epLsar; Fig 5B; Table 3). For reindeer, this decrease in anisotropy is supported at a larger time scale by the significant difference between Ante-Quina and Quina (p-value <0.05; Fig 4). The fact that this trend occurs for both cervids reinforces its relevance, despite being only “partially significant” (not significant at the thinner scale of blocks, but significant between Ante-Quina and Quina for reindeer). Despite previously mentioned paleoclimatic variations, paleoenvironmental changes and faunal turnovers (Figs 1 and 2), our DMTA results imply that the same types of resources were available and consumed by the four preyed ungulates throughout the sequence, at least during the last days or weeks prior to their death.

These seemingly contradictory findings are the result of a combination of the following two factors: (i) from a technical perspective, the apparent discrepancy between paleoenvironmental data acquired by DMTA analysis and from other proxies including isotopic analyses reflect the different periods of life and time scales that these approaches document [85]. While most other analyses document the feeding ecology, climatic and environmental conditions experienced by animal populations during few months to few years, DMTA on the other hand is an open window over the feeding preferences and resource availability for the animals shortly before their death ([15], the last few days or weeks, also referred as the « last supper effect »; [16, 17]). (ii) from an archaeological perspective, DMTA results might be explained, for a large part, by the impact of past hunting practices by Neanderthal groups on the sampled ungulate individuals. This latter aspect is further developed in the next section.

## Neanderthal subsistence strategies and habitat uses

The DMTA of ecologically plastic herbivores constitutes an efficient proxy for paleoenvironmental reconstructions [5, 58]. Based on this assumption and on Combe-Grenal ungulate textures (Fig 3), it is possible to infer that *Cervus elaphus*, *Rangifer tarandus* and *Bos/Bison* from Combe-Grenal were, at the time when they were killed, occupying an open, tundra-like habitat, characterized by the predominance of an herbaceous layer associated with shrubs that explain the grazing diet of *Cervus elaphus* and the mixed-feeding diet of both bovids and reindeer. No major change in the resource availability in the habitat(s) where these animals have been hunted, as witnessed by the absence of significant dental texture variations through time for the three ungulate taxa (Fig 5), is in favor of a certain habitat and dietary homogeneity throughout millennia. As these dental textures belong to animals that were hunted by successive Neanderthal groups that occupied Combe-Grenal rock-shelter, they likely represent a non-random sample of the wild population(s) from which they were selected, potentially biased in the direction of particular hunting seasons or preferential hunting territories. By
selectively hunting in particular landscapes of their territory, hunters can thus be seen as a
“buffering” agent of the paleoenvironmental signal. Tracking and killing ungulates in open set-
tings was probably easier for Neanderthal hunters, and this preference for open settings poten-
tially made environmental changes invisible in a DMTA analysis of archaeofaunas. Under
these conditions, DMTA provide us with a very important information: the continuous occur-
rence, throughout the sequence of Combe-Grenal, of grassland areas that were sufficiently
extensive to constitute a hunting ground for Neanderthal hunters and a major part of the diet
of the hunted animals shortly before their death. However, these results do not document
changes, expected as a consequence of climatic changes, in the location and/or variation in the
extent of these grassland areas in the hunting territories exploited by Neanderthal groups over
time. This potential “buffering” effect is undoubtedly an element to be taken into account
when aiming to infer the palaeoenvironmental changes occurring throughout an archeological
sequence using DMTA.

The regional climatic changes through time, notably the progressive habitat opening and
rapid climatic changes of MIS 4–3 induced rapid cycles of forest cover expansion and contrac-
tion, that impacted the food resources that were available for large ungulates [27–30], (Fig 1).
The modification towards peri-optimal habitats might have forced ecological adaptations for
several ungulate species, then followed by modifications of their range of distribution, that are
reflected in the archeological record by changes in faunal associations. In this context, Nean-
derthal hunting strategies seem to involve the use of open landscapes over the entire sequence,
despite clear paleo-environmental changes in the region. This is how we interpret the homoge-
neous DMTA results that we obtain, despite varying proportions of large ungulates across the
sequence. It is the combination of ungulate ecological plasticity and hunting selection by
Neanderthal groups that produces such an apparent “discrepancy” between proxies (e.g. pro-
portions of hunted ungulates on one hand, DMTA data on the other). MIS 5 to 3 paleoenvir-
onmental changes (such as decreasing tree cover and temperature, cf. Fig 1) might have
induced demographic changes in herbivore populations, yet some of these ungulates persisted
in open grasslands and were preferentially targeted by Neanderthal hunters.

Alternatively, the possibility of a specific hunting of ungulates at a given season during
which a large proportion of herbaceous monocotyledons was consumed also remains to be
explored. Data on season of death for cervids and bovids in layers 4–36 of Combe-Grenal is
extremely scarce: an analysis of tooth eruption sequences and wear allowed Guadelli [21;
p. 453–457] to propose reindeer hunting throughout the year in layers 22 and 23. Steele [86]
analysis of red deer teeth did not provide detailed information on layers 1 to 36 due to small
sample size. Finally, Binford [19] data is not available by layer and insufficiently detailed, mak-
ing its integration difficult. Acquiring further seasonal data would contribute to our under-
standing of Neanderthal hunting strategies by understanding whether these animals were
culled by Neanderthal populations during seasonal hunts or over several seasons. This aspect
needs supplementary investigation and will benefit from the comparison of our DMTA results
with those of cementum analyses [see for example 87], currently underway in the framework
of the ANR DeerPal project. Acquiring DMTA data from lower parts of the sequence (below
layer 36) would also be of interest to test if hunted *Cervus elaphus* always presented a grazing
diet at Combe-Grenal.

**Conclusion**

In this study, DMTA analysis of preyed *Cervus elaphus*, *Rangifer tarandus* and *Bos/Bison* from
Combe-Grenal, that reflects ungulate dietary preferences shortly before their death, support
the fact that these animals were, at the time when they were killed, occupying an open, tundra-
like habitat. Despite evidence for climatic and environmental changes in the vicinity of the site from MIS 5 to MIS 3, dental textures show relatively little variation through time. This apparent discrepancy between environmental changes and the feeding ecology of these dietary plastic ungulates highlights the continuous preference of Combe-Grenal Neanderthals for hunting in open landscapes, and the persistence of grassland areas in Combe-Grenal vicinity that were sufficiently extensive to constitute the principal hunting ground for Neanderthal hunters all along the sequence. Exploring the feeding ecology of preyed ungulates from other contemporaneous archaeological sites will allow us to decipher if these findings can be extrapolated at the more regional scale of southwestern France or if they reflect specific hunting strategies in the territories used by Neanderthal populations in the vicinity of Combe-Grenal.

DMTA permits to go further in our understanding of human-environment interactions: even when there are ample evidences for paleoenvironmental changes, zooarchaeologists still have to test whether or not these changes drastically altered hunting ground use or not. In the case of Combe-Grenal, it might be the case that they did not. The three taxa studied here, in addition to horses, formed most of Neanderthal subsistence for all of Combe-Grenal layers 36 to 4. DMTA analysis identifies only “open-living” individuals. This means that Combe-Grenal hunters “stayed in the open”, and were never forced to switch to hunting tactics adapted to close encounters in forested environments. As such, they never had, for example, to alter their lithic implements or mobility strategies for such use. Any discussion on the influence of paleoenvironmental change on material culture or human history should include such a precise understanding of local environmental conditions if misinterpretations are to be avoided.

**Supporting information**

**S1 Appendix. Individual information and values for texture parameters for all Combe-Grenal specimens.** Ref: reference for each specimen; c: layer; lm: lower molar; dex: right; sin: left; bd: bucco-distal facet; bm: bucco-mesial facet. (XLSX)

**S2 Appendix. Photosimulations and false color elevation maps of molar facets of Combe-Grenal Ungulates, scanned by E. Berlioz and E. Capdepon with the white light confocal microscope Leica DCM8 of TRACES Lab, Toulouse.**
(PDF)

**S3 Appendix. R markdown script including results, used to explore and analyze the dataset of the present study.**
(PDF)

**Acknowledgments**

Authors are grateful to the Musée national de Préhistoire (Les Eyzies), notably Stéphane Madeleine, Peggy Jacquement and Catherine Cretin, for their help and for giving access to the material. Comments by Ivan Calandra and Gregory James Smith greatly contributed to the improvement of the manuscript. Emilie Berlioz thanks Ghislain Thierry for his contribution to the R script.

**Author Contributions**

**Conceptualization:** Emmanuel Discamps.

**Data curation:** Emilie Berlioz.
**Formal analysis:** Emilie Berlioz.

**Funding acquisition:** Emmanuel Discamps.

**Investigation:** Emilie Berlioz, Eugénie Capdepon.

**Methodology:** Emilie Berlioz.

**Project administration:** Emmanuel Discamps.

**Supervision:** Emmanuel Discamps.

**Validation:** Emilie Berlioz.

**Visualization:** Emilie Berlioz, Eugénie Capdepon, Emmanuel Discamps.

**Writing – original draft:** Emilie Berlioz, Eugénie Capdepon, Emmanuel Discamps.

## References

1. 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. Quaternary International. 2017; 433:64–87.

2. Frid A, Dill L. Human-caused disturbance stimuli as a form of predation risk. Conservation Ecology. 2002; 6(1):1–16.

3. Tolon V, Dray S, Loison A, Zeileis A, Fischer C, Baubet E. Responding to spatial and temporal variations in predation risk: space use of a game species in a changing landscape of fear. Canadian Journal of Zoology. 2009; 87(12):1129–37.

4. Bonnot N, Morellet N, Verheyden H, Cargnelutti B, Lourtet B, Klein F, et al. Habitat use under predation risk: hunting, roads and human dwellings influence the spatial behaviour of roe deer. European journal of wildlife research. 2013; 59(2):185–93.

5. Berlioz E. Écologie alimentaire et paléoenvironnements des cervidés européens du Pléistocène inférieur: le message des textures de micro-usure dentaire [PhD Thesis]. Université de Poitiers, France; 2017.

6. Kerley GIH, Kowalczyk R, Cromsigt J. Conservation implications of the refugee species concept and the European bison: king of the forest or refugee in a marginal habitat? Ecography. 2012; 35(6):519–29.

7. Merceron G, Hofman-Kamińska E, Kowalczyk R. 3D dental microwear texture analysis of feeding habits of sympatric ruminants in the Białowieża Primeval Forest, Poland. Forest Ecology and Management. 2014; 328:262–9.

8. Bocherens H, Hofman-Kamińska E, Drucker DG, Schmölcke U, Kowalczyk R. European bison as a refugee species? Evidence from isotopic data on Early Holocene bison and other large herbivores in northern Europe. PLoS ONE. 2015; 10(2):e0115090. https://doi.org/10.1371/journal.pone.0115090 PMID: 25671634

9. Hofman-Kamińska E, Bocherens H, Drucker DG, Fyfe RM, Gumiński W, Makowiecki D, et al. Adapt or die—Response of large herbivores to environmental changes in Europe during the Holocene. Global Change Biology. 2019; 25(9):2915–30. https://doi.org/10.1111/gcb.14733 PMID: 31298814

10. Dietl GP, Kidwell SM, Brenner M, Burney DA, Flessa KW, Jackson ST, et al. Conservation paleobiology: leveraging knowledge of the past to inform conservation and restoration. Annual Review of Earth and Planetary Sciences. 2015; 43:79–103.

11. Louys J. Paleontology in ecology and conservation. Springer Science & Business Media; 2012.

12. Scott RS, Ungar PS, Bergstrom TS, Brown CA, Childs BE, Teaford MF, et al. Dental microwear texture analysis: technical considerations. Journal of Human Evolution. 2006; 51(4):339–49. https://doi.org/10.1016/j.jhevol.2006.04.006 PMID: 16908052

13. Scott JR. Dental microwear texture analysis of extant African Bovidae. Mammalia. 2012; 76(2):157–74.

14. Calandra I, Merceron G. Dental microwear texture analysis in mammalian ecology. Mammal Review. 2016; 46(3):215–28.

15. Teaford MF, Glander KE. Dental microwear in live, wild-trapped *Alouatta palliata* from Costa Rica. American Journal of Physical Anthropology. 1991; 85(3):313–9.

16. Merceron G, Escarguel G, Angibault JM, Verheyden-Tixier H. Can dental microwear textures record inter-individual dietary variations? PLoS ONE. 2010; 5(3):1–9. https://doi.org/10.1371/journal.pone.0009542 PMID: 20209051
17. Winkler DE, Schulz-Kornas E, Kaiser TM, Codron D, Leichliter J, Hummel J, et al. The turnover of dental microwear texture: Testing the "last supper" effect in small mammals in a controlled feeding experiment. *Palaeogeography, Palaeoclimatology, Palaeoecology*. 2020; 557:109930.

18. Laquay G. Recherches sur les faunes du Würm I en Périgord [PhD Thesis]. Verlag nicht ermittelbar; 1981.

19. Binford LR. The diet of early hominins: Some things we need to know before "reading" the menu from the archaeological record. Guts and brains: An integrative approach to the hominin record. 2007;185-222.

20. Guadelli JL. Les associations de grands mammifères de la séquence du pléistocène supérieur de Combe-Grenal (Dordogne, France) et la biochronologie des stades isotopiques 5 à 3. Large mammal assemblages in the upper pleistocene sequence of Combe-Grenal (Dordogne, France) and the biochronology of isotopic stages 5 to 3. In: 134e congrès du CTHS, Colloque 32 François Bordes. Les éditions du CTHS; 2012. p. 199-219.

21. Guadelli JL. Contribution à l'étude des zoocénoses préhistoriques en Aquitaine (Würm ancien et inter-stade würmiem [PhD Thesis]. Bordeaux 1; 1987.

22. Discamps E, Faivre JP. Substantial biases affecting Combe-Grenal faunal record cast doubts on previous models of Neanderthal subsistence and environmental context. *Journal of Archaeological Science*. 2017; 81:128-32.

23. Richards MP, Pellegrini M, Niven L, Nehlich O, Dibble H, Turq A, et al. Temporal variations in Equus tooth isotope values (C, N, O) from the Middle Paleolithic site of Combe Grenal, France (ca. 150,000 to 50,000 BP). *Journal of Archaeological Science: Reports*. 2017; 14:189-98.

24. Garralda MD, Vandermeersch B. Les Néandertaliens de la grotte de Combe-Grenal (Domme, Dordogne, France)/The Neanderthals from Combe-Grenal cave (Domme, Dordogne, France). Paléo, Revue d'Archéologie Préhistorique. 2000; 12(1):213-59.

25. Bordes F. A tale of two caves. *HarperCollins Publishers*; 1972.

26. Faivre JP, Discamps E, Gravina B, Turq A, Guadelli JL, Lenoir M. The contribution of lithic production systems to the interpretation of Mousterian industrial variability in south-western France: The example of Combe-Grenal (Dordogne, France). *Quaternary International*. 2014; 350:227-40.

27. Genty D, Blamart D, Ouahdi R, Gilmour M, Baker A, Jouzel J, et al. Precise dating of Dansgaard-Oeschger climate oscillations in western Europe from stalagmite data. *Nature*. 2003; 421(6925):833-7.
https://doi.org/10.1038/nature01391 PMID: 12594510

28. Wainer K, Genty D, Blamart D, Daëron M, Bar-Matthews M, Vonhof H, et al. Speleothem record of the last 180 ka in Villars cave (SW France): Investigation of a large δ18O shift between MIS6 and MIS5. *Quaternary Science Reviews*. 2011; 30(1-2):130-46.

29. Sanchez-Goñi MFS, Landais A, Fletcher WJ, Naughton F, Desprat S, Duprat J. Contrasting impacts of Dansgaard-Oeschger events over a western European latitudinal transect modulated by orbital parameters. *Quaternary Science Reviews*. 2008; 27(11-12):1136-51.

30. Genty D, Combourieu-Nebout N, Peyron O, Blamart D, Wainer K, Mansuri F, et al. Isotopic characterization of rapid climatic events during OIS3 and OIS4 in Villars Cave stalagmites (SW-France) and correlation with Atlantic and Mediterranean pollen records. *Quaternary Science Reviews*. 2010; 29 (19-20):2799-820.

31. 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. *Quaternary Science Reviews*. 2011; 30(19-20):2755-75.

32. Garralda MD, Giacobini G, Vandermeersch B. Neanderthal cutmarks: combe-grenal and marillac (France). A SEM analysis. *Anthropologie* (1962-). 2005; 43(2/3):189-98.

33. Dayet L, Faivre JP, Le Bourdonnec FX, Discamps E, Royer A, Claud É, et al. Manganese and iron oxide use at Combe-Grenal (Dordogne, France): A proxy for cultural change in Neanderthal communities. *Journal of Archaeological Science: Reports*. 2019; 25:239-56.

34. Bignon-Lau O, Catz N, Berlioz E, Veiberg V, Strand O, Merceron G. Dental microwear textural analyses to track feeding ecology of reindeer: a comparison of two contrasting populations in Norway. *Mammal Research*. 2017; 62(1):111-20.

35. Hermier R, Merceron G, Kostopoulos DS. The emblematic Eurasian Villafranchian antelope Gazellos-pira (Mammalia: Bovidae): New insights from the Lower Pleistocene Dafnero fossil sites (Northern Greece). *Geobios*. 2020; 61:11-29.

36. Merceron G, Berlioz E, Vonhof H, Green D, Garel M, Tütken T. Tooth tales told by dental diet proxies: An alpine community of sympatric ruminants as a model to decipher the ecology of fossil fauna. *Palaeogeography, Palaeoclimatology, Palaeoecology*. 2021; 562:110077.
37. Merceron G, Ramdarshan A, Blondel C, Boisserie JR, Brunetiere N, Francisco A, et al. Untangling the environmental from the dietary: dust does not matter. Proceeding of the Royal Society B. 2016; 283:1–8.

38. Discamps E. Hommes et hyènes face aux recompositions des communautés d'Ongulés (MIS 5–3): élé- ments pour un cadre paléoécologique des sociétés du Paléolithique moyen et supérieur ancien d'Europe de l'Ouest. [Bordeaux]: Université Bordeaux 1, France; 2012.

39. Smith C, Warren M. GLMs in R for ecology. Smith and Warren; 2019.

40. Zuur AF, Ieno EN, Elphick CS. A protocol for data exploration to avoid common statistical problems. Methods in Ecology and Evolution. 2010; 1(1):3–14.

41. Komsta L. outliers: Tests for Outliers [Internet]. 2022 [cité 7 oct 2022]. Disponible sur: https://CRAN.R-project.org/package=outliers

42. Wickham H, Chang W, Wickham MH. Package 'ggplot2'. Create elegant data visualisations using the grammar of graphics Version. 2016; 2(1):1–189.

43. Gastwirth JL, Gel YR, Hui WW, Lyubchich V, Miao W, Noguchi K, et al. Package 'lawstat'. R Founda- tion for Statistical Computing, Vienna, Austria; 2019.

44. Dag O, Dolgun A, Konar NM, Weerahandi S, Ananda M. onewaytests: One-Way Tests in Independent Groups Designs [Internet]. 2021 [cité 7 oct 2022]. Disponible sur: https://CRAN.R-project.org/package=onewaytests

45. Ripley B, Venables B, Bates DM, Hornik K (partial port ca 1998), Gebhardt A (partial port ca 1998), Firth D. MASS: Support Functions and Datasets for Venables and Ripley's MASS [Internet]. 2022 [cité 7 oct 2022]. Disponible sur: https://CRAN.R-project.org/package=MASS

46. Jr FEH, CD. (contributed several functions and maintains latex functions). Hmisc: Harrell Miscellaneous [Internet]. 2022 [cité 7 oct 2022]. Disponible sur: https://CRAN.R-project.org/package=Hmisc

47. Wei T, Simko V, Levy M, Xie Y, Jin Y, Zemla J, et al. corrplot: Visualization of a Correlation Matrix [Internet]. 2021 [cité 7 oct 2022]. Disponible sur: https://CRAN.R-project.org/package=corrplot

48. Marshall R. regplot: Enhanced Regression Nomogram Plot [Internet]. 2020 [cité 7 oct 2022]. Disponible sur: https://CRAN.R-project.org/package=regplot

49. Mazerolle MJ. AICcmodavg: Model Selection and Multimodel Inference Based on (Q)AIC(c) [Internet]. 2020 [cité 7 oct 2022]. Disponible sur: https://CRAN.R-project.org/package=AICcmodavg

50. Lenth RV, Buerkner P, Herve M, Jung M, Love J, Miguez F, et al. emmeans: Estimated Marginal Means, aka Least-Squares Means [Internet]. 2022 [cité 7 oct 2022]. Disponible sur: https://CRAN.R-project.org/package=emmeans

51. Team RC. R: A language and environment for statistical computing. R Foundation for Statistical Com- puting, Vienna, Austria. http://www/R-project.org/. 2013.

52. Sarkar D, Andrews F, Wright K (documentation), Klepeis N, Larsson J (colorkey title), Murrell P. lattice: Trellis Graphics for R [Internet]. 2021 [cité 7 oct 2022]. Disponible sur: https://CRAN.R-project.org/package=lattice

53. Ogle D [aut, cre], Doll J, Wheeler P, Dinno A (Provided base functionality of dunnTest()). FSA: Simple Fisheries Stock Assessment Methods [Internet]. 2022 [cité 7 oct 2022]. Disponible sur: https://CRAN.R-project.org/package=FSA

54. Gebert C, Verheyden-Tixier H. Variations of diet composition of red deer (Cervus elaphus L.) in Europe. Mammal Review. 2001; 31(3-4):189–201.

55. Azorit C, Tellado Ruiz MS, Oya A, Moro J. Seasonal and specific diet variations in sympatric red and fal- low deer of southern Spain: a preliminary approach to feeding behaviour. Animal Production Science. 2012; 52(8):720–7.

56. Geist V. Deer of the world: their evolution, behaviour, and ecology. Stackpole books; 1998.

57. Lovari S, Lorenzini R, Masseti M, Pereladova O, Carden RF, Brook SM. *Cervus elaphus*. The IUCN Red List of Threatened Species 2016. 2015.

58. Berlioz E, Kostopoulos DS, Blondel C, Merceron G. Feeding ecology of *Eucladoceros ctenoides* as a proxy to track regional climatic oscillations in Europe during the Early Pleistocene. Comptes Rendus Palevol. 2017; in press.

59. Skogland T. The effects of food and maternal conditions in fetal growth and size in wild reindeer. Rang- fer. 1984; 4(2):39–46.

60. Skogland T. Natural selection of wild reindeer life history traits by food limitation and predation. Oikos. 1989; 55(1):101–10.

61. Gunn A. *Rangifer tarandus*. The IUCN Red List of Threatened Species 2016. 2016.
62. Telfer ES, Kelsall JP. Adaptation of some large North American mammals for survival in snow. Ecology. 1984; 65(6):1828-34.

63. Klein DR. Comparative ecological and behavioral adaptations of Ovibos moschatus and Rangifer tarandus. Rangifer. 1992; 12(2):47-55.

64. Gaare E, Thomson BR, Kjos-Hanssen O. Reindeer activity on Hardangervidda. In: Fennoscandian Tundra Ecosystems. Springer; 1975. p. 206-15.

65. Nieminen M, Heiskari U. Diets of freely grazing and captive reindeer during summer and winter. Rangifer. 1989; 9(1):17-34.

66. Klein DR. The introduction, increase, and crash of reindeer on St. Matthew Island. The Journal of Wildlife Management. 1968; 32(2):350-67.

67. Geist V, Bayer M. Sexual dimorphism in the Cervidae and its relation to habitat. Journal of Zoology. 1988; 214(1):45-53.

68. Moen J, Andersen R, Illius AW. Living in a seasonal environment. In: Danell K, Duncan P, Bergstrom R, Pastor J, éditeurs. Large herbivore ecology, ecosystem dynamics and conservation. Cambridge University Press. Cambridge; 2006. p. 50-70.

69. Soubrier J, Gower G, Chen K, Richards SM, Llamas B, Mitchell KJ, et al. Early cave art and ancient DNA record the origin of European bison. Nature communications. 2016; 7(1):1-7. https://doi.org/10.1038/ncomms13158 PMID: 27754477

70. Markova AK, Kolfschoten T van, Puzachenko A. Evolution of European ecosystems during pleistocene-holocene transition (24-8 Kyr BP). 2019.

71. IUCN SSC Bison Specialist Group, European Bison Specialist Group. European Bison (Bison bonasus) Strategic Species Status Review 2020. Olech, W. and K. Perzanowski. 2022.

72. Guthrie RD. Frozen fauna of the mammoth steppe: the story of Blue Babe. In: Frozen Fauna of the Mammoth Steppe. University of Chicago Press; 1990.

73. Markova AK, Puzachenko AY, Van Kolfschoten T, Kosintsev PA, Kuznetsova TV, Tikhonov AN, et al. Changes in the Eurasian distribution of the musk ox (Ovibos moschatus) and the extinct bison (Bison priscus) during the last 50 ka BP. Quaternary International. 2015; 378:99-110.

74. Julien MA, Bocherens H, Burke A, Drucker DG, Patou-Mathis M, Krotova O, et al. Were European steppe bison migratory? 18O, 13C and Sr intra-tooth isotopic variations applied to a palaeoethological reconstruction. Quaternary International. 2012; 271:106-19.

75. Stuart AJ. Mammalian extinctions in the Late Pleistocene of northern Eurasia and North America. Biological Reviews. 1991; 66(4):453-562. https://doi.org/10.1111/j.1469-185x.1991.tb01149.x PMID: 1801948

76. Lorenzen ED, Nogués-Bravo D, Orlando L, Weinstock J, Binladen J, Marske KA, et al. Species-specific responses of Late Quaternary megafauna to climate and humans. Nature. 2011; 479(7373):359-64. https://doi.org/10.1038/nature10574 PMID: 22048313

77. Merceron G, Madelaine S. Molar microwear pattern and palaeoecology of ungulates from La Berbie (Dordogne, France): environment of Neanderthals and modern human populations of the Middle/Upper Palaeolithic. Boreas. 2006; 35(2):272-8.

78. Rivals F, Álvarez-Lao DJ. Ungulate dietary traits and plasticity in zones of ecological transition inferred from late Pleistocene assemblages at Jou Puerta and Rexidora in the Cantabrian Region of northern Spain. Palaeogeography, Palaeoclimatology, Palaeoecology. 2018; 499:123-30.

79. Wright E, Viner-Daniels S. Geographical variation in the size and shape of the European aurochs (Bos primigenius). Journal of Archaeological Science. 2015; 54:8-22.

80. Van Vuure C. Retracing the aurochs. In: Neue Modelle zu Masznahmen der landschaftsentwicklung mit groszen Pflanzenfressern und praktische Erfahrungen bei der Umsetzung. 2005. p. 261-5.

81. Helmer D, Monchot H. Un site mésolithique de chasse à l'aurochs (La Montagne, Sénas, Bouches-du-Rhône). Anthropozoologica. 2006; 41(2):215-28.

82. Wright E. The history of the European aurochs (Bos primigenius) from the Middle Pleistocene to its extinction: an archaeological investigation of its evolution, morphological variability and response to human exploitation. [PhD Thesis]. University of Sheffield; 2013.

83. Noe-Nygaard N, Price TD, Hede SU. Diet of aurochs and early cattle in southern Scandinavia: evidence from 15N and 13C stable isotopes. Journal of Archaeological Science. 2005; 32(6):855-71.

84. Marquet JC. Paléoenvironnement et chronologie des sites du domaine atlantique français d'âge Pléistocène moyen et supérieur d'après l'étude des rongeurs. Les Cahiers de la Claise, Supplément 2. 1993.

85. Davis M, Pineda-Munoz S. The temporal scale of diet and dietary proxies. Ecology and evolution. 2016; 6(6):1883-97. https://doi.org/10.1002/ece3.2054 PMID: 27087936
86. Steele TE. Red deer: their ecology and how they were hunted by Late Pleistocene hominids in western Europe. Stanford University; 2002.

87. Sánchez-Hernández C, Gourichon L, Pubert E, Rendu W, Montes R, Rivals F. Combined dental wear and cementum analyses in ungulates reveal the seasonality of Neanderthal occupations in Covalejos Cave (Northern Iberia). Scientific reports. 2019; 9(1):1–18.