Global Amphibian Extinction Risk Assessment for the Panzootic Chytrid Fungus

Species are being lost at increasing rates due to anthropogenic effects, leading to the recognition that we are witnessing the onset of a sixth mass extinction. Emerging infectious disease has been shown to increase species loss and any attempts to reduce extinction rates need to squarely confront this challenge. Here, we develop a procedure for identifying amphibian species that are most at risk from the effects of chytridiomycosis by 53 combining spatial analyses of key host life-history variables with the pathogen's predicted distribution. We apply our rule set to the known global diversity of amphibians in order to prioritize species that are most at risk of loss from disease emergence. This risk assessment shows where limited conservation funds are best deployed in order to prevent further loss of species by enabling ex situ amphibian salvage operations and focusing any potential disease mitigation projects.


Introduction
The IUCN Red List of Threatened Species reports one third of the about 6,500 extant amphibian species as threatened with extinction [1, for species list see ESI1].Chytridiomycosis, a disease caused by the amphibian chytrid fungus, Batrachochytrium dendrobatidis (Bd), plays a decisive role in this global biodiversity crisis [1][2][3][4][5][6] by driving rapid declines and species extinctions in pristine protected areas.The balance of evidence shows that Bd is spreading globally and, in response, this pathogen has been included as a notifiable disease in the Aquatic Animal Health Code of the World Organization for Animal Health (OIE) (Article 2.4.1.2;www.oie.int/fr/normes/fcode/fr_index.htm;latest access 8 July 2009).
An IUCN 'Amphibian Conservation Action Plan' has been developed [3,4] with its goal being the prevention of further loss of global amphibian biodiversity.This plan states that short-term ex situ breeding is the primary available conservation strategy for species under immediate threat from Bd to prevent further dramatic amphibian species loss owing to this panzootic.Antifungal bacteria reducing susceptibility of amphibians to Bd in nature may provide an additional solution to solve the amphibian crisis [7].Beyond this, any attempt to limit further spread of the pathogen requires intimate knowledge of the current distribution of Bd and identification of areas and amphibian species assemblages not currently affected by chytridiomycosis.There is a pressing need for tools to prioritize candidate species for targeted population monitoring, ex situ breeding as well as other treatments, to identify critical areas for the implementation of measures to limit the spread of the pathogen and, when possible, target disease mitigation.
Infection by Bd is expected to be particularly sensitive to environmental influences because the pathogen exclusively occurs in ectothermic hosts which are physiologically affected by ambient conditions [8].Furthermore, the pathogen is known to be particularly temperature and moisture dependent [e.g., 9,10].The in vitro growth optimum of Bd is at 17-25 °C, whereas temperatures higher than 29 °C, freezing and desiccation are lethal [11], findings that are supported by observations in the field [12][13][14][15][16]. Hence, the geographic extent of Bd's climatic niche can be readily assessed by Species Distribution Models (SDMs).These models give a prediction of potential distributions of species derived from their associated environmental parameters.In the case of pathogens, such models can provide predictions on impending epizootics in uninfected regions [e.g., 17,18].
Herein we undertake a worldwide risk assessment for the potential impact of Bd with the goal to identify amphibian species that are most at risk of future declines as a result of Bd invasion and infection.We did so by combining known amphibian species ranges with a prediction of the distribution of Bd based on a SDM integrated with the biological characteristics of host susceptibility.

Prediction of Bd Distribution
Climate information for Bd SDM building was obtained from Worldclim, version 1.4 [19], which is based on weather conditions recorded between 1950 and 2000 at spatial resolution of about 1 × 1 km².It was created by interpolation using a thin-plate smoothing spline of observed climate at weather stations using latitude, longitude and elevation as independent variables [20].The climate data set was obtained from the DIVA-GIS homepage (www.diva-gis.org;downloaded 15 May 2007), i.e., 36 monthly mean variables (each minimum and maximum temperature and precipitation, respectively).Based on these we calculated six 'bioclimate' variables for further processing with DIVA-GIS 5.4 [21]: 'annual mean temperature', 'maximum temperature of the warmest month', 'minimum temperature of the coldest month', 'annual precipitation', 'precipitation of the wettest month' and 'precipitation of the driest month'.Since it can be expected that Bd zoospores are able to survive in unfrozen microhabitat, water or on hosts during winter air freezing [11], subzero values in the 'minimum temperature of the coldest month' grids were pooled and set to 0 °C.More 'bioclimate' variables than the six used here can be obtained from Worldclim.However, to avoid model 'overfitting' [22] and multicollinearity of predictors, we restricted our study to the six 'bioclimate' variables mentioned which can be considered as biologically relevant parameters to Bd [10,11,14,16].Also, these variables have been suggested to perform well in SDMs [23][24][25] including those for Bd [26].We obtained 365 globally scattered Bd records (latitude/longitude) from www.spatialepidemiology.net/bd/ (accessed 15 August 2008).Records posted on this webpage were compiled from numerous scientific publications (a detailed list of references is also posted).Bd records were not randomly distributed over the world (Figure 1A), leaving the problem of possible sample selection bias which may violate SDM assumptions [27].To account for this, we extracted all 'bioclimatic' values at Bd records and performed a cluster analysis based on Euclidean distances, whereby resulting classes were blunted at a threshold leaving 200 classes.Only one record per class was used for further processing.This method reduces the amount of duplicate information in the feature space and thereby the impact of samples clumped in geographic space.Calculations were performed with XLSTAT 2008 (Addinsoft, www.xlstat.com;downloaded 1 July 2007).
GARP was used in an earlier model focusing on Bd's potential distribution in the New World [28].However, advances in SDM methodology [29] and current availability of Bd presence records enable modelling at a much finer spatial scale.Herein, Maxent 3.3.1 (www.cs.princeton.edu/~shapire/maxent;downloaded 25 May 2009) was used for Bd SDM calculation and mapping [30,31].Maxent has been developed within the machine learning community and implements a general purpose algorithm for making predictions and inferences from incomplete information.The Maxent algorithm estimates geographic distributions of species from locality point data and random background data by finding the maximum entropy distribution [30].Ideally, the area from which background data is obtained reflects those regions accessible to the target species [32].Therefore, we restricted the background samples to areas from which Bd was detected in the wild.Areas with climate conditions not analogous to those represented by background data may lead to uncertainties in model predictions.Maxent automatically allows an identification of the degree of uncertainty when projecting models ('clamping').In our SDM, the degree of 'clamping' was removed from the model prediction using the 'fade by clamping' option.In numerous comparative studies, Maxent has achieved better results than other presence only methods [summarized by 29,33].
It has been suggested that ensemble model predictions may enhance the reliability and robustness of SDM results [34].Therefore, we computed 100 models each trained with randomly chosen 75% of the 200 records for model training and subsequently integrated all results into a map indicating the average Maxent value per grid cell.The remaining 25% of the records were used for model evaluation through calculation of the Area Under the Curve (AUC) in Receiver Operating Characteristic curves [35], a threshold-independent index widely used in ecological modelling [35,36].In ROCs, the sensitivity values, the true-positive fraction against 1-specificity and the false positive fraction for all available probability thresholds are plotted [36,37].AUC values may range from 0.5 (random accuracy) to 1.0 (perfect discrimination).We received 'good' to 'excellent', following previously given definitions [37], AUC values, suggesting a high quality of our SDM output: mean AUC training data = 0.937; mean AUC test data = 0.910 (Figure 1A).Maxent allows to trace the relative contribution of each variable to the model.Herein, the 'annual mean temperature' had the highest explanatory power, followed by the 'maximum temperature of the warmest month' (Figure 1B).This pattern was consistent in all 100 models computed.The logistic output of Maxent chosen by us is a continuous, linear scaled map which allows fine distinctions to be made between the modelled probability of habitat suitability for Bd.Generally, the higher a logistic Maxent value the better the prediction and therefore the climatic suitability for a species under study.It has been proposed that this relationship can be directly related to a species' maximum possible abundance [38].Maxent calculates several threshold values at each run and values exceeding them may be interpreted as reasonable approximation of-in this case-Bd's potential distribution pending on the question at hand.We used the minimum training presence (mean = 0.049) as a strict threshold and the 10 percentile training presence (mean = 0.277; Figure 1C) as a more liberal threshold as suggested by Phillips et al. [30].In general, the mean Maxent value at the input records is typically 0.5 [30,31], which was therefore selected as third threshold as recommended by Liu et al. [39].Furthermore, the uppermost 25% of the logistic value was chosen as fourth threshold (logistic Maxent value of 0.75).

Identification of the Most Threatened Species
Geographic ranges of 6,156 amphibian species of all three orders were adopted from the IUCN Red List of Threatened Species, as categorized during the former IUCN/Conservation International/NatureServe 'Global Amphibian Assessment' ( [1], www.natureserve.org/getdata/amphibianmaps.jsp;accessed 7 July 2009; see ESI1, whereby taxonomy follows [40,41]).In order to assess Bd suitability within the geographic range of all amphibian species, we performed an overlay analysis of our Bd world map obtained from the SDM and the distributions of the 6,156 amphibian species.Higher logistic Maxent values at a given site are associated with a higher climatic suitability and most likely a higher maximum abundance for Bd [30,38].Since the mean minimum training presence and the lowest 10 percentile training presences in our model were 0.049 and 0.227 (Figure 1C) and in accordance with previous studies [39], we regard these as minimum thresholds for the environmental suitability to Bd.We calculated the total geographic range encompassed by each amphibian species and quantified the fraction suitable for Bd at four Maxent thresholds (>0.049, >0.227, >0.50, >0.75; see ESI1).Computations were performed with DIVA-GIS 5.4 and ESRI ArcMap 9.2.

Bd Risk Factor for Anuran Amphibians
In an approach to identify anuran species (from a total of 3,976) which exhibit highest probabilities of Bd-related decline or extinction, biological and life history information of Bd infected species showing rapid declines was used as explanative variables in a previously published study [42].Life history information comprised the degree of aquatic life-stage, the mean snout-vent length and the mean clutch size and geographic range size.Representative environmental values for each species were added using spatial datasets of altitude, annual actual evapotranspiration, net primary productivity, isothermality (a measure of annual temperature consistency), maximum temperature of the warmest month, precipitation seasonality, precipitation in the driest quarter and human population density [43][44][45].For model building, so called Generalized Linear Models specifying the link function as either 'logit' or complementary 'log-log' were applied, whereby the link with the lowest residual deviance was preferred [42].Expected effects of strong phylogenetic signals were removed using generalised estimating equations and Holm-Bonferroni corrections were made to account for Type I errors.According to the most predictive multi-predictor model (brier score = 0.06), species with aquatic life stages occupying small geographic ranges at high altitudes were most susceptible for Bd related declines.This final model was projected onto a larger data set, whereby a standardized value between 0 and 1 for each species was calculated, hereafter termed 'Biotic Index' (BI) (compare species list in ESI1 and see [42] for elaborate description of methods).
In order to generate a Bd 'Risk Factor' (RF) for anurans, we combined BI and the results of our quantification of species' distributions with high suitability to Bd (Maxent value >0.50, which is the mean score at the Bd records used for model building which was suggested as a suitable threshold [39]; ESI1) as: where area bd = the percentage of a species' distribution >0.50 Bd suitability.RF equally weights BI and area bd ranges between 0 and 1, whereby higher values reflect a higher threat balancing equally the fraction of the species area suitable to Bd and the BI.

Prediction of Bd Distribution
Our SDM suggests that highest suitability for Bd occurs in temperate and subtropical regions of both hemispheres, often near coasts; in the tropics, for instance, montane regions are identified to show a high suitability for Bd (Figure 2A).The potential distribution of Bd suggested by our SDM generally resembles that obtained through an earlier GARP approach [28] with fine scale patterns better depicted in our approach.
The currently known distribution of Bd records in the world is patchy (Figure 2A).Even when some of these gaps may be the result of limited collection efforts for Bd, as in Western Europe [4], it is remarkable that there are island-like high risk regions in which chytridiomycosis and its effects have not yet been observed [1,4].These include the Ethiopian Highlands, eastern Madagascar, the southern versant of the Himalaya, China's Yunnan province and considerable portions of tropical Asia (Figure 2A).The importance of these highly suitable locations being apparently free of the effects of the chytrid fungus is magnified considering that these regions harbour high levels of amphibian species richness and endemism (Figure 2B).
Evidence for extraordinary dispersal speed of Bd has been recorded in some regions [9] and apparently anthropogenic dissemination, especially through the international trade, has contributed decisively to the speed and extent of Bd's spread [46,47].Therefore, biosecurity measures and baseline surveys should be initiated at those areas highlighted through our analysis.

Figure 2. (A) Worldwide potential distribution of the amphibian chytrid fungus (Bd) and 200 out of 365 records (black dots) of this pathogen (see methods). The map is derived from a
Maxent SDM projected into geographic space based on six 'bioclimate' variables.Warmer colors are associated with higher Maxent values suggesting higher suitability for Bd, whereas grey areas are below the minimum training threshold and therefore considered to be unsuitable.This is equivalent to higher risk of invasion in regions from which Bd is currently absent.(B) Overlay highlighting worldwide regions of both high Bd suitability and high amphibian species richness.Warmer colors suggest higher overlap of Bd suitability (based on our SDM) and number of known amphibian species (see ESI1).(C) 'Risk Factor' (RF) for anuran amphibian species, calculated by combination of species life history traits linked with Bd-caused declines [42] and climatic suitability for Bd within the species' range (Appendices ESI1, ESI2).Warmer colors equal a higher proportion of species with a high RF per 0.5° grid cell.

Identification of the Most Threatened Species
The proportion of species within each amphibian order that exhibit a high overlap with areas that are environmentally suitable to Bd is summarized in Table 1 (for a detailed species list including results of the presented analysis, IUCN conservation status and Bd presence see ESI1).Approximately one sixth of all known amphibian species fall with their total distributions into regions potentially suitable to Bd (Maxent value >minimum observed training presence; Figure 1A) and more than 50% of all species exhibits an overlap of over half of their known geographic range with regions showing high Bd suitability.This number drastically increases to almost 6,000 (i.e., roughly all known amphibians) when over half of a species' range overlaps regions with a Maxent value >0.5 (Table 1).
Table 1.Number of species in the known amphibian orders which are most threatened by Bd through range overlap.Suitability for the pathogen is high at different percentages of their known geographic ranges.We here provide species numbers referable to environmentally suitable (minimum observed training presence) or highly suitable to Bd (values > 0.5) (Figure 2A).ESI1 provides detailed species-level information.

Bd Risk Factor for Anuran Amphibians
Most of the 833 anuran species which by their biology and life history show a high predicted susceptibility to Bd [42] occur in regions which at the same time are characterized by high suitability for Bd (listed in ESI1).Of these, we identified 379 species in which the entire geographic range is, in terms of climate, of high suitability to Bd (Maxent value > 0.5; for detailed list see ESI2).We consider these amphibians to be the most threatened by the emergence of Bd.As shown in Figure 2C, they are distributed all over the world including regions from where Bd is currently unknown (Figure 2A).
So far, little is known about the current infection or population status of most of the 'Top 379'.Perhaps due to the circumstance that many of them occur in regions from where Bd has not been recorded only seven of these species are reported to be infected with Bd in nature (list provided in ESI2).We suggest this is the result of limited surveillance for disease rather than the occurrence of healthy populations, as at least 42 species of the 'Top 379' (ESI2) have undergone so called 'rapid enigmatic declines' likely caused by the spread of Bd and the effect of chytridiomycosis [1,5].Generally, the 'Top 379' priority species outlined in ESI2 should be considered as the 'top of the iceberg'.The threshold of RF = 1 used for their identification was perhaps chosen in an arbitrary way since all species with RF > 0 may be affected in one or the other way.Available conservation funds should be invested for targeted population monitoring, compilation of life history and ecology information, and-if necessary-further in situ and ex situ efforts in descending order of RF, since the higher the RF value the higher is the potential threat caused by Bd.
The importance of emerging infectious diseases and anthropogenic dissemination of pathogens for biodiversity has in recent years been receiving increasing attention especially due to a suite of zoonotic diseases with impact on humans and economically important domestic animals.However, the integration of wildlife diseases in strategies for halting the loss of biodiversity is still remarkably premature [48,49].Of the 379 species identified under high risk of decline due to chytridiomycosis, 40% are categorized as 'Data Deficient' under the IUCN Red List of Threatened Species, whereas a compelling 94% of the species with sufficient information for proper assessment are categorized as threatened with extinction.However, it needs to be noted that species with restricted distributions and a combination of ongoing habitat threats would at the same time qualify as threatened based on IUCN's criteria B and as susceptible for Bd related declines according to the BI.Herein, we highlight the need to integrate novel approaches in the tool-set of conservation biology to mitigate the threat posed by pathogens.

Methodological Considerations
When interpreting our results, some issues related to the data and methods used herein need to be considered.Information on the distribution of the world's amphibians was adopted via the IUCN Red List of Threatened Species from the former 'Global Amphibian Assessment'.These data provide the most comprehensive and up to information currently available.Range information of each species is based on expert opinions; therefore the polygons describe in most cases the extent of occurrence and sometimes the area of occupancy of each species rather than providing an exact summary of existing populations [see 50].This may lead to an overestimation of the actually suitable areas in widely occurring species which may not be homogenously distributed throughout their general extent of occurrence (e.g., due to specific habitat requirements).Worldwide, these potential errors may not be homogenously distributed.Before conducting specific management actions, this should be acknowledged.However, since these polygons represent the best information available also used to derive the actual IUCN Red List status of amphibian species, they provide a suitable basis for our assessment.
Methodical inherent uncertainties can be assessed comparing multiple SDMs in ensembles of models providing more robust predictions [34].In our 100 models combined in the final model we observed rather low standard deviations of predictions per grid cells indicating low model variability and a generally good predictive ability (Figure 3).SDMs rely on specific assumptions including that environmental conditions are range limiting and that the range of the target species is in equilibrium with climate [51].Looking at the requirements of Bd regarding certain temperature and humidity conditions proven by both laboratory and field studies (see above), as well as the nearly world-wide distribution of Bd, these assumptions appear to be fairly met.Furthermore, SDMs commonly ignore biotic interactions which may lead to an exclusion of the modeled species from certain areas [52].Regarding Bd, the most obvious biotic interaction include the pathogens' relationship to amphibians.
(BI) and the 'Risk Factor' (RF) estimated in this paper.Furthermore species that were encompassed under 'rapid enigmatic declines' in the Global Amphibian Assessment 2004 (updated under the IUCN Red List 2009) and species in which Bd has been reported in the literature are pointed out.
ESM 2: List of most threatened species.The 379 species identified as most threatened through Bd, i.e.RF = 1.0 (BI = 1.0 and Maxent values > 0.50 at 100 % of their geographic range, see text).The 216 species listed as 'threatened with extinction' (i.e.IUCN Red List categories 'Vulnerable', 'Endangered' or 'Critically Endangered') are in bold while species that are too poorly known to asses the conservation status (i.e.categorised as 'Data Deficient' in the IUCN Red List) are italicised.A total number of 187 species that are assessed by the IUCN to be declining are marked with asterisks.Species in which Bd have been detected [55][56][57][58] are underlined.Entries are by family, following the taxonomy used by IUCN in the 2009 assessment [1, ESM 1].

Figure 1 .
Figure 1.Variation of AUC training and AUC test (A), and variable contribution (B), minimum and 10 percentile training presence (C) in 100 Maxent models.Abbreviations are: Bio1 = 'annual mean temperature'; Bio5 = 'maximum temperature of the warmest month'; Bio6 = 'minimum temperature of the coldest month'; Bio12 = 'mean annual precipitation'; Bio13 = 'maximum precipitation of the wettest month'; Bio14 = 'minimum precipitation of the driest month'.In boxplots, the minimum and maximum values are indicated as blue dots, 95% confidence intervals as short horizontal bars; the 1 st and 3 rd quartiles and the median are indicated as broad horizontal bars and means as red crosses.