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.