Identifying the unidentified fauna enhances insights into hominin subsistence strategies during the Middle to Upper Palaeolithic transition

Understanding Palaeolithic hominin subsistence strategies requires the comprehensive taxonomic identification of faunal remains. The high fragmentation of Late Pleistocene faunal assemblages often prevents proper taxonomic identification based on bone morphology. It has been assumed that the morphologically unidentifiable component of the faunal assemblage would reflect the taxonomic abundances of the morphologically identified portion. In this study, we analyse three faunal datasets covering the Middle to Upper Palaeolithic transition (MUPT) at Bacho Kiro Cave (Bulgaria) and Les Cottés and La Ferrassie (France) with the application of collagen type I peptide mass fingerprinting (ZooMS). Our results emphasise that the fragmented component of Palaeolithic bone assemblages can differ significantly from the morphologically identifiable component. We obtain contrasting identification rates between taxa resulting in an overrepresentation of morphologically identified reindeer (Rangifer tarandus) and an underrepresentation of aurochs/bison (Bos/Bison) and horse/European ass (Equus) at Les Cottés and La Ferrassie. Together with an increase in the relative diversity of the faunal composition, these results have implications for the interpretation of subsistence strategies during a period of possible interaction between Neanderthals and Homo sapiens in Europe. Furthermore, shifts in faunal community composition and in carnivore activity suggest a change in the interaction between humans and carnivores across the MUPT and indicate a possible difference in site use between Neanderthals and Homo sapiens. The combined use of traditional and biomolecular methods allows (zoo)archaeologists to tackle some of the methodological limits commonly faced during the morphological assessment of Palaeolithic bone assemblages.


Introduction
The investigation of behavioural shifts in prey selection and hunting strategies during phases of major changes in the material record is key to understanding the relationship between human behavioural evolution, cultural variation, and population dynamics (Delagnes and Rendu 2011;Discamps et al. 2011;Niven et al. 2012;Rendu et al. 2012;Steele 2004). Traditionally, such behavioural patterns have been approached through the analysis of the stone tools and faunal remains recovered from excavations at Palaeolithic sites. In particular, faunal specimens provide the opportunity to identify and document behaviour developed by human populations for the exploitation of their environment (Gaudzinski and Roebroeks 2000;Gaudzinski-Windheuser et al. 2014;Morin 2012;Pederzani et al. 2021;Smith et al. 2021;Stiner 1993). However, studying ancient fauna not only provides paleoenvironmental information, but when combined with the analysis of bone surface modifications related to human activity, it can finetune the timing of human occupations and helps to reconstruct human diet and interactions with other groups or even Extended author information available on the last page of the article species (Steele 2015). Indeed, faunal exploitation is related to a wide range of behaviours and cognitive aspects entwined with mobility, social organisation, technological development, and subsistence capacities (Marean and Assefa 1999).
However, Late Pleistocene bone assemblages are often highly fragmented, preventing proper taxonomic identification and anatomical attribution of many specimens based on morphology alone (Lyman 2002;Morin et al. 2017a). Several processes affect faunal remains, starting from decomposition, selective destruction in the sediment, post mortem transport, and burial, to preserve bone specimens that are potentially altered during excavation, cleaning treatment, and storage (Lyman 1994;Marean 1991). All these factors, geological, biological, and cultural, can lead to variability in faunal identification. Together with differential preservation, they can create a potential source of bias for the interpretation and quantification of relative abundances of taxa (Dirrigl 2002;Marean 1991;Marean and Kim 1998;Morin 2004;Morin et al. 2017aMorin et al. , 2017bPickering et al. 2006). Indeed, combined with the impact of human and carnivore activities at the site, these factors contribute to reduced proportions of taxonomically diagnostic bones resulting in a lower number of identifiable specimens. Such processes generate the potential to seriously distort various archaeological and ecological inferences (Faith 2007;Morin et al. 2017a).
Recent developments of biomolecular methods allow us to exploit the collagen preserved in these bone fragments to taxonomically identify faunal specimens (Buckley et al. 2009). The inclusion of the analysis of highly fragmented bone through proteomic screening using zooarchaeology by mass spectrometry (ZooMS) for the taxonomic assessment of Palaeolithic faunal assemblages has already demonstrated its great potential (Berto et al. 2021;Brown et al. 2021b, a;Buckley et al. 2017;Pothier Bouchard et al. 2020;Ruebens et al. 2022;Sinet-Mathiot et al. 2019;Welker et al. 2015) and highlighted the necessity to use a multi-methodological approach in studying human subsistence. Taxonomic identities from both the morphologically identified and the ZooMS-identified components can thus be correlated with bone surface modification analysis to address specimen surface preservation and bone accumulation agents through the reconstruction of taphonomic history. The analysis of the protein type I collagen provides taxonomic identity based on variations in amino acid sequences and allows the taxonomic identification of bone assemblages to be extended to the morphologically unidentifiable component. The previous application of ZooMS as a screening tool for faunal assemblages has provided variable results in terms of the comparability of the two components. Taxonomic abundances of the morphologically unidentifiable component of a faunal assemblage may not generally differ from the morphologically identifiable component (Berto et al. 2021;Buckley et al. 2017;Welker et al. 2016Welker et al. , 2017, but that does not necessarily indicate a pattern (Ruebens et al. 2022;Sinet-Mathiot et al. 2019). Moreover, such differences could reflect a specific human behavioural signature related to bone fragmentation and intensity of carcass processing (Sinet-Mathiot et al. 2019). A better understanding of the source(s) of variability will help in anticipating the potential differences that may occur within certain bone assemblages.
The zooarchaeological literature frequently contains body size class attributions of bone specimens that cannot be reliably assigned to a particular taxon or clade. It is generally assumed that these body size class attributions are reliable and reflect or contain taxonomic information about the bone assemblage as a whole. However, previous ZooMS research has highlighted that this is a potentially unreliable approach (Sinet-Mathiot et al. 2019). Here, we test the fragmentary component of bone assemblages of three Late Pleistocene sites: Bacho Kiro Cave (Bulgaria) and Les Cottés and La Ferrassie (France). They all show rich and well-preserved stratigraphic sequences spanning the Middle to Upper Palaeolithic transition (MUPT). These sites offer the opportunity to discuss diachronic changes in subsistence strategies during the period of possible interaction between Neanderthal and Late Pleistocene Homo sapiens populations (Hajdinjak et al. 2021;Higham et al. 2014;Hublin 2015;Hublin et al. 2020;Prüfer et al. 2021). This work explores the implications of incorporating the analysis of the morphologically unidentifiable bone component into the description of faunal assemblages in terms of both overall bone accumulation and aims to advance our interpretation of human subsistence strategies during the MUPT. We address methodological limits commonly faced during the morphological assessment of faunal assemblages and demonstrate how the addition of biomolecular methods, such as untargeted ZooMS screening, can complement our understanding of subsistence behaviour by providing a clearer picture of prey selection and site occupation. By including assemblages that span the MUPT in Europe, we are thereby able to demonstrate that the assessment of the fragmented component of bone assemblages through ZooMS can provide different patterns of species frequencies than previously interpreted based solely on the morphologically identifiable record.

Sample selection
This study includes the ZooMS analysis of bone material from three Late Pleistocene sites (Bacho Kiro Cave, Les Cottés, and La Ferrassie; SI Fig. 1, SI Table 1). All the material taxonomically identified through bone morphology by zooarchaeologists will be referred to as the morphology component. Similarly, all fragmentary specimens morphologically unidentifiable and taxonomically identified through ZooMS will be referred to as the ZooMS component. All three sites were recently excavated and have provided large, wellcontextualised, and highly fragmented bone assemblages of individually piece-provenienced faunal remains. Bone surface analyses of both the morphologically identified and the fragmentary unidentifiable bone assemblages were assessed using comparable zooarchaeological methods and protocols. All faunal data were derived from recent excavation campaigns, and specimens from both the morphology and ZooMS components show similar spatial distributions over the excavated areas. As the fragment size cutoff of 20 mm was used for piece-proveniencing artefacts during the excavation of the sites and in order to permit eventual additional biomolecular analysis, fragmentary and morphologically unidentifiable piece-provenienced specimens generally > 20 mm in length were selected for proteomic analysis. Indeed, the zooarchaeological analysis of fragments smaller than 2 cm is commonly limited as taphonomic attributes are size dependent. Bone material resulting from sediment sieving during the excavation of the archaeological sites is not included in this study. All morphologically unidentifiable piece-provenienced specimens from the La Ferrassie layer 6 faunal assemblage were selected for ZooMS analysis. In the case of Bacho Kiro Cave and Les Cottés, specimens were randomly selected by the zooarchaeologists among the morphologically unidentifiable material without any selection towards bone morphology nor specific surface modifications.

Bacho Kiro Cave
Bacho Kiro Cave (Dryanovo, Bulgaria) is located on the northern slope of the Balkan mountain range (Stara Planina) and about 70 km south of the Danube River. Previously investigated during the twentieth century (Garrod et al. 1939;Kozłowski and Ginter 1982), the site was reopened for excavation in 2015 by the National Archaeological Institute with Museum from the Bulgarian Academy of Sciences (Sofia, Bulgaria) and the Max Planck Institute for Evolutionary Anthropology (Leipzig, Germany). The archaeological sequence spans the Middle Palaeolithic (MP) through to the Upper Palaeolithic (UP). The archaeological material recovered from two sectors (Main Sector and Niche 1) from layers I and J was recognised as part of the Initial Upper Palaeolithic marked by the earliest occurrence of Late Pleistocene Homo sapiens in Europe . This starts around 45,990 cal BP in the upper part of layer J and considerably intensifies in layer I which is dated to 45,040-43,280 cal BP Pederzani et al. 2021). This material comprises the earliest and largest number of Homo sapiens bone tool and ornament assemblages in Europe, partly taxonomically identified through ZooMS (Martisius et al. 2022). The assemblage recovered from layer K was technologically associated with the MP and was deposited between 61 ± 6,000 and 51,000 year BP Pederzani et al. 2021). We investigated 1,595 faunal remains through ZooMS from layer I (814 specimens), layer J (438 specimens), and layer K (343 specimens) from both the niche 1 and the main sector . Zooarchaeological analysis was performed on 7,013 faunal remains from layers I, J, and K from both sectors following previously described methodology ) and including 1,453 specimens assigned to a taxonomic group (1,077 from layer I, 232 from layer J, and 143 from layer K).

Les Cottés
Les Cottés (Vienne, France) is a cave located in the corridor between the Parisian basin and the Poitou in westcentral France. The site was discovered in 1878 and was explored through several excavation campaigns (Bastin et al. 1976;Lévêque 1997;Pradel 1967), but the material included in this study is derived from an excavation initiated in 2006 by M. Soressi with support of the French Ministry of Culture and Communication and the Max Planck Institute for Evolutionary Anthropology (Soressi et al. 2010). Through ZooMS, we analysed 523 morphologically unidentifiable faunal specimens, which, together with the 152 presented in Welker et al. (2015) (137 undiagnostic fragments and 15 specimens analysed in a ZooMS blind test), means 675 specimens from Les Cottés were analysed with ZooMS. Of these, 220 are from the Mousterian (US08, dated between 46,051 and 42,034 cal BP using radiocarbon and between 55 and 48 ka according to the OSL measurements (Jacobs et al. 2015)), 217 are from the Châtelperronian (US06, dated between 42,961 and 40,344 cal BP), 168 are from the Protoaurignacian (US04 lower), and 70 are from the Early Aurignacian (US04 upper). The dates for the Aurignacian layers extend from 40,372 to 36,697 cal BP (Talamo et al. 2012) in radiocarbon years, or from 43 to 36 ka in OSL years (Jacobs et al. 2015). Interpretations coming from US04 upper are considered with caution due to the low number of specimens in comparison to the other layers. Bone surface analysis was standardised over the assemblage and was previously described elsewhere . Of a total of 5,169 bone remains assessed through traditional zooarchaeology, 1,922 bone and dental specimens were morphologically identified in the range of subfamily to species (397 specimens from US08, 166 from US06, 715 from US04 lower, and 629 from US04 upper).

La Ferrassie
The Grand Abri of La Ferrassie (Savignac-de-Miremont, France) is in the Dordogne region of south-western France in a tributary valley to the Vézère River and was first excavated during the twentieth century by Capitan and Peyrony and then by Delporte (Delporte and Delibrias 1984;Peyrony 1934). An excavation conducted from 2010 to 2015 by Turq and colleagues further refines the stratigraphic sequence spanning the MUPT (Guérin et al. 2015;Turq et al. 2012). The Châtelperronian layer (layer 6) was dated to between 45,100 and 39,520 cal BP ) marking the earliest appearance of this lithic industry in the region. The faunal material from this layer that was morphologically identifiable to taxon is limited to 17.5% of the bone assemblage (142 specimens) and is dominated by reindeer (Rangifer tarandus). All piece-plotted, morphologically indeterminate specimens were processed through ZooMS (527 specimens).

ZooMS methodology
ZooMS extraction protocols employed for this study were partially described previously (Buckley et al. 2009;van Doorn et al. 2011;Welker et al. 2016). All 2,645 specimens were sampled (10-30 mg) using pliers and placed into 96-well plates. Soluble collagen was extracted through incubation in 100 µL of 50 mM ammonium bicarbonate (AmBic) buffer at 65 °C for 1 h. In order to improve and verify the taxonomic identity obtained from soluble collagen, 440 specimens (70 for La Ferrassie, 369 for Les Cottés, and 1 for Bacho Kiro Cave) (SI Table 2) were demineralised in 130 µL 0.6 M HCl at 4 °C for 18-20 h, neutralised to pH 7, and solubilised again in AmBic. Then, 50 µL of the resulting supernatant was digested using trypsin (0.5 µg/µL, Promega) overnight at 37 °C, acidified using trifluoroacetic acid (20% TFA), and then cleaned on Hypersep C18 96-well plates (Thermo Scientific) using a vacuum manifold. In short, a 96-well deepwell plate (Eppendorf) is placed beneath the Hypersep plate to collect the solutions. C18 filter tips from the Hypersep plate were conditioned with 200 µL of 0.1% TFA in 50:50 acetonitrile and UHQ water (conditioning solution) and washed with 200 µL of 0.1% TFA and UHQ water (washing solution). Peptide extracts were then vacuumed through the filters slowly to ensure optimal binding efficiency. The obtained waste solution was discarded. Filters were then washed again with 200 µL of washing solution, and peptides were extracted in 100 µL of conditioning solution and transferred to a 96-well plate. Digested peptides were spotted in triplicate on a MALDI Bruker plate with the addition of α-cyano-4-hydroxycinnamic acid (CHCA, Sigma) matrix, using a multichannel pipette (Thermo Fisher).
MALDI-TOF-MS analysis was conducted at the Fraunhofer IZI in Leipzig (Germany), using an autoflex speed LRF MALDI-TOF (Bruker) in reflector mode, positive polarity, and matrix suppression up to 590 Da and collected in the mass-to-charge range 1000-3500 m/z. Triplicates were merged for each sample in R version 4.0.5 (R Core Team 2021) and MALDIquant v. 1.21 (Gibb and Strimmer 2012). First, we smooth the intensity using a moving average and remove the baseline using the TopHat approach. Then, for each sample, we align the replicate spectra using SuperSmoother and a signal to noise ratio of 3, sum the three replicates to obtain a single spectrum, and remove the baseline once more, again using TopHat. Spectra were exported as.msd files. Taxonomic identifications were made using mMass 5.5.0. (Strohalm et al. 2010) through manual peptide marker mass identification in comparison to a database of peptide marker series for all European Pleistocene mediumto large-sized mammals (Welker et al. 2016). To assess any potential contamination by non-endogenous peptides, we performed laboratory blanks alongside the samples. These remained empty of collagenous peptides, excluding the possibility of modern laboratory or storage contamination.
Peptide marker series can be similar for some closely related species, which is the case for the species belonging to the following taxonomic groups: Bos/Bison, Cervid/ Saiga, Equidae, and Ursidae. Cervid/Saiga can be attributed to either Cervus elaphus (red deer), Megaloceros giganteus (giant deer), Alces alces (elk), or Dama sp. (fallow deer). Equidae and Ursidae include, respectively, species from the genera Equus and Ursus, most likely Equus ferus and Equus hydruntinus or Ursus spelaeus and Ursus arctos. In order to facilitate the comparison between ZooMS and morphology components, the most common species and taxa were grouped into broader categories of Bos/Bison (Bos primigenius, Bison priscus, and Bos/Bison sp.), Cervid/Saiga (C. elaphus, D. dama, M. giganteus), Ursidae (U. arctos, U. spelaeus, and Ursus sp.), Capra sp. (C. ibex and Capra sp.), and Equidae (E. ferus, E. hydruntinus, and Equus sp.) (SI Table 3). At Bacho Kiro Cave, due to the high proportion of this taxonomic group, Cervid specimens from the morphology component were also included into the broader group Cervid/Saiga. Within the ZooMS component, the few specimens identified as Cervid/Saiga/Capreolus capreolus were included in the broader taxonomic group Cervid/Saiga, in order to allow the comparability of both components.
Suggested as an indicator of collagen preservation (Welker et al. 2017;Wilson et al. 2012), glutamine (Gln) deamidation ratios were calculated on all samples for peptide COL1α1 508-519 (Brown et al. 2021b, a), which is frequently observed in peptide fingerprints of collagen type I, following published protocols Wilson et al. 2012). The deamidation ratio ranges from %Gln = 1.0 with non-deamidated glutamines to %Gln = 0.0 indicating a full deamidation of the glutamines. The Gln deamidation ratios obtained during routine ZooMS screening have been previously suggested to assess bone assemblage homogeneity (spatial and temporal variability within a site), to detect stratigraphic outliers (intrusive material or differential bone preservation), to inform on the preservational quality of specific peptides and specimens, or to look at the taxonomic distribution from a biomolecular perspective (Sinet-Mathiot et al. 2019;Welker et al. 2017;Wilson et al. 2012), although with varying success (Brown et al. 2021b, a).

Zooarchaeological and taphonomic methodology
All taphonomic modifications were recorded on the morphology-and ZooMS-identified specimens by the respective zooarchaeologists, consistently within and between studied sites. Bone surfaces for both the morphologically identified and unidentified components were assessed through visual inspection, using magnification when needed (up to × 20 magnification) (Blumenschine et al. 1996). The maximum length of the bone specimens was measured individually with digital calipers.
Although traces of burning were recorded during taphonomic analyses using the scale proposed by Stiner and colleagues (Stiner and Kuhn 1995) (0: Unburnt to 6: Completely calcined), these burnt remains were excluded from subsequent ZooMS analyses due to poor collagen preservation. Weathering stages were recorded for all bones and provide a qualitative scale for understanding the exposure (short/long duration) of the bones prior to burial (Behrensmeyer 1978). A slightly modified scheme was used on Les Cottés bone assemblages where specific modifications were recorded related to weathering (see Rendu et al. (2019)). Specifically, weathering was recorded according to three variables: exfoliation (the peeling of bone surface), cracking (the emergence of longitudinal cracks on bone surface), and disintegration (the complete destruction of the bone). In addition, other recorded modifications included root etching and abrasion (expressed as a percentage of bone surface affected). The schemes range from 0% (no visible modification observed) through 100% (the whole bone surface covered) Behrensmeyer 1978;Blumenschine et al. 1996;Domínguez-Rodrigo et al. 2017;Fisher 1995;Lee Lyman 1994;Olsen and Shipman 1988;Soulier and Costamagno 2017).
For all three bone assemblages, human modifications included traces related to butchery and carcass processing (cut marks, scraping marks, chop marks, marrow bone breakage), and carnivore modifications included tooth marks, gnawing traces, and damage from bone breakage and digestion as well as rodent tooth traces. The number of identified specimens (NISP) represents the number of specimens assigned to a taxon.
When it was not possible to morphologically assign fragmentary bone specimens to a specific taxon, these were assigned to a specific body size class based on previous assignments (Morin 2012). The separation of specific taxa into different body size classes was normally done on the basis of both body and skeletal size (following (Morin 2012;Rendu et al. 2019;Smith et al. 2021)).
The combination of the ZooMS and morphology component allows for the assessment of skeletal element distributions and the possible identification of previously unrecognised skeletal elements, which has implications for our understanding of hunting strategies and carcass transport. As skeletal elements were identified, when possible, on taxonomically unidentifiable specimens, we aimed to correlate the skeletal part identifications with the ZooMS taxonomic identities in order to assess skeletal representation among both the morphology-and ZooMS-identified components. To assess skeletal element representation for the dominant taxa within each component, bone elements were categorised into body parts for each method of identification (cranial: cranium, mandible; axial: vertebrae, pelvis, rib; forelimb: humerus, radius, ulna; hindlimb: femur, tibia; distal limbs: carpals, metacarpal, tarsals, metatarsal, phalanges; LBN: long bone fragments, FBN: flat bone fragments) (based on Stiner (1991aStiner ( , 1991b). Within all three datasets, teeth and antler were categorised separately from the cranial body part, as their inclusion might bias the comparison between components. Indeed, antlers and horn cores tend to be rare and are more easily identified morphologically, reducing their representation in the ZooMS component. Anatomically unidentifiable specimens (NID) were excluded from the assessment of skeletal element representation as they did not provide substantial information.
Ecological diversity indices were calculated in order to investigate the effect of the addition of the ZooMS-identified specimens on the taxonomic diversity of the faunal community of each layer and site. We used the Shannon-Wiener index (H′) (Shannon 1948) and R package vegan v. 2.6-2 (Oksanen et al. 2019) to quantify the taxonomic diversity of our three faunal assemblages between the ZooMS and morphology component and to assess any variation in faunal diversity between layers of each site, taking into account the taxonomic richness and the distribution of their abundance. As the Shannon-Wiener index is sensitive to sample size, values should be considered with caution when the sample size is small. Along with species richness, Pielou's evenness (J′) measures taxonomic diversity by giving the count of individuals of each taxonomic group among each component and reflecting the evenness of the distributed abundances between taxa. The index value ranges from 0 (no evenness) to 1 (complete evenness).

ZooMS analysis
ZooMS analysis of all three datasets shows well-preserved collagen type I with a high success rate of taxonomic identification, up to the range of subfamily or genus, between 90 and 97% (SI Table 1). For 82% of the samples, the semi-destructive extraction protocol (AmBic) is sufficient to obtain a ZooMS identification. At Bacho Kiro Cave, collagen preservation is excellent (also noted by Fewlass et al. (2020)) resulting in a high proportion of discriminant taxonomic identities. All ZooMS samples could be extracted using only the AmBic protocol, while we extracted one specimen through acid demineralisation as well to verify its taxonomic identity (Castor fiber). At both Les Cottés and La Ferrassie, samples were processed using both AmBic and acid demineralisation protocols to improve and optimise taxonomic identifications (SI Table 2).

Deamidation between stratigraphic units and taxa
Glutamine deamidation ratios are calculated in order to detect potential intrusive material between archaeological layers or differential collagen preservation between taxa. Because the data is not normally distributed (Shapiro-Wilk normality test, p-value < 0.05), we used Wilcoxon-Mann-Whitney tests to compare the glutamine deamidation ratios between taxa and layers. The deamidation ratios presented here are derived from samples analysed using the acid-free protocol. At Bacho Kiro Cave, we observe that older samples from layer K show elevated levels of deamidation with values significantly different between layers (SI Table 5, SI Fig. 2). In contrast, we note overlapping deamidation values between layers at Les Cottés (SI Fig. 3) with the exception of US06 which showed values significantly different from US04 lower and US08 (SI Table 5). Glutamine deamidation ratios seem to overlap between dominant taxa which would suggest that they have undergone similar molecular diagenetic processes within each site (SI Fig. 4). However, a few exceptions could be identified. At Bacho Kiro Cave, deamidation ratios show similarities between taxa, particularly within layer K, but ursid specimens tend to have glutamine deamidation values significantly different from other taxa in layers I and J, notably in comparison with Bos/Bison, Capra sp., and Equidae (SI Table 6). At Les Cottés, all taxonomic groups show similar deamidation ratios within each layer, except for a few Rangifer tarandus specimens (n = 6) showing deamidation values significantly different from Equidae in US04 lower (Wilcoxon-Mann-Whitney tests, statistic = 64, p-values = 0.013, SI Table 6). At La Ferrassie, Rangifer tarandus and Cervid/Saiga specimens show deamidation values significantly different from Bos/Bison (SI Table 6). The statistical differences observed between some of the taxonomic groups and layers could be driven by discrepancies in sample sizes, taphonomic history, and site formation or butchery practices. However, further exploration is required in order to interpret these differences.

Taxonomic representation
Species representation among both ZooMS and morphology components are generally consistent within each site, but the addition of ZooMS permits the identification of taxa that were unrecognisable through morphology. At Les Cottés, ZooMS identified Felidae and Ursidae in the faunal community obtained from US06 but also resulted in the addition of Cervid/Saiga in US04 lower (SI Table 7). At Bacho Kiro Cave, the ZooMS analysis allowed for the identification of Elephantidae in layer J (SI Table 8). At La Ferrassie, the use of ZooMS results in a fourfold increase of the number of taxonomically identified specimens. Consequently, the taxonomic diversity for this layer was broadened, with the addition of Capra sp., Elephantidae, Rhinocerotidae, Ursidae, and several carnivores (SI Table 9).
Shannon-Wiener index calculations show that the diversity of the faunal community identified on a site can significantly change with the addition of ZooMS. More specifically, we observe an increase in the faunal diversity of the combined ZooMS and morphology components in the layers under study here at La Ferrassie, at Les Cottés, and at Bacho Kiro Cave (Fig. 1, SI Table 10). In contrast, the lower values of the Shannon-Wiener index, after the addition of ZooMS identities to the layer I faunal assemblage, at Bacho Kiro Cave indicate a lower taxonomic diversity. Such a pattern possibly emphasises a better identification rate within the morphology component related to a larger sample size or highlights a higher evenness of the ZooMS component due to the repeated identification of taxa showing a low abundance among the morphology component.
The occurrence of the dominant taxa, i.e. the taxa showing the highest proportions, among both components are consistent within each site (Bacho Kiro Cave: Ursidae, Equidae, Cervid/Saiga, Capra sp., and Bos/Bison; Les Cottés: Bos/Bison, Equidae, and Rangifer tarandus; La Ferrassie: Bos/Bison, Cervid/Saiga, Equidae, and Rangifer tarandus), but we observe differences in their relative contributions to the overall bone assemblage (Fig. 2). At Les Cottés and La Ferrassie, the ZooMS component indicates lower proportions of reindeer, offset by higher proportions of Bos/ Bison and Equidae (SI Table 7, 9 and 11). We note a ninefold increase in the proportion of Bos/Bison at La Ferrassie. At  SI Table 10 for details). Confidence intervals (2.5-97.5%) are given for each value Les Cottés, we observe an on average twofold increase of Bos/Bison and Equidae with the addition of ZooMS to the analysis of the faunal assemblage. At Bacho Kiro Cave, and similar to Les Cottés and La Ferrassie, Bos/Bison remains are slightly more abundant within the ZooMS component, particularly in layers I and J. Conversely, Ursidae show a similar pattern as reindeer at La Ferrassie and Les Cottés with slightly lower proportions notably in layers J and K. We note a large difference between methods of identification for Cervid/Saiga in layer I, but these differences are not consistent throughout the other layers. When comparing the faunal composition between layers to assess any changes or shifts in the NISP of different taxa across the MUPT, we note at Les Cottés a progressive decrease in the proportions of Bos/Bison offset by an increase of Equidae from US08 to US04, which is particularly clear through the use of ZooMS. Despite the low number of specimens analysed through ZooMS from the Early Aurignacian layer (US04 upper) of Les Cottés, the results obtained show a continuous pattern with those from the layers below in terms of taxonomic abundances between the dominant taxa.
While the categorisation of morphologically unidentifiable specimens into body size classes remains a useful tool when no other alternative is available for the interpretation of this component of the assemblage, the correlation between taxonomic identifications provided by ZooMS with the body size classes indicates inconsistencies. Therefore, the observations made previously at Fumane Cave therefore do not seem to be an exception, but rather the norm Sinet-Mathiot et al. 2019). We observe inconsistencies between body size class attributions, which are largely based on bone size and cortical thickness, and ZooMS taxonomic assignments (SI Table 4). For example, Ursidae specimens are present in most carnivore and ungulate body size class units (Fig. 3), several equid specimens are categorised among the large carnivore class and Caprinae and Capra sp. among the large ungulate class (Fig. 3). Although many zooarchaeologists are already using alternative nomenclatures (i.e. mammal classes or unknown instead of ungulate or carnivore classes (Castel 2011)) or standardisation tools (Discamps 2021), these results simply confirm that body size class attributions should be used with caution, especially when translating these classes to more specific taxonomic units and/or assessment of hominin subsistence strategies. When assigning bone specimens to generalised family attributions, one should cautiously avoid "taxonomic blindness" based on presumed abundance of cladistic assignments that are based on the thickness of the cortical bone.

Bone length distribution
As expected, larger bone fragments are generally more identifiable through comparative morphology as they often preserve more morphologically distinctive features. Smaller fragments tend to be identifiable only through ZooMS (Fig. 4). This pattern is particularly noticeable at Bacho Kiro Cave in layers I and K (SI Fig. 5). However, this is not the case for all taxa. We note a different specimen length distribution between both ZooMS and morphology components among dominant taxa. At Bacho Kiro Cave, Bos/Bison, Cervid/Saiga, and Equidae specimens show an opposite bipolar distribution of their specimen length whereas the two distributions are more similar for Capra sp. and Ursidae. Because the data is not normally distributed (Shapiro-Wilk normality test, p-value < 0.05), we used Wilcoxon-Mann-Whitney tests to compare the bone length distribution between taxa, layers, and method of identification. Bone specimens identified as Capra sp. and Ursidae through ZooMS show Carn, carnivore; Ung, ungulate. Each taxon is assigned a colour to help the visualisation of the graph. Numbers on the bars are the NISP per category a fragment length distribution significantly different from other taxonomic groups particularly in layer I (SI Table 12). Likewise, observations of the bone assemblage from Les Cottés indicate a similar trend with Bos/Bison and Equidae most often exhibiting opposite distributions, compared to the similar distributions of both ZooMS and morphology components for the reindeer specimens (SI Fig. 6). At Les Cottés, the bone length distribution of specimens identified morphologically as reindeer is significantly different from the Bos/Bison and of Equidae distributions in US04 and US08 (SI Table 13), but no differences are observed among the ZooMS component. When comparing the distribution between methods of identification, we also note significant differences for Bos/Bison and equid specimens in US04 and US08 (SI Table 14). The absence of metric measurements on the morphologically identified component from La Ferrassie prevents comparisons of bone length distribution between the ZooMS and morphological components. However, the ZooMS component represents 82.5% of the faunal assemblage, so a comparison of specimen length between dominant taxa for the ZooMS component is possible. Although specimens from the dominant taxa generally show similar length distributions, with a large proportion within the 2-3 cm range, equid bones tend to have fewer large fragments illustrated by a higher proportion of specimens within the smaller size classes (SI Fig. 7). Equid fragments identified through ZooMS present a length distribution significantly different from Bos/Bison, Cervid/Saiga, and reindeer (SI Table 15), most likely due to an overrepresentation of equid fragments of 2-3 cm counterbalanced by an underrepresentation of specimens of 3-4 cm. Nevertheless, it should be noted that Equidae is the taxa with the smallest sample size, which might influence these results.

Bone surface preservation
We investigated readability of the bone surfaces to rule out bone fragmentation related to environmental taphonomic factors. We find that, at Bacho Kiro Cave and La Ferrassie, the bone surfaces of specimens taxonomically identified both through ZooMS and morphology are only affected by low degrees of surface weathering, which cannot explain the differences in fragmentation between taxa and/or layers (SI Table 16 and 17).
Due to high stages of weathering at Les Cottés, many bone surfaces from the ZooMS component exhibit natural fractures. In particular, a large percentage of Bos/Bison fragments, from US06 and US08, and equid specimens in US04 indicate multiple types of surface damage (SI Fig. 8). These patterns are also recorded on reindeer at a high percentage (> 50% for US04 and US06) within the morphology component. The readability of the surfaces, which reflects how bones were affected by weathering or other factors possibly leading to fragmentation, is generally better for the reindeer specimens compared to bones from Bos/Bison and equids (SI Fig. 9).

Bone assemblage accumulator
We investigated bone modifications associated with carnivore and human activity to identify the accumulator agents of the bone assemblage. We find that, within all three sites, ZooMS analysis allows for improved association of taxonomic identity with taphonomic data, which in several cases provides additional behavioural information. Overall, the inclusion of ZooMS identifications within zooarchaeological analyses highlights a diverse range of taxa exhibiting bone modifications from carnivore and human activity (SI Figs. 10 and 11). These results are particularly informative for layers J and K at Bacho Kiro Cave and at La Ferrassie with the addition of three to four taxa previously unassociated with the modifying agents (humans).
With the addition of ZooMS, carnivore modifications were identified at Bacho Kiro Cave on Cervid/Saiga (layer J: 3% NISP, layer K: 23% NISP) and Capra sp. (layer J: 11% NISP, layer K: 11% NISP) and on Bos/Bison (31% NISP) and Equidae (31% NISP) in layer K. Carnivore modifications within the ZooMS component of layer K affected 21% of the remains from the dominant taxa, a considerably higher percentage than previously obtained through morphology (SI Fig. 12). At La Ferrassie, the proportion of carnivore activity within layer 6 appears relatively low compared to human activity as carnivore modifications were identified on only two Bos/Bison specimens within the ZooMS component (SI Fig. 12).
In addition to evidence of carnivore activity, anthropogenic modifications are also present on most taxa within all studied layers. Human modifications were recognised on equids (20% NISP) and Capra sp. (22% NISP) in layer J at Bacho Kiro Cave (Fig. 5), and we noted a relatively high proportion of percussion marks on Cervid/Saiga specimens from layer J (22% NISP) (SI Fig. 11 and SI Fig. 15). At La Ferrassie, human activity is identified on Cervid/Saiga (4% NISP) and Bos/Bison (6% NISP) but not on equid specimens, and percussion traces occur on a higher proportion of reindeer remains (9% NISP; Fig. 5 and SI Fig. 13). At Les Cottés, human modifications range between 10 and 20% among dominant taxa over all studied layers and occur at higher proportions on reindeer specimens (particularly in US04 and US06), mainly represented by cut marks and percussion traces ( Fig. 5 and SI Fig. 14). At Bacho Kiro Cave and Les Cottés, we note a progressive reduction of carnivore activity from the Late Middle Palaeolithic to the Upper Palaeolithic alongside an increase of human modifications at Bacho Kiro Cave, reinforcing patterns previously described Smith et al. 2021). In addition, we note the recurrent occurrence of anthropogenic modifications on carnivore remains (n = 93) from various taxa such as canids Fig. 5 Percentages of anthropogenic modifications within the ZooMS (orange) and morphology (blue) components on the dominant taxa at the sites of Bacho Kiro Cave, Les Cottés, and La Ferrassie. Numbers on the bars are the total NISP of specimens identified for the taxa (Canis lupus, Vulpes vulpes), felids (Panthera leo spelaea, Panthera pardus), cave hyaenas (Crocuta crocuta spelaea), and ursids within layers I and J at Bacho Kiro Cave while layer K exhibits only two carnivores remains with human modifications (SI Table 18). At Les Cottés, only two canid specimens show human modifications, and no human modifications were observed on carnivore remains within layer 6 at La Ferrassie (SI Table 18).

Skeletal representation
Due to their morphological specificities and as they are affected differently by taphonomic processes, teeth are largely represented in the morphology component and show the highest proportions among skeletal elements, particularly illustrated by the material from Bacho Kiro Cave and Les Cottés (Fig. 6 & SI Fig. 16). At Bacho Kiro Cave, the skeletal composition of the ZooMS component is mostly represented by long bones (LBN) and cranial and axial remains, with a higher proportion of axial elements within the ZooMS component explained by an overrepresentation of ribs (SI Table 19, Fig. 6). Rib elements are difficult to taxonomically identify as they do not retain many specific morphological features relative to their size and proportion in a skeleton. Long bone fragments (LBN) correspond to unidentified bone fragments from forelimbs, hindlimbs, and distal limbs (metacarpals and metatarsals). Bone specimens categorised as LBN within the ZooMS component are predominantly represented by diaphysis fragments (either from the midshaft or near the epiphysis of the bones) but rarely from the epiphysis, as illustrated by the example on the material of Bacho Kiro Cave (SI Fig. 17). Within the morphology component at Les Cottés, we observe relatively similar proportions of limb remains between the taxa, with the exception of the absence of hindlimb and distal limb remains recorded for Bos/Bison in US06 of both components, but higher proportions of cranial specimens from Bos/Bison and Equidae (SI Fig. 16). At La Ferrassie, the elemental representation of the ZooMS component only contributes to a small extent to the skeletal representation of the morphology component as most of the remains were unidentifiable and had not been assigned to a body part (SI Fig. 18).

Discussion
This study combined palaeoproteomic and zooarchaeological analysis of faunal material from three datasets covering the Middle to Upper Palaeolithic transition. It aims to overcome methodological limits in taxonomic identification Fig. 6 Skeletal distribution of the bone specimens identified through morphology (top) and ZooMS (bottom) from the dominant taxa at Bacho Kiro Cave. Numbers on the bars give the total NISP for each body part, layers, and ID method. Unidentified body parts (NID) were excluded from the plot. LBN, long bone fragment resulting from bone fragmentation and to address human subsistence and fauna processing behaviour during a period of possible interaction between Neanderthals and Late Pleistocene Homo sapiens groups in Europe. Together with a high success rate of taxonomic identification, the inclusion of ZooMS analysis of the fragmented, unidentifiable component of bone assemblages can identify species previously unrecognised through traditional morphological analysis and, furthermore, be integrated and correlated with traditional zooarchaeological, taphonomic, and ecological data (Berto et al. 2021;Sinet-Mathiot et al. 2019;Welker et al. 2015). In the case of highly fragmented bone assemblages, this addition can provide highly valuable information for the interpretation of human subsistence. This is notably exemplified in our study at La Ferrassie layer 6 with a fourfold increase in taxonomic identification through ZooMS compared to the morphologically identified component (NISP Morph = 142, NISP ZooMS = 518).

Prey selection and sampling bias
In the absence of alternative methods to address the fragmented component of Palaeolithic bone assemblages, previous studies of past human behaviour related to subsistence strategies have relied solely on morphologically identifiable fauna, excluding a vast majority of the available bone specimens. However, the fragmented component of Palaeolithic bone assemblages can differ significantly from the morphologically identifiable component, highlighted by differences in proportions of the dominant taxonomic groups between morphologically identified and ZooMS components. Our study does not reflect the pattern observed in several other ZooMS screening studies which found a similar taxonomic composition of dominant species between both components (Berto et al. 2021;Buckley et al. 2017;Welker et al. 2016Welker et al. , 2017. In this study, discrepancies in taxonomic abundances between both components are seen through an overrepresentation of reindeer and an underrepresentation of Bos/Bison and equids at the sites of Les Cottés and La Ferrassie. These differences seem to be related to differential identification rates between taxa, possibly creating a reporting bias in the representation of the dominant taxa depending on their ease of identification. Thus, taxa such as reindeer or Ursidae will be overrepresented in the morphologically identified component as they are easy to differentiate even when fragmented. On the contrary, Bos/Bison and Equidae are more difficult to distinguish when fragmented and are often categorised as unidentifiable remains.

Assemblage composition and identification rates
The uniformity in low weathering patterns on the bone material from Bacho Kiro Cave and La Ferrassie sites indicate that, throughout the stratigraphy, natural factors played a limited role in bone fragmentation. Overall, bone material was relatively quickly buried and suffered minimal re-exposure at these two sites. At Les Cottés, the degree of weathering was comparable among the dominant taxa, although reindeer showed slightly better bone surface readability. Further study is required to understand if this pattern could be explained through bone morphology or specific depositional conditions of the reindeer specimens (shorter exposition of the specimens prior to burying), especially knowing that glutamine deamidation ratios do not indicate a clear differential molecular preservation. Further, our analysis of collagen deamidation at each site does not provide a molecular diagenetic explanation for the differences in taxonomic proportion between the two bone components of each assemblage. When incorporating ZooMS identifications into the zooarchaeological analysis, we should keep in mind that both components, by definition, commonly show different bone length distributions, as larger fragments tend to be more easily identifiable morphologically. However, when comparing taxa, we note that this is not the case for all taxonomic groups (Pickering et al. 2006). Certain taxa, such as reindeer at Les Cottés and Capra sp. at Bacho Kiro Cave, can show a bone length distribution significantly different from other taxonomic groups (Bos/Bison and Equidae), possibly resulting from the size of the bone fragments most likely produced during marrow extraction and a different identification rate between these taxa. Indeed, because of the low cortical thickness relative to bone diameter and their smaller body size compared to Bos/Bison and equids, reindeer fragments will cover more of the reindeer bone proportionally, which would give it a better chance of preserving identifiable features. On the other hand, breaking open the bones of larger animals such as Bos/Bison or equids will produce larger fragments on average. Fragments of bovine bone specimens are often difficult to distinguish from homologous parts of Equidae or red deer as the skeletal elements of these taxa tend to overlap in size and morphology (Morin 2012). However, since reindeer are more easily identifiable, this results in increasing representation of this species within the morphological component alongside a limited proportion of identified Bos/Bison and equid specimens (Gobalet 2001).
The assessment of prey skeletal part distribution is often closely related to the taxonomic identification of the bone specimens. Small long bone shaft fragments tend to be difficult to identify due to a lack of diagnostic features on the bone diaphyses in combination with their high fragmentation rate due to marrow extraction (Morin et al. 2017a). Thus, the morphologically unidentifiable component analysed through ZooMS often contains a high proportion of long bone, particularly diaphysis portions, and rib fragments challenging the evaluation of skeletal distributions. Epiphysis portions tend to retain more specific morphological criteria facilitating the taxonomic identification of the remains. However, their representation within the long bone fraction of the ZooMS component do not strongly differ from the morphology component at Bacho Kiro Cave. Thus, an underrepresentation of epiphyses can also result from selective destruction due to various factors such as differential preservation and bone density, carnivore activity, specific butchering practices like extraction of bone grease, and post-depositional or sampling bias during the archaeological excavation (Binford 1981;Grayson and Delpech 2008;Morin 2010Morin , 2020Yravedra and Domínguez-Rodrigo 2009). Behavioural inferences such as carcass processing and the selected transport of different body parts are often made based on skeletal part representation and abundance (Bartram et al. 1999;Binford 1981;Klein et al. 1999;Marean and Assefa 1999). The integration of skeletal representation with the taxonomic identification obtained through ZooMS has the potential to add elements to the inventory of the faunal record, contributing to our understanding of the transport of articulated remains to the site.
However, further integrating ZooMS data with standard zooarchaeological investigations will require better comparability of the metrics used in both components. The number of identified specimens (NISP) (Grayson 1984) is commonly used as a proxy for species abundance among the ZooMS components, since the minimum number of skeletal elements (MNE) and minimum number of individuals (MNI) cannot be compared quantitatively to ZooMS data. If bone fragments are identified based on their taxonomy, regardless of morphology or surface preservation, experimental models may be built, which would help us better understand how bone fragmentation in Palaeolithic faunal assemblages varies according to specimen size, as well as identify skeletal elements. ZooMS can help refine alternative methods to calculate NME such as diagnostic landmarks on bone elements or refitting bone shafts (Marean et al. 2001;Morin et al. 2016;Stiner 1994).

Subsistence during the Middle to Upper Palaeolithic transition
The addition, through peptide mass fingerprinting, of taxonomically identified bone specimens to faunal assemblages spanning a transitional phase during human evolution contributes to our understanding of patterns of shifts observed during the MUPT. Our results contribute further detail to the general picture that, over this period, the hominin diet was dominated by a range of medium and large herbivores (Discamps et al. 2011;Gaudzinski-Windheuser and Niven 2009;Gaudzinski-Windheuser and Roebroeks 2011;Jaouen et al. 2019;Niven et al. 2012;Rendu et al. 2019;Richards et al. 2008;Smith 2015). Our work highlights the exploitation of a more diverse range of taxa by both hominins and carnivores, permitting the correlation of certain taxa with particular agents that were contributing to the bone accumulation on site, notably at Bacho Kiro Cave. Across dominant taxa, human modifications mainly consist of cut marks, with a low occurrence of percussion traces from marrow extraction, thus providing no suitable explanation for the difference in proportions between the components. The ZooMS analysis emphasises and refines shifts of proportions of taxa throughout the stratigraphy at Les Cottés, particularly between equids and Bos/Bison specimens . These shifts in the faunal composition could represent either a slow change in the prey availability in the environment around the site or human selection strategies paralleling the expansion of Late Pleistocene Homo sapiens over Europe. Nonetheless, while the morphologically identified fauna suggests a more specialised focus on hunting reindeer , our results suggest that this underestimates the exploitation of other species, in particular, Equidae. These results are particularly of interest within the framework of the debate about reindeer hunting specialisation (Grayson and Delpech 2002;Mellars 2004). Although the progressive increase of reindeer remains through the MUPT transition correlates with a progressive climatic degradation during MIS3 and can be explained by an adaptation of the human groups to environmental fluctuations (Banks et al. 2013;Discamps et al. 2011), the role of large ungulates in the human diet throughout the MUPT might have been underrepresented due to differential identification rates.
The incidence of carnivore modifications during late Neanderthal occupation of the sites suggests a context where both humans and carnivores were important in faunal accumulation and modification, still indicating frequent human occupation of the cave and sporadic carnivore visits, but the latter possibly more frequent than previously considered (Straus 1982). The progressive decrease of carnivore activity highlighted by the reduction of carnivore modifications from the MP to the UP at Bacho Kiro Cave and Les Cottés fits with the pattern previously detected in some other sites in Europe from this time period (Discamps 2014;Discamps et al. 2019;Rendu et al. 2019;Smith et al. 2021;Stiner and Kuhn 2006). This possible change of relationship of carnivores to humans from competitor to prey or source of raw material is emphasised by the appearance of human modifications on Ursidae remains during the IUP at Bacho Kiro Cave (alongside modifications on other carnivore species at Bacho Kiro Cave ). Homo sapiens started to exploit carnivore remains more intensively as a raw material, notably illustrated by the increase in bone artefacts made from cave bear bones and teeth at Bacho Kiro Cave and other sites in southeast Europe and southwest Asia during the IUP (Bosch et al. n.d.;Guadelli et al. 2011;Kuhn et al. 2009;Martisius et al. 2022;Stiner et al. 2013). Such specific needs in raw material can be investigated through skeletal part representation and carcass processing ). Furthermore, the higher percentage of carnivore traces in the Middle Palaeolithic layers at Bacho Kiro Cave and at Les Cottés attest to their repetitive use of the site correlated with possible short duration of human occupation Smith et al. 2021). The interaction between human groups and large carnivores seems to change during the MUPT and might indicate an increasing predatory pressure of human groups on their environment (Stiner and Kuhn 2006) and/or a shorter duration of site occupation by Neanderthals compared to Late Pleistocene Homo sapiens.
ZooMS screening of fragmentary components of Palaeolithic bone assemblages should be systematically undertaken alongside the taphonomic analysis of the taxonomically unidentifiable specimens (see, for example, Discamps (2021)). In addition, the integration of the faunal data obtained from aDNA retrieved from the sediment of an archaeological site with the zooarchaeological and ZooMS analysis of Palaeolithic faunal assemblages has the potential to provide a better understanding of the various episodes of occupation of a site or inform about the potential origin of the DNA preserved in the sediment.

Conclusion
The analysis of the morphologically unidentifiable component of Pleistocene bone assemblages offers an exciting new avenue for research. Our work on faunal assemblages from sites with occupational sequence that span the MUPT has highlighted inter-and intra-site differences between assemblages, taxa, layers, and identification methods. We emphasise that the morphologically unidentifiable component of faunal assemblages does not necessarily reflect the morphologically identified component. Certain taxa are more readily identifiable based on morphology compared to others. Their bone elements show particular features allowing for their recognition even when fragmented (Morin et al. 2017a). This results in a discrepancy in the identification rate of differing taxa during the analysis of bone material. Taxonomic abundances are influenced by these methodological limits, and any interpretation related to past human subsistence behaviour and hunting strategies can potentially be biased. Similar patterns might be expected in other monospecific faunal assemblages, and the assessment of morphologically unidentifiable bone fractions through ZooMS can reveal conditions that influence the variability of the results.
The integration of fragmentary bone components, identified through ZooMS or other biomolecular methods (Rüther et al. 2022), within a coherent zooarchaeological framework allows for a more exhaustive evaluation of the preserved bone assemblage, unlocking behavioural information based on skeletal part profiles, bone surface modifications, and ecological indices. Our large-scale, non-targeted ZooMS studies across the MUPT at Bacho Kiro Cave, Les Cottés, and La Ferrassie indicate an underestimated representation of the large ungulates such as Bos/Bison and Equidae, a progressive shift in prey selection from Bos/Bison to equids, a reduction in the frequency of site occupation by carnivores, and an increase in their exploitation by Upper Palaeolithic Homo sapiens over the course of their progressive dispersal across Europe. This approach provides complementary data for assessing preserved bone remains, contributes to our understanding of bone assemblage formation, and represents a future path for Palaeolithic zooarchaeology.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.