Epizootic to enzootic transition of a fungal disease in tropical Andean frogs: Are surviving species still susceptible?

The fungal pathogen Batrachochytrium dendrobatidis (Bd), which causes the disease chytridiomycosis, has been linked to catastrophic amphibian declines throughout the world. Amphibians differ in their vulnerability to chytridiomycosis; some species experience epizootics followed by collapse while others exhibit stable host/pathogen dynamics where most amphibian hosts survive in the presence of Bd (e.g., in the enzootic state). Little is known about the factors that drive the transition between the two disease states within a community, or whether populations of species that survived the initial epizootic are stable, yet this information is essential for conservation and theory. Our study focuses on a diverse Peruvian amphibian community that experienced a Bd-caused collapse. We explore host/Bd dynamics of eight surviving species a decade after the mass extinction by using population level disease metrics and Bd-susceptibility trials. We found that three of the eight species continue to be susceptible to Bd, and that their populations are declining. Only one species is growing in numbers and it was non-susceptible in our trials. Our study suggests that some species remain vulnerable to Bd and exhibit ongoing population declines in enzootic systems where Bd-host dynamics are assumed to be stable.


Introduction
The eastern slopes of the Andes have one of the highest diversity of amphibians in the world [1]. The recent epizootic (i.e., an epidemic in animal populations) of a pathogenic fungus, Batrachochytrium dendrobatidis (Bd; [2]) which causes the fatal disease chytridiomycosis, has caused the widespread collapse of montane Neotropical amphibian communities [3][4][5][6]. Although the pathways of disease transmission are still largely unknown, several studies have proposed that the disease may spread in epizootic waves throughout even the most remote PLOS  to the outbreak. We collected Bd infection prevalence and intensity estimates, important factors used to describe disease host/Bd dynamics [12,14,25]. As a first step towards understanding the impact of Bd on this community, we conducted Bd susceptibility trials in situ on the eight focal species. We compare our experimental results to population trends and infection dynamics in the wild, and explored whether there was any association between species susceptibility to Bd and change in population abundance before and after the epizootic.

Ethics statement
This work has been approved by the Animal Care and Use Committees of San Francisco State University (Protocol #A12-07) and Southern Illinois University (Protocol #13-027), and permits to carry on this research have been issued by the Peruvian Ministry of Agriculture. The Asociación para la Conservación de la Cuenca Amazónica authorized work at its Wayqecha Biological Station.

Study site
Our study took place in the Kosñipata Valley adjacent to Manu National Park near Cusco, southeastern Peru [26]. Manu National Park includes 17,163 km 2 of largely undisturbed Amazonian lowlands and Andean mountains from 300 m to 4020 m in elevation. In this study, we sampled habitats continuously from 1200 m to 3800 m, and captured frogs in montane forests, scrub, and high-Andean grasslands (Fig 1). Before the chytridiomycosis epizootic, 55 species occurred in this elevational range; after the outbreak, 19 species disappeared [5,24].

In situ Bd susceptibility trials
We conducted susceptibility trials in rustic conditions at the remote Wayqecha Biological Field station located in our study area [5,27]. We collected 122 frogs from the eight focal species for the susceptibility trials. We hand-captured frogs and placed them individually in plastic bags until they were placed in experimental chambers. All animals were swabbed on the skin for Bd, measured (SVL), weighed (g), and sexed. Because the study area contains Bd, we immersed all study animals in a 1% itraconazole solution for 5 minutes a day for seven consecutive days before susceptibility trials began [28]. After the anti-fungal procedure, each animal was randomly assigned into either the treatment group (Bd exposed) or the control group (no Bd exposure). We used 78 Bd-infected Telmatobius marmoratus, a highly Bd-infected species native to the study area, to expose our study species to Bd, because fungal cultures could not be maintained at the remote field station. Infected T. marmoratus were purchased from a live frog stand at the San Pedro market in Cusco, Peru. We purchased 16-20 live T. marmoratus, on four different dates (22 June, 3, 18 and 23 July 2012). In previous studies we showed that animals from this market stand were all highly infected with Bd [29,30]. Since we did our Bd susceptibility trials in situ, we analyzed Bd swabs after the end of the entire experiment. Animals in the susceptibility trials were either exposed to T. marmoratus Bd (treatment group, Bd exposed) or not exposed to T. marmoratus Bd (control group, not Bd exposed). The treatment group frogs (Bd exposed) were temporarily placed in a two-chambered exposure container with a single Bd-infected T. marmoratus. The exposure time was 5 hours per day for 5 consecutive days. Animals in the control group were not exposed to T. marmoratus because all T. marmoratus were Bd positive. The number of animals in the control group (not Bd exposed) ranged from 4-7 individuals depending on the species, while the number of animals in the treatment group (Bd exposed) ranged from 4-18 individuals, depending on the species. Exposure containers were 1.2 L plastic containers with two parallel mesh barriers affixed across the middle of the container. The mesh barriers created two equal-sized chambers within the exposure container. The mesh barriers were 1 cm apart from each other, preventing direct contact between the treatment frog and the Bd-infected T. marmoratus, but allowing water to flow freely between the two sides. We added 100 ml of water to each exposure container. We examined each exposure container visually every 10 min during each 5 hour exposure period to ensure that frogs stayed in contact with the water. Bd exposed individuals were exposed to the same T. marmoratus, unless the infected T. marmoratus died during the exposure period. When this occurred, the T. marmoratus was immediately replaced by another individual of the same species. Some T. marmoratus were used in multiple exposure sessions, but none were used in more than three sessions.
To determine Bd-infection status of our experimental animals, we collected skin swabs for each study animal at the time of capture, immediately following anti-fungal treatment, and once a week throughout the length of the experiment. For the treatment group (Bd exposed) and all T. marmoratus, we collected a skin swab immediately before and after the five-day exposure period. Skin swabs were analyzed using a standard qPCR protocol [31]. Body weight and SVL was recorded for each animal at the time of capture, immediately following the anti-Bd treatment, and at weekly intervals throughout the experiment.
Control frogs (not Bd exposed) were housed individually in 1.2 L plastic containers with vented isolation tops (lids). We could not expose them to uninfected T. marmoratus because all T. marmoratus were infected, but control frogs were kept in the same type of containers and under similar conditions as the Bd exposed individuals during the exposure period. The experiment began 14 June 2012 and ended 26 August 2012. All experiments were run under ambient conditions (they were outside) so that frogs experienced natural variations in temperature and photoperiod during the experiment. Two wet paper towels were placed in each container. Frogs were fed daily with locally captured invertebrates. Frogs that died during the anti-fungal treatment or during the Bd exposure period were excluded from the analyses. We formalin-preserved deceased individuals following standard museum protocols [32] and deposited the specimens in the herpetological collection of Centro de Ornitología y Biodiversidad (CORBIDI) in Lima, Peru. We previously reported survivorship data for G. excubitor and G. nebulanastes [27], but because this study did not include population abundance data, for sake of simplicity we include these data here. Data for all species are reported in S1 Table. Relative abundance, infection prevalence and intensity in wild populations We used 10x10 m 2 leaf-litter plots and visual encounter surveys to estimate frog species abundances on the basis of species-specific habitat preferences and life history traits. These data were collected during the wet seasons of 1998-1999 (before epizootic) and during the wet seasons of 2008-2009 (after epizootic, [5]) throughout the elevational range of each species, minimizing the importance of other seasonal and local factors affecting intra-and inter-annual variation in population abundance. Albeit climate may also have caused changes in abundance between the sampling years, our two recent sampling seasons differed sharply in precipitation [5], such that our after-epizootic data combine relative abundance from a range of climatic conditions. Leaf-litter plots (four plots for every 100 m change in elevation) were used to estimate relative densities (frogs/ha) for G. excubitor, Psychrophrynella usurpator, Pristimantis danae and P. pharangobates, while visual encounter surveys (covering the entire elevational range) estimated number of individuals per person-hour for G. nebulanastes, H. gladiator and the four species of Pristimantis. We calculated an index of population trends as the ratio between the average of the 2008-2009 estimates and the 1998-1999 estimates (1 = no change). For each species sampled in previous population surveys, we collected skin swabs to detect the presence and infection intensity of Bd from May-October 2012 (dry season) while conducting susceptibility trials. We calculated Bd infection prevalence by dividing the number of infected animals by the total number of animals assayed per species (for infection intensity, see below). Bd infection data are available online at the Amphibian Disease database (https://n2t.net/ark:/ 21547/AXY2).

Molecular methods
To measure Bd infection prevalence and intensity, the skin of each frog was swabbed with a sterile rayon-tipped swab (Medical Wire & Equipment MW113). Swabs were gently stroked across the skin a total of 30 times per frog: 5 strokes on each side of the abdominal midline, 5 strokes on the inner thighs of each hind leg, and 5 strokes on the foot webbing of each hind leg. This technique is widely used and does not harm the animals [33].
We used a quantitative Polymerase Chain Reaction (qPCR) assay to detect Bd and to quantify the intensity of Bd-infection [31]. We followed standard DNA extraction and real-time PCR methods [31], except that we analyzed single-swab extracts once instead of in triplicate [33]. Swabs were extracted using PrepMan Ultra and analyzed using the 7300 Real-Time PCR System (Life Technologies, Carlsbad, CA) at San Francisco State University. The qPCR assay estimates Bd genomic equivalents (GE) in each sample and we convert this to provide "zoospore equivalents" on each frog [12,14] to allow for comparison between studies.

Statistical analyses
We calculated Bd prevalence (proportion of swabbed frogs infected with Bd) by considering swabs as Bd-positive if zoospore equivalents > 0 and Bd-negative if zoospore equivalents = 0. Only species or groups with at least five swabbed individuals were considered for analyses of prevalence and only groups with at least five positive individuals were considered for analyses of infection intensity. We report prevalence data using Bayes credible intervals [34].
We used the Survival Package in R 2.15.2 to compare survivorship of control (not Bd exposed) and treatment frogs (Bd exposed). We used Cox's proportional hazards model with censoring [35,36] to assess the risk of dying for treatment frogs (Bd exposed). For all eight species we performed separate analyses of covariance with initial weight as the continuous explanatory variable and treatment group (control, treatment) as the categorical explanatory variable. We hypothesized that initial weight could have affected survival during the experiments. However, in none of the species initial weight affected survival, and thus weight was subsequently removed from the models.
We compared the relative abundances of the eight species of frogs by using generalized linear mixed effects model with Poisson errors. For leaf litter plots, we used the number of frogs found within each 100 m 2 plot. Relative abundances for visual encounter survey data were standardized by search effort, because duration and number of observers varied among surveys, and these data were rounded to integer prior to analysis. We conducted separate analyses for leaf-litter plots and visual encounter survey data. We only made within-species comparisons because abundance data were restricted to the unique elevational range of each species. Infection intensities were log-transformed (base 10) prior to analyses. We used non-parametric tests (Wilcoxon rank sum, Kruskal-Wallis) to compare maximum infection intensity between control (no Bd exposure) and treatment (Bd exposed) individuals, because there was substantial overdispersion in the data, and sample sizes were small. Averages are followed by the standard error.

Bd prevalence and intensity of infection in the field
Field sampling in 2012 found Bd infection in all species within our study area (Fig 1, Table 1), but prevalence varied among species (chi-square = 55.7, df = 7, p < 0.001) and exceeded 60% in the two riparian species, Hypsiboas gladiator and Pristimantis platidactylus. Infection intensities ranged from 0.01 to 26,208 (Figs 2-4), when compared to unexposed control frogs, but not in the other five species tested ( Infected frogs exhibited symptoms of chytridiomycosis (abnormal skin shedding, tetanic spasms, loss of righting reflex; see videos in S1 File) prior to death.
Phylogenetic relatedness and life-history characteristics did not explain susceptibility to elevated exposure to Bd (Table 2). For example, two closely related direct-developing marsupial species in the genus Gastrotheca showed different susceptibility to Bd; G. nebulanastes was highly susceptible while G. excubitor was not (Fig 2). Similarly, ground-dwelling direct-developing species in the genus Pristimatis also differed in susceptibility: P. platydactylus and P. toftae were susceptible (Fig 5), whereas P. danae and P. pharangobates were not (Fig 4).
There was no difference in survival between treatment and control animals in P. danae, but treatment individuals that died had higher infection intensity (Fig 4), consistent with Bd  (Fig 2), H. gladiator and P. usurpator (Fig 3,  Table 1), and in these species infection intensity in treatment animals (Bd exposed) remained below 1,000 ZE throughout the experiment (S1 Table).

Relative abundances before vs. after the epizootic
The ratio of relative abundances before the epizootic (1998-1999) compared to after the epizootic (2008-2009) was <1 for all species (except for P. usurpator) suggesting population declines (Fig 6; Table 1; S2 Table). These results indicate that the Bd epizootic caused the frog population declines, and that frog populations did not recover during the first decade after the epizootic.

Discussion
Amphibians are experiencing a mass extinction event [37] caused in part by the emergence of a lethal fungal pathogen (Batrachochytrium dendrobatidis, Bd). The disease-induced losses have been substantial (e.g. [8,38]), yet some species survive the epizootic and are assumed to be non-susceptible. However, because the pathogen remains in these systems after an epizootic, the level of susceptibility of surviving hosts has not been explicitly measured. Understanding host/pathogen disease dynamics is these systems is essential for the management and monitoring of surviving amphibian species following mass declines. Several studies have shown that infection intensities are directly tied to survival outcomes in wild populations [12,14]. Our findings support the hypothesis that hosts with high Bd infection intensities (>10,000 zoospore) have reduced survivorship [25]. Sustained, low Bd infection intensity has been found in host populations where host/Bd dynamics are stable [14], while high intensity infections often associate with continued mortality [12,39]. Comparisons Post-epizootic susceptibility to disease of temporal trends in Bd infection intensity of surviving hosts in the wild may help predict future host population stability. Rising Bd infection intensities may signal new outbreaks whereas stable or non-rising infection intensities may signal stable host populations. However, other factors may complicate interpretation of trends in infection intensities; for example occurrence of loss of infection has been shown to vary seasonally [39]. Furthermore, comparisons of Bd infection intensity data must be done cautiously because Bd standards used in quantitative PCR assays may vary in gene copy number [40]; yet, these data can be calibrated between studies based on the qPCR standards used and the gene copy number of the Bd strain infecting amphibians. Disease metrics (intensity and prevalence through time) coupled with population surveys can thus help predict how Bd will affect species in nature. Bd susceptibility trials improve our predictive capabilities regarding Bd effects on host species.
We assessed the effect of increased Bd exposure on survivorship of species that have endured an amphibian mass extinction caused by a Bd epizootic. The epizootic in the Post-epizootic susceptibility to disease Kosñipata Valley of southern Peru caused the loss of 19 of 55 species that occurred from 1200 to 3800m [5]. More than a decade has elapsed since the collapse of these amphibian communities. We examined focal species that continue to exist in the elevational range where Bd caused the highest loss of species. Since existing frogs (i.e., a subset of the 36 species that did not disappear) had survived the community collapse caused by Bd, we expected weak effects, if any, of increased Bd exposure on the focal species. Instead, we found that three of eight species we tested showed significant ongoing susceptibility to increased Bd exposure and our population trend estimates found that seven of our eight species were less abundant after the epizootic compared to our surveys conducted before the epizootic.
Our assessment of species susceptibility to Bd adds to our understanding of how the amphibian community is responding to Bd a decade after the epizootic, but susceptibility alone cannot describe the host/Bd dynamics. Other factors, especially those related to transmission risk and environmental exposure to Bd [5,41,42], likely influence persistence and recovery of amphibian populations at our study site. Field observations suggest that since 2014, populations of some species, such as P. toftae, are now recovering (A. Catenazzi, unpublished data). Although breeding cycles are poorly known for most species, if we assume continuous reproduction throughout the year and short time lag to sexual maturity (within two years) for all species, populations of surviving species could be recovering within 5-10 years unless host/Bd dynamics alternate between enzootic and epizootic stages. Future surveys will quantify the extent of these population recoveries, and more generally the status of populations known to have declined during the first decade after the epizootic. Post-epizootic susceptibility to disease An important finding of our study is that phylogenetic relationships were poor predictors of species susceptibility to Bd. Our results are consistent with findings from community-level assessments showing that the degree of population loss and decline among species is random with respect to the community phylogeny [3,5,8]. Differential susceptibility may reflect species level differences or perhaps differences in prior exposure intensity to Bd, although our experimental findings do not support the latter, immunoprotective hypothesis (see below). We have shown that skin defenses, particularly differences in the proportion of Bd-inhibiting skin bacteria, are associated with difference in susceptibility in the two species of Gastrotheca (Fig  2, [27]), but other factors such as thermal sensitivity of immune response may also affect disease outcome [43,44].
We found that exposure of field-caught frogs to Bd infected T. marmoratus caused mortality in patterns consistent with death by chytridiomycosis, as supported by symptoms such as frequent skin shedding and tetanic spams observed in highly infected individuals shortly before death (S1 File). Most T. marmoratus used for infection had high Bd-infection intensities and therefore presumably shed high numbers of zoospores (e.g., ZE ranging from 1,000-10,000). These values are similar to a proposed threshold that causes death in a number of amphibian species [7,12,25,45] (S1 Table). Our study suggests that frog-to-frog heterospecific transmission (with or without direct contact) can lead to chytridiomycosis and death even in species that were assumed to be resistant to Bd.
The itraconazole anti-fungal treatments did not completely remove Bd infections from study animals, as seen in other studies [46]. Thus, frogs in our Bd exposed treatment were not completely free of Bd when they were exposed to highly infected T. marmoratus. Our anti-fungal treatments did reduce Bd infection intensities on all individuals equally between treatment and control groups. We propose that exposure of animals that already have a low level infection may be more representative of actual dynamics in the wild where Bd is known to remain in systems after epizootics [15,47].
Our experimental design, and failure to clear infection in all individuals, indirectly provides support to previous suggestions that an immunoprotective effect of previous Bd exposure might not be a common response across all amphibian species [48][49][50][51]. Although we ignore the Bd exposure history of our experimental frogs, it is plausible to assume that many frogs at our study site, where infection prevalence can be as high as 90% [24], have previously been exposed to Bd. Among the three susceptible species in our study, P. platydactylus and P. toftae live at elevations with high Bd prevalence [24], and Bd prevalence in their populations is high (Table 1; [5]). Therefore, it is likely that the individuals of P. platydactylus and P. toftae that we used in our study had been infected with Bd prior to our experiment. Mitigation strategies emphasizing immunization should thus take into consideration inter-host differences in the effectiveness of protective response following immunization or prior exposure to live Bd.
Despite recent advances in our understanding of disease dynamics in chytridiomycosis [12,14,52,53], this disease remains a major threat to amphibian biodiversity [37]. Our study was conducted in a montane amphibian community 10 years after a Bd-epizootic eliminated 19 of 55 species. Some species fared worse than others. For example, over half of the stream-breeding species disappeared following the arrival of Bd [5]. We might have predicted that aquatic species would be at highest risk, but we would not necessarily have predicted that a terrestrial marsupial frog (Gastrotheca nebulanastes) would be susceptible (Table 1). Our a priori expectation was that phylogenetic relatedness would predict Bd susceptibility because two genera of frogs (Atelopus and Telmatobius) disappeared completely from Kosñipata Valley during the Bd epizootic, but phylogenetic relatedness does not appear to be an important predictor for susceptibility. Our Bd susceptibility trials demonstrate that chytridiomycosis continues to be a threat even for some species that survived an initial Bd-caused mass extinction and are now in environments that contain carrier hosts infected with Bd, similar to results from other systems [19,20,54]. More importantly, our results suggest that after an epizootic, Bd can still impact amphibian populations and may be regulating their populations at levels lower than pre-epizootic conditions. Our field-based empirical approach provides predictive information that can help determine priority species for conservation action even in regions where epizootics have previously decimated amphibian diversity and abundance.
Supporting information S1 Table. Raw data for susceptibility trials. This is the S1 Fig legend. (DOCX) S1 File. Videos showing symptoms of chytridiomycosis in highly infected frogs used in the susceptibility trials. Videos of diseased individuals used during the experiments. (DOCX) S2 Table. P-values from generalized linear models for change in relative abundances before (1998-1999) and after the epizootic (2008-2009) for leaf litter plots and visual surveys. "-"indicates that no data are available for a given survey and species; Table 2 column refers to data (population ratios and results of GLM analyses) reported in Table 2. When data from both sampling techniques are available, results from leaf litter plots are reported in Table 2 of the manuscript. (DOCX)