The wolf and the city: insights on wolves conservation in the anthropocene

After a long period of human persecution that led it to extinction in most of its distribution range, the wolf is undergoing a fast recovery. Despite being described as an elusive species only living in remote areas, wolves are recently occupying also human-dominated landscapes, increasing the frequency of direct contacts with humans. Nevertheless, it is unclear whether this situation is only caused by a numerical increase or partially facilitated by a higher tolerance of wolves with respect to human proximity. We focused on a European region offering an abundant, widespread and long-term monitored wolf population to analyse wolf pack distribution and its relationship with human presence across areas recolonized across four different time periods (1972, 1996, 2005 and 2016). Throughout areas recolonized during different periods, wolves were initially located in mountains and hills, occupying plains only in the recent past. Although they consistently tended to be located as far as possible from urban settlements, especially from those with over 5000 inhabitants, in 2016, 70% of packs included at least one urban settlement within the expected home range. Moreover, the distance kept by wolves from the nearest urban settlement was more constrained in areas of recent recolonization (2005 and 2016) and in the mountainous altitudinal range, likely due to a reduced availability of territories. We showed that wolves tend to keep as far away as possible from humans, but they can also permanently occupy human-dominated land-scapes to cope with the lower availability of space induced by their remarkably successful recolonization. Our results shed light on an upcoming scenario for the conservation of large carnivores.


Introduction
Wolves Canis lupus experienced historically one of the strongest persecutions ever recorded on a large carnivore by mankind. They were completely wiped out from a large portion of their original distribution by a steady and intense hunting, trapping and poisoning, in Western Europe and in North America (Boitani, 2000;Ripple et al., 2014). Differently from other carnivores, direct persecution was the actual cause of wolves' decline and extinction in several areas, rather than environmental modifications or strong reduction in prey availability (Breitenmoser, 1998;Boitani, 2000). Indeed, differently from large cats or bears, wolves could rely on their well-known ecological plasticity to survive in many different environmental contexts, even in those strongly modified by humans (Blanco & Cort es, 2007;Eggermann et al., 2011;Llaneza, L opez-Bao, & Sazatornil, 2012;Ahmadi, L opez-Bao, & Kaboli, 2014;Kuijper et al., 2016). Thus, their extinction was, in most cases, the outcome of a deliberate and constant effort aimed at removing a cause of strong nuisance for livestock breeding and sometimes a danger for humans themselves (Cayuela, 2004;Treves et al., 2004;Kusak, Skrbin sek, & Huber, 2005;Bisi et al., 2007). Even nowadays, human-induced mortality is the main source of mortality for wolves in large portions of their range, mainly in Western countries with variable human population density ( Alvares, Pereira, & Petrucci-Fonseca, 2000;Carreira & Petrucci-Fonseca, 2000;Lovari et al., 2007;Treves et al., 2017;Musto et al., 2021). As reported for other large carnivores (Klees van Bommel et al., 2020), the proximity to urban settlements may be dangerous for wolves, which may reduce movements (Ferreiro-Arias and Llaneza, submitted) or approach human settlements only at night  to avoid contacts with humans. This is consistent with the history of the species in Western Europe and North America where wolves were eradicated from all environments with the exception of forested remote areas characterized by limited human presence and scarcity of urban settlements, eventually represented by scattered farms or tiny villages (Promberger & Hofer, 1994;Boitani & Ciucci, 1995) making the wolf a typical representative of ecological refugees (Kerley, Kowalczyk, & Cromsigt, 2012). This phenomenon biased the perception of wolf habitat preferences, convincing researchers that wolf was a typical forest-dwelling species with a strong aversion for all human-related infrastructures, like roads and cities, and human-modified environments such as agricultural landscapes (Ciucci et al., 1997;Theuerkauf, Rouys, & Jedrzejewski, 2003;Cayuela, 2004;Włodzimierz Je z drzejewski et al., 2005;Oakleaf et al., 2006;Lesmerises, Dussault, & St-Laurent, 2012).
However, the recent history of this species in Western Europe is characterized by a strong change in the aforementioned patterns, as wolves are impressively increasing (Chapron et al., 2014;Mech, 2017). The decrease in direct persecution by man, which was partially linked to the decline of free-ranging livestock breeding as a major cause of conflict with humans, played a relevant role in allowing population recovery (Treves et al., 2004;Bisi et al., 2007;Ripple et al., 2014), together with array of national and European laws protecting large carnivores (Boitani, 2000). Moreover, rewilded areas increased due to the abandonment of mountain and hilly areas by people that moved to urban areas (Navarro & Pereira, 2015) mainly after Second World War. This allowed natural vegetation to increase (Navarro & Pereira, 2015) resulting in a cascade effect on the abundance of wild ungulates (Apollonio, Andersen, & Putman, 2010), which in turn led their predators to recover their former abundance and distribution (Bruskotter & Shelby, 2010;Chapron et al., 2014;Gippoliti et al., 2018).
In the face of this quickly changing situation and because of the ongoing expansion of this species in Europe, increasing knowledge on wolf presence and site selection in anthropogenic contexts is currently needed to improve effectiveness of conservation and management of this species.
To investigate this issue, we considered the proximity of wolf packs to urban settlements to test whether wolves still actively avoid human disturbance by settling as far as possible from urban settlements or if, as an alternative, wolves are increasing their tolerance towards humans and do settle without avoiding towns and cities.
We used data from pack distribution of a densely wolfpopulated region in Europe, Tuscany (Central Italy), that shows a widespread and long-term monitored wolf population as well as three well-characterized environments (mountain, hill and plain) differing in both magnitude of human presence and availability of natural habitats. To test the hypothesis that wolves avoid areas with higher density of urban settlements, we formulated the following predictions: 1 packs distribution would be unequally subdivided among the three macro environments present in the region following the density of urban settlements that are less abundant in mountains and more abundant in plains; moreover, the process of recolonization would be gradual, from least human-inhabited environments to most human-dominated ones; 2 packs would not be randomly distributed with respect to urban settlements as they prefer occupying locations as far as possible from humans.
To test the hypothesis that pack distribution is constrained by both the human presence and the intra-specific competition for spatial availability, we predicted that: 3 the distance of packs from the nearest urban settlement would be constrained towards stable values by higher wolf density (i.e. across the recolonization process) and by environmental and human factors (i.e. across the three altitudinal ranges).

Study area
Tuscany is a region of Central Italy (Lat. 43°25 0 N; Long. 11°00 0 E) that extends for 18,513 km 2 excluding the islands (where the wolf is absent).
Tuscany habitats range from mountain to Mediterranean ones, as its elevation on sea level varies from 0 to 2054 m a.s.l. Woods are mostly composed of beech Fagus sylvatica and white fir Abies alba at higher elevation followed by deciduous oak Quercus sp., and chestnut Castanea sativa and Mediterranean woods characterized by holm oak Quercus ilex and domestic pine Pinus pinea. Altogether, woods occupy almost 9000 km 2 and are progressively increasing due to abandonment by farmers of vast areas once cultivated. About 52% of the region is occupied by forests and semi-natural environment, 38% by agricultural areas, 9% by human-dominated areas and 1% by water bodies and wetlands (Tuscany Region, 2022). The protected areas cover almost 10% of the regional territory, for a total area of 2,270 km 2 . Human activity, as well as the dense road network (292 km/100 km 2 ), is widespread throughout the territory with a greater urban concentration along the Arno valley. The average human density of the study area is 181 inhabitants/km 2 and specifically of 15 inhabitants/km 2 in the mountain areas, 59 inhabitants/km 2 in hill areas and 328 inhabitants/km 2 in plain areas (Tuscany Region, 2022).
Agriculture and livestock farming are still relevant to the region's economy, free-ranging livestock breeding is practised, with 401,151 sheep and 20,165 goats raised within the region (Berzi, 2018).
Wolves never got extinct in Tuscany. Even during their minimal historical distribution in the 1970s, wolves were reported along the Thyrrenian coast on the metalliferous hills (between the provinces of Pisa, Livorno and Grosseto) and along the Apennine in the Casentino and Mugello areas (northernmost part of the province of Arezzo and Florence) and in a very small area of Massa Carrara province (Cagnolaro et al., 1974). Boitani & Ciucci (1996) confirmed the presence of wolves in these areas, with an expansion to the neighbouring territories, which included a large part of the Apennine chain and the province of Grosseto. Subsequent investigations Apollonio, 2013Apollonio, , 2014Apollonio, , 2015Apollonio, , 2016 have been monitoring wolves in this region since year 2005. In this time, the distribution of wolves increased from mountains to hilly and flat areas ( Fig. 1) as wolf presence was reported in 29 and 220 municipalities in 1974 and 2016, respectively.

Data collection
We identified wolf pack locations during a survey conducted between 2014 and 2016 in the whole region, by means of a multiple-method approach (as recommended, inter alia, by Duchamp et al., 2011 andAusband et al., 2014).
We used a combination of camera trapping Mattioli et al., 2018), wolf-howling recording and sonogram analysis (Gazzola et al., 2002;Palacios et al., 2016); direct observation and filming, taking advantage of a previous knowledge of the spatial habits of most packs and a network of local collaborators (i.e. hunters, shepherds and other volunteers). Each pack known in previous years or reported by the citizen science network was investigated using the techniques mentioned above, to ensure presence and current location of the territorial pair and the presence of offspring (Table S1).
Camera trapping was conducted on a year-round basis to detect wolf presence within opportunistically selected locations and to ascertain eventual reproduction. Territorial pairs were recognized by observation of marking behaviours (Llaneza, Garc ıa, & L opez-Bao, 2014). In cases of uncertainties on pack identity, this approach was integrated by genetic analyses through sampling of dead wolves collected (Scandura, Iacolina, & Apollonio, 2011;Canu et al., 2017). A detailed description of the camera trapping methodology and wolf individual recognition is reported in Canu et al. (2017) and Mattioli et al. (2018).
Wolf-howling surveys were conducted yearly from June to October, focusing on the previously known or supposed pack's home sites (through the collaborative group's prior knowledge), according to the methodology described by Gazzola et al. (2002) and Passilongo et al. (2010). For discriminating different packs based on howling responses, we adopted the criteria described in Apollonio et al. (2004). Wolf howling provided information on pack reproduction by sonographic analysis of chorus howls Palacios et al., 2016).
We considered 'pack' each social unit constituted by, at least, a territorial pair. When reproduction was ascertained, we identified the pack location with the rendezvous site, that approximately coincided with the site where pups were detected (i.e. the camera trapping site they were filmed or the site they emitted the recorded chorus howl). Conversely, when reproduction was not ascertained, we considered as pack location the site with the maximum number of detections of adults (i.e. the camera trapping site where a certain pack was filmed more times).
To analyse the differences in wolf distribution across recolonized areas during different periods, we assigned each pack location to a recolonization step of the Tuscan wolf population, by using the occupancy data at municipal level in 1972 (Cagnolaro et al., 1974), 1996(Boitani & Ciucci, 1996 and 2006 . These historical occupancy data were based on records of wolf presence reported by citizens and stakeholders. Only in the case of the 2006 survey, wolf occupancy was further verified by means of wolf howling, snow-tracking and genetic analyses . Analogously to our approach, these methods were only used to verify wolf presence where it had been reported by the citizen science network. Thus, the higher reliability and precision of the more modern methods made the results of 2006  and 2016 (this study) more conservative than the previous surveys. Based on the municipality in which each pack location fell, we assigned to each pack the reference year of the earliest record of occupancy of its municipality. The year of our monitoring (2016) was then assigned to those packs located in areas not previously reported as occupied by wolves (recent recolonization).

Data analysis
Following our predictions, analyses were performed by a three steps approach.

Distribution of observed pack locations
Within the first step, we aimed to evaluate pack distribution among three macro-habitats (with different urban settlements presence; Table S2) defined at municipal level by altitudinal ranges (mountain, hill and plain) by comparing through a chi-squared test the observed frequency of pack locations in each habitat, with respect to a frequency distribution proportional to habitat availability in Tuscany. We considered as mountainous area those ranges higher than 600 m a.s.l.; as hills those located between 200 and 600 m a.s.l., and as plain those lower than 200 m a.s.l. Mountains covered 19% of the region, while hills and plains cover 46% and 35%, respectively (Tuscany Region, 2022). We then assigned each observed pack location to its altitudinal range. This procedure was performed both on the whole set of pack locations and separately for each recolonization step. The available surface for each recolonization step area was calculated by subtracting the surface with the stable presence of wolves of the previous steps from that of the whole region (Fig. 1).

Effect of urban settlements on pack locations
In the second step, we tested whether pack location distribution was affected by the distance from urban settlements by comparing observed pack locations with randomly generated ones. We used a GIS software (QGIS 3.10 A CORUÑA) to generate 240 random points (from now on control pack locations) distributed among altitudinal ranges in the same proportions as the observed pack locations. The optimal amount of control pack locations was defined as the number of points needed to stabilize the variance of the distance from urban settlements. The value was found by a visual inspection of a line plot relating distance variance to the number of locations considered. We classified each control pack location by altitudinal range and recolonization step areas, as already described for observed pack location. We then generated the variable 'presence', assuming 0 and 1 in control and observed pack locations, respectively. 'Presence' was considered as the response variable of a generalized linear model (GLM) with a Logit link function and the predictors described in Table 1 (see Table S5 for their mean and standard error in the observed and control pack locations).
All those independent variables were tested for collinearity. First, we calculated the Pearson correlation coefficient between all possible pairs within the predictor variables considered (Zuur, Ieno, & Smith, 2007). We considered as nonnegligible correlations those with an r coefficient higher than 0.6 or lower than À0.6 (Zuur et al., 2009a). The surface covered by the nearest urban settlement was collinear with its number of inhabitants (for both any size-and large urban settlements). As the settlement surface had lower significance and score in a Random Forest rank (varImpPlot function from the 'randomForest' R package), we retained as predictors the number of inhabitants of the nearest urban settlement of any size and that of the nearest large urban settlement.
Then, once the least significant variables were removed, we recalculated Person's correlation coefficient and repeated this process until no residual correlated pairs remained. Subsequently, we performed a multicollinearity test using the corvif function of the 'AED' package (Zuur et al., 2009b) in R software to confirm the absence of multicollinearity among the remaining variables. All VIF values were less than three (see Chapter 26 of Zuur, Ieno, & Smith, 2007) indicating that there was no multicollinearity in the variables tested. All the residual variables were entered as predictors in the GLM stepwise models with backward elimination procedure. Different models were evaluated by combining unrelated predictors in all additive and multiplicative (interactions) ways. The best model was then selected according to Akaike's information criterion (Burnham & Anderson, 2004;AIC;Burnham, Anderson, & Huyvaert, 2011). We evaluated the relative importance of each predictor included in the best model by computing the Akaike weights from all the models whose AIC differed less than two points from the best one (Burnham & Anderson, 2004;Symonds & Moussalli, 2011). In order to avoid the issues arising from model averaging (see Banner & Higgs, 2017;Dormann et al., 2018), only the results of the best model were discussed. This approach was particularly appropriate in this analysis since our main goal was to describe the possible effect of urban settlements on the observed pack locations and not to predict their distribution. Nonetheless, we still performed a model averaging analysis and included the results in Appendix (Table S6).

Constraints of the variability of pack location distance from urban settlements
The limits to the free choice of pack location were tested by a comparison of the variability of the distance from the nearest urban settlement. We evaluated whether the coefficient of The distance (m) from the nearest urban settlement of any size. We considered as urban settlement any cluster of houses with a number of residents equal or greater than that of the smallest municipality in Tuscany (65 inhabitants) Distance from the nearest large urban settlement (dLU) The distance (m) from the nearest large urban settlement. We considered as large urban settlements the biggest urban areas that together accounted for 60% of region inhabitants (ISTAT, 2017). The smallest settlement of this category had 5046 inhabitants (ISTAT, 2017) Altitude (A) The altitude (m a.s.l.) of the pack location Altitudinal range (Ar) The altitudinal range (categorical): mountains (>600 m a.s.l.), hills (>200 m and <600 m a.s.l.) and plain (<200 m a.s.l.) Recolonization step (Cs) The reference year (1972,1996,2006  variation (CV) significantly differed among the four recolonization step areas (as proxies of the number of wolves and their regional density) and the three altitudinal ranges (as example of three different human densities). We computed CV via bootstrapping (with 100 resampling) by using the Boot package (Canty & Ripley, 2014) in R. The estimated means were then compared by ANOVA and Tukey honest significant difference post hoc test (R Core Team, 2022) to find actual differences among paired groups.

Distribution of observed pack locations
During the monitoring conducted between 2014 and 2016, 110 pack locations were identified. Distribution of pack location sites in the three altitudinal ranges was not proportional to their availability (v 2 = 71.747; P < 0.001) showing a preference for the mountains, where 54 pack location sites were observed over 21 expected. On the hills, 45 pack locations were observed compared with 51 expected. Finally, in the plain we observed 11 pack location instead of the 38 expected (Fig. 2). Accordingly, in the different recolonization step areas the number of packs for each area differed significantly from expectations (Table S4). The number of packs on the mountains was always higher than expected, while those in the plain were always lower than expected. Conversely, the number of packs in the hilly area until 1996 was less than expected, while after 1996, it was higher than expected.
With respect to the recolonization process, packs initially tended to occupy mountains and hills and later started to occupy plains (Table S4).

Effect of urban settlements on pack location
The best logistic GLM explaining 'presence' probability included as predictors: (1) the distance from the nearest urban settlement of any size, (2) the number of inhabitants of the nearest urban settlement of any size and (3) that of the nearest large urban settlement, (4) the altitude and (5) the interaction between the altitudinal range and distance from the nearest urban settlement of any size (Table 2). Considering the set of models with DAIC <2 (Table 2), the Akaike weight of the selected predictors averaged 1, 1, 0.72, 0.78 and 1, respectively (Table 3). Conversely, the predictors not included in the best model had much lower weights: 0.30, 0.62 and 0.62 for the distance from the nearest large settlement, the population density of the nearest settlement of any size and that of the nearest large settlement, respectively. The effects described by the best model were consistent with the outcome of the model averaging (Table 3,  Table S6).
The probability of pack presence was positively affected by increasing the distance from human settlements and by increasing altitude, while the number of inhabitants of the nearest urban settlement had a negative effect (Table 3). The intensity of the effect of the distance from the nearest urban settlement was lower in the hills with respect to mountain and plain. Wolf packs were located at an average distance of 2590 AE 128 m (mean AE standard error) and 1961 AE 103 m from urban settlements of any size, for observed and control pack locations, respectively. Analogously, observed packs were located at an average distance of 10913 AE 598 m from large settlements, while control points were at 9802 AE 442 m from them ( Fig. S1 and Table S5). Observed pack locations were located near to settlements with 565 AE 85 inhabitants, while control points were located near settlements with 14492 AE 5352 inhabitants. The 93% of packs were located at less than 5 km from urban settlements of any size. As for large urban settlements, the 18% of the observed pack locations were within 5 km from them (Fig. 3).

Constraints of the variability of pack location distance from urban settlements
The variability of the distance from the nearest urban settlement was significantly lower in areas recolonized recently (2016) in comparison with those recolonized in previous years (F = 138.6; d.f. = 3/396; P < 0.001; Fig. 4a and Table S7).
The comparison of the CV among the three macro-habitats (identified by the three altitudinal ranges) showed a statistically significant difference (F = 23.61; d.f. = 2/297; P < 0.001), with the distance from the nearest urban settlement being significantly less variable in mountains than in the hilly and plain ranges of the Region (Fig. 4b and Table S7). dU = distance of the wolf pack location from the nearest urban settlement; dLU = distance of the wolf pack location from the nearest large urban settlement; pU = number of inhabitants of the nearest urban settlement; pLU = number of inhabitants of the nearest large urban settlement; A = altitude of the pack location; DpU = density population of the nearest urban settlement; DpLU = density population of the nearest large urban settlement; Ar = altitudinal range (for more details on each predictor see Table 1); df = degrees of freedom; R 2 = coefficient of determination, the proportion of the dependent variable variability predictable by the model; loglik = log likelihood; AIC = Akaike's information criterion; DAIC = the difference between AIC values for two nested models; weight = Akaike weight. À0.0003 0.0001 À2.005 0.044 dU = distance of wolf pack location from the nearest urban settlement; pU = number of inhabitants of the nearest urban settlement; pLU = number of inhabitants of the nearest large urban settlement; A = altitude of pack location; Ar = altitudinal range (for more details on each predictor see Table 1); Coefficient estimate = estimated b coefficient of the predictor within the best model; Akaike weights = average Akaike weights of each predictor among the models with DAIC <2 (Table 2); se = standard error of estimated coefficient; z = z-ratio; P = Pvalue.
Animal Conservation (

Discussion
Wolves preferred mountains and tended to avoid plains across all the recolonization steps, while hills were avoided during the first recolonization step and selected during later steps. Proximity with an urban settlement and the number of its inhabitants reduced the likelihood of a wolf pack location. The distance of wolves to the nearest urban settlement was less variable in areas more recently recolonized, where the density of packs was the highest, and the same was true on the mountains with respect to plains and hills.

Distribution of observed pack locations
The results were in accordance with our first prediction. Mountainous areas were more occupied by wolf packs in comparison with plain areas, although their availability was lower in all the four recolonization steps. Moreover, the development of recolonization suggested that wolves preferred mountains and, once the latter were mostly occupied, started to locate at lower altitudes as well. Interestingly, this effect concerned hill areas until 1996 and then, as the pack density grew, started to spread in plain as well. Such phenomenon was likely mediated by a density-dependent dispersal, with young wolves being forced to move to lower altitude, more anthropized areas in order to establish new territories. The occupation of plains being lower than expected (particularly evident during the earlier recolonization steps) may be accounted to bottom-up and/or top-down processes. The low prey densities characterizing plains during 70s-90s (Apollonio, Andersen, & Putman, 2010) could have actually limited the potential expansion of wolves. Additionally, the  potential of human-caused mortality to strongly hinder wolf population recovery (Quevedo et al., 2019) may indeed be expected to be higher in the areas more densely inhabited by humans (i.e. plains). Nonetheless, the recent overall increase in ungulate communities even in low-altitude environments (Apollonio, Andersen, & Putman, 2010) and the humancaused mortality rates being nowadays unrelated to human density in Tuscany (Musto et al., 2021) may suggest that human-dominated plains are gradually losing their limiting effect on wolf populations. The tendency to avoid humans may have delayed the recolonization of human-dominated landscapes, leading many researchers to define suitable wolf habitats only forested areas with a low human impact (Mladenoff et al., 1995;Massolo & Meriggi, 1998;Mladenoff & Sickley, 1998;Salvatori et al., 2002;Gehring & Potter, 2005;Potvin et al., 2005;Włodzimierz Je z drzejewski et al., 2005;Karlsson & Sj€ ostr€ om, 2007;Je z drzejewski et al., 2008;Ahmadi, L opez-Bao, & Kaboli, 2014), where they were in fact confined as ecological refugees (Kerley, Kowalczyk, & Cromsigt, 2012). On the contrary, our data show that after 1996, wolves were increasingly present in anthropized environments (Table S3), which can thus be considered as suitable habitats for wolf presence. It is worth noting that these results describe the actual distribution of stable packs (i.e. social units composed by at least two individuals owing and defending a territory) in relation to human presence. Since our analytical approach did not distinguish between packs with and without ascertained reproduction, further studies are needed to investigate possible differences in resource selection patterns between reproductive and not reproductive packs. We are also aware that our results may have been influenced by not considering the human density changes from 1972 to 2016. Nonetheless, during this period the human density slightly increased in plains but decreased in mountains (Reynaud et al., 2020). This may not only represent a further explanation for the fast recolonization of mountains, but also highlight that wolves expanded in lowland, human-dominated environments when the human density was growing. As habitat selection depends on consumer density and/or resource availability (Avgar, Betini, & Fryxell, 2020), the approach of wild boar to urban areas (Cahill, Llimona, & Gr acia, 2003;Podg orski et al., 2013;Stillfried et al., 2017;Banti et al., 2021) may have favoured wolf expansion in that areas being this species its major prey item in the region (Bassi et al., 2012(Bassi et al., , 2017.

Effect of urban settlements on pack location
Among the wolf packs identified in 2016, 70% included at least one small and 18% at least one large urban settlement, respectively, within a 5 km radius, that is the radius of a hypothetical home range of 85-110 km 2 (Ciucci et al., 1997;Corsi, Dupre, & Boitani, 1999;Apollonio et al., 2004;Mattioli et al., 2018).
Although our control points were not of ascertained absence, but rather of undetected presence, the distances of observed pack locations to the nearest urban settlement differed significantly from the control, in contrast with  but in accordance with our second prediction. The tendency of wolves to locate farther from the nearest urban settlement in comparison with control points is consistent with previous studies showing that human infrastructure distribution negatively affects the likelihood of packs locating in a certain area (Capitani et al., 2006;W. Je z drzejewski et al., 2008;Bassi et al., 2015).
Moreover, the effect of the distance from the nearest urban settlement was lighter in the hilly range. This is probably the outcome of the presence of many wolves in an intensively anthropized environment with a greater availability of refuge areas (small wooded or bushy patches) than in the plain; however, they were highly localized in a patchy landscape (Table S2). This may have forced the wolf to settle where these specific conditions were found, with limited care for human presence. Indeed, the hillside is the range where human density is relatively high (see Section 2.1 and Table S2), but also the second altitude range most occupied by wolves with respect to its availability (Fig. 2).
The number of inhabitants in the nearest urban settlement also significantly influenced wolves' choices on where to locate, indicating that wolves preferred the outskirts of urban settlements with low population numbers (see also Table S5).

Constraints of the variability of pack location distance from urban settlements
As predicted, the variability of the distance from the nearest urban settlement differed significantly both across areas recolonized in different periods and across macro-habitats (i.e. Altitudinal ranges). Between 1972 and 2016, there was a gradual decrease in the possibility of choosing suitable locations due to the increase in pack density. Since there was a gradual increase in the amount of surface occupied by wolves and thus in the number of wolves themselves along the considered period (Tables S2 and S3), this further supports our hypothesis. That is, as the number of wolves increased, the availability of suitable areas decreased, leading wolves to select suboptimal areas. In contrast to what was predicted for the three altitudinal ranges, we found the least variability of the distance from the nearest urban settlement in the mountain range, considered as the most suitable macro-habitat. This is likely due to a higher density of wolves in this range, which was in fact more occupied than expected (Fig. 2). Thus, packs, having to keep a safe distance from humans, are forced to rearrange their territories to cope also with other packs and prey availability (densitydependent habitat selection, as explained by O'Neil et al., 2020).

Conclusions
In conclusion, we showed that wolves preferred locating as far as possible from humans but that they occupied locations relatively close to urban settlements in a densely inhabited region, likely to cope with the intra-specific spatial competition. As a consequence, the presence of wolves, even if conditioned by the presence of man, is spreading into the most anthropized areas without a saturation point being foreseen. It should be noticed that in the plain, we found only 10% of the total pack locations observed; therefore, a further expansion is probable in face of the increase in wolf packs. Karlsson & Sj€ ostr€ om (2007) showed how attitudes and perceptions of wolves varied with the distance of wolves from respondents. Probably, the impression that wolves are now closer to urban settlements than in the 1990 s is also because there are more and more packs in lowland areas, which are also the most anthropized. The increase in pack distribution and the expansion of their range can raise the possibility of encounters between wolves and humans; therefore, this progressive closeness between human beings and wolves, which it is already a major issue for large carnivore conservation, may increase over time. This is true both in Tuscany, where there is a ubiquitous presence of wolves, and in other European areas that are being quickly recolonized by wolves, as this species is able to live both in densely populated plains and in almost abandoned mountains, confirming its adaptability (Muhly et al., 2019). In such context, a proper management of human-wolf conflicts is pivotal in ensuring the conservation of wolves and in preventing a new wave of persecution of the species (Lute et al., 2018). Since risk-enhancing human behaviours are often the main cause of human-predator conflicts (Penteriani et al., 2016), where wolves are recovering after decades of absence it is necessary to encourage an appropriate behaviour of citizens, including a proper management of possible food items as garbage (Mohammadi et al., 2019) and domestic pets (Bassi et al., 2021).

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Table S1. Detailed effort and regional results of each pack known in years or reported by the citizen science network integrated survey. Table S2. Indicators of human presence within the considered altitudinal ranges. Table S3. Indicators of human presence within the areas available during each re-colonization step (see the text for more details). Table S4. Chi-square test results of the comparison between expected and observed wolf pack location frequency distribution in the three altitudinal ranges and in the four recolonization step areas. Table S5. Means and standard errors of all continuous variables used in the Generalized Linear Models used on wolf pack location. Table S6. Results of the model averaging among the set of equivalent Generalized Linear Models (DAIC <2) on wolf pack location. Table S7. Post hoc tests between the Coefficients of Variation of the distance from the nearest urban settlement, among the three different altitudinal ranges and the four different re-colonization steps. Figure S1. Violin plot representation of the distances of wolf pack locations from the nearest urban settlement; in black the observed wolf pack location and in grey the control wolf pack location.