Large carnivore expansion in Europe is associated with human population density and land cover changes

The recent recovery of large carnivores in Europe has been explained as resulting from a decrease in human persecution driven by widespread rural land abandonment, paralleled by forest cover increase and the consequent increase in availability of shelter and prey. We investigated whether land cover and human population density changes are related to the relative probability of occurrence of three European large carnivores: the grey wolf (Canis lupus), the Eurasian lynx (Lynx lynx) and the brown bear (Ursus arctos).


| INTRODUC TI ON
During the nineteenth and mid-twentieth centuries, large carnivores in Europe saw their numbers and distribution decline sharply, mainly due to human persecution, habitat loss and fragmentation (Chapron et al., 2014;Linnell et al., 2009). However, over the last few decades, most populations of Eurasian lynx (Lynx lynx), brown bear (Ursus arctos) and grey wolf (Canis lupus) have increased in Europe, returning to areas where they were previously locally extinct. In recent times, there have been records of wolf presence in all EU member states (excluding islands), and nearly all mainland states have at least one permanent and reproducing species of large carnivore (Chapron et al., 2014;Fernandéz-Gil et al., 2018).
In Europe, large carnivore species recolonized former range starting from source populations living in remote and mountainous regions (Linnell et al., 2010;Swenson et al., 1995), and also thanks to specific reintroduction and translocation programs for lynx and bears (Breitenmoser et al., 1998;Samojlik et al., 2018;Zedrosser et al., 2001). Wolves expanded their permanent range in the Italian and Iberian peninsulas (except for Southern Spain and Portugal; López-Bao et al., 2018) and recolonized France, Switzerland and Scandinavia Hindrikson et al., 2017).
Wolves that re-established within their previous distribution range in Eastern Germany and Western Poland originated from the Baltic population in northeastern Poland (Czarnomska et al., 2013;Nowak & Mysłajek, 2016). Those occupying central Spain dispersed from the Cantabrian Mountains (Silva et al., 2018), and in Sweden, the wolf population is descended from a few individuals that came from the Finnish-Russian border region (Linnell et al., 2005;Wabakken et al., 2001). Reintroduction and augmentation programmes have been developed to spur and support the recovery of lynx and bear populations in Central Europe (Breitenmoser et al., 1998;Samojlik 13 Callisto -Wildlife and Nature Conservation Society, Thessaloniki, Greece et al., 2018;Zedrosser et al., 2001), whereas wolf recolonization was entirely natural. Natural expansions of remnant populations were recorded in Fennoscandia for the lynx , and in the Baltics, Fennoscandia and the Dinaric mountains for the bear, even though other lynx populations in the rest of Europe did not increase in number nor expanded their range (e.g. Heurich et al., 2018).
Despite the different historical patterns in the expansion and re-establishment of wolves, lynx and brown bears, some common factors have been proposed that positively affect large carnivore recovery in Europe, including ecological, societal and cultural elements. After the end of World War II, Europe went through a series of major land use changes, with agricultural and pastoral land abandonment, decrease in natural grassland, increase in forest cover both due to afforestation and reforestation programmes and increased urbanization (Behr et al., 2017;Fuchs et al., 2013), with subsequent concentration of the human population in urban areas. Such changes are thought to be linked to increases in wild ungulate populations (Burbaite & Csányi, 2009;Burbaitė & Csányi, 2010) and increased availability of denning and refuge areas, as well as a decrease in human persecution Chapron et al., 2014). Socio-economic and land use trends further changed during the last thirty years. On average, cropland, grassland and forest areas remained stable or decreased slightly, while urbanization increased even further at the expense of cropland (Li et al., 2018). In parallel, legislation changes led by the Council of Europe (1979) and the EU Habitats Directive (1992) provided a basis for halting the persecution that had occurred in the previous decades, therefore reversing the trend towards population decline by introducing a fully protected status, strict regulations on hunting, and the ban on poison baits (Boitani, 1992;Díaz, 2010;Evans, 2012). Furthermore, increased protected area coverage may have also contributed to large carnivore recolonization, together with a shift in the societal perception towards these species (Dressel et al., 2015), although conflict and illegal persecution still persist in some areas (Liberg et al., 2012). Importantly, the aforementioned changes considerably varied between different parts of Europe, partly reflecting the fall of the Iron Curtain and its socio-economic consequences. In Western Europe, forest area remained largely stable with the exception of a mild decline in parts of the Alps, Apennines, Pyrenees and Southern Scandinavia, while an increase was observed in Eastern Europe and North Scandinavia (Nowosad et al., 2019). Conversely, the human population in Europe increased from around 480.5 million people in 1992 to approximatively 509.7 million in 2015 (The World Bank, 2019). However, while population density increased in nearly all Western and Southwestern regions of the EU, it decreased in most of Northeastern, Eastern and part of Southeastern Europe.
The relative effects of these changes in land cover structure on large carnivore distributions at the European scale remain unclear.
Unfortunately, ground data are scarce at such a large scale. For example, a complete dataset of spatio-temporal trends in ungulate distribution and abundance at the European scale is not available as a measure of prey abundance. Similarly, actual on-ground law enforcement and societal changes in the perception of large carnivores have been assessed at the local scale but are difficult to quantify or estimate at large spatial scales (van Heel et al., 2017;Røskaft et al., 2007).
However, high-resolution land cover time series between 1992 and 2015 together with spatially explicit estimates of human population density, which can be used as a proxy of human pressure on predators and prey populations (Basille et al., 2009;Milanesi et al., 2017;Woodroffe, 2000), are now available can be used to predict habitat suitability for large carnivores at the European scale.
In this study, we investigated if land cover, human population density and protection status are related to the occurrence of Eurasian lynx, brown bear or grey wolf across Europe and over the past 24 years. We compiled a dataset of >50,000 occurrence points and matched annual land cover data as well as spatially explicit human population density information covering the period between 1992 and 2015. We fitted multi-temporal species distribution models (SDMs)-also known as time-calibrated models-then predicted how the relative probability of occurrence of the species (i.e. henceforth defined as habitat suitability) changed through time by projecting our models at different time steps within the temporal frame considered. Unlike classic SDMs, this approach matches occurrence data collected at a given time with the respective temporal covariate (Ancillotto et al., 2016;Maiorano et al., 2013;Nogués-Bravo, 2009), thus reducing the niche underestimation due to limited or spatially biased observations, and preventing model extrapolation beyond the environmental space of the training dataset (Maiorano et al., 2013;Nogués-Bravo, 2009). Finally, we projected the model through time to assess which factors (land cover, landscape configuration, human population density and protection status) can better underlie the temporal change in the species presence over the last 24 years, thereby highlighting the regional and species-specific differences in the possible drivers of change. We expected an increase in suitability in recently recolonized areas in response to a decrease in the spatial cohesion of cropland (i.e. intensity), an increase in forest cover and a decrease in human population density in rural areas.
We also expected to find a positive effect of conservation status and hunting legislation but little or no effect of protected areas as shown in previous studies Tammeleht et al., 2020).

| Species distribution data
We integrated occurrence data for the three species from both published (see Appendix S1) and unpublished datasets provided by several coauthors (HA, CB, FC, DC, PC, MK, MM, YM, LP, NS, JLB, IR, IT, AZ, TZ-K). We filtered the data to comprise only those that unambiguously pertain to large carnivores, including direct observations, snow tracking, and Very High Frequency (VHF) or Global Positioning System (GPS) relocations, and also locations where animals were non-invasively detected (i.e. non-invasive genetic sampling, camera trapping) or where they were found dead. When possible, we discarded occurrence points clearly associated with dispersing animals possibly sampled in areas of sporadic occurrence (based on expert judgement) which could potentially have biased the model. We only retained occurrence points collected between 1992 and 2015 with a spatial precision ≤10 km (matching the environmental data timeframe and resolution, see below) and within European countries west of 64° longitude (including non-EU member states). From these, we only retained one occurrence point per 10 km grid cell per year to avoid including pseudo-replicates, resulting in a final dataset of 21,852 presence points for wolf, 17,943 for lynx and 5,224 for bear. The final distribution of occurrence records corresponds well to the European atlas maps for the three species presented in Chapron et al. (2014) (Figure 1, see Figures S1-S3 for distribution of data points over time).

| Land cover and human population density
We focused our study on dynamic covariates for which multiple temporal layers were available (i.e. land cover, human population density and protected areas designation) and included topography as the only static variable in order to control for terrain roughness, which may limit carnivore dispersal and influence refuge site selection.
We identified a total of 15 variables (Table 1), which can be divided into four main classes: land cover, landscape configuration, anthropogenic pressure and topography. We used a resolution of 10 km with the projection ETRS89-Lambert Equal Area, matching the grid used to produce the European atlas of large carnivores (Chapron et al., 2014). The same resolution has been used in previous distribution modelling studies on the same species in Europe (Chapron et al., 2014;Milanesi et al., 2017) and is considered to be a good compromise between information precision and the average home range of these species in Europe that can vary between 100 km 2 to more than 1,000 km 2 , depending on the species and region (Dahle & Swenson, 2003;Jedrzejewski et al., 1996Jedrzejewski et al., , 2007Mancinelli et al., 2018).
We retrieved land cover data from the ESA CCI Land Cover (CCI-LC) project (ESA, 2017; downloaded from https://www.esalandc over-cci.org/), which provides maps at 300 m spatial resolution from 1992 to 2015. We aggregated the 37 different classes of land cover into 9 macro-classes considered to be relevant for the ecology of the species (see Table S3 for reclassification rules) and resampled each layer at 10 km as the percentage of land cover per grid cell.
Because landscape composition and configuration may influence large carnivores' spatial ecology and predator-prey dynamics, thus influencing predator habitat selection (Gorini et al., 2012), we also considered the diversity, interspersion and cohesion of land cover types (McGarigal et al., 2012). We calculated the Shannon diversity index to describe landscape heterogeneity based on the different proportions of land cover type per grid cell. We calculated the Contagion index (Li & Reynolds, 1993)  We considered three proxies of anthropogenic pressure: distance from urban areas, human population density and its coefficient of variation within grid cells as a measure of heterogeneity. We measured the Euclidean distance from urban agglomeration using the urban layer from ESA CCI Land Cover (CCI-LC) maps. We obtained the human population density data at 30 arc-second horizontal res- To measure the heterogeneity in human population density, we calculated the coefficient of variation of human population density values within the 10 km cells. The effect of heterogeneity of human population density may be conditional to the level of population density. For example, cells with low population density and high heterogeneity may be less suitable than areas with low population density and low heterogeneity, but the opposite may be true for high population density cells. Consequently, we have also included the interaction between the human population density and its coefficient of variation in the model. As for topography, we generated slope and a terrain ruggedness index (TRI) layers using the EarthEnv-DEM90 digital elevation model available at https://www.earth env. org/DEM (Robinson et al., 2014). Because there have been substantial changes in the protection status of all three species and in the coverage of protected areas during the time period covered, we also included three covariates: law implementation, hunting regulations and protected area coverage. We retrieved the year of protection per species per country based on a literature search (see Table S1).
For EU countries, we considered the year 1992 as baseline, unless a country entered the EU later (e.g. Slovenia) or signed the Habitat Directive subsequently (e.g. Romania). If a species was included in Annex V of the Habitat Directive (Appendix S2, Table S1), it was considered protected but hunted. When the only information available concerned hunting regulation, we considered the species not protected and hunted. Implementation of the law and hunting variables were included in the model as binary variables. To consider the role of protected areas on the conservation of large carnivores, we calculated the percentage of protected area coverage per 10 km grid cell per year. We used the dataset WDPA (IUCN & UNEP-WCMC, 2020) and overlaid the 10 km grid to its shapefile, considering all types of protected areas. We resampled all environmental variables to a resolution of 10 × 10 km. Changes to each predictive variable over the 24-year period is presented in Figure S4.
To meet normality assumptions, all the predictors were inspected and, if necessary, transformed with either a logit or log 10 -transformation (Table 1) and then standardized (i.e. mean = 0; standard deviation = 1). For each model, we then discarded collinear variables having variance inflation factors (VIF) >3 and Pearson correlation > 0.7 (Table S4, Figures S5-S7) (Dormann et al., 2013;Zuur et al., 2010), retaining variables deemed more biologically relevant for the species. Specifically, we gave priority to forest cover over cropland cover, to human population density over Euclidean distance from urban areas (whose effect can change over shorter distances than 10 km resolution) and to Contagion index over the  Table  S3 for reclassification rules and concise description of the land cover variables) Shannon index. However, since cropland expansion is an important driver of biodiversity loss (Kehoe et al., 2017) and may be associated with the abundance of large carnivore prey (Bunnefeld et al., 2006), we also repeated the analyses with cropland instead of forest cover (see Appendices S1-S2). To further assess whether low correlation values could influence our conclusions, we ran a simulation analysis to estimate the potential biases in model coefficients. This analysis indicated that potential errors in parameter estimates derived of multicollinearity were negligible (<0.01) (Appendix S3).

| Modelling
We generated pseudo-absences by sampling points from NUTS2 regions (Nomenclature of Territorial Units for Statistics, https:// ec.europa.eu/euros tat/web/nuts/backg round, Figure S8) where presence data have been recorded. Because NUTS2 only includes EU member countries and associated countries (e.g. Norway and Switzerland), we integrated the missing countries (Belarus, Bosnia and Herzegovina, Kosovo, Moldavia, Russia and Ukraine) using the regional map from https://gadm.org ( Figure S9). To compensate for the geographic and temporal bias in the presence data, we sampled as many pseudo-absences as the number of occurrence points present in the same NUTS2 region per year. As we had <10,000 occurrence points for brown bears, we sampled 10,000 background points as suggested in Barbet-Massin et al. (2012) proportionally to the presence points in the same NUTS2 regions.
We built the multi-temporal SDMs using generalized linear mixed models (GLMM) with a LASSO (Least Absolute Shrinkage and Selection Operator) penalizing algorithm (GLMMLasso, Schelldorfer et al., 2014).
Since the average habitat suitability for the species in different regions in Europe can differ due to factors beyond the predictors considered, such as different local management practices, diverse perceptions of these species by local communities and possible context-specific responses to the different predictors across Europe, we used the NUTS2 regions as random effects in order to allow the intercept of the model to vary across these regions. Unlike a standard GLMM, the GLMMLasso does not require a stepwise regression for model selection, which presents several drawbacks (Whittingham et al., 2006).
LASSO is an alternative to the multiple stepwise selection which enhances predictive accuracy and model interpretability (Tibshirani, 1996). The L-1 penalty shrinks regression coefficients towards zero in the attempt to maximize predictive ability, while allowing for model regularization and variable selection, thus avoiding overfitting (Groll & Tutz, 2014;Schelldorfer et al., 2014). The amount of shrinkage is de-

| Spatio-temporal trend in habitat suitability
For each species, we predicted the species habitat suitability for the years 1992 and 2015. To represent the change in habitat suitability at the regional level, we then averaged the difference in habitat suitability using a focal statistic of the values with a 5 by 5 matrix of weights (which is equivalent to a 50 km radius), re-assigning the new averaged value to the 10 × 10 km grid. Furthermore, to highlight areas of recent or possible future recolonization, we masked the maps of change in habitat suitability by creating a buffer of 100 km radius around the grid cells of both sporadic and permanent occurrence (Chapron et al., 2014). Finally, to examine the role of different drivers in influencing the spatio-temporal changes in habitat suitability, we estimated a spatial partial response by reprojecting each model at 1992 and 2015 by using the average of all predictors but one along the whole time period. We then repeated the procedure for all the most important predictors.

| Softwares and packages
All GIS operations for land cover data were conducted in GRASS GIS Euclidean distance from urban areas was calculated in R using functions of the GIS software SAGA (Conrad et al., 2015).

| RE SULTS
The optimal selected lambdas (the tuning parameter for the GLMM Lasso models) had a value of 1 for the wolf, 11 for the Eurasian lynx and 3 for the brown bear. The cross-validation demonstrated good predictive performances of the three models, with the Boyce index ranging between 0.878 and 1 for the grey wolf, 0.922 and 1 for the Eurasian lynx, and 0.909 and 1 for the brown bear ( Figure   S10).
The two most important variables contributing the most to the model fit were forest cover and human population density (Figure 2).
Forest cover and human population density showed a positive and negative correlation with the habitat suitability of all three species, respectively ( Figure 3, Table S5). However, forest cover was more important for the lynx, and human population density was more important for the wolf and the bear (Figure 2). The interaction between human population density and its heterogeneity within cells was also important and showed a positive effect in all species, indicating that areas with high population density are more suitable when heterogeneity is also high ( Figure 2, Table S5). Lynx habitat suitability was also slightly negatively correlated with shrubland ( Figure 3, Table S5). Landscape configura-

tion variables (Contagion index, Crop Cohesion and Habitat Diversity)
had little effect, except for habitat diversity, which showed a positive relationship with habitat suitability of the lynx and bear (Figure 3, Table   S5). Implementation of the law showed no effect for wolf, a slightly negative effect for lynx and a slightly positive effect for bear. Hunting regulations had a weak negative effect on the wolf model, but positive for the lynx and the bear models. Finally, the percentage of protected areas coverage had a small but positive impact for the wolf and the bear but negative for the lynx (Figures 2 and 3). The models including cropland instead of forest showed a very similar pattern (Figures S13-S16, Table S6). As expected, cropland had a negative effect on the three species ( Figure S13, Table S6), and however, its effect was only significant for the bear and the lynx ( Figure S14).

| D ISCUSS I ON
In this study, we assessed whether land cover, landscape configuration, human density, protection status and protected area coverage were associated with changes in habitat suitability for large carnivore across Europe during the period between 1992 and 2015. We found that the factors we evaluated had similar effects across the three large carnivore species: habitat suitability is consistently correlated with increasing forest cover, decreasing human population density and cropland cover. As a result of the changes in forest cover and human population density throughout Europe, our model predicts increasing suitability in Eastern Europe and the Balkans, and mixed trends in Scandinavia and Central and Western Europe.
These changes largely match the recent recovery of large carnivores in Europe, but also resulted in many areas inside the known distribution of all three species experiencing a decrease in habitat suitability (e.g. Alps, Pyrenees, Northern Balkans; although not noticeably altering the geographic pattern of habitat suitability for the three species; Figure 4). This contradicts our original expectation of a generalized increase in habitat suitability F I G U R E 4 Map of predicted habitat suitability for the grey wolf, Eurasian lynx and brown bear in Europe in 1992 (a, d, g), in 2015 (b, e, h) and the difference of predicted suitability per cell (c, f, i) between 1992 and 2015. We averaged the value of the difference using a focal statistic of the values within a 50 km radius for facilitating the visualization of the change at the European scale. Masked areas indicate areas beyond 100 km from the species occurrence points collected and permanent distribution areas reported in Chapron et al. (2014) where large carnivore populations have recently recovered. This leads to two important considerations. First, while land cover and human population density are associated with the habitat suitability of the three species, other factors that are not captured by these variables must have played a crucial role in the recent expansion of large carnivores in Europe (Behr et al., 2017;Llaneza et al., 2012). Second, the recent recolonization of areas where we estimated a decrease in habitat suitability suggests that large carnivore species have a high adaptability to habitat, and therefore, habitat availability should not be necessarily considered the limiting factor for the first order of habitat selection (sensu Johnson, 1980) and the subsequent distribution and expansion of these species in Europe (López-Bao et al., 2015. Overall, this would support the idea that these species have a remarkable ability to persist and thrive in human-modified landscapes as long as human tolerance and policy are favourable (Basille et al., 2009;Chapron et al., 2014;Reinhardt et al., 2019). Nevertheless, at finer scales (Johnson, 1980), other factors that do not emerge at the scale of our analyses may be more important for habitat selection (Ciucci et al., 2020;Mancinelli et al., 2018).
Our results suggest that many areas in Europe have less suitable habitat for these species compared to 1992. While there has been an increase in forest cover of about 25% between 1950 and 2010 (Fuchs et al., 2013), the timeframe analysed in our study  was characterized by a slight reversal in this trend, with an overall loss of 0.2% of forest cover in Europe as a whole (Li et al., 2018;Nowosad et al., 2019). In particular, forest cover slightly decreased in some of the important areas for the conservation of large carnivore species, such as the Alps, Pyrenees and Estonia, with several patches losing between 10% and 30% of forest cover (Nowosad et al., 2019). In parallel, human population density has increased overall in Europe (The World Bank, 2019), with decreasing population in many rural areas and an increase in urban and suburban areas, and a shift towards northwestern countries of Europe (EuroStat, 2016). The estimated increase in habitat suitability in Eastern Europe may have fostered the recolonization of Central Europe with more dispersing animals able to settle in rural areas where human population density has decreased and natural vegetation has increased. Furthermore, recent improvement in the suitability of rural areas may have also increased connectivity among large carnivore populations, thus facilitating long-distance dispersal events (Bartoń et al., 2019).

F I G U R E 5
Spatially explicit partial response depicting the temporal trend in habitat suitability for the grey wolf in response to change in human population density (a), forest cover (b), mosaic of forest and shrubland cover (c) and sparse vegetation cover (d); for the Eurasian lynx in response to change in human population density (e) forest cover (f), shrubland cover (g) and Shannon Index (h); for the brown bear in response to change in human population density (i), forest cover (j), shrubland cover (k) and mosaic of cropland cover and natural vegetation (l) Our analysis reveals a lack of association between changes in large carnivore habitat suitability and protected areas, in line with previous findings (Tammeleht et al., 2020). This result can likely be explained by the fact that, on average, European protected areas are smaller than the area requirements of a single individual large carnivore, and are thus expected to play a small role (Linnell et al., 2001;. Furthermore, global analyses have shown that protected areas are typically established in areas previously not impacted by humans (Joppa & Pfaff, 2009;Pressey et al., 1993), which in this context would make them unlikely to contribute to recolonization success. Similarly, the implementation of international conservation laws in Europe did not show any clear effect in our models.
This lack of a consistent pattern could be explained by mismatches between law implementation (i.e. our data) and their actual enforcement, which is likely to vary within and between countries, and temporally. A clear example is the Habitats Directive, which impose the protection of the three large carnivores among many other species of conservation concern (Evans, 2012 Third, large-scale models aimed at estimating habitat suitability through relative probability of presence fall short in accounting for density-dependent habitat selection, which we can expect to vary between areas of stable presence and areas that have been recently recolonized. According to the principle of ideal free distribution (Fretwell & Lucas, 1969), the most suitable habitat is expected to be filled first, and less suitable habitat is expected to be occupied after density in the most suitable habitat increases above a certain threshold. To further complicate this picture, according to ideal despotic distribution, dominant territorial animals may monopolize the most suitable habitat so that low suitable habitat may host a higher density of low ranked individuals (Fretwell, 1972). Finally, in sourcesink dynamics, species may also occur in unsuitable habitat (Pulliam & Danielson, 1991 (Reinhardt et al., 2019). Human tolerance towards large carnivore differs across countries, regions and individuals (Dressel et al., 2015;Gangaas et al., 2013;Glikman et al., 2019;Piédallu et al., 2016). The effort invested in the conservation of these species (e.g. through reintroduction programs) has contributed to their recovery Chapron et al., 2014), sometimes in spite of poor public acceptance (Clark et al., 2002;Tosi et al., 2015). Indeed, recent recolonization of areas historically occupied by large carnivores has also created new conflicts that may compromise conservation outcomes if not adequately managed. This is particularly the case in countries where the species returned after a long period of absence and preventive measures were abandoned, and thus, compensation is favoured over prevention (Bautista et al., 2019).
As a note of caution, the data used as predictors in our models also come with some degree of uncertainty. First, as all land cover products, individual cells in the CCI-LC dataset may be misclas- Despite fine-scale errors, these products are deemed suitable for assessing spatio-temporal changes in land cover over large spatial scales (Huang et al., 2020;Li et al., 2018;Liu et al., 2018;Mousivand & Arsanjani, 2019;Nowosad et al., 2019). Second, human population density data are modelled products, which are deemed suitable for large-scale analyses of human impact (e.g. Venter et al., 2016). However, these layers may hold a level of error when interpreted at fine scale. For example, rural population change and land use change are very much related, and population decrease and rural-urban migration resulted in demographic ageing and land abandonment from the 60s and throughout the late 20th century (Prokopová et al., 2018), first in Eastern Europe and the Mediterranean regions, then in Scandinavia. This widespread land abandonment likely benefitted large carnivore population recovery . Thus, while land cover and human population variables are expected to be collinear at fine scales, at a 10-km resolution and European scale the correlations were low and showed a negligible effect on the results (Appendix S3).
We have shown that recent changes in land cover and human to change. It is likely that habitat availability does not currently limit large carnivore distribution in Europe, and this has important implications for their conservation and management, pointing to scenarios of human-carnivore co-existence and a limited role of the current network of protected areas. While human population growth has slowed in recent decades, the European population is expected to keep increasing in the near future (Eurostat, 2017).
The resulting pattern of redistribution, and its cultural and socio-economic consequences, may have important implications for the occurrence of large carnivores in European landscapes, especially considering that the EU-28 has seen an increase in the rural population in the 2010-2015 period (Eurostat, 2017). An important conservation challenge in Europe will be to capitalize on the socio-economic and landscape changes to create new opportunities for species to recover (Ceausu et al., 2015), as well as active education, legislation and management to mitigate human-wildlife conflicts in newly recolonized areas Rigg et al., 2011;Treves & Karanth, 2003).

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13219.

DATA AVA I L A B I L I T Y S TAT E M E N T
Codes and datasets used for the analyses are openly available on