Between the mining trench and the excavation of the shaft (cross-section J/K), 17 stratigraphic units (I–XVII) have been documented with a thickness of more than 4 m, varying in thickness, color, and texture (Figure 2B). The average thickness is 18 cm (σ = 9.3) and ranges between 2 and 50 cm, with the lower levels having the greatest thickness. They exhibit frequent erosive contacts and a sub-horizontal or northeastward dipping arrangement (levels VII–XV), with a bulging morphology that suggests post-depositional collapse processes. The colors (Figure 2C) range from various shades of brown, gray, pink, and reddish hues, with a dominant texture of silts and sands and a notable granular structure. The consistency is highly heterogeneous, varying from soft to hard depending on the extent of secondary calcite precipitation affecting each level. Porosity is generally absent or minimal (when present, it is always vacuolar in type), and the contacts between the layers are mostly sharp (Barroso Ruíz et al., 2011). So far, no hiatus has been detected in the stratigraphy.

Based on sedimentological, geomorphological, and geochemical analyses (Barroso Ruíz et al., 2011), it has been determined that this area (profile J/K) of the site was heavily exposed to anthropogenic fire activity. The stratigraphic units have been grouped into three major structural complexes (Monge et al., 2014): I (I–VII), II (VIII–XII), and III (XIII–XVII), with evidence suggesting a possible increase in thermo-alteration intensity from complex III to I, as indicated by the increase in organic carbon from complex III to I (Barroso Ruíz et al., 2011). In contrast, moving southward from this area, in trench 8, evidence of fire activity appears to be significantly lower (Monge et al., 2020). These data support the hypothesis of the existence of a recurrently used fireplace at the same location for thousands of years (Barroso Ruíz et al., 2011; Solano García, 2024). The presence of lithic industry remains and burned bone in all levels (Barroso Ruíz et al., 2011), along with low carnivore activity (Solano García, 2024) and the absence of guano (Monge et al., 2016), indicates a recurrent use of the site by human groups.

The chronology of the site was initially estimated to be over 130,000 years based on ²³⁰Th/²³⁴U dating of a sample from the upper part of the archaeological fill (Botella Ortega et al., 2006). The analysis of lithic industry and faunal spectrum suggested a broad chronological range from MIS 11 to MIS 5 (Barroso Ruíz et al., 2011; Moigne et al., 2016). Shortly afterward, a new series of ESR/U-series dates were published (Falguères et al., 2019) which, despite some limitations inherent to the method, identified ages ranging from MIS 9 to MIS 7, suggesting continuous human occupation between 320–180 ka. Recently, the discovery of a human tooth was reported (Bermúdez et al., 2023) which was dated using amino acid racemization and yielded a date of 104,300 years. However, the tooth came from a loose block from surface cleaning, possibly originating from previously eroded upper layers. Nevertheless, this date allows us to estimate an upper limit for the original fill, providing a timeframe from MIS-9/10 to MIS-5d, though the chronology of the lower part of the stratigraphic sequence remains uncertain (Figure 2B).

In this study, remains from grids J7, K6, K7, and K8, which make up the J/K stratigraphic profile (Figure 1D), were used. We analyzed 3,700 remains, of which 2,000 are indeterminate fragments smaller than 2 cm that cannot even be classified as long bone or non-long bone (Pickering et al., 2003). These small fragments fall outside the scope of the zooarchaeological analysis of this study but were used to assess the overall fragmentation of the assemblage (Solano García, 2024). Stratigraphic units IV–XVII were selected from the entire sequence, while levels I to III were not considered due to the low number of remains, as these levels could not have been systematically excavated in this area. The analysis focuses on the exploitation of ungulates throughout the stratigraphic sequence. The presence of carnivores is attested by a low proportion of remains from Ursus arctos, Lynx pardinus, Canis lupus, and Felis silvestris across the entire sequence (Barroso Ruíz et al., 2011; Solano García, 2024) (Supplementary Table S1), which will be referred to as a general carnivore group (Table 1). Small fauna, such as birds, lagomorphs, tortoises, amphibians, rodents, etc. (Barroso Ruíz et al., 2011), were also excluded from the analysis, as they fall outside the objectives of this study (see Solano García, 2024 for more details). The sample used consists of 3,514 remains, including 1,394 identified herbivore remains (NISP), 145 carnivore remains, 1,569 remains classified by size category, and 406 indeterminate remains.

The preservation state of the fossil remains is generally good, as they exhibit a high degree of mineralization. Between 65 and 80% of the remains, depending on the level, show no cortical modifications, and when they do, they appear as small dissolution pits. There are no evident signs of weathering (Behrensmeyer, 1978) throughout the sequence, although longitudinal and transverse fissures (25%) filled with precipitated calcite are common, sometimes contributing to the compaction of the fossil due to the weight of the sediments. Most of the fissures appear to be a direct result of fire action (Cain, 2005; Stiner, 2005; Théry-Parisot et al., 2010) and their subsequent development through fossil-diagenetic karstification processes. The most frequent alterations are calcite concretions and the presence of manganese oxide distributed as small dendritic patterns on the cortical surface. In contrast, biochemical modifications are scarcely described. Moisture and carbonated waters have been the main agents of alteration, aside from the anthropogenic effects of fire (Solano García, 2024).

The accumulation is anthropogenic in origin, as evidenced by the large volume of lithic tools, retouchers (Moigne et al., 2016), and anthropogenic activities related to butchery and food cooking (Barroso Ruíz et al., 2011). Carnivore remains are scarce throughout the assemblage (3.9%), and their impact is minimal (0.8%), except in unit XIII (1%), where there is a slight increase but with low activity (0.5%). Additionally, the skeletal representation and the dominant adult mortality profile (Stiner, 1994; Steele, 2004), along with the abundance of cut marks, intentional fracturing, and extensive use of fire, also suggest an anthropogenic accumulation. Access to carcasses was likely primary for large ungulates (Bos, Bison, Equus, and Cervus) and secondary for larger animals like elephants or rhinoceroses and probably Sus (Solano García, 2024).

It is important to note that the excavation area of each grid, except for J7, does not reach a full square meter in surface area. The excavation of this zone corresponds to the edge of the mining shaft dug along the K band (Figure 1D). Consequently, depending on the width of the sides and the depth, the excavated surface area will vary. On the other hand, grid J7, which covers the entire surface, began to be excavated in 1995 from level V at level XIII. Therefore, the material density is biased and uneven across different levels and may influence the frequencies of the number of identified remains. Nevertheless, the sample presented allows for the first diachronic zooarchaeological approach to nearly the entire stratigraphic sequence that we have so far.

Taxonomic and anatomical identifications were carried out using the osteological collection of the Center for Prehistoric Research of Tautavel (CERPT) and specialized reference literature (Heintz, 1970; Pales et al., 1971; Schmid, 1972; Hillson, 2005; Barone, 2010). At the taxonomic level, some remains could be identified as Bos primigenius and Bos bison (Brugal, 1984; Gee, 1993), but for the purposes of this study, they have been unified under Bos/Bison (BB). When taxonomic identification was not possible, size-based criteria (adapted from Bunn, 1981; Bunn et al., 1986)Bunn, 1981; Bunn et al., 1986): were used to categorize animal mass: (1) small, 5–100 kg (e.g., Capra, Dama, Canis, Lynx); (2) medium, 100–250 kg (e.g., Sus scrofa, Cervus elaphus, Equus hydruntinus); (3) medium-large, 250–500 kg (e.g., Equus ferus, Ursidae); (4) large, 500–800 kg (e.g., Bos primigenius, Bison priscus); (5) very large, >800 kg (e.g., Stephanorhinus, Palaeoloxodon). To simplify the analysis, unidentified herbivore remains were grouped into two sets: large herbivores (LH: categories 3, 4, 5), mainly composed of categories 3 and 4; and small herbivores (SH: categories 1, 2). Anatomically unidentified remains were classified into one of the following categories: long bone, flat bone, and articular bone (Saladié et al., 2011).

The determination of age at death, when possible, is based on the fusion of epiphyses and dental eruption/replacement patterns and wear (Mariezkurrena, 1983; Wegrzyn and Serwatka, 1984; Stiner, 1990, 1998; Forsten and Moigne, 1998; Guadelli, 1998; Hillson, 2005; Barone, 2010) and has been classified into infant, juvenile, adult, and senile categories.

Abundance measures for all ungulates were recorded following the number of identified specimens (NISP), the minimum number of elements (MNE), the minimum number of individuals (MNI), and the minimum animal unit (MAU). The skeletal survival index per individual (SI) was calculated using the following formula: MAU × 100/MNI (Binford, 1981; Grayson, 1984; Klein and Cruz Uribe, 1984; Lyman, 1994a, 2008). To synthesize the information for each ungulate by stratigraphic units, the MAU for each taxon was calculated based on the total MNE (excluding teeth) and the complete skeleton, along with its corresponding SI. The %ΣMAU represents the proportion of MAU for each taxon within each stratigraphic unit as a summary value, and not Binford (1981) index, which has been replaced by the SI equivalent. However, it has been verified that the results are equivalent. The carnivore/herbivore percentage ratio (Klein and Cruz Uribe, 1984) has also been included. Additionally, the total MAU for each stratigraphic unit was calculated by summing all the complete skeletons present for each taxon.

For the analysis of the anatomical representation of the stratigraphic sequence, given the large number of skeletal elements and with the aim of identifying general exploitation patterns and synthesizing the results, the elements have been grouped into six main categories: skull, axial, long bones, girdles (scapula + coxal), articular bones, and phalanges. However, individual elements will also be used for analysis, and they will be referenced when necessary.

To analyze the fragmentation of the ungulate assemblage by stratigraphic unit, the MNE/NISP ratio was used, applied also to long bones. The intensity of fragmentation was calculated by counting incomplete elements (fragments/NISP × 100) and comparing it with the overall set of fragmented and complete remains. The degree of fragmentation was also measured based on the integrity of the remains and the circumference size of long bones (Villa and Mahieu, 1991). These data were complemented by a descriptive statistical analysis of the length of the remains and a count of all small indeterminate splinters not assigned to size categories.

Additionally, to verify the intensity of fracturing, all fragmented and complete long bones were counted by stratigraphic unit, taking into account the diaphysis and epiphysis ratios (Todd and Rapson, 1988; Lyman, 1994b, 2008). The count of fractured phalanges was also used as an indicator of intensified marrow exploitation (Starkovich, 2017) along with the count of mandibles. The longitudinal opening of these types of elements is a clear distinctive sign of marrow extraction.

Anthropogenic activity has been synthesized based on the count of cut marks (Binford, 1981; Blumenschine et al., 1996; Domínguez-Rodrigo et al., 2009), percussion marks (Blumenschine and Selvaggio, 1988; Blumenschine et al., 1996; Pickering and Egeland, 2006; De Juana and Domínguez-Rodrigo, 2011), fresh bone fractures (Villa and Mahieu, 1991) and burned bones (Shipman et al., 1984; Stiner et al., 1995). Additionally, the number of carnivore tooth marks (Bunn, 1981; Fisher, 1995; Blumenschine et al., 1996) has been counted as contrasting information, highlighting carnivores as modifying agents in the assemblage.

The observation of surface bone modifications was carried out using a 10× handheld magnifying lens and a Dinolite device. The graphical analysis of the anatomical representation of each taxon by stratigraphic unit was performed using ArcGIS^® software by Esri (ESRI Inc, 2023), where the MAU normalized by NISP was utilized (Orton, 2010). All statistical analysis and graphical output were conducted using the R programming language (R Core Team, 2024).

This study analyzed 1,394 remains (39.67%) identified as belonging to ungulates. As complementary information, we also present the count of unidentified material (56.2%), classified into the categories of Large Herbivores and Small Herbivores (44.65%), as well as the number of carnivore remains (4.13%) and indeterminate remains (11.55%) for each stratigraphic unit (3,514 NR; Table 1).

The number of remains is quite disproportionate throughout the entire sequence, as is the NISP (CV = 124.2 and 86.2, respectively), with a marked bias in the upper layers (IV–VI) and layer XVII, while a higher number is observed in the central layers. This bias is attributed to the smaller volume excavated in these levels, although level V may also show a potential research bias with a lower identification rate of around 0.2%, compared to the 3.1% average for the total NR. Despite the marked discrepancies in NISP throughout the sequence, we have confirmed that the use of MAU and SI as analytical variables does not affect the results and is correlatable with the material density per unit of analysis.

A total of nine ungulate taxa have been identified throughout the sequence: Bos primigenius, Bison priscus, Cervus elaphus, Dama dama, Sus scrofa, Equus ferus, Stephanorhinus hemitoechus, Palaeoloxodon antiquus, and Capra sp. (Barroso Ruíz et al., 2011). The first six taxa, counting B. priscus and B. primigenius as a single group (BB), are present in almost the entire sequence, with the exception of Dama, which is absent in levels V, VII, and XIV; Stephanorhinus in levels V and XVII; and BB in level XVII. On the other hand, Palaeoloxodon is only found in levels VII, XII, and XVII, and Capra in level XVII. The absence of Dama and Stephanorhinus in level V may be due to the possible excavation or research biases mentioned earlier.

Carnivores are present in all levels of the analyzed sequence (Table 1, Supplementary Table S1), although there is considerable variability in the number of remains, ranging from just one in level VIII to 39 remains in level XIII. Unit VI has the highest carnivore-herbivore ratio (17.14%), followed by unit XIII (11.34%), within an overall carnivore-herbivore ratio of 4.89%. Among the carnivores, lynx remains are the most abundant (46.1%), followed by bear (39.2%), wolf (11.8%), and wildcat (2.9%). Bears have the widest distribution across the entire J/K stratigraphic sequence, except in levels V to VIII, followed by lynx, present in eight levels, wolves in six, and wildcats in three. However, according to Barroso Ruíz et al. (2011), brown bears are the most represented in NR across the entire site.

The predominant taxon among herbivores throughout the stratigraphic sequence (Figure 3), based on representation indices (NISP, x¯ = 47.2, σ = 46.8; MNE, x¯ = 28.1, σ = 25.5; MAU, x¯ = 0.16, σ = 0.11), is the horse, followed by BB (NISP, x¯ = 19.8, σ = 19.5; MNE, x¯ = 12.85, σ = 9.95; MAU, x¯ = 0.08, σ = 0.06) and red deer, with very similar values (NISP, x¯ = 20.2, σ = 15.2; MNE, x¯ = 12.8, σ = 7.55; MAU, x¯ = 0.10, σ = 0.06), although the red deer's MAU is slightly higher or comparable. The rest of the taxa, with lower representativity in terms of the number of remains and elements, are represented by wild boar, fallow deer, rhinoceros, elephant, and ibex. The number of individuals also confirms this taxonomic hierarchy, with horse (39), BB (26), red deer (24), wild boar (20), but more rhinoceros individuals (13) than fallow deer (11), and lastly, elephant (3) and ibex (1).

The total MNI (Minimum Number of Individuals) is 137 for the ungulate species (Table 1). The vast majority correspond to adult individuals (74.45%), compared to 16.8% juveniles, 6.6% subadults, and 2.2% elderly. The central levels, between VIII and XV, have the highest proportion of individuals, with a peak in level X, which contains 20 individuals, including two of the three elderly ones. Levels IX and XV have the highest proportion of elements relative to those expected (MNE/MNI; 12.3 and 11.2, respectively), followed by VII (8.2), XI, and XIII (7.1). The taxon with the greatest anatomical representation in this case is Equus (10.1), followed by Cervus (7.5) and BB (6.4), on one hand, and Dama (4), Stephanorhinus (3), and Sus (2.9), on the other. The species with the largest number of juveniles are Equus and Cervus (25% and 25.6%, respectively), followed by wild boar (20%), BB, and rhinoceros (7.7%). Subadults are present in BB (11.5%), Equus (10.3%), Sus (5%), and Cervus (4.2%), while elderly individuals are found in Equus (5.1%) and Cervus (4.2%). The remainder belong to adults, with proportions above 60%. In the case of fallow deer, all individuals are adults.

It is clear that the means of the main ungulate taxa show variability (Figure 3). After comparing the assumptions of normality and homogeneity of variances (S-W and Levene tests) for the indices of the six main ungulates (NISP: D = 0.26, df = 80, p < 0.001; MNE: D = 0.24, df = 80, p < 0.001; MAU: D = 0.207, df = 80, p < 0.001), a K–W test, effect size (η²; Wolverton et al., 2016), and post-hoc tests were conducted. The results indicate statistically significant differences between the means of the six taxa (NISP: H = 29.35, df = 5, p < 0.001, η² = 0.338; MNE: H = 31.29, df = 5, p < 0.001, η² = 0.365; MAU: H = 38.48, df = 5, p < 0.001, η² = 0.465; SI: H = 35.795, df = 5, p < 0.001, η² = 0.428). These results suggest the differences are reliable. Post-hoc comparisons (Wilcox & Conover-Iman test using Holm's method) identified a greater similarity between BB, Cervus, and Equus (Group 1) (NISP: H = 2.97, df = 2, p = 0.226, η² = 0.0129; MNE: H = 3.68, df = 2, p = 0.159, η² = 0.0223; MAU: H = 3.81, df = 2, p = 0.149, η² = 0.0242) and between Dama, Sus, and Rhino (Group 2) (NISP: H = 0.34, df = 2, p = 0.843, η² = −0.0221; MNE: H = 0.62, df = 2, p = 0.734, η² = −0.0184; MAU: H = 2.87, df = 2, p = 0.238, η² = 0.0116). In summary, we can assume two main distinct groups based on their representativity in the stratigraphic sequence, with the first group showing significantly higher representativity compared to the second (Figure 4).

The results indicate that Group 1 forms the core trio of species with the highest taxonomic and anatomical representation within the assemblage and throughout the entire stratigraphic sequence. This does not imply that there are no differences between stratigraphic levels, as it is clear that in some levels the representation of these taxa is higher than others, but they always maintain a similar presence per level, comparable to the rest of the layers. Group 2 behaves similarly according to the data, but with a much lower representation. When we applied the same test for the survival index (%SI), which takes into account the MNI, it confirmed the similarity within the first group (%SI: H = 4.73, df = 2, p = 0.94, η² = 0.0364), and within the second group (%SI: H = 8.06, df = 2, p = 0.0177, η² = 0.0808). However, there is greater similarity between Rhino and Sus. This is because the skeletal survival of Dama is more similar to that of BB.

The overall skeletal representation based on MNE is dominated by long bones (37.51%), cranial remains including isolated teeth and antlers (31.81%), and axial skeleton elements (12.2%), followed by articular bones, girdles, and phalanges. However, the MAU (Figure 5) indicates a slightly different distribution, where long bones still dominate (57.57%) along with cranial remains (24.70%), but these are followed by girdles (9.66%), articular bones (3.14%), phalanges (2.84%), and axial elements (2.10%). This same distribution is reflected in the skeletal survival index, although the representation of phalanges (2.8%) is slightly higher than that of articular bones (2.4%).

When comparing the MAU by taxa (Figure 5), the same general distribution is observed, with long bones predominating, followed by cranial elements, phalanges, or articular and axial bones, except for Sus scrofa. The wild boar shows a different pattern, where cranial elements dominate (75%) compared to long bones (20.3%), followed by articular bones (2.5%), phalanges (1.9%), and axial elements (0.2%).

According to the MAU by anatomical elements (Figure 4, Supplementary Figure S1), metapodials are the most represented (16.7%), followed by tibiae (10.40%), humeri (8.73%), femora (8.36%), radii (7.61%), and ulnae (5.57%). However, if we break down all the elements, tibiae (10.28%) are the best represented, followed by cranial elements (9.92%), mandibles (9.00%), metatarsals (8.88%), humeri (8.63%), femora (8.26%), radii (7.53%), scapulae (5.51%), ulnae (5.51%), metacarpals (4.34%), and coxae (4.04%). The remaining elements are below 2%.

At the stratigraphic level, most layers contain skewed (non-symmetrical) data. Only levels X, XI, and XV show data symmetry. However, when we analyse the MAU (Figure 5) as a contrast index, we confirm that there is no statistically significant difference between them, based on the K–W test (H = 14.904, df = 13, p = 0.3134, η² = 0.00362) and the two-sided Jonckheere-Terpstra test (JT = 65,792, p = 0.8461). In contrast, the %SI (Figure 5) shows significant differences (H = 29.878, df = 13, p = 0.005, η² = 0.0321), although the effect size is very close. This is because the differences are solely related to layer X. Thus, the JT test, which is more robust and takes into account the order of the stratigraphic sequence, indicates similarity (JT = 68,622, p = 0.24). This result supports the feasibility of conducting analyses to identify differences and similarities between species groups and anatomical categories for the entire assemblage. In fact, the correspondence analysis between anatomical categories throughout the sequence shows how most levels cluster around the long bones and near the origin of the coordinates, indicating little variation (Figure 6B).

The species with the highest number of elements are Equus, Cervus, BB, Sus, and Dama. However, as we have seen, Sus presents a different skeletal configuration (Supplementary Figure S1). When comparing the total anatomical representation of the four taxa based on the MAU, there are also significant differences (H = 229.56, df = 5, p < 2.2e-16, η² = 0.488), although post-hoc tests indicate a potential similarity between BB, Cervus, and Equus, and Dama with Cervus. However, when pairwise tests are conducted, no correspondence is found. This is due to a certain variability in combinations of categories by elements. When analyzing these skeletal categories for the four mentioned taxa, there is statistically significant similarity between skulls, long bones, and girdles (H = 2.6531, df = 2, p = 0.2654, η² = 0.00197), with the distribution of skulls and girdles being much more similar between BB, Cervus, and Equus (H = 2.0769, df = 2, p = 0.354, η² = 0.000254). These results indicate a similar configuration among the species in Group 1, with a predominance of long bones, followed by skulls and girdles and mandibles in similar proportions (Figure 6A). As we have seen, Dama presents a distribution approaching that of Group 1, but its lower representation, for now, does not allow it to be fully included in that group and instead places it in Group 2, as noted earlier. Group 2 shows statistical similarity between articular bones and phalanges (H = 0.15342, df = 1, p = 0.6953, -η² = 0.0470), long bones and girdles (H = 0.32613, df = 1, p = 0.5679, η² = −0.0157), and skulls and girdles (H = 1.3651, df = 1, p = 0.2427, η² = 0.00936).

The survival index for the entire assemblage (Table 1) shows that the red deer has the highest skeletal survival (34%), followed by the horse (24.3%) and BB (17.5%). In Group 2, with lower indices, the best represented is the fallow deer (11.2%), followed by the rhinoceros (7%) and the wild boar (6.1%). At the skeletal level (Figure 5), long bones are the best represented (58.1%) for all taxa except for the wild boar, where cranial remains predominate. Cranial remains are the second most represented (25.5%) for the other taxa, followed by girdles (9.2%), phalanges (2.8%), and articular bones (2.4%). Among the long bones, metapodials have the highest survival (15.5%), followed by tibiae and humeri (10.5%). These are followed by femora (8.3%), radii (7.4%), and ulnae (5.4%). Scapulae (5.9%) are better represented than pelves (3.3%), and phalanges (2.8%) are better represented than articular bones (2.4%), although the difference is minimal. As can be seen, the survival of the girdles is consistent with that of humeri and femora, considering their anatomical connection. Additionally, it indicates a higher representativity of anterior stylopodia and posterior zeugopodia, although the MAU is similar for the stylopodia, and the tibia stands out more than the radius among the zeugopodia.

Lastly, we conducted an analysis of the distribution of forelimbs and hindlimbs based on the MNE. At first glance, the frequency of elements is markedly uneven between Groups 1 and 2 (Figure 7A). As a preliminary approach, we calculated the Log (Forelimb/Hindlimb) ratio for all species with at least two elements within the same level. The overall average for the entire dataset is x¯ = −0.16, σ = 0.54, min = −1.39, max = 0.85, indicating a near balance between forelimbs and hindlimbs, though with a slight tendency toward the hindlimb. Some taxa, such as Cervus, lean toward the forelimb (Figures 7B, C), while others, like Equus and Sus, show a preference for the hindlimb.

This distribution is particularly apparent in the density plots (Figures 7C, D), where a clear balance can be observed in taxa such as BB and Stephanorhinus, though with a predominance toward the hindlimb, especially in the case of Equus. Stratigraphically, although there is variability and extreme values in some species like BB in Levels IV or VII, the differences are not statistically significant, even though Equus and Cervus exhibit the highest variability throughout the sequence. Taxa in Group 2 are dominated by the hindlimb at all levels, except for Dama, which shows a predominance of the forelimb in Level XV.

Since Group 2 exceeds 50% missing values, we do not consider it suitable for further analyses, as these cannot be interpreted as absences (value 0) given the current sample size and excavated area. Likewise, data imputation in this case would not provide reliable results, so we focused on Group 1 for a more detailed exploration. Missing values (BB: 15%; Cervus and Equus: 3.8%) were imputed using the “mice” algorithm (m = 5, method = “pmm”, maxit = 50) (van Buuren and Groothuis-Oudshoorn, 2011), which is robust and flexible for such distributions (BB: W = 0.94, p = 0.45; Cervus: W = 0.94, p = 0.46; Equus: W = 0.95, p = 0.545). Levene's test [Pr(>F) = 0.19, p = 0.83] indicates homogeneity of variances among species. The analysis (H = 2.29, df = 2, p = 0.32) suggests no statistically significant differences between the ratios of the three species. However, the correlations between species are low, indicating little direct relationship in terms of the forelimb-hindlimb ratio. This could suggest that their representation across stratigraphic levels is independent or influenced by different factors.

The analysis of ratios across stratigraphic levels (H = 10.82, df = 12, p = 0.54) also indicates no statistically significant differences between levels. The distributions are evenly spread and independent of stratigraphic level. However, while the variability is not statistically significant, it may still provide some stratigraphic insights.

The principal component analysis between species (Figure 8A) reveals a distinct behavior of Equus compared to BB and Cervus. Component 1 explains forelimb representation. As a result, Equus contributes very little to this component, appearing mainly in Levels VII, VIII, and XII. The distribution of BB and Cervus is the most similar throughout the sequence, showing a tendency toward the forelimb but with opposing differences in certain levels, such as XI and XVI. In this context, BB exhibits the greatest balance in forelimb-hindlimb proportions throughout the sequence, while Cervus shows a stronger tendency toward the forelimb (Figure 7C). Component 2, which explains hindlimb representation, is strongly dominated by Equus, contributing hindlimb elements across all levels except the three mentioned earlier. Except for Level VIII, entirely dominated by forelimb elements, the remaining levels tend, to varying degrees, toward the hindlimb.

To better understand the association between species elements and stratigraphic levels, we performed a correspondence analysis (Figure 8B). Most levels cluster near the origin, suggesting a balance between forelimb and hindlimb, though with a tendency toward the hindlimb. Some levels are outliers, such as Level VIII, composed solely of forelimb elements dominated by Equus, and its opposite, Level XVII (imputed data), composed entirely of hindlimb elements dominated by BB.

There does not appear to be a clear pattern in the representation of forelimbs and hindlimbs throughout the sequence. What we can deduce is that Equus is the dominant species, with a strong representation of hindlimbs, while Cervus shows a slight tendency toward forelimbs, often contrasting with Equus at various levels. Considering that BB tends toward anatomical balance, while Equus and Cervus show opposing tendencies, the overall dataset generally reflects a balanced representation between forelimbs and hindlimbs, as indicated by the global mean and statistical tests performed.

A total of 98.6% of the remains included in this study are fragmented, and 97.3% retain less than half of their original size (Table 2). Moreover, the vast majority (91.9%) retain between 0 and 25% of their size. Their distribution across the sequence is very homogeneous (x¯ = 90.3, σ = 6.4, CV = 7.1), with hardly any differences in terms of integrity between stratigraphic levels. However, when analyzing the length of ungulate remains, we do find significant variability (min = 2 mm, max = 283 mm, x¯ = 49.9, σ = 30.2, CV = 60.6). The overall average for all levels is between 4.5 and 6 cm, except for layer V, which has an average length of 23 mm. These results already point to a high degree of fragmentation. In fact, the circumference of long bones is <180° in 94% of the cases (x¯ = 93.2, σ = 8.1, CV = 8.7), and in 60% for those <90° (x¯ = 55.2, σ = 13.6, CV = 8.7). The ratio of fractured remains is 19.9 for each complete element in the entire assemblage. An increase in this ratio is observed from the lower levels to the upper ones, but with some variability in the central levels.

The highest degree of fragmentation is observed in Group 1, with MNE/NISP averages (Table 3) around 0.7, in the following order: Equus, BB, and Cervus. In contrast, Group 2 has values above 0.9: Rhino, Dama, and Sus. Among skeletal categories, cranial elements and long bones are the most fragmented (0.76), followed by girdles (0.79), axial elements (0.8), and phalanges (0.9). Horns and antlers (0.49), along with femora (0.59), exhibit the highest degree of fragmentation, followed by tibiae (0.7) and radii (0.72). Next are skull remains, humeri, scapulae, ribs, and vertebrae (0.74–0.82). A third group consists of pelvises, teeth, metapodials, ulnae, and phalanges (0.83–0.90), and finally, articular bones, which are mostly complete (0.97). If we compare the index solely for long bones, again the horse, BB, and red deer show values between 0.7 and 0.76, followed by the fallow deer (0.86), with wild boar and rhinoceros above 0.96. However, when we analyse the intensity of fragmentation using the %IFOin (Table 2), which considers whole elements, we observe that red deer and fallow deer have the highest values, followed by the horse, BB, and rhinoceros. The wild boar shows the lowest fragmentation intensity, as its anatomical configuration, dominated by teeth and largely complete articular bones, reduces the proportion of fractured elements.

Long bones are the most represented in the assemblage (n = 1,817). Of these, 35.88% have been identified (NISP = 652), with an MNE of 328. Within this group, we found only two complete elements (metatarsals), two complete shafts (humerus and metatarsal), and 26 complete epiphyses (1.6%). The remainder consists of diaphyses (83.1%) and fragmented epiphyses (14.9%). The overall ratio between shafts and epiphyses is 5.56. Among the identified elements, tibiae and femora are the best represented, followed by radii-ulnae and humeri. Metapodials, as a group, are well represented, with metatarsals being the most abundant. More than half of the long bone remains (56.8%) could not be identified. At the stratigraphic level, there is considerable variability in the number of long bone remains, but there is a strong positive correlation (ρ = 0.774, p = 0.001; τ = 0.641, p = 0.001) between the number of shafts and epiphyses throughout the sequence.

Fragmentation was also analyzed for phalanges and mandibles (Table 3). A total of 79 ungulate phalanges were recorded, with first phalanges being the most abundant (50.6%). Of these, 77% are fractured, and 6.3% are longitudinally fractured, primarily first phalanges. Although the index of longitudinally fractured phalanges suggests similarities between first and third phalanges, complete phalanges account for 16.4%, with third phalanges being the most common. The taxa with the highest incidence of longitudinally fractured phalanges are the red deer, BB, and the horse. As for mandibles, the NR is slightly higher than that of phalanges. All mandibles are fractured, and of these, 39.4% are longitudinally fractured, significantly raising the longitudinal fracture index (64.9%). Like the phalanges, the taxa from Group 1 show the highest proportion of this fracture type. At the stratigraphic level, they exhibit the same variability as the other indices.

Regarding skeletal preservation, there is a marked bias against the axial skeleton, basipodia, and acropodia. To a lesser extent, but still with significant bias, the girdles and even the metapodia are affected. In contrast, the regions with the highest survival are the stylopodia, zeugopodia, and cranial skeleton. Long bones are the elements that have been best preserved, particularly the diaphyses. The relationship between diaphyses and epiphyses shows a considerable difference, with a significant bias against epiphyses, especially the distal ones. Differential preservation is evident, but it is not related to bone density and therefore not to post-depositional processes. However, there are other significant taphonomic factors, such as fire, that have undoubtedly contributed to the differential preservation of the assemblage (Solano García, 2024).

The anatomical regions with the highest frequency of cut marks in the assemblage are the axial skeleton and the stylopodia, followed by the zeugopodia and girdles. Cut marks on metapodia and acropodia are present but less frequent, and are almost absent on basipodia. The distribution of cut marks and their various attributes on the represented elements are very similar across all stratigraphic units. Virtually the entire butchery processing sequence is present in all stratigraphic units, except for periosteum scraping, which shows greater variability.

Intentional fracturing evidence is present across all units and species, but occurs most frequently in Group 1. Throughout the sequence, the predominant fracture profiles are longitudinal, with a prevalence of pointed and oblique types. These are followed by spiral and pointed spiral fractures. The predominant fracture angle at all levels is right, except in level IX, where acute angles prevail. Acute angles are the second most common in all units, while mixed angles are less frequent (Solano García, 2024).

The proportion of burnt remains throughout the J/K profile is remarkably high (82.51; Table 4). The analysis of their distribution along the sequence reveals a clear difference between Complex III (x¯ = 82.3%, σ = 10.3%) and Complexes II (x¯ = 96%, σ = 3.3%) and I (x¯ = 96.1%, σ = 3%), which share the same proportion. Therefore, from Level XIII onwards, an increase in thermally altered elements is observed, persisting up to the top of the sequence.

In this analysis, all bone remains showing cut marks, green fractures, percussion marks, and evidence of burning were recorded (Table 4). Additionally, remains with carnivore tooth marks were included as complementary information. From this count, we observe that 7.62% of the remains bear cut marks, 12.9% show green fractures, 1.58% display percussion marks, and 82.5% of the remains exhibit evidence of fire. Carnivore tooth marks were recorded at a relatively low rate (0.79%) across the entire J/K assemblage. At the stratigraphic level, as in previous analyses, there is variability in the number of remains with anthropogenic modifications per level, due to the density of remains per excavated surface. However, there is a strong correlation between the various types of modification and the total number of remains without teeth for each stratigraphic level (cutmarks: ρ = 0.877, p < 0.01; percussion: ρ = 0.822, p = 0.002; burnt: ρ = 0.972, p < 0.01; fractures: r_s = 0.832, p < 0.01), except for the few carnivore tooth marks (r_s = 0.557, p < 0.06). Therefore, the results are consistent and proportional throughout the stratigraphic sequence. In other words, where there is a higher number of remains, it is reasonable to expect a higher number of anthropogenic modifications. Thus, given that the highest number of elements is found in the central levels, it is also here where the greatest number of modifications is located.

Nearly 90% of the remains are burnt, with greater representation in levels XVII–X and level IX. There is no preferential burning pattern by taxon, as around 90% of their remains also show signs of burning. However, this is not the case for the rest of the surface modifications. As expected, the taxa from Group 1 show the highest incidence of anthropogenic modifications. When we analyse the relationship between the six main taxa and cut marks, percussion marks, fractures, and carnivore tooth marks, there is a statistically significant association (X² = 25.124, df = 15, p = 0.04831), given the strong association within Group 1 (Figure 9A). In fact, Cervus and Equus show a strong association with cut marks and fractures, while BB is more associated with fracturing (19.8%). The red deer has a stronger relationship with percussion (24%), cut marks (20.5%), and fracturing (18.4%). The horse shows an almost equal association (~60%) with percussion, fracturing, and cut marks. Interestingly, fallow deer shows a low incidence of inflicted modifications, distancing it from the main group, while the rhinoceros shows some association primarily with fractures. In the case of the wild boar, with only two recorded cut marks and one percussion mark, it stands out as an atypical species within the assemblage. Carnivore activity is mostly associated with Cervus (50%) and Equus (30%).

Anatomical categories show a strong association with surface modifications (X² = 92.0623, df = 21, p = 7.0971e-11). Long bones are highly associated with percussion activities (92%), fracturing (90.3%), and cut marks (59.8%), while axial elements (18.2%), girdles (10.6%), and mandibles (6.8%) are primarily associated with cut marks (Figure 9B). Phalanges show a stronger correlation with percussion (8%) and fractures (3.7%) than with other modifications. Finally, cranial elements are distanced from the main group, with a low connection to cut marks (0.8%) and fractures (0.5%). Carnivore tooth marks are predominantly associated with long bones (50%) and girdles (20%) from red deer and horse.

In summary, the presence of the three main species throughout the sequence reflects the exploitation of meat, marrow, and fat, as indicated by the fracturing and its correlation with utility indices, and the systematic distribution of predominant cut marks on ribs and long bone diaphyses (Solano García, 2024).