The Eurasian beaver range expansion reveals uneven future trends and possible conservation issues: an European assessment

The Eurasian beaver is a keystone species and landscape-capable ecosystem engineer, which went close to extinction until the 19th century. Recently, thanks to legal protection and reintroduction programs, the species has recolonized much of its past range. However, in some countries this process did not occur. Objectives. Our objective is to model the potential distribution of the Eurasian beaver for current and future conditions, on a continental scale, at river and sub-basin level. We focus on the protected areas of Italy and Portugal for possible reintroductions. Methods. The study area is Europe, with a subset focusing on Italy and Portugal. We produce species distribution models for current and future conditions using climate change scenarios and predicting changes in river flow, including topographic and human disturbance variables. We then deepen suitability-related issues within Italian and Portuguese protected areas. Results. We find that the Eurasian beaver current suitability is comparable to its known distribution, although some potentially-suitable spots occur in Italy (where there are two occurrences), while the Iberian Peninsula and the Balkan countries host scattered suitable spots. Future scenarios predict a general lowering of suitability in Central and Northern Europe. Portuguese protected areas generally host unsuitable territories, while the Italian ones have reported a tangled scenario, depending on the biogeographical sector. Conclusions. Our results may support the large-scale management of the beaver, both for countries already hosting this species and those planning a reintroduction. The framework used may be applied to other species, and for different topics, from biogeography to conservation.


Introduction
A global decline in biodiversity is underway: an average demographic decline of 68% was found in monitored populations of vertebrates between 1970 and 2016 (Living Planet Index 2022), including more than 1,300 mammal species (over the about 6,000 assessed) considered as currently threatened (IUCN 2021).
Among the biogeographic regions, the Palearctic shows the smallest decrease, with an average decline of the monitored populations of 24%, owing to several conservation measures which avoided severe declines (Pereira et al. 2010;Living Planet Index 2022). Indeed, in the last 20 years, this area hosted several species' comeback, thanks to focused management actions and international regulations: reintroduction programs succeeded, such as for the European bison (Bison bonasus), the red deer (Cervus elaphus), the alpine ibex (Capra ibex) (Saether et al. 2002;Pucek et al. 2004;Apollonio et al. 2010) and for some European carnivores (Chapron et al. 2014). When properly planned, reintroductions have thus proven to give positive results (Frey and Walter 1989;Johnson and Cushman 2007;Linnell et al. 2009), especially when applied to keystone or umbrella species (Linnell et al. 2009;Ripple and Beschta 2012). Among these, an example of the successful reintroduction of a keystone species is the Eurasian beaver (Castor fiber) (Gaywood 2018). Reintroductions of this species have taken place in 25 European countries, most of which have been successful (Nolet & Rosell, 1998;Halley & Rosell, 2002). However, in some cases reintroductions have failed due to low numbers of individuals released, or mortality caused by the stress of moving individuals (Myrberget 1967;Nolet & Rosell, 1998;Dewas et al., 2012). Previous analysis revealed that of 87 past reintroductions, 53% were successful (Macdonald et al. 1995).
The Eurasian beaver, an herbivore with a semi-aquatic lifestyle, is the second largest rodent in the world and is recognized as an ecosystem engineer due to its ability to modify rivers and streams through the creation of dams (Wright et al. 2002;Rosell et al. 2005). It is also defined as a keystone species (Mills et al. 1993), considering its role in maintaining the stability and functionality of ecological communities it lives in, as well as its unmatched influences on the ecosystems (Mills et al. 1993;Rosell et al. 2005).
This species was widely distributed in Palearctic, but populations started declining during the Middle Ages, up to the 19th century, mainly due to hunting for fur, meat and castoreum (a secretion produced for scent marking and historically used by humans for medical purposes) (Halley and Rosell 2002;Halley et al. 2021). Some estimates account for only 1200 individuals surviving, divided into eight isolated populations in the Palearctic, of which five of them in its Westernmost part (Nolet and Rosell 1998;Halley et al. 2021). In recent times, legal protection and reintroductions (started since 1922), together with natural recolonization, allowing the species to recover (Halley and Rosell 2002;Halley et al. 2021).
These reintroductions played an important role in restoring populations of the Eurasian beaver within its old range, but they did not involve all the European countries: in fact, Italy, Portugal, Bulgaria, Liechtenstein, Slovenia, Luxemburg, and Moldova were not part of these programs (Halley et al. 2021). Still, beavers naturally recolonized some of them (Fasel 2018;Herr et al. 2018), except for Portugal (Halley et al. 2021). About Italy, although Halley et al. (2021) report a dubitative recolonization, some very recent findings of beavers' signs prove the occurrence of this species in two different Italian locations: in the north-east of the country (near the town of Tarvisio, bordering Slovenia) (Pontarini et al. 2018;Loy et al. 2019), and in Central Italy (Tuscany and Umbria regions) Viviano et al. 2022), even with a possible reproduction event Pucci et al. 2021). Thus, we are now maybe witnessing a natural recolonization in north-eastern Italy, which could continue spontaneously, or which could be assisted by direct or indirect human mediation (i.e., reintroduction and/or area protection). The presence of the species in central Italy can instead be traced to an illegal reintroduction Pucci et al. 2021). On the contrary, Portugal has no similar chances to date but could be involved (as done for other European countries) in reintroduction programs, especially after the very recent findings of beavers' signs at the Spain-Portugal border (Calderón et al. 2022).
Moreover, one must simultaneously consider the biodiversity crisis that species are facing (and will face in the next future) considering the ongoing global change, which threatens the Eurasian beaver itself as well as the ecosystems it contributes to shape. One of the main solutions to handle these issues is to support Eurasian beaver conservation through a suitability distribution model, which could both support reintroductions and predict current and future suitable territories for this species. In fact, the IUCN "Guidelines for Reintroductions and Other Conservation Translocations" stress the importance of verifying the availability of suitable habitats as a key part of planning reintroductions (IUCN/SSC 2013). Also, the knowledge of areas' suitability in time was shown to strongly aid protection measures, for instance finding gaps or fallacies in protected areas' (PAs) networks (Iannella et al. 2018;Cerasoli et al. 2019).
In this research, we aim at supporting both possible reintroduction measures and protection/management of areas, through a zonal-focused modeling approach. Starting from a climate-based Ecological Niche Modelling, we build finer Species Distribution Models sensu Soberon and Peterson (2005) through the weighted overlay framework presented by Iannella et al. (2021), assessing the suitability for each European river. We integrate into the weighted process the environmental predictors which are reported to be ecologically important for the Eurasian beaver. In fact, we encompass not only climate change scenarios but also rivers' flow changes inferred for the future (Blöschl et al. 2019), as these are one of the main predictors of the beavers' occurrence (Alakoski et al. 2020). Also, we consider habitats used and rivers' human disturbance level, to support landscape-scale management, projecting suitability from a river scale to basins, thus accounting for a higher territorial level. Finally, we couple the weighted models' predictions and the post-modelling analyses performed with a country-scale protected areas network analysis. Indeed, we assess PAs' coverage focusing both on Italy and Portugal, to identify and support the management of reintroductions or natural recolonizations.

Target species and study area
The target species is the Eurasian beaver, a semi-aquatic rodent currently distributed in Central, Eastern and Northern Europe, with a recent increase of its population (Halley et al. 2021). Though, from the eight refugial populations where the species managed to survive from its initial decline to between the 19th and 20th centuries (Nolet and Rosell 1998), strong bottleneck deriving from genetic drift is observed (Halley et al. 2021). Notwithstanding natural recolonization and reintroduction programs, Italy, Portugal and southern Balkans remain uncolonized, with an exception for Italy. In fact, signs of the presence of the beaver have recently been found in two distinct areas in Italy as mentioned above: in the north-east of the country (near the town of Tarvisio, bordering Slovenia, and in Val Pusteria, Trentino Alto-Adige region), the species' occurrence is due to natural recolonization (Pontarini et al. 2018), while the central Italian records (Tuscany and Umbria regions) are attributable to illegal reintroductions Mori et al. 2022;Viviano et al. 2022).
We gathered occurrence localities for the Eurasian beaver from the online dataset GBIF (GBIF 2022); we applied a temporal (year > 1970) and spatial precision (coordinate uncertainty < 100 m) filter to discard old and low spatial-quality data. Then, we processed the resulting occurrences through the 'Spatially Rarefy Occurrence Data for SDMs' tool (set at = 10 km) of the SDMtoolbox 0.9.1 (Brown et al. 2017) for ArcGIS Pro 2.9.3 (ESRI Inc. 2022) to make the spatial resolution of both predictors and occurrences comparable (Sillero and Barbosa 2020).
We selected two different study areas' scales for the purposes of the present work. The first one is the European continent, considering the beaver's increase of populations occurred there and the threats that global change poses on these. To support the management of natural recolonization and/or reintroduction projects, we then focus on two sub-areas, targeting the biogeographic regions referring to Italy and Portugal, and not including Balkans considering that a detailed paper is already available (Smeraldo et al. 2017).

Environmental predictors
For the purposes of the species distribution modelling, we took into consideration several environmental variables. To encompass the climatic aspect, we downloaded nineteen bioclimatic variables from WorldClim (version 2.1, 2.5 arc-min) (Fick and Hijmans 2017) for both present and future (2050 and 2070) conditions. For each future scenario, we selected the Shared Socioeconomic Pathways (SSPs) 2.45, 3.70, and 5.85 to involve all but one (the 1.26, the most optimistic) different possible trajectories (Riahi et al. 2017).
We collected topography data from ETOPO1 (~ 1 km resolution, available at https:// www.ngdc.noaa.gov/mgg/global/) and processed them through the 'Surface Parameters' tool in ArcGIS Pro. Specifically, considering that the speed of running waters is an important parameter for the Eurasian beaver to settle (Hartman and Törnlöv 2006;Dittbrenner et al. 2018), the Profile curvature was preferred to the Slope alone. This predictor is applied when one wants to differentiate the acceleration or deceleration of flow at a specific surface (ESRI Inc. 2022). This analysis is then useful to identify areas where slow-(or fast-) flowing water-adapted species may occur. We thus calculated the Profile curvature raster for the whole study area.
We obtained rivers' spatial data and their corresponding connectivity and annual mean flow from Grill et al. (2019). Then, for the "human impact" information, we performed a raster calculation ('Weighted sum' tool in ArcGIS Pro) by averaging three different parameters: the degree of regulation of the rivers, the road density, and the degree of urbanization of the territory, obtaining this information from the HydroATLAS dataset (Lehner and Grill 2013;Linke et al. 2019). Also, we used sub-basin data from the HydroBASIN dataset (Linke et al. 2019), and starting from the data reported in (Blöschl et al. 2019) regarding the projected changes in flow rate in European rivers in the future, we created a map for future changes of rivers-related spatial information.
About the habitat predictors, we downloaded the 100 m-precision map of global terrestrial habitats from the Copernicus repository (https://land.copernicus.eu/global/products/lc), containing 47 different terrestrial habitat types, as defined by the IUCN habitat classification scheme (Jung et al. 2020).

Ecological niche modelling
The nineteen Worldclim bioclimatic variables were checked for multicollinearity through the Variance Inflation Factor (using the 'vifstep' algorithm of the 'usdm' R package (Naimi 2015), threshold = 10). To reduce the variability possibly caused by the use of individual General Circulation Models (GCMs) in future projections (Stralberg et al. 2015), we selected and managed three different GCMs, namely the BCC-CSM2-MR (Wu et al. 2019), the IPSL-CM6A-LR (Boucher et al. 2020) and the MIROC6 (Tatebe et al. 2019).
We used the 'biomod2' package (Thuiller et al. 2016) to build the ENMs. This package allows to create ensemble models, which are the result of the combination of individual algorithms' outcomes, corresponding to different modeling techniques. We first generated five sets of 2,000 pseudo-absences each using the Surface Range Envelope algorithm (quantile = 0.05) ( Barbet-Massin et al. 2012). We parametrized the single models through the 'BIO-MOD_ModelingOptions' as follows: GLM: type = "quadratic", interaction.level = 3; GBM: n.trees = 10,000, interaction.depth = 3, cv.folds = 10; GAM: type="s_smoother", interaction. level = 3; RF: n.trees = 10,000, bag.fraction = 0.75, shrinkage = 0.001. Then, the 'BIOMOD_ Modeling' function was used to calibrate the single models, and the 'BIOMOD_Ensemble-Modeling' algorithm was subsequently run to obtain the Ensemble models, using as input only single models with high performance (AUC > 0.8 and TSS > 0.7) (Cerasoli et al. De Simone et al. 2020;Di Gregorio et al. 2021). We chose the 'wmean' algorithm to merge single models, as it averages models into the Ensemble one based on a weighted mean function, in which weights are the discrimination scores (Thuiller et al. 2016). Finally, we obtained the projections to future scenarios by using the 'BIOMOD_EnsembleForecasting' function.
To buffer the uncertainties deriving from future projections calibrated on each GCM and concurrently merge them into one year/SSP single model, we first checked for model extrapolation (i.e., the dissimilarity from the calibration conditions) by assessing the Multivariate Environmental Surface Similarity (MESS) (Elith et al. 2010), computed through the function 'mess' of the 'dismo' package (Hijmans and Elith 2016). Then, the MESS maps were included in the Multivariate Environmental Dissimilarity Index (MEDI) algorithm (Iannella et al. 2017), a form of weighted raster average which down-weights extrapolation based on projections and MESS, reducing the uncertainties in future projections.

Data preparation and weighted overlay process
To merge all the information into the weighted maps and the corresponding weighted workflow (Fig. 1a), we transformed in a 1-to-10 scale all the previously mentioned predictors to make them comparable each other through the 'Reclassify' tool in ArcGIS Pro.
The connectivity value available from the Grill et al. (2019) dataset, which divides rivers based on free flowing status into three categories, was reclassified to new values (i.e., 1, 5, and 10) for each study area's river. Moreover, we took into account the ecological needs and hydrographic preferences of the Eurasian beaver (Hartman and Törnlöv 2006;Dittbrenner et al. 2018) and reclassified the rivers' data accordingly, on the basis of their reported average annual flow (Grill et al. 2019). Also, the "human impact" information deriving from the processing of Road Density, Urban Extent and Degree of Regulation parameters of the RIVERAtlas dataset, was converted to the new scale, where the lowest human impact has the highest value.
To reclassify the profile curvature and habitat needs to the common scale, we calculated the frequencies of occurrence in each value of these two predictors by extracting this information in GIS ('Extract multi values to points' tool), following the procedure of .
The climate-related Ensemble Models' output maps were reclassified as well, converted from their 0-to-1,000 output to the 1-to-10 scale.
Finally, we used the 'Weighted Overlay' tool in ArcGIS Pro to merge the predictors (free flowing status, annual average flow, human impact, profile curvature, terrestrial habitat type, and bioclimatic-based ensemble models) into a single map. This tool merges a given set of rasters, sharing a common evaluation scale, through a weighted averaging in which each input raster is also assigned a specific percentage set by the operator. This workflow has been shown to increase the models' discrimination power if compared to the Ecological Niche Models alone , deepening the species' distribution model to a more realistic one.
About the future scenarios, we kept all variables as fixed except two, namely the climatebased Ensemble models and the annual average flow. For this latter, we performed a raster algebra to increase or decrease the values based on the indications of published predicted future flow change rates in the study area (Blöschl et al. 2019).

Landscape modelling and conservation applications
We applied the Zonal Statistics tool (which allows to perform stats on a geometry based on the values or attributes of another one) to the weighted models obtained and to the HydroBASINS dataset, which reports the boundaries of the European sub-basins (at Sub-Basin level = 10) (Lehner and Grill 2013). We thus linked the values of the rivers' averaged weighted suitability to the hydrographic basins, leading to an explicit spatial information at a landscape scale.
As one of the purposes of the present work is to support possible reintroductions of the Eurasian beaver, we focused on Italy and Portugal, countries where the species was present in historical times but where stable populations have still not re-established after its eradication. We thus took into consideration the PAs occurring in these countries and measured the environmental suitability within them. Subsequently, to summarize the vast amount of information, we divided Italy and Portugal based on their biogeographic sectors: Italy has therefore been divided into the Alps, the Po plain, the Northern, Central and Southern Apennines, and the Apulian province, excluding the two main islands. Similarly, we separated Portugal into four areas, namely the Luso-Extremadurense, the Carpetano-Iberico-Leonesa, the Cantabro-Atlantica and the Gaditano-Onubo-Algarviense sectors, based on the information from Costa et al. (1998).
We finally assessed the environmental suitability within the PAs for each biogeographic sector, thus obtaining these values at all the landscape scales: from the single river to the sub-basin level, for the whole Western Palearctic area. Aside from the numeric results, the maps we obtained permit to deepen the suitability of even a single river, which may be helpful if one wants to properly assess conservation and/or reintroduction measures within a certain territory (e.g., a protected area).

Results
After the filtering process, we gathered 2,848 occurrence localities for the whole European continent (Fig. 1b); also, taking into account the VIF results, we used for the modelling process nine bioclimatic variables, namely BIO2 (mean diurnal range), BIO6 (minimum temperature of the coldest month), BIO8 (mean temperature of the wettest quarter), BIO9 (mean temperature of the driest quarter), BIO11 (mean temperature of the coldest quarter), BIO15 (precipitation seasonality), BIO17 (precipitation of the driest quarter), BIO18 (precipitation of the warmest quarter) and BIO19 (precipitation of the coldest quarter). About the Eurasian beaver's habitat preferences (obtained through the extraction in GIS) used in the weighted process, we found the highest frequencies in the "Arable land", "Boreal forests", "Temperate forests", and "Boreal shrubland" categories (Appendix S1a). Also, when considering the Profile raster, we found the highest frequency in flat (or nearly-flat) areas, with a decreasing frequency as the slope increases (Appendix S1b).
The weighted map we obtained for the current conditions (Fig. 2) depict a suitability that is mostly in line with the species' current occurrence (Fig. 1b), with a plethora of small-and medium-sized rivers being more suitable to the larger ones ( Fig. 2 and its example zoomed insets).
Indeed, some territories not currently hosting the species are predicted as suitable, such as for a large part of Belarus, Estonia, Latvia, Lithuania, and western Russia (even though this latter shows slightly lower values), Finland, Great Britain, Croatia, and Italy (Fig. 2). On the contrary, the Iberian Peninsula, the Alps, some major eastern Balkans territories, and Turkey report a general low-suitable rivers' trend. This is also observed at a "surface" spa-tial scale: when projecting suitability to sub-basins, they spread the rivers' suitability values to a territorial scale, mostly highlighting the areas where low suitability occurs (Fig. 3a). Differences between the very high and very low patches become more marked, highlighting a wide eastern European area (from eastern Poland to southern Finland), but limited to the westernmost part of Russia (Fig. 3a). Similarly, a continuous suitable area is found in the southern Scandinavian peninsula (mostly occurring in Sweden), and not extending further north of the 10th parallel (Fig. 3a). Some medium-extent suitable areas appear in central Europe; even if interspersed with continuous major rivers' basins, they show high suitability values (Fig. 3a). About southern Europe, scattered and small-sized highly-suitable areas occur in the Balkans and, to a lesser extent, in northern and central Italy; the Iberian Peninsula reports some more scattered and smaller patches, with lower suitability values (Fig. 3a).
In the future, a general deterioration of suitability in many wide areas is predicted, with 2070 being worse than 2050 (Fig. 3b,c). The central and northern Europe, which currently host highly suitable areas, are predicted to lower towards medium values (Fig. 3b,c). Specifically, the areas falling into Poland, Belarus, and Ukraine are predicted to shifting towards medium or low suitability, while Finland e Scandinavian peninsula remain suitable, with high suitability bypassing the 10th parallel. Also, some northern Russian territories will become more suitable with respect to the current conditions, even though this suitability increase would be negligible (Fig. 3b,c).
When applying this information to the Portuguese and Italian PAs' network, a complex framework emerges (Fig. 4). In Portugal, most of the PAs are predicted to currently cover medium-low suitable territories, with some very low ones in the inland (Carpetano-Iberico-Leonesa and Luso-Extramadurense, respectively PO2 and PO4) except for the southern part of the Luso-Extramadurense (the Beja region, Fig. 4a). On the contrary, some medium and medium-high patches occur in the coastal biogeographic sectors (PO3, the Gaditano- Onubo-Algarviense), as found, for instance, in the Leiria area (Fig. 4a). About Italy, we obtained a very diverse asset, even within a single biogeographic sector. For instance, the PAs hosting a high suitability occur at the westernmost and easternmost extremes of the Alps (IT1), while its central part results instead in very low values (Fig. 4a). The Padano-Venetian plain (IT2) is rather poor of PAs, but some of them report high or medium suitability; the Northern and Central Apennines (IT3 and IT4) score high suitability, especially in PAs occurring in the Apennine chain and in the southern coastal part of IT3 (the Maremma and Metalliferous Hills territories) (Fig. 4a). About the Apulian and Southern Apennines sectors (IT5 and IT6), we find a general medium-low and low patches, except for the karstic semi-mountainous territory of the Murge plateau, at the center of the IT5 (Fig. 4a).
When projecting these results to future, some changes throughout the 2050 and 2070 scenarios may be noted in Italy, while minor ones resulted for Portugal (Fig. 4b,c). In fact, some suitable patches in the westernmost and easternmost extremes of the Alps, in the Maremma and Metalliferous Hills (IT3), in the Central Apennines, and in the Murge Plateau (IT5) are predicted to progressively increase from 2050 to 2070 (Fig. 4b,c). We also found a similar future, but with a lower magnitude, for the medium-suitability patches occurring in Leiria and Beja areas (PO3 and PO4, respectively). To give an overview at a landscapescale of these specific trends, we report the averaged future changes (Fig. 5). As mentioned before, we predicted minor suitability changes within Portuguese PAs, with an overall trend of slight increases or decreases in suitability for all the four biogeographic sectors (Fig. 5a). Differently, we observed more pronounced variations for Italian PAs, with considerable losses for Padano-Venetian plain, Northern and Central Apennines sectors and slight gains for the other ones (Fig. 5b). Although the predicted negative changes or stability, it must be noted that the starting values are not less than a medium-high suitability.

Discussion
The recolonization of territories by the Eurasian beaver is a today a fact; however, one must remember that this event basis its success on reintroduction programs (Halley et al. 2012). Reflecting the IUCN guidelines, where a strong indication of ceased pressures and threats on the target species is recommended to proceed with the reintroduction (IUCN/SSC 2013), these in fact succeeded. Starting from the legal protection of the species (which, for instance, cannot be hunted as in the past), and following genetic and territorial assessment, country-scale actions were accomplished (Halley et al. 2021). Indeed, the local dimension of the activities carried out since now could represent a limitation for future conservation strategies aiming at halting possible range variations triggered by the ongoing global change.
Notwithstanding the great importance of the Eurasian beaver in terms of ecosystem services-related conservation issues Thompson et al. 2021), an overview of suitability projections for future scenarios was currently missing for the Europe. Moreover, considering that Eurasian beaver populations develop on a watershed-scale (Hartman 1994), we managed to combine the general view with this spatial data.
In our research, the model calibrated on current conditions mirrors the Eurasian beaver known distribution, indicating not only its climatic suitability, but also representing its other ecological requirements (topography and habitat), thus narrowing the climate-based predictions to more specific territories. Also, the current model highlights some patches where the species is not reported, where it is likely to spontaneously recolonize, or where reintroduction plans are proposed (Treves et al. 2020;Halley et al. 2021). For this last case, the Italian peninsula and Wales are the two areas where our results could support specific planning for future management, especially for the future scenarios, in which the currently-suitable areas are predicted to shrink by our models. If considering the Balkans' countries, a growth of suitable territories is instead predicted (even though with a slight high-suitable patches decrease), thus Bosnia and Herzegovina and Serbian stakeholders should pay special attention in the ongoing 'establishment' phase of the former country and expected population increase for the latter (Halley et al. 2021). About the main European River, the Danube, we found its basin showing a low environmental suitability, since the river itself is unsuitable for the creation of dams and burrows for beavers, considering its size and flow, and manmade modifications. However, high values of environmental suitability are found in the tributaries of the Danube almost along its entire length. Thus, we can consider the Danube itself as unsuitable for creating core areas for the Eurasian beaver, but still important for its movement and connectivity, confirming the spatial and populational observations recently published (Halley et al. 2021). In contrast, Spain, where beavers' reintroduction in the past led to a medium-sized population currently expanding (Halley et al. 2021), should take preventive conservation actions, as a suitability decrease is predicted by our models for the areas currently inhabited by the Eurasian beaver population.
An important remark must be done for the territories which acted as refugia (following Halley et al. 2021) during the years when the Eurasian beaver suffered declines, where we found a decrease of the watershed-based suitability. The Russian and Belarusian patches resulted in the highest suitability decrease. The territories where the Norwegian refugia occurs are predicted to keep high suitability in future scenarios, while the German and French ones are predicted to reduce both in size and suitability. This is an alarming result if one considers the importance of these areas in terms of their historical role, as well as for the possible ecological functionality persisted (Thompson et al. 2021) which may be rapidly lost, as well as for the genetic bottlenecks there occurring (Halley et al. 2021).
A separate mention concerns Finland, where the Eurasian beaver coexists with the North American beaver (Castor canadensis) which was erroneously reintroduced in the past when the taxonomic status was still uncertain (Parker et al. 2012). Finland has high values of environmental suitability, but in this case the expansion of the Eurasian beaver is hindered by the presence of the North American beaver. This latter often seems to come out as the winner of ecological competition (Parker et al. 2012), even though there are opposite cases, such as for Austrian and Polish ones (Halley and Rosell 2002;Danilov 2009). The eradication of the North American beaver, considered an invasive species in Finland, would support the expansion process of the Eurasian beaver, in a territory bearing large highly-suitable patches.
Recent research state that, generally, Eurasian beaver population are expected to slow its increase in the next future and become stable in the next decades (Wróbel 2020;Halley et al. 2021). The rapid populations' growth was observed to be a consequence of natural recolonization or reintroductions (Halley et al. 2021), a phenomenon which deserves attention for countries which want to successfully reintroduce this species. As studied in our case, Italy and Portugal (which are the only two southwestern European countries where the species occurred in the past) should address robust plans if they want to achieve a proper recolonization process. Protected areas play a key role in the management of the territories they cover, and usually support conservation at the forefront; indeed, in the case of Italian and Portuguese ones, the general status is rather complex.
For instance, Italy currently hosts the overall highest suitability within the Padano-Venetian plain protected areas, where these are scattered and suffer from a severe fragmentation caused by the greatest urbanization of the whole Italian peninsula. The only beaver occurrence known site falls within the Northern Apennine biogeographic sector, for which a future decrease in suitability is predicted, and where protected areas are dispersed in wide landscapes. A slightly lower suitability results for Central Apennines, which instead are not currently recolonized, but where highly-suitable and interconnected patches occur. Going southwards, the other sectors (Southern Apennines and Apulia) seem not to be colonizable, as the basins' suitability drastically lowers, with corresponding lowering of the protected areas occurring there. A more complex scenario is predicted for Alps' sector, where a general low suitability resulted for its protected sites. If analyzing in detail, one can note an eastern and western "polarization": in fact, the central Alps basins result not suitable, while the ones at the farthest are not only suitable, but near the other European populations (French, and Austrian-Slovenian ones). This may support a natural recolonization, or reintroduction programs, such as the one studied for the Piedmont region (Treves et al. 2020).
On the other hand, Portugal still does not have any Eurasian beaver occurrences within its borders, thus possible reintroduction programs could take advantage of the present research. Environmental suitability is generally low, except for some protected areas where average suitability values are found (in the Gaditano-Onubo-Alsarviense sector). However, these areas appear fragmented and poorly connected to each other. Conversely, greater proximity between protected areas is present in the Luso-Extremadurense sector, but only a few protected areas have good suitability values. The Carpetano-Iberico-Leonesa sector has low suitability values but is a relevant region because it borders the Spanish populations of Eurasian beaver. This leads to believe that, as things stand at present, a natural recolonization of Portugal is unlikely in the next years, and a highly-informed reintroduction is the only possible option.
The Eurasian beaver is a key species primarily due to its role as an ecosystem engineer, and can produce diverse ecosystem services, capable of services worth millions to hundreds of millions of US dollars (USD) per year. Concrete examples concern the provision of habi-tat and biodiversity (estimated at USD 133 million), and the sequestration of greenhouse gases (estimated at USD 75 million) (Thompson et al. 2021). Other ecosystem services provided by the Eurasian beaver include water purification, moderation of extreme events, nutrient cycling, recreational hunting and fishing, water supply, and non-consumptive recreation (Brazier et al. 2021;Thompson et al. 2021).
The ability of the beaver to restore habitat can help improve the ecological status of rivers, especially through some actions that the beaver carries out in rivers: maintenance of levees through a dynamic of their creation / destruction (Campbell-Palmer et al. 2016); spreading plants through foraging (Campbell-Palmer et al. 2016); replication of small dams in addition to the construction of dams (Pollock et al. 2014); increase habitat heterogeneity, forming different microhabitats, such as ponds and wetlands (Law et al. 2017). For Portugal in particular, the restoration actions of the beaver in rivers could reduce the costs of post-fire human interventions, carried out between 2018 and 2020 to reduce flood events and their risks, estimated at around 11.43 million € (https://sniambgeoviewer.apambiente. pt/GeoDocs/geoportaldocs/Docs/APA_Relatorio_Incendios2017_RCM11A2018.pdf). Some of the central areas of Portugal, where some protected areas with medium suitability have emerged, could benefit from the presence of the beaver. Also, further fire-related impacts may be partially mitigated by the Eurasian beaver, as dams formed by beaver (in this case, Castor canadensis) improve fire resistance of riparian vegetation (Fairfax and Whittle 2020). In Portugal, in 2017, fires in the central regions of the country caused the loss of 40-60% of the forested areas of some districts (DGAPPF, 2017); in fact, it should be remembered that Portugal is among the countries of the European Union most affected by fires, which burned an annual average of 147,000 hectares of land between 2000 and 2018 (San-Miguel-Ayanz et al. 2019). Reintroduction projects in these areas could contribute to improve the ecological status of rivers, partially solving the problem of fires, and adding further ecosystem services, not only for Portugal (as reported in this specific case) but also for any historically-colonized territory.
The beaver is a species that can also increase wildlife tourism, as demonstrated by the case study in Otter, England, where the beaver has been reintroduced (Auster et al. 2020). The reintroduction in this river has resulted in an advantage for local businesses, with the magnitude of the benefit depending on the commercial initiatives that have exploited the opportunity (Auster et al. 2020). Other forms of beaver-centered tourism are "beaver safaris" (Halley et al. 2012).
In view of a possible reintroduction of the Eurasian beaver in Italy or Portugal, it is necessary to consider also the possible presence of other species that may interact with it. The beaver is in fact a prey for many mammals widely present in the western Palearctic, such as the grey wolf (Canis lupus) or the fox (Vulpes vulpes) (Kile et al. 1996). In addition to predation, the presence of the beaver can also be the cause of competition events with a species partially occupying a similar niche, such as the coypu (Myocastor coypus), an invasive species in Europe. Studies in central Italy reported an intermediate overlap in the activity pattern between the newly-found beavers ) and coypu, with the latter having more diurnal habits . However, slight ecological differences occur between the two species, in terms of feeding behavior ) and habitat needs (Ruys et al. 2011).
Another important factor halting beaver expansion is the presence of man-made dams, which hinder river connectivity. In fact, given that a natural water regime is necessary for the survival of beavers (Nolet & Rosell, 1998), man-made dams can alter it, preventing individuals from dispersing. Also, beavers' food availability, supplied by intact riparian forests (Nolet & Rosell, 1998), diminishes when man-made dams occur. Thus, where populations are present in areas affected by artificial dams, it is important to evaluate technical measures to guarantee the connectivity of the river system, such as adopting naturalistic engineering alternatives.
Finally, in a reintroduction context, one should consider the topics of diseases and parasites. In fact, the relationship between coypu and beaver is concerning also considering the possible transmission of pathogens. Indeed, coypu is an intermediate host for Echinococcus multilocularis, one of the most pathogenic parasitic zoonoses in central Europe. Several cases of Echinococcus multilocularis in beaver have been reported in Switzerland, Austria, and Great Britain (Janovsky et al. 2002;Cronstedt-Fell et al. 2010;Barlow et al. 2011). Therefore, the choice of possible reintroduction sites also depends on the co-occurrence of coypu, which thus represents, among other species, both as a direct competitor and as a possible source of pathogens. About Italy, where the coypu already occurs, appropriate measures must be taken; in Portugal and Spain, the suitable habitats of coypu are predicted to expand in the future (Schertler et al. 2020), thus specific actions should be addressed to achieve successful reintroductions (for the former country) or to preserve existing populations (for the latter).

Conclusion
Our study modelled for the first time the European beaver suitable habitats at a continental scale, allowing to draw important considerations on the future of a species in the past close to extinction, but which today is once again spreading. Although the Eurasian beaver is a generalist and adaptable species, a difficult future is predicted in the southernmost areas of the western Palearctic, while in central Europe and in Scandinavia, populations may persist. In Italy (even in the future) large areas suitable for reintroduction will remain, especially in the Central Apennines, while for Portugal only a few protected areas would be suitable. Further studies are needed to understand the feasibility of reintroductions of the Eurasian beaver in Portugal, or the persistence/colonization/reintroduction of this species in Italy. Despite the small population occurring in Central Italy being the result of illegal reintroduction, we recommend keeping and protecting it, following the Spanish case, where a similar situation occurred recently and the beavers have been protected, adopting the provisions of the Habitat Directive. Moreover, given the evidence of reproduction, it is possible that in the near future this population could spread more in Central Italy. Therefore, we recommend maintaining populations after carefully monitoring their spatial range and concurrently examining and sensitizing public opinion. It is also important to evaluate the potential ecological impact of the species in the basins where it occurs, taking into account the fauna currently present and its possible interactions.
Given the continuous populations' growth and the expansion of the range, it is likely that the beaver will be able to naturally recolonize all its historical areal. Therefore, reintroductions would favor this process, which in the case of Italy is already taking place. In the case of reintroductions, monitoring programs will be needed, as the beaver tends to move widely when colonizing a watercourse to find its ideal habitat. The benefits linked to the beaver activity are different and would allow direct gains (wildlife tourism) and savings from the beaver activity in waterways (construction of dams, creation of wetlands, resistance against fires).