En route to the North: modelling crested porcupine habitat suitability and dispersal flows across a highly anthropized area in northern Italy

The crested porcupine (Hystrix cristata) underwent a rapid and widespread range expansion in Italy. Nowadays the species is moving towards the northernmost regions of the country and its occurrence is increasing in the highly anthropized Po Plain. Our objectives were to evaluate the suitability of the Po Plain for the species, as well as to identify dispersal corridors connecting the northern Apennines occurrence areas and the Prealps. We modelled the species home-range scale habitat suitability based on an ensemble modelling approach. Additionally, a habitat suitability prediction carried out at a finer scale was used to parametrize the landscape resistance, based on which we modelled the potential dispersal corridors for the species using a factorial least-cost path approach. The ensemble prediction estimated a potential occurrence of the crested porcupine in 27.4% of the study area. The species occurrence probability was mainly driven by the distribution of extensive cultivations, woodlands and shrublands, and water courses and by the annual mean temperature. Conversely, the movements of the species resulted mainly sustained by woodlands and shrublands and highly hindered by simple arable lands and rice paddies. The connectivity prediction showed that three main dispersal routes are likely to connect crested porcupine occurrence areas in the northern Apennines to currently unoccupied but highly suitable areas in the Prealps. The study allowed us to identify the areas in the Prealps with the highest probability to be colonized by the crested porcupine in the near future and provided important insights for the conservation of a strictly protected species in a human-dominated landscape.


Introduction
The crested porcupine (Hystrix cristata) is a large rodent (up to 12-15 kg; Mori and Lovari 2014) native of North Africa which inhabits the Italian peninsula possibly since the early medieval times (Bertolino et al. 2015;Trucchi et al. 2016). Like many other generalist mammal species, during the last decades, the crested porcupine underwent a rapid and widespread range expansion in Italy: since the 1970s, the species almost doubled its original range once confined in central and southern Italy (Mori et al. , 2021. The expansion of the species can be mainly explained by legal protection [the crested porcupine is listed in Berne Convention (all. II), Habitat Directive (all. IV) and National Italian Law n. 157/1992 (strictly protected species)], and by environmental changes. In particular, the woodland expansion due to the abandonment of pastures and other traditional types of agriculture on mountain and hilly areas (Falcucci et al. 2007) greatly drove the successful colonization of new areas across the Apennine ridge in both northern (Mori et al. , 2018 and southern (Mori et al. 2021) Italian regions. Indeed, the species may adapt to several habitat types, from woodlands to suburban areas passing through agroecosystems (Mori et al. 2014a;Lovari et al. 2017), but it always requires dense vegetation cover for denning and feeding (Monetti et al. 2005;Mori et al. 2014aMori et al. , b, 2017Lovari et al. 2017). Recently, different studies pointed out the importance of changes in climatic factors such as temperature, drought and precipitation regimes, in promoting the species range expansion (Mori et al. 2018(Mori et al. , 2021. Specifically, during the 1 3 last decades, climate changes caused a decrease in the number of days of snow cover in mountain areas leading to the colonization of higher elevations by shrublands and woodlands (Gehrig-Fasel et al. 2007;Bani et al. 2019) triggering a plausible increase of the suitability of areas once considered inhospitable for the species (i.e., high mountain areas covered in snow and ice for long periods during the year).
Predictive distribution models developed for the crested porcupine in Italy (e.g., Mori et al. 2018;Milanesi et al. 2020) generally limited the suitable range for the species south of the Po Plain, the intensively cultivated lowland area interposed between the Apennines and the Alps (but see Mori et al. 2013). The colonization of the Po Plain and the northernmost western region of Italy is expected within the next 30-40 years under severe scenarios of climate changes (Mori et al. 2018). Despite these long-term predictions, crested porcupine occasional observations are rapidly increasing in the Po Plain area: given the peculiar morphology and the easily identifiable signs of presence, the expansion of the species in new areas is frequently documented on social networks, local newspapers, and citizen-science projects.
A recent study (Mori et al. 2021) predicted some areas of the Po Plain and of the Prealps (the mountain belt that borders the southern limit of the Alpine chain) as suitable for the crested porcupine. The Po Plain is characterized by the presence of scattered natural (e.g., riparian vegetation, broad-leaved woodlands, and shrublands) and semi-natural (e.g., hedgerows and arboreal cultivations) vegetated patches, which can represent a sub-optimal habitat supporting the occurrence of many forest-dwelling mammal species (Balestrieri et al. 2015(Balestrieri et al. , 2016Chiatante et al. 2017;Dondina et al. 2019). These landscape features can also mitigate the effects of fragmentation by increasing connectivity and supporting long-range dispersal movements of individuals, as observed for many mammal species both in other geographic contexts [e.g., carnivores (Hilty and Merenlender 2004;Šálek et al. 2009;Virgós 2001) and rodents (Bani et al. 2018;Silva and Prince 2008)] and in the Po Plain [e.g., the hazel dormouse Muscardinus avellanarius (Dondina et al. 2018a); the European badger Meles meles and the roe deer Capreolus capreolus (Dondina et al. 2019); the Italian wolf Canis lupus italicus ]. The possibility of crossing the Po Plain by the crested porcupine would have the important consequence of leading the species towards the Prealps, an area potentially suitable for colonization (Mori et al. 2021), where occurrences start to be detected.
Predicting areas of plausible future colonization is crucial for the conservation of wildlife, especially for species directly persecuted by humans, such as the crested porcupine, that, despite the high level of legal protection, is still considerably poached in Italy (Amori et al. 2008;Laurenzi et al. 2016;Lovari et al. 2017; but see also Cerri et al. 2017). Therefore, the aims of this study were (1) to predict the potential species distribution (i.e., habitat suitability) in the Po Plain and in the Prealps, (2) to model potential connections through the Po Plain able to support the species dispersal flows from the northern Apennines to the Prealps, (3) to identify the areas of the Prealps with the highest probability to be colonized by the species in the near future.
Although multiple future scenarios of expansion for this species have been already developed for the whole country (Mori et al. 2018(Mori et al. , 2021Milanesi et al. 2020), focusing on a particular area can provide important details useful for enhancing the effectiveness of conservation actions. In other words, species protection and mitigation measures aimed to reduce human-wildlife conflicts should include the peculiarities of the landscape at the regional scale. For these reasons, we focused on habitat suitability and ecological connectivity for the crested porcupine of the highly anthropized Po Plain, a key area for the species expansion towards the northernmost regions of Italy. We hypothesized that the expansion of the crested porcupine in the Po Plain would be driven by the presence of natural vegetation cover (i.e., woodlands and shrublands), mainly restricted along river and stream banks, and by the presence of semi-natural elements typical of extensive cultivated areas. Accordingly, we expected that the Prealps, with their extended natural woodland cover, would offer large and continuous areas suitable for the persistence of the species.

Study area
The study area encompassed the Po Plain as well as the area of the northern Apennines and the Prealps under 800 m a.s.l. (Fig. 1); it covers about 65,000 km 2 within the Italian regions of Piedmont, Lombardy, Veneto, Friuli-Venezia Giulia, and Emilia Romagna (Supplementary material: S1- Figure S1.1).
The conformation of the Po Plain, closed between the Alps and the Apennines, causes high levels of relative humidity throughout the year. Thus, the climate of the study area can be considered humid continental with cool-humid foggy winter (mean temperatures during coldest months: 0-5 °C) and warm-moist summer (mean temperatures during hottest months: 22-25 °C).

Habitat suitability analysis
We removed from the dataset all the occurrence data related to carcasses because they likely belonged to dispersal individuals (they were often detected along roads away from areas permanently occupied by the species) and the environmental preferences which drive habitat selection are known to be significantly different between dispersers and resident individuals (Mateo-Sánchez et al. 2015a;Dondina et al. 2018a). We checked the absence of spatial autocorrelation of data using Moran's I test with 999 permutations (Cliff and Ord 1981).
To assess the habitat suitability of the study area for the crested porcupine, we statistically compared the environmental characteristics of locations used by the species with those of random locations (Manly 2010) using an ensemble modelling approach.
A set of pseudo-absence locations were randomly generated across the study area adjusted to overcame the sampling bias (Phillips et al. 2009): we produced a probability density map using a kernel density estimator based on the occurrence locations and we used this bias map to generate the Fig. 1 The study area (red line) pseudo-absence locations (Clements et al. 2012;Fourcade et al. 2014). The number of pseudo-absence locations was set equals to the number of locations used by the species to maximize the pooled performance of the set of models used in the ensemble modelling procedure (see below) (Barbet-Massin et al. 2012).
According to the published literature (Mori et al. , 2018(Mori et al. , 2021, we considered 13 covariates describing both the land-cover composition and the topographic and climatic characteristics at each presence and pseudo-absence location. The land-cover variables were measured within a circular buffer of 2.5 km designed around each presence and pseudo-absence location (home-range scale); this spatial resolution is considered optimal for the species since crested porcupines range for over 3-5 km per night to reach the feeding areas (Mori et al. 2014a). The land-cover composition within the buffer was defined by measuring the fractional cover of eight land cover-based variables: urban areas, simple arable lands, rice paddies, extensive cultivations [i.e. permanent crops (vineyards and fruit orchards), complex cultivations (small cultivated land parcels with different cultivation types), and cultivations with natural elements (areas principally occupied by agriculture interspersed with significant natural or semi-natural areas)], meadows and pastures, broad-leaved and mixed woodlands and shrublands, coniferous woodlands, and water courses. We jointly considered woodlands and shrublands as other studies carried out on habitat use by the crested porcupine highlighted the importance of both these habitats (Mori et al. 2018). The land use cartography adopted to calculate the fractional cover of the land-cover variables was the Corine Land Cover 2018 (http:// www. coper nicus. eu) with a spatial resolution of 100 m. We considered two topographic variables derived from the digital elevation model with a spatial resolution of 20 m (http:// www. sinan et. ispra mbien te. it), which were altitude and aspect. As climatic variables, we considered the annual mean temperature (BIO1), the temperature seasonality (BIO4) and the annual precipitation (BIO12). The climatic maps were obtained from the World-Clim dataset (http:// www. world clim. org) for the current period (1970-2000 average) at a 30-s resolution (≈ 1 km).
We calculated the pairwise Pearson's correlation coefficient between covariates and we removed three covariates (altitude, temperature seasonality and annual precipitation) from the analysis because they were highly correlated to one or more other covariates (correlation coefficients > 0.75).
To investigate the species habitat suitability, we developed an ensemble model using the BIOMOD2 package (Thuiller et al. 2009) in R version 4.0.3 (R Core Team 2020). The ensemble model was implemented using a regression method (GLM), a classification method (CTA), and two machine-learning methods (GBM and RF). For each model we performed 10 runs. In each run, we randomly selected 75% of the data as a training set, while the remaining 25% as a test set. The discriminating power of each model was evaluated through the area under the curve (AUC) of the receiver operating characteristic (ROC; Swets 1988) and the true skill statistic (TSS). The ensemble prediction was developed by integrating the output of the single models adopting a weighted-averaging approach which weighted the single models according to their discriminating power (Thuiller et al. 2009;Rodríguez-Soto et al. 2011). The performance of the ensemble model was evaluated using the AUC of the ROC (Scherrer et al. 2019;Fernandes et al. 2019;Mohammadi et al. 2019).
We finally used the output of the ensemble model to produce an ensemble predictive map (resolution of 100 m) for the whole study area showing a gradient of suitability ranging from 0 to 1000. To carry out the prediction, we used 100-m resolution maps of aspect and annual mean temperature. Additionally, for each land cover-based variable, we produced a raster map by assigning to each pixel (sized 100 m × 100 m) the fractional cover of the variable in a 2.5km buffer area around the pixel.

Dispersal flows modelling
To model the predicted crested porcupine dispersal flow from the northern Apennines to the Prealps, we used the universal corridor network simulator software UNICOR (Landguth et al. 2012). UNICOR adopts a modified Dijkstra's algorithm (Dijkstra 1959) within a factorial least-cost path approach to overcome the limitation of the classic leastcost path analysis associated with the number of sources and target locations and to generate a synoptic measure of landscape connectivity .
Two data layers were used as input: (1) a resistance surface depicting the cost of movement for the crested porcupine throughout the landscape, and (2) a point file containing the geographic coordinate of each source location.
To obtain the resistance surface, we developed a fine scale ensemble model. In this model, we considered as presence and pseudo-absence locations the same locations considered in the habitat suitability ensemble model developed at the home-range scale, while we considered as covariates the land-cover variables measured within a circular buffer of 250 m designed around each presence and pseudoabsence location. Specifically, we considered as variables the fractional cover of seven land-cover types (urban areas, simple arable lands, rice paddies, extensive cultivations, broad-leaved and mixed woodlands, shrublands, and water courses). The fractional cover of meadows and pastures and coniferous woodlands were excluded from the analysis because they were poorly represented at the fine scale (Fattorini et al. 2014). For this analysis, we separately considered broad-leaved and mixed woodlands and shrublands because there is no published study regarding the role that these habitats may play in sustaining the crested porcupine dispersal. We considered a fine spatial scale to modelling dispersal flows because habitat suitability models have proved to be a useful alternative to movement data to parameterize the landscape resistance to animal movement provided that environmental variables are calculated at a finer spatial scale compared to the scale generally used to parametrize habitat suitability for a species, i.e. the home-range scale (Ziółkowska et al. 2016;Dondina et al. 2020;Torretta et al. 2020). Dispersers typically depend on the availability of local resources during their movements (Zeller et al. 2012), while resident individuals usually have broadscale ecological requirements (Mateo Sánchez et al. 2014). The ensemble modelling procedure was the same as that described above for the habitat suitability analysis performed at the homerange scale. Based on the output of the ensemble model, we carried out a fine-scale prediction. Specifically, for each land cover-based variable, we produced a raster map by assigning to each pixel (sized 100 m × 100 m) the fractional cover of the class in a 250-m buffer area around the pixel. To obtain the landscape resistance map, we converted the suitability scores derived from the prediction of the ensemble model developed at the fine spatial scale using a negative exponential function (Mateo-Sánchez et al. 2015a, b).
To obtain the source locations, we considered the ensemble predictive map developed through the ensemble model performed at the home-range spatial scale (2.5 km), i.e., the suitability map obtained for the whole study area through the habitat suitability analysis. We converted the ensemble prediction continuous map into a binary one (suitable/unsuitable) considering a threshold value estimated by maximizing the AUC (Thuiller et al. 2009); values higher and lower than the threshold represented suitable and unsuitable areas, respectively. We then superimposed a grid of 5 km × 5 km cells [(the optimal spatial resolution for the species; Mori et al. (2018)] to this map and selected the cells covered by suitable habitat for at least 75% of their extent. The centroids of these cells were used as source locations for connectivity modelling. While in the Apennines and lowland area we retained all source locations, they were fictitiously rarefied in the Prealps (forcing a minimum distance of 10 km) to simulate a very sporadic presence of the crested porcupine in this area .
To represent the potential of the entire study area in supporting the crested porcupine movements, the estimated least-cost paths were smoothed as a function of cumulative cost using a kernel density function and a Gaussian probability distribution (Shahnaseri et al. 2019). The smoothed least-cost paths were then overlapped to obtain a map of corridor density (Cushman et al. 2013a) where each pixel had a value corresponding to the frequency of least-cost paths crossing that pixel. To model the potential dispersal flows for the crested porcupine in the whole study area, we set no threshold for maximum dispersal distance when modelling the least-cost paths (e.g., Cushman et al. 2013b;Shahnaseri et al. 2019).
The accuracy of the estimated corridor density was evaluated through the Boyce index (Boyce et al. 2002) by splitting the pixel values of the corridor density map into ten equal continuous classes using quantiles. The Boyce index varies between − 1 (i.e., counter prediction) and + 1 (i.e., consistent prediction; Hirzel et al. 2006). As validation dataset, we used all the occurrence data related to carcasses (n = 44) because they likely belong to dispersal individuals which preferentially move along areas characterized by high corridor density values.

Results
We used a total of 122 species locations to run the habitat suitability and connectivity analyses. The data were not affected by spatial autocorrelation (Moran's I = − 0.004; P = 0.777).

Habitat suitability analysis
Among the single models developed within the ensemble modelling procedure performed at the home-range spatial scale, the GLM had the lowest performance in predicting habitat suitability ( The ensemble model performed at the home-range scale showed an excellent discriminating power (AUC = 0.952). The most important variables in predicting the crested porcupine occurrence were the fractional cover of extensive cultivations, broad-leaved and mixed woodlands and shrublands, water courses and the annual mean temperature (Table 1).
Considering the resulting threshold value (547), the ensemble prediction estimated the occurrence of the crested porcupine in 27.4% of the study area. The prediction map showed that the most suitable areas for the species were located in the hilly and mountain areas of both the Apennines (in the southern part of the study area) and the Prealps (in the northern part of the study area). Highly suitable areas were also located along the main rivers (e.g., Po and Ticino rivers) and some minor affluents. Other two particularly suitable areas for the species were the lowland area that borders the Apennine hilly area in the southern part of the Lombardy region and the hilly areas of the Piedmont region.
Contrariwise, the most unsuitable areas for the crested porcupine were located in the intensively cultivated plain, where the lowest suitability values are found in areas characterized by the presence of rice paddies, i.e., the areas east and west of the Ticino River, and a small area south of the Po River in the easternmost part of the study area (Fig. 2).

Dispersal flow modelling
The ensemble model developed at the fine spatial scale showed a good discriminating power (AUC = 0.809). The model highlighted that the most permeable land-cover types for the movement of the crested porcupine were broadleaved and mixed woodlands and, second, shrublands and extensive cultivations, while the most hostile land-cover types are simple arable lands and rice paddies (see Supplementary material: S4 for more details regarding the output of both the single models and the ensemble model performed at the fine spatial scale).
The dispersal flows analysis was performed on 149 source locations (62 in the Apennines; 42 in the lowland and hilly areas in the Piedmont region; 45 in the Prealps). The cumulative density of the buffered least-cost paths across the study area showed the potential connections (i.e., dispersal flows) between the Apennines and the Prealps (Fig. 3).
As expected, the suitable areas for the species in the hilly and mountain regions of the Apennines were all well connected to each other. The Apennines area was also well connected to the central part of the Po River by a well-structured network of potential dispersal flows characterized by a high least-cost path density. The Po River corresponded to the  there were three main routes to reach the Prealps. The first route was composed of two parallel flows that radiated to the western part of the study area. The northern flow followed the Po River reaching the Prealps area located northwest of Turin, while the southern ran parallel to the first about 10 km away and reached the Prealps area located south-west of Turin. The second main route to the Prealps perfectly followed the course of the Ticino River, from the city of Pavia to Lake Maggiore. The third route was composed of two parallel flows about 15 km apart which originated in the central region of the Po River and moved toward northeast reaching the Prealps north to the cities of Verona and Vicenza. This area was well connected to the eastern Prealps by relatively strong connections. Conversely, the south-easternmost part of the study area was connected to the species population inhabiting the Apennines areas by a unique quite weak connection. Boyce index suggested that the estimated least-cost paths density was highly robust (Spearman rank correlation coefficient = 0.81), thus able to predict the most suitable dispersal routes for the crested porcupine between the Apennines and the Prealps (Fig. 4). The predicted-to-expected ratio had an increasing trend from the first to the second to last classinterval, showing a lack of presences (i.e., carcasses) within the last class-interval (i.e., the strongest paths).

Discussion
The crested porcupine is one of the most studied mammal species in Italy; among the most recent published researches, particular attention was given to the observed rapid spread of the species throughout Italy (Mori et al. , 2021 and  (Mori et al. 2018(Mori et al. , 2021Milanesi et al. 2020). In this study, we provided the first assessment of habitat suitability and ecological connectivity for the crested porcupine in the highly anthropized Po Plain, a key area for the species expansion towards the northernmost regions of Italy.
Similarly to other species (e.g., the Italian wolf; Fabbri et al. 2007), the crested porcupine mainly exploited the Apennine chain for its range expansion, taking advantage of the continuous woodland cover, mild climate, and gentle slopes. Once crossed the northern Apennines ridge, the crested porcupine faced a completely different landscape, characterized by wide unsuitable areas for the species, as predicted by the home-range scale ensemble model. Specifically, the intensive agricultural lands that dominate the Po Plain cannot be permanently exploited by the species. To persist throughout the year, the crested porcupine relies on the presence of densely vegetated areas (Monetti et al. 2005;Lovari et al. 2013) and this makes the mono-specific cereals cultivations, completely bare or covered with little vegetation during the winter, scarcely able to support the species over time.
In our study area, the most suitable locations for the species were the hilly and mountain areas of the Apennines, where the crested porcupine presence is stable and regular, and the Prealps, where the presence of the species is generally still occasional and reproduction was rarely documented (Spada et al. 2008). In these areas, woodlands and shrublands, one of the most important explanatory variables in predicting the occurrence of the crested porcupine in the study area, have a continuous and dense cover. As stated above, the importance of woodlands and shrublands for the species has been widely recognized: the crested porcupine needs these habitats, even in the form of small patches (e.g., bramble bushes), for denning and feeding. In detail, the species prefers areas covered with thick vegetation to increase protection and stability of the dens and it mainly feeds on underground storage vegetal organs and wild fruits, which are abundant in woodlands and shrublands (Monetti et al. 2005;Mori et al. 2014aMori et al. , b, 2017Lovari et al. 2017;Mori and Fattorini 2019). Other suitable areas for the species were located along several water courses within the Po Plain. Here the presence of the species is still scattered and restricted; the most suitable water courses are those characterized by thick natural and semi-natural vegetation cover along banks, for example the Po and the Ticino (Lombardy-Piedmont regions) rivers, but also the Adige and Piave (Veneto region) and Tagliamento (Friuli-Venezia Giulia region) rivers, and some of their affluents like the Trebbia, Nure and Taro streams in Emilia Romagna. Conversely, the areas surrounding water courses characterized by the lack of dense arboreal or shrub-like cover along banks (e.g., the Adda River) resulted inhospitable for the species.
Another important explanatory variable positively related to the occurrence of the crested porcupine was the fractional cover of extensive cultivations. The areas characterized by this land-cover type are usually highly heterogeneous (e.g., small land parcels with different cultivation types) and frequently interspersed with natural or semi-natural vegetated areas. This landscape context represents a suitable environment for the species because it provides seasonal abundant food resources. For example, permanent crops (i.e., vineyards and fruit orchards) are typical of extensive cultivations areas and provide concentrated, predictable, and easily accessible sources of food. Indeed, fruits growing on the lowest branches and those fallen on the ground represent an important component of the diet of crested porcupine Lovari et al. 2017). Accordingly, those areas particularly suitable for permanent crops [i.e., the hilly areas of Pavia province in Lombardy and the Prealps in Veneto (vineyards) and the hilly areas of Cuneo province in Piedmont (fruit orchards)] resulted highly suitable for the occurrence of the species.
The last important variable in driving species occurrence probability at the home-range scale was the annual mean temperature. Specifically, the estimated relationship with the species occurrence was mainly negative probably because suitable land covers for the species (i.e., woodlands and shrublands and extensive cultivations) are mainly located in mountain and hilly areas, where temperatures are usually lower than those characterizing the lowland areas.
In this study, we developed the first analysis of ecological connectivity for the species in Italy. The fine scale ensemble model showed that the land-cover types likely used by the species during dispersal movements are the same used by resident individuals: woodlands, shrublands and extensive cultivations. Nevertheless, the home-range ensemble model revealed that woodlands and shrublands and extensive cultivations are equally suitable for the species, which jointly exploits these habitats to find shelter and food. Conversely, woodlands, and second shrublands, resulted to be more important for species movements compared to extensive cultivations. This result may indicate that movement behaviour of the crested porcupine may be likely mainly driven by the need to find shelter, rather than by the availability of food resources. The crested porcupine is a generalist species that feed on a variety of vegetal food resources in relation to their availability (Lovari et al. 2017); it is thus plausible that the availability of food is not the limiting factor for the species during dispersal movements. The fine scale ensemble model also revealed that species movements are highly hindered by simple arable lands and rice paddies. These land-cover types have no element able to provide shelter to the species. Specifically, rice paddies, annually flooded during spring and summer, represent a very hostile matrix feature for movements of ground-dwelling mammals (Dondina et al. 2016).
As expected, the predicted suitable areas in the Apennines resulted all well connected to each other by strong dispersal connections. The obtained results suggested that more and more individuals may move from the hilly areas of the northern Apennines to the banks of the Po River following a complex network of potential dispersal flows. Once reached the Po River individuals coming from the Apennines can easily move westward along the river and reach unoccupied areas in hilly and mountain areas of western Piedmont region. Individuals can also move north-westward through the dispersal route which follows the Ticino River or, on the other side, move north-eastward following the suitable routes towards Veneto and Friuli-Venezia Giulia regions. The crested porcupine will plausibly follow one of these three main alternative ways to reach the Prealps; we thus expect that the Prealps areas with the highest probability to be colonized by the crested porcupine in the near future correspond to the ending points of these routes, i.e. the hilly and mountain areas around Turin in Piedmont region, the areas around the Lake Maggiore at the boundaries between Piedmont and Lombardy regions, and the hilly area north to the city of Verona and Vicenza at the boundaries between the Veneto and the Trentino regions. We are quite confident regarding the likelihood of the identified dispersal routes based on the result obtained through the validation analysis carried out by means of an external dataset.
Despite the existence of multiple potential dispersal routes, considering both the landscape characteristics of the Po Plain (i.e., the limited and fragmented woodland cover) and the dispersal ability of the species (Mori and Fattorini 2019), we expect that the colonization process of the Prealps by the crested porcupine will probably occur in multiple generations. This prediction gains even more reliability considering that the Po Plain is the most densely populated area in Italy (average population density of 450 inhabitants/km 2 ) and it is crossed by several high traffic roads and highways, which could significantly slow the species expansion process through the phenomenon of road killing.
On the other hand, human-mediated actions, as the illegal introduction occurred in the Varese province possibly around 2007 , may furtherly support the species range expansion across the Po Plain.
Notwithstanding the limited extension of the study area and, consequently, of the number of species occurrences considered for developing the habitat suitability model and the potential dispersal corridors, compared with other nation-wide researches (Mori et al. , 2018(Mori et al. , 2021Milanesi et al. 2020), we obtained highly reliable results, which provide important insights for the conservation of the strictly protected crested porcupine in northern Italy, especially in light of the ongoing expansion process across the Po River Plain. In this human-dominated area, where the fragmentation of suitable habitats increases species vulnerability, conservation measures should be aimed at strengthening the effectiveness of the existing ecological network (McRae et al. 2008;Dondina et al. 2018b) to support the colonization of unoccupied but highly suitable areas. Moreover, the prediction of the Prealps areas with the highest probability to be colonized by the species in the near future can help to develop future conservation actions aimed at mitigating the conflict between the species and human activities typical of areas of new colonization. Indeed, the tolerance of farmers towards the crested porcupine is very low. Despite the fact that the actual damage to agricultural activities is limited, this large rodent is still considered one of the worst agricultural pests and is still extensively poached in Italy (Laurenzi et al. 2016).
Author contributions ET and OD conceived the study, conducted all data analyses and drafted the manuscript, while VO and LB contributed to the data analyses and reviewed the manuscript. ET and SM collected most of the data. All authors read and approved the final manuscript.
Funding Open access funding provided by Università degli Studi di Pavia within the CRUI-CARE Agreement. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.
Availability of data and materials All data analysed during this study are included in the supplementary materials to this published article. The datasets used and analysed during the current study are available from the corresponding author on reasonable request.
Code availability Not applicable.

Conflict of interest
The authors declare that they have no competing interests.
Ethics approval This study complies with the guidelines or rules for animal care and use for scientific research.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will 1 3 need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.