Bioclimatic and anthropogenic variables shape the occurrence of Batrachochytrium dendrobatidis over a large latitudinal gradient

Amphibian chytridiomycosis, caused by the fungus Batrachochytrium dendrobatidis (Bd), has caused the greatest known loss of biodiversity due to an infectious disease. We used Bd infection data from quantitative real-time PCR (qPCR) assays of amphibian skin swabs collected across Chile during 2008–2018 to model Bd occurrence with the aim to determine bioclimatic and anthropogenic variables associated with Bd infection. Also, we used Bd presence/absence records to identify geographical Bd high-risk areas and compare Bd prevalence and infection loads between amphibian families, ecoregions, and host ecology. Data comprised 4155 Bd-specific qPCR assays from 162 locations across a latitudinal gradient of 3700 km (18º to 51ºS). Results showed a significant clustering of Bd associated with urban centres and anthropogenically highly disturbed ecosystems in central-south Chile. Both Bd prevalence and Bd infection loads were higher in aquatic than terrestrial amphibian species. Our model indicated positive associations of Bd prevalence with altitude, temperature, precipitation and human-modified landscapes. Also, we found that macroscale drivers, such as land use change and climate, shape the occurrence of Bd at the landscape level. Our study provides with new evidence that can improve the effectiveness of strategies to mitigate biodiversity loss due to amphibian chytridiomycosis.

For each species sampled, the proportion of swabbed animals that were Bd-positive and the median infection load are summarized in Table 1. The infection load in Bd-positive amphibians ranged from 0.1 to 630,720 ZE (Zoospore equivalents; median = 234; the highest load was obtained from a P. thaul individual). Despite the presence of some high infection loads, swabbed animals did not exhibit clinical signs consistent with chytridiomycosis. Of the total number of infected frogs, only 6.5% (774 individuals) had more than 10,000 ZE/swab.    Fig. 3A-C). Also, Bd infection was positively associated with anthropogenic biomes (95% CI 0.9-1; Fig. 3D), which means a higher Bd prevalence was observed in highly anthropogenic impacted ecosystems. Using model averaging, the anthropogenic biomes variable better explained Bd prevalence, and it was present in almost all the candidate models with low AICc scores (delta AICc < 4; Fig. 3D; Table S2). Also, the effect of ecoregions was significant, explaining differ-   Table 2). Half of the clusters, all of similar size, were located in the contiguous Metropolitan and Valparaíso administrative regions, the most densely populated area in the country. A large single cluster was identified further south, covering an area highly impacted by land-use change to exotic pine and eucalypt monocultures. The remaining three clusters were smaller in size and located near the cities of Angol, Valdivia and Puerto Montt, all in southern Chile. In the four spatial clusters in central Chile, the observed Bd prevalence was almost three times higher than would be expected by chance (Loglikelihood ratio range = 16.31-117. 73 and relative risk range = 2-3.7). The Global Moran's I index was statistically non-significant (p = 0.2).

Discussion
Amphibian chytridiomycosis is well recognized as a main causative factor of the current global amphibian extinction crisis 10,57 . Therefore, it is essential to identify risk factors facilitating the occurrence of Bd, as well as high risk areas of infection, which can provide guidance for effective conservation management 58 . Our results show that, at the landscape level across a large latitudinal gradient in Chile, Bd occurrence is: (i) biased towards certain families and species that use aquatic habitats; (ii) largely determined by bioclimate and human associated risk factors such as altitude, annual mean temperature, annual precipitation and anthropogenic biomes; and (iii) grouped in spatial clusters associated with urban centres. Overall, our data shows that Bd is widely distributed in Chile, infecting an extensive number of amphibian species and over broad altitudinal (0 to 3460 m) and latitudinal (18° to 48°S, comprising 3400 km) gradients. This included the finding of Bd infection in a Puerto Eden frog (Chaltenobatrachus grandisonae) near Villa O'Higgins (Fig. 1A), extending the previously known southernmost global record of Bd by 588 km further   59,60 . From the total of 40 sampled anuran species in Chile (64% of total richness) 49 , we found that 24 species showed evidence of Bd infection (see Table 1). All 11 species which had a Bd infection prevalence > 30% were exclusively aquatic amphibians, likely due to a higher contact rate with the infective stage of Bd in the aquatic environment 53,61 . Of these aquatic species, high Bd prevalences were found at the extreme north of the Andes in the altiplano frogs T. pefauri (3308 m) and T. chusmisensis (3365 m; see Table 1). Chytridiomycosis-related mortality has been described in the related species, T. pisanoi and T. atacamensis from northern Argentina 62 and T. marmoratus from Peru 17,63 . In addition, Bd-implicated disappearances of Telmatobius populations have been described in Peru 17 and Bolivia 64 .
The effects of Bd infection on anurans are highly variable and species-, population-or context-specific. For example, some species exhibit high disease-induced mortality, while others experience no detrimental individual or population effects of disease while maintaining enzootic infections 57,65 . In some cases, populations with cryptic, enzootic infections can experience chytridiomycosis-related mortality under certain environmental conditions, such as drought 41 . Vredenburg et al. 65 proposed a 10,000 ZE infection load rule, above which lethal disease invariably occurs. In our study, the species with highest Bd loads were C. gayi (median 436 ZE and a maximum of 409,440 ZE) and R. darwinii (median 1170 ZE and a maximum of 84,709 ZE; Table 1), possibly indicating impacts at the population level in these species. We recommend, therefore, that longitudinal population monitoring and Bd surveillance programmes be initiated or continued for these threatened species at Bd-positive sites (e.g., see Binational Conservation Strategy for Darwin's Frog 66 ). We found statistical differences in ZE between aquatic and terrestrial species, which emphasizes the role of aquatic species and aquatic environments in the maintenance and spread of Bd. Population declines due to Bd infection can occur in the absence of evident mass mortality or other obvious signs of disease, as evidenced by the Bd-driven declines of R. darwinii in Chile 53 and the possible extinction of its sister species, R. rufum, which has not been observed since 1981 48,52,66 .
Our best ranked model included the effects of altitude, annual mean temperature, annual precipitation, anthropogenic biomes and ecoregions, as predictors of Bd infection (Fig. 3). Anthropogenically-disturbed ecosystems proved to be one of the most important predictors that explain Bd infection. Also, altitude, annual mean temperature and annual precipitation were positively associated with Bd prevalence. Amphibian chytridiomycosis has been reported at high elevations, for example in the Rocky Mountains 67 , the Sierra Nevada 65 , the Pyrenees 68 and the high Andes 63 , suggesting that cold high-altitude environments do not necessary limit Bd spread and subsequent impacts on wild populations. The arrival of Bd in high altitude areas has been facilitated by human movement that has spread the fungus among isolated water bodies, but also climate change can facilitate such spread modifying the environment that anurans inhabit 63 and further force the severity of infection 35 . In addition to altitude, temperature and precipitation appear to be relevant climatic variables shaping the occurrence of Bd in Chile, as previously reported for other world regions 21,27,36,43,[69][70][71][72][73] . However, other studies have found sometimes different patterns (e.g., 39,41,74 ), suggesting that the mechanisms between these climatic factors and Bd occurrence are complex.
In our study, the highest prevalences of Bd infection were detected near to densely populated human settlements. Our results show that high human perturbation (anthropogenic biomes) is correlated with an increase in Bd infection probability, highlighting the importance of human activities on the epidemiology of Bd, possibly due to human-assisted pathogen introduction and spread (e.g., through the transport of Bd contaminated water and sediment) between anuran populations 24,57,75 . Additionally, human activities can spread Bd through the movement of infected amphibians, as has been shown to occur with amphibian trade, including the introduction of exotic amphibians 24,46,47,[76][77][78] . Also, it has been proposed that the reduced connectivity among amphibian populations resulting from human perturbation of the environment might impact host skin microbiome, affecting the innate immunity in amphibian skin against pathogens 79 . Habitat fragmentation can also affect ecological (e.g., colonization/extinction, host physiology, etc.) and evolutionary (e.g., local adaptation, evolution of resistance/ tolerance mechanisms, etc.) processes that can affect host-parasite interactions 19,79 . The development of increased susceptibility to infection through amphibian immunosuppression as a result of environmental contamination (e.g., pesticides) and habitat perturbation 80 are also possible impacts of human activities that increase the persistence and spread of Bd.
Our models showed that anthropogenic impacts and climate variables could synergistically interact and exacerbate infection risk (Table S2). The mechanisms enabling such synergy remain unclear, but our results support anthropogenic disturbance as a driver of Bd infection risk. In this context, anthropogenic biomes are generalizations for the restructured terrestrial biosphere due to agriculture, forestry and urbanization 81 . Elevated risk of Bd infection in areas close to human activities and settlements has been described previously in both temperate and tropical regions 23,25,36,43,75,82 . Most studies based species distribution models of Bd in the Americas have found an association of Bd occurrence with several climatic variables, notably precipitation, temperature and seasonality 13,21,32,33,44 , although few incorporate explicitly the effect of human impact, such as urban centres, in the analyses 43 . Interestingly, Zumbado-Ulate et al. 33 found a higher Bd occurrence in undisturbed ecosystems or protected areas, highlighting the fact that Bd occurrence is context specific as can be influenced by many factors including time of Bd introduction, species and population susceptibility, among others.
As with prevalence, the same association has been shown with intensity of infection, namely a positive association of Bd loads with anthropogenic disturbance 71 . The highest observed Bd prevalence was in the Chilean Matorral ecoregion, an area considered as a priority for global biodiversity conservation 83 . This region harbours a high level of anuran endemism 49 , yet contains the highest human population density in the country (almost 90% of the Chilean population). Consequently, increasing urbanization is resulting in deforestation and habitat loss which is negatively impacting amphibian populations 84 . Such environmental changes could favour the invasion of alien species, including Bd 24 .
All clusters of Bd-infection were located in central-south Chile, suggesting that amphibians in this region are at a higher risk of Bd-infection than elsewhere in the country. Although these findings are similar to those  43 , by using a different approach in a much wider area of Chile, this strengthens the hypothesis of urban centres playing an important role in the epidemiology of Bd. The four clusters located in the Metropolitan and Valparaíso regions could result from a potential initial introduction of Bd in this part of Chile and its subsequent spread 46 with fewer clusters in the south of Chile. When Bd is introduced to a new geographic location, first foci of introductions represented by narrow spatiotemporal cluster(s) occur followed by subsequent spread over time 68,85 . This has been seen in Spain, where Bd infection shows a pattern of introduction and spread along the Pyrenees with narrow spatial clusters, indicating recent introductions into Iberian biomes 35 . Our results are consistent with such a pattern having occurred in our study area. An alternative hypothesis is that instead of high-risk areas we might be capturing oversampled regions 86 . Therefore, we recommend considering other methods such as species distribution models 43 or kriging interpolation 33 to have a more accurate picture of the identified high-risk areas. Batrachochytrium dendrobatidis has been recognized as causing the greatest recorded loss of biodiversity due to a single pathogen 10 . Therefore, having a better understanding of the factors that shape the occurrence of Bd at the landscape level is valuable for conservation strategies and actions at the regional or country levels to halt the loss of biodiversity due to chytridiomycosis 14,87,88 . In this study, we provide evidence linking climate and anthropogenic factors as macroscale drivers of Bd occurrence at the landscape level. In particular the interaction of altitude, temperature, precipitation and human modified landscapes, appears to be the most relevant factors facilitating establishment and spread of Bd over large areas. Although we did not observe any mass mortalities or obvious signs of disease, we found high infection loads (> 10,000 ZE) which have been shown to lead to fatal chytridiomycosis in other amphibian species, with potential impacts at the population level. Although some species can coexist initially with high loads of Bd, host defence mechanisms such as anti-Bd skin microbiota, immunity, and skin peptides could reduce Bd infection loads to allow host-pathogen co-existence 89,90 . Strategies are required to prevent the further spread of Bd and to mitigate its impacts where it already occurs 10 . Our study provides valuable information for the design of such conservation strategies as long-term monitoring and control of biotic and abiotic factors (e.g., environmental disinfection, anti-Bd microbiome bioaugmentation or environmental management to reduce Bd exposure) 91,92 in areas with high occurrence of Bd and especially in areas where amphibian populations are under a high degree of anthropogenic pressure. Animal capture, biosecurity and sampling. Each amphibian was located through diurnal and nocturnal captures by direct observation and caught by hand or, in the case of aquatic species (i.e., Calyptocephalella gayi, Telmatobius spp. and X. laevis), caught using herpetology nets or funnel traps baited with chicken liver 94,95 . Following capture, each individual was handled for sampling with the use of clean disposable nitrile gloves and then released back to the exact point of capture. To minimize any false positive results and to avoid pathogen cross-contamination within or between study sites, a strict field sampling and disinfection protocol was followed 20 . For Bd detection, a non-invasive skin swab (MW100, Medical & Wire Equipment Co.) was obtained from each amphibian by firmly running it five times each over the ventral abdomen, the pelvis, both ventral hind limbs (femur and tibia), and the plantar surface of both hind feet, to complete a total of 35 strokes 52 . Swabs were kept in a refrigerated box until being stored frozen at -80 °C once back at the laboratory until they were analysed.

Ethic statement. This study was approved by Bioethics
Bd detection assay. Extraction of DNA from skin swabs and subsequent detection of Bd DNA using a specific real time qPCR assay was done 20 . For each sample, diagnostic assays were performed in duplicate, and standards of known zoospore concentration were included within each PCR plate as positive controls. We assumed that a Bd-positive swab indicated Bd infection. By including known concentrations of Bd DNA in serial diluted positive control (four standards of 100, 10, 1 and 0.1 Bd genomic equivalent) wells on each PCR plate, we were able to quantify infection intensity, which we defined as the number of zoospore equivalents/swab (ZE). To quantify and correct the infection intensity per swab, each genomic value was multiplied by 120 following Hudson et al. 96 .
Bd prevalence and infection loads by family, host ecology and ecoregion. We first calculated prevalence by counting the number of positive animals in a particular taxonomic family, host ecology (aquatic vs. terrestrial), or ecoregion divided by the total number of samples within that category. Host ecology was www.nature.com/scientificreports/ defined by considering whether the adult frogs of each species spent most of their time in or out of the water (see Table 1). We estimated 95% binomial confidence interval (95% CI) with a logistic (logit) parameterization for each category using the binom.confint function (R package 'binom') in the statistical software R v.3.1.3 97 . We evaluated whether there was a trend in Bd prevalence over time using a binomial generalized linear model (GLM) using year as an explanatory variable. The deviance of a null GLM model was estimated to explore the contribution of time (year) as an explanatory variable. An autocorrelation function 'acf ' was used to explore a potential temporal autocorrelation of the residuals. Finally, we applied odds ratio (OR) statistics to estimate the probability the amphibian to having Bd at each site in different years.
Modelling Bd prevalence across the landscape. We employed an information-theoretic modelling approach to contrast the adequacy of different working hypotheses explaining the geographic occurrence of Bd infection in our landscape gradient. In order to model Bd infection, we used bioclimatic and anthropogenic factors as explanatory variables and Bd prevalence from each of the 162 surveyed sites as response variable. Eight variables derived from landscape-scale geographic layers were used as predictors in the statistical modelling. Explanatory variables included annual mean temperature, temperature seasonality, annual precipitation, altitude, human footprint, anthropogenic biomes, ecoregions, and amphibian species richness (see Table S1 for a full description of each variable and data sources 81,[98][99][100][101][102][103] . We extracted all data for each sampled anuran to GPS coordinates using raster layers of 30 s (~ 1 km 2 ) spatial resolution with QuantumGIS v.3.8.2 104 . The assumption of normality of the data, the presence of outliers, the number of absolute zeros in the response variable and the collinearity between environmental variables and anthropogenic factors were explored. We only retained those variables with a correlation coefficient < 0.7 by the Pearson and Kendall tests. The homogeneity of variance was verified using the residuals of the model, plotting residuals vs. fitted values, and making a similar set of conditional boxplots for the residuals 105 . Twenty-one GLMs were constructed based on biological hypotheses using a binomial error structure (link = logit) for Bd infection model to evaluate factors that influence the occurrence of Bd infection at the individual level (i.e., Bd prevalence as response variable by study site). Then, we evaluated the degree of support for each competing model based on Akaike information criterion (AIC), which takes into consideration the likelihood of each model while penalizing for the number of included parameters to obtain the best-fitted model, and we used Akaike weights (AICw) to quantify the relative support of each model in the set (R package 'MuMIn'). As a result, a suitable model that describes the factors associated with Bd prevalence in anurans was obtained. To validate the GLM final model, the residuals were plotted against fitted values to assess homogeneity. In addition, we plotted a histogram and q-q plot of the residuals to verify normality of the data ( Figure S1). The existence of patterns in the residuals due to the assumption of independence was verified with a plot of the residuals against each of the explanatory variables. Finally, the effect of each significant variable in the GLM was displayed with R package 'effects' 106 . Bd spatial cluster analysis. Each sampled site was categorized as Bd-positive or Bd-negative according to the results of the Bd qPCR assays: a site was considered positive if at least one individual swab sample tested positive. Visualization of the sample sites was carried out using QuantumGIS and projected for analysis using the WGS 1984 datum as a coordinate system. Spatial distribution was characterized by the Moran's I spatial autocorrelation 107 , to identify spatial autocorrelation globally. We used Kulldorf´s clustering algorithm 108 under Bernoulli probability model, using the software SatScan v.9.4.4 109 to identify any cluster of Bd-positive samples across space with the proportion of infection at a given sample site. The model was run using Bd locations under the null hypothesis that cases were randomly distributed in space. The model was set to scan for areas with high Bd-positives numbers to test for clusters with a spatial occurrence higher than that outside the cluster. Briefly, the number of observed and expected Bd-infected amphibians is counted by a scanning window that moves across space for each location and variable window sizes 108 . Scan statistics allows the detection of the most "unusual" excess of observed Bd-positives and therefore provides georeferenced high-risk areas of Bd infection. Distributions of the likelihood ratio and its corresponding P-value were obtained using Monte Carlo simulation by generating 999 replications of the data set under the null hypothesis. The test statistics were computed for each replication and the test was deemed significant at P < 0.05.