Aim: 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. Location: Germany. Methods: Using telemetry data as presence points, we compared the predictive performance of five modelling approaches based on combinations of distribution modelling algorithms—GLMM, MaxEnt and ensemble modelling—and two background point selection strategies. We used a homogeneous Poisson point process to draw 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 non-stationarity in habitat selection pointed to wolves settling in the best areas first, and filling in nearby lower-quality habitat as the population increases. Our results caution against spatial extrapolation and space-for-time substitutions in habitat models, at least with expanding species.