Humans and climate change drove the Holocene decline of the brown bear

The current debate about megafaunal extinctions during the Quaternary focuses on the extent to which they were driven by humans, climate change, or both. These two factors may have interacted in a complex and unexpected manner, leaving the exact pathways to prehistoric extinctions unresolved. Here we quantify, with unprecedented detail, the contribution of humans and climate change to the Holocene decline of the largest living terrestrial carnivore, the brown bear (Ursus arctos), on a continental scale. We inform a spatially explicit metapopulation model for the species by combining life-history data and an extensive archaeofaunal record from excavations across Europe with reconstructed climate and land-use data reaching back 12,000 years. The model reveals that, despite the broad climatic niche of the brown bear, increasing winter temperatures contributed substantially to its Holocene decline — both directly by reducing the species’ reproductive rate and indirectly by facilitating human land use. The first local extinctions occurred during the Mid-Holocene warming period, but the rise of the Roman Empire 2,000 years ago marked the onset of large-scale extinctions, followed by increasingly rapid range loss and fragmentation. These findings strongly support the hypothesis that complex interactions between climate and humans may have accelerated megafaunal extinctions.


Supplementary Methods
Evaluation of potential collinearity problems in the meta-analysis. To evaluate potential collinearity problems that may arise from linear relationships between model covariates, we calculated variance inflation factors (VIFs) 1 . VIFs for all covariates were lower than 8.16, less than the threshold of 10 above which collinearity may adversely affect regression results 1 (Supplementary Table 1). To further assess how robust our analysis was with respect to collinearity among predictors, we calculated two alternative variables to characterize the growing season. First, we used growing degree-days (GDD; base 5 °C; data available at http://www.sage.wisc.edu/atlas/) 2 , as a proxy for the length of the growing season, because the length of the growing season is indirectly related to net primary productivity via plant growth 3 . Second, we used the mean monthly temperature of the growing season (T gs ), because, similar to the length of the growing season, the mean temperature of the growing season is indirectly related to net primary productivity via plant growth 3 . We calculated T gs as the mean temperature across the months, in which the minimum monthly temperature was above 5 °C. We used a threshold of 5 °C for the growing season to be consistent with the base temperature for growing degree-days. Using these alternative variables in the analysis reduced VIFs from 8 to 3 but led to identical conclusions, suggesting that collinearity does not affect the conclusions from this analysis (Supplementary Table 1).

Uncertainty in assumptions about Holocene changes in land-use intensification and colonization of Europe by the brown bear from Asia.
To account for uncertainties due to different assumptions about historical changes in land-use intensification (i.e., per capita landuse intensity) 4,5 , we considered two different scenarios in our analysis ( Supplementary Fig. 4).
The first scenario (HYDE 3.1 baseline) effectively omitted land-use intensification by assuming that per capita land-use remained approximately constant over time and that the proportion of land used for agriculture, pastures and urban areas increased linearly with increasing population density. In this scenario, land-use intensity closely resembled the pattern in 1961 CE, the first year for which global Food and Agriculture Organization of the United Nations statistics are available [4][5][6] . The second scenario (HYDE 3.1 concave) assumed a non-linear, concave relationship between population density and per capita land use through time. That is, low-density populations with high per capita land use first expanded to fill all usable land and then intensified land use (i.e., use less land per capita) as population densities increased over time 4,5 . This assumption is similar to the assumption of the KK10 model 7,8 .
The global estimates of land use for the early Holocene from the HYDE 3.1 concave scenario and the KK10 model 4,7 closely resembled each other and provided an upper-bound estimate of land use in the early Holocene, whereas the HYDE 3.1 baseline scenario provided a lowerbound estimate of land use during the same period.
To account for uncertainty regarding colonization of the European continent by the brown bear from Asia 9-11 we considered two scenarios. The first scenario assumed that there was no permanent source population at the eastern border of the study area, while the second scenario allowed for colonization from Asia from a permanent source population at the eastern border. To implement these two scenarios, we fixed the state of the eastern most cells in the occupancy matrix Z s,t in the second scenario to permanently occupied, i.e., we set z s,t = 1, whereas these cells were not informed a priori in the first scenario, so that the occupancy state of these cells was estimated by the model (Supplementary Fig. 5).
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 and two scenarios for colonization of the European continent by the brown bear from Asia (yes versus no), respectively.
Compilation of the archaeofaunal database. Archaeofaunal records were obtained from the Holocene vertebrate database 9,[12][13][14] . In most cases, subfossil bones were context dated from their assignment to an archaeological layer. These dates were often radiocarbon supported and obtained from other materials (e.g., bones or charcoal) from the layer in which subfossil bones were recovered 9,12,13,15 . In a minority of cases (<1% of records) the absolute age of subfossil bones was available from 14 C dating (AMS dating and conventional radiocarbon dating). In the few cases where absolute dates were available these were preferred over relative dates, because relative dates are generally less precise and can sometimes be erroneous due to complicated stratigraphy 14 .
In addition, we applied the following protocol to all records of the archeofaunal database: All archaeofaunal records were mapped in Google Earth (https://earth.google.de/) to validate location information from the database and correct assignment of geographic coordinates (latitude, longitude). Where location information was incorrect, the correct site or location was searched for to obtain accurate geographic coordinate information. Where information for records was missing, the original reference was searched and information filled in. If crucial information on the date or location of the record could not be found, it was rejected.
Implementation of least cost path distances in the metapopulation model. To account for coastline shape and elevation differences (terrain shape), we calculated dispersal distances d between cells as Least Cost Path (LCP) distances. We used a transition function that yields full permeability for flat and downhill (negative) slopes over land, and decreases towards zero for uphill (positive) slope. Water was given extremely low permeability (0.001) to allow rare events of crossing water. The exact shape of the permeability function followed the equation:  To assess the robustness of the relationships, we performed this analysis with three variables describing the growing season (a, net primary productivity; b, growing degree days; and c, temperature during the growing season, T gs ). The model was fitted using a gamma-distribution with an identity-link function for reproductive rate. Given are standardized parameter estimates with their standard error (s.e.m.). For the random effect and the residuals the standard deviations (SD) are given. The R 2 marginal was obtained by comparing the fit of the full model with the fit of a reduced model without the term in question. We assessed whether the analysis was affected by collinearity among predictors, by calculating variance inflation factors (VIFs) and by using alternative variables to characterize the growing season (see Supplementary Methods). Generally VIF-values above a threshold of 10 indicate that collinearity might adversely affect regression results 1 . Using alternative variables to characterize the growing season in the analysis reduced VIFs from 8 to 3 but led to identical conclusions, suggesting that collinearity does not affect the conclusions from this analysis. x   lon and lat, longitude and latitude (decimal degree); LI, inter-birth interval (years); LS, litter size (cubs female −1 ); RR, reproductive rate (cubs female −1 year −1 ); FM, female mass (kg); nLI, nLS, nFM, sample sizes for the three life-history variables; Tgs, temperature of the growing season; Tws, temperature of the winter season; NPP, net primary productivity (kg m −2 a −1 ); GDD, growing degree days (base 5 °C); MAT, mean annual temperature # means parameterization for site types