**1. 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 managemen<sup>t</sup> efforts [11,12].

Costa Rica, with an area of only 51,100 km2, is home to a grea<sup>t</sup> diversity of amphibians [13]. More than 200 of the approximately 8000 described amphibian species are present in Costa Rica [14], and new species continue to be described. The vast species richness confined to a relatively small

area is due to complex biogeographic events and climatic conditions throughout the country, and a long history of work has been done by in-country taxonomic specialists [13,15]. Costa Rica is also an example of a country where numerous amphibian population declines have occurred in response to multiple environmental threats [16], highlighted by the enigmatic disappearance of the golden toad (*Incilius periglenes*) [17]. However, several species that catastrophically declined in the last thirty years, such as the harlequin frog, the Golfito Robber frog, and the Holdridge's toad, have been recently rediscovered in viable populations [18–20]. These findings sugges<sup>t</sup> that highly susceptible species can recover from or at least persist when faced with deadly threats [21]. Thus, Costa Rica is an excellent location to study not only how amphibian communities are a ffected by environmental threats but also their resistance and resilience from declines [22].

One widespread cause of amphibian declines is the introduction of the pathogen *Batrachochytrium dendrobatidis* (*Bd*) [23]. This fungus causes chytridiomycosis [24], a potentially deadly skin disease that has contributed to the decline of at least 500 amphibian species globally [10]. In Central America, amphibian declines peaked during the 1980s and 1990s and have been linked to the introduction of *Bd*, which caused deadly outbreaks of chytridiomycosis (i.e., epizootics) [16,25,26]. It has been suggested that *Bd*-driven epizootic declines mostly a ffected species in highland lotic environments because moisture and temperature in these sites matches the optimal conditions for *Bd* growth in the lab [27,28]. However, it is also well known that some amphibian species su ffered unexpected and unexplained declines in lowland environments (<700 m above sea level) during the 1980s and 1990s, likely due to chytridiomycosis [16,19,29–31]. After the declines, the evolution of resistance and tolerance mechanisms in amphibian communities [32], and/or the evolution of less-pathogenic strains of *Bd* [33], might have allowed susceptible amphibians to persist with endemic *Bd* infection (i.e., enzootics) [22,34–36]. However, susceptible species are still at a high risk of extinction under endemic infection if conditions shift in favor of the pathogen. For example, the introduction of an invasive species that is also a competent reservoir might amplify infection in the environment to epizootic levels [37–39]. Thus, examining the life history traits and conditions that may favor outbreaks of *Bd* is the key to understanding the underlying mechanisms behind why some infected species declined more severely than others and which species are most vulnerable to future outbreaks [40,41].

In this study, we present an updated list of all the amphibian species of Costa Rica, quantifying species diversity in each herpetological province and describing their conservation status. We also identified the e ffect of geography (herpetological province and altitudinal belt) and life-history traits associated with foraging and reproduction on current infection with *Bd*. For this, we sampled for *Bd* at nine tropical localities across Costa Rica from 2015–2018. In addition, we built a robust dataset by adding records from studies that detected *Bd* across Costa Rica from 2005–2018 in multi-species amphibian assemblages. We hypothesized that *Bd* is widespread across herpetological provinces and altitudinal belts in Costa Rica and would exhibit an infection intensity below theoretical thresholds associated with mass mortalities [42]. To compare across life-history traits, we developed an index that combines foraging habitat, reproductive habitat, and type of development. We hypothesized that *Bd* infection would vary across habitats, with the highest prevalence and infection intensity found in species with the greatest use of cool and humid environments [43,44]. The knowledge from this work will aid policy-makers in identifying the most threatened regions and taxonomic clades to develop better conservation strategies in Costa Rica [22,45,46].

## **2. Materials and Methods**

## *2.1. Species Assessment*

We updated the last o fficial list of amphibian species in Costa Rica published in 2011 [47] by consulting the Herpetological Database ("Herp Database") of the Museo de Zoología at Universidad de Costa Rica (http://museo.biologia.ucr.ac.cr/) and taxonomists' lists [48,49]. In addition, we georeferenced 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).

For every species, we showed their status in Costa Rica [50] according to the International Union of Conservation of Nature (IUCN) [51] as follows: NA = "not applicable," DD = "data deficient," LC = "least concerned," NT = "near threatened," VU = "vulnerable," EN = "endangered," CR = "critically endangered," and EX = "extinct in the wild" (for additional details see http://www.iucnredlist.org/). We also included environmental vulnerability scores (EVS) [52], a regional vulnerability index that classifies amphibians and reptiles into four levels of risk: "no immediate risk" (EVS < 3), "low vulnerability" (EVS of 3–9), "medium vulnerability" (EVS of 10–13), and "high vulnerability" (EVS of 14–17). A high EVS indicates species that are restricted in distribution, occur in a single life zone, and have a highly derived reproductive mode. The EVS for Costa Rican amphibians reported here were extracted from Sasa et al. [49]. Finally, we compiled a list of all the species that have been screened for *Bd* and the methods used for detection: histology or polymerase chain reaction (PCR).

## *2.2. 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.

## *2.3. 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–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–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).

**Figure 1.** Map of 20 survey sites across Costa Rica. Sites are color-coded by herpetological province.

**Table 1.** Summary of studies where *Batrachochytrium dendrobatidis* (*Bd*) was detected in multi-species amphibian assemblages using conventional PCR and quantitative PCR (qPCR) in Costa Rica between 2005–2018. The table shows surveyed localities, herpetological province, sampling period, percentage of infection, and Holdridge's altitudinal belt. Symbology: CL—Caribbean Lowlands, MSCC—Montane Slopes and Cordillera Central, PN—Pacific Northwest, and PS—Pacific Southwest.

