Modeling the Distribution of the Chytrid Fungus Batrachochytrium dendrobatidis with Special Reference to Ukraine

Amphibians are the most threatened group of vertebrates. While habitat loss poses the greatest threat to amphibians, a spreading fungal disease caused by Batrachochytrium dendrobatidis Longcore, Pessier & D.K. Nichols 1999 (Bd) is seriously affecting an increasing number of species. Although Bd is widely prevalent, there are identifiable heterogeneities in the pathogen’s distribution that are linked to environmental parameters. Our objective was to identify conditions that affect the geographic distribution of this pathogen using species distribution models (SDMs) with a special focus on Eastern Europe. SDMs can help identify hotspots for future outbreaks of Bd but perhaps more importantly identify locations that may be environmental refuges (“coldspots”) from infection. In general, climate is considered a major factor driving amphibian disease dynamics, but temperature in particular has received increased attention. Here, 42 environmental raster layers containing data on climate, soil, and human impact were used. The mean annual temperature range (or ‘continentality’) was found to have the strongest constraint on the geographic distribution of this pathogen. The modeling allowed to distinguish presumable locations that may be environmental refuges from infection and set up a framework to guide future search (sampling) of chytridiomycosis in Eastern Europe.


Introduction
Amphibians are the most sensitive group of vertebrates-41% of currently known species are threatened with extinction [1][2][3]. Although habitat loss clearly poses the greatest threat to amphibians, a fungal disease threat (Batrachochytrium dendrobatidis Longcore, Pessier & D.K. Nichols 1999) has affected an increasing number of species [4]. This disease caused by the chytrid fungus B. dendrobatidis (Bd) has been linked to the declines in amphibian species globally and represents the greatest documented loss in biodiversity attributable to a pathogen [5]. A hypervirulent and globally distributed pandemic lineage (B. dendrobatidis-GPL) is considered the leading cause of population reduction in amphibians [6].
The Bd fungus has a wide geographical and host range in Europe; however, there currently is still a lack of complete understanding regarding the potential impacts of the amphibian communities. As of now, it has been registered in almost all Anura species in 26 European countries [7]. The susceptibility of various amphibian species across 2 of 18 Europe is still not fully understood, but declines in populations of species such as the common midwife toad (Alytes obstetricans (Laurenti, 1768)) and the common toad (Bufo bufo (Linnaeus, 1758)) [8]. In other species in which the fungus presence has been confirmed, infection with Bd does not always result in disease development, which could explain why mass die-offs of amphibians such as those observed in other parts of the world have not been seen [9]. Some species of anurans have limited immunity to Bd, while others such as water frogs (Pelophylax spp.) exhibit tolerance [10,11]. Some research also showed that populations of the yellow-bellied toad (Bombina variegata (Linnaeus 1758)) can coexist with Bd. However, the future safety of these toads may be compromised due to the potential effects of climate change and other environmental factors. The absence of widespread die-offs in many areas of Europe may be attributed to the fact that Bd has multiple strains with varying levels of virulence [12].
Although Bd is widely prevalent, there are identifiable heterogeneities in the pathogen's distribution that have been linked to environmental parameters [12]. In this respect, species distribution models (SDMs) have proven to be useful tools for predicting Bd distribution and elucidating the importance of a wide range of environmental covariates considered to affect Bd occurrence. Based on published data on B. dendrobatidis occurrence, SDMs can be used to evaluate the factors affecting its occurrence and predict its distribution [13][14][15]. The first developed Bd SDMs were global in scope [16,17].
It is well known that one major assumption of SDMs is that the species being modeled is in equilibrium with its environment. Considering invasive species, their distribution today may not be in equilibrium with current environmental conditions, and it may take centuries or millennia for invasive species to stabilize. However, global or other large-scale SDMs may provide overarching clues to Bd climate niche characteristics, and downscaled models such as those at the continental scale show that different climate metrics can be important predictors of Bd occurrence at finer spatial scales [17]. By considering geographically more restricted areas such as Ukraine or Eastern Europe, we anticipate obtaining a better view on Bd niche requirements that will aid us in reaching our objectives. Nevertheless, despite this reservation, habitat suitability models have been shown to be highly predictive [18,19].
Early detection and rapid response to incoming aliens are required for a successful management response [20]. Early-response strategies involve surveying risk areas under threat of invasion. SDMs help to identify such areas at risk that are suitable for the establishment of alien species, for instance by matching suitable climate conditions [21]. Identifying areas where a species is more likely to occur can also be used to guide sampling protocols and prioritizing areas of study [22]. Our objectives were to: (1) Identify priority survey areas in Eastern Europe (with a special focus on Ukraine) where future outbreaks of Bd could occur ("hotspots"), but perhaps what was equally important was the recognition of locations that may be environmental refuges ("coldspots") from infection (especially for amphibians that have certain sensitive conservation status) [23]; (2) Identify bioclimatic and other environmental conditions that constrain the geographic distribution of this pathogen in the study area.
It is necessary to mention that global or other large-scale SDMs may provide overarching clues to Bd climate niche characteristics, but downscaled models such as those at the continental scale show that different climate metrics can be important predictors of Bd occurrence at finer spatial scales. Hence, the scale of Ukraine or Eastern Europe is a much finer spatial scale than previously addressed by those recent large models, which makes this study exactly the next needed step.

Materials and Methods
As input, standard distribution models (SDMs) require georeferenced biodiversity observations. The globally dominant lineage is the pandemic one, while other lineages are of limited distribution or known transmission routes [24]. Therefore, this study did not take into account the lineage identity of positive samples for which there is yet insufficient data. Additionally, the lineages of Bd are not well mapped, resulting in a lack of necessary data for analysis. Localities for Bd were gathered from GBIF [25] and the literature [26][27][28][29][30][31][32][33][34][35]. Using quantitative PCR as described in [36], 20 new records of Bd from Belarus were added to the unpublished database of (A. A. Kulikova. To account for sampling bias, we used the nearest neighbor distance method ('ntbox' package in R; Osorio-Olvera et al., 2020) for thinning the data [37].
Because many limitations are associated with SDM projections, particularly when it comes to building an SDM for a species expanding its home range in a new area [38], we employed for the analysis only records of European localities. In addition, a global description of the niche does not account for the specificities of local invasive ranges (local environment, local biotic interactions, and specific human uses) [39]. Studies using records solely from the invasive range have demonstrated the ability to make adequate predictions regarding the expansion of invasive alien species [40][41][42].
SDMs commonly utilize associations between environmental variables and known species-occurrence records to identify environmental conditions within which populations can be maintained. SDMs extrapolate in situ habitats to identify geographic regions that have similar combinations of these values [43]. SDMs may be primarily climate-driven, meaning that the variables used to develop them typically portray climatic factors [44]. In this study, several types of environmental variables besides climatic at a geodetic resolution of 2.5 arc minutes were used as proxies for the fundamental niche [45].
First, 19 bioclimatic variables were downloaded from the WorldClim website (http://www.worldclim.com/version2 (accessed on 12 November 2022)) indicating a general trend of precipitation and temperature as well as extremity and seasonality of temperature [46]. However, we excluded four variables (bio8, bio9, bio18, and bio19) owing to their known spatial artifacts when following the protocol implemented in previous similar studies [47,48]. Other "quarter" variables were removed as well because they were highly correlated with monthly values and largely carried redundant information. Second, our model included a set of 16 climate and two topographic variables (the ENVIREM dataset downloaded from http://envirem.github.io (accessed on 20 January 2023); Table 1), many of which were likely to have direct relevance to ecological or physiological processes determining species distributions [49][50][51][52]. The included topographic variables were potentially important as well because they could modify the effect of the climate descriptors. Aspect cosine (min, max), aspect sine (min, max), eastness (min, max), eastness (min, max), northness (min, max), roughness (max, mean), slope (max, mean *), topographic position index (min, max), terrain ruggedness index (max, mean) Third, topography as measured by elevation and its derived variables (e.g., slope and aspect) was key for characterizing spatial heterogeneity and the abiotic environment in a given area, subsequently driving hydrological, geomorphological, and biological processes. All developed variables were available for download and visualization at the EarthEnv project site (http://www.earthenv.org/topography (accessed on 21 January 2023)) [52].
Fourth, metrics of land cover were included in our models [53] (Table 1). Land cover information offers a powerful first-order proxy for locally expected biodiversity and ecological processes [54]. Land cover is also considered relevant in models aimed at predicting species distributions because it adds realistic information on habitat fragmentation and human influence, which are not represented in more commonly used sets of climatic variables [55]. Consensus information on the prevalence of 12 land cover classes across the globe was downloaded from the EarthEnv project site (https://www.earthenv.org/landcover (accessed on 25 January 2023)) [56].
Finally, soil variables were included in our analyses. We assumed that the integration of soil factors in SDMs could help improve our understanding of factors that limit the Bd distribution range. These comprised 9 layers of soil physical and chemical properties recommended or previously used for building SDMs [57][58][59]. Soil grids were obtained from the Land-Atmosphere Interaction Research Group at Sun Yatsen University [60]. The values of the first five layers (0-49 cm) in the profile were derived and averaged for subsequent modeling.
Predictors often show high collinearity, and most SDM approaches require the selection of one among those strongly correlated [61]. In order to carry out such selection, the 'removeCollinearity' function in the 'virtualspecies' R (v. 3.3.3) package was employed [62]. This function analyzes the correlation among predictors, and using a 0.7 cut-off [63] returned a vector containing the names of those that were not colinear. It also grouped predictor variables according to their degree of collinearity, so from each such group consisting of two or more variables, those putatively most relevant to Bd could be selected. Because predictor variables are commonly skewed or have outliers, the Spearman correlation method was applied.
Today, a number of modeling approaches are at hand [64] depending on the environmental questions and available data; in particular, presences and absences. However, in the case of an emerging disease, the use of presence-only data is crucial to avoid possible complications due to false negatives and temporal variations in pathogen prevalence [65]. Therefore, a presence-only approach was chosen.
SDMs were generated by employing Bayesian additive regression trees (BARTs), a cutting-edge technique in this field. Running SDMs with BARTs has been substantially facilitated by the development of the R (v. 3.3.3) package 'embarcadero' [66], which is highly effective at identifying informative subsets of predictors. The package includes methods for generating and plotting response curves, illustrating the effect of selected variables on habitat suitability. The algorithm computes habitat suitability values ranging from 0 for a fully non-suitable habitat to 1 for a fully suitable habitat.
Then, 'Boruta', one of the most effective predictor selection algorithms implemented in R (v. 3.3.3) [67], was used to create a custom-made set of predictors for building a consensus model.
To prepare the input, we used pre-modeling functions from the 'flexsdm' R (v. 3.3.3) package [68]. The calibration area was defined using buffers around presence points. Because Bd is effective at dispersing, relatively large 500 km buffers were chosen [27]. Filtering the occurrence data was used to reduce sample bias by randomly removing points where they were dense (oversampling) in the environmental and geographical spaces.
The models were evaluated using the area under the receiver operating characteristic curve (AUC) [44] and the true skill statistic (TSS) [69]. AUC scores range from 0 to 1 (with 0 for systematically wrong model predictions and 1 for systematically perfect model predictions); AUC values 0.7 to 0.8 are considered acceptable, while values > 0.8 are considered to be good to excellent. TSS values range from −1 to +1 (with −1 corresponding to systematically wrong predictions and +1 to systematically correct predictions) [70].
TSS values <0.4 are considered poor, 0.4-0.8 useful, and >0.8 are good to excellent. Because AUC has its drawbacks [71], we employed the continuous Boyce index, which only requires presences and measures how much model predictions differ from random distribution of the observed presences across the prediction gradients [72]. Thus, it is an appropriate metric in the case of presence-only models. It is continuous and varies between −1 and +1. Positive values indicate a model that presents predictions that are consistent with the distribution of presences in the evaluation dataset, values close to zero mean that the model is not different from a random model, and negative values indicate counter predictions [73]. Evaluation metrics were calculated using the script posted by A.M. Barbosa on R-bloggers (https://www.r-bloggers.com/2022/05/model-evaluationwith-presence-points-and-raster-predictions/ (accessed on 7 January 2023)).
For guidance regarding the search for Bd and distinguishing areas in Ukraine where the pathogen most possibly might occur, the consensus distribution model was categorized into three frequency distribution classes (i.e., low, medium, and high) using Jenks' natural breaks classification, which (like k-means clustering) maximizes the variance between classes while minimizing the variance within classes [74].
We used a 50% habitat suitability threshold [75] as a cut-off above which responses of the most contributive variables could be analyzed in terms of their impact on habitat suitability.
Maps of habitat suitability in the GeoTIFF format were processed and visualized in SAGA GIS (v.2.14) [75]; statistical data was analyzed using the PAST software (v. 4.03) package [76] and/or the R environment (v. 3.3.3) [77].

Results and Discussion
The update of published and unpublished Bd-occurrence data yielded a total of 234 non-duplicate georeferenced records across Europe. The results of grouping predictor variables according to their degree of collinearity at a cutoff of 0.7 yielded a subset of metrics included in the analyses (Table 1). Variables selected for checking their relevance for Bd are marked with an asterisk.

SDMs Based on Selected WorldClim v.2 Predictors
Using the pre-modeling functions from the 'flexsdm' R package, the number of presence records was reduced to 114. The BART model showed a very good performance with a Boyce index of 0.951 and an AUC and TSS of 0.824 and 0.525, respectively.
The BART model indicated the importance of Temperature Seasonality, Annual Precipitation, and the Min. Temperature of Coldest Month (Figure 1) in predicting the occurrence of Bd. Temperature Seasonality is a measure (in units of standard deviation (SD)) of temperature variability over the course of a year. Figure 1A depicts that increasing the variability in temperature first established a hump-shaped relationship with habitat suitability within a certain interval (using a 50% threshold, that would be approximately between 4000 and 8000 units) ( Figure 1A). Further increasing variability produced a sharp negative effect. Notably, Temperature Seasonality was highly correlated with Temperature Annual Range, a surrogate for 'continentality' [78], and in the study area, both exhibited a strong longitudinal gradient with values increasing toward the east. Regarding Annual Precipitation ( Figure 1B), rising values were found as low precipitation values increased until reaching the mark of around 600 mm, after which the curve steadily approached precipitation values for which habitat suitability was the highest. Previously, Annual Precipitation was found along with Annual Mean Temperature to highly influence the distribution of Bd [24]. Other studies, in addition to stressing the importance of annual precipitation, reported the optimum rainfall for Bd to be between 1500 and 2500 mm a year [79,80] or at an initial modal maximum of around 1200-1400 mm realistic for the European study area [81]. In any event, the graph in Figure 1B leads to similar conclusions.

R PEER REVIEW 7 of 18
Annual Precipitation was found along with Annual Mean Temperature to highly influence the distribution of Bd [24]. Other studies, in addition to stressing the importance of annual precipitation, reported the optimum rainfall for Bd to be between 1500 and 2500 mm a year [79,80] or at an initial modal maximum of around 1200-1400 mm realistic for the European study area [81]. In any event, the graph in Figure 1B leads to similar conclusions. In terms of the Min. Temperature of Coldest Month (Figure 2A), habitat suitability remained low before reaching a 50% threshold in a period after −9 °C, which was well below the critical thermal minima of +4 °C [82]. In other places (for instance, the USA (WY, ME, CO, and CA), Bd localities reach the lowest coldest temperatures-down to −19.6 °C [74]. However, the pathogen is hardly exposed to such temperatures because Bd is strongly associated with aquatic habitats and host species that hibernate in aquatic microhabitats rather than terrestrial [83]. Therefore, the fungus is largely buffered from such extreme external conditions. In terms of the Min. Temperature of Coldest Month (Figure 2A), habitat suitability remained low before reaching a 50% threshold in a period after −9 • C, which was well below the critical thermal minima of +4 • C [82]. In other places (for instance, the USA (WY, ME, CO, and CA), Bd localities reach the lowest coldest temperatures-down to −19.6 • C [74]. However, the pathogen is hardly exposed to such temperatures because Bd is strongly associated with aquatic habitats and host species that hibernate in aquatic microhabitats rather than terrestrial [83]. Therefore, the fungus is largely buffered from such extreme external conditions.

SDMs Based on Selected ENVIREM Predictors
After performing the pre-modeling functions, the number of presence records was reduced to 161. The model showed a very good performance with a Boyce index of 0.971, AUC = 0.793, and TSS = 0.471.
The BART algorithm distinguished 'Continentality' and 'Annual potential evapotranspiration' as important predictors. 'Continentality' (or the same as Temperature Annual Range from the WorldClim v.2 dataset) is influenced by distance from oceans; i.e., it is a proxy of maritimity and actual Continentality [84]. The response curve for this variable presented in Figure 2B shows a general hump-shaped relationship with habitat suitability ( Figure 2B). Using the 50% habitat suitability threshold, suitable conditions for Bd were theoretically found between +17 and +24 °C, after which the response curve sharply declined to a negligible level. Maritimity below +17 °C, as found in much of the Atlantic biogeographical region of Europe [85], appeared unfavorable for the pathogen (likely because of the June-July temperatures), whereas toward the east (roughly beyond the longitude of 26° E) the limiting factor was freezing winter temperatures; this exhibited itself in a profound way as shown by the steeply declining response curve.
'Annual potential evapotranspiration' relates to the ability of the atmosphere to remove water through evapotranspiration processes and is strongly influenced by temperature [86]. It is regarded as an index meant to represent available environmental energies and ecosystem productivity [87]. Based on Figure 3A, it can be deduced that suitable conditions for Bd were apparently found (once again using the 50% threshold) between 460 and 900 mm/year, after which the response curve continued to rapidly decline ( Figure 3A). Places with high rates of 'Annual potential evapotranspiration' appeared to be less suitable for the pathogen. Using the 'contour lines' module in SAGA GIS, it could be shown that these areas were located primarily in Southern Europe (below a latitude of approxi-

SDMs Based on Selected ENVIREM Predictors
After performing the pre-modeling functions, the number of presence records was reduced to 161. The model showed a very good performance with a Boyce index of 0.971, AUC = 0.793, and TSS = 0.471.
The BART algorithm distinguished 'Continentality' and 'Annual potential evapotranspiration' as important predictors. 'Continentality' (or the same as Temperature Annual Range from the WorldClim v.2 dataset) is influenced by distance from oceans; i.e., it is a proxy of maritimity and actual Continentality [84]. The response curve for this variable presented in Figure 2B shows a general hump-shaped relationship with habitat suitability ( Figure 2B). Using the 50% habitat suitability threshold, suitable conditions for Bd were theoretically found between +17 and +24 • C, after which the response curve sharply declined to a negligible level. Maritimity below +17 • C, as found in much of the Atlantic biogeographical region of Europe [85], appeared unfavorable for the pathogen (likely because of the June-July temperatures), whereas toward the east (roughly beyond the longitude of 26 • E) the limiting factor was freezing winter temperatures; this exhibited itself in a profound way as shown by the steeply declining response curve.
'Annual potential evapotranspiration' relates to the ability of the atmosphere to remove water through evapotranspiration processes and is strongly influenced by temperature [86]. It is regarded as an index meant to represent available environmental energies and ecosystem productivity [87]. Based on Figure 3A, it can be deduced that suitable conditions for Bd were apparently found (once again using the 50% threshold) between 460 and 900 mm/year, after which the response curve continued to rapidly decline ( Figure 3A). Places with high rates of 'Annual potential evapotranspiration' appeared to be less suitable for the pathogen. Using the 'contour lines' module in SAGA GIS, it could be shown that these areas were located primarily in Southern Europe (below a latitude of approximately 47 • N). In a related manner, the minimum monthly potential evapotranspiration was found to be an important driving factor of spatial patterns of amphibian chytridiomycosis [88].

SDMs Based on Selected Topographical Variables from the EarthEnv Dataset
Filtering the number of presence records reduced them to 138. The BART mod showed an acceptable performance with the continuous Boyce index reaching 0.778 a an AUC and TSS of 0.766 and 0.396, respectively. The BART algorithm highlighted 'Slop as the most important predictor. It is worth noting here that elevation is amongst the top graphical variables repeatedly considered important for shaping the distribution of B Findings in this respect ranged from establishing evidence of a positive correlation of t pathogen occurrence with elevation [89,90] to the negation of any such correlation [91]. our case, it was notable that the BART algorithm fully excluded 'elevation' from t model.
K.M. Kriger and J.-M. Hero (2007) found a greater prevalence of Bd infection amo stream-breeding amphibians in Australia and suggested that dissemination of the path gen is greatly assisted by flowing water linked to slope gradient [92]. This could occ because flagellated zoospores of Bd rarely swim more than 2 cm prior to encysting [8 Besides that, the fungus prefers cool temperatures [85,92] and thus should grow better streams rather than in stagnant water bodies. In another study, the odds of being thre ened by Bd were found to be five times higher in stream microhabitats [93]. In our ca habitat suitability noticeably increased with 'slope' and reached a maximum at arou 25-26% ( Figure 3B). Interestingly, the results of a special study showed that runoff creased as slope gradient reached a critical point of 25%, then runoff decreased [94].

SDMs Based on Selected Topographical Variables from the EarthEnv Dataset
Filtering the number of presence records reduced them to 138. The BART model showed an acceptable performance with the continuous Boyce index reaching 0.778 and an AUC and TSS of 0.766 and 0.396, respectively. The BART algorithm highlighted 'Slope' as the most important predictor. It is worth noting here that elevation is amongst the topographical variables repeatedly considered important for shaping the distribution of Bd. Findings in this respect ranged from establishing evidence of a positive correlation of the pathogen occurrence with elevation [89,90] to the negation of any such correlation [91]. In our case, it was notable that the BART algorithm fully excluded 'elevation' from the model.  (2007) found a greater prevalence of Bd infection among stream-breeding amphibians in Australia and suggested that dissemination of the pathogen is greatly assisted by flowing water linked to slope gradient [92]. This could occur because flagellated zoospores of Bd rarely swim more than 2 cm prior to encysting [85]. Besides that, the fungus prefers cool temperatures [85,92] and thus should grow better in streams rather than in stagnant water bodies. In another study, the odds of being threatened by Bd were found to be five times higher in stream microhabitats [93]. In our case, habitat suitability noticeably increased with 'slope' and reached a maximum at around 25-26% ( Figure 3B). Interestingly, the results of a special study showed that runoff increased as slope gradient reached a critical point of 25%, then runoff decreased [94].

SDMs Based on Selected Land Cover Variables from the EarthEnv Dataset
After filtering, the number of presence records was reduced to 184. The BART model showed a good performance with the continuous Boyce index reaching 0.937 and an AUC and TSS of 0.832 and 0.536, respectively. The BART model strongly highlighted 'Open water' and 'Cultivated and Managed Vegetation' as important predictors (Figure 4). tween the percentage of land covered with cultivated and managed vegetation; in other words, crop-and farmland and habitat suitability ( Figure 4B). Impacts of land-use changes from increased agricultural production are commonly considered negative because they usually alter the habitat physically or chemically such that survival of resident organisms is questionable [95]. This can apply to both Bd and its amphibian hosts. For instance, pesticides can inhibit the immune response in amphibians, making them more susceptible to disease [96], but on the other hand certain fungicides are capable of reducing the number of Bd zoospores on frogs [97]. Nevertheless, initially the increase in the percentage of crop-and farmland in the landscape favored Bd habitat suitability, which apparently was mediated by the hosts while there was still a substantial amount of natural habitat. Amphibians have been found breeding in a variety of habitats that are substantially different from their former natural breeding habitats. In this respect, agricultural landscapes are of no exception [98], and the corresponding human infrastructure has been shown to provide beneficial environments to amphibian species [99,100]. Even so, a further increase of the percentage of crop-and farmland in the landscape (beyond approximately 60%) leads to a sharp drop in habitat suitability (assuming that heavily exploited agricultural areas are not fit for the pathogen).

SDMs Based on Selected Soil Feature Variables from the Land-Atmosphere Interaction Research Group Dataset
Filtering the number of presence records reduced them to 181. The BART model showed a good performance with the continuous Boyce index reaching 0.975 and an AUC and TSS of 0.795 and 0.468, respectively. In terms of top importance, the algorithm pointed toward the concentration of hydrogen ions (pH) ( Figure 5). Because life cycles of both the pathogen and its hosts are dependent on the availability of water, it was no surprise that the percentage of open water in the landscape was an influential factor. However, in excess it negatively impacted Bd habitat suitability. In this respect, an optimum was reached at around 28%, after which the response curve underwent a gradual decline ( Figure 4A). A similar hump-shaped relationship was found between the percentage of land covered with cultivated and managed vegetation; in other words, crop-and farmland and habitat suitability ( Figure 4B). Impacts of land-use changes from increased agricultural production are commonly considered negative because they usually alter the habitat physically or chemically such that survival of resident organisms is questionable [95]. This can apply to both Bd and its amphibian hosts. For instance, pesticides can inhibit the immune response in amphibians, making them more susceptible to disease [96], but on the other hand certain fungicides are capable of reducing the number of Bd zoospores on frogs [97].
Nevertheless, initially the increase in the percentage of crop-and farmland in the landscape favored Bd habitat suitability, which apparently was mediated by the hosts while there was still a substantial amount of natural habitat. Amphibians have been found breeding in a variety of habitats that are substantially different from their former natural breeding habitats. In this respect, agricultural landscapes are of no exception [98], and the corresponding human infrastructure has been shown to provide beneficial environments to amphibian species [99,100]. Even so, a further increase of the percentage of crop-and farmland in the landscape (beyond approximately 60%) leads to a sharp drop in habitat suitability (assuming that heavily exploited agricultural areas are not fit for the pathogen).

SDMs Based on Selected Soil Feature Variables from the Land-Atmosphere Interaction Research Group Dataset
Filtering the number of presence records reduced them to 181. The BART model showed a good performance with the continuous Boyce index reaching 0.975 and an AUC and TSS of 0.795 and 0.468, respectively. In terms of top importance, the algorithm pointed toward the concentration of hydrogen ions (pH) ( Figure 5). The partial response curve illustrating the relationship between acidity/basic habitat suitability ( Figure 5) indicated the best conditions were in the range of pH b 5.5 and 6.5 (with an optimum found around pH 6). This does not mean the pathog not occur in the field beyond this range. If it does, apparently this occurs under less able conditions. Experimentally isolates of Bd have been shown to have the most g at pH 6-7, less growth at pH 8, and minimal growth at pH 4 and 5 [82].
In summary, evaluation metrics for SDMs built on separate environmental d showed satisfactory results, meaning the applied predictor variables more or les captured habitat characteristics of the fungus species. Perhaps only the model ba topographic variables showed a relatively reduced performance; its evaluation m were all lower compared to the rest but nevertheless pointed out the importance of ' Influential predictors as assessed by the BART algorithm were pooled and su to a selection procedure using the R program 'Boruta'. In the end, the algorithm s the following variables: Annual Precipitation, Max. Temperature of Warmest Month tinentality, Gravel content, Organic carbon, PET seasonality, Evergreen/Dec Needleleaf Trees, Open Water, Deciduous Broadleaf Trees, Cultivated and Manage etation, and Urban/Built-up. Interestingly, roughly half of these were land cover va from the EarthEnv dataset. The final BART algorithm with the combined metrics re in two top predictors: Continentality and Cultivated and Managed Vegetation. The sponding SDM ( Figure 6) showed a pattern of greater Bd habitat suitability to th The partial response curve illustrating the relationship between acidity/basicity and habitat suitability ( Figure 5) indicated the best conditions were in the range of pH between 5.5 and 6.5 (with an optimum found around pH 6). This does not mean the pathogen cannot occur in the field beyond this range. If it does, apparently this occurs under less favorable conditions. Experimentally isolates of Bd have been shown to have the most growth at pH 6-7, less growth at pH 8, and minimal growth at pH 4 and 5 [82].
In summary, evaluation metrics for SDMs built on separate environmental datasets showed satisfactory results, meaning the applied predictor variables more or less fully captured habitat characteristics of the fungus species. Perhaps only the model based on topographic variables showed a relatively reduced performance; its evaluation metrics were all lower compared to the rest but nevertheless pointed out the importance of 'Slope'. Influential predictors as assessed by the BART algorithm were pooled and subjected to a selection procedure using the R program 'Boruta'. In the end, the algorithm selected the following variables: Annual Precipitation, Max. Temperature of Warmest Month, Continentality, Gravel content, Organic carbon, PET seasonality, Evergreen/Deciduous Needleleaf Trees, Open Water, Deciduous Broadleaf Trees, Cultivated and Managed Vegetation, and Urban/Built-up. Interestingly, roughly half of these were land cover variables from the EarthEnv dataset. The final BART algorithm with the combined metrics resulted in two top predictors: Continentality and Cultivated and Managed Vegetation. The corresponding SDM ( Figure 6) showed a pattern of greater Bd habitat suitability to the west and south of the area modeled.  A consensus distribution model (Figure 8) was categorized using Jenks' natura breaks classification into three frequency distribution classes: low (habitat suitability be The downscaled Bd habitat map for Ukraine (Figure 7) similarly showed the greatest Bd habitat suitability to the west with some suitable patches along the Dnipro River, in the Lower Danube area, and in the Crimea. Since the greatest diversity of amphibian species in Ukraine is observed in the Carpathians and forest regions [101], the confirmed presence of Bd tends to be a potential threat to a large proportion of the country's batrachofauna (17 of 19 species; excluding Salamandra salamandra (Linnaeus, 1758)). In addition, suitable patches of habitat for Bd were found along the Dnipro River and in wetlands in the Lower Danube area.
Sporadic areas in the Crimea could accommodate the fungus as well if it ever reaches the peninsula (perhaps via the Northern Crimean irrigation canal). A consensus distribution model (Figure 8) was categorized using Jenks' natural breaks classification into three frequency distribution classes: low (habitat suitability between 0.04 and 0.27), medium (between 0.27 and 0.58), and high (between 0.58 and 0.93). Maximum values of habitat suitability (>0.9) were observed in the following administrative regions (oblasts): Ivano-Frankivsk, Transcarpathia, Lviv, Chernivtsi, Ternopil, Volyn, and Rivne ( Figure 8). In summary, regarding the search for Bd, these are areas in Ukraine where the pathogen most possibly might occur, and focusing on them would reduce the areas for direct field samplings and facilitate the identification of potential "hotspots" as well as "coldspots".   Maximum values of habitat suitability (>0.9) were observed in the following administrative regions (oblasts): Ivano-Frankivsk, Transcarpathia, Lviv, Chernivtsi, Ternopil, Volyn, and Rivne ( Figure 8). In summary, regarding the search for Bd, these are areas in Ukraine where the pathogen most possibly might occur, and focusing on them would reduce the areas for direct field samplings and facilitate the identification of potential "hotspots" as well as "coldspots".

Conclusions
The continuing spread of Bd suggests that the geographic distribution of this pathogen is greater than currently known. Despite the apparent global invasion of Bd and a corresponding spate of past amphibian losses [84], there are many locations where this disease-causing pathogen has not yet been detected [102]. The use of broad environmental data to model the distribution of such a small fungal organism may cause some uncertainty depending on the geographical and study-variable scales; however, diverse research groups are using these tools to model pathogens' distributions [22,65,103]. Accordingly, ecological niche modeling (ENM) is an effective way to evaluate how these environmental factors affect current species distributions [104]. In this study, we focused on

Conclusions
The continuing spread of Bd suggests that the geographic distribution of this pathogen is greater than currently known. Despite the apparent global invasion of Bd and a corresponding spate of past amphibian losses [84], there are many locations where this diseasecausing pathogen has not yet been detected [102]. The use of broad environmental data to model the distribution of such a small fungal organism may cause some uncertainty depending on the geographical and study-variable scales; however, diverse research groups are using these tools to model pathogens' distributions [22,65,103]. Accordingly, ecological niche modeling (ENM) is an effective way to evaluate how these environmental factors affect current species distributions [104]. In this study, we focused on an attempt to highlight important variables shaping the current niche of a pathogenic organism using a variety of sets of bioclimatic and environmental predictors. The modeling algorithm in this case pointed toward 'Continentality' and 'Cultivated and Managed Vegetation' as comprehensive predictors of Bd distribution, therefore binding in such a way bioclimatic and human-induced factors. In this respect, unveiling macroscale environmental drivers of Bd and their interactions is crucial for proper conservation management of amphibians in the wake of an expanding disease [105,106]. In addition to the need to identify risk factors facilitating the occurrence of Bd, from a management perspective it is profoundly important to point out high risk areas of infection that potentially can favor the fungus. This is particularly substantial for Ukraine, where surveys for detecting the pathogen are yet to be undertaken.
The results of this study can inform the development of a strategic surveillance and monitoring program for Ukraine amphibian populations and associated threats (including Bd) as well as the development of biosecurity priorities to safeguard unique regions and taxa. Funding: This study was partly supported by Latvian Council of Science project "Ecological and socioeconomic thresholds as a basis for defining adaptive management triggers in Latvian pond aquaculture" (lzp-2021/1-0247) and the project "Modern risks of degradation of ecosystems of Ukraine on the example of model zoocenoses: analysis of factors in terms of biological safety" (project #0122U000708).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.