Isotopic signature in isolated south-western populations of European brown bear (Ursus arctos)

Stable isotope analysis of animal tissue samples is increasingly used to study the trophic ecology of target species. The isotopic signatures respond to the type of diet, but also to the environmental conditions of their habitat. In the case of omnivorous, seasonal or opportunistic feeding species, the interpretation of isotopic values is more complex, as it is largely determined by food selection, either due to individual choice or because of availability. We analysed C and N isotopes in brown bear (Ursus arctos) hair from four isolated populations of south-western Europe (Cantabrian, Pyrenees, Central Apennines and Alpine) accounting for the geographical and climatic differences among the four areas. We found inter-population differences in isotopic signatures that cannot be attributed to climatic differences alone, indicating that at least some bears from relatively higher altitude populations experiencing higher precipitation (Pyrenees) show a greater consumption of animal foods than those from lower altitudes (Cantabrian and Apennines). The quantification of isotopic niche space using Layman’s metrics identified significant similarities between the Cantabrian and Central Apennine samples that markedly differ from the Pyrenean and Alpine. Our study provides a baseline to allow further comparisons in isotopic niche spaces in a broad ranged omnivorous mammal, whose European distribution requires further conservation attention especially for southern isolated populations.


Introduction
The brown bear (Ursus arctos) is the largest carnivore of Europe. Its geographical range includes the whole Holarctic realm and the IUCN currently considers this species as "Least Concern". Within Europe, populations of brown bear are increasing in number although a great disparity in population size occurs between the south and the north of the continent (Chapron et al. 2014). Due to their small size, highly fragmented and human-impacted landscapes, bear populations from Italy and Spain/France are of a particular conservation concern ). In the first half of the twentieth century, the bear population in the Cantabrian mountain range of northern Spain split into two subpopulations (western and eastern). Estimates of population size in the 1990s were no more than 60 individuals for the western subpopulation (Wiegand et al. 1998) and about 14 individuals in the eastern subpopulation (Clevenger and Purroy 1991). The most recent estimates from genetic sampling give a figure of 223 bears, most in the western subpopulation and only 19 in the eastern subpopulation of the Cantabrian Mountains (Pérez et al. 2014). In the Pyrenees, bears were almost extinct in 1995 with an alarmingly reduced population of only 5 individuals (Camarra and Dubarry 1997). After the introduction of 2 males and 6 females from Slovenia in 1996-1997 and 2006, 30 individuals were counted in 2015 (Piédallu et al. 2016) and a minimum of 46 detected in 2017 (Parres et al. 2020). Italian brown bears persist with two different and isolated populations. One, in the central Apennines, coincides with the 1 3 relict population of Apennine brown bears (Ursus arctos marsicanus), numbering around 50 individuals (Ciucci et al. 2015) and critically endangered according to the regional IUCN criteria (Rondinini et al. 2013). The second, in the central Alps, originates from the successful re-introduction in Trentino and currently comprises a minimum of 43 individuals (Tosi et al. 2015) with no connectivity with the neighbouring Slovenian bear population.
Concerns related to the long-term future of these populations remain and major conservation efforts have been put in place, especially to protect the endemism of the Apennine (Thomsen et al. 2021) and Cantabrian Mountain populations (Rodríguez et al. 2007). Within this context, tools are now needed to understand differences in ecology and behaviour of these isolated populations to ensure their functional role within south-European landscapes. In addition, the particularities of each population must be known in order to adapt specific conservation policies (Martínez-Abraín et al. 2021).
Previous genetic evaluations suggested that these southwestern bear populations reflect glacial refugia (Taberlet and Bouvet 1994); however, a recent genomic approach confirmed, at least for the Apennine population, an origin of about 1500 years which in turn resulted in the remarkable evolution of specific phenotypic and behavioural traits (Benazzo et al. 2017). Studies of ancient DNA also do not support the hypothesis of glacial refugia (Ersmark et al. 2019). In the north of the Iberian Peninsula, the Iberian Pleistocene lineage of bears was replaced by one coming from a cryptic Atlantic refuge that entered through the Pyrenees after the Late Glacial Maximum ).
Available information on the feeding ecology of all these bear populations confirms a strong element of seasonality and opportunism, which is typical of the species (Swenson et al. 2021). Cantabrian and Pyrenean brown bears have a diet generally characterised by the major consumption of vegetable compared to animal matter (Couturier 1954;Clevenger et al. 1992;Braña et al. 1993;Naves et al. 2006;Rodríguez et al. 2007). In the Cantabrian Mountains, the food most frequently eaten in spring is grass, although carrion consumption also increases compared to the rest of the year. In summer, the dominant food is fleshy fruits, and in autumn and winter, hard mast (nuts). The season in which animal protein is most important is summer, although the volume of insects (ants, for example) in their diet is almost double that of mammals (Braña et al. 1993). Pyrenean brown bears have been reported to consume more animal-derived food in addition to tubers (such as Conopodium majus), roots and, occasionally, aquatic animals such as frogs and trout (Couturier 1954). The diet of the Apennine brown bears was also described as being characterised by major consumption of a seasonally diversified array of vegetables (Zunino and Herrero 1972;Ciucci et al. 2014). In terms of annual energetic contribution, the ranking for the Apennine bear population is as follows: (i) hard mast (beechnuts and acorns) and fleshy fruits, (ii) herbaceous vegetation and insects (mostly ants), (iii) wild ungulates and livestock, and (iv) roots (mainly carrots). Hard mast is particularly important in autumn during hyperphagia, with fleshy fruits (especially buckthorn berries) contributing more during late summer. Summer data support higher consumption of herbs, forbs and animal proteins (wild ungulates, particularly red deer as expected by their availability, Ciucci et al. 2014;Careddu et al. 2021).
For the re-introduced Alpine population, preliminary data based on DNA metabarcoding support a great consumption of vegetable matter (De Barba et al. 2014) with high frequency of occurrence for plants belonging to the Asteraceae, Apiaceae and Maleae families. Damage reports by individual bears also suggest that livestock consumption can be relatively high although a quantitative incidence of this phenomenon in the diet of the whole population is lacking (Groff et al. 2014. Altogether, these data on the southwestern populations of brown bear support a highly flexible ecology in relation to food availability with predominant consumption of vegetable matter, as expected from brown bear geographic diet variation (Bojarska and Selva 2012).
Using methods proposed by Turner et al. (2010) and Jackson et al. (2011), we analyse the isotopic niche of these populations testing the hypothesis that isotopic niche should differ between geographically distinct populations. Stable isotope analysis (SIA) is an important complementary method to assess diet, feeding ecology and ecophysiology of mammals (Crawford et al. 2008;Ben-David and Flaherty 2012). For dietary studies, carbon and nitrogen isotopes are most often used. In terrestrial ecosystems, the carbon isotopic signature is primarily related to the proportion of C3 or C4 plants in the diet, as the two types of plants differ in the way they take up CO 2 during photosynthesis and have divergent δ 13 C values (O'Leary 1988;Farquhar et al. 1989). In Europe, C4 plants are not common. Its percentage ranges from 1.1 in central Europe to 5.6 in south-western Europe, and is even lower in mountainous areas with temperate vegetation, where the arid conditions that favour its development do not exist (Collins and Jones 1986;Pyankov et al. 2010). However, in some regions there is extensive cultivation of maize (Zea mays), a C4 plant, which can be consumed by wildlife directly or indirectly (through the intake of livestock fed on this type of crop). If ingested regularly, its isotopic signature will be appreciably recorded in the tissues of the consumer. Inputs of marine protein in diet produce also a distinctive δ 13 C signature (Chisholm et al. 1982) but, unlike their North American relatives that consume seasonally migrating salmon, there is no evidence of marine protein intake in European brown bears (Swenson et al. 2007). Other sources of δ 13 C variation are related to climatic and geographical variables, such as rainfall, temperature, insolation, altitude or tree cover, with sometimes cumulative, sometimes opposing effects on isotopic values (Dawson et al. 2002;Körner et al. 1991;Diefendorf et al. 2010). In herbivores, the mean observed diet-hair offset in δ 13 C was + 3.2‰, with a range of + 2.7 to + 3.5‰ (Sponheimer et al. 2003a). Each step in the food chain, from the primary consumer onwards, produces an additional increase in δ 13 C up to 1‰ (Bocherens and Drucker, 2003).
The nitrogen isotopic signature in plants also depends on geographic and climatic factors related to temperature that influences the microbial activity in soils and thus the fixation of atmospheric N. To a lesser extent, it depends also on altitude and rainfall (Amundson et al. 2003). The dietary isotopic offset is more pronounced than for carbon. δ 15 N increases by 3.5 to 5‰ at each step of the trophic chain from primary up to tertiary consumers (Bocherens and Drucker 2003), although in various types of herbivores there are notable differences between species (Sponheimer et al. 2003b).
There have been a number of studies exploring isotopic variation in bears as a proxy to discriminate diet between species and populations (Hildebrand et al. 1996(Hildebrand et al. , 1999Hobson et al. 2000). Isotopes vary considerably between different tissues, with hair being demonstrated to record bear's diet during the hair growth period, without undergoing subsequent turnover, unlike, for example bone collagen (Hilderbrand et al. 1996). Meta-analyses of the North American grizzly (U. arctos horribilis) based on isotopes extracted from hair demonstrated the expected high incidence of salmon diet in populations closer to the coastline, while terrestrial prey (wild ungulates) represent a larger portion of grizzly diet from mountainous and more internal regions (Mowat and Heard 2006). For European brown bear, with the exception of Apennine bears (Careddu et al. 2021), there are currently no extensive isotopic data published based on hair even if this is relatively available in the field and has been targeted for non-invasive genetic analyses (Pérez et al. 2009;De Barba et al. 2014;Gervasi et al. 2012). Vulla et al. (2009) proposed an ecogeographical gradient for representative omnivorous carnivores (including bears), suggesting that at higher latitudes omnivorous species should consume more animal matter. Accordingly, we hypothesised that the nitrogen isotope ratio should increase following a geographical gradient from the Apennines, through Pyrenees, then Cantabrian Mountains and finally central Alps. Other factors that may have an influence are geographical and climatological, as there are differences in altitude between the populations studied, and subtle climatic variations. Considering the high dependency of all these populations on plant food, we predict a geographical effect both in nitrogen and carbon isotopic signatures reflecting those variations.

Specimen collection
The data set comprises hair (n = 32) from 20 distinct locations opportunistically collected during genetic non-invasive surveys or from museum specimens (Fig. 1). Hairs from the Cantabrian population (n = 14) were collected in hair traps or rub trees between 2007 and 2014, aside from one female sample collected in 1991. The sex of nine individuals was known.
The sample from the Pyrenees included three specimens from the Museu de Ciènces Naturals de Barcelona of the early twentieth century (of which one was male, one female and one juvenile), and one bear re-introduced from Slovenia in 2011. The Italian samples were all collected from hair snag or rub trees sites both in the Apennines (2001, n = 7) and central Alps (2013, n = 7). Specimens from the Apennines were not individually recognised while all the hairs from Alpine population belonged to 6 adult males and 1 female that were constantly monitored in the Adamello Brenta National Park.

Sample preparation and isotopic extraction
We collected at least 0.5 g of guard hair for each specimen, excluding the scalp. The use of whole hairs makes it possible to obtain an isotopic signal averaged over the entire period of hair growth, avoiding the difficulty of assigning a specific time interval for each section of hair, since hairs from different parts of the body do not grow at the same rate or reach the same length.
SIA was carried out in three laboratories: Institute of Geology of the University of A Coruña, Spain (IUX-UDC; n = 27), NERC National Environmental Isotope Facility, British Geological Survey (Nottingham, UK, NEIF; n = 4, © UK Research and Innovation 2022) and the University of Liverpool, UK (n = 2).
At IUX-UDC, dirt, natural oils and museum preservatives were removed from the samples following a standard protocol: the hair was sonicated for 1 h in a 2:1 methanol:dichloromethane solution (O'Connell et al. 2001), after which the procedure was repeated. This was followed by two 20-min rinse in Millipore water, also in a sonicator. The hair was then left to dry in a desiccator chamber. Since the hairs grow over several months during which the diet may vary, the samples were homogenised by cutting the entire length into segments of less than 1 mm and mixing them for analysis. Two sub-samples of approx. 200 mg of guard hair from each animal were weighed in tin capsules using a UMX-2 balance (Mettler Toledo). The determination of δ 15 N and δ 13 C was carried out by combustion in an EA1108 elemental analyser (Carlo Erba Instruments) linked via a ConFlo III interface to a MAT253 isotope ratio mass spectrometer (ThermoFinnigan). In each analytical sequence, USGS 40 (− 4.52‰), USGS41a (+ 47.55‰), IAEA-N-1 (+ 0.4‰), IAEA-N-2 (+ 20.3‰) and USGS-25 (− 30.4‰) were used as secondary standards for δ 15 N. For δ 13 C, USGS 40 (− 26.39‰), USGS41a (+ 36.55‰), NBS 22 (− 30.031‰) and USGS 24 (− 16.049‰) were used. Replicates (n = 10) of an internal standard (acetanilide) evaluated the precision of the measurement (standard deviation), resulting in ± 0.15‰ for both C and N. The preservation of hair keratin was measured by its atomic C:N ratio, which should be between 2.9 and 3.8 (O'Connell et al. 2001).
At the NEIF, each analysis was performed on a composite of one to six whole hairs (equivalent to 0.6 mg for combined δ 13 C and δ 15 N isotope ratio measurement). Isotope ratios of carbon and nitrogen were measured by continuous flowelemental analysis-isotope ratio mass spectrometry (CF-EA-IRMS). The instrumentation was a ThermoFinnigan continuous flow system, comprising a Delta Plus XL isotope ratio mass spectrometer interfaced with a "Flash/EA" elemental analyser via a ConFlo III interface. All samples were analysed in triplicate. Hair carbon and nitrogen isotope ratios were calibrated using an in-house reference material Hair-2 (calibrated against CH7, USGS40, USGS41, N-1 and N-2, IAEA) and using additional check standards of nylon and Hair-1 (calibrated against CH7, USGS40, USGS41, N-1 and N-2, IAEA). Replicates (n = 17) of an internal standard  evaluated the precision of the measurement (standard deviation), resulting in ± 0.03‰ for carbon and ± 0.12‰ for nitrogen.
At University of Liverpool, washed (Milli-Q water 18 MΩ cm −1 ) de-fatted (2:1 methanol:dichloromethane) dry hair samples were ground (liquid N 2 ) and weighed into tin cups. They were then analysed in duplicate using an elemental analyser (Costech) coupled to a Delta V Advance mass spectrometer (Thermoquest). Samples were combusted at 1000 °C and diluted with N 2 prior to entering ECS 4010 CH/N/CN reaction tube and the mass spectrometer. USGS24, USGS40 and USGS41 were used to determine the accuracy (δ13C < 0.1‰, δ15N < 0.5‰) and precision (δ13C < 0.05‰, δ15N < 0.7‰) of the EA/IRMS and to correct the isotopic values obtained from the bear samples.
Isotopic values were extracted at least twice for each specimen and to ensure consistency across the labs, six hairs collected from the same locations were analysed by at least two labs separately and values compared. For carbon values, standard deviation ranged for the same specimens between 0.065 and 0.71 (mean st.dev = 0.41) while for nitrogen the range was between 0.18 and 0.63 (mean st.dev = 0.30). These data were within the range of intra-population variation and were smaller than inter-population standard deviation (mean st.dev δ13C = 0.91; mean st.dev δ15N = 1.65). The whole raw data is available in the Electronic Supplementary Material.

Calibration curve
Due to samples being collected in different time periods, we accounted for the Suess effect caused by the contribution of CO 2 depleted in δ 13 C due to the burning of fossil fuels (Keeling 1979) using a calibration curve. We used data from Friedli et al. (1986) (covering years 1744 to 1953 aD), Keeling et al. (1989Keeling et al. ( ) (years 1978Keeling et al. ( to 1988, Leuenberger et al. (1992) (40,000 BP to 1270 aD), Francey et al. (1999) (years 1006 to 1997 aD) and Graven et al. (2017) (years 1851 to 2015 aD). All these data cover the last 40,000 years, and specially the last 100 years, across the range of our samples. The data were fitted to an exponential curve with the statistical programme PAST (Hammer et al. 2001). The line obtained was as follows: where x is the year of collection before 2019, being 2019 = 0. This curve showed a good fit (R 2 = 0.97) and allowed us to calibrate all the carbon isotopic values to the year 2014, which corresponds to the most recent date of our samples.

Data analysis
The final dataset consisted of 32 individual samples from 16 distinct locations for which latitude and longitude were available (Supporting Information 1). To test the hypothesis that isotopic signature differs among populations, we first computed Layman metrics in the R environment using packages "SIBER" and "siar" (Layman et al. 2007;Turner et al. 2010). These include total isotopic niche area (TA, used as a measure of the total foraging width of a population, Layman and Allgeier 2012), mean distance to the centroid of each population (CD), and the eccentricity of the scatter of isotopic values in the isoscape (E). Eccentricity is used to determine whether the points scatter similarly in all direction, in which case it will be close to 0.0, or if the points scatter in a linear fashion, with this value approaching 1.0 (Turner et al. 2010).
To compare niche area and overlap among the brown bear populations, a Bayesian statistical model [Stable Isotope Bayesian Ellipses in R (SIBER)] was applied to the isotopic data to compute and delimit stable isotope Bayesian ellipses, correcting ellipses areas (SEAB) for sample size in each group, and measure the level of overlap between populations (Turner et al. 2010).
To test the hypothesis of geographical/environmental gradient in the isotopic values, each geographic location was characterised by latitude, longitude, altitude and additional climatic data (precipitation and temperature) combined into 19 bioclimatic variables as described by Fick and Hijmans (2017). These climatic parameters are representative of historical climatic conditions that covers 30 years between 1970 and 2000 and have been proved to effectively characterise eco-physiological adaptations of small (e.g. Moreno-García and Baiser 2021) and large (e.g. Di Marco et al. 2021) mammals. Localities with the same bioclimatic parameters were combined into one data point and this reduced the sample to n = 14 distinct geographic samples. Carbon and nitrogen values were equally averaged for each unique location.
The 19 bioclimatic variables together with altitude were first reduced through a principal component analysis (named ENV PCA) of the correlation matrix (recommended when variables have different magnitudes, Hair 2010), in order to account for the high level of correlation and reduce data dimensionality (Jiang 2018) using PAST (Hammer et al. 2001). Latitude and longitude data were transformed into a truncated geographic distance matrix subjected to a principal y = −6.46 − 2.22 * 0.98 ⋀ x coordinates of neighbourhood matrix (PCNM) using the package vegan (Oksanen et al. 2022). This procedure generates PCNM eigenfunctions that represent a spectral decomposition of the spatial relationships among our 14 sampled locations and accounts for spatial bias in statistical models (Borcard and Legendre 2002). Linear models were built in order to identify the best predictors (ENV PC scores and PCNMs) of carbon and nitrogen values across the 14 unique geographical samples. Dependent variables were initially identified using a forward selection procedure through the function "forward. sel", package ade.spatial (version 03-16, Dray et al. 2022). This was applied in two separate instances using carbon or nitrogen as separate independent variables. Multiple linear models of several complexities were compared using Akaike information criteria for small samples (AICc) and delta AICc (Brewer et al. 2016) to identify the best combination of ENV PCs and PCNMs (independent variables) that predicts variation into carbon or nitrogen isotopic data (dependent variables).

Results
The preservation of hair keratin was adequate in all samples, even those from museum preserved pelts, with a mean atomic C/N ratio of 3.4 (sd = 0.1, range between 3.2 and 3.6). The δ 13 C values extracted from the bear hair sample ranged from − 23.15 to − 20.15‰ (mean ± SD: − 22.33 ± 0.89‰, n = 32) while δ 15 N values for the bears ranged from 1.1 to 9.2‰ (mean ± SD: 5.56 ± 1.65‰, n = 32).
The isotopic values averaged by populations show a large disparity especially for the Alpine and Pyrenean sets. The standard ellipses of these two populations are also large. This is to be expected in the Pyrenean population due to the scarcity of data and their different chronology (Fig. 2, Table 1).
The Bayesian estimate for SEA (= SEAb) differed significantly between populations (one-way ANOVA, F = 32.678, df = 3, 431, p = 0001). However, post hoc multiple comparisons, based on Tukey, identified overlap between the Pyrenean and Alpine populations and between the Cantabrian and the Apennine (Table 2). These results concur with the density plot (Fig. 3), which shows the modes for standard ellipse area of the Pyrenean and Alpine populations at similar levels and the Cantabrian and Apennine populations at similar levels on the y axis.
The polar histograms and the polar density plots (Fig. 4) show that the Cantabrian and Apennine populations occupy the most similar isotopic niche space. This matches the isotope biplot, which showed the Apennine standard ellipse and convex hull fitting within those of the Cantabrian. The Pyrenean population has the highest position on the δ 15 N axis of all the populations, and along the δ 13 C axis was higher than the Cantabrian and the Apennines. Generally, the Alpine population has slightly lower δ 15 N and larger δ 13 C than the Cantabrian and Apennine populations.
The estimated area of overlap between each pair of ellipses was calculated independently for 1000 draws. When comparing the Pyrenean and Alpine populations and the Apennine and Cantabrian populations, most of the estimates for the overlaps occupy the 20-30% range; however, some estimates range up to 80% (Table 3).
The standard ellipses for the Cantabrian and Alpine populations are likely distinct and occupy different niches. Similarly, for the Pyrenean and Apennine populations, the highest frequencies occupy the 5-10% range, and none of the estimates reaches higher than 50% (Table 3).
The smallest overlap estimates were between the Apennine and Alpine populations with the highest frequency ranging between 0 and 2%, supporting significantly different isotopic niches. The Cantabrian and Pyrenean populations' histogram had a frequency range of 0-60% and the highest frequencies were for 15-20%, showing some cross-over (Table 3). Overall, these overlap estimates suggest that there is some level of distinction between the standard ellipses of all the populations, the least being between Pyrenean and Alpine and Cantabrian and Apennine and the highest being between Apennine vs Alpine. The 14 unique isotopic sampling locations could be clearly separated through ENV PCA that extracted 13 PC vectors of which the first two explained cumulatively 81.80% of variance (Fig. 5A, Supplementary Information 1). ENV PC1 (45.49% var.) was positively loaded on altitude (r = 0.78) as well as multiple bioclimatic parameters associated with precipitation such as BIO12 (r = 0.93), BIO13 (r = 0.90), BIO14 (r = 0.86) and BIO16 (r = 0.96) (see also Supplementary Information 1 for additional loading values). This axis separates alpine locations (negative ENV PC1 scores) from the rest due to their relatively lower rainy precipitations and some samples from low altitudes (Fig. 5A). ENV PC2 exhibit an opposite trend being negatively loaded on altitude (r = − 0.54) and positively on temperature parameters (BIO1 r = 0.76, BIO2 r = 0.76, BIO3 r = 0.81, BIO6 r = 0.96, BIO9 r = 0.82, BIO11 r = 0.96, see also Supplementary Information 1) and it distinguished Alpine and Pyrenean locations (negative scores, lower temperatures) from the Apennine and Cantabrian mountains range (positive scores).

Isotopic niche analysis
We identified significant geographical differences in isotopic signature of south-western European brown bear populations. The geographical isolation appears to have impacted ecological habits in the brown bear that Cantabrian brown bears inhabit steep highlands with sparse tree cover and feed predominantly on vegetable matter . Similarly, Apennine brown bears feed on predominantly vegetable matter such as herbs and fleshy fruits. They have an overall low consumption of large mammals, which peaks in the spring/early summer. During summer and autumn, males have higher δ 15 N values than females (Careddu et al. 2021).
We found that the Apennine and Cantabrian populations are likely to have similar isotopic niches. This may be due to their predominantly herbivorous diets. Field data has shown that both Cantabrian and Apennine brown bear faeces contain predominantly plant material with scavenged animal protein supplementing their diet (Clevenger et al. 1992;Ciucci et al. 2014). A more recent isotopic analysis also noted a relatively lower degree of meat consumption for the Apennine bears compared to other populations (Careddu et al. 2021).
Pyrenean brown bears have been found to have δ 13 C values slightly more positive than ungulates inhabiting the same area, and δ 15 N values more positive than ancient Mont Ventoux brown bears, suggesting a greater consumption of meat (Bocherens et al. 2004). The Pyrenean and Alpine populations both had similar isotopic values, which is likely due to the higher prevalence of animal protein throughout the year. Recent studies indicate that the Pyrenean reintroduced bear population uses the same habitat as the extinct native bears (Palazón 2017), which had already been described as feeding on more animal protein than the Cantabrian bears (Couturier 1954). Still, the results from the Pyrenean bears should be taken with caution, for several reasons: they are only 4 individuals, 2 of them adults from the beginning of the twentieth century. A third individual is a juvenile, also from the beginning of the twentieth century, whose high δ 15 N value seems to indicate that its isotopic signature could be affected by lactation. The fourth bear is current, reintroduced from Slovenia and a larger sample size is required to better interpret the feeding ecology of the current population.

Influence of geography and bioclimatic parameters
For mammalian herbivores, variation in soil δ 15 N, root depth of dietary plant and the composition of woody plants in the diet can all influence the δ 15 N value. Generally, δ 13 C in herbivorous mammals is influenced by the composition of C3 Table 3 Frequencies of the Bayesian estimates for the overlap between the 95% isotopic ellipses for four populations of brown bear (Ursus arctos) expressed as a proportion of the non-overlapping area of the two ellipses (%). Above the ranges with the highest frequency and below the entire range of frequencies . In B the plot shows the positive association between averaged δ15N and ENV PC1 scores plants in their diet. However, the canopy effect also has an influence (Cormie and Schwartz, 1994;Vogel et al. 1990;Ambrose and DeNiro 1986). Other studies have also shown an interaction between δ 13 C and climatic variables, including hours of sunshine, precipitation amount, humidity and temperature (Van Klinken et al. 1994). Our data for the first time allowed to disentangle the geographical signal (that was strong especially on the brown bear hair δ 13 C values) from the climatic one that influence δ 15 N values through temperature and altitudinal parameters summarised by ENV PC1. In previous research, the impact of altitude on isotopic signatures was demonstrated for plants, and subsequently for domestic and wild ungulates, with a negative altitudinal gradient for δ 15 N (Männel et al. 2007;Hofman-Kamińska et al. 2018). The same negative correlation was observed in European Pleistocene cave bears (Krajcarz et al. 2016), whose diet was shown to be primarily plant-based (Bocherens 2019). This is mainly due to the differences in precipitation and temperature. The relationship found between δ 15 N and ENV PC1 is consistent with the observation of a more carnivorous diet in the Pyrenean and some individual Alpine bears. Higher precipitation parameters imply more snowfall that has been demonstrated to impact significantly diet of brown bear at large spatial scale (Bojarska and Selva 2012). It is important to note that the relationship we identified between δ 15 N and altitude + precipitation parameters is not impacted by geographical bias that still occurs in this data as supported by the inclusion of PCNM5 in the best predictive model. This is possibly a limitation related to our sample and such relationship should be better tested on a larger database.
The δ 13 C values in our sample showed a stronger bias in geographical sampling with none of the ENV PCs being selected as best predictors. Generally, for plants, there is a positive relationship between temperature and δ13C (Liu et al. 2017); however, this association has never been tested accounting for geographical bias. Plant δ 13 C are expected to show latitudinal and altitudinal trends (Körner et al. 1991;Kohn 2010), related to climatic and insolation conditions in each area. This is mirrored by some herbivorous species such as white-tailed deer (Cormie and Schwartz, 1994), but this trend does not extend to all herbivores. For instance, there was no systematic association between δ 13 C values in bone collagen and Pleistocene cave bears from various European regions (Krajcarz et al. 2016) or different chronology during MIS 3 .
Our data suggests that a possible strong spatial fidelity in levels of herbivory occurs in the diet of the sampled brown bears; however, this would again require a larger spatial variation. Mowat and Heard (2006) analysed isotopic variation in 81 distinct populations of North American grizzly bears and identified a strong link between level of salmon in the diet and carbon/nitrogen isotopic variation. Still, no spatial filter was included in their diet and  Careddu et al. (2021) on Apennine bears revealed for δ 13 C variation to be linked with individual preferences. Managed bears generally exhibited 0.9‰ higher in δ 13 C than non-managed bears and estimated herb consumption has been detected to change from 15.5 up to 73% within the same season. Consequently, it is possible that the climatic pattern observed in our data reflects proportional differences within the sampled populations of individuals that consume more anthropogenic sources of food or show highly distinct preference in herb consumption.

Bear impact on livestock and crops
Data on livestock predation suggest a relatively higher level of damage for the Alpine population than for the Cantabrian, with livestock damage ratios of 0.26 ± 0.045 for Western Cantabrian population, 0.070 ± 0.043 for the Eastern one and 0.47 ± 0.23 for Trento (Bautista et al. 2017). Across the Pyrenees, there is a strong dissymmetry (Catalonia 0.47 ± 0.23 and France 6.8 ± 1.8) in the damage ratio to livestock, and our Pyrenean bears, with the highest values of δ 15 N being on the Catalonian side. Tosi et al. (2015) found that for bears inhabiting Trentino between 2000 and 2012, their number was correlated with the number of damage events. Considering that we sampled mostly males some of which were directly responsible of livestock damages in 2013 (i.e. M11, M6 and MJ2G1, Groff et al. 2014Groff et al. , 2015, we would expect to find elevated δ 15 N values in these bears, which is not always the case. In summary, the data seem to indicate that although the damage to livestock is apparently significant, it is not systematically reflected in the isotopic values, calling into question the true impact of this type of feeding on the individuals involved. The impact on crops is more difficult to identify from isotopic values, except in the case of maize. Maize is a C4 plant with a distinctive high δ 13 C value. It is cultivated in both the northern Iberian and Italian peninsula, including the highland regions, where, on the other hand, wild C4 plants are very rare (Collins and Jones 1986). Specifically in Apennine bears, no direct or indirect contribution of C4 plants in the diet was detected (Careddu et al. 2021). The δ 13 C values in the four populations studied are consistent with a food chain based on C3 plants. Some alpine specimens show carbon signatures somewhat biased towards positive values, which could have been acquired by feeding on maize or by consumption of maize-fed cattle. However, the δ 13 C values (which do not exceed − 20‰ in any case) are not high enough to point to an appreciable influence of this type of crop on their diet.
Evidently, variation within and between the population sampled can be particularly high in isotopic values. This suggests that livestock/crop consumption can be difficult to detect in brown bears and more controlled experimental data might be needed to elucidate this issue. Our data suggest relatively similar (lower) rate of livestock and crop consumption in Cantabrian and Apennine bears compared to (higher) rate in the sampled populations of central Alps and the Pyreneans.

Conclusion
Isotopic analyses might provide additional insights on the ecological adaptations of brown bear populations. We identified potential overlap for the Cantabrian and Apennine bears while individuals sampled in central Alps and Pyreneans were possibly impacted by altitudinal gradients and greater consumption in animal protein. Our data set the scene for merging databases from different locations across Europe to obtain further information on the variability in brown bear feeding behaviour that might become critical to plan for the future of the species.
at the British Geological Survey were funded by the NERC National Environmental Isotope Facility.
Data availability All data obtained during the development of the work are attached as supplementary information along with the submission of the manuscript.

Declarations
Ethics approval The authors declare that the study followed the institutional and national ethical guidelines for scientific research in the sites where data were collected.
Consent to participate All the authors consent to participate in the development of the investigation.

Consent for publication
All the authors consent for the publication.

Conflict of interest 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/.