Heavy reliance on plants for Romanian cave bears evidenced by amino acid nitrogen isotope analysis

Heavy reliance on plants is rare in Carnivora and mostly limited to relatively small species in subtropical settings. The feeding behaviors of extinct cave bears living during Pleistocene cold periods at middle latitudes have been intensely studied using various approaches including isotopic analyses of fossil collagen. In contrast to cave bears from all other regions in Europe, some individuals from Romania show exceptionally high δ15N values that might be indicative of meat consumption. Herbivory on plants with high δ15N values cannot be ruled out based on this method, however. Here we apply an approach using the δ15N values of individual amino acids from collagen that offsets the baseline δ15N variation among environments. The analysis yielded strong signals of reliance on plants for Romanian cave bears based on the δ15N values of glutamate and phenylalanine. These results could suggest that the high variability in bulk collagen δ15N values observed among cave bears in Romania reflects niche partitioning but in a general trophic context of herbivory.


Results
Judging from the C content (26.0-46.1%), N content (9.1-16.2%) and C/N atomic ratio (3.3-3.5), all the cave bear as well as the cave lion and horse bone collagen are well preserved (Table 1) 55,56 . Large variations in δ 13 C and δ 15 N values of bulk collagen for the cave bears were observed: −19.9 to −22.3‰ for δ 13 C and 5.2 to 9.8‰ for δ 15 N values, respectively (Fig. 2). The bears sampled from each cave corresponded to different time horizons. Out of six cave bear samples, two have been molecular dated (USR10 from Măgura cave and USR67 from Răsuflătoarei cave), and the other four have been radiocarbon dated (Table 1). Based on calibrated 14 C age and the mean of the posterior samples for molecular dated samples, ages of the samples from Măgura fell within the interval 28-31 ka BP, while the oldest bears were from Cioclovina Uscată with 14 C ages of c. 43 and 48 cal ka BP. The samples from Răsuflătoarei cover a wider time-span (20-49 ka BP), due to the wide confidence interval depicted based on molecular tip dating, given that USR67 is the only representative within its clade. In each cave and related timespan, both high (above 8‰) and lower δ 15 N values of bulk collagen have been recorded (see more details in Table 1).
Strong correlations between the δ 15 N values of several AAs and those of bulk collagen were observed for the investigated cave bears, among which phenylalanine showed the strongest one (R 2 = 0.926, P < 0.01). Other AA such as hydroxyproline (R 2 = 0.704, P < 0.05) also showed strong correlations with bulk collagen δ 15 N values. The scatterplot presenting δ 15 N Phe vs. δ 15 N Glx values shows that some cave bear individuals from the Romanian caves exhibited isotope values similar to those analyzed in previous studies for cave bears from Late Pleistocene Belgian sites, though three individuals with bulk collagen δ 15 N values above +8‰ analyzed in our study exhibited higher www.nature.com/scientificreports www.nature.com/scientificreports/ δ 15 N Phe values than the others without any overlaps. Most relevant for this study are the canonical trophic position (TP; 1 = primary producers, 2 = primary consumers, 3 = secondary consumers, etc…) estimates based on eqn. 1 (see Methods section) indicating that all cave bears, including those with high δ 15

Discussion
Feeding behavior of the cave bears based on AA δ 15 n values. The bulk collagen δ 13 C and δ 15 N values are within the range of published data for cave bears from Romanian caves in general ( Fig. 3) 24,26 . Thus, the isotopic results that were obtained on single AAs in the present study are likely to apply also to most other cave bears from Romania, including the ones claimed to have been omnivores, carnivores or piscivores, based on the high δ 15 N values of their bulk collagen. Our results rule out at least the following possibilities; (i) significant contribution of aquatic resources such as fish to their diets; (ii) carnivory as the main feeding behaviour 54 . The TP estimates for all individuals based on δ 15 N of glutamate and phenylalanine were around 2, which rather indicates a highly plant-dependent feeding behavior for all the analyzed cave bears, including those with a high δ 15 N value of their bulk collagen. The strongest correlation was observed between the δ 15 N value of bulk bone collagen and the δ 15 N value of phenylalanine that represents δ 15 N of the nitrogen source in a food web. This indicates that the bulk collagen δ 15 N values of the Romanian cave bears reported here can be explained by the following two possibilities, that are not mutually exclusive: (i) a δ 15 N shift of the baseline in this local ecosystem, and/or (ii) the consumption of some specific plants with high δ 15 N values by these cave bear individuals. This means that the reported high δ 15 N values in the Romanian cave bears' bulk collagen can be explained by the exclusive consumption of plants rather than by the trophic level effect caused by consumption of animal meat. As shown in Fig. 2, δ 15 N Glx and δ 15 N Phe values for some individuals, especially the ones with high bulk δ 15 N values, overlap with those of strictly herbivorous species, such as late Pleistocene mammoths from other regions, suggesting that a purely herbivorous species could exhibit such carbon and nitrogen isotopic values of bulk collagen in this palaeoenvironmental context 57,58 .
It was shown that there is currently a potential error of ~0.3 units for the TP estimation calculated by the propagation of errors for each parameter in eqn.1 59 . Considering this error together with a general offset of 0.1-0.2 (accuracy) between calculated TP based on eqn.1 and known TP that were estimated based on well-controlled field observations or laboratory feeding experiments for extant animals, our estimates for all individuals are well within the range of acceptable values for herbivorous mammals. In addition to predatory carnivores, fossils that so far yielded TP-values above 2.3 originate from brown bears, wild boars, foxes and badgers 54,60 , in contrast to the cave bears with TPs less than 2.2 analyzed in the present study.
The question remains which kind of plants could have been consumed by the cave bears with bulk collagen high δ 15 N values. Mammoth and some Pleistocene fallow deer have been found to exhibit higher δ 15 N value and lower δ 13 C values relative to the other coeval herbivorous species 43,44 . Some horses from the Late Pleistocene of western Germany also exhibited similar high δ 15 N and low δ 13 C values and are interpreted as having the same ecological niche as mammoths 61,62 . It is conceivable that the cave bears with high bone collagen δ 15 N values fed on plants, since the high values were also found in species with indisputable herbivorous dietary behavior. Possible candidates of such plants with high bulk δ 15 N value consumed by cave bears were those consumed by some fallow deer in the Pleistocene 44 , while those consumed by mammoth such as dry, mature grasses and sedges were less likely based on the lower resistance of cave bear teeth against abrasion. Mushrooms, another plant food with high δ 15 N values, can be also excluded as an explanation for the high δ 15 N values of some cave bears for the two following reasons: (1) mushrooms typically have high δ 15 N and high δ 13 C values (not low δ 13 C values), (2) mushrooms are consumers, like animals, and they display a trophic position value higher than 2 52 , therefore their consumption would lead to a TP significantly higher than 2. Ellipses indicate 90% and 50% confidence intervals for isotopic compositions of the cave bears analyzed in this study; except for a few individuals like the one with values plotted at the right-upper corner, many cave bears with δ 15 N value of collagen around 7-8‰ that might be interpreted as an omnivorous/carnivorous signal can be categorized to highly-plant-dependent feeders.
Scientific RepoRtS | (2020) 10:6612 | https://doi.org/10.1038/s41598-020-62990-0 www.nature.com/scientificreports www.nature.com/scientificreports/ One additional factor specific to bears that has to be considered is hibernation, since unlike the other mammalian taxa discussed above, including Homo sapiens, bears hibernate. Indeed, it has been demonstrated that a longer hibernation period caused by climate cooling could result in higher δ 15 N of cave bear bulk collagen 63,64 . Although we did not see any statistically significant chronological trend in δ 15 N values of AAs, only a limited number of individuals have been analyzed so far with this approach. Further investigations on the possible influence of this factor on the δ 15 N values of individual AAs of Ursids should be performed on more extinct cave bears as well as extant bears in order to clarify how nitrogen metabolism such as urea recycling during hibernation possibly affects the δ 15 N of AAs in bone collagen 65,66 . Yet, our data currently does not support a major physiological control of the high bulk collagen δ 15 N value in bone that remodels intensively during the warm season 30 . evolution of feeding behavior and niche partitioning among cave bear populations. Despite a predominantly herbivorous diet, Romanian cave bears exhibit a very large variation of their isotopic values that probably reflects ecological differences. Such isotopic differences in bone collagen reflect long-term dietary differences, since the turnover of collagen is slow in bone, averaging diet isotopic composition over many years 67 . It is possible that, in contrast to other regions of Europe, the absence of mammoths played a role in allowing cave bears to have access to a more diverse feeding niche and that competition was stronger among cave bears in Late Pleistocene Romania than in other regions, although an exact contemporaneity between the populations with diverse bulk collagen δ 15 N values must be examined. If the observed variability in δ 15 N values of bulk collagen among the Romanian cave bears was the result of high competition among cave bears in a context of a feeding niche restricted to herbivory, it is possible that some individuals specialized on other plants rather than turned to omnivory including animal food resources. Such a dietary strategy is a very unusual case among bears. Assuming that the herbivory for individual cave bears observed in our study was a general trait across all Late Pleistocene cave bear taxa, the most parsimonious interpretation would be the loss of omnivorous feeding behavior after the divergence of cave bears from the brown/polar bear lineage and before the basal divergence of cave bear lineages (e.g., U. ingressus and U. spelaeus; see also low bulk δ 15 N values for U. kudarensis 68 ). However, it still remains an open question why this extinct animal that could live over a wide geographical range reaching Transbaikalia in Eastern Siberia with diverse ecological conditions within the mammoth steppe ecosystem finally became extinct 69,70 . It is possible that a dietary adaptation similar to that seen in extant giant panda that consumes large quantities of bamboos but has a macronutrient profile, namely the percent of energy derived from protein, similar to that of carnivores, might have facilitated the full shift for the cave bear lineage towards highly plant-dependent feeding behavior while keeping some ecological plasticity in the mammoth steppe 71 .
Our findings have interesting implications regarding the evolution of herbivory among carnivorans. The foraging strategy of an animal is largely determined by body mass via energetic needs 72 . Cave bear populations were likely adapted to their local habitats, with different altitudes being associated with different features, such as body size, tooth morphology, and tooth wear patterns (e.g., ref. 21,33,73 ). Although the tooth microwear and the isotopic composition of bone collagen reflect the diet of the analyzed individuals at different pre-mortem time periods (e.g., several weeks vs several years or more before the death 74 ), the observed variation in the bulk collagen δ 15 N values for cave bears in the Romanian region possibly reflected a difference in foraging strategy on plants associated with different body sizes 24 . Extant Carnivora species with herbivorous feeding behavior (giant panda, red panda and binturong) have smaller body size and live only in temperate and tropical/subtropical Asia, highlighting the paradox of cave bears that had much larger body sizes and lived in a colder and drier environment. Besides, these three extant species are all on very long evolutionary branches, namely they diverged more than 10 million years ago from their closest relatives. Cave bears, in contrast, share a much more recent ancestor within an omnivorous/carnivorous clade, which is another ecological paradox of this extinct taxon. It is possible that the modern herbivorous carnivorans are relicts of a dietary trend that was much more widespread in the past, and that their rarity is rather due to the recent extinction of the large mammalian species rather than to the impossibility for Carnivora to have a diet restricted to plant food. Our results emphasize the need to investigate more in detail the diversity in food sources, especially plants, consumed by cave bears, as a unique case among Carnivora to further understand the benefits, costs and limitations of herbivory. conclusions A highly plant-dependent feeding behavior for cave bears was demonstrated by the TP estimates based on δ 15 N Glx and δ 15 N Phe values. The high δ 15 N of bulk bone collagen of some cave bears can be explained by a higher δ 15 N baseline of the exploited food chain that is reflected in the δ 15 N Phe values, namely plants with high δ 15 N values. Therefore, the high δ 15 N value of bulk collagen of some cave bears in the Romanian region, a unique feature compared to cave bears from the other regions of Europe analyzed so far, can be explained by the consumption of plants. Further research is needed to clarify which plant species were eaten by those cave bears by using other approaches such as δ 13 C analysis of individual AAs to identify protein sources 75,76 . It is also of importance to see if this trend towards herbivory was common for other cave bear populations from broader regions including Asia. www.nature.com/scientificreports www.nature.com/scientificreports/ Bone collagen from 6 adult cave bears from the three caves was subjected to δ 13 C and δ 15 N analysis as well as AA-specific δ 15 N analysis. One cave lion collagen sample from the North Sea, an area that was terrestrial during the Late Pleistocene and one horse sample from the site of Lommersum (Germany) from the late Pleistocene, have been analyzed together with the cave bear collagen using the same analytical protocol.

Collagen extraction and bulk isotopic analysis. The samples were washed in an ultrasonic bath in
acetone, rinsed several times with demineralized water, dried at 35 °C for 72 h and crushed to powder of 0.7 mm grain size. Collagen extraction was performed after ref. 78 with a minor modification as described in ref. 79 . In summary, this acid-base-acid method includes a first step of demineralization with 1 M HCl, a second step in 0.125 M NaOH to remove humic acids and lipids, and a final step of gelatinization at pH 2 during 17 h at 100 °C. Isotopic measurements were done using an NC2500 elemental analyzer connected to a Thermo Quest Delta þ XL mass spectrometer. Using the 'δ' (delta) value, the isotopic compositions are expressed as follows: δ (‰) ≡ 10 3 [R sample /R standard − 1], where the R denotes the 13 C/ 12 C ratio for carbon and the 15 N/ 14 N ratio for nitrogen, with the international reference (standards) being V-PDB for δ 13 C values and atmospheric nitrogen (AIR) for δ 15  The GC oven temperature was programmed as follows: isothermal hold at 40 °C for 3 min; temperature ramp to 110 °C at 15 °C min −1 ; ramp to 150 °C at 3 °C min −1 ; ramp to 240 °C at 2.5 °C min −1 ; ramp to 270 °C at 6 °C min −1 ; and subsequent holding isothermally at 270 °C for 4 min. Carrier gas (He) flow rate through the GC column was 1.4 ml min −1 . The CO 2 generated in the combustion furnace was eliminated by a liquid nitrogen trap. Standard mixtures of ten AAs with known δ 15 N values were injected into the GC/IRMS every two to five runs to confirm reproducibility of the isotope measurements. The precision of the reference mixtures was 0.65-0.90‰. Nitrogen isotopic composition of maximally the following eight AAs was determined as their N 2 gas peaks from derivatives showed fine resolution on the GC/IRMS chromatogram: valine (Val), leucine (Leu), isoleucine (Ile), proline (Pro), serine (Ser), glutamate (Glx), phenylalanine (Phe) and hydroxyproline (Hyp). All reported δ 15 N values for glutamate included the contributions from the α-amino group of glutamic acid and glutamine, as glutamine is converted to glutamic acid during acid hydrolysis. Data for alanine and glycine were not obtained for all samples due to relatively large peak sizes in this study.
Radiocarbon dating. Cave bear collagen was dated directly using AMS radiocarbon dating. Radiocarbon dating was performed by the Labor für Ionenstrahlphysik at Eidgenössische Technische Hochschule-ETH in Zurich (Switzerland). The obtained radiocarbon ages were calibrated using the OxCal v4.2.4 software using the IntCal13 atmospheric curve 82,83 . All dates were calibrated to BP dates with 2σ (95.4%) probability.

inferring molecular ages
Ancient DNA extraction and library building. For two samples lacking radiocarbon dates (USR10 and USR67), molecular tip-dating was performed. DNA was extracted from 50 mg of cave bear petrous bones following the protocol of ref. 84 , with reduced centrifugation speeds as described in ref. 85 . Prior to library building, each DNA extract has been quantified with the Qubit 2.0 fluorometer with high sensitivity reagents (Thermo Fisher Scientific) by using one microlitre out of 25 μL extract. To support the validity of the measurements, positive and negative controls were also measured. Illumina single-stranded libraries were then prepared following the protocol of ref. 86 , after the removal of uracil residues and abasic sites. The optimal number of library amplification PCR cycles has been assessed with qPCR as described in ref. 85 . A volume of 80 μL including 20 μL template library, Accuprime Pfx DNA polymerase and eight base-pair tailed-primers was used for Indexing PCR to generate dual-indexed library molecules. Amplified libraries were purified using commercial silica spin-columns (Qiagen MinElute). Prior to single-end sequencing on Illumina platforms, library concentration and length distribution were quantified using Qubit 2.0 and 2200 TapeStation (Agilent Technologies), respectively.
Sequencing and bioinformatic analyses. Shotgun sequencing was performed on an Illumina NextSeq 500 sequencing platform. We obtained 75 bp single-end reads using the NextSeq 500/550 High Output Kit v2 (75 cycles) sequencing kit, following the procedures described in ref. 87 . Raw reads were trimmed with minimum overlap of one nucleotide and reads shorter than 30 bp were discarded using cutadapt v1.12 88 . The processed reads Scientific RepoRtS | (2020) 10:6612 | https://doi.org/10.1038/s41598-020-62990-0 www.nature.com/scientificreports www.nature.com/scientificreports/ were mapped to the reference mitogenome sequence of Ursus spelaeus (Genbank Acc. No. EU327344, ref. 89 ) with default parameters using the bwa v0.7.15 "aln" algorithm 90 . The alignment was then filtered for mapping quality (Q > 30), sorted by read position and potential PCR duplicates were removed using SAMtools v0.1.19 90 . Mitochondrial genome reconstruction and alignment. ANGSDv0.920 was used to generate consensus sequences for both individuals 91 . Using the option -doFasta 3 that takes into account bases with highest effective depth 92 , the mitochondrial consensus has been constructed only considering read mapping and Phred base quality scores >30. The two mitogenomes were aligned to a published alignment of cave bear mitochondrial sequences 15 . The repetitive section of the d-loop was removed from the alignment because it could not be aligned unambiguously. The final alignment comprised 19 sequences of 16,468 nucleotides, publicly available (Dryad database and GenBank). The GenBank accession numbers of the mitogenomes generated in this study are: MN311249 for USR10 and MN311250 for USR67. Molecular tip-dating. Ages of the two samples were estimated using Bayesian phylogenetic tip dating 93 , a method shown to provide accurate age estimates for Late Pleistocene European cave bears 15 . PartitionFinder v.2.1.1 has been used for data partitioning and substitution model selection 94 . All substitution models available in BEAST v.1.8.2 were considered and the best partitioning scheme has been computed under the Bayesian Information Criterion using linked branch lengths and the greedy search algorithm (Supplementary Table 2) 95 .
Tip-dating was performed separately for each sample by analyzing its mitochondrial sequence in BEAST v1.8.2, using median 14 C dates as calibration from the dataset in ref. 15 (Supplementary Tables 3). For each tip date its posterior distribution has been sampled using a uniform, uninformative prior of zero to one million yBP. To accommodate changes in effective population size during the time-span of the tree, a piecewise-constant coalescent Bayesian Skyline population model was selected. Unlinked strict molecular clocks were set for each data partition with the posterior distribution of the substitution rate of each sample with a uniform, uninformative prior of zero to 20% mutations per million years. The software Tracer v1.6 was used to assess if the MCMC chain had run for a sufficient length to achieve burn-in and adequate sampling of all parameters (ESS > 200) 96 . The Maximum Clade Credibility Tree was selected from the posterior sample using TreeAnnotator and visualised in FigTree ( Supplementary Fig. 1).