Infection with Batrachochytrium dendrobatidis is common in tropical lowland habitats: Implications for amphibian conservation

Abstract Numerous species of amphibians declined in Central America during the 1980s and 1990s. These declines mostly affected highland stream amphibians and have been primarily linked to chytridiomycosis, a deadly disease caused by the chytrid fungus Batrachochytrium dendrobatidis (Bd). Since then, the majority of field studies on Bd in the Tropics have been conducted in midland and highland environments (>800 m) mainly because the environmental conditions of mountain ranges match the range of ideal abiotic conditions for Bd in the laboratory. This unbalanced sampling has led researchers to largely overlook host–pathogen dynamics in lowlands, where other amphibian species declined during the same period. We conducted a survey testing for Bd in 47 species (n = 348) in four lowland sites in Costa Rica to identify local host–pathogen dynamics and to describe the abiotic environment of these sites. We detected Bd in three sampling sites and 70% of the surveyed species. We found evidence that lowland study sites exhibit enzootic dynamics with low infection intensity and moderate to high prevalence (55% overall prevalence). Additionally, we found evidence that every study site represents an independent climatic zone, where local climatic differences may explain variations in Bd disease dynamics. We recommend more detection surveys across lowlands and other sites that have been historically considered unsuitable for Bd occurrence. These data can be used to identify sites for potential disease outbreaks and amphibian rediscoveries.

efforts should prioritize actions on endangered taxa that are rapidly declining and the habitats that protect these species (Brooks et al., 2006;Foden et al., 2013;Giraudo & Arzamendia, 2018). However, there is often incomplete information on which populations are suffering the greatest declines and which locations provide them with the best chances of long-term persistence. For example, for several endangered species or clades, the majority of conservation actions have been designed based on opportunistic field studies conducted in sites where historic declines occurred (Kriger & Hero, 2007a). The potential bias caused by this unbalanced sampling might lead researchers to overestimate the rate of decline or to miss less dramatic declines and environmental threats across the range of the declining species. Therefore, extending the sampling to heterogeneous habitats across the entire geographic distribution of threatened species is crucial to detect and quantify potential threats as well as to establish suitable and more effective conservation actions (Hitchman, Mather, Smith, & Fencl, 2018;Miller et al., 2018;Olson et al., 2013).
Historic research on global amphibian population declines provides numerous examples of conservation actions in response to environmental threats in specific ecosystems. During the last four decades, at least 43% of described amphibian species declined or became extinct worldwide from multiple causes (Collins, 2010;Monastersky, 2014;Stuart et al., 2004;Wake & Vredenburg, 2008;Young et al., 2001). One widespread cause of amphibian population declines is the introduction of infectious pathogens. For example, Batrachochytrium dendrobatidis (Longcore, Pessier, & Nichols, 1999) (hereafter Bd) is a fungus that causes chytridiomycosis, a deadly cutaneous disease that affects amphibians in all continents where amphibians occur (Berger et al., 1998;Fisher, Garner, & Walker, 2009). Global assessments conservatively estimate that chytridiomycosis has caused the severe decline or extinction of over 200 species (Skerratt et al., 2007). Highland stream-dwelling amphibians have been hypothesized to be more prone to massive Bd-related die-offs than amphibians in other habitats (Hero, Williams, & Magnusson, 2005;Hirschfeld et al., 2016;Lips, 1998;Lips, Reeve, & Witters, 2003). Evidence suggests that tropical highland stream environments match the range of ideal abiotic conditions where Bd reproduces best in the laboratory (Berger et al., 2004;Longcore et al., 1999;Piotrowski, Annis, & Longcore, 2004). However, the spatial dynamics of Bd are intricate and still poorly understood. It is known that the intensity and occurrence of epizootic outbreaks and length of negative effects upon amphibian communities have varied globally (Catenazzi, 2015). In addition, numerous field studies show that prevalence and intensity of Bd infection vary with host species, microhabitat, temperature, humidity, seasonality, and geographic location (Kinney, Heemeyer, Pessier, & Lannoo, 2011;Kriger & Hero, 2007b;Kriger, Pereoglou, & Hero, 2007;Phillott et al., 2013;Searle, Gervasi et al., 2011). Thus, identifying conditions that constrain the geographic distribution of this pathogen will help elucidate why some species and populations suffer declines from Bd and identify locations that may be environmental refuges from infection (Murray et al., 2011;Rödder, Veith, & Lötters, 2008;Rosenblum et al., 2013).
The strong elevational gradients in the mountain ranges of Central America (Savage, 2002) create habitat heterogeneity and high endemism of amphibians in midlands and highlands (>800 m elevation). The cool and moist environments in tropical highlands provide suitable conditions for the Bd epizootic that occurred in Central America during the 1980s and 1990s, causing the extinction of an unknown number of amphibian species, especially highland stream-breeding species (Cheng, Rovito, Wake, & Vredenburg, 2011;Lips, Diffendorfer, Mendelson, & Sears, 2008;Pounds et al., 2006;Pounds & Crump, 1994;Rovito, Parra-Olea, Vasquez-Almazan, Papenfuss, & Wake, 2009). Historical declines in montane amphibian species reflect why most studies on amphibian host-Bd dynamics in the tropics have been conducted in premontane and upper elevation localities (Lips, 1999(Lips, ,1998Puschendorf, Bolaños, & Chaves, 2006;Ryan, Lips, & Eichholz, 2008). For example, a considerable amount of Bd infection data has been opportunistically collected from montane ecosystems, increasing the focus of conservation actions on highlands while overlooking other potential environments where amphibians may also be impacted by Bd (Puschendorf, Hodgson, Alford, Skerratt, & VanDerWal, 2013). For example, the suitability of lowland ecosystems for the spread of Bd has been frequently disregarded (Puschendorf et al., 2009) even though is known that some amphibian species (Figure 1) and clades have suffered dramatic unexplained declines in these zones (Chaves et al., 2014;La Marca et al., 2005;Puschendorf et al., 2009;Ryan et al., 2008;Whitfield et al., 2007;Zumbado-Ulate, Bolaños, Gutiérrez-Espeleta, & Puschendorf, 2014).
Thus, combining information on infection prevalence and abiotic conditions (e.g., from the WorldClim dataset) across the entire geographic distribution of a host can provide a more informative distribution of both the host and pathogen to identify potential hotspots of future disease outbreaks and potential environmental refuges from disease (Green, 2017;James et al., 2015;Rödder et al., 2010).
In this study, we sampled for Bd at four tropical lowland locations in Costa Rica and contrasted Bd prevalence and intensity of infection across study sites. We hypothesized that different host-pathogen dynamics occur across study sites because they exhibit latitudinal and altitudinal variation (Kriger & Hero, 2008;Kriger et al., 2007).
We extracted all 19 bioclimatic variables of the WorldClim to describe the different ranges of temperature and precipitation across study sites, which are the main environmental variables that affect Bd growth and dispersal (Nowakowski et al., 2016;Savage, Zamudio, & Sredl, 2011). Additionally, we hypothesized that all study sites would exhibit low levels of Bd prevalence and intensity of infection suggesting stable enzootic infections of Bd (Retallick, McCallum, & Speare, 2004;Scheele, Hunter, Brannelly, Skerratt, & Driscoll, 2017;Woodhams et al., 2008). Finally, we also expected a higher prevalence of Bd in amphibian assemblages occurring in permanent streams than in ephemeral ponds and terrestrial assemblages, as has been found in previous studies (Kriger & Hero, 2007a;Lips et al., 2003).

| Lowland sampling sites
We sampled four assemblages of amphibians between November and December 2011, at four tropical lowland locations in Costa  (Holdridge, 1967). Study sites consisted mostly of tropical moist forest and tropical wet forest with transitional ecosystems including semi-deciduous and evergreen forests, with temperature and precipitation ranges characteristic of these life zones. Our four sampling sites grouped into two main zones:

| Pacific zone
Here, we focused on the areas surrounding the small towns of

| Pathogen detection
At each site, four people systematically searched for amphibians for 36-48 hr during the day and night (9-12 hr/person). Within each site, we conducted visual encounter surveys of amphibians (Heyer, Donnelly, McDiarmid, Hayek, & Foster, 1994) and classified them by the habitat where they were captured: stream-dwellers (permanent flowing water), pond-dwellers (standing ephemeral waterbodies such as swamps, pools, and ditches), and forest-dwellers (leaf-litter, tree holes, or bromeliad plants in the understory, and canopy).
Caught amphibians were stored individually in clean, unused plastic bags. Each individual was inspected for visible signs of chytridiomycosis, such as hyperplasia, hyperkeratosis, abnormal shedding, depigmentation, and lethargic behavior (Berger et al., 1998;Voyles et al., 2009) and swabbed to detect Bd with a cotton swab (Medical Wire and Equipment, MW-113) using nitrile gloves. To swab, we ran a total of 20 strokes on every individual as follows: five strokes on one hand, five strokes on the ventral patch, five strokes on one foot, and five strokes along inner thigh. Swabs were stored dry in 1.5 ml Eppendorf tubes and frozen at −20°C until DNA extraction. All amphibians were immediately released after sampling. During this study, we followed field protocols (Kriger, Hines, Hyatt, Boyle, & Hero, 2006;Skerratt et al., 2008) which were approved by the National System of Conservation Areas of Costa Rica (SINAC, research permit 001-2012-SINAC) which ensures that animals are being cared for in accordance with standard protocols and treated in an ethical manner.
We extracted DNA from swabs using PrepMan Ultra (Boyle, Boyle, Olsen, Morgan, & Hyatt, 2004). All extractions were diluted 1:10 in 0.25X TE buffer and run in singlicate (Kriger, Hero, & Ashton, 2006) following diagnostic quantitative PCR (qPCR) standard protocols (Boyle et al., 2004) using an Applied BioSystems Prism 7300 Sequence Detection System to test for the presence and quantity of Bd genome equivalents. All Bd-positive samples were run again in singlicate confirmatory assay. Negative controls (DNase/RNase-free distilled water) were run in triplicate on every 96-well PCR plate.

| Data analysis
We were interested in understanding how Bd prevalence and intensity varied among our study sites and habitats (predictor variables).
For our analyses, we pooled all species together instead of using species as predictor or running independent tests for each species because the samples sizes per species were highly variable (from 1-44). This high variance in the sample size could produce significant models that may be an artifact of opportunistic sampling instead of a real pattern. Therefore, we analyzed habitat as a proxy of amphibian community composition, because the species variable was 100% correlated to habitat. To contrast Bd prevalence, we used fix-effects generalized linear models (GLMs) to find the most suitable model using binomial response variables (infected or not infected). Candidate models were ranked according to the Akaike's information criterion (AIC) to determine the relative importance of predictor variables within each model set. The model with the lowest AIC was considered the most robust (Burnham & Anderson, 2004). To compare infection intensity among locations and habitats (predictors), we generated fix-effects general linear models (LMs) with data only from infected individuals. We built our models using the log-transformed Bd load (estimated number of genomic equivalents) as a response variable and included site and habitat as predictors. Candidate models were ranked according to the coefficient of regression (R 2 ), with the model with the highest R 2 considered the most robust (Zar, 2013). For the most robust GLM, we tested the significance of the predictors using an ANOVA with a chi-square approximation to find the probabilities of predictor variables within the most suitable models, and for the most robust LM we used an ANOVA. Finally, we conducted post hoc, pairwise comparisons (Tukey test) to confirm where the differences occurred between significant predictors.
To describe the local abiotic environment for the sampled lowland sites, we generated buffers (radius = 10km) around each one of our four study sites. Because we wanted to achieve a full description of the abiotic environment, we extracted values for all the cells occurring within each buffer (mean = 355 cells/site, Table 1) from all 19 bioclimatic variables of WorldClim (version 1.4; www.worldclim.org) at a spatial resolution of 30 arc-s (Hijmans et al., 2005).
We compared the abiotic environment among sites using a principal component analysis (PCA). To contrast climatic dissimilarities between lowland study sites, we also generated a pairwise matrix of Euclidean distances between the centroids of climatic envelopes. All analyses were conducted in R v.3.5.1 (R Core Team, 2014).

| RE SULTS
We screened a total of 348 adult amphibians from 47 species for Bd (346 frogs and two salamanders, Table 2). From this list, a total of 44 species are classified as least concern and three are categorized as threatened: Oophaga granulifera is classified as vulnerable (VU), tested positive for Bd and total prevalence of Bd was 54.6%. We did not detect Bd on three of the amphibian families sampled, including Plethodontidae, the only family of Salamanders in the Neotropics; however, the sample size for these families was very small.
Prevalence of infection showed high heterogeneity among sites with values ranging from 0.0% in Rincon de Osa to 68.6% in Punta Banco (Table 3). This variation in Bd prevalence was best explained by the interaction effects model (Table 4), which showed significant effects of locality (p < 0.01), and that the variation of Bd prevalence by site depends on the habitat (p < 0.001; Figure 3a). Despite being close in proximity, amphibian assemblages from Sarapiqui showed significant higher prevalence of Bd than assemblages from Siquirres (p < 0.01, Figure 3a, Tables 3 and 5). There was a nonsignificant trend for Bd prevalence to be lower in Siquirres than Punta Banco (p = 0.06). We also found high prevalence of Bd across habitats (Table 3), but no significant differences between habitats in our Similarly, the differences in the infection intensity across study sites ( Figure 3b, Table 3) were best explained by the interaction model (R 2 = 0.19, Table 4), which also showed significant effects of location (F 2,166 = 15.5, p < 0.001) and the interaction between habitat and location (F 3,166 = 3.6, p < 0.01). Levels of infection intensity were significantly lower in Sarapiqui (Figure 3b, Tables 3 and 5) when compared to Punta Banco (p < 0.001) and Siquirres (p < 0.01). Overall, the infection intensity ranged from 0.1 to 63,861 genome equivalents and four individuals had more than 10,000 zoospore genomic equivalents, a theoretical number that is considered a threshold that results in mass mortality and rapid population decline . However, none of the sampled individuals including the TA B L E 1 Mean values (standard deviation) of the 19 bioclimatic variables from the WorldClim dataset and loads (coordinates) for PCA axes 1 and 2 showing the specific contribution of each of the bioclimatic variables used in the environmental analysis of four lowland sites in Costa Rica   four that were heavily infected showed any evident signs of disease.
Remarkably, three of these heavily infected individuals belong to the Critically Endangered species Craugastor taurus.
In our PCA analysis of the 19 bioclimatic variables, we retained the first two axes (Table 1) because they accounted for 98% of the total variance of our data. A tridimensional representation of PCA axes 1 and 2 (PCA 3 included as reference) shows four separated clusters of points, each one representing a study site ( Figure 4). As expected, we found the highest similarity in climatic conditions occurred among sites in each zone (Appendix A). We found that bioclimatic variables associated with precipitation (Annual Precipitation, Precipitation of Wettest Quarter, Precipitation of Driest Quarter) make a higher contribution in the variance of our climatic data than other variables (Table 1).

| D ISCUSS I ON
We found Bd infections at three of the four lowland sites sampled and in 70.2% of the 47 sampled species for an overall Bd prevalence of 54.6% (Tables 2 and 3). Furthermore, we did not detect signs of disease in heavily infected individuals during the study and found low levels of infection in most of our samples. Similar community composition and population dynamics observed during our study and later visits (unpublished data) suggest that host-pathogen dynamics in surveyed lowlands are exhibiting enzootic dynamics, rather than epizootic dynamics (Brem & Lips, 2008;Perez et al., 2014;Woodhams et al., 2008). Our findings also suggest that the distribution of Bd in Costa Rica is wider than historically considered (Puschendorf et al., 2009) and that the population declines during the 1980s and 1990s may not have been restricted to highlands.
Comparable results were found in lowlands of Panama where Bd has been detected in multiple lowland sites (Kilburn et al., 2010;Perez et al., 2014;Woodhams et al., 2008). We suggest that future studies should include replicated sampling across seasons and sites that are outside the optimal environmental conditions for Bd growth, especially since most of these optimal conditions been estimated from lab studies.
The three endangered species sampled (Craugastor taurus, Agalychnis lemur, and Oophaga granulifera) tested positive for Bd (Table 2). The populations of C. taurus and A. lemur that we surveyed also tested positive in past surveys Whitfield et al., 2017). The continuous occurrence of these endangered species and the lack of clinical signs of chytridiomycosis in Bd-infected individuals (Berger et al., 1998;Voyles et al., 2009) Puschendorf et al., 2011). Further studies on these endangered TA B L E 4 Candidacy generalized linear models (GLMs) and linear models (LMs) used to determine the best predictors of prevalence of Batrachochytrium dendrobatidis and infection intensity in amphibian assemblages from four lowland sites and three lowland reproductive habitats in four lowland sites of Costa Rica

Model AIC (GLMs) R 2 (LMs)
Site*habitat ( Note. The most robust models were selected according the highest values for the Akaike information criteria (AIC) for the generalized linear models (GLMs) and the coefficient of regression (R 2 ) for the linear models (LMs). lowland populations can lead to management plans that protect and stabilize these relict populations.

F I G U R E 3 Prevalence and intensity of infection of
The absence of Bd in fourteen surveyed species could be an artifact of the small sample sizes (1-10 individuals,  (Thorpe et al., 2018;Whitfield et al., 2017), to describe host-pathogen population dynamics, or preferably survey multiple species across seasons to obtain more accurate estimates of prevalence and infection intensity for all species within the amphibian community Kinney et al., 2011;Vredenburg et al., 2010).
Our results showed that Bd was widespread across lowlands during the time of study, but Bd prevalence and intensity might exhibit seasonal dynamics. However, to detect a seasonality effect, multi-season studies collecting samples from a variety of amphibian assemblages must be conducted (Kinney et al., 2011;Phillott et al., 2013;Savage et al., 2011). Similar studies conducted in lowlands of  Zumbado-Ulate et al., 2014). Similarly, prevalence of Bd varied from <5% to around 35% in an amphibian assemblage in tropical lowland forest across 1-year period (Whitfield et al., 2012). Therefore, follow-up studies across lowlands in Costa Rica are needed to identify seasonal dynamics of Bd in Costa Rica, which may help design more suitable conservation strategies for lowlands endangered populations.
We did not find Bd in our samples from Rincon de Osa, and a similar study also reported a very low prevalence of Bd in the same study sites and nearby zones across the Osa Peninsula (Goldberg, Hawley, & Waits, 2009). Although our detected prevalence in Rincon de Osa was 0%, our binomial confidence interval (0%-95%) overlaps with the prevalence value presented in this study. Therefore, our result for Rincon de Osa might be an artifact of our low sample size (n = 24) which is not large enough to achieve 95% certainty of detecting 1 positive individual, based on the minimum disease prevalence of ≥5% in infected amphibian assemblages (Skerratt et al., 2008). Climatic conditions at Rincon de Osa might constrain the dispersal and growth of Bd allowing coexistence between susceptible frogs and Bd (i.e., environmental refuge from chytridiomycosis, Puschendorf et al., 2011). However, the extirpation of the Golfito robber frog in this area, where it was abundant before the 1980s and 1990s (Chaves et al., 2014), suggests this may not be the case.
We also found the highest levels of Bd prevalence in the Caribbean sites which coincide with studies conducted in the nearby locations within the same geographic zone (Whitfield et al., , 2013(Whitfield et al., , 2012.
Thus, even within lowland zones, there is large variation in Bd prevalence across zones and sites.
Our statistical models showed no differences among habitats in relation to prevalence and infection intensity. However, there was a trend for higher prevalence of Bd in forest assemblages than in aquatic assemblages (lotic and lentic, Table 3), which differs from other similar studies (Brem & Lips, 2008;Kriger & Hero, 2007a;Lips et al., 2003). Some of the sampled species (e.g., Craugastor fitzingeri,  (Kriger & Hero, 2007a;Lips et al., 2003). Lentic environments are more exposed to sunlight, resulting in temperatures >30°C (Adams et al., 2017), which in laboratory conditions is unsuitable for Bd (Piotrowski et al., 2004). Lentic environments also sustain invertebrates that feed on zoospores reducing the proportion of infected individuals (e.g., Daphnia spp., Searle, Mendelson, Green, & Duffy, 2013).
However, our findings suggest that the role of terrestrial lowland ecosystems in the dispersal of Bd might have been underestimated.
(but see Whitfield et al., 2012;Whitfield et al., 2013). Therefore, multiseason studies contrasting Bd dynamics across habitats are needed to elucidate the role of microhabitats in sustaining Bd.
We found significant evidence that every site of study represents an independent local abiotic environment according to the 19 environmental predictors that we used in our analysis ( Figure 4).
This climatic independence was consistent with the heterogeneous prevalence of Bd, which suggests that every site exhibits a different host-pathogen dynamic in response to local environmental conditions. However, irregularity in elevation gradient across our study sites (Kriger & Hero, 2008), especially in the study site of Siquirres, where elevations varied from 400 to 600 m, could have influenced the differential prevalence we found across lowlands. We recommend controlling for elevational gradients (Kilburn et al., 2010) in follow-up studies. Seasonality and particularly differences in pre-  (Herrera, 1985). Other studies conducted at larger scale have also shown seasonal and latitudinal variation of Bd prevalence and infection Kinney et al., 2011;Kriger et al., 2007;Phillott et al., 2013).
Future studies should evaluate the effect of elevational gradients on the amphibian host-Bd dynamics.
Our results suggest that researchers should expand their sampling across the entire distribution of focal species and communities instead of only focusing on sites of historical declines. An adequate seasonal description of the suitable abiotic environment of pathogens across the host amphibian home range may help identify disease-free sites for effective repatriation or to determine instances were more technical strategies are needed to secure maintenance of declined populations (e.g., antifungal treatments to clear infection, bioaugmentation with commensal bacteria, habitat manipulation, ex-situ conservation) (Garner et al., 2016;Scheele et al., 2014). Furthermore, conducting more seasonal sampling in lowlands will increase the record of presence-absence datasets on Bd and can be used to generate more robust species distribution models (SDMs) from nonopportunistically collected data (Puschendorf et al., 2013). SDMs can help identify hotspots for future outbreaks of Bd and can be used to predict potential locations for amphibian rediscoveries (García-Rodríguez et al., 2012;Puschendorf et al., 2009). Recent validation surveys have led to the discovery of relict peripheral populations that occur in potential environmental refuges from disease Raffel & Fox, 2018;Scheele, Hunter, Skerratt, Brannelly, & Driscoll, 2015), validating increased surveys outside the boundaries of core geographic distributions (Abarca, Chaves, García-Rodríguez, & Vargas, 2010;Chaves et al., 2014;González-Maya et al., 2013;Jiménez & Alvarado, 2017;Nishida, 2006). A comprehensive assessment of a pathogen's distribution, prevalence, and infection intensity can lead to more effective disease-management strategies based on specific locations, habitats, and species.

ACK N OWLED G M ENTS
We thank the InterAmerican Network of Academies of Science (IANAS)

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest exists.

AUTH O R CO NTR I B UTI O N S
HZ-U, AG-R, CS, and VV developed the ideas and designed methodology; HZ-U and AG-R conducted data collection and data analysis; HZ-U, AG-R, CS, and VV led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

DATA ACCE SS I B I LIT Y
The data associated with this publication (Data files title: Bd_ lowlands_Costa_Rica) are deposited at Dryad data repository.