Endemic Infection of Batrachochytrium dendrobatidis in Costa Rica: Implications for Amphibian Conservation at Regional and Species Level

: Batrachochytrium dendrobatidis ( Bd ) has been associated with the severe declines and extinctions of amphibians in Costa Rica that primarily occurred during the 1980s and 1990s. However, the current impact of Bd infection on amphibian species in Costa Rica is unknown. We aimed to update the list of amphibian species in Costa Rica and evaluate the prevalence and infection intensity of Bd infection across the country to aid in the development of e ﬀ ective conservation strategies for amphibians. We reviewed taxonomic lists and included new species descriptions and records for a total of 215 amphibian species in Costa Rica. We also sampled for Bd at nine localities from 2015–2018 and combined these data with additional Bd occurrence data from multiple studies conducted in amphibian communities across Costa Rica from 2005–2018. With this combined dataset, we found that Bd was common (overall infection rate of 23%) across regions and elevations, but infection intensity was below theoretical thresholds associated with mortality. Bd was also more prevalent in Caribbean lowlands and in terrestrial amphibians with an aquatic larval stage; meanwhile, infection load was the highest in direct-developing species (forest and stream-dwellers). Our ﬁndings can be used to prioritize regions and taxonomic groups for conservation strategies.


Introduction
Anthropogenic threats including habitat destruction, pollution, climate change, introduction of invasive species, and pathogens are causing a rapid and severe decline in global biodiversity [1]. Scientific consensus states that we are in the midst of a sixth mass extinction event [2,3]. Within vertebrates, amphibians are the most endangered taxonomic class with approximately 41% of described species classified as "globally threatened" [4,5]. The majority of the amphibian declines have occurred in the tropics of Australia, Central America, and South America [6,7], and have been observed even in seemingly pristine and protected environments [8,9]. However, information is still lacking regarding which species are suffering the greatest declines and which abiotic and biotic factors are contributing the most [10]. Identifying threatened species and factors contributing to global amphibian declines is vital for effective conservation and management efforts [11,12].
Costa Rica, with an area of only 51,100 km 2 , is home to a great diversity of amphibians [13]. More than 200 of the approximately 8000 described amphibian species are present in Costa Rica [14], the distribution of all amphibian species within the five Costa Rican herpetological provinces (see Section 2.3.1). For this, we extracted all collection points for each species from the "Herp Database" (Datum WGS1984) and mapped them using a shapefile of the Costa Rican herpetological provinces and QGIS software 3.8.1 (QGIS Development Team, http://qgis.osgeo.org).

Field Dataset
To add to existing datasets of amphibian distribution and Bd infection, we surveyed nine amphibian assemblages across Costa Rica in both versants (Caribbean and Pacific) and at elevations ranging from sea level to 1385 m ( Figure S1). All surveys were conducted during the months of June and July between 2016-2018, except in the locality of Alto Lari, which was sampled in March 2015. At each site, we conducted visual and acoustic encounter surveys searching for amphibians in streams, ponds/puddles, and forest (leaflitter and canopy), and then caught individuals to screen them for Bd (see below). In total, we screened for Bd from 267 amphibians from 33 species (see Tables S1 and S2, Figure S2). Four of those species were classified in threatened categories: Oophaga granulifera (VU), Ptychohyla legleri (EN), Craugastor ranoides (CR), and C. taurus (CR).
All observed amphibians were collected using nitrile gloves and temporally placed 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 [24,53]. We swabbed (using MW-113 swabs) each individual's skin to detect Bd as follows: five strokes on one hand, five strokes on the ventral patch, five strokes on one foot, and five strokes along the inner thigh. The swabs were placed into 1.5 mL screw-cap tubes and stored dry at −20 • C until fungal DNA extraction. Once swabbed, all animals were released back to the site of their collection. During this study we followed field protocols [54,55] approved by the National System of Conservation Areas of Costa Rica (SINAC), the Comisión Nacional para la Gestión de la Biodiversidad (CONAGEBIO), and animal care protocols from the Purdue Animal Care and Use Committee (PACUC 1604001392), ensuring that all animals were being cared for in accordance with standard protocols and treated in an ethical manner (research permits 001-2012-SINAC, R-019-2016-OT-CONAGEBIO, R-023-2016-OT-CONAGEBIO, R-057-2016-OT-CONAGEBIO, R-060-2016-OT-CONAGEBIO).
We conducted diagnostic quantitative polymerase chain reaction (qPCR) on each swab to quantify Bd infection load following standard protocols [56], with the following modifications: (1) the fungal DNA was extracted using 60 µL of PrepMan Ultra, and (2) an internal positive control (IPC) was used to detect inhibitors [57]. Fungal DNA was diluted 1:10 in 0.25× TE buffer and run in singlicate [58] using a Step One Plus (Applied Biosystems, Woburn, MA, USA). Negative controls (DNase/RNase-free water) were run in triplicate on every 96-well qPCR plate. We classified samples as positive when both dyes (Bd probe and IPC) amplified in each well. Samples absent of IPC amplification were considered inhibited. In order to eliminate inhibitors, we diluted 5 µL of a new dilution in 0.25×TE buffer in a proportion of 1:100. Ten samples were classified as inhibited and then determined to be negative after dilution. Quantification curves for genomic equivalents were constructed using 1000, 100, 10, and 1 zoospore quantification standards derived from a gBlock ® Gene fragment (Integrated DNA Technologies, Coralville, IA, USA). In order to calculate the zoospore genomic equivalents in the original sample, we multiplied the qPCR score by the dilution factor of 120 (dilution factor = 60 × 20 × 1/10). We estimated prevalence with 95% binomial confidence intervals (CIs) by locality.

Combined Dataset
We generated a dataset from multiple studies that screened for Bd in multiple amphibian assemblages in Costa Rica after the year 2000 using conventional PCR and qPCR methods [59][60][61][62][63][64][65] ( Figure 1, Table 1) (including the 267 individuals from the 33 species we tested in the "field dataset" (see Section 2.2 and supporting data). In total, this "combined dataset" consisted of 1750 individual records from 79 species and 20 localities at elevations ranging from sea level to 2000 m. We identified the year 2000 as the starting of post-decline because most epizootic outbreaks of Bd occurred during the 1980s and early 1990s [26,30,66]. We also assumed that Bd expanded its range across Costa Rica by 2000 due to the rapid rate of spreading that this pathogen exhibits in tropical locations [67,68]. Although Bd was detected in 405 swabs in this dataset, quantification through qPCR was conducted only in 351 Bd-positive swabs (from the "field dataset" and three of the seven reviewed studies [61][62][63]). We did not consider studies that used histology as a method of detection because most of these studies evaluated samples that were taken before 2000. We also excluded records of individuals that were identified only at the genus level, e.g., Craugastor spp. [64] and Agalychnis spp. [61], and cases where only one species was screened for Bd (e.g., Atelopus varius in the locality of Uvita [65]). Finally, we classified all sampled amphibians according to herpetological province, altitudinal belt, and life history traits (foraging habitat, reproductive habitat, and type of development).  (Integrated DNA Technologies, Coralville, IA, USA). In order to calculate the zoospore genomic equivalents in the original sample, we multiplied the qPCR score by the dilution factor of 120 (dilution factor = 60 × 20 × 1/10). We estimated prevalence with 95% binomial confidence intervals (CIs) by locality.

Combined Dataset
We generated a dataset from multiple studies that screened for Bd in multiple amphibian assemblages in Costa Rica after the year 2000 using conventional PCR and qPCR methods [59][60][61][62][63][64][65] ( Figure 1, Table 1) (including the 267 individuals from the 33 species we tested in the "field dataset" (see Section 2.2 and supporting data). In total, this "combined dataset" consisted of 1750 individual records from 79 species and 20 localities at elevations ranging from sea level to 2000 m. We identified the year 2000 as the starting of post-decline because most epizootic outbreaks of Bd occurred during the 1980s and early 1990s [26,30,66]. We also assumed that Bd expanded its range across Costa Rica by 2000 due to the rapid rate of spreading that this pathogen exhibits in tropical locations [67,68]. Although Bd was detected in 405 swabs in this dataset, quantification through qPCR was conducted only in 351 Bd-positive swabs (from the "field dataset" and three of the seven reviewed studies [61][62][63]). We did not consider studies that used histology as a method of detection because most of these studies evaluated samples that were taken before 2000. We also excluded records of individuals that were identified only at the genus level, e.g., Craugastor spp. [64] and Agalychnis spp. [61], and cases where only one species was screened for Bd (e.g., Atelopus varius in the locality of Uvita [65]). Finally, we classified all sampled amphibians according to herpetological province, altitudinal belt, and life history traits (foraging habitat, reproductive habitat, and type of development).   We classified all surveyed assemblages within the five herpetological provinces proposed by Savage [13] and modified by Sasa and colleagues [49].
Caribbean Lowlands: This faunal area represents 30% of Costa Rica and includes the lowlands of the Caribbean versant and the northernmost region of the country, predominantly consisting of lowland wet forest. Sampling for Bd through PCR has been conducted in the localities of La Virgen de Sarapiquí, Puerto Viejo de Sarapiquí, Tres Equis de Turrialba, Guayacán de Siquirres, Kekoldi, and the remote Alto Lari. Alto Lari was surveyed as part of a recent expedition following an enigmatic path that connects the Caribbean Lowlands with the highlands of Cordillera de Talamanca and is known as "the Gabb's route" [69].
Pacific Northwest: This herpetological province includes the lowlands of the Pacific Northwest and extends into the western side of the Central Valley, in the Meseta Central Occidental (Central Valley) up to the base of Cerros de Ochomogo. The Pacific Northwest consists of predominantly Lowland Dry Forest vegetation and constitutes 24% of Costa Rica's area. This province contains a distinctive dry season that lasts five to six months. Within the Pacific Northwest, sampling for Bd has been conducted in the tropical dry forest at Guanacaste National Park (Santa Rosa and Santa Elena Peninsula stations).
Pacific Southwest: Encompassing the lowlands of the Pacific central and south, the herpetological provinces consist primarily of lowland wet forest and lowland moist forest and accounts for 15% of the country's area. This herpetological province is biogeographically related to the Caribbean Lowlands and species have more recently been differentiated between these herpetological provinces due to isolation caused by the uplifting of the Cordillera de Talamanca. Within the Pacific Southwest, sampling for Bd through PCR has been conducted in the localities of Punta Banco-Burica, Golfito, Rincón de Osa, Piro, Corcovado, and Uvita.

Montane Slopes and Cordillera Central:
This area represents 23% of Costa Rica and occurs along all of Costa Rica's mountain ranges from 500-2100 m elevation in Cordillera de Guanacaste, 500-3400 in Cordillera Central, and 500-1600 in Cordillera de Talamanca (Lower Talamanca). The Montane Slopes and Cordillera Central province includes regions that receive the highest annual precipitation in the country. Sampling for Bd through PCR has been conducted in the localities of Monteverde, San Vito de Coto Brus, Las Tablas, Tinamastes de Pérez Zeledón, San Rafael de Heredia, and Santo Domingo de Heredia.
Cordillera de Talamanca: Found at elevations above 1600 m (Upper Talamanca). This is the smallest faunal area (8% of Costa Rica) and consists primarily of montane rainforest and subalpine pluvial paramo. This faunal area is the least explored herpetological province of Costa Rica and there is no published data for Bd detection through PCR or qPCR in this faunal province.

Foraging-Reproduction Habitat Index
To compare Bd infection dynamics across taxonomic groups, we developed a foraging-reproduction habitat index (FRHI). The FRHI was created to classify species with a system of three letters that represented life history traits associated with foraging and reproduction (Table 2). First, we classified species according to their development into indirect-(I) or direct-developing amphibians (D). Second, we classified species according to their foraging habitat into terrestrial (T), arboreal (A), pond/puddle-dwellers (P), stream-dwellers (R), and phytotelma (F). Finally, species were classified according to their reproductive habitat into terrestrial (T), arboreal (A), pond/puddle-breeders (P), stream-breeders (R), and phytotelma (F). Table 2. Categories and taxonomic examples of the foraging-reproduction habitat index (FRHI) that we developed for this study to analyze current prevalence of Batrachochytrium dendrobatidis across taxonomic groups. Symbology: First letter represents development: (I) indirect-or (D) direct-developing amphibians. Second letter represents foraging habitat: terrestrial (T), arboreal (A) pond/puddle-dwellers (P), stream-breeders (R), and phytotelma (F). Third letter represents reproductive habitat: terrestrial (T), arboreal (A), pond/puddle-dwellers (P), stream-breeders (R), and phytotelma (F).

Statistical Analysis
We reduced our "combined dataset" to 1741 samples from 74 species and 20 localities because there was insufficient information to accurately classify nine records of the species Diasporus tigrillo, D. vocator, Hyloscirtus palmeri, Triprion spinosus, and Cruziohyla calcarifer in the FRHI. For our analyses, we pooled all species together instead of using species as a predictor because the samples sizes per species were highly variable (from 1-177), which could produce significant models that may be an artifact of opportunistic sampling instead of a real pattern. Instead, we used the FRHI, which is highly correlated with taxonomic group. We were unable to include time of sampling as a predictor in our analyses because these data were missing in several of the amphibian assemblages sampled. All our analyses were performed with the R package "stats" [71].
We analyzed Bd prevalence with fixed-effects generalized linear models (GLMs) using infected status as a binomial response variable (uninfected or infected) and herpetological province, altitudinal belt, and the FRHI as predictors. Ranking of the candidate GLMs followed the Akaike's information criterion (AIC) where the model with the lowest AIC was considered the most robust [72]. To analyze Bd infection intensity (estimated as the number of genomic equivalents), we analyzed the 351 Bd-positive swabs where Bd infection intensity was quantified through qPCR (see Section 2.3). We used linear models (LMs) to compare Bd infection intensity (response variable) across herpetological provinces, altitudinal belts, and FRHI (predictor variables). We log-transformed the Bd infection intensity to reduce skewness. Statistical significance of models was tested with ANOVA. For both, GLMs and LMs, we conducted post hoc, pairwise comparisons (Tukey's honestly significant difference; HSD-test) to confirm where the differences occurred between significant predictors. We were unable to run mix-effects models or fixed-effects interaction models because some combinations of predictors presented missing or low values, causing convergence difficulties.
Overall, 200 species have been classified into IUCN categories at the country level and 15 species need future assessment. A total of 155 species do not fulfill the criteria to be considered in the threatened categories, including 136 species listed as LC, 18 as DD, and one as NA (a category for taxa that occur in the region but have been excluded from the regional Red List for a specific reason). Within threatened categories, two species are listed as EX, 24 as CR, ten as EN, seven as VU, and two as NT. Regionally, lowlands exhibited the lowest percentage (0-2%) of DD species (Figure 2b). Similarly, about 75-80% of species occurring in lowlands are listed as LC (Figure 2c). In highlands, 6-10% of species are categorized as DD (Figure 2b) and 26% of species in Cordillera de Talamanca are classified within threatened categories (Figure 2d). According to EVS, a total of 81 species were classified as "no immediate risk," four species at "low vulnerability," 50 species at "medium vulnerability," and 48 species at "high vulnerability" (Table S3).
In our review, we found a total of 105 amphibian species (49%) that have been screened for Bd in Costa Rica (103 anurans and only 2 species of salamanders) (Table S4). In the field, the most common method used to detect Bd was qPCR, especially after 2005. Conventional PCR was used only in one study in the Caribbean Lowlands [64]. Histology and qPCR have also been used in retrospective studies on preserved specimens from declined and extinct species.

Endemic Dynamics
Overall, Bd prevalence in Costa Rica was estimated to be 0.23 (60% of species tested positive for Bd) ( Table S5). The most robust GLM found both herpetological province and the FRHI as significant predictors of Bd prevalence (AIC = 1700, p < 0.001, Table S6). Among herpetological provinces ( Figure  3), the highest percentage of infected individuals was found in the Caribbean Lowlands (34%) and the lowest in the Pacific Northwest (4%). The Mountain Slopes, Cordillera Central, and Pacific Southwest had a similar percentage of infected individuals (≈23%). Furthermore, Bd was proportionally more prevalent in amphibians with terrestrial foraging and larval stage in phytotelma (ITF), pond-breeding treefrogs (IAP), and direct-developing species that breed in the forest (leaf-litter frogs DTT, rain frogs, DAT) (Figure 4a).
The species Craugastor taurus (the Golfito robber frog) was the species that had the highest average infection load (average Bd load of 11632.4 versus 571.6 genomic equivalents or 2.51 versus 1.18 after log transformation) ( Table 4). We found an effect of the FRHI on infection load (F8342 = 7.91, p < 0.01, Table S6). Direct-developing frogs with terrestrial reproduction (robber frogs and leaf-litter frogs; DTR and DTT respectively) had the highest Bd loads (Figure 4b, Table 4).

Endemic Dynamics
Overall, Bd prevalence in Costa Rica was estimated to be 0.23 (60% of species tested positive for Bd) ( Table S5). The most robust GLM found both herpetological province and the FRHI as significant predictors of Bd prevalence (AIC = 1700, p < 0.001, Table S6). Among herpetological provinces (Figure 3), the highest percentage of infected individuals was found in the Caribbean Lowlands (34%) and the lowest in the Pacific Northwest (4%). The Mountain Slopes, Cordillera Central, and Pacific Southwest had a similar percentage of infected individuals (≈23%). Furthermore, Bd was proportionally more prevalent in amphibians with terrestrial foraging and larval stage in phytotelma (ITF), pond-breeding treefrogs (IAP), and direct-developing species that breed in the forest (leaf-litter frogs DTT, rain frogs, DAT) (Figure 4a).
The species Craugastor taurus (the Golfito robber frog) was the species that had the highest average infection load (average Bd load of 11632.4 versus 571.6 genomic equivalents or 2.51 versus 1.18 after log transformation) ( Table 4). We found an effect of the FRHI on infection load (F 8342 = 7.91, p < 0.01, Table S6). Direct-developing frogs with terrestrial reproduction (robber frogs and leaf-litter frogs; DTR and DTT respectively) had the highest Bd loads (Figure 4b, Table 4).

Species Assessment
We presented the first updated list of Costa Rican amphibians since 2011 [47]. Compared to the last list, we added 10 anurans, 9 salamanders, and 1 caecilian (Table 3) for a total of 215 species (Table  S3). As is common throughout the world, anurans exhibited the highest amphibian species richness in Costa Rica, with 72% of listed species. However, the richness of salamander species is also high (25%). In Costa Rica, the diversity and endemism of amphibians (especially salamanders) increase with elevation and complex mountain topography [13]. Proportionally, in terms of number of species per unit area (km 2 ), the richest herpetological province was the Cordillera de Talamanca. In this herpetological province, the number of species continued to increase and most of the newly described species in our report came from this remote and almost inaccessible province [69]. The Montane Slopes and Cordillera Central presented the highest number of species (158 species). Within this province, numerous mountain ranges provide multiple microhabitats for niche differentiation and further speciation [13,49]. In lowlands, the highest number of amphibian species occurred in the Caribbean Lowlands with 101 species. However, the Pacific Southwest presented more species per unit area. The Pacific Northwest only had 66 species, which was also the lowest number of species per unit area. This pattern may be attributed to the warm and dry conditions that occur in most part of this herpetological province [13,49].
According IUCN, the species Craugastor escoces and Incilius periglenes are classified as EX in Costa Rica. However, C. escoces was recently rediscovered [89]. Similarly, several species that remained undetected after the 1980s and 1990s such as Incilius holdridgei [20], Craugastor taurus [19], and Atelopus varius [18] have been rediscovered in peripheral populations during the last few years. However, the number of extinct species could be higher because multiple threatened species still remain undetected in the field (e.g., Craugastor andi, Incilius fastidiosus, Atelopus senex). We recommend expedition surveys to find populations of declined and data-deficient species [90] and captive-breeding for species where ex situ reproduction has been successful (e.g., harlequin frogs.) [91]. Although we acknowledge that there are limited funds available for these types of conservation efforts, knowledge from these sites is essential to be able to identify conditions that favor the persistence of threatened species and identify species that should be targeted for future conservation efforts [92].
Lowlands of Costa Rica exhibited the lowest proportion of DD species (0-2%; Figure 2b) and the highest proportion of LC species (75-80%; Figure 2c). On the other hand, highlands exhibited the highest percentage of DD species (6-10%; Figure 2b). Similarly, Cordillera de Talamanca had the highest percentage of threatened species (26%; Figure 2d). Based on these findings, we strongly recommend increasing the sampling effort in the montane and subalpine altitudinal belts (>2800 m) that exclusively occur in Cordillera de Talamanca and Montane Slopes and Cordillera Central. These herpetological provinces present the highest rate of endemism (especially for salamanders) and contain several of the recently described species [85,86]. Conducting expeditions and long-term studies in highlands will aid in monitoring threatened species and reducing information gaps, allowing for more accurate assessments of amphibian species.
To better evaluate the vulnerability of amphibian species, we utilized EVS (Table S3). This index relies on ecological information for categorizing threat levels, which makes application easy for most species in a specific region [52]. Unlike the IUCN Red List of Threatened Species, previous evaluations of threats are not considered by this index [49]. For that reason, species that are classified as LC by IUCN can be classified as highly vulnerable in this index (e.g., Duellmanohyla rufioculis). We categorized 48 species (24 salamanders, 21 anurans, and 3 caecilians) in "high vulnerability" (e.g., Atelopus senex, Craugastor andi, Bolitoglossa pesrubra, Nototriton guanacaste, and Oscaecilia osae). These species exhibited the highest EVS values because their habitats are restricted and because they exhibit complex reproductive modes. Quantifying environmental threats and combining information from both indexes will help policy-makers to prioritize conservation actions for threatened species.
In Central America, habitat destruction is the most important threat impacting amphibian populations [49,52]. Although approximately 30% of Costa Rica remains forested and protected, rapid urbanization, extensive agriculture, excessive pesticide use, illegal traffic, and inappropriate waste management negatively affect numerous amphibian populations. However, even in seemingly pristine locations, amphibian declines have occurred [93]. Additionally, climate change has been associated with the decline of several amphibian species in Costa Rica by affecting their reproduction and likely increasing susceptibility to pathogens [28]. Although it has not been found in Central America, we recommend screening for the recently emerged fungus Batrachochytrium salamandrivorans [94], which causes chytridiomycosis in salamanders. Conveniently, swabbing methods and qPCR allow for accurate detection of both fungi species in the same assay, which may facilitate rapid population assessments. For a fully detailed review of the environmental threats for amphibian communities in Costa Rica, we recommend the work of Sasa et al. [49].

Post-Decline Dynamics
In this study, we found strong evidence that Bd is widespread in Costa Rica. Our results also suggest that post-decline Bd exhibited enzootic dynamics, characterized by high prevalence of infection across regions and pathogenic loads below thresholds associated with mass mortalities [39,42]. We found Bd in all the herpetological provinces and altitudinal belts surveyed (Figures 3 and 4, Tables S5 and S6), for a total infection rate of 23%. We also found that Bd was more prevalent in terrestrial amphibians with an aquatic larval stage and direct-developing frogs exhibited the highest pathogenic loads.
We found the lowest infection rate (<5%) in the Pacific Northwest and the highest (33%) in the Caribbean Lowlands ( Figure 3). However, these values may be an effect of sampling periods. A study conducted at La Selva Biological Station in the Caribbean Lowlands [95] reported infection rates varying from <5% during the warmest months (May to early November) to 35% in the coolest months (mid-November to January) in three common amphibian species. In addition, the gradual population declines observed at La Selva over several decades [31] and opportunistic observations of small-scale mortality events during cold periods [96] suggest that Bd may be causing mortality in amphibians long after its initial invasion. Similar mortality events in response to seasonality could be occurring in amphibian communities in "refuges from decline" in the Pacific Northwest [62,97] and Pacific Southwest [63] of Costa Rica; however, they have not yet been documented. There, seasonal changes in precipitation and temperature caused Bd prevalence to vary from >5% in the peak of the dry season (March and April) to 80% in the coldest months (November-December). Therefore, we recommend follow-up studies at these sites to identify whether seasonal disease dynamics are causing mortality events in regions that have been considered unsuitable for Bd [98].
Bd was found across all altitudinal belts across Costa Rica. Similar results of high Bd prevalence across all elevations has been found in Panama [99][100][101]. These findings suggest that current environmental conditions are suitable for Bd at most elevations in Central America [63]. It is also plausible that Bd-driven declines during the 1980s and 1990s were not exclusively restricted to highlands [30] but were relatively undetected at lower elevations. Another hypothesis is that species with high susceptibility historically occupied high elevations sites, but severely declined or went extinct after Bd was introduced, leaving only species with mid-to-low susceptibility across elevations [102]. On the other hand, the absence of samples from montane and subalpine belts and uncontrolled variables (e.g., changes in species composition, climatic disturbances) could have reduced the statistical power to determine changes in Bd prevalence across altitudinal belts. We suggest that future studies increase sampling at high elevations (>2700 m) to better understand the local spatial dynamics of Bd across elevations.
Our results showed that infection with Bd was common in amphibians across all life-history traits evaluated in the FRHI (Table 2). However, Bd was significantly more prevalent in terrestrial amphibians with a larval stage (Figure 4a), especially those that complete metamorphosis in phytotelma (ITF). All the species within the ITF category belonged to the family Dendrobatidae, which have previously been shown to easily acquire Bd infection [62,95]. The high susceptibly of the Dendrobatidae family is likely due to their preferred habitat (e.g., water-filled bromeliads for many species), as it offers suitable conditions for Bd infection [103]. In addition, dendrobatid adults forage in the tropical forest floor and stream-associated low vegetation, which are environments that can sustain Bd [104]. Regarding infection intensity (Figure 4b), the FRHI showed similar results to studies that have used the aquatic index [43]. We found that direct-developing species with terrestrial reproduction had significantly higher infection load than other species with different life-history traits (Figure 4b, Table 4). This life-history trait is exhibited by leaf-litter frogs and all the species within the Craugastor punctariolus clade (robber frogs), which is one of the clades most affected by chytridiomycosis in Central America [29,105]. Robber frogs spend a majority of their life cycle along fast-flowing streams [106], an aquatic environment that seems highly suitable for Bd in Central America [107]. In addition, these frogs appear to be highly susceptible to Bd-driven mass mortalities outside warm and dry ecosystems [19,30,108]. Our results suggest the FRHI is particularly useful for identifying taxonomic units that are more susceptible to Bd [92].

Conclusions
Our results demonstrated that the number of identified amphibian species in Costa Rica is still growing, and there may be potential future additions (e.g., Bolitoglossa anthracina [109] and B. indio [110]). A continuous assessment of species and regions is needed to identify continuing threats to amphibian biodiversity. We found that Bd was widespread across species, herpetological provinces, and altitudinal belts in samples collected since 2000. Conducting more studies in remote regions, such as Cordillera de Talamanca, may help to better describe spatial dynamics of both amphibian hosts and Bd. In addition, future studies should test whether seasonal disease dynamics are causing mortality events in regions that are considered unsuitable for Bd. Under potential scenarios of climate change, environmental conditions may shift to ideal ranges for Bd infection [28] and seasonal regions that sustain critically endangered species (e.g., tropical dry forest) may experience future outbreaks of chytridiomycosis [111]. We also recommend continuous surveillance of invasive species, which might amplify Bd in the environment, causing future epizootics [37]. This vital information will aid in the development of more effective conservation strategies for amphibians across a broader range of habitats [46,[112][113][114].
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/11/8/129/s1. Figure S1: Map of Costa Rica showing elevational gradient and nine study sites surveyed for Batrachochytrium dendrobatidis, Figure S2: Prevalence and intensity of infection of Batrachochytrium dendrobatidis in amphibian assemblages from the nine surveyed sites in Costa Rica, Table S1: List of species where Batrachochytrium dendrobatidis was surveyed in Costa Rica for our "field dataset", Table S2: Prevalence and infection intensity of Batrachochytrium dendrobatidis at nine sites in Costa Rica, Table S3: Distribution of amphibians in Costa Rica according to herpetological province and elevation, Table S4: List of species that have been screened for Batrachochytrium dendrobatidis in Costa Rica, Table S5: List of species where Batrachochytrium dendrobatidis was surveyed in Costa Rica in our combined dataset, Table S6: Candidacy generalized linear models (GLMs) and linear models (LMs) used to determine the best predictors of prevalence and infection intensity of Batrachochytrium dendrobatidis in amphibian assemblages from Costa Rica.