Alternating carnivore and Neanderthal activities at Escoural Cave: insights from the taphonomic and machine learning analysis of leporid remains

Exploring the varied subsistence strategies and cave occupation patterns of Neanderthals is key to understanding their complex behaviors and ecological adaptations. Small game consumption, in particular, is considered a relevant indicator of their behavioral complexity. Rabbit assemblages from Pleistocene cave sites provide valuable insights into Neanderthal interactions with small prey and potential competition with carnivores. Here, we present the first detailed taphonomic analysis of faunal remains from Escoural Cave (Portugal), where a European rabbit (Oryctolagus cuniculus) assemblage was found alongside Middle Paleolithic stone tools and some macromammal remains. This study combines traditional zooarchaeological and taphonomic analysis of the rabbit remains with multivariate statistics and machine learning methods to establish the origin of the accumulation, and the implications for Neanderthal subsistence and cave use. Results from the taphonomic analysis show no evidence of human consumption but abundant evidence of small terrestrial carnivore activity, primarily from lynxes. This could indicate a sequential occurrence of Neanderthal and carnivore activities in the cave, with Neanderthal activities likely related to something other than rabbit consumption. Our study contributes to characterizing Iberian carnivore fossil accumulations and differentiating between faunal assemblages accumulated by carnivores and those by hominins. Additionally, we show that the use of machine learning analysis provides a robust and objective method for identifying and classifying taphonomic signatures, enhancing the accuracy and reliability of our interpretations. Future work will focus on analyzing additional faunal collections from both past and new excavations at Escoural, to test whether carnivores and Neanderthals focused on different prey types and how they shared the cave space.

1 Introduction

In recent years, our understanding of Neanderthal complex behavior has significantly increased. Studies have demonstrated diverse and context-dependent strategies, including advanced planning, social structures (e.g., Duveau et al., 2019; Skov et al., 2022), varying landscape use and mobility (e.g., Marín et al., 2019; Valde-Nowak and Cieśla, 2020; Moncel et al., 2019), and campsite organization (e.g. Gabucio et al., 2023; Vaquero, 2022; Roebroeks and Soressi, 2016). Regarding their subsistence behavior, while Neanderthals often focused on hunting ungulates (e.g. Daujeard et al., 2019; Daujeard and Moncel, 2010; Gaudzinski-Windheuser and Kindler, 2012; Livraghi et al., 2021; Marín et al., 2017, 2020; Richards et al., 2000; Smith, 2014; Yravedra and Cobo-Sánchez, 2015; Yravedra et al., 2024), they also occasionally exploited large carnivores, such as lions, bears and wolves, as suggested by evidence of hunting and butchery marks (Russo et al., 2023; Auguste, 1995; Herranz-Rodrigo et al., 2024; Buccheri et al., 2016; Hussain et al., 2022). Moreover, the hunting and consumption of small game, such as leporids, have also been documented and recognized as evidence of sophisticated foraging strategies and adaptive behaviors (Pérez Ripoll, 2001; Blasco et al., 2022). Research has shown that small game, including leporids, birds, tortoises, small terrestrial carnivores, and marine resources, also played a role in Neanderthal diet (Real, 2020; Blasco and Fernández Peris, 2012; Brown et al., 2011; Carvalho et al., 2018; Cochard et al., 2012; Gabucio et al., 2014; Gómez-Olivencia et al., 2018; Gutiérrez-Zugasti et al., 2018; Hardy et al., 2013; Nabais and Zilhão, 2019; Pelletier et al., 2019; Romandini et al., 2023; Rufà et al., 2014; Stiner et al., 1999, 2000).

Rabbits occur abundantly at archaeological sites and are present across various biotopes in the Iberian Peninsula given their high reproductive rates and rapid adaptability (Delibes and Hiraldo, 1981; Pérez Ripoll, 1992; Rodríguez-Hidalgo et al., 2020; Rufà et al., 2014). Understanding the origins of leporid accumulations has long been an important focus for taphonomists, not only because they were a part of human subsistence during the late Pleistocene and early Holocene but also due to the complexities they present, as these accumulations can result from many different processes (e.g. Pérez Ripoll, 1992; Cochard, 2004; Yravedra et al., 2019). The most significant cause of death in rabbits recovered from archaeological Middle Paleolithic contexts is predator activity (e.g. Pérez Ripoll, 1977; Rufà et al., 2014). Rabbits are prey to around 40 predators in Western Europe (Delibes and Hiraldo, 1981; Cochard, 2004), including terrestrial carnivores such as the Iberian lynx (Lynx pardinus) (Lloveras et al., 2008a; Rodríguez-Hidalgo et al., 2013, 2015, 2020), foxes (Vulpes vulpes) (Cochard, 2004; Hockett and Haws, 2002; Krajcarz and Krajcarz, 2012; Lloveras et al., 2012a,b; Sanchis Serra and Pascual Benito, 2011), wildcats (Felis silvestris) (Cohen and Kibii, 2015; Lloveras et al., 2016; Scheifler et al., 2020), wolves (Canis lupus) (Schmitt and Juell, 1994; Armstrong, 2016; Lloveras et al., 2020), as well as birds of prey, including nocturnal (e.g., Lloveras et al., 2009, 2012a,b; Sanchis, 1999; Yravedra, 2004) and diurnal raptors (e.g. Hockett, 1996; Lloveras et al., 2008b).

Rabbit bone accumulations can result from a variety of processes that are not mutually exclusive, with several agents contributing to leporid assemblages within the same deposit. According to Callou's (1997) classification, these accumulations can be categorized as intrusive (e.g., accidental deaths within the site due to rabbits entering or burrowing), exogenous (e.g., remains brought into the site by predators), or anthropic (i.e., associated with human activities such as hunting, butchering or consumption). This classification also encompasses natural causes, such as accidental deaths in natural traps or burrows, where the accumulation occurs without direct involvement from other predators (Pelletier et al., 2016; Seuru et al., 2024). Thus, rabbit bones in fossil assemblages may reflect contributions from non-human predators, humans, natural processes, or a combination of these factors (e.g. Callou, 1997; Yravedra, 2008; Yravedra and Andrés, 2013).

Identifying the specific accumulation origins is essential for understanding broader ecological interactions between Neanderthals and carnivores. These interactions can shed light on the evolutionary pressures that shaped hominin behavior and reveal critical aspects of Neanderthal diet, competition for resources, and foraging strategies, while also providing insights into Neanderthal use of shelters, risk mitigation, and responses to carnivore presence (e.g. Camarós et al., 2017; Hussain et al., 2022; Rivals et al., 2022; Romandini et al., 2023; Rosell et al., 2017; Linares-Matás and Yravedra, 2024).

Faunal accumulations are often palimpsests of multiple overlaid occupation events produced by several bone-modifying agents, which is challenging for archaeological interpretation. Additionally, complex site formation processes, particularly in cave sites, can obscure the original context and sequence of events (e.g., Rufà et al., 2014; Saladié et al., 2018; Linares-Matás and Yravedra, 2024). Despite these complexities, researchers attempt to classify the different types of dynamics or co-occurrences of hominin and carnivore activity at sites along a continuum of scenarios, which considers both the timing (temporality) and the nature (proportion and types) of the activities carried out by each agent: (a) sequential succession of hominin and carnivore occupations, (b) concurrent events of carnivore and human activity (with the presence of stone tools and fireplaces), (c) human settlement with sporadic carnivore visits, and (d) carnivore den with occasional human occupation (Picin et al., 2020; Villa et al., 2004; Linares-Matás and Yravedra, 2024; Pinto-Llona et al., 2023; Yravedra, 2006). Taphonomy has the potential to enhance the accuracy and reliability of these classifications by incorporating advanced statistical methodologies. However, current approaches often lack the robustness needed to clearly differentiate the origins of rabbit accumulations, especially at Pleistocene and early Holocene sites where evidence is often ambiguous. This highlights the need for new, more rigorous methods to better identify human involvement in complex assemblages. To address this gap, our study applies advanced statistical methods and contributes to developing a more comprehensive methodological framework for future research.

Southwestern Iberia's Middle Paleolithic sites are critical for understanding Neanderthal flexibility and resilience in response to changing environmental conditions and carnivore presence, as this region may have served as a refugium for Neanderthals. At some of these sites, which present the challenges of containing palimpsests and complex deposits, there is a need for re-excavation and re-analysis using modern analytical techniques to obtain more accurate and comprehensive data. One significant gap in current research in the region is the lack of faunal taphonomic analysis, which is essential for understanding the interactions between Neanderthals, their environment, and other fauna. The present study is part of a broader project (FINISTERRA: ERC-2021-COG-101045506) aimed at excavating and reanalyzing several such sites.

Here we present a detailed taphonomic and zooarchaeological analysis of a sample of rabbit bones recovered from the Middle Paleolithic site of Escoural Cave to assess the role Neanderthals or other agents played in the formation of the assemblage. Although this analysis is preliminary and based on a subset of materials, it represents the first comprehensive taphonomic study conducted for Escoural cave and offers essential contextual information into the role of small prey species in Neanderthal-carnivore interactions and the use of cave spaces. Additionallly, the methods developed here will provide a foundational framework for future studies at Escoural.

1.1 Site overview

The Escoural Cave or Gruta do Escoural (38°32′55″ N, −8°8′18″E, Alentejo, Portugal) is a dolomitic cave located about 15 km from the city of Montemor-o-Novo, near the village of Santiago do Escoural (3 km) (Lattao et al., 2023). Discovered in 1963 during mining activities for extracting marble, the cave holds significant potential for investigating Neanderthal-carnivore interactions during the Middle Paleolithic. The karstic system is composed of a main gallery, which spans ~12 × 15 m, connected with ~17 large karstic channels of around one meter in diameter (Figure 1). The original entrance collapsed after the Middle Paleolithic occupations, although investigations are ongoing to establish the exact timing of this event.

Figure 1

Figure 1. (A) Geographic location of Escoural Cave in Portugal; (B) cave outline and location of the test profiles excavated in 2020 inside (P1) and outside (P2); (C) stratigraphic section of P1 showing the defined geological levels; (D) spatial vertical distribution of the recovered archaeological materials in each level.

Excavations were conducted during the 1960s and 1990s and unearthed abundant archaeological findings, including Chalcolithic occupations outside the cave, Neolithic burials, Upper Paleolithic engravings and paintings, and archaeological remains from the Upper and Middle Paleolithic. The Middle Paleolithic levels were characterized by abundant quartz industries and faunal remains. Close to 600 lithic artifacts were reported from the original excavations inside the cave (Otte and da Silva, 1996), while our work in 2021 recovered 206 additional lithics from P1. Preliminary analysis of this assemblage aligns largely with previous findings, indicating a predominantly expedient technology, alongside the presence of blanks obtained through discoidal and Levallois methods. U-Th dating of a horse tooth suggested an age of ~48 Ka for the Middle Paleolithic levels identified inside the cave (Otte and da Silva, 1996). Current sampling and dating efforts aim to further refine the chronology.

1.2 Previous studies

Earlier research on the faunal remains recovered from the Middle and Upper Paleolithic levels at the cave indicate it was occupied by humans and other carnivores (Otte and da Silva, 1996). Although contributions would have come from both groups, carnivores were likely the primary accumulating agents. The reported identified assemblage represents 27.8% of the total, increasing to 73.4% when including rabbit remains. Despite being very fragmented, the assemblage showed good preservation. The macrofaunal collection is dominated by several species of equids and cervids but also includes a diverse taxonomic list including bovids, suids, ibex, and carnivores (see Supplementary Table S1). Smaller fauna include rabbits, birds, mollusks, gastropods, reptiles, and microfauna (see Otte and da Silva, 1996 for the exhaustive taxonomic list). Carnivore activity on macromammals was mainly attributed to hyenas, alongside digested bones and coprolites. However, a portion of the assemblage, including gastropods, frogs, lizards, bats, most carnivores, and rodents, would have been mainly accumulated by nocturnal birds of prey. Birds were likely accumulated by wildcats or lynxes, while toads and rabbits could have been collected by either carnivores or humans. Hyenas are believed to have been the primary accumulating agent, and potentially responsible for much of the assemblage, although lynxes and wolves would also have contributed by hunting various animals, including rabbits, hares, deer, and young suids (Otte and da Silva, 1996).

The authors of the study also observed that the ratio of large animal remains relative to carnivores increased from the Middle Paleolithic, suggesting changes in environmental conditions or selective hunting by carnivores. Equid remains also appear to increase in the upper levels, potentially indicating a drier and colder environment. However, in the Upper Paleolithic levels, the number of identified specimens decreases due to higher fragmentation or increased mixing and reshaping of the deposits, which suggests that these observations should be taken with caution. Additionally, their findings were based on the Number of Identified Specimens (NISP), including various ratios used in the analysis. The reliance on NISP requires careful interpretation, as estimates may not be accurate due to fragmentation of the archaeofaunal collection. Moreover, the absence of detailed taphonomic analyses or reported modification frequencies underscores the need for further studies, including both the re-examination of previously excavated materials and the analysis of new finds.

2 Materials and methods

2.1 Materials

This study focuses on materials recovered during recent excavations from the Middle Paleolithic layers of Escoural cave. The analysis of additional materials housed in Lisbon is planned but currently on hold due to museum renovations. Meanwhile, the FINISTERRA project (https://finisterra.icarehb.com/) has initiated new excavations at the site opening and excavating larger areas to understand the connection between the exterior and interior of the cave and reconstruct site formation processes. Excavations began in 2020 with the expansion of the profiles left by previous excavation areas in one of the SE chambers (Profile 1 in Gallery 17—see Figure 1B) and outside the cave (Profile 2), and included sampling for dating and micromorphological analyses, the latter of which are currently ongoing (Alzate-Casallas et al., in prep.). This study focuses on remains from Profile 1, which yielded more abundant and better preserved materials than Profile 2. The vertical spatial distribution of the remains in Profile 1 revealed an overlap between bones and lithics, primarily in Levels 2, 3, and 4 (Figures 1C, D). Faunal remains decreased drastically in Level 5 and were absent in Levels 6 and 7. Preliminary microstratigraphic analyses suggest potential reworking of the sediments; however, distinct depositional episodes are evident, as multiple geological levels could be identified. The deposition of these levels was not continuous but alternated with breaks during which calcite crusts formed.

Profile 1 inside Escoural Cave yielded macrofaunal remains; yet, the sample size (N = 88) was insufficient for a comprehensive zooarchaeological and taphonomic analysis. However, rabbit remains were abundant in Geological Levels (GL) 2, 3, and 4. This study focuses on the rabbit remains from GL 2, where the sample size was N = 1,313. In this level, macromammal remains amounted to N = 37, and plotted rabbit remains to N = 27. These 27 plotted remains, recovered due to their larger size and identifiability, are each associated with precise spatial coordinates, whereas the majority of specimens were recovered from sieve bags, resulting in only approximate spatial locations corresponding to each 10 L bucket of excavated sediment. Each rabbit bone specimen recovered from the sieve bags was assigned a unique identification number in the lab to ensure the reproducibility of the analysis.

2.2 Methods

Zooarchaeological and taphonomic attributes of the bones were recorded using E5 data entry program (McPherron, 2022), and analyzed according to anatomical representation, breakage, bone surface modifications, and digestion. Bone breakage planes and cortical surfaces were inspected using hand lenses and a Euromex stereomicroscope (Nexius Zoom NZ 1902-P) with magnification up to 45 × . Wherever possible, the resulting data were treated statistically by comparing them to other archaeological and actualistic rabbit assemblages to enhance the robustness of the analysis and help identify the accumulating agent. Assemblages were selected based on the availability of detailed taphonomic and zooarchaeological data presented in a consistent format to enable reliable comparison and analysis. We used the dataset and site selection criteria from Rodríguez-Hidalgo et al. (2020) as a baseline for our analysis, adopting their format of taphonomic variables. The variables included in our analyses are MAU% or RA% data, as well as percentages of cut marks, burned remains, digested remains, tubes/cylinders, tooth/beak marks, complete long bones, and adults. Sites with very incomplete data on these variables were excluded. Additionally, a minimum threshold of 9 MNI was established to ensure that the percentages were representative of the assemblage, and sites or actualistic studies with smaller rabbit assemblages were not included. The following sections outline the specific methodologies used.

The entirety of the R code used for the analysis, datasets, and visual representations contained in this paper can be accessed through our online research compendium at https://osf.io/zyu4p/. We used the rrtools package by Marwick et al. (2018) to create a research compendium and write a reproducible journal article. The provided files include the complete set of raw data used in the analysis, along with a custom R project (Wickham, 2015) containing the code required to generate all tables and figures. To enable maximum reuse, the code is made available under the MIT license, data under CC-0, and figures under CC-BY (additional details can be found in Marwick, 2017).

2.2.1 Taxonomic identification and age profiles

The zooarchaeological analysis involved identifying animal bones to the genus or species level where possible, and the selection of rabbit remains for detailed analysis. The taxonomic identification of leporid remains was based on comparative collections housed at ICArEHB laboratory and reference bibliography (e.g. Barone et al., 1973; Callou, 1997). Following established zooarchaeological methods (e.g. Reitz and Wing, 1999), we determined the skeletal element, the preserved portion of the bone, and bone laterality. We used the stage of epiphyseal fusion to establish the approximate age of rabbits at death (e.g. Jones, 2006; Callou, 2003; Cochard, 2004; Hockett and Bicho, 2000). Specifically, we used the distal epiphysis of the humerus, which begins to fuse at about 2 months of age, and the distal epiphysis of the tibia, which fuses at around 3 months, to identify juveniles (<3 months) (Jones, 2006). The remaining epiphyses of the long bones typically fuse between 9 and 10 months (Pelletier et al., 2019), which marks the threshold between subadults and adult individuals. The latter are characterized by the complete fusion of the proximal epiphysis of the humerus and tibiae (Pelletier et al., 2019).

2.2.2 Anatomical representation

We calculated the Number of Identified Specimens (NISP) (Lyman, 1994), the Minimum Number of Elements (MNE) (Bunn, 1982; Marean et al., 2001), the Minimum Number of Individuals (MNI) (White, 1953; Lyman, 2008), as well as the Minimum Number of Animal Units (%MAU) (Binford, 1978). The MNE was estimated by counting the number of times the most represented portion of a certain element was present. This estimation involved including the consideration of diaphyses as well as epiphyses (Marean et al., 2001). Long bone shafts were identified when possible to a specific element using anatomical landmarks like muscle insertions or foramina, as well as diaphyseal cross-section and medullary cavity shape and size (Yravedra and Domínguez-Rodrigo, 2009). MNI was determined by considering the most abundant anatomical element, bone laterality and animal age. MAUs were calculated from the MNE values, as established by Binford (1978, 1981, 1984), by dividing the MNEs of each skeletal part by the number of times that such element occurs in the animal's skeleton. Each resulting MAU was subsequently transformed into %MAU by standardizing each MAU with reference to the largest original MAU, also as described originally by Binford (1978, 1981). Following the established methodology for the study of leporid assemblages, we calculated the relative abundance (%RA) (Dodson and Wexlar, 1979; Andrews, 1990). In addition, the proportions of skeletal elements were evaluated using the indices outlined in Table 1 and described and applied in studies such as Andrews (1990), Lloveras et al. (2008a) or Rodríguez-Hidalgo et al. (2020).

Table 1

Table 1. Indices used to evaluate the proportions of leporid skeletal elements.

Moreover, to evaluate whether biases in skeletal representation were consequence of bone density-mediated attrition or influenced by other post-depositional or biostratinomic processes affecting the preservation of specific skeletal parts, we calculated Spearman's correlation coefficient (Pavao and Stahl, 1999). This coefficient determines if there is a significant correlation between the bone density values of rabbit bones and the relative abundance of each skeletal element.

2.2.3 Principal component analysis (PCA) and machine learning (ML) analysis of skeletal part abundances

Subsequently, we performed a PCA analysis using skeletal part abundance data (%RA values) following the approach in Rodríguez-Hidalgo et al. (2020) to visualize the distribution of a set of assemblages and sites generated by different predators by reducing the dimensionality of the dataset and to identify potential patterns and clusters among the samples. We chose relative abundance data for the statistical analysis, because it provides standardized and comparable values, unlike the indices of skeletal part representation calculated from NISP and MNE, and thus ensures more reliable outcomes. Assemblages were included in the reference dataset only if their minimum number of individuals (MNI) was higher than 9. The reference dataset included assemblages generated by terrestrial predators (lynx, wildcat, fox, and coyote), diurnal predators (eagles and vultures), nocturnal predators (owls), and humans (Neanderthals and modern humans), as well as a combination of non-ingested and scat or pellet remains (see Supplementary Table S2 for complete dataset and references). We combined species into broader groups, as individual classifications would often comprise only one or two samples per group. This results in a dataset with increased variability, but enabled the implementation of this analysis and subsequent ML models. The PCA was performed using the factoextra library in R (Kassambara and Mundt, 2020). The dataset was preprocessed to handle missing values using the “missForest” library (Stekhoven and Buehlmann, 2012) and normalized.

Next, we explored the use of machine learning (ML) techniques to identify the agent of accumulation of rabbits at Escoural based on skeletal representation data. Machine learning models generally perform better with larger datasets, because small sample sizes can lead to overfitting, where the model learns noise and specific details of the training data rather than general trends. However, there are several techniques that can help mitigate this effect, especially feature selection, cross-validation, and hyperparameter tuning (Vabalas et al., 2019). To ensure robust models and minimize the risk of overfitting, we used a combination of these techniques. We first used Random Forest to select the 10 most important features of the dataset. Next, using a correlation matrix, we identified and removed the highest correlated features (>0.7) to reduce multicollinearity. Then, we split the data into training (70%) and test (30%) sets using stratified sampling. Stratified sampling is similar to k-fold, but ensures that each fold has approximately the same percentages of samples of each target class as the entire dataset, which can help in maintaining the balance of class distributions (Kuhn and Johnson, 2013). Subsequently we bootstrapped the training set with 100 iterations, involving repeatedly resampling the training data with replacement, coupled with internal cross-validation (Kuhn and Johnson, 2013). We trained three different ML models, Random Forest (RF), Support Vector Machine (SVM), and Gradient Boosting Machine (GBM), with hyperparameter tuning and cross-validation to ensure that they generalized well to unknown data. We calculated precision, recall, F1-score and other performance metrics across multiple resamples, and compared the outcomes of the three models. We used the “caret” (Kuhn, 2008) “e1071” (Meyer et al., 2023), and “randomForest” (Liaw and Wiener, 2002) libraries in R.

2.2.4 Breakage

All identified bone specimens were classified into the breakage categories defined and described by Lloveras et al. (2008b). Breakage categories vary depending on bone element, and percentages are also recorded for complete elements and isolated and in-situ teeth. Long bones and metapodial breakage categories include complete, proximal epiphysis (PE), proximal epiphysis and shaft (PES), shaft (S), shaft and distal epiphysis (SDE), and distal epiphysis (DE) portions. Scapulae, innominates, mandibles, crania, and ribs are classified into four to six breakage patterns, whereas patellae, carpals, tarsals, and teeth were only classified as complete or fragmented. Phalanges are recorded as complete, proximal, distal or fragmented, vertebrae are divided into complete, vertebral body, vertebral epiphysis or spinous processes. In the case of immature individuals, the long bone shafts preserving unfused epiphyses were considered complete elements (Rodríguez-Hidalgo et al., 2020).

Additionally, we documented the type of fracture (green, dry, modern or indeterminate) (Villa and Mahieu, 1991), shaft circumference (Bunn, 1982) and length, along with fracture outlines (Bunn, 1983; Villa and Mahieu, 1991). We also recorded bone cylinders, or “tubes,” which are shafts of long bones that preserve the complete circumference and have broken ends resulting from consumption. Following the criteria established by Rodríguez-Hidalgo et al. (2020), we considered only the humerus, tibia, and femur as potential “tubes” due to their significant marrow content, and distinguished them from “fake tubes,” long bone cylinders with sign of dry breakage (jagged or right-angle-smooth ends; see also Gabucio et al., 2024). We also identified notches according to the classifications proposed by Pickering and Egeland (2006). Helical and V-shaped fractures, typically indicative of carnivore breakage, were also recorded (Villa and Mahieu, 1991; Sanchis, 1999).

2.2.5 Bone surface modifications and digestion

We recorded the presence, number, location and distribution (isolated, dispersed or grouped) of pits, punctures, gnawing, digestion, helical fractures, and crenulated edges (Andrés et al., 2012; Binford, 1981; Capaldo and Blumenschine, 1994; Haynes, 1980; Lloveras et al., 2014; Sanchis et al., 2011). We also recorded measurements of the pits, punctures and scores using a microscope camera. We considered the possibility that these modifications were caused by hominins, although the criteria for differentiating them from marks made by other carnivores remain ambiguous (Fernandez-Jalvo and Andrews, 2016; Landt, 2007; Saladié et al., 2024). We compared the widths of pits and punctures on the cortical surfaces of long bone shafts with samples of marks made by small carnivores, based on studies by Andrés et al. (2012), Massigoge et al. (2014), Rodríguez-Hidalgo et al. (2013, 2015), and Saladié et al. (2024). In the comparative dataset, the animals consumed were primarily small, except for the human-made marks from Andrés et al. (2012). We included two samples of tooth marks made by lynxes, two by foxes, three by humans, one by hyena cubs, and one by wildcats (see Supplementary Table S5). We also carefully examined bones for cut marks (Pérez Ripoll, 1992; Binford, 1981; Blumenschine et al., 1996; Fernández Jalvo et al., 1999; Fisher, 1995) and several postdepositional modifications, including the presence of weathering, root etching, chemical corrosion, marks produced by bacteria colonies, manganese oxide staining, polishing, rounding, trampling, and microabrasion (Fernandez-Jalvo and Andrews, 2016). We also considered potential damage caused by fire, following stages described by Fernandez-Jalvo and Andrews (2016). Conspicuous marks were photographed using a microscope camera [Model: MicrosCopiaDigital XM Full HD Camera ALPHA B Series, resolution: 1080p (Full HD)].

2.2.6 Hierarchical Clustering on Principal Components (HCPC) and K-means cluster analysis on taphonomic data

We conducted a Hierarchical Clustering on Principal Components (HCPC) analysis on a combination of breakage and bone surface modification data, following a similar approach in Rodríguez-Hidalgo et al. (2020) and using the dataset shared by these authors, which includes reference assemblages generated by lynx, fox, eagles, owls, and humans (see Supplementary Table S3). Before performing the HCPC analysis, we handled missing data through imputation using the “missForest” package (Stekhoven and Buehlmann, 2012). The HCPC analysis was conducted using the “FactoMineR” package in R (Le et al., 2008), which performs PCA on the data, followed by hierarchical clustering on the principal components. We evaluated the results by calculating the silhouette scores and using the gap statistic method. Following this, we implemented the K-means clustering algorithm with six clusters. For this, we used the “cluster” (Maechler et al., 2023) and “factoextra” (Kassambara and Mundt, 2020) packages in R.

3 Results

3.1 Anatomical representation

A total of 1,189 European rabbit (Oryctolagus cuniculus) specimens (NISP) were identified to a specific skeletal element. Based on left calcanei, a total of 19 individuals (MNI) was estimated. Based on unfused distal tibiae and distal humeri, we estimated a minimum of three infantile individuals (younger than 2 or 3 months) and a minimum of six juveniles. Adult rabbit individuals therefore comprise 52.6% of the total MNI.

The total MNE is 775. Phalanges and metatarsals are the highest represented elements based on MNE, followed by isolated upper molars and metacarpals. Scapulae and ribs appear in much lower numbers (Table 2). The relative abundance (%RA) values show a predominance of calcanei and metatarsals (more than 60% of the RA), followed by tibiae, innominates, and cranial elements (up to 60% of the RA). The lowest represented elements in relative terms are axial remains like ribs and vertebrae (<5% of the RA; Figure 2).

Table 2

Table 2. Skeletal part frequencies of rabbit remains expressed in number of identified specimens (NISP), %NISP, minimum number of elements (MNE), %MNE, relative abundance (RA), %RA, and minimum animal units (MAU), %MAU.

Figure 2

Figure 2. Barplot showing the relative abundance (%RA) of rabbit skeletal elements.

The indices for the relative proportions of skeletal elements indicate a significant overrepresentation of postcranial elements compared to cranial elements (Table 3). This overrepresentation remains evident even when using indices that give considerable weight to the mineral density of skeletal elements, such as the TA/MD index. However, it is less pronounced when considering the MNEs (146.67; Table 3). Among limb bones, posterior (hind limb) elements are better represented than anterior (front limb) bones across the autopodium, zygopodium and stylopodium regions, and this relationship is consistent from proximal to distal elements (Table 3). The Spearman correlation coefficient of −0.02 (Figure 3) suggests a very weak negative correlation between relative abundance and bone density values. The p-value (0.93) indicates the correlation is not statistically significant. This implies that there is no meaningful relationship between skeletal representation and bone density, which suggests that density-mediated attrition does not appear to have biased the preservation of the assemblage.

Table 3

Table 3. Rabbit skeletal parts relative proportions.

Figure 3

Figure 3. Scatterplot illustrating the relationship between %RA and bone density values. The linear regression line shows no evident trend, supporting the lack of significant correlation.

3.2 Principal component analysis (PCA) and machine learning (ML) analysis of skeletal part abundances

The PCA analysis revealed a plot with considerable overlap between all predator groups. The variance explained by the principal components was relatively low, accounting for only 41.9% and 16% of the variance, respectively (Figure 4). This indicates high variability within the dataset. The PCA analysis also revealed a large distribution of samples across the principal component axes. While some samples cluster closely according to their respective groups (e.g. several non-ingested lynx samples at the top of the factor map and several owl assemblages underneath), indicating similarity in their anatomical representation characteristics, others are dispersed and located far from their group centers. This wide distribution suggests substantial variability within and between groups. A significant contributing factor to this high variability is the inclusion of both ingested (e.g., predator scats and pellets) and non-ingested samples in the dataset. These two types of assemblages often exhibit opposing characteristics, as the elements missing from the non-ingested samples are typically those that are digested, and vice versa. The groups from the diurnal raptors, nocturnal raptors, and terrestrial carnivores overlap significantly, making it impossible to distinguish these groups based solely on PCA. While the human sample forms a somewhat distinct cluster, it almost completely overlaps with the other groups. These results complicate a straightforward classification of the Escoural assemblage, which is positioned near the human cluster but also inside the raptors and terrestrial carnivores clusters. This could possibly be suggesting a complex taphonomic history involving multiple accumulating agents, however, it is impossible to determine whether the position of the sample indicates such complexity or simply reflects the ambiguity of the groups' distributions. The observed variability in the PCA results is to be expected given the diverse nature of the dataset and the fact that predators generally exhibit variable behaviors. It is also possible that the PCA is insufficient to discern clear patterns of difference between groups, and machine learning techniques might be better suited to uncover patterns.

Figure 4

Figure 4. PCA biplot showing the distribution of samples across the two principal components, with samples colored according to their respective groups. Confidence ellipses illustrate the variability within groups. See similar plot with the labels of the assemblages in Supplementary Figure S1.

Following the results of feature selection and variable correlation analysis, we excluded the skeletal abundance values for cranium, humerus, and metacarpals (see Supplementary Figure S2 and Supplementary Table S4). Tables 4, 5 summarize the results of the ML models. The RF model showed moderate performance and balanced outcomes in precision and recall across most groups. However, it struggled with identifying terrestrial carnivores. Compared to RF, the SVM underperformed, showing lower accuracy and kappa values, and was not good in distinguishing between classes, especially nocturnal and diurnal raptors. The best results were obtained using GBM. This model showed better performance in identifying humans and terrestrial carnivores and yielded higher accuracy, kappa, precision, recall and F1-scores across the different predator groups (Tables 4, 5).

Table 4

Table 4. Selected hyperparameters and performance metrics for random forest (RF), support vector machine (SVM), and gradient boosting machine (GBM) models, including accuracy, 95% confidence intervals (CI), and Kappa values.

Table 5

Table 5. Sensitivity (recall), specificity, positive predictive value, negative predictive value, balanced accuracy, precision, and F1-score for each predator group across the random forest (RF), support vector machine (SVM), and gradient boosting machine (GBM) models.

The classification results of the GBM indicated terrestrial carnivores as the primary accumulating agents of the Escoural sample with a probability of 0.99. In contrast, the RF model classified Escoural as a human-made accumulation, while the SVM model also grouped Escoural with terrestrial carnivores. The high confidence of the GBM model in classifying the Escoural sample, despite the variability of the dataset, raises concerns about overfitting. However, given the comprehensive cross-validation and the model's consistent performance (see Supplementary Figure S3), the prediction appears reliable. However, it is important to interpret these results within the broader context of taphonomic and archaeological evidence. Escoural could have been influenced by multiple agents, in this case, several terrestrial carnivores and humans, and this classification should be viewed as identifying the predominant agent rather than excluding other possible contributors.

3.3 Breakage

Approximately 26.8% of the rabbit specimens are complete. Table 6 details the preservation status of the different bone portions. Smaller bones, such as phalanges, metapodials, and carpals/tarsals, are typically found intact. In contrast, longer bones, including limb bones, girdles, and cranial elements, were predominantly recovered in a broken state. The size of the rabbit remains ranges from 3.98 to 60.36 mm, with an average length of 14.23 mm, and over 85% of the specimens exceeding 10 mm in length (Figure 5). Upper and intermediate limb bones were highly fragmented, with shafts being the most represented bone portion. Specifically, shafts accounted for 48% of radii-ulnae specimens, 66% of humeri, 76% of tibiae, and 80% of femora. Conversely, lower limb bones were often found complete, with 70% of metacarpals and 22% of metatarsals recovered intact. With the exception of the ulnae, which preserved the proximal end significantly more than the distal end, the remaining upper and intermediate long bones showed somewhat varied preservation rates between their proximal and distal ends. For instance, the humeri preserved 14% of their proximal ends compared to 18% of their distal ends, with similar variations observed in other long bones. The radii and femora exhibited more preservation of the proximal ends, while the tibia had greater distal preservation (Table 6). Additionally, metapodials preserved their proximal ends around 36% of the time and their distal ends 14%. Cranial portions were primarily represented by M and NC fragments, while mandibles were mostly represented by MBI portions. No complete crania or mandibles were recovered. Innominates were mainly represented by fragments retaining the acetabulum (AIL, AIS, AISIL), scapulae retained the GC portion, and ribs were primarily found as shaft fragments (S). Phalanges, carpals, tarsals, and isolated teeth were mostly complete (Table 6).

Table 6

Table 6. Representation of different breakage categories.

Figure 5

Figure 5. Specimen length distribution of the recovered rabbit remains.

Green fractures were recorded on 34% of the rabbit assemblage, while dry fractures were identified on ~6% of the specimens. Additionally, 34% of the specimens exhibited indeterminate breakage patterns, and another 6% had recent fractures. Table 7 summarizes the frequencies of green breakages on the different portions of marrow-rich long bones (humerus, femur and tibia) and the number of green broken shaft cylinders or “tubes.” In total, around 70% of the marrow-rich long bone sample presented green breakages. Cylinders occurred most often on tibiae. Fifteen marrow-rich long limb bones (one humerus, five femora, and nine tibiae) retained the epiphysis and part of the shaft presenting a helical or V-shaped fracture, which has been described as a morphotype of carnivore breakage (Rodríguez-Hidalgo et al., 2020) (Figure 6). While shaft fragments were more abundant, epiphyses were still fairly represented across all long bones. This pattern suggests that felid activity may have played a more significant role in the assemblage compared to canids, as felids tend to leave more epiphyses unconsumed. Only 9% of the marrow-rich long bone specimens presented clear dry breakages, which suggests limited postdepositional damage and an overall high integrity of the rabbit assemblage. A total of 22 notches were recorded on marrow-rich long bones. Most notches were “complete or type A” (60%), although there were some micro-notches (type E) (36%) and overlapping notches (type C).

Table 7

Table 7. Counts and frequencies (parentheses) of green breakage across the different marrow rich long bone portions and tube counts.

Figure 6

Figure 6. Examples of typically preserved marrow-rich long bone breakage portions and tubes in the Escoural rabbit sample. (A) examples of humeri specimens, represented mostly by distal epiphysis plus a small portion of the shaft, (B) femora specimens generally represented by proximal epiphysis plus shaft, (C) tibiae specimens preserving distal epiphysis plus shaft, (D) examples of tubes or cylindrical femora and tibia shafts.

3.4 Bone surface modifications and digestion

Approximately 75% of the bone cortical surfaces exhibited good or moderate preservation, indicating an overall high preservation state of the assemblage. Fissures resulting from changes in humidity and concretions, typical modifications associated with karstic environments, were present in 2% of the assemblage. Weathering on the rabbit bones was minimal, with 98% of the bones classified as stage 0, indicating no visible weathering, and the remaining 2% at stage 1, which is characterized by slight splitting of bone parallel to fiber structure (Andrews, 1990). These weathering stages suggest that the remains were exposed for 0–2 years (stage 0) or 1–5 years (stage 1) in a wet temperate climate (Andrews, 1990). This further supports the overall high preservation and limited exposure of the assemblage. Traces of trampling and microabrasion were observed also on ~2% of the remains. Biochemical marks were present on 8% of the specimens, while root etching occurred on 23%. Manganese staining affected 60% of the bones, with 35% of these showing generalized or complete coverage. Rounding was observed in 16% of the bone sample, with 15% classified as R1 (rounding affects the anatomical edges and fractures at the microscopic level) and only 1% as R2 (rounding affects certain areas of the bone and can be observed at the macroscopic level). Polished surfaces or edges were found on <5% of the specimens. Signs of burning were present on only two bone fragments, characterized by a black, homogeneous coloration on the bone surface. Rodent gnawing was observed on 11 rabbit bone specimens.

Carnivore damage in the form of tooth marks was observed on 96 specimens, accounting for 8% of the NISP of the rabbit bone assemblage. Digestive damage affected 9% of the sample. Specific types of damage include pits on 55 specimens (4.6%), punctures on 27 specimens (2.3%), scoring on 35 specimens (2.9%), and crenulated edges on 85 specimens (7.2%). Tooth marks appeared mainly on long bones and innominates. Metatarsals and tibiae are the most frequently tooth-marked elements (18.8 and 17.7% respectively), followed by femora (13.5%) and innominates (11.5%). Table 8 shows the counts and percentages of tooth-marked and digested specimens across elements. Hind limbs (femora, tibiae, metatarsals, and innominates) exhibit a significantly higher amount of carnivore damage in absolute terms. However, in relative terms, ulnae and innominates stand out as the elements with the highest percentages of tooth marks. Innominates also show significantly higher percentages of punctures and scores compared to other elements (Figure 7). Figure 8 provides examples of damage inflicted by carnivores on the rabbit bone sample. All three types of tooth marks (pits, punctures, and scores) occur in relatively similar percentages across long bones.

Table 8

Table 8. Counts and percentages per NISP of tooth-marked and digested bone specimens per element.

Figure 7

Figure 7. Heatmap of types of carnivore damage across skeletal elements. TM refers to the sum of the remains with any kind of tooth mark modification (pits, punctures and scores).

Figure 8

Figure 8. Examples of tooth marked and digested specimens. (A) Specimens showing pits and crenulated edges, (B) pelvis fragment with carnivore scores, (C) examples of punctured specimens, (D) several digested calcanei.

The comparison between tooth mark sizes from Escoural and those documented for various terrestrial carnivores and humans indicates that the marks fall within the size ranges associated with foxes, wildcats, and humans (Figure 9). Indeed, mean tooth mark size is relatively small and several v-shaped scores, typical of cubs, were identified. The upper range of the 95% confidence interval for the Escoural sample reaches the range of one of the lynx samples. This analysis indicates that rabbits were likely consumed by small carnivores, although the overlap makes it unclear whether multiple predators were involved. Lynx and especially lynx cubs should not be completely ruled out as contributors, and humans could have been involved as well. Despite these observations, the variability within groups, as well as the overlap between carnivore types, makes it challenging to draw definitive conclusions about the specific predator involved. This analysis alone is insufficient to clearly identify the accumulating agent.

Figure 9

Figure 9. Mean values and 95% confidence intervals of tooth pit and puncture widths on cortical surfaces of long bone shafts on bone assemblages modified by fox, humans, hyena cubs, lynxes and wildcats. Comparative samples from Andrés et al. (2012), Massigoge et al. (2014), Rodríguez-Hidalgo et al. (2013, 2015), and Saladié et al. (2024). The values can be found in Supplementary Table 4.

Regarding digestion, among all elements, tarsals exhibited the highest amount of digestion (41%). This was followed by ulnae at 20%, humeri at 12%, and phalanges at 8% (Table 8). Approximately 66% of the digested remains were classified as light, around 27% as moderate and 7% as heavy. Extreme digestive damage was not documented. Anthropogenic marks in the form of cut marks were not observed. Only two burned remains (a metatarsal and an indeterminate long bone shaft) were observed.

3.5 Hierarchical Clustering on Principal Components (HCPC) and K-means cluster analysis on taphonomic data

The HCPC analysis reveals several clusters based on the selected taphonomic variables that more or less align with the known predator groups. The discrepancies could be due to high variability and overlap in taphonomic modification patterns within and between the predator groups, or to limitations in the dataset and the resolution of the method used. The first two principal components explain 63.6% of the variance. Clusters 1, 2, and 3 include assemblages where humans were involved (Figure 10; see also Supplementary Figure S4). Cluster 4 contains assemblages of mixed origins, including one with human modifications as well as others associated with fox and lynx activity. Assemblages with direct anthropogenic evidence are clearly separated along the first dimension from the remaining clusters, which include only carnivore-made assemblages. Cluster 5 contains a mixture of carnivore-made assemblages. The Escoural assemblage clustered with an assemblage made by fox and one made by wildcat (Álvarez et al., 2012), and lies very close to Navalmaíllo Rockshelter, an assemblage interpreted as generated by lynxes (Arriaza et al., 2017). It also falls near to other assemblages made by lynxes and foxes, further supporting the interpretation that the assemblage may have mostly originated from terrestrial carnivore activity (Figure 10). The average silhouette width showed that clusters had varying quality (see Supplementary Figure S5a). Whereas clusters 2, 3, and 7 seemed to be well-defined, clusters 1, 4, 5, and 6 had some misclassified samples and moderate quality (as shown in Supplementary Figure S5a and Supplementary Table S6a).

Figure 10

Figure 10. HCPC cluster plot of taphonomic data.

The K-means clustering analysis with six clusters yielded a slightly higher average silhouette width of 0.37, suggesting the analysis is more reliable (Supplementary Figure S5b and Supplementary Table S6b). In this case, the clusters showed similar alignment with the known groups. For example, cluster 1 seemed to be mostly related with non-ingested assemblages, cluster 3 with ingested samples, and clusters 5 and 6 contained only assemblages with human input (Figure 11). In fact, the coordinates of the centroids for each cluster, which describe the main characteristics of each cluster, show that an important feature in describing cluster 3 are indeed digested remains, and clusters 2, 3, and 4 are characterized by the amount of adult individuals (Table 9). Interestingly, Escoural showed taphonomic similarities to Cova del Gegant IIIa, which has been interpreted as an assemblage made primarily by lynxes (Rodríguez-Hidalgo, 2022), and Navalmaíllo, which has a similar interpretation (Arriaza et al., 2017). This cluster also contains an eagle nest assemblage and an Layer III at Estebanvela, which has a higher proportion of evidence for non-human accumulation agents, such as carnivores and raptors than Estebanvela I and II (Yravedra et al., 2019). This groups seems to be defined best by the proportion of adult individuals and tubes or cylinders (Table 9). Clusters 1 and 2 lie close and are also mainly characterized by the amount of adults and complete long bones, respectively. Cut marks, burned remains and tooth marks do not appear to be very important in this classification.

Figure 11

Figure 11. Factor map of the K means clustering analysis of taphonomic data.

Table 9

Table 9. Centroid values for each cluster obtained from k-means clustering, showing the mean values of each taphonomic feature within each cluster.

4 Discussion

Faunal studies have revealed that Neanderthals displayed complex and flexible subsistence strategies and variability in site use patterns, and that rabbit hunting and consumption were sometimes part of them (Cochard et al., 2012; Carvalho et al., 2018; Pérez et al., 2017; Rufà et al., 2014; Pelletier et al., 2019). At Escoural Cave, clear evidence for human consumption of rabbits in the form of cut marks is absent. Only two burned specimens were recorded, suggesting that cremation occurred most likely post-depositionally, given the small number of remains (Pérez et al., 2017; Rodríguez-Hidalgo et al., 2020; Téllez et al., 2022). This indicates that Neanderthals may have focused on other food sources or different subsistence strategies during their occupation of the site.

The results suggest that small and medium-sized terrestrial carnivores, primarily lynxes, and possibly also red foxes and wildcats, played the main role in the accumulation and modification of the remains. These three carnivores are represented in the faunal assemblage from Escoural, analyzed in previous studies (Otte and da Silva, 1996), and Iberian lynx and a smaller carnivore have been documented in the remaining assemblage recovered in 2020 from Profile 1. Substantial carnivore activity is reflected in the frequencies of carnivore damage and digested elements. The anatomical representation showed a predominance of hind limb elements like tibiae and metatarsals. Phalanges and tarsals were also highly represented. This pattern is consistent with accumulations by Iberian lynxes and by red foxes, which show a high prevalence of distal limb elements and often the preservation of complete autopodial bones (Rodríguez-Hidalgo et al., 2013, 2015, 2020; Rodríguez-Hidalgo, 2022). Similarly, wildcat accumulations show a high representation of cranial elements, metatarsals, and femora (Lloveras et al., 2016). In contrast, diurnal and nocturnal raptors exhibit a different pattern, with a higher proportion of proximal upper limbs (humerus, femur) (Rodríguez-Hidalgo et al., 2013).

The high confidence of the GBM model in classifying the Escoural sample as predominantly created by terrestrial carnivores supports this interpretation. However, pinpointing the exact predator species remains challenging. Mortality profiles could potentially help in this differentiation, as foxes and wildcats generally prey on juvenile rabbits due to their smaller size and relative vulnerability, which results in age profiles skewed toward younger individuals. Lynxes, in turn, as well as larger predators in general, including humans, tend to prey on prime adults (Pelletier et al., 2019). The Escoural sample is dominated by adults, which would suggest lynx as the primary accumulators. Additionally, the presence of infant rabbits in the assemblage could indicate individuals that likely died within their burrows, as infant rabbits typically remain in the burrow until they reach the juvenile stage.

Regarding breakage patterns, high rates of green fractures and minimal dry fractures suggested limited postdepositional damage and high assemblage integrity, and typical carnivore breakage patterns, including helical fractures, were predominant. The analyzed sample also showed extensive breakage of long bones, which has been typically observed in lynx-made rabbit assemblages. However, red foxes and wildcats and humans can also cause significant breakage, often associated with gnawing and scavenging behaviors. Notably, foxes tend to leave a higher number of tooth marks on skeletal elements compared to lynxes (Rodríguez-Hidalgo et al., 2013). This characteristic aligns with the findings from Escoural, which exhibits a relatively high frequency of tooth marks compared to other assemblages attributed to lynxes, such as Cova del Gegant IIIa and Navalmaíllo (Rodríguez-Hidalgo et al., 2020; Arriaza et al., 2017). The comparison of tooth mark sizes indicates that the damage observed in the Escoural assemblage could be consistent with various agents, including small carnivores and potentially humans, but the overlap and variability in tooth mark sizes prevent a definitive identification of the contributors. Moreover, there is no additional evidence of anthropogenic modification on the bones to support human involvement conclusively.

However, equifinality in rabbit assemblages complicates the identification of human activities. Studies have shown that even when cut marks and burning damage, which are exclusively attributed to humans, are considered, detecting human involvement as predators of rabbit assemblages remains challenging. Most experiments aimed at describing the human taphonomic signature on rabbit remains have focused on exclusive products of human activity, including butchering processes and cooking, which leave cut marks and burned remains (Lloveras et al., 2009; Sanchis et al., 2011). However, these modifications are variable and sometimes even absent in archaeological assemblages (Hockett and Bicho, 2000; Hockett and Haws, 2002; Ripoll, 2004). Given the fact that rabbits can be processed without the use of tools, additional taphonomic features can be caused by human activity, namely tooth marks, shafts with notches on fractured edges or cylinders resulting from marrow extraction (Saladié et al., 2024; Hockett and Haws, 2002; Cochard et al., 2012). The relatively low proportion of cylinders observed at Escoural aligns with patterns typically associated with lynx activity, while humans generally produce a much more extensive pattern of cylinders along with the disappearance of epiphyses. Additionally, the presence of tooth marks on metapodials and phalanges, which are not commonly found in human accumulations as these elements are usually not consumed and often remain articulated, further supports the likelihood of carnivore activity at the site (Yravedra et al., 2019). In the case of Escoural, no metapodials were found with articulated phalanges, which may provide additional evidence favoring carnivore involvement over human agency. Yet, some taphonomic evidence at Escoural suggest that Neanderthal activity could potentially be underrepresented or masked, like the predominance of adult individuals, extensive breakage or some tooth mark sizes. These and other typical characteristics of Neanderthal consumption of rabbit, which include a relatively balanced representation of cranial and postcranial remains and a significant presence of metapodials, scapulae and humeri, have been observed in several rabbit accumulations with Neanderthal input, including Bolomor Cave (Blasco and Fernández Peris, 2012), Arbreda Cave (Lloveras et al., 2010), Les Canalettes (Cochard et al., 2012), Columbeira Cave (Carvalho et al., 2018), El Salt (Pérez et al., 2017), or Teixoneres (Rufà et al., 2014). As Saladié et al. (2024) have pointed out, misinterpreting human marks could be as problematic as overestimating carnivore involvement in the formation of assemblages. Given the lack of clear evidence, we cannot argue that humans contributed to the Escoural assemblage, but we also cannot entirely exclude the possibility of some human involvement until more data and comprehensive analyses are available.

Taphonomic attributes left by the same predator species exhibit significant variability. The same type of predator, including humans, can create bone accumulations with particular characteristics depending on how the site is used, e.g. as a consumption area, defecation and abandonment zone, breeding den, transitional area, or prey acquisition zone (Cochard, 2004). At the same time, very similar modification patterns can be found on scats, pellets or non-ingested remains in nests and dens of different predators. This means that this variability must be considered when accurately characterizing archaeological assemblages. On the one hand, expanding frames of reference is crucial to address the challenge of identifying human activity in rabbit assemblages beyond cut marks and burning, as well as to distinguish between other non-human predators. On the other hand, several authors have pointed out that only a comprehensive approach combining multiple variables can aid in identifying the agent(s) responsible for rabbit archaeological assemblages (Lloveras et al., 2012a,b; Rodríguez-Hidalgo et al., 2020). Thus, taphonomic studies have introduced the use of multivariate analyses that consider multiple criteria, such as anatomical representation, breakage patterns, and the frequency and type of surface modifications to classify sites among a set of archaeological and reference assemblages to help reduce the effects of equifinality (e.g. Arriaza et al., 2017; Rodríguez-Hidalgo et al., 2020). However, since assemblages can have a complex history involving multiple accumulating agents and processes, studies have to consider that finding a single perfect prediction may not accurately reflect the true nature of the assemblage. For instance, the Escoural sample is composed of both ingested and non-ingested remains and very likely has contributions from various predators. While these analyses can provide valuable insights, they must be interpreted in the context of the potential for multiple agents and processes influencing the assemblage.

Similarly, the application of machine learning models presents its own set of challenges. Although they can handle high-dimensional data and complex interactions between variables more effectively and provide a probabilistic framework that accounts for the possibility of multiple accumulating agents, the inherent variability in predator behavior, the limitations of current reference frameworks and comparative collections, and the fact that datasets are generally small, can impact the accuracy and reliability of the predictions. Thus, while machine learning offers powerful tools for taphonomic analysis, its results should be integrated with traditional taphonomic evidence and expert interpretation to form a holistic understanding of the assemblage formation. Nevertheless, the use of ML models, can reduce the subjectivity of manual classification, as well as provide a more standardized approach in identifying taphonomic agents. Even on small datasets, models can establish a framework that can be applied to larger datasets in the future. As more data becomes available, the models can be retrained and improved so that their predictive power is increased.

This study provides an example of how multiple analytical methods (PCA, HCPC, ML) can contribute to interpreting complex rabbit assemblages by complementing traditional taphonomic analysis. Additionally, it expands the reference frameworks for carnivore fossil accumulations, which can be useful for future studies.

Our findings also offer insights into Neanderthal behavioral complexity. The absence of evidence for rabbit consumption by Neanderthals at Escoural suggests a context-dependent variability in their subsistence strategies, supporting the idea that Neanderthals were highly adaptable, selecting alternative food resources or exploiting different ecological niches according to specific environmental conditions. Moreover, the lack of anthropic consumption indicates that human occupations were likely of short duration, allowing carnivores to occupy the site during periods of human absence. This suggests potential behavioral influences between Neanderthals and carnivores, affecting the duration of occupations, settlement choices, mobility patterns, and other aspects of subsistence strategies (Picin et al., 2020; Daujeard et al., 2012; Linares-Matás and Yravedra, 2024; Pinto-Llona et al., 2023). Future research will include the analysis of additional faunal remains from both previous and recent excavations and will explore whether humans focused on other prey types at Escoural. Ongoing dating and microstratigraphic analyses will also help refine the site's occupation history.

5 Conclusions

The rabbit assemblage from Escoural Cave resulted from the activity and accumulation of a mixture of terrestrial carnivores, among which the Iberian lynx stands out. The classification of the Escoural sample as predominantly influenced by terrestrial carnivores is supported by consistent taphonomic evidence across all the analyses included here, with the lynx identified as the primary agent. However, the variability and overlap observed in the PCA and HCPC analyses between different small and medium-sized terrestrial carnivore species, indicate that multiple agents may have contributed to the assemblage. The results also highlight the challenges of distinguishing between different accumulation agents in highly variable datasets, but underscore the potential of advanced statistical and ML techniques to enhance taphonomic analyses.

Complexity in subsistence behavior includes the ability to adapt to different environmental and ecological conditions. At Escoural, even though Neanderthals were using the cave or its surroundings, they were not consuming the rabbits accumulated there. This shows that Neanderthal subsistence behavior was variable and dependent on the context. While Neanderthals are known to have consumed small game in other contexts, the absence of rabbit consumption in this assemblage suggests they may have had alternative food sources available during their occupation of the cave.

Future studies will further explore whether the cave could have functioned mainly as a carnivore den, and have been occupied by Neanderthals only intermittently. The cave has yielded a very diverse faunal spectrum, which includes multiple carnivore species, suggesting varied interactions between different species and a dynamic use of the site over time. Thus, Escoural holds significant importance in the southern Portugal Mousterian context. Research on Middle Paleolithic occupation in central Portugal has provided significant insights into Neanderthal behavior and adaptations. Key sites such as Figueira Brava, Gruta da Oliviera, Gruta da Companheira, Foz do Enxarrique, and Almonda Cave have revealed a complex picture of Neanderthal life, yielding evidence of diverse subsistence strategies, including the hunting of large game and the exploitation of aquatic resources. Despite these advances, more comprehensive studies and excavations using modern analytical techniques are still needed. Continued research at Escoural will soon contribute to a deeper understanding of Neanderthal behavior and the interplay between human and carnivore activities.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

Author contributions

LC-S: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing. AR: Formal analysis, Methodology, Writing – review & editing. JC: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was funded by the European Union (ERC-CoG-2021, FINISTERRA, 101045506). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. AR is a beneficiary of the Individual Call to Scientific Employment Stimulus—3rd Edition, promoted by the Portuguese FCT (reference: 2020.00877.CEECIND). JC was funded by Fundação para a Ciência e a Tecnologia, contract reference DL57/2016/CP1361/CT0026. LC-S is funded by the European Union's Horizon Europe programme under the Marie Skłodowska-Curie Postdoctoral Fellowship grant agreement No. 101152206.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fearc.2024.1473266/full#supplementary-material

References

Álvarez, M. C., Kaufmann, C. A., Massigoge, A., Gutiérrez, M. A., Rafuse, D. J., Scheifler, N. A., et al. (2012). Bone modification and destruction patterns of leporid carcasses by Geoffroy's cat (Leopardus geoffroyi): an experimental study. Quat. Int. 278, 71–80. doi: 10.1016/j.quaint.2011.11.029

Crossref Full Text | Google Scholar

Andrés, M., Gidna, A. O., Yravedra, J., and Domínguez-Rodrigo, M. (2012). A study of dimensional differences of tooth marks (pits and scores) on bones modified by small and large carnivores. Archaeol. Anthropol. Sci. 4, 209–219. doi: 10.1007/s12520-012-0093-4

Crossref Full Text | Google Scholar

Andrews, P. (1990). Owls, Caves and Fossils. London: Natural History Museum Publications.

PubMed Abstract | Google Scholar

Armstrong, A. (2016). Eagles, owls, and coyotes (oh my!): taphonomic analysis of rabbits and guinea pigs fed to captive raptors and coyotes. J. Archaeol. Sci. Rep. 5, 135–155 doi: 10.1016/j.jasrep.2015.10.039

Crossref Full Text | Google Scholar

Arriaza, M. D. C., Huguet, R., Laplana, C., Pérez-González, A., Márquez, B., Arsuaga, J. L., et al. (2017). Lagomorph predation represented in a middle Palaeolithic level of the Navalmaíllo Rock Shelter site (Pinilla del Valle, Spain), as inferred via a new use of classical taphonomic criteria. Quatern. Int. 436, 294–306. doi: 10.1016/j.quaint.2015.03.040

Crossref Full Text | Google Scholar

Auguste, P. (1995). Chasse et charognage au Paléolithique moyen: l'apport du gisement de BiacheSaint-Vaast (Pas-de-Calais). Bull. Soc. Prehist. Fr. 92, 155–168. doi: 10.3406/bspf.1995.10003

Crossref Full Text | Google Scholar

Barone, R., Paavaux, C., Blin, P. C., and Cuq, P. (1973). Atlas of Paris: Rabbit Anatomy. Issy les Moulineaux Cedex: Masson and Co. Editeurs.

Google Scholar

Binford, L. R. (1978). Nunamiut Ethnoarchaeology. Clinton Corners, NY: Eliot Werner Publications/Percheron Press.

Google Scholar

Binford, L. R. (1981). Bones: Ancient Men and Modern Myths. Cambridge, MA: Academic Press.

Google Scholar

Binford, L. R. (1984). Butchering, sharing, and the archaeological record. J. Anthropol. Archaeol. 3, 235–257. doi: 10.1016/0278-4165(84)90003-5

PubMed Abstract | Crossref Full Text | Google Scholar

Blasco, R., Cochard, D., Colonese, A. C., Laroulandie, V., Meier, J., Morin, E., et al. (2022). “Small animal use by Neanderthals,” in Updating Neanderthals (Academic Press), 123–143. doi: 10.1016/B978-0-12-821428-2.00010-X

Crossref Full Text | Google Scholar

Blasco, R., and Fernández Peris, J. (2012). A uniquely broad spectrum diet during the Middle Pleistocene at Bolomor Cave (Valencia, Spain). Quat. Int. 252, 16–31. doi: 10.1016/j.quaint.2011.03.019

Crossref Full Text | Google Scholar

Blumenschine, R. J., Marean, C. W., and Capaldo, S. D. (1996). Blind test of interanalyst correspondence and accuracy in the identification of cut marks. Percussion marks and carnivore tooth marks on bone surface. J. Archaeol. Sci. 23, 493–505. doi: 10.1006/jasc.1996.0047

Crossref Full Text | Google Scholar

Brown, K., Fa, D. A., Finlayson, G., and Finlayson, C. (2011). “Small game and marine resource exploitation by Neanderthals: the evidence from Gibraltar,” in Trekking the Shore: Changing Coastlines and the Antiquity of Coastal Settlement, ed. N. F. Bicho, J. A. Haws, and L. G. Davis (Cham: Springer), 247–272. doi: 10.1007/978-1-4419-8219-3_10

Crossref Full Text | Google Scholar

Buccheri, F., Bertè, D. F., Berruti, G. L. F., Cáceres, I., Volpe, L., and y Arzarello, M. (2016). Taphonomic analysis on fossil remains from the Ciota Ciara Cave (Piedmont, Italy) and new evidence of cave bear and wolf exploitation with simple quartz flakes by Neanderthal. Riv. Ital. Paleontol. Stratigr. 122, 41–54. doi: 10.13130/2039-4942/7674

Crossref Full Text | Google Scholar

Bunn, H. T. (1982). Meat-eating and human evolution: Studies in the diet and subsistence patterns of Plio-Pleistocene hominids in East Africa (Ph.D. Dissertation). Berkeley, CA: University of California, Berkeley.

Google Scholar

Bunn, H. T. (1983). Comparative analysis of modern bone assemblages from a San hunter-gatherer camp in the Kalahari Desert, Botsuwana, and from a spotted hyena den near Nairobi, Kenya. Anim. Archaeol. 163, 141–148.

Google Scholar

Callou, C. (1997). “Diagnose différntielle des principaux éléments squelettiques du lapin (Genre Oryctolagus) et du Lèvres (Genre Lepus) en Europe Occidentale,” in Fiches d'ostéologie animale pour l'archéologie. Series B: Mammifères. Centre de recherches Archéologiques du CNRS, eds. N. Desse-Berset, and J. Desse (Valbonne-Sophia Antipolis), 8–21.

Google Scholar

Callou, C. (2003). De la garenne au clapier: Étude archéozoologique du Lapin en Europe occidentale. Mus. Natl. Hist. Nat. 189:356.

Google Scholar

Camarós, E., Cueto, M., Rosell, J., Díez, J. C., Blasco, R., Duhig, C., et al. (2017). Hunted or scavenged Neanderthals? Taphonomic approach to hominin fossils with carnivore damage. Int. J. Osteoarchaeol. 27, 606–620. doi: 10.1002/oa.2584

Crossref Full Text | Google Scholar

Capaldo, S. D., and Blumenschine, R. J. (1994). A quantitative diagnosis of notches made by hammerstone percussion and carnivore gnawing on bovid long bones. Am. Antiq. 59, 724–748. doi: 10.2307/282345

Crossref Full Text | Google Scholar

Carvalho, M., Pereira, T., and Manso, C. (2018). Rabbit exploitation in the middle Paleolithic at Gruta Nova da Columbeira, Portugal. J. Archaeol. Sci. Rep. 21, 821–832. doi: 10.1016/j.jasrep.2018.09.003

Crossref Full Text | Google Scholar

Cochard, D. (2004). Etude taphonomique des léporidés d'une tanière de renard actuelle: apport d'un référentiel à la reconnaissance des accumulations anthropiques. Rev. Paléobiol. 23, 659−673.

Google Scholar

Cochard, D., Brugal, J.-P., Morin, E., and Meignen, L. (2012). Evidence of small fast game exploitation in the Middle Paleolithic of Les Canalettes Aveyron, France. Quat. Int. 264, 32–51. doi: 10.1016/j.quaint.2012.02.014

Crossref Full Text | Google Scholar

Cohen, B. F., and Kibii, J. M. (2015). Experimental study of bone modification by captive caracal (Caracal caracal); a model for fossil assemblage analysis. J. Archaeol. Sci. 62, 92–104. doi: 10.1016/j.jas.2015.07.008

Crossref Full Text | Google Scholar

Daujeard, C., Fernandes, P., Guadelli, J. L., Moncel, M. H., Santagata, C., Raynal, J. P., et al. (2012). Neanderthal subsistence strategies in Southeastern France between the plains of the Rhone Valley and the mid-mountains of the Massif Central (MIS 7 to MIS 3). Quat. Int. 252, 32–47. doi: 10.1016/j.quaint.2011.01.047

Crossref Full Text | Google Scholar

Daujeard, C., and Moncel, M.-H. (2010). On neanderthal subsistence strategies and land use: a regional focus on the Rhone Valley area in southeastern France. J. Anthropol. Archaeol. 29, 368–391. doi: 10.1016/j.jaa.2010.05.002

Crossref Full Text | Google Scholar

Daujeard, C., Vettese, D., Britton, K., Béarez, P., Boulbes, N., Crégut-Bonnoure, E., et al. (2019). Neanderthal selective hunting of reindeer? The case study of Abri du Maras (south-eastern France). Archaeol. Anthropol. Sci. 11, 985–1011. doi: 10.1007/s12520-017-0580-8

Crossref Full Text | Google Scholar

Delibes, M., and Hiraldo, F. (1981). “The rabbit as prey in the Iberian Mediterranean ecosystem,” in Proceedings of the world lagomorph conference, Vol. 1979 (Guelph: University of Guelph), 614–622.

Google Scholar

Dodson, P., and Wexlar, D. (1979). Taphonomic investigations of owl pellets. Paleobiology 5, 275–284. doi: 10.1017/S0094837300006564

Crossref Full Text | Google Scholar

Duveau, J., Berillon, G., Verna, C., Laisné, G., and Cliquet, D. (2019). The composition of a Neandertal social group revealed by the hominin footprints at Le Rozel (Normandy, France). Proc. Nat. Acad. Sci. 116, 19409–19414. doi: 10.1073/pnas.1901789116

PubMed Abstract | Crossref Full Text | Google Scholar

Fernández Jalvo, Y., Andrews, P., and Denys, S. C. (1999). Cut marks on small mammals at Olduvay Gorge Bed I. J. Hum. Evol. 36, 587–589. doi: 10.1006/jhev.1999.0292

PubMed Abstract | Crossref Full Text | Google Scholar

Fernandez-Jalvo, Y., and Andrews, P. (2016). Atlas of Taphonomic Identifications: 1001+ Images of Fossil and Recent Mammal Bone Modification. Cham: Springer. doi: 10.1007/978-94-017-7432-1

Crossref Full Text | Google Scholar

Fisher, D. C. (1995). Bone surface modifications in zooarchaeology. J. Archaeol. Method Theory 2, 7–68. doi: 10.1007/BF02228434

Crossref Full Text | Google Scholar

Gabucio, M. J., Bargalló, A., Saladié, P., Romagnoli, F., Chacón, M. G., Vallverdú, J., et al. (2023). Using GIS and geostatistical techniques to identify Neanderthal Campsites at archaeolevel Ob at Abric Romaní. Archaeol. Anthropol. Sci. 15:24. doi: 10.1007/s12520-023-01715-6

Crossref Full Text | Google Scholar

Gabucio, M. J., Cáceres, I., Rodríguez-Hidalgo, A., Rosell, J., and Saladié, P. (2014). A wildcat (Felis silvestris) butchered by Neanderthals in Level O of the Abric Romaní site (Capellades, Barcelona, Spain). Quat. Int. 326–327, 307–318. doi: 10.1016/j.quaint.2013.10.051

Crossref Full Text | Google Scholar

Gabucio, M. J., Sanz, M., and Daura, J. (2024). Lynxes and foxes among hyenas, bears and humans. Unit 3 of the Cova del coll Verdaguer Upper Pleistocene site (Barcelona, Iberian Peninsula). Quat. Sci. Rev. 333:108671. doi: 10.1016/j.quascirev.2024.108671

Crossref Full Text | Google Scholar

Gaudzinski-Windheuser, S., and Kindler, L. (2012). Research perspectives for the study of Neandertal subsistence strategies based on the analysis of archaeozoological assemblages. Quat. Int. 247, 59–68. doi: 10.1016/j.quaint.2010.11.029

Crossref Full Text | Google Scholar

Gómez-Olivencia, A., Sala, N., Núñez-Lahuerta, C., Sanchis, A., Arlegi, M., Rios-Garaizar, J., et al. (2018). First data of Neandertal bird and carnivore exploitation in the Cantabrian Region (Axlor; Barandiaran excavations; Dima, Biscay, Northern Iberian Peninsula). Sci. Rep. 8:10551. doi: 10.1038/s41598-018-28377-y

PubMed Abstract | Crossref Full Text | Google Scholar

Gutiérrez-Zugasti, I., Rios-Garaizar, J., Marín-Arroyo, A. B., Rasines del Río, P., Maroto, J., Jones, J. R., et al. (2018). A chrono-cultural reassessment of the levels VI–XIV from El Cuco rock-shelter: a new sequence for the Late Middle Paleolithic in the Cantabrian region (northern Iberia). Quat. Int. 474, 44–55. doi: 10.1016/j.quaint.2017.06.059

Crossref Full Text | Google Scholar

Hardy, B. L., Moncel, M.-H., Daujeard, C., Fernandes, P., Béarez, P., Desclaux, E., et al. (2013). Impossible Neanderthals? Making string, throwing projectiles and catching small game during Marine Isotope Stage 4 (Abri du Maras, France). Quat. Sci. Rev. 82, 23–40. doi: 10.1016/j.quascirev.2013.09.028

Crossref Full Text | Google Scholar

Haynes, G. (1980). Evidence of carnivore gnawing on Pleistocene and recent mammalian bones. Paleobiology 6, 341–351. doi: 10.1017/S0094837300006849

Crossref Full Text | Google Scholar

Herranz-Rodrigo, D., Estaca-Gómez, V., Torres, T., de Baena, J., and Yravedra, J. (2024). El aprovechamiento de Ursus spelaeus en el Paleolítico Medio. Novedades desde la cueva de El Esquilléu (Cillórigo de Liébana, Cantabria). ZEPHYRVS. XCII, 161–173.

Google Scholar

Hockett, B., and Haws, J. (2002). Taphonomic and methodological perspectives of leporid hunting during the upper Paleolithic of the Western Mediterranean basin. J. Archaeol. Method Theory. 9, 269–302. doi: 10.1023/A:1019503030246

Crossref Full Text | Google Scholar

Hockett, B. S. (1996). Corroded, thinned and polished bones created by golden eagles (Aquila chrysaetos): taphonomic implications for archaeological interpretations. J. Archaeol. Sci. 23, 587–591. doi: 10.1006/jasc.1996.0055

Crossref Full Text | Google Scholar

Hockett, B. S., and Bicho, N. F. (2000). The rabbits of picareiro cave: small mammal hunting during the late upper palaeolithic in the Portuguese estremadura. J. Archaeol. Sci. 27, 715–723. doi: 10.1006/jasc.1999.0496

Crossref Full Text | Google Scholar

Hussain, S. T., Weiss, M., and Nielsen, T. K. (2022). Being-with other predators: cultural negotiations of Neanderthal-carnivore relationships in Late Pleistocene Europe. J. Anthropol. Archaeol. 66:101409. doi: 10.1016/j.jaa.2022.101409

Crossref Full Text | Google Scholar

Jones, E. L. (2006). Prey choice, mass collecting, and the wild European rabbit (Oryctolagus cuniculus). J. Anthropol. Archaeol. 25, 275–289. doi: 10.1016/j.jaa.2005.11.002

Crossref Full Text | Google Scholar

Kassambara, A., and Mundt, F. (2020). _factoextra: Extract and Visualize the Results of Multivariate Data Analyses_. R package version 1.0.7. Available at: https://CRAN.R-project.org/package=factoextra (accessed July, 2024).

Google Scholar

Krajcarz, M., and Krajcarz, M. T. (2012). The Red Fox (Vulpes vulpes) as an accumulator of bones in cave-like environments. Int. J. Osteoarchaeol. 24, 459–475. doi: 10.1002/oa.2233

Crossref Full Text | Google Scholar

Kuhn, M. (2008). Building predictive models in R using the caret package. J. Stat. Softw. 28, 1–26. doi: 10.18637/jss.v028.i05

Crossref Full Text | Google Scholar

Kuhn, M., and Johnson, K. (2013). Applied Predictive Modeling, Vol. 26. New York, NY: Springer, 13. doi: 10.1007/978-1-4614-6849-3

Crossref Full Text | Google Scholar

Landt, M. J. (2007). Tooth marks and human consumption: ethnoarchaeological mastication research among foragers of the Central African Republic. J. Archaeol. Sci. 34, 1629–1640. doi: 10.1016/j.jas.2006.12.001

Crossref Full Text | Google Scholar

Lattao, V., Hipólito, C. G., Hugo, G. S. G., Maria, H., Elena, M., Pierluigi, R., et al. (2023). Cave sediments stable isotopic analysis to attempt an environmental reconstruction: the case of Cadaval and Escoural (Portugal) and Postes (Spain) caves XI. Coimbra: Congresso Nacional de Geologia.

Google Scholar

Le, S., Josse, J., and Husson, F. (2008). FactoMineR: an R package for multivariate analysis. J. Stat. Softw. 25, 1–18. doi: 10.18637/jss.v025.i01

Crossref Full Text | Google Scholar

Liaw, A., and Wiener, M. (2002). Classification and regression by random forest. R News 2, 18–22.

Google Scholar

Linares-Matás, G. J., and Yravedra, J. (2024). Competing forces: Subsistence strategies and human-carnivore interactions during the middle to Upper Palaeolithic transition in Northern Iberia. Quat. Sci. Rev. 334:108703. doi: 10.1016/j.quascirev.2024.108703

Crossref Full Text | Google Scholar

Livraghi, A., Fanfarillo, G., Colle, M. D., Romandini, M., and Peresani, M. (2021). Neanderthal ecology and the exploitation of cervids and bovids at the onset of MIS4: a study on De Nadale cave, Italy. Quat. Int. 586, 24–41. doi: 10.1016/j.quaint.2019.11.024

Crossref Full Text | Google Scholar

Lloveras, L., Moreno García, M., and Soler i Masferrer, N. (2010). The application of actualistic studies to assess the taphonomic origin of Musterian rabbit accumulations from Arbreda Cave (North-East Iberia). Archaeofauna Int. J. Archaeozool. 19, 99–119.

Google Scholar

Lloveras, L., Moreno-Garcia, M., and Nadal, J. (2008a). Taphonomic study of leporid remains accumulated by the Spanish Imperial Eagle (Aquila adalberti). Geobios 41, 91–100. doi: 10.1016/j.geobios.2006.11.009

Crossref Full Text | Google Scholar

Lloveras, L., Moreno-García, M., and Nadal, J. (2008b). Taphonomic analysis of leporid remains obtained from modern Iberian lynx (Lynx pardinus) scats. J. Archaeol. Sci. 35, 1–13. doi: 10.1016/j.jas.2007.02.005

Crossref Full Text | Google Scholar

Lloveras, L., Moreno-García, M., and Nadal, J. (2009). The eagle owl (Bubo bubo) as a leporid remains accumulator: taphonomic analysis of modern rabbit remains recovered from nests of this predator. Int. J. Osteoarchaeol. 19, 573–592. doi: 10.1002/oa.995

Crossref Full Text | Google Scholar

Lloveras, L., Moreno-García, M., and Nadal, J. (2012a). Assessing the variability in taphonomic studies of modern leporid remains from Eagle Owl (Bubo bubo) nest assemblages: the importance of age of prey. J. Archaeol. Sci. 39, 3754–3764. doi: 10.1016/j.jas.2012.06.033

Crossref Full Text | Google Scholar

Lloveras, L., Moreno-García, M., and Nadal, J. (2012b). Feeding the foxes: an experimental study to assess their taphonomic signature on leporid remains. Int. J. Osteoarchaeol. 22, 577–590. doi: 10.1002/oa.1280

Crossref Full Text | Google Scholar

Lloveras, L., Moreno-García, M., Nadal, J., and Thomas, R. (2014). Blind test evaluation of accuracy in the identification and quantification of digestion corrosion damage on leporid bones. Quat. Int. 330, 150–155. doi: 10.1016/j.quaint.2013.07.033

Crossref Full Text | Google Scholar

Lloveras, L., Nadal, J., and Fullola, J. M. (2020). Distinguishing the taphonomic signature of wolves from humans and other predators on small prey assemblages. Sci. Rep. 10:8030. doi: 10.1038/s41598-020-64716-8

PubMed Abstract | Crossref Full Text | Google Scholar

Lloveras, L., Thomas, R., Cosso, A., Pinyol, C., and Nadal, J. (2016). When wildcats feed on rabbits: an experimental study to understand the taphonomic signature of European wildcats (Felis silvestris silvestris). Archaeol. Anthropol. Sci. 10, 449–464. doi: 10.1007/s12520-016-0364-6

Crossref Full Text | Google Scholar

Lyman, R. L. (1994). Relative abundance of skeletal specimens and taphonomic analysis of vertebrate remains. Palaios 9, 288–298. doi: 10.2307/3515203

PubMed Abstract | Crossref Full Text | Google Scholar

Lyman, R. L. (2008). Quantitative Paleozoology. Cambridge, MA: Cambridge University Press. doi: 10.1017/CBO9780511813863

Crossref Full Text | Google Scholar

Maechler, M., Rousseeuw, P., Struyf, A., Hubert, M., and Hornik, K. (2023). cluster: Cluster Analysis Basics and Extensions. R package version 2.1.6.

Google Scholar

Marean, C. W., Abe, Y., Nilssen, P. J., and Stone, E. C. (2001). Estimating the minimum number of skeletal elements (MNE) in zooarchaeology: a review and a new image-analysis GIS approach. Am. Antiq. 66, 333–348. doi: 10.2307/2694612

PubMed Abstract | Crossref Full Text | Google Scholar

Marín, J., Daujeard, C., Saladié, P., Rodríguez-Hidalgo, A., Vettese, D., Rivals, F., et al. (2020). Neanderthal faunal exploitation and settlement dynamics at the Abri du Maras, level 5 (south-eastern France). Quat. Sci. Rev. 243:106472. doi: 10.1016/j.quascirev.2020.106472

Crossref Full Text | Google Scholar

Marín, J., Rodríguez-Hidalgo, A., Vallverdú, J., de Soler, B. G., Rivals, F., Rabuñal, J. R., et al. (2019). Neanderthal logistic mobility during MIS3: zooarchaeological perspective of Abric Romaní level P (Spain). Quat. Sci. Rev. 225:106033. doi: 10.1016/j.quascirev.2019.106033

Crossref Full Text | Google Scholar

Marín, J., Saladié, P., Rodríguez-Hidalgo, A., and Carbonell, E. (2017). Neanderthal hunting strategies inferred from mortality profiles within the Abric Romaní sequence. PLoS ONE 12:e0186970. doi: 10.1371/journal.pone.0186970

PubMed Abstract | Crossref Full Text | Google Scholar

Marwick, B. (2017). Computational reproducibility in archaeological research: basic principles and a case study of their implementation. J. Archaeol. Method Theory, 24, 424–450. doi: 10.1007/s10816-015-9272-9

Crossref Full Text | Google Scholar

Marwick, B., Boettiger, C., and Mullen, L. (2018). Packaging data analytical work reproducibly using R (and friends). Am. Stat. 72, 80–88. doi: 10.1080/00031305.2017.1375986

Crossref Full Text | Google Scholar

Massigoge, A., Gutierrez, M. A., Alvarez, M. C., Kaufmann, C. A., Rafuse, D. J., Gonzalez, M. E., et al. (2014). Estudio comparativo de las marcas de dientes producidas por dos pequeños carnívoros sudamericanos. Rev. Chil. Antropol. doi: 10.5354/0719-1472.2014.36264

Crossref Full Text | Google Scholar

McPherron, S. P. (2022). E5 [Computer software]. Available at: https://github.com/surf3s/E5 (accessed November, 2023).

Google Scholar

Meyer, D., Dimitriadou, E., Hornik, K., Weingessel, A., and Leisch, F. (2023). _e1071: Misc Functions of the Department of Statistics, Probability Theory Group (Formerly: E1071), TU Wien_. R package version 1.7-14. Available at: https://CRAN.R-project.org/package=e1071 (accessed July, 2024).

Google Scholar

Moncel, M. H., Fernandes, P., Willmes, M., James, H., and Grün, R. (2019). Rocks, teeth, and tools: new insights into early Neanderthal mobility strategies in South-Eastern France from lithic reconstructions and strontium isotope analysis. PLoS ONE 14:e0214925. doi: 10.1371/journal.pone.0214925

PubMed Abstract | Crossref Full Text | Google Scholar

Nabais, M., and Zilhão, J. (2019). The consumption of tortoise among last interglacial Iberian Neanderthals. Quat. Sci. Rev. 217, 225–246. doi: 10.1016/j.quascirev.2019.03.024

Crossref Full Text | Google Scholar

Otte, M., and da Silva, A. C. (Eds.). (1996). Recherches préhistoriques à la grotte d'Escoural, Portugal. Liège: Service de Préhistoire, Université de Liège.

Google Scholar

Pavao, B., and Stahl, P. (1999). Structural density assays of leporid skeletal elements with implications for taphonomic, actualistic and archaeological research. J. Archaeol. Sci. 26, 53–66. doi: 10.1006/jasc.1998.0299

Crossref Full Text | Google Scholar

Pelletier, M., Brugal, J.-P., Cochard, D., Lenoble, A., Mallye, J.-B., Royer, A., et al. (2016). Identifying fossil rabbit warrens: insights from a taphonomical analysis of a modern warren. J. Archaeol. Sci. Rep. 10, 331–344. doi: 10.1016/j.jasrep.2016.10.016

Crossref Full Text | Google Scholar

Pelletier, M., Desclaux, E., Brugal, J.-P., and Texier, P.-J. (2019). The exploitation of rabbits for food and pelts by last interglacial Neandertals. Quat. Sci. Rev. 2019:105972. doi: 10.1016/j.quascirev.2019.105972

Crossref Full Text | Google Scholar

Pérez Ripoll, M. (1977). Los mamíferos del yacimiento musteriense de Cova Negra (Játiva, Valencia), Serie Trabajos Varios S.I.P. 23. Valencia: Diputación de Valencia.

Google Scholar

Pérez Ripoll, M. (1992). Marcas de Carnívoros, fracturas incipientes y mordeduras de carnívoros en huesos prehistóricos del Mediterráneo español. Alicante: Instituto de Cultura, Juan Gil-Albert (Colección Patrimonio).

Google Scholar

Pérez Ripoll, M. (2001). “Les marques antròpiques en ossos de conill,” in De Neandertals a Cromanyons: l'inici del Poblament Humà a les Terres Valencianes, ed. V. Villaverde (València: Universitat de València), 119–124.

Google Scholar

Pérez, L., Sanchis, A., Hernández, C. M., Galván, B., Sala, R., Mallol, C., et al. (2017). Hearths and bones: an experimental study to explore temporality in archaeological contexts based on taphonomical changes in burnt bones. J. Archaeol. Sci. Rep. 11, 287–309. doi: 10.1016/j.jasrep.2016.11.036

Crossref Full Text | Google Scholar

Picin, A., Blasco, R., Arilla, M., Rivals, F., Chacón, M. G., Gómez de Soler, B., et al. (2020). “Short-term Neanderthal occupations and carnivores in the Northeast of Iberian Peninsula,” in Short-term Occupations in Paleolithic Archaeology: Definition and Interpretation, eds. J. Cascalheira, and A. Picin (Cham: Springer), 183–213. doi: 10.1007/978-3-030-27403-0_8

Crossref Full Text | Google Scholar

Pickering, T. R., and Egeland, C. P. (2006). Experimental patterns of hammerstone percussion damage on bones: implications for inferences of carcass processing by humans. J. Archaeol. Sci. 33, 459–469. doi: 10.1016/j.jas.2005.09.001

Crossref Full Text | Google Scholar

Pinto-Llona, A. C., Estaca, V., Grandal-d'Anglade, A., Romero, A. J., and Yravedra, J. (2023). Alternation between humans and carnivores in the occupations of the Mousterian site of Sopeña rock-shelter (Asturias, Spain). Quat. Sci. Rev. 108468. doi: 10.1016/j.quascirev.2023.108468

Crossref Full Text | Google Scholar

Real, C. (2020). Rabbit: more than the Magdalenian main dish in the Iberian Mediterranean region. new data from Cova de les Cendres (Alicante, Spain). J. Archaeol. Sci. Rep. 32:102388. doi: 10.1016/j.jasrep.2020.102388

Crossref Full Text | Google Scholar

Reitz, E. J., and Wing, E. S. (1999). Zooarchaeology. Cambridge University Press.

Google Scholar

Richards, M., Pettitt, P., Trinkaus, E., Smith, F., Paunović, M., Karavanić, I., et al. (2000). Neanderthal diet at Vindija and Neanderthal predation: the evidence from stable isotopes. Proc. Natl. Acad. Sci. USA. 97, 7663–7666. doi: 10.1073/pnas.120178997

PubMed Abstract | Crossref Full Text | Google Scholar

Ripoll, M. P. (2004). “La consommation humaine des lapins pendant le Paléolithique dans la région de Valencia (Espagne) et l'étude des niveaux gravétiens de La Cova de les Cendres (Alicante),” in Actes des XXIV rencontres internationales d'Archéologie et d'Histoire d'Antibes (Antibes: Editions APDCA), 191–206.

Google Scholar

Rivals, F., Bocherens, H., Camarós, E., and Rosell, J. (2022). “Diet and ecological interactions in the Middle and Late Pleistocene,” in Updating Neanderthals (Academic Press), 39–54.

Google Scholar

Rodríguez-Hidalgo, A. (2022). Comments on ‘Hyaenids, felids and canids as bone accumulators: Does the natural history of extant species support zooarchaeological inferences?', by Palomares et al. [Quat. Sci. Rev. 284 (2022) 107459]. Quatern. Sci. Rev. 295:107632. doi: 10.1016/j.quascirev.2022.107632

Crossref Full Text | Google Scholar

Rodríguez-Hidalgo, A., Lloveras, L., Moreno-García, M., Saladié, P., Canals, A., Nadal, J., et al. (2013). Feeding behaviour and taphonomic characterization of non-ingested rabbit remains produced by the Iberian lynx (Lynx pardinus). J. Archaeol. Sci. 40, 3031–3045. doi: 10.1016/j.jas.2013.03.006

Crossref Full Text | Google Scholar

Rodríguez-Hidalgo, A., Saladié, P., Marín, J., and Canals, A. (2015). Expansion of the referential framework for the rabbit fossil accumulations generated by Iberian lynx. Palaeogeogr. Palaeoclimatol. Palaeoecol. 418, 1–11. doi: 10.1016/j.palaeo.2014.11.010

Crossref Full Text | Google Scholar

Rodríguez-Hidalgo, A., Sanz, M., Daura, J., and Sánchez-Marco, A. (2020). Taphonomic criteria for identifying Iberian lynx dens in quaternary deposits. Sci. Rep. 10:7225. doi: 10.1038/s41598-020-63908-6

PubMed Abstract | Crossref Full Text | Google Scholar

Roebroeks, W., and Soressi, M. (2016). Neandertals revised. Proc. Nat. Acad. Sci. 113, 6372–6379. doi: 10.1073/pnas.1521269113

PubMed Abstract | Crossref Full Text | Google Scholar

Romandini, M., Silvestrini, S., Real, C., Lugli, F., Tassoni, L., Carrera, L., et al. (2023). Late Neanderthal “menu” from northern to southern Italy: freshwater and terrestrial animal resources. Quat. Sci. Rev. 2023:108233. doi: 10.1016/j.quascirev.2023.108233

Crossref Full Text | Google Scholar

Rosell, J., Blasco, R., Rivals, F., Chacón, M. G., Arilla, M., Camarós, E., et al. (2017). A resilient landscape at Teixoneres Cave (MIS 3; Moià, Barcelona, Spain): the Neanderthals as disrupting agent. Quatern. Int. 435, 195–210. doi: 10.1016/j.quaint.2015.11.077

Crossref Full Text | Google Scholar

Rufà, A., Blasco, R., Rivals, F., and Rosell, J. (2014). Leporids as a potential resource for predators (hominins, mammalian carnivores, raptors): an example of mixed contribution from level III of Teixoneres Cave (MIS 3, Barcelona, Spain). Comptes Rendus Palevol. 13, 665–680. doi: 10.1016/j.crpv.2014.06.001

Crossref Full Text | Google Scholar

Russo, G., Milks, A., Leder, D., Koddenberg, T., Starkovich, B. M., Duval, M., y Terberger, T., et al. (2023). First direct evidence of lion hunting and the early use of a lion pelt by Neanderthals. Sci. Rep. 13:16405. doi: 10.1038/s41598-023-42764-0

PubMed Abstract | Crossref Full Text | Google Scholar

Saladié, P., Rodríguez-Hidalgo, A., Marín, J., i Poch, J. V., and Carbonell, E. (2018). The top of the Gran Dolina (Atapuerca, Spain) sequence: a zooarchaeological and occupational perspective. Quatern. Sci. Rev. 195, 48–71. doi: 10.1016/j.quascirev.2018.07.010

Crossref Full Text | Google Scholar

Saladié, P., Rosas, A., García-Tabernero, A., Fidalgo, D., Fero Meñe, M., Ebana Ebana, C., et al. (2024). An actualistic taphonomic model of human tooth marks on bone remains: a sample recovered in villages of continental Equatorial Guinea. J. Archaeol. Sci. Rep. 55:104514. doi: 10.1016/j.jasrep.2024.104514

Crossref Full Text | Google Scholar

Sanchis Serra, A., and Pascual Benito, J. L. (2011). Análisis de las acumulaciones óseas de una guarida de pequeños mamíferos carnívoros (Sitjar Baix, Onda, Castellón): implicaciones arqueológicas. Archaeofauna 20, 47–71. doi: 10.15366/archaeofauna2011.20.003

Crossref Full Text | Google Scholar

Sanchis, A. (1999). Análisis tafonómico de los restos de Oryctolagus cuniculus a partir de la alimentación de Bubo bubo y Vulpes vulpes y su comparación con materiales antrópicos (Degree Thesis). València: Universitat de València, 260.

Google Scholar

Sanchis, A., Morales Pérez, J., and PérEZ-riPoll, M. (2011). “Creación de un referente experimental para el estudio de las alteraciones causadas por dientes humanos sobre huesos de conejo,” in Actas Del Segundo Congreso Internacional de Arqueología Experimental, Ronda, Málaga, eds. A. Morgado, J. Baena, and D. García (Granada: Universidad de Granada), 343–349.

Google Scholar

Scheifler, N. A., Álvarez, M. C., Rafuse, D. J., Kaufmann, C. A., Massigoge, A., González, M. E., et al. (2020). Taphonomic signature of Geoffroy's cat (Leopardus geoffroyi) on small sized preys: a comparative study of ingested and non-ingested leporid bones. J. Archaeol. Sci. Rep. 31:102340. doi: 10.1016/j.jasrep.2020.102340

Crossref Full Text | Google Scholar

Schmitt, D. N., and Juell, K. E. (1994). Toward the identification of coyote scatological faunal accumulations in archaeological contexts. J. Archaeol. Sci. 21, 249–262. doi: 10.1006/jasc.1994.1025

Crossref Full Text | Google Scholar

Seuru, S., Burke, A., and Perez, L. (2024). Evidence of an age and/or gender-based division of labor during the Last Glacial Maximum in Iberia through rabbit hunting. J. Archaeol. Sci. Rep. 56:104560. doi: 10.1016/j.jasrep.2024.104560

Crossref Full Text | Google Scholar

Skov, L., Peyrégne, S., Popli, D., Iasi, L. N. M., Devièse, T., Slon, V., et al. (2022). Genetic insights into the social organization of Neanderthals. Nature 610, 519–525. doi: 10.1038/s41586-022-05283-y

PubMed Abstract | Crossref Full Text | Google Scholar

Smith, G. M. (2014). Neanderthal megafaunal exploitation in Western Europe and its dietary implications: a contextual reassessment of La Cotte de St Brelade (Jersey). J. Hum. Evol. 10:007. doi: 10.1016/j.jhevol.2014.10.007

PubMed Abstract | Crossref Full Text | Google Scholar

Stekhoven, D. J., and Buehlmann, P. (2012). MissForest - non-parametric missing value imputation for mixed-type data. Bioinformatics 28, 112–118. doi: 10.1093/bioinformatics/btr597

PubMed Abstract | Crossref Full Text | Google Scholar

Stiner, M. C., Munro, N. D., and Surovell, T. A. (2000). The tortoise and hare. Small-game use, the broad-spectrum revolution, and paleolithic demography. Curr. Anthropol. 41, 39–73. doi: 10.1086/300102

PubMed Abstract | Crossref Full Text | Google Scholar

Stiner, M. C., Munro, N. D., Surovell, T. A., Tchernov, E., and Bar-Yosef, O. (1999). Paleolithic population growth pulses evidenced by small animal explotation. Science 283, 190–193. doi: 10.1126/science.283.5399.190

PubMed Abstract | Crossref Full Text | Google Scholar

Téllez, E., Saladié, P., Pineda, A., Marín, J., Vallverdú, J., Chacón, M. G., et al. (2022). Incidental burning on bones by Neanderthals: the role of fire in the Qa level of Abric Romaní rock-shelter (Spain). Archaeol. Anthropol. Sci. 14:119. doi: 10.1007/s12520-022-01577-4

Crossref Full Text | Google Scholar

Vabalas, A., Gowen, E., Poliakoff, E., and Casson, A. J. (2019). Machine learning algorithm validation with a limited sample size. PLoS ONE 14:e0224365. doi: 10.1371/journal.pone.0224365

PubMed Abstract | Crossref Full Text | Google Scholar

Valde-Nowak, P., and Cieśla, M. (2020). “Models of raw material exploitation as an indicator of middle paleolithic mobility: case studies from uplands of Northern Central Europe,” in Short-term occupations in Paleolithic archaeology: Definition and Interpretation, eds. J. Cascalheira, and A. Picin (Cham: Springer), 105–120. doi: 10.1007/978-3-030-27403-0_5

Crossref Full Text | Google Scholar

Vaquero, M. (2022). “The organisation of living spaces in Neanderthal campsites,” in Updating Neanderthals (New York, NY: Academic Press), 207–225. doi: 10.1016/B978-0-12-821428-2.00001-9

Crossref Full Text | Google Scholar

Villa, P., Castel, J. C., Beauval, C., Bourdillat, V., and Goldberg, P. (2004). Human and carnivore sites in the European Middle and Upper Paleolithic: similarities and differences in bone modification and fragmentation. Rev. Paléobiol. 23, 705–730.

Google Scholar

Villa, P., and Mahieu, E. (1991). Breakage patterns of human long bones. J. Hum. Evol. 21, 27–48. doi: 10.1016/0047-2484(91)90034-S

Crossref Full Text | Google Scholar

White, T. E. (1953). A method of calculating the dietary percentage of various food animals utilized by aboriginal peoples. Am. Antiq. 18, 396–398. doi: 10.2307/277116

Crossref Full Text | Google Scholar

Wickham, H. (2015). R Packages: Organize, Test, Document, and Share Your Code. Sebastopol, CA: O'ReillyMedia. Inc.

Google Scholar

Yravedra, J. (2004). “Implications taphonomiques des modifications osseuses faites par les vrais Hiboux (Bubo bubo) sur les lapins,” in Pettits animaux et sociétés. Du complement alimentaire aux ressources utilitaires XXIV. Reencontres Internationales d'Archeologie et d'histoire d'Antibes, eds. J. P. Brugal and J. Desse (Antibes: APDCA), 321–323. I.S.B.N. 2-904110-38-0

Google Scholar

Yravedra, J. (2006). Acumulaciones biológicas en yacimientos arqueológicos: amalda VII y esquilleu III-IV. Trabajos Prehist. 63, 55–78.

Google Scholar

Yravedra, J. (2008). Los lagomorfos como recurso alimenticio en Cueva de Ambrosio (Almería, España). Zephirus 61, 63–80.

Google Scholar

Yravedra, J., and Andrés, M. (2013). “Estrategias de subsistencia entre los grupos magdalenienses de La Peña de Estebanvela (Ayllón, Segovia),” in Ocupaciones Magdaleniense en el Interior de la Península Ibérica, ed. C. Cacho (Ayllón: La Peña de Estebanvela), 213–244. ISBN: 978-84-616-5992-0.

Google Scholar

Yravedra, J., and Cobo-Sánchez, L. (2015). Neanderthal exploitation of ibex and chamois in southwestern Europe. J. Hum. Evol. 78, 12–32. doi: 10.1016/j.jhevol.2014.10.002

PubMed Abstract | Crossref Full Text | Google Scholar

Yravedra, J., and Domínguez-Rodrigo, M. (2009). The shaft-based methodological approach to the quantification of long limb bones and its relevance to understanding hominid subsistence in the Pleistocene: application to four Palaeolithic sites. J. Quat. Sci. 24, 85–96. doi: 10.1002/jqs.1164

Crossref Full Text | Google Scholar

Yravedra, J., Estaca-Gómez, V., Grandal-d'Anglade, A., and Pinto-Llona, A. C. (2024). Neanderthal subsistence strategies: new evidence from the Mousterian Level XV of the Sopeña rock shelter (Asturias, northern Spain). Archaeol. Anthropol. Sci. 16:6. doi: 10.1007/s12520-023-01914-1

Crossref Full Text | Google Scholar

Yravedra, J., Herranz, D., Sesé, C., López-Cisneros, P., Linares-Matás, G. J., Pernas-Hernández, M., et al. (2019). Lagomorph exploitation during the upper Palaeolithic in the Northern Iberian Peninsula. new evidence from Coímbre Cave (Asturias, Spain). Quat. Int. 506, 59–68. doi: 10.1016/j.quaint.2018.06.016

Crossref Full Text | Google Scholar