The amphibian-killing fungus in a biodiversity hotspot: identifying and validating high-risk areas and refugia

Amphibian chytridiomycosis, due to infection with the fungus Batrachochytrium dendrobatidis (Bd), has been associated with the alarming decline and extinction crisis of amphibians worldwide. It is essential for conservation management to identify regions with high or low suitability for Bd. We use a species distribution model to estimate the environmental suitability of Bd in the Chilean Winter Rainfall–Valdivian Forest biodiversity hotspot. Fourteen environmental variables were used as predictors in the statistical modeling (Maxent, generalized linear models, random forest) which also included 56 independent Bd+ localities. High-risk areas (i.e., suitability above a defined threshold) were validated through prospective field surveys conducted in 2017. As results from Maxent, which only uses presence data, were the only results retained, refugia (i.e., suitability below a defined threshold) were validated with the Bd absences (N = 12) used in the GLM and RF modeling. Our results showed that (1) suitability for Bd increased with human footprint and with shorter distances to urban centers and water bodies and decreased with elevation; (2) climate was not a major factor shaping the current distribution of Bd; and (3) the model predicted high-risk and refugia areas fairly well. Surveys of 24 new localities in high-risk areas confirmed that 23 were Bd+; hence, these areas warrant consideration for long-term Bd surveillance, population monitoring, and disease mitigation. In addition, five localities with apparent Bd absence were found in the predicted high-risk areas. Our models showed that refugia can exist near high-risk areas and Bd+ sites. Four localities with apparent Bd absence were located within the refugia predicted by the model. Preventing Bd transmission to such refugia is of paramount importance for persistence of Bd-susceptible amphibian populations. The identification and validation through prospective field surveys of high-risk areas and refugia are imperative to develop strategies to prevent further arrival and establishment of Bd and also, by identifying amphibian species or populations of conservation concern in such areas, will help to guide specific actions to reduce the biodiversity loss caused by chytridiomycosis.


INTRODUCTION
Emerging infectious diseases (EIDs) are increasingly recognized as a major threat to biodiversity (Daszak et al. 2000. Furthermore, given the current scenario of rapid environmental change driven by anthropogenic activity, there is growing interest in understanding the role of global change in the emergence and spread of infectious diseases (Jones et al. 2008. This is particularly true for amphibians, where population declines and extinctions at the global scale reported in Australia, Costa Rica, and Chile have been linked to the EID chytridiomycosis, caused by the fungus Batrachochytrium dendrobatidis (Bd; Lips et al. 2006, Schloegel et al. 2006, Skerratt et al. 2007, Wake and Vredenburg 2008, Soto-Azat et al. 2013a, b, Valenzuela-S anchez et al. 2017. At the local scale, Bd prevalence (i.e. number of positive animals/number sampled) and the host response to Bd exposure are likely the result of the complex interaction between the host, the pathogen, and the environmental conditions they encounter (Murray and Skerratt 2012, Spitzen-Van Der Sluijs et al. 2014, Blaustein et al. 2018. Overall, our understanding of the factors underlying the distribution and prevalence of Bd has greatly increased over recent years (Berger et al. 2016). Globally, climate, elevation, and species richness have been shown to be associated with Bd occurrences (Olson et al. 2013, Xie et al. 2016). At regional levels, it is known that Bd prevalence varies with latitude (Kriger et al. 2007), climate (Murray et al. 2011), human activities (James et al. 2009, Adams et al. 2010, Liu et al. 2013a, and with the interaction of these factors (Bacigalupe et al. 2017). By analyzing Bd distribution and spread at regional scales, it is possible to identify localities that are likely to be highly suitable (high-risk) and those that are unsuitable for Bd. The latter might represent potential disease refugia for amphibians (Keppel and Wardell-Johnson 2012). Such information regarding Bd habitat suitability would be useful to conservation managers, especially if areas can be confidently classified as high-risk for, or as refugia from, this important threat. For example, Bd refugia may protect susceptible amphibians from the constant pressure of Bd spread, or potentially serve as source or destiny areas for amphibian conservation-based translocations.
Correlative species distribution models (SDMs) are extensively used tools that evaluate the association between the presence of a species in a set of locations and environmental data to predict its potential distribution in terms of habitat suitability across landscapes (Guisan and Thuiller 2005, Elith and Leathwick 2009, Estay et al. 2014. In this context, SDMs have been used to predict the global (R€ odder et al. 2009, Liu et al. 2013a,b, Xie et al. 2016, continental (Ron 2005, Murray et al. 2011, James et al. 2015, and regional (Puschendorf et al. 2009, Ghirardi et al. 2011, Moriguchi et al. 2015, Flechas et al. 2017 distribution of Bd. Although the value of Bd environmental suitability to local conservation managers has been disputed (Riley et al. 2013), its utility has also been demonstrated (Murray et al. 2011). For example, Murray et al. (2011) showed that areas of high-predicted environmental suitability for Bd in Australia (i.e., Queensland and New South Wales) were indeed areas where Bdrelated declines have been documented. Nevertheless, Riley et al. (2013) reported that infected populations of Crinia georgiana in southwestern Australia inhabiting areas with high-predicted environmental suitability for Bd show no evidence of declines. However, it should be noted that conclusions from Riley et al. (2013) came from studying Crinia georgiana, a non-threatened, non-declining species, which is also not susceptible to chytridiomycosis. Although SDMs are not without their shortcomings (Meineri et al. 2015, Qiao et al. 2015, they provide an important baseline of key information to guide future sampling efforts or to test specific hypotheses regarding disease risk (Murray et al. 2011).
Here, we use an SDM to estimate the environmental suitability for Bd of the Chilean Winter Rainfall-Valdivian Forest biodiversity hotspot (Mittermeier et al. 2011), which has high amphibian species richness and endemism, including a large number of amphibian species threatened with extinction (Vidal and D ıaz-P aez 2012). This hotspot incorporates the most urbanized, dense, and economically productive areas of Chile, as well as natural environments (Barbosa and Villagra 2015). We have recently shown species-specific patterns of Bd prevalence and a non-random geographic distribution of Bd throughout this hotspot; prevalence decreases with latitude and increases with economic development, interacting with some components of climate (Bacigalupe et al. 2017). Given the inherent complexity of the relationship between Bd, its host, and the environment, we used field surveys to validate the SDM-modeled refugia and the high-risk areas (i.e., suitability below and above defined thresholds, respectively). To the best of our knowledge, such ground-truthing of a Bd SDM has not been previously conducted. This novel approach for Bd-refugia identification offers a powerful management tool to sustain Bd-vulnerable amphibian host populations.

Bd occurrence
In order to construct an SDM for Bd in the Chilean biodiversity hotspot, Bd occurrences were compiled through an extensive survey (2008( , Bacigalupe et al. 2017) of the Mediterranean and Valdivian forest ecoregions, which together form the Chilean biodiversity hotspot (30°-46°S; Fig. 1, Mittermeier et al. 2011). In addition, an extensive literature search was carried out in www.bd-maps.net (Olson et al. 2013), Google Scholar, and ISI Web of Knowledge using the following search term combinations: Batrachochytrium dendrobatidis OR chytrid Ã AND Chile. Cited literature in all retrieved articles was screened for further relevant publications. Only studies that used molecular methods (PCR and qPCR) as a means of Bd detection were included. From these resources, we constructed a database of known Bd occurrence for the biodiversity hotspot (Appendix S1: Table S1).

Predictor variables
Fourteen environmental variables derived from landscape-scale geographic layers that could potentially influence the occurrence of Bd were used as predictors in the statistical modeling. These comprised bioclimate, elevation, distance to nearest water body, topographic wetness index, net primary production of biomass, amphibian species richness, and several human developmental variables (see Table 1 for a full description and for data sources). Layers of all predictor variables were cropped to our study area and resampled to 1-km resolution using ArcGIS 10.1. A principal component analysis of the bioclimate variables was conducted to reduce their high collinearity (i.e., the first 4 principal components accounted for 97.7% of the bioclimatic variability and thus were used as the final bioclimatic predictors; Appendix S1: Table S2).

Modeling methods
In order to increase the spatial independence between sites and reduce the effect of the heterogeneous distribution of the sampling effort, sites that were within 3-km of each other were grouped, with the centroid of each grouped site assigned as the new georeferenced location (Muths et al. 2008, Rohr et al. 2011. We considered a grouped site to be Bd positive if at least one individual showed evidence of Bd infection at any of the individual sites. The SDM was conducted using Maxent, generalized linear models (GLM, logit link), and random forest (RF) under default settings, with prevalence (i.e. the proportion of sampled sites where a species is present, Santika 2011) set to 0.8, as we have 56 Bd+ localities out of a total of 66 (see Results). Bd absences (N = 10) for GLM and RF were based in localities where more than 21 individuals were sampled and all were Bdnegative by qPCR. This minimum sample size was based on an expected infection prevalence of 13.3% (Bacigalupe et al. 2017), an assumed test sensitivity of 100%, and a confidence interval of 95%. For each modeling method, random sampling of 25% of the data was used to evaluate the model performance, while the remaining 75% was used to fit the SDM. This procedure was repeated 10 times, and the AUC and jackknife values reported represent the mean of the 10 tests (Phillips et al. 2017). In Maxent, a geographical sampling bias correction was applied, including a 50-km buffer around each surveyed site. The area outside this buffer was automatically excluded (Syfert et al. 2013). In Maxent, the gain is defined as the improvement in penalized average log-likelihood compared to a null model (Elith et al. 2011) and we used it to evaluate the relative importance of each predictor variable by (1) jackknifing the observed value of gain; (2) comparing the decrease in gain when fitting a model using all but one variable with the gain of ❖ www.esajournals.org  the full model (i.e., including all variables); and (3) comparing the increase in gain when fitting a model using only one variable with the gain of the full model. Models were evaluated using the area under the curve (AUC) of the receiver operating characteristic, which measures the ability of the model to discriminate between sites with Bd presence and background sites. All modeling analyses were performed using Maxent 3.3.4 (Phillips et al. 2017) and the R package biomod2 (version 3.3-7; Thuiller et al. 2016).

Identification of thresholds to determine high-risk areas and refugia
Most SDMs give the results of the occurrence probability or environmental suitability for the target species as a continuous prediction in the study area. However, for many conservation management applications, discrete or presence/ absence information is more useful (Liu et al. 2013a,b). In the context of our work, we considered the identification of two areas: (1) those with a high probability of Bd occurrence, which  (2009) Human footprint A composite index which characterizes human influence on the land based on accessibility, anthropogenic land use, population density, and infrastructure Sanderson et al. (2002) ❖ www.esajournals.org we defined as high-risk; and (2) those with a low probability of Bd occurrence, which we defined as refugia (Puschendorf et al. 2009). Thus, threshold values were needed to transform the continuous output of each model into a discrete representation. For high-risk areas, we tested the performance of maximizing the sum of specificity and sensitivity (mSSS) and the prevalence of the training model, which maximizes the discrimination between presences and absences (Phillips et al. 2017). These thresholds have been shown to perform very well for presence-only models (Liu et al. 2013a,b) as only results from Maxent were retained (see Results).
For refugia, we tested a very conservative threshold in order to be confident that an area identified with low suitability had a minimal probability of Bd occurrence. Therefore, we used the fixed cumulative value 1, which rejected only the lowest 1% of possible predicted values (Phillips et al. 2017).

Validation of SDM predictions
In addition to the SDM results, high-risk areas were validated via prospective field surveys conducted in 2017. These surveys collected 692 new Bd samples from 24 randomly chosen independent localities that were not included in the model training but were identified as high risk for Bd by our model. As we only retained the results from Maxent (see Results) which only uses presence data, we evaluated whether localities with apparent Bd absence (N = 12, 620 individual samples) were located within the refugia or high-risk areas predicted by the model. Two of these localities were found during the field surveys in 2017.

RESULTS
The compiled dataset for Bd occurrence comprised samples from 3881 individual amphibians (3539 from our field sampling, and 342 obtained from the literature) and represented 29 species (50% of Chilean amphibian richness), 20 of which had at least one individual positive for Bd (Appendix S1: Table S1). The dataset comprised 143 sites surveyed between 2005 and 2017, 85 (59%) of which had confirmed Bd presence. We used 56 independent Bd+ occurrence localities to carry out environmental niche models (Fig. 1).

Environmental niche modeling
Maxent had the highest average AUC compared with RF and GLM models (Maxent = 0.82; RF = 0.64; GLM = 0.61). As the predictive capacity of both RF and GLM was just above a random model (i.e., AUC = 0.5), the information they provide regarding the coefficients and the relative importance of the independent variables is unreliable (Qiao et al. 2015). Therefore, Maxent was the only modeling procedure retained for further analyses. The exclusion of human footprint caused the greatest reduction in both the regularized training gain and the test gain (Appendix S1: Fig. S1). The environmental variables with the highest gains when used in isolation were human footprint, elevation, distance to nearest water body, and distance to nearest urban center (Appendix S1: Fig. S1). Climatic variables were not primary drivers of Bd environmental suitability. The response curves for the predictors with the highest gains showed that environmental suitability for Bd increased with human footprint and shorter distances to urban centers and water bodies. In addition, environmental suitability for Bd decreased with elevation.

High-risk areas and refugia: Validation of SDM predictions
The threshold performance in the classification of training presences (true positives and false negatives), testing absences (true negatives, false positives), and testing presences (true positives, false negatives) for Bd is summarized in Fig. 2; Appendix S1: Table S3. In brief, a higher proportion of Bd+ sites were found from the prevalence of the training model threshold than from mSSS. In addition, compared to what was found with mSSS, prevalence of the training model threshold yielded a higher accuracy and a substantially lower omission error when classifying presences surveyed during 2017 ( Fig. 2; Appendix S1: Table S3). Therefore, the prevalence of the training model threshold was the one we used to define high-risk areas. Likewise, the fixed 1 threshold had zero omission error (Appendix S1: Table S3), thus providing a reliable indicator for Bd absence. Overall, the areas identified as refugia and high-risk for Bd covered 22.6 and 27.1% of the study area, respectively. Additionally, mean amphibian species richness in the high-risk ❖ www.esajournals.org areas was 6.1 AE 2.6 SD while in the refugia was 4.3 AE 3.1 SD. Despite this slight difference, both types of areas show a relatively high diversity of amphibians.
Twenty-three out of 24 sites surveyed during 2017 in the predicted high-risk areas were positive for Bd, which represents a 96% validation of this aspect of the model predictions (Fig. 1). The AUC from the Maxent model for these new occurrences was 0.815. Of the 12 sites with apparent Bd absences, four were located within the refugia predicted by the model (i.e., Chilo e Island) and five were located within the predicted high-risk areas (Fig. 1). The remaining three sites were located between the defined thresholds.

DISCUSSION
Amphibian chytridiomycosis is recognized as a primary driver factor of the global amphibian decline crisis (Skerratt et al. 2007, Wake and Vredenburg 2008, Collins and Crump 2009, Fisher et al. 2012. For effective conservation management, therefore, it is prudent to identify areas with high or low suitability for Bd. This knowledge can also inform where to allocate amphibian population monitoring and disease surveillance efforts (Bosch et al. 2015, Langwig et al. 2015. For the Chilean biodiversity hotspot, our results show that (1) suitability for Bd increased with human footprint and shorter distances to urban centers and water bodies and decreased with elevation; (2) climate is not a major factor shaping the current distribution of Bd, and (3) the model predicted high-risk and refugia areas fairly well, as they were validated with new data. A recent regional analysis of Bd prevalence in the Chilean biodiversity hotspot indicated that, once temporal and spatial effects are taken into account, Bd occurrence tends to be higher in regions with higher gross domestic product, particularly near large urban areas (Bacigalupe et al. 2017). Our SDM model found a positive association between Bd occurrence and proximity to urban centers and, more general, with increasing human footprint. At least two different processes might produce this pattern. First, it might be indicative of human-mediated transmission of the fungus (e.g., intentional or accidental introduction of Bd with infected animals or via the transportation of contaminated water or damp substrates; e.g., Garner 2007, Soto-Azat et al. 2016). This is consistent with other studies that have also found a positive relationship between Bd occurrence and human influence at both local (Adams et al. 2010) and regional (James et al. 2015) scales, but not at the global scale (Liu et al. 2013a,b). Second, this pattern might occur if the human footprint translates into reduced habitat quality for amphibians or increases other stressors, which may increase their susceptibility to Bd infections (Daszak et al. 2003, Gray et al. 2007, Ramsey et al. 2010, Becker et al. 2016, Blaustein et al. 2018, potentially increasing pathogen occurrence. Distribution patterns associated with populated centers and roads are common in the first stages of many biological invasions (Pauchard and Alaback 2004, Christen and Matlack 2009, Skultety and Matthews 2017. Therefore, our results might indicate that Bd invasion in Chile is fairly recent, a hypothesis supported by two additional lines of evidence. First, a retrospective epidemiological study based on museum specimens showed that the earliest record of Bd in Chile is from a four-eyed toad collected from A recent invasion, or invasions, of Bd into Chile would suggest that the distribution of Bd has not yet reached equilibrium with its environment. That is, an organism is considered to be at equilibrium if it occurs in all suitable locations and is absent from all unsuitable ones (Ara ujo and Pearson 2005). Thus, under equilibrium settings it would be reasonable to expect that climate (or other non-human environmental factors) were the main drivers of a species distribution. It is important mentioning that non-equilibrium settings violate one of the key assumptions of environmental niche modeling, because they usually involve extrapolating to novel environments, well beyond where the species is found (Jeschke andStrayer 2008, Elith andLeathwick 2009). However, this might not represent a problem in our study, as our data set has Bd occurrences fairly well distributed throughout the biodiversity hotspot (Fig. 1), encompassing the wide range of climates, environments, and human impact factors along its latitudinal range. Therefore, the extrapolated area is located within the same latitudinal range of the biodiversity hotspot, although beyond the 50-km buffer that we constructed around Bd occurrence locations. The thresholds for fixed accumulative value 1 (F1), prevalence of the training model (Prev), and for maximizing the sum of sensitivity and specificity (mSSS) used for refugia and high-risk delimitation are highlighted in blue (refugia) and red (high-risk) horizontal lines. Training presences were the Bd+ occurrences used to fit the model. Validation presences were the Bd+ occurrences obtained during the prospective fields surveys in 2017. Twenty-three of those occurrences are above the Prevalence threshold (high-risk areas), while the remaining 5 occurred in areas between thresholds. Absences with an Ã represent the Bd-localities obtained during the 2017 field surveys.
A recent Bd introduction could explain why climatic or other environmental factors do not play a dominant role predicting the current distribution of Bd in our study area. This result does not agree to what similar studies have shown for other regions (Puschendorf et al. 2009, Murray et al. 2011, Flechas et al. 2017. However, those studies have only used climatic variables as predictors, possibly because the Bd distribution in those regions is in equilibrium with the environment. In our study, only the PC2, which was positively associated with mean temperatures (i.e., annual mean, minimum temperature of coldest month, and mean temperature of coldest and wettest quarter, Appendix S1: Table S2), showed some test gain when used as the only predictor in fitting the model (Appendix S1: Fig. S1). Nevertheless, once the effect of the invasion stabilizes and Bd reaches an equilibrium with the environment, climatic factors would take more prominence in explaining its distribution and therefore in predicting high-risk and refugia areas. It was no completely surprising that precipitation variables, incorporated mostly within PC1, were not important in predicting suitability for Bd. This is because distance to water bodies encompasses a much more reliable proxy of the presence of permanent water, certainly an important component that facilitates Bd persistence and transmission. Our results also indicate that suitability for Bd decreases with elevation, which is opposite to what other authors have found regarding the best conditions for Bd growth Alford 2005, Scheele et al. 2017). Nevertheless, the high altitude areas where suitability is indeed very low (Andes mountains in blue in Fig. 1) are actually areas where human impact is very low, strengthening the fact of a recent invasion of Bd in the hotspot.

Management implications
Although correlative SDMs have been applied to a broad range of ecological questions, there has been considerable interest in predicting possible ranges of invasive species (Peterson 2003). Therefore, predicting the occurrence of Bd with accuracy represents a key tool for spatially explicit, landscape-based conservation planning.
Our surveys of 24 new localities in high-risk areas predicted by the model confirmed Bd occurrence in 23 of them. Identified high-risk areas were mainly situated along the coast and central valley of Chile (Fig. 1). While the prediction of a site being high-risk (i.e., high suitability for Bd) does not necessarily mean that Bd is already present, it identifies areas where Bd surveillance and amphibian population monitoring efforts are needed (Langwig et al. 2015). In addition, five localities with apparent Bd absence were found in the predicted high-risk areas (Fig. 1). The identification of na€ ıve populations of Bd-susceptible species in these areas should be a priority, either to develop management strategies to prevent Bd arrival or to rapidly identify novel Bd occurrences and consider Bd mitigations strategies. This is particularly relevant given that most of the identified high-risk areas are regions of high endemic amphibian diversity, such as the Nahuelbuta range (Vidal and D ıaz-P aez 2012). Although refugia are areas with low Bd suitability accordingly to our model, our results indicate that suitability in our system is more related with the potential of Bd to disperse into an area, mediated by human impact, rather than with environmental variables such as climate. This would suggest that the identified refugia are such because Bd has not yet reached them and therefore, preventing its eventual arrival is of paramount importance. In this context, refugia provide a critical line of defense because, at least in the Chilean biodiversity hotspot, they can exist not only nearby high-risk areas but also close to sites with known Bd presence (Fig. 1). Four localities with apparent Bd absence were located within refugia predicted by the model, all within southern Chilo e Island (Fig. 1). Low Bd suitability areas, therefore, represent potentially valuable pathogen refuge for native amphibians, as well as candidate sites for translocating animals, in the framework of reintroduction programs. Predicted refugia warrant systematic monitoring to identify potential Bd incursion and to enable rapid responses to prevent pathogen establishment (Garner et al. 2016, Hudson et al. 2016. However, as only 24.4% of the areas predicted to be refugia are under some formal aspect of land protection, education and outreach are key to engage the public about the importance of these sites and to ensure the protection of these important amphibian habitats. Furthermore, the lack of legislative and organizational governance to address emerging wildlife ❖ www.esajournals.org diseases might limit the range of potential actions that can be implemented. For example, a response strategy currently exists in the United States for the potential invasion of Batrachochytrium salamandrivorans (Bsal) aimed at identifying, among other things, the monitoring needs that could better inform management and conservation strategies for the country (Gray et al. 2015, Grant et al. 2016. In this context, it would be highly relevant to hold meetings and workshops with relevant stakeholders, including government agencies, private land owners (protected and non-protected areas), and local communities in order to design a monitoring plan and implement national preventive biosecurity protocols and campaigns such as the one employed with the invasion of the freshwater algae Didymosphenia geminata (Didymo) in south Chile (e.g., http://didymo.sernapesca.cl).
In this study, we developed a feasible approach for guiding conservation management to help mitigate the impact of the EID amphibian chytridiomycosis. Using Chile as a model system, we present key baseline information to enable the prioritization of efforts to detect and manage Bd in the Chilean biodiversity hotspot. We clearly identify which areas are priorities for monitoring for conservation purposes (e.g., refugia in southern Chilo e, Bd-negative high-risk areas) by using a fairly simple and well-known statistical methodology that can be readily applied elsewhere. The statistical methodology we used (i.e., SDM) is standard and has been already proposed as a tool to guide conservation efforts against chytridiomycosis (e.g., Puschendorf et al. 2009, Murray et al. 2011, Flechas et al. 2017. Although environmental suitability for Bd is context-dependent, SDMs could be used to identify areas of concern in other regions (particularly in Central and South America, two areas severely impacted by chytridiomycosis) where prioritized Bd surveillance efforts and amphibian population monitoring could be implemented. New information about Bd is rapidly emerging; hence, adaptive management of response efforts will likely be needed. Results from SDMs should be regularly updated once new data on Bd presence becomes available or as more precise and relevant information about environmental factors that enable or suppress Bd establishment is obtained.