Modelling the amphibian chytrid fungus spread by connectivity analysis: towards a national monitoring network in Italy

The emerging amphibian disease, Batrachochytrium dendrobatidis (Bd), is driving population declines worldwide and even species extinctions in Australia, South and Central America. In order to mitigate effects of Bd on amphibian populations, high-exposed areas should be identified at the local scale and effective conservation measures should be planned at the national level. This assessment is actually lacking in the Mediterranean basin, and in particular in Italy, one of the most relevant amphibian diversity hotspots in the entire region. In this study, we reviewed the available information on Bd in Italy, and conducted a 5-year molecular screening on 1274 individual skin swabs belonging to 18 species. Overall, we found presence of Bd in 13 species and in a total of 56 known occurrence locations for peninsular Italy and Sardinia. We used these occurrence locations and climate data to model habitat suitability of Bd for current and future climatic scenarios. We then employed electric circuit theory to model landscape permeability to the diffusion of Bd, using a resistance map. With this procedure, we were able to model, for the first time, the diffusion pathways of Bd at the landscape scale, characterising the main future pathways towards areas with a high probability of Bd occurrence. Thus, we identified six national protected areas that will become pivotal for a nationally-based strategic plan in order to monitor, mitigate and possibly contrast Bd diffusion in Italy.


Introduction
As part of the recent "biodiversity crisis," many amphibian populations are declining worldwide (e.g., Blaustein et al. 1994;Wake and Vredenburg 2008;Catenazzi 2015;Scheele et al. 2019). Indeed, the comprehensive analysis by Stuart et al. (2004) indicated that about one-third of all amphibian species is threatened with extinction, while almost half are experiencing regional or local declines due to often interacting causes, such as habitat reduction, pollution, climate change and emerging diseases. Among pathogens, the recently described chytrid fungus Batrachochytrium dendrobatidis  is considered the main cause of population declines in different continents, driving species to extinction in Australia, South and Central America Scheele et al. 2019). Batrachochytrium dendrobatidis (Bd) is a highly virulent pathogen that infects the skin of all the three orders of amphibians [i.e. Anura (frogs), Caudata (salamanders) and Apoda (caecilians)], causing chytridiomycosis, a frequently lethal disease that produces immunosuppression, depletion of plasma electrolytes and cardiac electric dysfunctions in amphibians (Berger et al. 2004(Berger et al. , 2016. The global Bd occurrence and its effects on amphibian populations have been recently reviewed and mapped (Olson et al. 2013;Lötters et al. 2009;Scheele et al. 2019). To date, six different Bd lineages have been identified by multilocus sequence typing, but only the global pandemic lineage (GPL) seems associated to widespread chytridiomycosis outbreaks that caused populations declines . The main factor spreading the pandemic Bd lineage in different parts of the world and in different time periods is the international trade of amphibians and other aquatic animals for food, research, collection or company (Olson et al. 2013).
Various studies, in the last fifteen years, mapped the known distribution of Bd-infected amphibian populations and analysed the possible future consequences of this disease at the global or continental level (e.g. Ron 2005;Becker and Zamudio 2011;Doherty-Bone et al. 2020;Ribeiro et al. 2020). However, this kind of assessment is lacking at the regional level in the Mediterranean biogeographic area, a well-known hotspot of natural and human-adapted ecosystems (Blondel and Aronson 1999;Myers et al. 2000). In particular, the central Mediterranean region has been shown as a potential suitable area for Bd by a variety of global studies (e.g. Fig. 4 in Ron 2005;Fig. 3 in Lötters et al. 2009;Fig. 2 in Liu et al. 2013). At the centre of the Mediterranean basin, Italy represents one of the most relevant hotspots of biodiversity with a high concentration of amphibian species and in particular of endemics (Sindaco et al. 2006). In fact, due to its central geographic position within the Mediterranean, complex geological history, contrasted geomorphology, variable climates and the long-lasting coevolution between rural landscapes and wildlife (Cevasco et al. 2015), Italy hosts a highly diverse and unique amphibian fauna, comprising about half of all amphibian species described in Europe (Table 1; Temple and Cox 2009;Rondinini et al. 2013). Various local studies have been published on Bd occurrence in Italy (see Table 1), but up to now only one mass-mortality event has been observed in Sardinia Tessa et al. 2013). In this island the pandemic lineage GPL has been recorded . Besides, in some areas along the Apennine mountain range the observed decline of the Apennine yellow-bellied toad, Bombina pachypus, was attributed to Bd infection (Stagni et al. 2004;Canestrelli et al. 2013). However, declines are also observed in other areas where Bd has been screened for, but is apparently absent (Canessa et al. 2013a(Canessa et al. , 2019. In these latter areas, the Apennine yellow-bellied toads' declines were attributed to major habitat changes rather than pathogens (Canessa et al. 2013b). Therefore, the ecological effects of chytridiomycosis on the conservation status of Italian amphibians remain poorly understood, and in some cases enigmatic. For this reason, is of primary importance to delineate a nationwide monitoring and intervention plan, and model Bd occurrence to forecast possible future outbreaks. In this research, we reviewed all the available information on Bd occurrence in Italian amphibian populations and we added original data, obtained over a five-year screening on peninsular Italy (Grasselli et al. unpublished data). Then we used both bibliographic and original data to model the present and near-future pathways of Bd diffusion in Italy, by both building habitat suitability models and ecological connectivity models at the landscape scale (McRae et al. 2008;Dickson et al 2019). This approach allowed us to forecast the most probable pathways of Bd diffusion and to better understand the impact of this disease on amphibian populations, within and outside protected areas, and to possibly identify the best mitigation measures to be implemented at the national and regional scale.

Study rationale and framework
The principal aim of this study was to employ occurrence data of Bd in Italy, both from published research and original data, to model its current and future habitat suitability, to identify possible pathways in order to define a national monitoring network and to drive conservation measures. For that purpose, we collected all the available information on Bd occurrence data from published research that confirmed the presence of Bd by PCR and that possessed a reliable geographic locality. Moreover, we collected new data from 1274 skin swabs obtained between 2015 and 2019. Occurrence data were used to build habitat suitability models for both current and near future climatic conditions (year 2050), considering that climate has been found to have a predominant effect on Bd occurrence (e.g. Liu et al. 2013;Xie et al. 2016;Flechas et al. 2017). These habitat suitability models, in turn, have been used as a basis for building connectivity models by means of electric circuit theory, in order to delineate landscape permeability and to identify actual and future channels of Bd diffusion in Italy. Our aim is to identify which national protected areas will be most likely interested by massive diffusion of Bd in the near future. In this way, we will be able to design a rapid response monitoring network and to close possible knowledge gaps on the diffusion of Bd, concentrating the monitoring effort in those areas that will be interested by a rapid Bd spread (Fig. 1).

Bd screening in Italy: field sampling and occurrence data
In the present study, we used data from multiple sources: published records and original data of Bd in peninsular and insular Italy. Bibliographic data were obtained from published research spanning from 2004 to 2015 (Table 1) and yielded 33 occurrence locations ( Fig. 1; Supplementary Material Appendix 1, Table 1). Original data belong to a nationwide screening, conducted between 2015 and 2019 which mainly occurred in four national protected areas (Grasselli et al. unpublished data). We analysed 1274 amphibian swabs, belonging to 18 species, from 114 sampling locations. Swabbing procedures were standardised ) and field operators were formed during a 2-days field course held by the University of Genova, dedicated to the formation of national Parks' staff and 1 3 provided all participants with written instructions. Each swab was obtained by firmly rubbing its tip on the abdomen, pelvic area and leg pits of the captured individual (Boyle et al. 2004). All locations of Bd presence, both from published research and national screening, were georeferenced on a 1 × 1 km grid.
The extraction of nucleic acids from skin swabs for qPCR was performed as described by Boyle et al. (2004). Briefly, nucleic acids were extracted from swabs using 200 μl Prep-Man Ultra along with 30-40 mg of Zirconium/silica beads (Biospec Products) and then incubated in a bead-beater. After centrifugation (1 min at 13 × 10 3 × g) and incubation in a boiling water bath for 10 min, supernatant was recovered and stored at-80 °C. The Bd SYBR green assay (Canessa et al. 2017)

Habitat suitability models
Habitat suitability models were fitted using Maxent (version 3.4.0), a machine-learning algorithm that creates a model of habitat suitability (or probability of presence) for a given species, based on gridded predictors at observed species' occurrence locations (Phillips  Guisan et al. 2017). The obtained model of habitat suitability can then be projected onto a different set of predictors, representing a different landscape or a future condition, thus obtaining a suitability model for another region or for a future scenario for the species of interest (Elith et al. 2011;Merow et al. 2013). This approach has been proved to outperform other methods based on presence-only (or presence background) data (Elith et al. 2006). Moreover, it also may have some advantages for the specific application on Bd data: (i) by relying on presence-only data there is no need for absence data that could be difficult to collect in the case of Bd, given the non-exhaustive sampling, the mixed origin of the data (i.e. bibliographic and original), (ii) Maxent algorithm proved to be reliable with a small amount of presence locations (Pearson et al. 2007;van Proosdij et al. 2016;Bacigalupe et al. 2019), which is the case of Bd occurrence data in the present study.
As environmental predictors, 19 bioclimatic variables representative of current and future (2050) Fig. S1). As a result, we retained 9 bioclimatic variables: BIO1-Annual mean temperature, BIO2-Mean temperature diurnal range, BIO5-Max temperature of warmest month, BIO6-Min temperature of coldest month, BIO8-Mean temperature of wettest quarter, BIO9-Mean temperature of driest quarter, BIO12-Annual precipitation, BIO13-Precipitation of wettest month, BIO14-Precipitation of driest month. In addition to bioclimatic variables, we also included in the predictor set a variable representing the density of the hydrographic network (HYD; Ribeiro et al. 2020). For projecting our habitat suitability model on future conditions, we considered two Representative Concentration Pathways (RCP): RCP2.6 and RCP8.5, representing the minimum and the maximum emission pathways for the year 2050 respectively. For each RCP we used projections from four CMIP5 Global Circulation Models (GCM): CCSM4, HadGEM2-CC, IPSL-CM5A-LR and MPI-ESM-LR. Since our main interest was to project Bd habitat suitability and spread to future climatic conditions, we avoided to include land cover, vegetation and anthropogenic variables in the models (despite these variables may contribute to determine Bd occurrence: e.g. Liu et al. 2013;James et al. 2015;Bacigalupe et al., 2019), because they may experience an unpredictable and radical change in a Mediterranean environment by year 2050, and therefore their reliability for future projections is doubtful.
For what concerns habitat suitability model building, availability of environmental conditions was drawn from 10,000 random background points. Model performance was determined using the area under the receiver operating characteristic curve (AUC), which ranges from 0.5 for completely random models to 1.0, for perfectly predictive models, and considering that models with an AUC > 0.8 have a good predictive ability (Merow et al. 2013). Random sampling of 75% of the data were used to fit the model while the remaining 25% were used to evaluate model performance. This procedure has been repeated 20 times and AUC values and habitat suitability model predictions were averaged from these repeated runs (Merow et al. 2013;Bacigalupe et al. 2019). Finally, we projected our habitat suitability model to each one of the four potential future climatic conditions forecasted by selected GCMs, and then we averaged GCMs projection within the same RCP, hence obtaining two habitat suitability model projections for year 2050, for RCP2.6 and RCP8.5 1 3 respectively (Guisan et al. 2017;Salas et al. 2017). Finally, we compared the gains and losses of suitability between the current and future bioclimatic conditions.

Landscape connectivity models
Connectivity models for Bd were built using electric circuit theory (McRae et al. 2008). Electric circuit theory calculates electrical flow between nodes on a resistance surface, with the flow representing dispersal/diffusion pathways across various types of resistances across the landscape (McRae et al. 2008). In particular, animal/pathogen movements are modelled via random walk (i.e. a random process describing a path that consist of subsequent random steps on the space) across all possible movement pathways, assuming that the intensity of current flow between a couple of nodes is proportional to the number of times an individual walks the path under consideration (McRae et al. 2008). In other words, current flow equals the likelihood of movement across the landscape, and therefore it represents landscape connectivity (McRae et al. 2008;Dickson et al. 2019). Specifically, circuit theory analysis produces a map of electric current density, equivalent to the likelihood of use for the given path (McRae et al. 2008), by considering all possible pathways connecting these nodes on a resistance/conductance map.
Circuit-theory connectivity models have been typically built by using occurrence locations or habitat patches as electrical nodes: however, this approach may be restrictive if the objective of the study is to estimate landscape permeability in the whole study area (Koen et al. 2014;Pitman et al. 2017). To overcome this problem, a "wall-to-wall" or "omnidirectional" approach has been applied: by placing current source and sink nodes outside of the study area (Koen et al. 2014;Pelletier et al. 2014;Pitman et al. 2017) or by using a moving-window algorithm (McRae et al. 2016), thus obtaining a landscape-scale connectivity map which is independent of nodes placement. Given the sparseness of Bd data in Italy, we opted to use a "wall-to-wall" method, independent of nodes' placement.
We built circuit theory connectivity models using Omniscape connectivity algorithm implemented within Julia programming environment (McRae et al. 2016;Landau et al. 2021). This algorithm applies Circuitscape (McRae et al. 2013;Hall et al. 2021) iteratively across the entire map extent using a circular moving window of given radius, and using two raster maps as input: (i) the resistance map, defining the cost of crossing each pixel, and (ii) the source map, defining the amount of current to be injected in every pixel (McRae et al. 2016). Omniscape algorithm evaluates connectivity between every possible pair of pixels within the moving window, and then the resulting maps are summed to obtain a final map of cumulative current flow across the whole landscape.
We converted our habitat suitability models of Bd, for present and future conditions, to resistance maps (1-suitability) in order to model landscape permeability, an approach that has been successfully employed, among others, for frogs (Falaschi et al. 2018), martens (Balestrieri et al. 2019), leopards (Pitman et al. 2017), bears (Zeller et al. 2020) and elephants (Buchholtz et al. 2020). Given the sparseness of Bd data and considering that habitat suitability models also represent probability of presence of the modelled species in a given cell (Guisan et al. 2017), we also employed habitat suitability model values of every cell of the map as source values for the amount of current to be injected in the node. That is, if a cell has a suitability = 0.9, then a corresponding amount of current will be injected trough the corresponding node. We considered this approach to be more informative, rather than using few occurrence locations that are under-representative of the actual distribution of Bd across the study region (McRae et al. 2016), or by arbitrarily placing 1 3 nodes outside (or randomly inside) of the study area (Koen et al. 2014;Pitman et al. 2017). In our analysis, we set a radius for the moving window equal to 50 km and a block size of 3 (i.e. a parameter coarsening the source surface: it represents the size, in pixels, of each block composing the moving window, at the centre of which the source point is located).
By running Omniscape algorithm, three different outputs for every habitat suitability model (one for current conditions and two for future conditions, RCP2.6 and RCP8.5, respectively) were produced: (i) cumulative current flow (CCF), which is the total current flowing through the landscape, (ii) flow potential (FP), which represents a null model of current flow, assuming movement unconstrained by resistance (i.e. using a constant resistance map), and (iii) normalized current flow (NCF; calculated as CCF/FP) which helps to distinguish between areas where current is impeded (i.e. movement is contrasted by barriers; NCF < 1), diffuse (i.e. no resistance to movement; NCF = 1), intensified (NFC > 1) or channelled (i.e. movement is constrained between barriers and is greater than expected by the null model; NFC > > 1; McRae et al. 2016). Finally, by overlying national protected areas' borders to CCF and NFC maps, we calculated the mean current flow for each national protected area and identified those areas where the likelihood of Bd diffusion is higher for both present and future conditions (Littlefield et al. 2017;Choe and Thorne 2019).

Bd screening in Italy: field sampling and occurrence data
We found presence of Bd in 71 out of 1274 skin swabs (corresponding to 5.5%), and Bd infection occurred in 13 species out of 18 tested by qPCR (Table 1). The overall Bd prevalence in Italian amphibians was 6% (Bayesian 95% credible intervals 4-7%), with a high variation among species (Table 1). Although, the present research does not focus on individual Bd load, in general the Bd load obtained by qPCR ranged from 1 to 140 Bd GEs per swab, which may be considered as a relatively low infection load. From our screening, we obtained 23 original occurrence locations in different parts of peninsular Italy to be added to bibliographic data, thus yielding to a total of 56 Bd occurrence location, that were employed for building habitat suitability model (Supplementary Material Appendix 1, Table 1).

Habitat suitability models
The habitat suitability model built for Bd had a good predictive performance (AUC = 0.879 ± 0.021; mean and standard deviation, calculated over 20 replicated runs). Out of the 10 predictor variables considered, BIO6 (Min temperature of coldest month), BIO8 (Mean temperature of wettest quarter) and BIO13 (Precipitation of wettest month) had the highest relative contribution to the maxent model (39.4%; 14.1%; 8.3% of variable contribution, respectively). At the same time BIO6 is the variable with the highest gain when used in isolation and which determines the highest decrease in gain when omitted from the model. As regards the current bioclimatic conditions, Maxent habitat suitability model shows that maximum habitat suitability for Bd is predicted in the northernmost and southernmost part of the Apennine mountain chain, in Sicily and Sardinia islands, while in the central part of Italy and in the Alps the suitability seems markedly lower (Fig. 2).
Habitat suitability model projections on future conditions predict a loss of suitability for the most suitable areas actually identified with current climatic conditions, and conversely an increase of suitability in North-Western Alps and central Italy. This predicted change of suitability is markedly higher for the RCP8.5 with respect to RCP2.6 scenario (Fig. 2).

Landscape connectivity models
Cumulative current flow (CCF) and Normalized current flow (NCF), obtained from the application of Omniscape algorithm, are presented in Figs. 3 and 4, respectively, for both current and future conditions. Cumulative current flow, which is the total current flowing through the landscape, shows a higher connectivity in the northernmost and southernmost portions of the Apennine mountain chain, and for Sardinia and Sicily islands, while landscape connectivity seems to be reduced in the central part of peninsular Italy. For future conditions, CCF seems to be reduced in the northern and southern extremes of the Apennine mountain chain, with respect to the current condition model, but landscape connectivity increases in central Italy, in particular when considering the most severe climate change scenario (RCP8.5- Fig. 3). For what concerns NCF, the flow appears as unconstrained in the major part of the landscape (Fig. 4), while in some areas the flow is channelled by barriers that create some pinch points: e.g. in the southernmost Apennines, Sicily and Sardinia islands and in some Alpine valleys as well. For future conditions NCF seems to be subject to smaller changes with respect to CCF, and as a general rule the amount of landscape portions with channelled flow is reduced in both RCP2.6 and RCP2.8 scenarios, while diffused and intensified flow increases overall (Fig. 4). At current conditions, there are four national protected areas experiencing a potentially high level of Bd diffusion: three in the  Fig. 5). The current flow in the remaining national protected areas is lower or even negligible. Conversely, for what concerns future bioclimatic conditions, connectivity models for both RCP2.6 and RCP8.5 predict a decrease of potential diffusion in these same protected areas (Sila, Pollino, Aspromonte,  Golfo di Orosei/Gennargentu National Parks), while an increase of Bd flow for the Circeo National Park and Dolomiti Bellunesi National Park is expected (Fig. 5).

Discussion
In peninsular Italy, Bd is present all along the Apennine mountain chain and also in the Po plain, with an overall prevalence of 6%. When this value was compared with other European countries, with a reported sample larger than 1000 individuals, a similar prevalence was found in Germany (7%, in 3064 individuals) but a much higher one was observed in Spain (20%, in 1149 individuals; data from supplementary materials in Baláž et al. 2013). However, this latter high prevalence may be explained by an over-sampling of amphibians in areas with well-known Bd outbreaks, such as Central Spain and the island of Majorca (Bosch et al. 2001;Bosch and Martinez-Solano 2006;Walker et al. 2008;Baláž et al. 2013). Concerning the individual Bd load, the Italian data were in the lowest range of those reported by Baláž et al. (2013) that measured up to 4067 GEs, but usually much lower.
The application of species distribution and habitat suitability models to the chytrid amphibian fungus has been widespread in the last 15 years. Several studies modelled the current distribution and habitat suitability of Bd at a global (e.g. Ron 2005;Rödder et al. 2009;Lötters et al. 2009;Liu et al. 2013;Xie et al. 2016) or continental scale (e.g. James  Rahman et al. 2020;Zimkus et al. 2020); while other studies investigated Bd distribution at a smaller, regional scale, focusing on specific amphibian diversity hot-spots (e.g. Puschendorf et al. 2009;Seimon et al. 2015;Flechas et al. 2017;Miller et al. 2018;Bacigalupe et al. 2019). However, ecological relationships among hosts and parasites are complex and the outcomes of their interactions vary in association with global environmental variables but also with the complexity of the biological community at the local scale (Benício et al. 2019. Indeed, this relationship seems to be generally non-linear and that high biodiversity may dilute parasite occurrence (Halliday and Rohr 2019). This implies that highly diverse and rich ecosystems could inhibit the diffusion of wildlife diseases at the local scale. Therefore, the role of areas with complex amphibian communities could act as ecological barriers to Bd spread, and should become important areas for the monitoring of the chytrid fungus. Among the studies modelling Bd occurrence at a global scale, several identify Italy as a high-suitability/high-risk area for Bd (Rödder et al. 2009;Lötters et al. 2009;Liu et al. 2013). Within our study, we identified some specific area of high suitability for Bd in Italy (e.g. Sila, Pollino, Aspromonte National Parks), while other areas were predicted to be less suitable at the current climatic conditions (e.g. Foreste Casentinesi National Park). The major part of the studies inferring habitat suitability and distribution of Bd, highlighted an overwhelming effect of climate. For instance: Puschendorf et al. (2009) found that high temperature seems to constrain the distribution of the pathogen at small scale in Costa Rica, Flechas et al. (2017) identified mean temperature and precipitation seasonality as main drivers of Bd in Colombia, while Liu et al. (2013) observed a relationship with annual temperature range at a global scale. In this study, among the bioclimatic variables included in the habitat suitability model, we observed that extremely high or low temperatures, in particular in the wettest quarter of the year, (BIO6 and BIO8) were main predictors of Bd suitability and acted as limiting conditions for his occurrence. These results can be explained by both an increased rate of epidermal renewal, driven by higher temperatures, which may in turn reduce Bd infection (Piotrowski et al. 2004), or alternatively producing physiological stress, which may limit Bd reproductive success (Piotrowski et al. 2004). Besides climate, also vegetation (Liu et al. 2013), land cover (James et al. 2015) and anthropogenic factors have been found to shape the distribution of Bd (Liu et al. 2013;Bacigalupe et al. 2019).
Habitat suitability and distribution models also confirmed the link between epizootic chytridiomycosis and amphibian worldwide decline, highlighting how areas of rapid amphibian decline overlaps with those of higher suitability for Bd at a global scale (Rödder et al. 2009;Lötters et al. 2009;James et al. 2015). The majority of studies modelling Bd distribution only focused on current climatic conditions, while few studies also projected distribution and suitability models on future climatic scenarios, predicting that Bd may decrease globally in some regions by 2100, but with a shift towards higher latitude and altitudes (Xie et al. 2016). This trend has also been confirmed at a smaller spatial scale (e.g. Seimon et al. 2015;Xie et al. 2016;Miller et al. 2018). In the present study we also observed a general contraction of Bd suitable areas by 2050, considering both circulation pathways. In particular, we observed that the suitability loss will mainly occur in southern and coastal Italy, while suitability gain will be observed in mountainous areas, such as central Apennines, western Alps and eastern pre-Alps and Alps. Despite the variety of approaches employed, local studies identified two main issues: (i) the identification of high-risk and refuge area will be of primary importance, (ii) the inadequacy of local strategies to monitor and mitigate Bd expansion (e.g. Flechas et al. 2017;Bacigalupe et al. 2019;Rahman et al. 2020). In Italy, while we addressed the first issue by both collecting new occurrence data and modelling habitat suitability, we also tried to overcome the second problem by modelling diffusion pathways at the landscape scale. Applying this procedure, we identified the national protected areas that should be more involved in monitoring and prevention actions to mitigate Bd diffusion.
Connectivity models based on electric circuit theory, have been widely adopted to model the spread of diseases or pathogens in wildlife populations, such as deer chronic wasting disease (Nobert et al. 2016), or rabies in raccoon populations (Algeo et al. 2017). Notably, among the almost 300 studies using circuit theory biology, ecology and conservation science (Dickson et al. 2019), only one involved the amphibian chytrid pathogen Bd (i.e. Becker et al. 2017). By modelling forest connectivity between populations of an Hylid frog in Brazil, Becker et al. (2017) found that skin microbiome similarity and Bd load are related to landscape connectivity and natural vegetation gradient, but a landscape connectivity model of Bd spread is actually lacking. Therefore, our study, at least to our knowledge, is the first one employing circuit theory to model the diffusion of Bd at the landscape scale. By applying this method, we were able to identify four national protected areas (Pollino, Sila, Aspromonte and Golfo di Orosei and Gennargentu National Parks) that may experience maximum Bd diffusion for the current climatic conditions (i.e. Cumulative Current Flow). These areas of maximum diffusion are also identified by the presence of channelled or intensified movement as predicted by normalized current flow, meaning that in these areas the spread rate is potentially high (McRae et al. 2016).
For these areas, we suggest the development of an intensive molecular screening plan in order to track any possible change in Bd presence, population prevalence and individual load of resident amphibians. Furthermore, three out of four of these areas are in geographic and ecological connectivity (Pollino, Sila and Aspromonte National Parks), as resulted from our model, while the fourth area (Maddalena national Park) is in ecological connection with populations interested by a Bd induced mass-mortality that occurred in Sardinia Tessa et al. 2013). Indeed, according to the Italian national biodiversity strategy, national parks should become "focal points for research and monitoring networks … in terms of biodiversity" (MATTM 2010, page 31). Therefore, national parks identified as areas of high Bd diffusion will represent a fundamental tool to implement monitoring, awareness and mitigation strategies. Moreover, we suggest that a particular monitoring effort should be spent in those area where current and predicted Bd flow is channelled between climatic suitability barriers, within a national park or protected area, such as Dolomiti Bellunesi in the Alpine region and Circeo National Parks in the central region of the Apennine mountain chain. Finally, despite the fact that our study allowed the identification of high priority areas, in order to detect early Bd spread and diffusion; our results are obtained from an incomplete sampling of the national area. Therefore, the identification of National Parks (or other priority areas) involved in the monitoring network could benefit from a more intensive sampling effort or from the application of different techniques: such as the implementation of a national monitoring programme, relying upon environmental DNA sampling and Occupancy modelling .