2 | METHODS
2.1 | Study system
Our study area encompassed ~13,000 km2in Gävleborg and Dalarna counties in southcentral Sweden with low human density and heavily managed forests of Scots pine (Pinus sylvestris ) and Norway spruce (Picea abies ). Bears were captured via remote drug delivery from a helicopter (Arnemo & Evans, 2017) in spring soon after den emergence (March–May). All capture procedures were conducted in accordance with the Swedish Environmental Protection Agency, Swedish Board of Agriculture, and Swedish Ethical Committee on Animal Research.
2.2 | Sample collection
Bear hair samples were collected from between the shoulders of brown bears during spring captures 1995–2020. After collection, hair samples were placed in individual paper envelopes, labeled accordingly and stored dry at room temperature. To estimate brown bear diet proportions, we also collected brown bear foods within the study area in 2014 and 2015. We collected hair from between the shoulder blades of local, wild moose harvested during the regular hunting season in 2014. Because the moose harvest in Sweden includes all demographic classes of the population, our moose hair sample set includes adults and sub-adults of both sexes. Wild berries (fruits) of the three primary species consumed (bilberry, lingonberry and crowberry) were collected from random locations within the study area in summer 2014. Ants of the generaFormica and Camponotus were also collected within the study area in the summer of 2014. Specimens were mostly adult workers collected from ant hills (Camponotus ) and by sampling coarse woody debris and tree stumps (Formica ) in clearcuts of different age classes at random locations in the study area.
2.3 | Stable isotope analysis
When processing brown bear hair, we separated as much underfur as possible out of the sample, removed large surface contaminants, and weighed hairs to the nearest milligram. Preparatory procedures followed the protocol for cortisol concentration measurement (Macbeth et al., 2010) and samples for stable isotope estimation were portioned out after the grinding step (Sergiel et al., 2017; Appendix S2). We washed each sample three times with 40 μl HPLC grade methanol per mg hair for three minutes per wash to remove other external contaminants (Sergiel et al., 2020). After the hair had dried for at least 24 hours, it was ground to a fine powder in a mixer mill (Retsch MM4000; Retsch GmbH, Germany) at 30 Hz and then put into plastic vials. To measureδ 13C and δ 15N in hair, we followed Koehler at al. (2019) and weighed 1 mg of powdered hair into precombusted tin capsules. Encapsulated hair was combusted at 1030°C in a Carlo Erba NA1500 or Eurovector 3000 elemental analyser. The resulting N2 and CO2 were separated chromatographically and introduced to an Elementar Isoprime or a Nu Instruments Horizon isotope ratio mass spectrometer. We used two reference materials to normalize the results to VPDB and AIR: BWBIII keratin (δ 13C = -20.18,δ 15N = +14.31 per mil, respectively) and PRCgel (δ 13C = -13.64, δ 15N= +5.07 per mil, respectively). Within run (n = 5) precisions as determined from both reference and sample duplicate analyses and from QA/QC controls were ± 0.1 per mil for both δ 13C and δ 15N.
We corrected δ 13C values for the anthropogenic depletion of 13C in the atmosphere (i.e., the Suess Effect) by applying a -0.022‰ correction per year (Chamberlain et al., 2005) and used results from published feeding experiments on ursids (Felicetti et al., 2003; Hilderbrand et al., 1996; Rode et al., 2016) to estimate the isotopic discrimination factors (TDFs) between bear hair and bear serum (Appendix S1.a). We used a similar procedure to estimate TDFs between moose hair and moose meat and offal (Appendix S1.b).
2.4 | Statistical analysis
To answer our hypotheses stated above, we did two separate analyses. First, we used stable isotopes and Bayesian mixing models to estimate annual dietary proportions among three different demographic classes to estimate annual diet over the 25-year study period. Second, we used linear regression models to explain annual variation in dietary proportions and stable isotopes relative to indices of annual food availability.
2.4.1 | Dietary proportion estimation
Based on previous research, we expected brown bear dietary proportions to vary among the demographic classes of independent (no longer dependent on their mother) bears in our population (Steyaert et al., 2013; Swenson et al., 2007).We subset the stable isotope data by demographic classes (females with dependent offspring, solitary females, and solitary males) and ran three separate diet estimation models using year and bearID as random effects. We used previous diet estimates for our population (Mikkelsen et al., 2023) to derive informative priors for our models. We removed six outliers that had particularly highδ 13C values and fell outside the mixing polygon. Each model was run with three chains with 3000000 iterations, a burn-in of 1500000 and a thin rate of 500. We used graphical output as well as fit statistics to determine if each model had run for sufficient time to converge and to ensure proper chain mixing (Semmens et al., 2009; Stock et al., 2018). All analysis was completed in R (R Core Team, 2024) using package MixSIAR (Stock & Semmens, 2013).
Model estimates of the dietary proportions of moose in males had a distinct bimodal distribution, which may arise from the model failing to converge on a single estimate, or from the population having two different diets among males in our sample (i.e., two different possible solutions to the equation). Larger, older males may be more predatory than younger bears (Welch et al ., 1997), thus the bimodal distribution may represent the proportion of moose in the diet for subadult males vs. adult males. To test this, we included an additionala posteriori model for males with an adult and subadult categorical variable as a fixed effect to determine whether this resolved the bimodal distribution.
2.4.2 | Annual variability in food availability
Berries. We used berry inventory information from the Siljansfors Experimental Forest, which is adjacent to the bear monitoring area to estimate the annual productivity of bilberry and lingonberry. Each year, berry production on 54-60 0.25 m2 circular plots was inventoried 2006– 2020. The number of ripe bilberries were counted between July and end of August and the number of ripe lingonberries were counted between end of August and mid-September. Following Hertel et al. (2018), we calculated an annual berry production index as the annual deviation of berry abundance from the 14-year average for each plot. We then created a model to predict berry production with year as a fixed effect and predicted the annual deviation of berries produced. This was scaled between 0 and 1 with indices approaching 0 denoting years of lower-than-average berry production and indices approaching 1 denoting years of higher-than-average berry production (Hertel et al., 2018).
Moose. Moose harvest and observation data was downloaded from Statistik älgdata (https://algdata-apps.lansstyrelsen.se/algdata-apps-stat; Singh et al., 2014) for the counties of Gävleborg and Dalarna. Data in this system are citizen-reported moose observations in the first seven days of the hunting season (October) adjusted by observer effort/observation hours. In Sweden, reporting harvested moose is required by law (Singh et al., 2014). The moose observation database also records the sex, age (calves and adults) and the number of calves with an observed female (singles vs twins), which indicates the overall moose population size, as well as the annual recruitment rate of calves surviving from birth in spring to the fall (Kalén et al., 2022). We used the annual number of total moose observed and harvested after accounting for hunter effort (Singh et al., 2014) as indicator of moose population size. We also used the total number of calves observed as an indicator of annual calf production because calves represent the age class most commonly preyed on by bears in the study area (Swenson et al, 2007).
2.4.3 | Linear regression modelling and model selection
We used mixed-effects linear regression models with the lme4 package (Bates et al., 2014) in R (R Core Team, 2024) to explain variation in diet proportions, δ 13C, andδ 15N values within our population based ona priori hypotheses (Table S1). We used individual ID and demographic class (females with dependent offspring, solitary males, and solitary females) as random effects and used a build-up modeling strategy in which we began by determining the best relationship for each covariate considered (linear, log-linear, or quadratic), and then retained that structure throughout modelling. We used bear age, sex, and annual indices of food availability as fixed effects. All models were compared to the null model to determine whether fixed effects explained more variation than the intercept only, and variables that performed better than the intercept only were used to build more complex models that included additive effects and two-way interactions. We used an information-theoretic criterion for small sample sizes (AICc ) and the relative differences between models (ΔAICc ) when determining the model with the best fit given the data for final inferences (Burnham & Anderson, 2002). For models with similar AICc values, we compared beta estimates, the 95% confidence intervals around the beta estimates, and model variance to select the most parsimonious model (Burnham & Anderson, 2002).
Although we had bear data from 1995–2020, our annual food availability data was limited to 2006–2020, restricting our inference regarding the drivers of variation in δ 13C andδ 15N to this 14-year period. Prior to analysis, we compared the demographic composition, means, medians, and standard errors of our subset data (2006–2020) to the full dataset (1995–2020) to ensure there were no obvious differences in the data (Supplementary Materials S3).
Diet proportions were estimated using the full data set and each demographic class separately, and we used the subset data (2006–2020) and all demographic classes together for linear regression modelling explaining variation in diet proportion and stable isotope values.