Understanding habitat selection of range‐expanding populations of large carnivores: 20 years of grey wolves (Canis lupus) recolonizing Germany

The non‐stationarity in habitat selection of expanding populations poses a significant challenge for spatial forecasting. Focusing on the grey wolf (Canis lupus) natural recolonization of Germany, we compared the performance of different distribution modelling approaches for predicting habitat suitability in unoccupied areas. Furthermore, we analysed whether grey wolf showed non‐stationarity in habitat selection in newly colonized areas, which will impact the predictions for potential habitat.

background points from either the minimum convex polygons derived from telemetry or the whole area known to be occupied by wolves.Models were fit to the data of the first years and validated against independent data representing the expansion of the species.The best-performing approach was then used to further investigate non-stationarity in the species' response in spatiotemporal restricted datasets that represented different colonization steps.
Results: While all approaches performed similarly when evaluated against a subset of the data used to fit the models, the ensemble model based on integrated data performed best when predicting range expansion.Models for subsequent colonization steps differed substantially from the global model, highlighting the non-stationarity of wolf habitat selection towards human disturbance during the colonization process.
Main Conclusions: While telemetry-only data overfitted the models, using all available datasets increased the reliability of the range expansion forecasts.The

| INTRODUC TI ON
There is an increasing need to restore wildlife populations and ecosystems in the 21st century (United Nations, https:// www.decad eonre stora tion.org), given the ongoing defaunation in the Anthropocene (Dirzo et al., 2014;Schmitz et al., 2023;Young et al., 2016).This holds particularly true for the world's large carnivores, who play key functional roles in ecosystems (Ripple et al., 2014), yet have been extirpated over much of their historical ranges in the past (Chapron et al., 2014).Restoring these species can bring major benefits in the context of rewilding and to human societies (Fernández et al., 2017;Perino et al., 2019), but also lead to increasing human-carnivore conflicts.As a result, predicting where and how large carnivore populations expand is important for their management and conservation, and ultimately for fostering coexistence between people and carnivores (Treves et al., 2004).Europe, where after centuries of persecution and local extirpation large carnivores have expanded since the 1990s-and particularly the case of the grey wolf (Canis lupus, hereafter 'wolf')-provides a unique, but largely unexplored natural experiment to understand how large carnivore recolonization processes take place (Chapron et al., 2014;Cimatti et al., 2021;Perino et al., 2019).
While the range expansion of wolves in Europe is a recent phenomenon that has been studied at broad scales (Nowak et al., 2017;Reinhardt et al., 2019), wolf recolonizations in North America have been studied in more detail scales (Gantchoff et al., 2022;Lesmerises et al., 2012;e.g. Mladenoff et al., 2009;Ripple & Beschta, 2012).For example, low road density and natural cover determined range expansion and habitat selection in the Midwest of the United States (Gantchoff et al., 2022;Mladenoff et al., 2009), whereas forest cover and food resources have been important in boreal forests (Lesmerises et al., 2012).Avoidance of high human-density areas is also a common pattern (van den Bosch et al., 2022).Analyses of North American wolf habitat selection patterns over time and space have thus shown considerable variability, possibly related to the increasing wolf population size (Mladenoff et al., 1997(Mladenoff et al., , 2009)).When the initial wolf population is establishing in an area, it may select the most suitable habitat first, with lower-quality habitat subsequently filled as the number of packs increases (O'Neil et al., 2020).However, the landscapes which are colonized by wolves in North America and Europe differ hugely, with higher human population density, human pressure and habitat fragmentation in Europe, and it thus remains unclear whether the insights gained for the wolf expansion in North America can be transferred to Europe and other situations.For Europe, different attempts have been made to assess potential wolf habitat during the expansion of the populations (Fechter & Storch, 2014;Jedrzejewski et al., 2008;Louvrier, Chambert, et al., 2018;Marucco & McIntire, 2010;Nowak et al., 2017), but there is a huge variability in habitat selection factors found across these studies, from road avoidance and forest selection to altitude or farmland areas.It remains unclear whether this variation reflects wolf preferences, differences in environmental conditions or food availability, or whether it is an artefact of the modelling process, especially due to data obtained at different phases of the expansion process.
Forecasting species range expansion is not easy as it requires accounting for the potentially changing behaviour of the species in space and over time (O'Neil et al., 2020).While species distribution models (SDMs) have been shown to be very useful for delineating potentially suitable habitat for species based on environmental predictors (Elith & Leathwick, 2009;Guillera-Arroita et al., 2015), these models rely on the assumption of stationary processes to predict potential future areas for range expansion.
In the context of habitat selection, stationarity refers to a constant response of the species to environmental variables across time and space.Making reliable forecasts under non-stationary conditions-in other words, when the individuals react differently to certain environmental variables as may be the case of populations expanding beyond their current distribution-, may risk yielding spurious results due to poorly specified models as well as inappropriate spatial inference (Osborne et al., 2007;Rollinson et al., 2021).Non-stationarity is important in two ways.First, non-stationarity in environmental conditions occurs when the environmental variables defining the currently occupied habitat of a species differ from those defining habitat that the species can expand into.This makes the transferability of models a difficult task (Yates et al., 2018), leading to a potential underestimation of range expansion.Second, the habitat selection of individuals during a process of population expansion can be non-stationary, leading to changing functional responses of individuals towards environmental factors over time (Matthiopoulos et al., 2011;Venne & Currie, 2021).These functional responses subsume several mechanisms behind selection (such as predator escape or foraging habitat selection), and hence are correlative simplifications, difficult to transfer from one situation to the next.Recently, it has been suggested that habitat models incorporating data from many populations are often better suited for assessing range expansion potential than models based on individual populations (Kuemmerle et al., 2018;Zurell et al., 2018).Robust modelling techniques and careful cross-validation in space and time are needed to deal with non-stationarity.Data integration tools (Zipkin et al., 2017(Zipkin et al., , 2021) ) and ensemble modelling approaches (Liu et al., 2019;Thuiller et al., 2009) corroborate considerable promise in this regard.
Here, we analyse such models' predictive ability to capture range expansions, using the example of a rapidly expanding grey wolf population in Germany.The growth of the wolf population from one pack in the year 2000 to approximately 130 packs over a 20-year period (DBBW 2021; www.dbb-wolf.de/ home) is a major conservation success regarding the natural recolonization of the former range of the species, but also a challenge for co-existence of people and carnivores in shared landscapes.On one hand, better wolf protection laws effectively restored wolf populations, bringing back key ecological functions given their role as top predators, and might have positive biodiversity effects by controlling herbivore populations or providing resources to scavenger species (Ripple & Beschta, 2012;Wilmers et al., 2003).On the other hand, the presence of wolves may increase conflicts with farmers and hunters (König et al., 2020;Treves et al., 2011).These positive and negative aspects of returning wolves typically result in a highly polarized public attitude towards this iconic species (Arbieu et al., 2019;Treves et al., 2004).Although wolves are protected in most European countries including Germany (Chapron et al., 2014), illegal killing as retaliation for livestock depredation, in fear of competition for game and fear of wolves became a major source of wolf mortality (Ansorge et al., 2006;Bautista et al., 2019;Linnell & Alleau, 2016;Nowak et al., 2021;Sunde et al., 2021;Suutarinen & Kojola, 2017).Understanding where wolves are likely to expand is therefore urgently needed to proactively manage potential conflicts and to foster co-existence.
In this study, we aimed at determining the best modelling strategy to obtain a reliable distribution map for potential habitat outside the current species range in an expanding population.To answer this, we (1) tested different modelling approaches-combinations of datasets and modelling algorithms-to predict the expansion of the population and validated the model fit with independent data and (2) analysed whether data from newly colonized areas showed non-stationarity in habitat selection that could affect the predictions for the species' range expansion forecasts.Specifically, we used telemetry data and monitoring locations of territories documenting the wolf population expansion in Germany for the first 20 years of population establishment, while addressing two main questions: (1) How do data from different data sources affect model performance in an expanding population?and (2) Can we detect non-stationarity in the species' response to the environmental variables at different colonization steps, and how does this affect predictions of habitat suitability?

| Study area
Our study region was the highly human-dominated landscape of Germany in Central Europe.Germany has a population of circa 83 million inhabitants (2018) on an area of about 358,000 km 2 , resulting in an average human population density of 232 people/km 2 .About 14.2% of the country's surface is occupied by settlements and roads, with a total length of 13,000 km of motorways and 125,000 km of interstate and main roads, while forest areas represent about 30%.Although Germany has 16 National Parks, covering a terrestrial area of approximately 2150 km 2 (Federal Agency for Nature Conservation-BfN, www.bfn.de/ en), this only constitutes 0.6% of Germany's land surface.

| Wolf data
The grey wolf is currently recolonizing Germany, with the first territories with pups established in Saxony, eastern Germany, in 2000(Reinhardt et al., 2021).Since then, the population has spread across mainly the northern part of Germany, with a few individuals dispersing to the south (Appendix S2 in Supporting Information, Figure S2.1).
Our wolf data comprised two complementary datasets.Our first dataset consisted of GPS telemetry locations of 20 collared resident wolves from 2009 to 2018 (hereafter 'telemetry data').
Manipulation of individuals in the field was done following a strict ethical approach and within approved projects (Appendix S2: Table S2.1).We did not consider dispersing individuals for the analyses as our focus was on the establishment of territories.
Three individuals (FT2, MT3, WF5) dispersed after collaring from their natal to new breeding territories, in which case only data from the newly established territories were used, as natal territory data were already available from a parent wolf.For ID1, data from its natal (ID1) and breeding territory (ID1p) were included, as no parent data were available.Additional information on the collared individuals is provided in Appendix S2: Table S2.2.The final telemetry dataset consisted of 3841 locations from 21 home ranges (12 adults and 9 subadults including one home range shift, 183 ± 104 locations per home range) followed between 2 and 23 months (Appendix S2: Figures S2.2 and S2.3).The average home range size, defined as the 95% Minimum Convex Polygon (MCP), was 234.4 ± 129.0 SD km 2 (range 45.8-461.8km 2 ).GPS collars provided locations every 1-4 h.Thus, to avoid temporal and spatial autocorrelation, we selected one location randomly per day (Hiller et al., 2015) and further filtered the remaining data to one location per grid cell of 100 × 100 m (resolution of the environmental variables, see below).
Our second dataset consisted of centroids of known wolf territories monitored annually since 2000 (www.dbb-wolf.de/ Wolfs vorko mmen/ terri torien/ karte -der-terri torien; hereafter 'monitoring data') (Appendix S2: Figure S2.1).The centroids of the territories were assessed after the end of the monitoring year as the central point of all activity signs (scats; camera trap images; telemetry data when available; opportunistic sightings; hunt remains) that were assigned to the same territory; however, GPS locations of these activity signs were not available.Territory areas for our study were delineated with a radius of 8 km around their centroids, resulting in territory sizes of approx.200 km 2 .This refers to the standard procedure of the Germany monitoring protocol (Reinhardt et al., 2015).

| Environmental variables
Wolves in Central and Northern Europe occur mainly in areas with high forest cover, low density of roads and low human population density (Cimatti et al., 2021;Jedrzejewski et al., 2008;Louvrier, Chambert, et al., 2018), although they can also use open areas such as mineral extraction sites (Ansorge et al., 2006).These variables are consistent with the species preferences in other areas of their range, where they favour forested and low-density road areas (Mladenoff et al., 1999(Mladenoff et al., , 2009)).Therefore, we used five environmental layers to represent wolf habitat suitability, with four of them capturing human disturbance factors: Human population density (Pop_den), Human footprint (HFP), Distance to settlements (Dist_settl) and Distance to roads (Dist_roads) (Appendix S2: Table S2.3).Natural environment was represented by Corine Land Cover data (CLC2012), which were reclassified into seven classes based on habitat features relevant for wolves (Appendix S2: Table S2.4):forest, green area with cover (i.e.forest), green area without cover, agriculture, urban, water bodies and bare areas without cover.
We rasterized and resampled all environmental layers to a 100 × 100 m resolution for the analyses and tested for multicollinearity using Pearson's correlation coefficient (r < |.7| in all cases, Appendix S2: Figure S2.4).After the habitat suitability modelling, we aggregated the resulting maps to a coarse resolution of 10 × 10 km 2 following the EU Reference Grid (https:// esdac.jrc.ec.europa.eu).
We did this for visualization and safety reasons, to avoid showing exact locations of high suitability to avoid illegal killings.We kept a 10-km buffer around the German border to account for suitable areas in neighbouring countries that are connected to suitable areas within Germany.

| Predicting habitat suitability for the expanding wolf population
We applied a 'use versus availability' approach, where the GPS locations of the individuals represent the presence of the species that were compared with background points representing the available range of values for the explanatory variables (Elith et al., 2011).As we were modelling the distribution of a species that is expanding its range, we used a two-step approach (Figure 1).First, we evaluated different modelling methods and data selection procedures to find the best-performing approach in predicting habitat suitability for new areas by validating our model with the most recent data not used during the model fitting procedure (Step 1, see below).Then, we used that modelling approach to obtain the final habitat suitability map for the species using all available data (Step 2, see below).

| Step 1: Modelling habitat suitability for an expanding species
Predicting the habitat suitability of an expanding species comes with multiple challenges.For starters, predictions outside the range used to create the distribution model might be unreliable (Guillera-Arroita et al., 2015).Additionally, the area selected for the background point extraction can influence model outcomes (Kramer-Schadt et al., 2013;Merow et al., 2013;Phillips et al., 2009).With expanding wolf populations being such a sensitive topic due to human-wildlife conflicts and the strong effect of management measures for the species conservation (e.g.Cimatti et al., 2021), our main goal was to obtain a model as reliable as possible.
For that reason, we tested five modelling approaches, representing a gradient in the background selection from an individual-specific to a species-broad approach and compared single model algorithms with an ensemble model approach (Figure 1 4. an ensemble modelling approach (Liu et al., 2019;Marmion et al., 2009), as it has been proposed to outperform the individual algorithms (Breiner et al., 2015), based on the presence and background data from the individual home ranges ('telemetry model-Ensemble'); and 5. an ensemble model using the same approach as before, but with the background points obtained from the areas occupied by the In all cases, the presences were obtained from the telemetry data as the GPS locations of the collared individuals, temporally thinned to one location per day to avoid the effects of temporal autocorrelation.For the individual-id analyses, background points were selected within the 100% Minimum Convex Polygon (MCP100) derived for each individual.For both telemetry models, background points were extracted following a homogeneous Poisson point process within the MCP100.sThe integrated models combined the GPS and monitoring data by using the thinned collared data as presences and the known territories from the monitoring scheme as available background areas to randomly extract the background points also following a homogeneous Poisson point process (details in Appendix S2: Section S2.1).In all cases, we extracted 10 times more background points than observations.Both presences and background points were spatially filtered to retain no more than one location per grid cell to avoid overrepresentation of environmental variables (Table 1).
To assess model performance for the range expansion, we set aside the most recent data (after year 2018) from all datasets and used two different evaluation approaches.First, we ran the models with data from before 2018, separated into training (80%) and test (20%) sets.This test set was used to evaluate the model performance within the observed range of the species ('test1-evaluation data', Table 1).
Second, we used the most recent data (data from 2018 to 2020) to validate the model's performance including those areas where the species expanded to ('test2-validation data', Table 1).We used receiver operating characteristic area under the curve (AUC, Elith  3) delineate the area of data origin for the regional models that correspond with step 2 of our modelling approach to assess non-stationarity in habitat selection during the range expansion (cf. Figure 2).AUC, area under the curve; TSS, true skill statistic.et al., 2006), true skill statistic (TSS) based on maximized sensitivity and specificity (Allouche et al., 2006;Liu et al., 2013) and continuous Boyce index values (Hirzel et al., 2006) on the evaluation and validation test data sets.All SDMs were run in R v. 4.0 (R Core Team, 2020), using the packages lme4 (Bates et al., 2015) for GLMMs and biomod2 for MaxEnt and ensemble distribution models (Thuiller et al., 2023).

| Step 2: Predicting suitable wolf habitat across Germany ('global model')
After the validation, we used the best-performing approach (integrated model-ensemble) to obtain the final wolf distribution in the study area, now using all the available data (until 2020) to obtain a model as reliable as possible (Figure 1.2).In this case, given the high number of presence points in the full dataset, we used the same number of random background points as presences to optimize algorithm performance (Liu et al., 2019).We ran each model 10 times using a random selection of points in each repetition, where 80% were assigned as training set and 20% of the points were assigned as test set.The final ensemble model (hereafter 'global model') was built as a weighted average based on AUC, with the threshold for inclusion greater or equal to AUC = 0.7 to account only for the bestperforming models.

| Analysis of non-stationarity in the responses during colonization steps ('regional models')
We tested for non-stationarity in wolf habitat selection by parameterizing models representing different colonization steps of the expanding population (Figure 1.3).We divided the data into three different regions, representing different temporal colonization  Telemetry data filtered to one random observation per day and per grid cell.
b Random points within 100% MCP, 10× presences, filtered to one per grid cell.
c Random points with 0.99 probability in 100% MCP and 0.01 outside MCPs, filtered one point per grid cell.
d Random points with 0.99 probability in 10 km buffer territories and 0.01 outside buffer, filtered one point per grid cell.
steps (Jarausch et al., 2021).Region 1 refers to the area that was To analyse if regional differences were caused by individuals' breeding status rather than region, we additionally ran two binomial GLMMs with individual-id as a random effect and the environmental variables as predictors in interaction with breeding status (breeder, non-breeder) or colonization step (region 1, region 2, region 3).We compared both models using Akaike Information Criteria for small sample size (AICc), which supported the use of region instead of breeding status of the individuals (ΔAICc = 4344.56,Appendix S2: Tables S2.5 and S2.6; Figure S2.6).

| Habitat suitability effects on reproduction and assessment of territory numbers
To assess how habitat suitability might be associated with the reproductive success of the recolonizing population, we ran an exploratory analysis based on the number of reproductive events recorded in each territory monitored, based on opportunistic data obtained from the monitoring surveys.For this analysis, we fitted a GLM using the number of reproduction event per territory as the response variable and the mean habitat suitability of the territory as an explanatory variable, with a log link and a Poisson error distribution.
Finally, we transformed predictions into a binary map to estimate the potential number of wolf territories in Germany, assuming maximum occupancy of the study area and territory sizes of ~200 km 2 (Reinhardt & Kluth, 2016).For this, we aggregated the predictions of the original distribution model consisting of a 100 × 100 m resolution to a coarser resolution of 10 × 10 km, from which it was possible to assign wolf territories by selecting two consecutive cells.We used the existing monitoring territories of the wolf population and calculated mean habitat suitability within these territories and their associated variance as a reference value.
We then generated two binary maps, a conservative map (mean habitat suitability minus one standard deviation), and a more optimistic map based on the lower 95% confidence interval of the mean habitat suitability value.Isolated cells that were too small to contain a wolf territory (i.e., single, unconnected 10 × 10 km 2 cells) were removed.

| RE SULTS
Almost all modelling approaches performed adequately to predict wolf habitat suitability for the period used to fit the models, although their performance dropped when validated against the expansion data (i.e. the most recent data of the expanding population) (  The exploratory analysis of the number of reproductive events showed a clear positive association between the total number of reproductive events and higher habitat quality (β = 2.60 ± 0.54, p < .001, Figure 3).The mean habitat suitability within the existing monitoring territory was 0.386 ± 0.107.Finally, the binary territory maps from the global ensemble model estimated the number of potential wolf territories in Germany to range between 652 and 1407 (Figure 4), depending on the threshold variance used for habitat suitability values.These numbers are similar to a prior exploratory study with only one modelling algorithm and dataset (Kramer-Schadt et al., 2020).

| DISCUSS ION
Large carnivores provide key ecosystem functions but are also associated with human-wildlife conflicts.Forecasting range expansions of large carnivores into human-dominated landscapes can therefore help to proactively devise strategies that foster co-existence.
However, predicting range expansions robustly is a very challenging task as SDMs typically assume species' responses to environmental factors to be constant over time and in space, which may not be the case for range-expanding species.This can make predictions a tightrope walk towards accuracy and predictive power.Using the natural experiment of an expanding wolf population in Germany, we showcased the best modelling approach to obtain reliable distribution models for an expanding population.Furthermore, we revealed non-stationary habitat selection in the expanding wolf population, which could hide potential habitats if not accounted for during the data-collecting process.
The non-stationarity in the response was likely explained by the increasing number of wolves in regions with long-time presence of the species, as region 1 was the longest-occupied one and contained the highest density of individuals (www.dbb-wolf.de).This might have resulted in the use of lower-quality habitat once high-quality habitat is occupied as the population increases (Mladenoff et al., 2009;O'Neil et al., 2020).Specifically, wolves appear to establish packs first in high-quality areas characterized by forest cover and lower human disturbance, which was shown in the models for the recently colonized areas (regions 2 and 3), filling areas with increasing human disturbance when the species density increases and no more high-quality habitat is available.This last case is clearly shown in the models based on region 1, which was occupied for 20 years and shows territories in areas with higher human disturbance values.
The observed differences in habitat suitability between the regional models might point to a flexible territory selection of the species depending on the available habitat that is not already occupied by conspecifics.
Differences in habitat quality are expected to translate into differences in fitness (Mosser et al., 2009;O'Neil et al., 2017).
Furthermore, the flexible habitat selection of the species from high-to low-quality habitats as the number of individuals increases may also create source-sink dynamics in the population, with packs in low-quality areas maintained by immigration from high-quality areas (Mladenoff et al., 2009;O'Neil et al., 2017).
Both, the difference in fitness and potential source-sink dynamics are partially supported by our exploratory analysis on reproduction, as the more suitable habitats were associated with a higher number of reproductive events.However, this result should be interpreted with caution, as we do not have information about pup survival and the opportunistic nature of our reproduction data might bias results.
Methodologically, we showed that the combination of different data sources is highly beneficial, as models using both datatypes (telemetry data for the presence of the species and territory monitoring Additionally, the ensemble modelling approach outperformed the single MaxEnt algorithm, especially when validated against the data from the population expanding to new areas.Although ensemble approaches are not necessarily better than other models (Hao et al., 2019), in the particular case of the expanding population, single algorithms like MaxEnt might again overfit the model (Fourcade et al., 2014;Merckx et al., 2011) and play against the ability of this approach to extrapolate the predictions to new territories.By combining multiple algorithms, the ensemble model might overcome some of the limitations that can appear in single algorithm models (Norberg et al., 2019), especially for cases like large carnivores, where usually data are scarce (Breiner et al., 2015).
We showed how validating models against quasi-independent data, collected from periods after the training and test data, helps to properly assess the predictive power of the model, as suggested by Treves et al. (2011) and Treves and Rabenhorst (2017).Models that relied only on one dataset performed overall poorly in predicting range expansion.The model based only on telemetry data underestimated the potential distribution of wolves as it was highly dependent on a single individual's habitat selection and thus, a limited representation of the population (Milakovic et al., 2011).Thus, monitoring programs will benefit from purposefully gathering different types of data, and we advise agencies to use a combination of extensive and intensive methods (e.g., territory mapping and telemetry data, respectively) and to accurately assign a GPS location to each activity sign.
Our results highlight the importance of assessing the population status of expanding carnivore populations before attempting to forecast range expansion.As already shown in the case of hunting site selection by wolves, the species has flexible behaviour, thus predictive modelling should be done very carefully (Treves et al., 2011).This also means that the commonly used Representation of the number of reproductive events registered at each monitoring territory related to mean habitat quality (x-axis) and the standard deviation or the habitat quality (y-axis).Each pie chart represents one monitoring territory, with the size representing the total number of years that the territory was monitored.The histograms over the axis represent the number of events per value of habitat suitability.
space-for-time approach could be problematic in range-expansion models, as habitat selection might vary across different regions and/or at different time steps.Although differences in reproductive status could be partially responsible for the observed non-stationarity in wolf habitat selection, our suggestion that non-stationarity was due to colonization time steps was supported by the fact that the models accounting for regions performed much better than the model including the breeding status of individuals (Appendix S2: Tables S2.4 and S2.5).
A number of conservation and management implications derive from our work.Wolf populations show a leap-frogging colonization pattern, based on long-distance dispersal events and subsequent spreading and filling-in of habitat in the surrounding, initially occupied, prime habitat (Jarausch et al., 2021;Louvrier, Duchamp, et al., 2018;Reinhardt et al., 2019).The historical expansion of wolves in Germany followed an east to north-west path, driven by wolves' preference for forested areas away from human disturbance (Cimatti et al., 2021;Louvrier, Chambert, et al., 2018;e.g. Mladenoff et al., 2009).The first individuals arrived from Poland, following the expansion of the local population (Czarnomska et al., 2013;Reinhardt et al., 2021) and quickly expanded over the available habitat.In accordance with other European analysis of wolf habitat selection (Fabbri et al., 2007;Jedrzejewski et al., 2008;Louvrier, Chambert, et al., 2018;Ražen et al., 2016), our results show that wolf-preferred natural areas and avoided high human disturbance.We also identified suitable habitat in remote, forested,  According to our binary territory maps, the number of potential territories (assuming a territory size of 200 km 2 ) ranges from ~650 to 1400 in Germany, compared to a current number of 158 packs, 38 pairs and 22 solitary territorial wolves (monitoring year 2020/21; www.dbbw-wolf.de).Prey density could be a limiting factor, although we did not explicitly include it in our models.The main prey of wolves in Germany are wild ungulates, and especially roe deer, which accounts for more than 50% of wolf diet (Wagner et al., 2012).There were three reasons why this variable was left out.First, to our knowledge, there is no reliable, fine-scale information on ungulate abundance across Germany.Second, we expect ungulate densities in the study area to be highly correlated with land cover, especially forest cover (Heurich et al., 2015).
Third, prey density is very high across much of Europe, including Germany (Massei et al., 2015;Reinhardt et al., 2021).In a regional analysis of prey abundance versus wolf abundance (Appendix S3), we did not find any effect of wolf pack numbers on the population trends of any prey species.We furthermore note that prey numbers are currently increasing in Germany since the 1980s (Burbaite & Csányi, 2009, German hunters association, 2022).

| CON CLUS IONS
We highlighted the importance of explicitly accounting for nonstationary habitat selection and, when modelling range-expanding carnivores, that species are not in equilibrium with their environment.Forecasting range expansions of species, and particularly of generalist species such as many large carnivores, comes with high levels of uncertainty.As limiting as this may seem, analysing range expansions can help tremendously to improve our understanding of species' ecology, and more generally about the usefulness and limitations of SDMs.Analysing the natural experiment of an expanding wolf population in Germany, we here demonstrated the advantage of integrating diverse occurrence datasets to increase the predictive power of models to forecast range expansions.Likewise, we uncover major non-stationarity in wolf habitat selection over time (i.e.early colonization phase vs. established populations), resulting in considerable bias in habitat predictions if non-stationarity is not considered.Our results highlight that habitat models based on data from recent colonized areas will underestimate habitat for rangeexpanding species, and instead all available data, and particularly those from areas where populations are already established should be used.Furthermore, we showed how the combination of multiple data sources and the use of ensemble modelling approaches improve the performance of the distribution models with the aim to predict new areas of expansion, therefore we recommend monitoring agencies to collect data on all possible signs of the species presence to increase prediction accuracy.

ACK N O WLE D G E M ENTS
We are thankful to Hermann Ansorge, Markus Ritz (both SGN Görlitz) and Sandra Balzer (BfN) for their support during the study.
For providing the wolf territory data, we thank the authorities in the

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest.

PEER R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/ api/ gatew ay/ wos/ peer-review/ 10. 1111/ ddi.13789 .

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in this manuscript is available in the Dryad reposi-

O RCI D
Aimara Planillo https://orcid.org/0000-0001-6763-9923 .1).Specifically, we tested: 1. a generalized linear mixed-effects model (GLMM) accounting for the individual's identity ('individual-id model') with background points obtained from each individual's home range; 2. a machine learning algorithm (MaxEnt; Elith et al., 2011; Phillips et al., 2006) as a common approach for SDMs favoured by the scientific community (Fourcade et al., 2014) using GPS locations as presences and background points obtained from the individual home ranges, but without individual identity ('telemetry model-MaxEnt'); 3. a similar MaxEnt model, but using as background points locations extracted from the areas occupied by the species as determined by the German monitoring program ('integrated model-MaxEnt'); species, as determined by the monitoring program ('integrated model-Ensemble').For the ensemble model, we used seven modelling algorithms including regression (Boosted Regression Trees, Multiple Adaptive Regression Splines, Generalized Linear Models and Generalized Additive Models) and machine-learning methods (Random Forest, MaxEnt and Artificial Neural Network).All models used a fine spatial resolution of 100 × 100 m, representing the resolution of the environmental data associated with the GPS locations used as presences.

F
I G U R E 1 Modelling framework highlighting the main analyses and modelling steps.Panel (1) shows step 1 of the modelling approaches related to first finding the best data structure and modelling algorithm.Panel (2) shows the final global model.Red ellipses in panel ( a first colonized by wolves, with a population established in 2005 and present in the region since then, and which holds the highest density of individuals.Region 2 groups individuals from a second, intermediate colonization phase that was the result of a long dispersal event from region 1 around 2009.Finally, region 3 groups individuals from the most recently established western territories, which were colonized after 2012 following another long dispersal event (Appendix S2: Figure S2.5).Wolf occurrence data in these re-gions covered areas of 1582 km 2 in region 1, 957 km 2 in region 2, and 810 km 2 in region 3. Background data were subsampled as random points in a 50-km buffer around the 95% MCPs in each region, and an ensemble model was run with the same approach as the global model described above.Each of the resulting regional models was then subtracted from the global model to highlight discrepancies between these models.We assessed the extend of the environmental differences between model training and the predictions outside the observed range using multivariate environmental similarity surfaces(MESS, Elith et al., 2010), highlighting areas of higher uncertainty in the predictions.
the model that accounted for individual-id('individual-id model')    and the MaxEnt model based on telemetry data ('telemetry model-MaxEnt').On the other hand, the ensemble models had the highest scores, with the highest AUC value for the ensemble model based on the combined data from telemetry and monitoring territory locations ('integrated model-Ensemble', Table2).The global model based on all available data and the approach 'integrated model-Ensemble' identified the Land Cover (class Forest) and the distance to roads as the most relevant variables, both with a positive effect on wolf habitat suitability.This model predicted areas of high habitat suitability in the north-east and south of the existing wolf population, with lower habitat quality towards the western part of Germany (Figure2, top left).Considerable suitable habitat areas that are currently not occupied occurred in southern Bavaria and the forested areas in central Germany (the Harz and Odenwald-Spessart-Rhön lower mountain ranges).Regarding the non-stationarity in the responses, the habitat suitability models built from data exclusively from each of the areas related to the colonization steps (regions 1-3) varied widely in their habitat suitability values (Figure 2, second row), showing a clear change in the species selection across the expansion.While the general patterns remained similar, the model based on data from region 1 (earliest colonization) was more optimistic in its predictions of habitat suitability than the global model, especially in the northern part of Germany, assigning higher values to areas closer to roads.The MESS analysis showed that the ranges of environmental values represented in the data from region 1 were similar to the ranges of values in the predicted area, suggesting it to be a good representation of the country (Figure 2, third and fourth rows).The habitat prediction based on wolf data from region 2 (intermediate colonization step) showed generally lower habitat suitability than the global model, although with higher suitability values in the southern part of Germany.Finally, the habitat suitability map derived from the model with wolf occurrence data from region 3-the most recently colonized area-underestimated habitat suitability for wolves across most of Germany.This model also accounted for the least environmental variability throughout the country, as environmental values associated with data from region 3 did not represent the full range of possible values across the country, especially for higher human disturbance areas (Figure 2, bottom right plot).

F
I G U R E 2 Wolf ensemble habitat suitability models based on the combined datasets.Top left, model based on the global dataset.Top right, definition of the three different regions based on the colonization steps.Rows two to four show the habitat suitability maps (second row), the differences with the global map (third row) and the multivariate environmental similarity surfaces (MESS, fourth row) for the regional ensemble models based on region 1 (left), region 2 (middle) and region 3 (right).Differences to global model: orange colours show areas with habitat quality over-predicted by the regional model.MESS: negative values in red show environmental values outside of the range used in training the model.data as a representation of the background or availability) clearly outperformed models based only on one dataset in terms of predictive performance(Renner et al., 2019;Zipkin et al., 2021).While models using background points obtained from the MCP of the individuals performed very well for the current conditions, their accuracy when extrapolated to new areas clearly decreased, pointing to an overfitting of the model to the specific individuals used during the fitting process.In contrast, using background from all areas with the known presence of the species (represented by the territories in the monitoring data) provided the model with enough flexibility to predict the expansion.

F
Binary maps showing suitable wolf territory areas (orange) in Germany.Two thresholds were applied for binary conversion: (a) mean minus one standard deviation of habitat suitability values in the monitored territories and (b) lower 95% CI of these values.

Federal
States responsible for wolf monitoring.SKS is associated with the DFG RTG 2118 BioMove.We thank Jamie Burton for proofreading the manuscript.Open Access funding enabled and organized by Projekt DEAL.FU N D I N G I N FO R M ATI O N This research was partially funded by the Federal Agency for Nature Conservation (BfN) and the Federal Ministry for the Environment, Nature Conservation and Nuclear Safety (BMU) in the frame of the project 'Federal Documentation and Consultation Centre on Wolves' (FKZ 3515 82 4100).Telemetry projects were supported by multiple institutions: own funds from HNEE, WZI (Wolfskompetenzzentrum Iden-Gluecksburger Heide), SUNK-Stiftung (Stiftung Umwelt, Natur-und Klimaschutz des Landes Sachsen-Anhalt-Oranienbaumer Heide), Federal Forest Services, Bundesforstbetrieb Mittelelbe and Primigenius gGmbH, LAU (state authority for environmental protection), BfN, BMU, SMUL (Saxony State Ministry of Environment), WWF, NABU (Naturschutzbund Deutschland), GzSdW (Gesellschaft zum Schutz der Wölfe) and IFAW (International Fund for Animal Welfare), NLWKN (Lower Saxony Water Management, Coastal Defence and Nature Conservation Agency), Ministry of Agriculture and Environment in Mecklenburg-Western Pomerania, The German Hunting Association and The Circle of Friends of wild Wolves.AP was partly supported by the Bridging in Biodiversity Science-BIBS-project funded by the German Federal Ministry of Education and Research BMBF (funding number 01LC1501).The publication of this article was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)project number 491292795.
tory: https:// doi.org/ 10. 5061/ dryad.m63xs j461.Telemetry data of wolf individuals is provided at 50 km resolution due to concerns of illegal killing of the species, which is protected in Europe under Bern Convention Annex II (strictly protected species) and Habitats Directive, and in Germany under the Federal Nature Conservation Act ( §44 BNatschG 29.07.2009).Additionally, the r code used in the analysis can be found in GitHub repository: https:// github.com/ aplan illo/ 2023_ WolfH abita tGermany.

TA B L E 1
Number of points and environmental values used in each modelling approach comparing different datasets.

Data type n points Main land use pop_den HFP Dist_settl [m] Dist_roads [m] Mean SD Mean SD Mean SD Mean SD Individual-id model
Note: Data type refers to each of the data subsets used during the modelling process: before 2018 training (80%) and test (20%) evaluation set ('test1Eval'), and after 2018 expansion validation set ('test2Val'), with the second part of the name '_pres' and '_bg' identifying presence and background points, respectively.The superscripts refer to the origin of each subset of data.Environmental variables are explained in Appendix S2: TablesS2.3andS2.4.Main land use 104 = forest, 102 = agricultural areas.

Table 2
Evaluation metrics for the five modelling approaches based on different background selection and algorithms.
).While Boyce index and TSS values showed little variation in the evaluation of the different datasets, AUC values showed larger differences, pointing to this metric as the most sensitive score for model performance and thus, we base the following results on the AUC values.When evaluated against new independent data (test2-validation data), the largest performance drops occurred inTA B L E 2Note: Test 1-evaluation data (before 2018) corresponds to the data used to fit the model.Test 2-validation data (after 2018) corresponds with the most recent data from the expanding population that includes new areas and was not used for model fitting.For MaxEnt models, AUC and TSS values correspond with the average value of the 10 model replicates.For ensemble models, the ensemble predictions were evaluated.Abbreviations: AUC, area under the curve; TSS, true skill statistic.