The “Bear” Essentials: Actualistic Research on Ursus arctos arctos in the Spanish Pyrenees and Its Implications for Paleontology and Archaeology

Neotaphonomic studies of large carnivores are used to create models in order to explain the formation of terrestrial vertebrate fossil faunas. The research reported here adds to the growing body of knowledge on the taphonomic consequences of large carnivore behavior in temperate habitats and has important implications for paleontology and archaeology. Using photo- and videotrap data, we were able to describe the consumption of 17 ungulate carcasses by wild brown bears (Ursus arctos arctos) ranging the Spanish Pyrenees. Further, we analyzed the taphonomic impact of these feeding bouts on the bones recovered from those carcasses. The general sequence of consumption that we charted starts with separation of a carcass’s trunk; viscera are generally eaten first, followed by musculature of the humerus and femur. Long limb bones are not broken open for marrow extraction. Bears did not transport carcasses or carcass parts from points of feeding and did not disperse bones appreciably (if at all) from their anatomical positions. The general pattern of damage that resulted from bear feeding includes fracturing, peeling, crenulation, tooth pitting and scoring of axial and girdle elements and furrowing of the upper long limb bones. As predicted from observational data, the taphonomic consequences of bear feeding resemble those of other non-durophagus carnivores, such as felids, and are distinct from those of durophagus carnivores, such as hyenids. Our results have paleontological and archaeological relevance. Specifically, they may prove useful in building analogical models for interpreting the formation of fossil faunas for which bears are suspected bone accumulators and/or modifiers. More generally, our comparative statistical analyses draw precise quantitative distinctions between bone damage patterns imparted respectively by durophagus (modelled here primarily by spotted hyenas [Crocuta crocuta] and wolves [Canis lupus]) and non-durophagus (modelled here by brown bears and lions [Panthera leo]) carnivorans.


Introduction
A major focus of neotaphonomic research is to understand the cause-and-effect relationships that are involved when modern carnivores act as bone-modifying agents. From paleontological and archaeological perspectives, a major objective of this kind of research is to use what we learn about the taphonomic impact of carnivores in order to, in stepwise fashion, isolate evidence of noncarnivore, presumably hominin-induced effects in fossil fauna palimpsests. It is from that point that hypotheses of carcass foraging by ancient human ancestors can proceed.
Fossils of the genus Ursus (including the extinct cave bear species Ursus deningeri and Ursus spelaeus) are very common in assemblages from European Pleistocene sites-including those that also contain hominin remains and hominin-produced stone artifacts. These regular paleontological associations of cave bear fossils and especially Neanderthal (Homo neanderthalensis) fossil and archaeological remains formed the basis of dramatic and influential hypotheses that Neandertals both hunted and worshipped cave bears (e.g., [22][23][24][25][26][27]).
In contrast to their prominent role in scenarios of ancient European life, from a neotaphonomic perspective, ursids are among the least studied of all common, large carnivores. Given that most experts recognize that bears were very likely significant taphonomic agents in past (especially in karst contexts, where their hibernation-related activities and probable cannibalistic inclinations held great potential to modify the spatial distribution and condition of the bones of their conspecifics and other animals (e.g., [28][29][30][31][32][33][34][35][36][37][38][39][40]), this imbalance is doubly unfortunate.
We acknowledge that Ursus arctos arctos is obviously not an exact match for its prehistoric congeners (e.g., [41][42][43]). For instance, unlike Pleistocene bears, extant Pyrenean brown bears face little to no feeding competition for large carcasses (see below) and, depending on their exact distribution on the landscape, have access to different arrays and quantities of other foods. These differences probably impact the relative intensity of bear feeding behaviors, as well as condition their decision making processes about whether to transport (or not transport) carcass parts to sheltered feeding sites. We do not believe, however, that these differences disqualify modern brown bears as suitable taphonomic referents for extinct cave bears. Instead, we highlight relevant continuities between the two types of bear: modern brown bears are about the same size as were cave bears [41]; the skulls of both taxa have ''domed'' frontals and prominent, similarly placed sites for the attachments of large temporalis and masseter muscles [44][45][46], imbuing both with powerful bite forces [43,45]; both-based on behavioral observations for modern brown bears and on isotopic [37,47,48] and occlusal microwear [49] evidence for cave bears-are/were omnivorous.
Couturier's [50] work in the Cantabrian Mountains (northern Spain) was the first study to describe the sequence by which modern brown bears consumed large vertebrate carcasses, but his observations were also unsystematic and lacked a contemporary taphonomic focus. It was almost thirty years later that Haynes [51][52][53] became the first researcher to describe modifications of selected skeletal elements of cattle (Bos taurus) fed to captive bears; Haynes then used these results to interpret taphonomic damage on carcass remains of wild North American ungulates. Haynes's observations and conclusions-that bears usually abandon feeding on a carcass once it has been defleshed and eviscerated, and that this behavior produces only minimal taphonomic damage, including minor tooth pitting, scoring, furrowing, crushing and perforating of cortical bone (see below)-have often been recruited in analogical models that assert ursid involvement in the accumulation and/or modification of some Pleistocene fossil faunas (e.g., [54,32,33,36,37]).
More recently, Saladié et al. [55] fed fleshed, but disarticulated, limb segments of young cattle, pigs (Sus scrofa) and sheep (Ovis aries) to captive brown bears in the Barcelona Zoo and in the Hosquillo Natural Park (Spain). In contrast to the general patterns of bear-induced taphonomic damage documented by Haynes [51][52][53], Saladié et al. [55] found that bears generated bone surface destruction and modification patterns very comparable to those produced by other, better-studied large carnivores. Such patterns include: a high degree of bone breakage; the abundant production of ''diaphyseal cylinders'' (i.e., long limb bone [LLB] specimens that lack epiphyses, which were chewed-off, but that retain intact, full-circumferenced shafts); the furrowing, and scooping-out of trabeculae from LLB epiphyses. Saladié et al. [55] also showed that the metric range of individual bear-produced tooth marks overlaps with that previously documented for African lions (Panthera leo).
Confusing these apparently clear taphonomic patterns were the conclusions of Sala and Arsuaga [56], who tracked the modification of nine domestic equid carcasses consumed by free-ranging brown bears in the Cantabrian Mountains. In addition to the interesting behavioral data that Sala and Arsuaga [56] phototrap study yielded (e.g., the Cantabrian bears often moved carcasses before feeding on them, but never into caves or other secluded refuges; observations that are obviously relevant to models that prehistoric bears collected carcasses and bony residues in shelter sites), their taphonomic results, in contrast to those of Saladié et al. [55], included much reduced intensities of tooth-marking and fracturing of bones and clustering of most tooth marks on axial, rather than on appendicular, skeletal elements. Sensibly, Sala and Arsuaga [56] attribute the stark difference in the taphonomic patterns that emerged from their study and those documented by Saladié et al. [55] to the fact that the former was conducted with free-ranging subjects and whole carcasses and the latter with captive animals and selected carcass portions of much smaller ungulates.
In summary, each study in the development of research on brown bears as taphonomic agents has been an incremental improvement to its predecessor(s); each has also proven of great value in the interpretation of fossil faunas suspected to have been (at least partially) formed and/or modified by prehistoric bears. By presenting new data and conclusions on the taphonomic behaviors and consequences of free-ranging brown bears in the Pyrenees of Lleida (Spain), here we contribute to the evolving research program on this critically important but understudied large carnivore. Importantly, our study advances the state of knowledge about brown bears as taphonomic agents because our study group is composed of free-ranging bears who were granted unfettered access to complete, fresh carcasses of seventeen ungulates.
In addition, we contextualize our results by presenting a crossspecies, multivariant analysis, comparing our taphonomic data to those generated in other studies on other major large carnivore taphonomic agents, including spotted hyenas (Crocuta crocuta), lions and wolves (Canis lupus). (Taphonomic data on brown hyenas [Parahyaena brunnea] and striped hyenas [Hyaena hyaena] are utilized more restrictively in selected analyses.) In doing so, this analysis achieves a broader goal of further refining those taphonomic signatures that are respectively diagnostic of durophagus (modelled here by spotted hyenas and wolves) and non-durophagus (modelled here by brown bears and lions) carnivorans. These results will prove useful to paleontologists and archaeologists investigating questions about Pleistocene large mammal terrestrial ecosystems.

Ethics Statement
All the research images analysed in this study were obtained by the monitoring teams from the Conselh Generau d'Aran and DAAM (Departament d'Agricultura, Ramaderia, Pesca, Alimentació i Medi Natural) of the Generalitat de Catalunya and funded by the Ministerio de Agricultura, Alimentación y Medio Ambiente from Spanish Government. All the carcasses were recovered under the supervision of these teams, who are the responsible of the protection of wildlife in the studied areas (Pallars Sobirà and Valh d'Aran, Spain), and managed the administrative permissions required, including the owner's permission in the case of domestic animals. No animals were damaged or sacrificed by the researchers and/or rangers during the present project.

Study area and subject population
The Pyrenees mountain system extends westward from the Mediterranean Sea to the Atlantic Ocean, separating the Iberian Peninsula from the rest of Europe ( Figure 1). The formation of the range, related to the Alpine orogeny, resulted in metamorphic peaks (some of which attain altitudes of 3,000 m above sea level), with deep valleys in between. The region is characterized by an Alpine climate and concordantly graded vegetation, dependent on relative altitude. The Pyrenees are sparsely populated by humans, who subsist there largely on a mixed economy, based mostly on animal husbandry and tourism.
Brown bears were naturally distributed in the Pyrenees until the recent historical period, but were all but extirpated by the 1950s. Transporting wild brown bears from Eastern Europe (mostly Slovenia), a reintroduction program, under the banner of Project Life (funded by the European Union and the French, Catalonian and Aragonian governments), began in 1996. Since that time, the brown bear population has been continuously monitored (originally with external satellite-connected transmitters and other intraperitoneal devices and now by use of photo-and videotraps, direct observation, hair traps and fecal collection), as has the transformations of their range. Today, the Pyrenean brown bears comprise a controlled population of 25-30 individuals, with 12-14 adults of a balanced sex ratio.
In order to assure the bears' recurrence there, current photoand videotrap locales are baited with fruits, ungulate carcasses and glandular odors. It is at these locales that research/monitoring teams are able to regularly observe the bears' health and behavior. Based on longitudinal data generated at observation points, adult males have annual ranges of 700-2,000 km 2 , while females have annual ranges of ,300 km 2 . Individuals of both sexes range between 1,600-2,300 m above sea level, across the intramountain regions of Pallars Jussà and Pallars Sobirà (in the Mediterranean watershed) and Val d'Aran (in the Atlantic watershed). Bear diet in these regions is primarily herbivorous but is supplemental by some animal product gained through predation and scavenging of ungulate ( ) and smaller vertebrate carcasses, as well as some insectivory. Because of the absence of other large carnivores in the monitored regions, bears face a distinct lack of interspecific competition for edible resources from vertebrate carcasses. This factor probably accounts, at least in part, for the observation that bears generally feed recurrently on a single large carcass for anywhere from only one to just less than eight weeks.

Bear-derived sample
We analyzed the remains of seventeen variously sized ungulate carcasses fed on by bears in Pallars Sobirà and Val d'Aran during spring and summer months between 2010 and 2013 ( Figure 1). For this sample, Table 1 summarizes the number and identity (age, sex, weight) of bear consumer(s) per carcass, carcass taxon, estimated ontogenetic age of carcass and state of carcass completeness when collected by our research team. These carcasses were made available to the bears by: (1) bear predation; (2) bear scavenging of natural deaths; (3) road kills transported by forest rangers to observation points along routes used by bears. We were able to document bear consumption of carcasses in the case of predation on domestic animals because livestock owners notified bear monitoring teams of those kills. In cases of bear scavenging, monitoring teams installed photo-and videotraps in the vicinity of carcasses. Thus, human intervention with carcasses was minimal until cessation of bear feeding, which was signaled by prolonged lack of visits by bears to the carcass. Exact location and spatial dispersion of the remains of each carcass were documented using a GPS and camera.

Taphonomic analyses
Based on body masses of living adult animals, the sample of carcass remains was divided into three categories: large (.300 kg; includes cattle and horses), medium (300-100 kg; includes red deer), and small (,100 kg, includes roe deer, sheep and goats). Similarly, the bear consumers were divided into large (.100 kg) and small (,100 kg) individuals.
We used strong incandescent light and 10x power hand lens, as well as an Olympus SZ 11 stereoscopic (magnification up to 110x), together with an ESEM (FEI QUANTA 600 Environmental Scanning Electron Microscope), to identify and classify bone surface damage on the collected ungulate bones, following established standards (e.g., [51,[57][58][59][60]). The damage we quantified included: (1) Tooth pits and punctures are caused by the penetration of the tooth (or teeth) of taphonomic agent into (pitting) or through (puncturing) a bone's cortex. Both types of damage results from static loading of a bone's surface between the agent's upper and lower teeth [57,58]. (2) Tooth scores are ''essentially elongated tooth pits that are usually U-shaped in cross-section, where a consumer's tooth cusp maintained brief contact with […] the bone surface as the tooth cusp dragged away from the pit basin'' [61, p.124]. (3) Furrowing is the chewing-driven deletion (e.g., ''scooping out'') of trabeculae, usually from an LLB's epiphysis [51,53,58,62]. (4) Crushing is produced by the collapse of localized areas of cortical tissue under the static loading of a consumer's jaws and teeth, forming cracks and splitting in these regions [58]. (5) Crenulated edges often occur on flat bones in the form of ragged, notched terminations imparted by a consumer in areas where its bite force overcame the structural strength and density of bone tissue [57,58]. (6) Classic and general peeling of bone [63] are regions of bone where cortical layers have been stripped away [64], usually by the action of a consumer's anterior teeth; with incipient peeling, cortical layers are simply frayed and perhaps splayed or bent but not removed from the bone [63].
For quantitative analyses of bone damage distributions and frequencies, we divided carcasses into anatomical regions, including: the skull (cranium and mandible); axial skeleton (ribs and vertebrae); girdles (scapulae and pelves); LLBs (including humeri, radioulnae, femora and tibiae, but excluding metapodials, which are essentially meatless and thus rarely intensively modified by non-durophagus meat-eaters). Teeth and podials were excluded from most of our quantitative analyses. For overall anatomical region distribution and frequency patterns, we display our statistical analytical results on ternary graphs [65], which are based on a likelihood method that is, in turn, an improvement over traditional bootstrapping methods. With ternary graphing, the likelihood of confidence intervals does not depend on sample size [65]. An additional advantage of the likelihood method is that it can accommodate counts of zero, displayed on ternary graphs at vertices. For length and breadth (mm) of bear tooth pits, punctures and scores, we performed a Mann-Whitney test between all possible pairs of samples to determine which differed significantly, including variables such as tissue bone (cancellous and thin cortical bones), size of bear consumers and size of ungulate bones.

Analysis of LLB modification patterns
Our more specific analytical treatment of the taphonomic alteration of LLB specimens was explicitly comparative. It is obvious that the more alike the processes and fundamental variables (e.g., environment/habitat, community structure, season,  Table 3. Frequencies of long limb bone (LLB) portions (P = proximal; D = distal) with tooth mark damage, tooth pit/tooth score ratios and frequencies of LLB cylinders and complete LLBs in comparative faunal samples (data sources listed in Table 2; data for brown bears are from this study). etc.) of actualistically derived taphonomic models and the paleontological/archaeological subject materials they attempt to model, the more robust is the heuristic utility of the former (e.g., [14,66,67]). Accordingly, we strove to control, to the greatest extent possible, as many variables as possible across our three comparative samples. First, the composition of our comparative samples is uniform, consisting only of LLB specimens. Second, all samples formed under wild, free-ranging situations, through the taphonomic actions of large, gregarious carnivores, including spotted hyenas, lions and wolves (samples formed by brown and striped hyenas are also included in our analyses, but in more minor ways). Third, we were able to largely control for carcass taxon and size by utilizing observations of samples that are composed in large parts by the bones of equids modified by our comparative carnivore subjects. Table 2 summarizes key aspects of the comparative samples we utilized in this study [4,10,11,15,[67][68][69][70][71][72].
For the comparative, interspecific analyses of the LLB subsamples, we focused on the following variables: (1) Operating from the hypothetical assumption that, all things being equal, durophagus carnivores are likely to produce much greater frequencies of furrowing on LLB epiphyses, we scored the simple presence or absence of furrowing on the ends of LLBs. More specifically, we conducted multivariate statistical analyses of the occurrence of furrowing on distal humeri, proximal ulnae, proximal and distal femora, and proximal tibiae.
(2) We quantified tooth mark (including pits, punctures and scores) frequencies per bone specimen. A range of measures for pits and punctures has been established to allow comparison with those provided by other researchers [55,56,73] and on other carnivores (e.g., [36,[72][73][74][75][76][77]). (3) All things being equal, durophagus carnivores are predicted to destroy bone epiphyses completely in greater frequency than do non-durophagus carnivores (e.g., [78]). Therefore, for each sample, we quantified LLB completeness thusly: complete; cylinder (a LLB specimen lacking its epiphyses, represented by a portion of diaphysis that for, at least some of its length, maintains its complete, original circumference); LLB diaphysis fragment (a shaft fragment that maintains ,100% of its complete, original circumference). (4) All things being equal, durophagus carnivores are predicted to impart a higher proportion of tooth pits to tooth scores on dense cortical bone than are non-durophagus, who primarily focus only on meat-stripping rather than on meat-stripping and bone-cracking, as do durophagus carnivores (e.g., [6]). Thus, we calculated pit: score ratios for all comparative samples.

Statistical analysis of LLB modification patterns
In order to document relationships among analytical components of the comparative LLB samples, we conducted principle component analyses (PCA), which produce factors that result from the reduction of dimensionality caused by multiple variables. With exploratory PCA, the analyst aims to improve prediction and variance accountability by detecting those variables that do not contribute significantly to sample variance; with confirmatory PCA, the analyst uses selected variables to maximize sample variability and sample component relationships. The use of continuous numerical variables may result in bias of PCA solutions due to the heterogeneity of these values and overemphasis of the weight of variables displaying high numerical values. For this reason, variables are usually centered and scaled prior to their statistical analysis. However, in the present analysis, all variables involved the use of percentage values; they have similar scale, so there was no need to center and scale variables.
As opposed to dimension reduction by orthogonal projection as performed in PCA, in multidimensional scaling (MDS) points are chosen so that stress (the sum of the squared differences between the inter-sample disparities and the inter-point distances) is minimized [79]. The MDS option we selected for our analyses is the identity transformation, which consists of taking the intersample disparities as the inter-sample dissimilarities themselves. This metric MDS approach uses a Pythagorean metric analysis of inter-point distances, which includes an iterative majorization algorithm to find the MDS solution [79]. This algorithm was considered to have converged as soon as the relative decrease in stress was less than 10 26 . The algorithm was also stopped once greater than 5,000 iterations were performed. In MDS, points are   Table 4. Brown bear induced-damage by ungulate body size (see text for definitions) and skeletal element. related in a low-dimensional Euclidean space [79], with data spatially projected by regression methods admitting non-linearity. The use of MDS in the present study therefore complements the PCA test. We also employed a canonical variate analysis (CVA). CVA focuses on data grouped into K classes by transforming original variables into canonical variables defined by square distances between the means of the groups obtained by Mahalanobiss D 2 . This is scale invariant. CVA produces a higher degree of separation between the group means than does either PCA or MDS [80]. The biplot axes for CVA were determined by the group means.
In a separate database, we entered LLB completeness data (see above, point 3) and tooth score/tooth pit ratio (see above, point 4) data. Both data sets were analyzed together using a Redundancy Analysis (RDA), which combines a multiple regression with a PCA. Data were previously normalized using a Hellinger transformation method [81]. Missing data were converted using group average values.
Because results of the multivariate statistical analyses summarized above differentiated our comparative carnivore groups (see below Results), our next step was to create thresholds from which frequencies of different modification types could be used as lines to demarcate the different types of carnivoran taphonomic agents. To this end, we submitted data for all taxa and all variables to a tree-based analysis.
Regression and classification trees usually suffer from high variance, but a procedure with a repeated sequence of data sets derived from the same original sample will decrease that variance. This is the general purpose of bootstrap aggregation, more commonly known as bagging, which splits the original training data set (TDS) into multiple data sets derived from bootstrapping the original TDS. Similarly, the powerful procedure of random forests (RF) performs a bootstrapping approach, but samples only a subset of variables for each tree. Thus, RF produces a final solution that includes a selection of variables that are important for correct classification of the analytical set. RF produces hundreds of trees that are repeatedly fitted to bootstrapped sets of data. The results are contrasted against a validation test, from the observations (about one third) not used for the training data set (these observations are referred to as out-of-bag (OOB) observations). RF produce estimates on how many iterations are needed to minimize the OOB error. The importance of each response variable is determined by mean decreased error (MDE) for regression trees (RT), whereas the Gini index is more useful for classification trees (CT). The following step is the creation of a classification tree using the most useful variables selected by the RF test.
We conducted an RF test by selecting three-variable sets in each aggregated bootstrap iteration [10,15,69,70]. A total of 500 trees was generated. The test was applied to all variables in Table 3. Subsequently, it was applied only to the variables dealing with LLB portion modification. Variable selection for CT analysis was carried out when variables showed a mean decreased Gini index of .0.5.
The variables with higher MDE values were then selected to obtain single regression trees. CT tests were carried out by selecting a complexity parameter (cp) of 0.005.
doi:10.1371/journal.pone.0102457.t006 Table 1 summarizes many relevant details of the 17 ungulate carcass occurrences modified by bears. In general, feeding bears rarely moved carcasses any appreciable distances, if at all. Most carcasses maintained near natural anatomical position and connectivity of skeletal elements upon recovery by our team. In the rare instances of distarticulation and disassociation of parts (between 50 to 250 m 2 from a main carcass cluster), we note that much of this separation might be explained the actions of nonbear, secondary consumers, such as foxes (Vulpes vulpes), martens (Martes martes) and vultures (Gyps fulvus, Gypaetus barbatus)-all of which are active in the region (Observations [OB] 4-7). Minor water flows could have also carried remains away from carcass concentrations (OB 17). The very incidental movement of carcasses by bears most probably relate to the minor agitations caused by infrared lights on the photo-and videotraps and human scents in the vicinities of carcasses. But, in agreement with the observations of Sala and Arsuaga [56] in Cantabria, the Pyrenean brown bears did not transport carcasses or carcass parts to dens or other shelter sites. Although Cantabria is home to wolves, like us, Sala and Arsuaga [56] did not observe large, non-ursid carnivores influencing bear feeding behavior. In more intensive competitive situations bears might be more likely to transport and secure carcasses in shelter refuges. Thus, as it now stands, neotaphonomic observations of extant brown bears hold little relevance for modelling potential carcass transport by extinct ursids, which are sometimes implicated in the substantial accumulation of European Pleistocene cave faunas.

Patterns of carcass consumption by brown bears
In nine instances (OB 4,(8)(9)(10)(11)(12)(13)(14)17), bears skinned carcasses before consumption (see also, [50,82] (Figure 2). All cases of skinning involved small sized carcasses; no large carcasses were completely skinned. An example of partial skinning by a bear is shown in Video S1, in which a young bear consumes the viscera of an adult male deer without removing the carcass's whole skin. Skinning seems to have facilitated separation of a carcass's trunk and limbs, after which the trunk's viscera was removed and consumed. We do not associate any subsequent superficial modifications on recovered bones to skinning by bears. Across carcass sizes, removal and consumption of viscera usually commenced at the abdomen of the carcass, with bears employing their forepaws to wrench back ribs dorsally, away from their sternal articulations. In only one case (OB 7), involving the consumption of a large (.500 kg) cow carcass, did a bear initiate removal of viscera differently, by moving from the carcass's throat in a caudal direction toward the lower abdomen and viceversa. To deal with this large carcass, the bear consumer took a multistep tact to opening the cow's ribcage, as is viewable in Video S2. The bear initiated the process by compressing the cow's ribcage under its forepaws and body weight, fracturing ribs. Next, the bear pried apart the ribcage with its forepaws, and then bit the sternal ends of ribs, separating them with brisk jerking of its head and the assistance of its forepaws (some examples of these actions are shown in Figure 3). In parallel, some portions of musculature covering vertebrae and pelves, and to a lesser extent upper LLBs, are defleshed and consumed (OB 1,2,6,8,11,(13)(14)(15)(16).
Following the consumption of viscera and meat, bears focused on chewing-directed extraction of intrabone nutrients from especially the proximal epiphyses of humeri and tibiae, both epiphyses of femora, from the innominates and from the bodies of vertebrae. Sala and Arsuaga [56] observed similar behaviors among the Cantabrian bears they studied. However, unlike the observations of Sala and Arsuaga [56], as well as those of Haynes [53] and Pinto and Andrews [83], the Pyrenean bears did not fracture open LLB diaphyses for access to yellow marrow. Unique among bear studies, we did observe an instance of a bear breaking open a deer spine along the lumbar section in order to gain access to edible spinal cord tissues (Video S3).

Overview of the taphonomy of the brown bear-derived sample
The aggregated osteological sample from all 17 ungulate carcasses consumed by bears shows a variety of consumptionrelated damage, including tooth punctures, pits and scores, crushing, the furrowing of LLB epiphyses, crenulation, fracturing and peeling (Table 4). Of the total recovered number of identified specimens (NISP) of 1,173, 508 (43.3%) bear at least one occurrence of taphonomic modification, with most of the damaged specimens deriving from bones of the axial (NISP = 459) and girdle (n = 22) regions (Tables 5, 6; Figure 4). Only five skull specimens and 21 LLB specimens (as well a single patella) preserve bearinduced modifications. Ungulate carcass size seems to have had only minimal influence on divergences from these general patterns of assemblage-level modification, with slightly more intensive damage observed on scapulae and pelves derived from smaller carcasses (Table 4).
Here we summarize the general patterns of bone surface modification (BSM) frequencies and anatomical distributions:  Most fractures are observed on vertebrae (NISP = 124), including on especially lumbar vertebrae and specifically on lumbar transverse and spinous processes, and on ribs (NISP = 90) ( Figure 5G). Fractures resulted primarily from evisceration-related activities, when ribcages were detached from sterna and pried apart by bears.
As with fractures, we associate most peeling damage with evisceration and wrenching of ribcages by bears. One hundred thirty-four specimens (11.4% of the total aggregated NISP) are peeled. The most common type of peeling is classic peeling (NISP = 76), which is seen especially prominently on lumbar vertebrae (NISP = 31), with much damage on transverse processes and much more rarely on spinous processes ( Figure 5F). Fewer (NISP = 12) thoracic vertebrae show classic peeling, with most of it on spinous processes. Only two goat scapulae show this type of modification -on portions of the coracoid process (OB 12) ( Figure 6D). A NISP of 30 ribs are also classically peeled, with most of that damage clustered at rib sternal ends and midshafts ( Figure 5B). A total NISP of 58, including 50 ribs of small and medium sized ungulates and  Table S1 for descriptive statistics seven lumbar vertebrae of small ungulates are incipiently peeled ( Figure 5C). Our observations of peeling damage in the bear-modified sample are especially significant, given that this is type of taphonomic damage was heretofore primarily associated with the feeding activities of hominins (e.g., [63,64,84]) and chimpanzees (Pan troglodytes) (e.g., [85][86][87]). Anecdotal evidence suggested that other taphonomic agents, including spotted hyenas (Gary Haynes, personal communication, 2013, in [63]), also occasionally peel bones, but to our knowledge it is still relatively uncommon in non-primate-produced faunas (e.g., [63]). Thus, the high frequency of peeling produced by bears was unexpected, especially considering its apparent lack of occurrence in previously studied bear-created faunas, including that from Cantabria, described in great detail by Sala and Arsuaga [56]. We note that peeling is distributed primarily on bone specimens from small ungulates (NISP = 101) in our sample, and that, in contrast, Sala and Arsuaga [56]'s sample comprises large carcass remains; this distinction might be reason for the disparity in inter-assemblage observations of peeling damage. 2. Crushing damage is distributed broadly across skeletal parts, including on 25 cervical, 20 thoracic and 10 lumbar vertebrae of large, small and medium ungulates, on the caudal vertebra of a sheep (OB 10), and on three scapulae, seven innominates and four sacra of small and medium carcasses. Crushing is also observed on the nasal bones of a goat (OB 9), on the olecranon process of a red deer ulna (OB 4). Among the damaged vertebrae, crushing is identified mainly on the spinous (NISP = 31; 55.3%) and transverse (NISP = 16; 28.6%) processes. Occasionally, this modification is detected on mammillary apophyses (NISP = 5; 8.9%) and on the vertebral body (NISP = 4; 7.1%). Several ribs (NISP = 43) are crushed, including on both ends, and on necks and midshafts. Given its wide anatomical distribution, crushing is associated with both ribcage opening activities and general feeding/chewing.  3. Longitudinal cracking of specimens is relatively rare in the aggregated bear-derived sample (NISP = 11), but clusters in the axial region and is thus probably related to prying and wrenching of ungulate thoraxes when bears sought access to viscera ( Figure 5A). Seven rib specimens from small ungulates and one from a large ungulate are cracked. The other cracked specimens include two thoracic vertebrae (one of a small ungulate, one of a large ungulate) and a lumbar vertebra of a medium ungulate. 4. Fifty specimens (4.3% of the total NISP) are crenulated, resulting mainly from directed chewing on exposed bone surfaces. As expected, based on their (at least partial) flat morphology and thinly distributed cortical bone, most crenulation is preserved on axial and girdle elements, concentrated (by definition) on specimen terminations ( Figure 5D; Figure 6A-C; Figure 6H). Only two LLB specimens-a cow femur (OB 7) and a red deer humerus (OB 6) -have crenulated ends, which are, in these cases, essentially the diaphysis-ward extension of intensive furrowing that commenced on trabeculae-filled epiphyses. Proper furrowing of LLB (excluding meatless metapodials) epiphyses (resulting from bear chewing that was directed in order to gain access to intrabone grease and red marrow) occurs on 17 of 64 (26.6%) of the recovered non-metapodial LLB specimens ( Figure 7A-D). In contrast, only 2.9% (NISP = 19) of the total axial element NISP of 646 is furrowed (the basicranium of a goat  ] also show damage that is best characterized as furrowing). Together, humeri (NISP = 6) and femora (NISP = 6) account for 70.6% of the furrowed LLB sample, with all six humeri displaying furrowing on their proximal epiphyses only. For damaged femora, furrowing is distributed more evenly anatomically, with some femora furrowed only proximally, others furrowed only distally and still others on both epiphyses. Only one femur of goat (OB 12) shows a more intensive consumption on its proximal epiphysis, leading to an almost incipient scooping-out; however, the head is still preserved and this alteration therefore cannot be considered as a classic scooping-out. Four tibiae and one ulna are also furrowed.  Figure 5E). One hundred thirty-three axial specimens and 14 girdle specimens are punctured, pitted and/or scored ( Figure 6E-G). Very few skull specimens (NISP = 3) are punctured, pitted and/or scored. Likewise-and in conspicuous contrast to faunas modified by most other large carnivores (e.g., [6,14,15,51,58,59,61,67,72,78,88])-few of the recovered LLB specimens are punctured, pitted and/or scored. Moreover, none of the punctured/pitted/scored LLB specimens in the bearderived sample is damaged on its diaphysis; instead all puncturing/pitting/scoring is restricted to LLB epiphyses. The  (Figure 8; Table S1; Figure S1). This range is consistent with that documented in other bear studies by Saladié et al. [55] and Sala and Arsuaga [56]. Small bears created the largest marks in our sample, including those on the remains of a horse (OB 1) and those on red deer bones (OB 15). The only other occurrence of bones that comes near to preserving similarly large marks was created by a large bear that fed on small roe deer. On the face of it, it seems counterintuitive that it was small bears who created the largest tooth marks in our aggregated sample. We suspect that small, young bears, recently on their own, might face more intraspecific competition for carcass resources and thus, when able, extract carcass  nutrients more intensively and completely than larger, older bears.

Modification patterns on LLB: Comparative results
PCA yielded a two-dimensional solution accounting for 61.48% of the sample variance. The Eigen value for the first dimension is 45.13% (explaining most of the variance) and the second is 16.35%. Most variables display a similar value, with proximal radius and proximal tibia having the least impact on the final solution (Table 7). Considering both components, these portions, together with both femoral epiphyses, show a reduced impact  compared to other LLB portions. Both ends of the humerus (especially the distal epiphysis) display a higher discriminatory value for differentiating taphonomic agents (Figure 9). This may have to do with taxon-specific styles by which different carnivores consume ungulate forelimbs.
CVA yields a two-dimensional solution that explains 69.5% of the sample variance, with the first dimension explaining 53.6% of that variance. This first dimension was mostly concerned with intergroup differences based on comparisons of the distal humerus, distal radius and distal tibia, and, to a lesser degree, the distal femur. In distinction, the proximal humeral, radial, tibial, and femoral epiphyses were substantially less useful for differentiating groups in the first dimension. The proximal humerus was the most useful discriminator in the second dimension. Thus, with regard to modification of ungulate humeri (especially the distal epiphyses of humeri), the CVA supports the PCA in stressing distinct patterns created by brown bears, spotted hyenas, lions and wolves-the four major taphonomic agents investigated here. In addition, the ungulate distal radius clearly separates the taphonomic impacts of the comparative carnivores in Euclidean space. In sum, the proximal and (especially) distal epiphyses of ungulate humeri and the distal epiphyses of ungulate radii are modified most intensively by spotted hyenas, next most intensively by wolves, and significantly less intensively by lions and brown bears. These comparative conclusions are apparent in the mean relative absolute error for each bone portion, with the proximal (PCA) and distal (CVA) humerus and the distal radius displaying the lowest values.
Additionally, a two-dimensional, MDS-yielded solution shows that the proximal and distal epiphyes of humeri are the most influential first-dimensional variables ( Figure 10). The distal humerus also had the highest score for the second dimension. The distal radius and distal femur were also highly influential in the first dimension (Table 8). These LLB portions also show the smallest mean relative absolute error (Table 9).
Collectively, the results of all three statistical analyses converge to show that the four major taphonomic agents investigated here are readily separated based on the relative intensity that each modifies LLB epiphyses. As predicted by the general, observationally based understanding that spotted hyenas and wolves are (in contrast to brown bears and lions) durophagus carnivorans, these two subjects generated the most intense taphonomic damage observed on LLB specimens. On the intrasample scale, the wolfgenerated sample evinces the greatest variability of all four samples. We propose this high intrasample variability is, at least in part, a product of the aggregated nature of the wolf-generated sample, with some ungulate bones being those of young animals, modified over the course of several feeding/bone chewing episodes, while others derive from adult individuals [67]. However, we also note that the spotted hyena-generated sample is similarly variable to that generated by wolves; together, the considerable intrasample variability of the spotted hyena-and wolf-generated samples stands in contrast to the reduced intrasample variability of the brown bear-and lion-generated samples. This clear separation between the two groups-spotted hyenas and wolves on one hand, brown bears and lions on the other-probably reflects a more generalizable separation between durophagus and non-durophagus carnivores: the latter, acting primarily as carcass eviscerators and defleshers-and not as bonecrackers-avoid or minimize bringing their teeth in contact with bone surfaces; the former feed with less compunction-not only consuming innards and flesh but also readily demarrowing bones, and, in the course, inflicting a high degree of damage to LLBs. This broad contrast between durophagus and non-durophagus carnivorans is readily apparently in our PCA (and to a less extent in our CVA), which shows obvious separation between spotted hyenas and the other taphonomic agents that we investigated.
The RDA that combined the LLB epiphysis-modification dataset and the dataset composed of LLB completeness frequencies and tooth pit:tooth score ratios yielded a two-dimensional solution whose accumulated constrained eigenvalues (total explained variance) amounted to 96.9%. The first dimension explained 73.6% of the sample variance and the second dimension accounted for 23.3% of sample variance. The most influential variables for the first dimension are the distal humerus, distal radius and distal tibia ( Table 10). The two most important constraining variables are LLB completeness frequencies, which explain most of the variance in the first dimension. The tooth pit:tooth score ratio accounts for a substantial part of the variance of the second dimension ( Figure 11).
The RF test that employed all variables produced a solution in which the OOB error was stabilized after 230 trees, although data on brown bear-and spotted hyena-induced modifications (representing, respectively both extremes of LLB modification) required fewer than 100 trees ( Figure 12). The result produced an OOB estimate of the error rate of 30%. The solution is correct 70% of the time. The mean decreased Gini value showed that the most important variables were the tooth pit:tooth score ratio, distal humerus, LLB completeness frequencies, distal radius and proximal femur. Data for spotted hyenas and lions produced no error of classification, but, until the OOB was stabilized, data for wolves and brown bears were confounded up to 60% of the time (see also the proximity and overlap in PCA and CVA tests).
The RF analysis of LLB portion data only yielded a less stable OOB error (Figure 13), although it was similar to the previous analysis (30%). The mean decreased Gini value showed that the distal humerus and distal radius are the most discriminant variables, followed by the other portions. Brown bears signatures overlapped with those produced by wolves wolves .50% of the time until OOB stabilization was reached.
A CT with the selected set of variables from the first RF test shows that tooth pit:tooth score ratios of ,0.48 correctly identify lion-generated samples 100% of the time (Figure 14). If modifications on the distal humerus are .46.5%, this correctly classifies hyenas in 100% of the observations. If the damage on the distal humerus is ,46.5%, then only a combination of damage on the proximal femur and the distal radius can correctly differentiate bears from wolves in 100% of the observations. Further, our CT analysis of LLB portion variables reveals that modification distributions of the distal humerus differentiate spotted hyena-derived samples from those produced by the other three taphonomic agents that we investigated. Last, consideration of damage frequencies on the distal humerus, proximal radius and proximal femur, in combination, differentiates samples modified by brown bears wolves in 100% of the observations (Figure 15).

Conclusions
Previous neotaphonomic studies demonstrated that bears are efficient consumers of small-sized ungulate carcasses, resulting taphonomically in modest degrees of bone damage and a lack of medium-or long-distance transport of carcasses from points of initial discovery [36,53,55,56]. These results, along with those on the morphology, metrics and anatomical distribution of their tooth marks, prompted some workers (e.g., [55]) to characterize the taphonomic potential of bears as consistent with that of betterstudied lions. Data and results presented here disagree with this assessment. Most of the bear-imparted damage in our sample occurs not on LLB specimens, but instead on axial and girdle specimens.
That said, compared to previous neotaphonomic studies of bears (e.g., [53,55,56]), the Pyrenean brown bears did inflict a relatively high frequency of tooth punctures/pits/scores -on 14.4% of the total aggregated NISP. Consistent with the results presented by Haynes [53] in a different sample of bear-modified bones, most of the pits and scores in our sample are light and are positioned on thin cortical bone; there is a very modest degree of furrowing into the trabeculae of LLB epiphyses, but no substantial punctures through the thick cortices of LLB diaphyses (the only perforated bone in the Pyrenean bear-modified sample is a red deer rib [OB 15]). Based on the complementary results of these independent studies of bears, we feel comfortable generalizing that bears-unlike spotted hyenas and other durophagus carnivores-do not regularly use their teeth to puncture LLB shafts. Following from this conclusion, we remain cautious about well-known interpretations that fossil LLB diaphyseal specimens from the Middle Pleistocene cave sites of Yarimburgaz (Turkey) [31][32][33], Troskaeta (Spain) [36,89] and Divje Babe (Slovenia) [54] were punctured by extinct cave bears during the course of cannibalistic feeding. Although it is possible that bears at these sites fed with great intensity (perhaps under stressful environmental and/or individual conditions; see, e.g., our discussion of OB 15, above) and did, in fact, impart punctures on LLB diaphysis, we believe that, additional detailed analyses are required, since the punctured bear LLB diaphyses may have also been damaged by the taphonomic actions of durophagus carnivores with robust, conical tooth cusps, such as hyenas or wolves. Indeed, by refining our understanding of the taphonomic behaviors and consequences of extant brown bears, this study, in stepwise manner, contributes to ongoing efforts to model the possible impact of extinct ursids on conditioning paleontological and archaeological faunas in (especially) European Pleistocene contexts.
More generally, the results of our comparative statistical analyses draw the most precise distinctions between damage patterns imparted respectively by durophagus and non-durophagus carnivorans. By in large, these results corroborate commonsense and earlier qualitative observations that durophagus carnivorans can inflict much more intensive damage to bone and are more effective bone destroyers than are non-durophagus carnivorans. But the novel contribution of this research is that it provides quantitatively based, taphonomic thresholds for several commonly preserved paleontological variables-thresholds that can be used, with great accuracy and nearly 100% confidence, to separate faunas formed predominantly by durophagus carnivorans from those formed predominantly by non-durophagus carnivorans.
In summary, the observations and data that emerged from this study provide a clear pattern of ungulate carcass consumption by Pyrenean brown bears, regardless of bear age, sex or body size. But, considering the impact of local environment on conditioning taphonomic behavior (e.g., [14]), additional investigations of extant ursids in other locales and habitats is certainly merited. Ultimately, results from all these studies should be combined in order more completely characterize the taphonomic potential of extant brown bears and to better define their role as referents for modelling the taphonomic behaviors of their extinct cogeners.