Meta-analysis

Reproduction imposes one of the most expensive energetic costs on organisms with parental care, such as mammals34. Life-history theory predicts that resource availability and metabolic demands cause a trade-off between life-history traits related to adult survival and reproduction. Low resource availability or high metabolic demands may both force an organism to invest in self-survival, but at the cost of reproduction34. Because life-history strategies are typically under strong environmental selection34, comparative approaches that investigate how investment in reproduction varies across populations in contrasting environments can offer insights into the factors that cause life-history trade-offs.

We used a meta-analysis to test the hypothesis that energy availability during the growing season and energy expenditure during the wintering season impose a life-history trade-off on the reproductive rate of the brown bear (Ursus arctos). To do so, we compiled information on litter size (LS, mean number of cubs per litter), inter-birth interval (LI, mean number of years between litters), and female body mass (FM, kg) for 43 brown bear populations from North America (n = 33), Europe (n = 6) and Asia (n = 4) from 25 published studies (Supplementary Table 3). We excluded five populations where data was not available for all variables from the analysis. Based on the reported litter sizes and inter-birth intervals, we calculated the reproductive rate as RR = LS/LI (cubs per female per year) for each of the remaining 38 populations.

We compiled the minimum, maximum, and average monthly temperatures and average monthly precipitation at the locations of the study populations (spatial resolution of 2° × 2°; data available at http://www.worldclim.org/)35. Brown bears enter and exit the den at a mean temperature of approximately 0 °C36. Therefore, we calculated the mean temperature of the wintering season (T ws ) as the mean temperature across the months, in which the maximum monthly temperature was below 0 °C. Moreover, we estimated net primary productivity (NPP; g dry matter m−2 a−1), based on the minimum of the temperature- and precipitation-limited rates of NPP using the Miami model37. The two limiting rates were calculated independently as simple empirical functions of mean annual temperature (T) and total annual precipitation (P), following the equations given in Adams et al.37:

$$\begin{array}{c}NPP=\,{\rm{\min }}(NP{P}_{T},\,NP{P}_{P}),\,{\rm{with}}\,\\ NP{P}_{T}=3000{(1+\exp (1.315-0.119T))}^{-1}\,{\rm{and}}\\ NP{P}_{P}=3000(1-\exp (-0.000664P)).\end{array}$$ (1)

The model does not treat the effects of seasonality, CO 2 , humidity, light capture, or vegetation and soil types explicitly37. However, we chose this method, because it allowed for a direct comparison with the results of the metapopulation model. Given the lack of detailed data about NPP during the Holocene, the Miami model provides a good first-order characterization of the major climatic controls on ecosystem productivity during the past 12,000 years.

We fitted a generalized linear mixed effects model to analyse variation in reproductive rate using a gamma distribution with an identity-link function. In this model, we included a random effect for continent (North America, Europe and Asia) to account for differences in the phylogeographic histories of brown bear populations across continents38,39,40. Moreover, we included female body mass as a fixed explanatory variable to account for local adaptations to differences in environmental conditions among the brown bear populations41, 42. We included NPP and T ws into the model to test the hypothesis that energy availability during the growing season and energy expenditure during the wintering season constrain the reproductive rate of the brown bear. As the sample size varied among the study populations, we performed the analysis with and without weighting of the observations by their sample size. We calculated the weight of observations as w = log 10 (N FM ) + log 10 (N LI ) + log 10 (N LS ), where N FM , N LI and N LS were the sample sizes for female body mass, inter-birth interval and litter size, respectively. In the case of observations with missing information about sample size, we assumed a sample size of one. Weighted and non-weighted analyses led to identical conclusions. In the main text we report results of the weighted analysis, because it incorporates variability in uncertainty in the estimates of life-history traits among populations. We assessed whether the analysis was affected by collinearity among predictors, by calculating variance inflation factors and by using alternative variables to characterize the growing season and found the results of our analysis to be robust (Supplementary Methods and Supplementary Table 1).

Metapopulation model

To estimate the spatiotemporal population dynamics of the European brown bear during the Holocene, we fitted a Bayesian hierarchical metapopulation model24, 43 to detection/non-detection data based on archaeofaunal and recent occurrence records. Below we give detailed information about the spatial and temporal resolution of the analysis, the data used, and the structure of the model.

Spatial and temporal resolution of analysis

For our analysis, we defined seven periods: (1) 10000–7000 BCE, (2) 7000–3000 BCE, (3) 3000 BCE – 0 CE, (4) 0–1500 CE, (5) 1500–1800 CE, (6) 1950–1970 CE, and (7) 2010–2015 CE. We used the midpoint of each of the seven periods to determine the length of the time intervals Δt between the respective periods in the metapopulation model (e.g., Δt 12 = 8500–5000 = 3500 years). The time intervals Δt between periods were 3500, 3500, 2250, 900, 310, and 52.5 years, respectively. Moreover, we used a 100 × 100 km2 grid (100-km grid, hereafter) including 789 terrestrial cells across Europe for the analysis (Fig. 2c). By using a coarse grain of 100 km in our analysis, we minimized the problem of missing data in the archeofaunal occurrence records and avoided inflating data quality due to inappropriate downscaling of reconstructed past land use and climate data44.

Occurrence data

For the first five periods (covering 10000 BCE – 1800 CE) we compiled spatially explicit detection histories of the brown bear across Europe using a database containing published and unpublished archeofaunal records from dated archaeological and paleontological excavation sites19, 25,26,27 (Supplementary Methods and Supplementary Fig. 2). The excavation sites included settlements (n = 2,735), castles (n = 371), caves (n = 274), burial sites (n = 264), middens and rubbish pits (n = 174), cult and religious sites (n = 76), moors and riverbeds (n = 11), as well as 272 sites without known designation. For each record, we extracted data on the presence and abundance of subfossil bone remains of brown bear, red deer (Cervus elaphus), and roe deer (Capreolus capreolus) (total abundance, median: 7, range: 1–26,428). Records containing red or roe deer, but not brown bear, served as negative controls, because both deer species were typical elements of the European fauna during the Holocene and the most frequent prey of European Stone Age hunters25, 26. The rationale behind this approach is that if remains of red or roe deer were present at a given site, it was suitable for preservation of bone material and, thus, for detection of brown bears, if the species had been present. This allowed us to account for imperfect detection in the archaeofaunal record, while simultaneously modelling the spatiotemporal variation in brown bear extinction and colonization rates across the European continent. In total, our dataset comprised 4,177 archeofaunal records from 3,461 excavation sites (records from archaeological layers with distinct temporal context were treated separately). We used the above records of brown bear presence and absence in each cell and period to obtain detection/non-detection data (i.e., detection histories) for the first five periods (i.e., 10000 BCE – 1800 CE). To inform the model about the distribution of the brown bear in the last two periods, we used recent occurrence data from Chapron et al.20 and Chestin45 (Supplementary Fig. 2). These data were resampled to the 100-km grid.

Land-use data

We compiled land-use data from the HYDE 3.1 database46. These data represent spatially explicit reconstructions of fractional land use covering the period 10000 BCE to 2000 CE. The data have a spatial resolution of 5′ × 5′ and a temporal resolution of 1,000-year time intervals for the BCE period, then 100-year time intervals until 1700 CE, and 10-year time intervals until 2000 CE. For analysis, we extracted three variables reflecting human impact on the environment, namely the proportion of a given cell covered by agriculture, pastures, and urban areas. We used the sum of these variables as a proxy for human land use, that is, the total fraction of land appropriated by humans. Then we averaged the fractional land use for each cell across the time intervals that fell in each of the seven periods used for the metapopulation model. The time intervals that occurred at the margins of the periods received only half the weight in the calculations compared to those time intervals in the core of the periods. We resampled the land-use data to the 100-km grid for analysis. To account for uncertainties due to different assumptions about historical changes in land-use intensification (i.e., per capita land-use intensity)47, 48, we considered two different scenarios of the HYDE 3.1 database in our analysis (Supplementary Methods and Supplementary Fig. 4).

We decided to use land use as a proxy for human impact, because unlike other measures, such as human population density, land use integrates the various effects of humans on brown bear populations. First, previous work has shown that habitat loss due to land use severely constrains the occurrence and breeding range of the European brown bear49, 50. Second, the potential for human-bear conflicts due to predation on livestock or wild ungulates increases in intensively used landscapes51. Therefore, land-use intensification historically coincided with the persecution and local extirpation of the brown bear and other large carnivores in Europe21, 52. This is why persecution and habitat loss are generally the most important drivers of extinction risk in mammals, and particularly in large carnivores22, 53. While the behaviour of humans in relation to bears might well have changed through time, there is evidence that humans have hunted bears since at least 9000 BP54 and that bear hunting and the associated rituals and beliefs were deeply rooted in the prehistoric culture of humans in the northern hemisphere55. Thus, we expect land use to be a reliable proxy of the combined effects of hunting, persecution and human-induced habitat loss.

Climate, net primary productivity and elevation data

We compiled climate data from Luterbacher et al.56, Pauling et al.57, and Mauri et al.30. The data from Luterbacher et al.56 and Pauling et al.57 represent spatially explicit seasonal temperature and precipitation reconstructions covering the period 1500–2000 CE. The data have a spatial resolution of 0.5° × 0.5° and a temporal resolution of 1-year. The data from Mauri et al.30 represent spatially explicit pollen-based climate reconstructions covering the period 12000 BP (before present) to 100 BP. The reconstructions are given as anomalies relative to late preindustrial time, centred approximately around 100 BP (i.e., 1850 CE), with a dating uncertainty of ± 250 years30. The data have a spatial resolution of 1° × 1° and a temporal resolution of 1,000-year time intervals until 1000 BP and a single 900-year time interval until 100 BP, respectively. From these datasets, we extracted mean annual temperature, total annual precipitation (to estimate NPP see below) and mean winter temperature for analysis.

We first calculated 31-year running averages for each of the three variables from Luterbacher et al.56 and Pauling et al.57 for each cell for the period 1500 to 2000 CE. Then we used the 31-year running averages for the year 1850 CE as the baseline to calibrate the anomalies reported in Mauri et al.30. We estimated net primary productivity (g dry matter m−2 a−1) from mean annual temperature and total annual precipitation using the Miami model37, 58, 59 in equation (1). For analysis, we used the same time intervals and weighting scheme as for the land-use data to calculate mean net primary productivity and mean winter temperature for each cell of the 100-km grid in each of the seven periods.

We extracted high-resolution elevation data from WorldClim35 using an initial resolution of 0.5′ × 0.5′ (~1 km2 at the equator). For analysis, we resampled the elevation data to the 100-km grid.

Model structure

To decompose the contribution of humans and climate change to the Holocene decline of the European brown bear, we combined a spatially explicit metapopulation model with a path analysis in a Bayesian hierarchical framework24, 60. In particular, we based our analysis on a hierarchical state-space formulation of a spatially explicit metapopulation model that accounts for imperfect detection of brown bear occurrence in the archeofaunal record24. To describe detection/non-detection data y s,t in grid cell s in period t, the hierarchical state-space model is based on two linked stochastic processes24: a sub-model for the observations conditional on the unobserved state process, i.e., y s,t | z s,t (observation model) and a sub-model for the unobserved or partially observed state process z s,t (state model).

State model

The state model has a simple formulation in terms of initial occurrence probability, i.e., at t = 1, which we will designate μ z,s,t = 1 . The initial occurrence states are assumed to be independent identically distributed Bernoulli random variables, denoted as z s,t = 1 ~ Bernoulli(μ z,s,t = 1 ). In all subsequent periods, the occurrence state z s,t + 1 is a Bernoulli trial with a success parameter μ z,s,t + 1 that depends on whether cell s was occupied in period t and on the value of two parameters that indicate either the probability of extinction, p φ , or colonization, p γ , so that:

$$\begin{array}{c}{\mu }_{z,s,t+1}={z}_{s,t}(1-{p}_{\phi ,s,t})+(1-{z}_{s,t}){p}_{\gamma ,s,t}\,{\rm{with}}\\ {z}_{s,t+1}|{z}_{s,t} \sim {\rm{Bernoulli}}({\mu }_{z,s,t+1}).\end{array}$$ (2)

Because the lengths of the periods Δt in our dataset were unequal (i.e., they ranged from 3500 to 52.5 years), the assumptions of the discrete-time Markov process were violated61, preventing direct estimation of transition probabilities from our data. Therefore, we used Kolmogorov’s forward equations61 (equivalent to matrix exponentiation) to map the colonization (r γ ) and extinction (r φ ) rates for each time interval Δt on the colonization (p γ ) and extinction (p φ ) probabilities, respectively:

$${p}_{\gamma ,s,t}=1-({r}_{\phi ,s,t}+{r}_{\gamma ,s,t}\exp (-({r}_{\phi ,s,t}+{r}_{\gamma ,s,t}){\rm{\Delta }}t))/({r}_{\phi ,s,t}+{r}_{\gamma ,s,t}),$$ (3)

$${p}_{\phi ,s,t}=1-({r}_{\gamma ,s,t}+{r}_{\phi ,s,t}\exp (-({r}_{\phi ,s,t}+{r}_{\gamma ,s,t}){\rm{\Delta }}t))/({r}_{\phi ,s,t}+{r}_{\gamma ,s,t}),$$ (4)

We modelled the extinction rate (r φ,s,t ) as a function of land use (LU), winter temperature (WT), net primary productivity (NPP), and elevation (ELE) with a logarithmic-link function:

$$\mathrm{log}({r}_{\phi ,s,t})={\beta }_{\phi ,0}+{\beta }_{\phi ,LU}L{U}_{s,t}+{\beta }_{\phi ,NPP}NP{P}_{s,t}+{\beta }_{\phi ,WT}W{T}_{s,t}+{\beta }_{\phi ,ELE}EL{E}_{s,t},$$ (5)

We modelled the colonization rate r γ,s,t of the species as a function of spatial connectivity to extant cells with a logarithmic-link function. Connectivity describes the distance-dependent influence of all potential neighbouring source populations s' via a negative exponential dispersal kernel62:

$$\mathrm{log}({r}_{\gamma ,s,t})=\lambda \sum _{s

e {s}^{^{\prime} }}\exp (-\alpha {d}_{s,{s}^{^{\prime} }}){z}_{{s}^{^{\prime} },t}$$ (6)

where z s′,t is the occupancy of cell s′ in period t, d s,s′ is the distance between cells s and s′, 1⁄α is the average dispersal distance, and λ is a constant. We did not estimate the additional scaling constant λ from our data62, because direct estimation of this parameter from our data would cause parameter redundancy with the intercept β φ,0 in equation (5). To avoid problems with parameter identifiability, we estimated β φ,0 and set λ = 143. Using alternative values of λ did not affect any of our conclusions, because it causes a proportional shift of the estimate for the intercept β φ,0 . To account for coastline shape and elevation differences (terrain shape), we calculated dispersal distances d between cells as Least Cost Path (LCP) distances (Supplementary Methods). We extracted high-resolution elevation data from WorldClim35 (0.5′ × 0.5′ resolution) and resampled these data to a resolution of 5 × 5 km2 for calculation of Least Cost Path (LCP) distances. To account for uncertainty regarding colonization of the European continent by the brown bear from Asia19, 38, 63 we considered two scenarios with and without colonization of Europe from Asia, respectively (Supplementary Methods and Supplementary Fig. 5).

Observation model

We linked the latent occurrence state z s,t of the state model through the observation model to each occurrence record i associated with cell s in period t, y i,s[i],t[i] . Therefore, we specified the observation model conditional on the latent process z s,t as:

$${y}_{i,s[i],t[i]}|{z}_{s,t} \sim {\rm{Bernoulli}}({z}_{s,t}{p}_{i,s[i],t[i]}),$$ (7)

Thus, if cell s is occupied in period t, the data are Bernoulli trials with a parameter for the detection probability, p i,s[i],t[i] . If cell s is unoccupied in period t, then the data are Bernoulli trials with P(y s,t = 1) = 0. For the first five periods that were informed by archeofaunal records, we modelled the detection probability p i,s[i],t[i] associated with the ith record as a function of the total number of bone remains (NSP) of the three species in a given excavation record and the type of excavation site (TYP) from which the records originated using a logit-link function as:

$${\rm{logit}}({p}_{i,s[i],t[i]})={\beta }_{p,0}+{\beta }_{p,NSP}\,\mathrm{log}(NS{P}_{i})+{\beta }_{p,TYP}TY{P}_{i},$$ (8)

The number of bone remains reflects both sampling effort and the potential of an excavation site to preserve bone material4, 28 and also may reflect the intensity of hunting activities and thus the hunting pressure on the regional fauna64. The type of excavation site may reflect cultural habits that could affect the prevalence of brown bear remains in certain types of excavation sites (e.g., burial sites associated with rituals)65. For the 272 of the 4,177 excavation records for which we did not have information about the designation of the site, we randomly assigned a site type TYP i based on the relative occurrence of the seven site types in the different periods following a multinomial distribution: TYP i,s[i],t[i] ~ Multinomial(P n,t[i] ), where P is a matrix with n rows indicating the seven site types and t columns indicating the five periods for which we had archeofaunal records, and each cell contains the frequency of occurrence of each site type in each period. Therefore, the uncertainty originating from these records was fully propagated through all steps of model estimation. For the last two periods (i.e., 1950–70 s and recent) that were informed by species distribution maps, we assumed perfect detection, so that y i,s[i],t[i] = z s,t . We are confident that this is a reasonable assumption, given the coarse grain of our analysis (i.e., 100 × 100 km) and the huge effort spent in the monitoring of European brown bear populations since the 1950s20.

Path analysis

Besides this state-space model describing the direct effects of the explanatory variables on the extinction rate, the path model contained three regression equations to describe the potential causal relationships among the predictor variables. To account for spatial autocorrelation, these regressions included a random effect for grid cell. We modelled the effect of elevation on winter temperature as:

$$\begin{array}{c}{\mu }_{WT,s,t}={\beta }_{WT,0}+{\beta }_{WT,ELE}EL{E}_{s,t}+cel{l}_{WT,s},\,{\rm{with}}\\ W{T}_{s,t} \sim {\rm{Normal}}({\mu }_{WT,s,t},{\tau }_{WT}),{\rm{and}}\,cel{l}_{WT,s} \sim {\rm{Normal}}(0,{\tau }_{WT,cell}).\end{array}$$ (9)

We modelled the effect of winter temperature and elevation on net primary productivity as:

$$\begin{array}{c}{\mu }_{NPP,s,t}={\beta }_{NPP,0}+{\beta }_{NPP,WT}W{T}_{s,t}+{\beta }_{NPP,ELE}EL{E}_{s,t}+cel{l}_{NPP,s},\,{\rm{with}}\\ NP{P}_{s,t} \sim {\rm{Normal}}({\mu }_{NPP,s,t},{\tau }_{NPP})\,{\rm{and}}\,cel{l}_{NPP,s} \sim {\rm{Normal}}(0,{\tau }_{NPP,cell}).\end{array}$$ (10)

Finally, we modelled the effect of the aforementioned variables on human land use as:

$$\begin{array}{c}{\mu }_{LU,s,t}={\beta }_{LU,0}+{\beta }_{LU,NPP}NP{P}_{s,t}+{\beta }_{LU,WT}W{T}_{s,t}+{\beta }_{LU,ELE}EL{E}_{s,t}+cel{l}_{LU,s},\,{\rm{with}}\\ L{U}_{s,t} \sim {\rm{Normal}}({\mu }_{LU,s,t},{\tau }_{LU})\,{\rm{and}}\,cel{l}_{LU,s} \sim {\rm{Normal}}(0,{\tau }_{LU,cell}).\end{array}$$ (11)

Model selection

To select informative paths, we used Stochastic Search Variable Selection (SSVS) with global adaptation66. The variable selection procedure can be seen as one of deciding which of the regression parameters β j in a given model are equal to zero. To denote whether the variable j is present in the model, we used an indicator variable I j (where I j = 1 indicates presence, and I j = 0 absence of variable j in the model). The prior inclusion probability of all explanatory variables was set to P = 0.5. The SSVS approach uses a mixture prior for each regression parameter β j following the equation:

$$P({\beta }_{j}|{I}_{j})=(1-{I}_{j})\cdot {\rm{Normal}}(0,{\tau }_{\beta })+{I}_{j}\cdot {\rm{Normal}}(0,g{\tau }_{\beta }),$$ (12)

where the first density (the spike) is centred around zero and has a small variance66. Here we use a random effects variant of SSVS66, in which in the spike part of the prior τ β is fixed to a constant (τ β = 3,600), and in the slab part of the prior the product gτ β is estimated by the model. The variance σ 2 β = gτ β −0.5 was drawn from a uniform prior between 0 and 20. This form of global adaptation has the advantage of facilitating the tuning of the variable selection, because the distribution of each coefficient is shrunk towards the correct region of the parameter space by the other coefficients in the model. In order to infer which of the variables should be included in the models, we used Bayes factors66.

Model implementation

The models were fitted using an Markov Chain Monte Carlo (MCMC) algorithm with a Gibbs sampler in JAGS 67 called from R 68 by using package rjags 69 (for the JAGS code see Supplementary Computer Code 1). We used uninformative prior distributions for most parameters. Running the models with a uniform prior for the dispersal parameter α suggested that the information contained in our dataset did not suffice to inform this parameter, because the estimate for this parameter was poorly defined. Therefore, we specified an informative prior following a gamma distribution (shape = 3, rate = 1) with highest mass around 2 for the dispersal parameter α, which corresponds to an average dispersal distance of 1⁄α = 50 km (given the 100-km grid). This value roughly corresponds to the average natal dispersal distance for Scandinavian brown bears (males: 108.3 km ± 27.4 km; females: 15.7 km ± 2.4 km)70. The decision of whether to use an informative or uninformative prior for the dispersal parameter did not affect any of the results or conclusions.

We ran 10 parallel chains for each model. The initial values for the chains were drawn randomly from uniform distributions. The models were run for 4,500 iterations with an adaptive burn-in phase of 2,000 iterations and a thinning interval of 25 iterations, resulting in 100 samples per chain, corresponding to 1,000 samples from the posterior distribution. The chains were checked for convergence, temporal autocorrelation, and effective sample size using the R package coda 71. In the main text, we report the parameter estimates after pooling the posterior samples from the two-factorial design with two scenarios for changes in per capita land-use intensity during the Holocene (constant versus decreasing) and two scenarios for colonization of the European continent by the brown bear from Asia (yes versus no), respectively (Supplementary Methods and Supplementary Figs 4 and 5). Therefore, these estimates account for potential biases due to the uncertainty in these model assumptions. The models were run on the high-performance computing cluster of the PL-Grid Infrastructure.

Data availability

The datasets generated and analysed during the current study are available from the corresponding author upon reasonable request.