Amphibian survival compromised by long-term effects of chytrid fungus

Chytridiomycosis, the disease caused by the fungal pathogen Batrachochytrium dendrobatidis (Bd), has been unambiguously implicated in the decline of amphibian populations worldwide. However, the impact of this devastating infectious disease can be difficult to gauge without empirical data on the population-level effects of Bd. Often, assessments of the amphibian chytridiomycosis panzootic are based primarily on expert opinions; as a result, declines in tropical areas are promptly attributed to Bd while its impact on temperate species not suffering from adult mass mortalities is frequently overlooked. Here, we investigated the survival probability in an amphibian species from a temperate area that until now has not been considered to be severely impacted by the disease. Specifically, we related individual survival to Bd infection status using long-term capture-mark-recapture data of male spiny common toads (Bufo spinosus) in Sierra de Guadarrama National Park in central Spain. Even though the study population has demonstrated potential for adaptation to Bd and die-offs of adult individuals have not been recorded, our results clearly indicated that the probability of survival was lower for Bd-positive individuals. Moreover, the probability of becoming Bd-positive was higher than the probability of clearance, driving the population to a slow but certain decline. These results are consistent with other indicators of a negative population trend and suggest that the impact of Bd on temperate species of less concern may be greater than previously thought.


Introduction
Emerging diseases are one of the most dangerous threats to biodiversity, with a fourfold increase in the last 60 years (Daszak et al. 2000;Jones et al. 2008). In particular, fungal diseases have become a major concern for many taxa (Fisher et al. 2012). A good example of this is chytridiomycosis, which is caused by the chytrid fungus Batrachochytrium dendrobatidis (Bd). This disease has caused mass mortality and population declines in amphibians around the world, jeopardizing hundreds of species (Scheele et al. 2019). In the last few decades, numerous studies have reported amphibian declines linked with this emerging disease (e.g., Bosch et al. 2001;Schloegel et al. 2006); however, many have focused only on short-term impacts (e.g., mass mortalities) on susceptible species in wet climates of Central America, South America, and Australia (Scheele et al. 2019). The long-term impacts of chytridiomycosis (e.g., slow population declines due to reduced survival), and its effects on species and areas traditionally considered to be less vulnerable, remain poorly known (but see Briggs et al. 2010;Seimon et al. 2017;Valenzuela-Sánchez et al. 2017;Bosch et al. 2018).
In temperate areas of Europe, chytridiomycosis has impacted many species (Garner et al. 2005). At the Peñalara Massif in the Sierra de Guadarrama National Park (Central Spain), mass mortalities of metamorphs of midwife toads (Alytes obstetricans), fire salamanders (Salamandra salamandra), and spiny common toads (Bufo spinosus) were first reported more than twenty years ago (i.e., 1997-1999) and represented the first report of chytridiomycosis in Europe (Bosch et al. 2001(Bosch et al. , 2014Bosch and Martínez-Solano 2006;Garner et al. 2009). Initially, common toads benefited from the almost complete disappearance of midwife toads and colonized several ponds in this area (Bosch and Rincón 2008). However, a few years later, a demographic study of the common toad population of Laguna Grande, the biggest pond in Peñalara, revealed an incipient decline and a low recruitment (Bosch et al. 2014). Indeed, from 1999 to 2016 there was a detectable decrease in the number of clutches of B. spinosus in the few permanent ponds located in different drainages of the massif . At the same time, a quantitative genetic study carried out in 2013 revealed that this toad population might possess the ability to adapt to chytridiomycosis, since the heritable basis of the Bd load has a predominantly additive component (Palomar et al. 2016). This means that parents exhibiting low Bd load will have offspring exhibiting low Bd load, if infected. If the Bd load is indeed heritable in these toads, the population might be able to evolve some resistance or tolerance to the disease, adapt to it, and counteract its decline. However, recent population trends have not been characterized and it remains unknown if the decline persists or if this species has managed to adapt to chytridiomycosis.
Capture-mark-recapture (CMR) studies enable the investigation of population dynamics of wildlife species on a long-term scale, and are considered an effective tool in clarifying the impact of diseases on natural populations (e.g., Murray et al. 2009;Phillott et al. 2013). In particular, the multi-event modeling approach, an extension of classical CMR models, allows researchers to test whether and how survival and recapture probability are affected by individual biological states, such as breeding site (Fernández-Chacón et al. 2013), size class (Fernández-Chacón et al. 2015) or infection status (Conn and Cooch 2009;Muths et al. 2020). Previous research that has applied CMR models to the epidemiological study of amphibian populations affected by chytridiomycosis has reported a decrease in the survival probability of infected individuals and populations (e.g., Murray et al. 2009;Pilliod et al. 2010;Phillott et al. 2013;Brannelly et al. 2015;Russell et al. 2019). Furthermore, 1 3 these studies have revealed that this lethal disease may continue to impair amphibian survival even after decades of coexistence (e.g., 30 years in Murray et al. 2009;20 years in Phillott et al. 2013).
In the present study, we took advantage of 13 years of CMR monitoring to evaluate, for the first time, how chytridiomycosis affects adult survival in an amphibian population in decline but with proven evolutionary potential for adaptation to the disease. We hypothesized that Bd plays a key role in shaping the demographic trend in the population of spiny common toads in Laguna Grande, even though no dead adults have been reported and common species from temperate zones are traditionally considered less susceptible to chytridiomycosis.

Study area
The Peñalara Massif is a protected alpine area (768 ha) in the Sierra de Guadarrama National Park (N 40.839877, W -3.957451) composed of several small glacial valleys located 1800-2430 m a.s.l. The landscape consists mainly of granitic outcrops, alpine meadows, heathlands dominated by Cytissus oromediterraneus and Juniperus communis nana, and forests of Pinus sylvestris below the timberline. Within the massif there are approximately 250 catalogued ponds where up to nine amphibian species breed: B. spinosus, Rana iberica, Pelophylax perezi, Epidalea calamita, A. obstetricans, Hyla molleri, S. salamandra, Ichthyosaura alpestris, and Triturus marmoratus. Laguna Grande (elevation 2018 m a.s.l., surface area 5452 m 2 , maximum depth 4.7 m) is one of the largest ponds in the massif and harbors the majority of the spiny common toad population.

Sampling
Individual spiny common toads were captured, marked, and recaptured from 2006 to 2018 at Laguna Grande. Multiple nighttime capture sessions (ca. 15) were conducted during the breeding season of each year (from April to June). During each session, all adult toads on the shore of the pond were captured by dip net and their sex and PIT (passive integrated transponder; Microplus 2 × 12 mm, Insvet Inc., Huesca, Spain) tag number were recorded. If an animal lacked a tag, a new one was inserted subcutaneously using a syringe. All animals were released immediately after the procedure. Furthermore, from 2008 to 2016, a randomly chosen sub-sample of animals (588 of 1369 captures) was tested for Bd. A sterile swab was rubbed 20 times on each animal following standard protocols (Hyatt et al. 2007). Swabs were stored dry at 4 °C and analyzed less than two weeks after collection for the presence of Bd with the Prepman extraction and TaqMan real-time PCR protocol (following Boyle et al. 2004). Individual swabs were analyzed in duplicate and with Bd genomic equivalent (GE) standards of 100, 10, 1, and 0.1, as well as a negative control.

Data analysis and modeling approach
The field records of marked toads were used to build a 13-year capture-recapture data set. Although it would have been desirable to include females, the dataset contained insufficient 1 3 data on recaptures for females to allow the construction of a model (246 marked in 13 years and just 69 recaptured); for this reason, only males were used for the following analyses. The male data set contained, for each sampling occasion, information on whether the individual was captured or not and, if captured, its infection status. This information was formatted and analyzed using a multi-event modeling approach (Pradel 2005), an extension of classical capture-recapture models that links field observations to a series of underlying individual states defined in the model structure (see below and Online Appendix A). By defining parameters that allow for changes in biological states between captures, multievent models provide an ideal framework for epidemiological studies by estimating the probability of an animal acquiring an infection and recovering from it (Conn and Cooch 2009).
Our field data consisted of four types of observations or "events" that were codified in the dataset as follows: "not seen" (0), "seen alive and Bd-negative" (1), "seen alive and Bd-positive" (2), and "seen alive with infection status unknown" (3). From this set of events, we estimated annual toad survival and infection rates by constructing a model pattern based on transition matrices that linked the observed events to transitions between the underlying states in which individuals could be found at a given sampling occasion ( Fig. 1). We modeled transitions among three states: alive Bd-positive (infected, "I"), alive Bd-negative (uninfected, "U"), and dead ( †); between sampling occasions, toads could change state according to the transitions shown in Fig. 1. The probabilities associated with each change of state are defined in the full transition matrix (Φ), which can be written as: where S: the annual apparent survival probability or the probability of surviving while remaining available for sampling within the study area. As sampling occurred only within the study area, it was not possible to distinguish between permanent emigration and mortality, so the survival estimates reported for our system correspond to apparent survival

Biological processes Observational processes
Bd positive toad  Bd negative (U) as well as the shifts between states, for instance from Bd positive to negative I»U 1 3 (hereafter "survival"). Ψ: the probability of transitioning from one state to another, conditional on survival. Two types of transitions are possible: Ψ U→I : infection transitions (from Bd-negative to Bd-positive state). Ψ I→U : clearance transitions (from Bd-positive to Bd-negative state). These model parameters could be estimated separately by splitting the full transition matrix into a two-step series of transition matrices representing survival and state transition processes. Our model pattern assumes that ecological processes occur before the observational ones, with survival being the first step in our sequence of transition matrices. The last part of the sequence corresponds to the observational process, which is also split into two matrices, mirroring how observations are made in the field. This was necessary in our case due to the existence of ambiguous observations (i.e., infection status unknown) in our data (for a similar approach see Conn and Cooch 2009). Matrices B1 and B2 (below) represent the observational process.
where p: the probability that a marked toad, alive and in the study area, is recaptured. δ: the probability of state assignment of recaptured toads, with two possibilities: δ I : the probability that a toad in state "I" is recorded as Bd-positive. δ U : the probability that a toad in state "U" is recorded as Bd-negative.
In matrix B1, individuals from each alive state (I and U) can either be encountered with probability p or not encountered with probability 1−p. Dead individuals ( †) are not recorded in the field, so they remain always unobservable with probability = 1. Thus, we have three outcomes: "not encountered" ( −), "encountered" for Bd-positive individual (I), and "encountered" for Bd-negative individual (U). These outputs from matrix B1 (columns) become the input (rows) of matrix B2. For an "encountered" Bd-positive individual (row 2), the status may then be assessed (δ I ), in which case the final event will be "seen alive and Bd-positive" (event 1, second column of B2 matrix) or not assessed (1−δ I ), in which case the event will be "seen alive with infection status unknown" (event 3, last column of B2 matrix).
Our multi-event models were built and fitted to the data using the program E-SURGE (Choquet and Nogue 2010; see Online Appendix A), but prior to the model selection process, a Goodness-of-fit (GOF) test was conducted to check if our data met the assumptions of a departure model that considers all parameters to be state-and time-dependent, namely the Arnason-Schwarz (AS) model (Pradel et al. 2003). GOF tests were performed using U-CARE (Choquet et al. 2009), a statistical program that helps users to detect sources of lack of fit in their capture-recapture data (mainly caused by differences in survival and recapture probabilities among individuals) and to redefine the structure of the departure model to accommodate these heterogeneities. The stratification of data in different individual states can already control for some sources of heterogeneity, but in order to scale model deviances and correct for the remaining sources of lack of fit, we introduced the overdispersion coefficient, or ĉ (calculated as the sum of chi-square results for each test divided by the total number of degrees of freedom), when performing the analyses in E-SURGE.
Model selection was based on Akaike's information criterion corrected for overdispersion (QAIC) and we retained as good candidate models those with the lowest QAIC values (Burnham and Anderson 2002). The model selection process began with a general model considering full time (year) and state effects on survival (S), recapture (p), and transition probabilities (Ψ I→U and Ψ U→I ) and constancy in state assignment probabilities (δ U and δ I ).
The modeling process consisted of removing state and time interactions ("*") from S, Ψ, and p parameters and of testing constancy (".") and alternative additive ( +) combinations of time and state until the most parsimonious model structure was found. We began testing such effects on recapture probability (p) first, and retained the most parsimonious structure for p in the subsequent modeling step. Then, we followed the same modeling procedure for the S and Ψ parameters until a final model was obtained.
In order to confirm the decline and the low recruitment found previously in the studied population with 5-year data (Bosch et al. 2014), we also performed a time-symmetric open model or Pradel model (Pradel 1996) to the male toad CR dataset using the program MARK (White and Burnham 1999). This model allowed the estimation of λ (population growth rate), φ (apparent survival), and f (recruitment). We conducted two analyses. In the first one, we estimated survival and recruitment considering survival as time depending (as found in the above multi-event model) and testing two models: constant recruitment and time dependent recruitment. In the second analysis, we estimated survival and population growth considering both survival and growth rate as time depending. More details about the models are provided in supplementary material (Online Appendix A).

Results
Over the 13 years of the study (2006-2018), we captured a total of 879 toads, 632 of which were males (72%). Per breeding season, an average of 150 individuals (from 86 to 292) were found, some of which were recaptured several times. Males were recaptured a maximum of nine times (nine years) and females a maximum of 13. However, the number of captures per year decreased with time (Fig. 2). Each year a subsample of captured toads (ranging from 46 individuals, 21% of captured, in 2010 to 152 individuals, 86%, in 2009; average of 65 individuals, 43%) was tested for Bd, which was found to have an overall prevalence of 27.79% (see Daversa et al. 2018 for more detailed information about infection dynamics over the years).

Capture-recapture dataset, goodness-of-fit testing, and model selection in multi-event model
During our 13-year study period, 632 male toads were marked and released in the sampling area and 438 animals were subsequently reencountered. The GOF tests performed for this encounter dataset yielded significant results, indicating some lack of fit (see Online Appendix B); we therefore applied the obtained ĉ value of 1.474 to scale model deviances in the subsequent modeling step.
In the modeling process, we began with a multi-event model structure that considered the effect of full time and state interactions on survival, recapture, and transition 1 3 parameters (Model 4, Table 1; see all tested models in this table). This multi-event departure model included extra observational parameters that accounted for uncertainty in state assignment, thus incorporating potential sources of lack of fit. State assignment probabilities were always kept state-dependent (see Methods). Removing state effects on recapture probability (p) increased model parsimony (Model 3 vs. Model 4), so we used this structure to start modeling the survival (S) and transition (ψ) parameters. Removing time and/ Year

Survival, transition, recapture, and state assignment probabilities in multi-event model
The annual survival of male toads varied over time, with similar trends for Bd-positive and Bd-negative individuals. Survival peaked during the central years of the study period and fell in later years (Fig. 3a). Mean annual survival (estimate ± standard error) was lower for Bd-positive (S I = 0.754 ± 0.016) than for Bd-negative individuals (S U = 0.883 ± 0.027). Transition probabilities also varied over time, but unlike state-dependent survival, infection and clearance transitions showed different temporal patterns (Fig. 3b). On average, the probability of becoming Bd-positive (Ψ U→I = 0.583 ± 0.039) was higher than the probability of becoming Bd-negative (Ψ I→U = 0.225 ± 0.019). Annual recapture probabilities (p) ranged from 0.375 ± 0.056 to 0.794 ± 0.042, with a mean value of 0.593 ± 0.017 (Fig. 3c).
Estimates of state assignment probabilities (δ) indicated that captured Bd-negative toads were all correctly identified as Bd-negative (δ U = 1.000; standard error could not be computed by the model) and none was classified as "unknown". In the case of captured toads belonging to the Bd-positive state, the probability of them being identified as Bd-positive was low (δ I = 0.119 ± 0.011), indicating that most encountered Bd-positive individuals were classified as "unknown".

Discussion
In this study, we investigated the probability of adult survival in a population of spiny common toads infected with Bd. We found that the impact of chytridiomycosis is still detectable despite the number of years that has passed since the initial outbreak (~ 25 years) and the adaptive potential of the population (i.e., the heritability associated with Bd load). Our estimates of survival probability were lower in Bd-positive individuals, and the probability of becoming Bd-positive was higher than the probability of clearance. Assuming that this pattern, identified in males, also applies to females and taking into account the low estimated recruitment and population growth rate, the population seems destined to suffer a slow but steady decline which will likely increase its vulnerability to environmental stochasticity. This work highlights the uncertainty that chytridiomycosis can induce even in populations of temperate-climate species of least concern.

Effect of chytridiomycosis in adults
We found that the estimated probability of survival was lower for Bd-positive individuals in all years of our study except one (Fig. 3a). These results resemble the findings of Year Fig. 4 Recruitment, survival, and population growth estimates of male toads during the study period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). Bars represent the 95% confident intervals other CMR studies in which chytridiomycosis affected adult survival even when decades had passed since the disease outbreak (Murray et al. 2009;Pilliod et al. 2010;Phillott et al. 2013;Brannelly et al. 2015;Valenzuela-Sánchez et al. 2017;Russell et al. 2019). This lower survival probability might be due to Bd-induced increase in adult overwintering mortality, which has been suggested as a possible cause for decline in populations without recorded annual mass mortalities (Wetsch et al. 2022;Valenzuela-Sánchez et al. 2017). Mean estimate of growth rate showed a declining population. Indeed, if we consider the number of captures per year as a proxy of the abundance of toads, we detected a significant negative trend over the studied period (Fig. 2). This was not the result of a decline in sampling effort, since we did not see the same trend in recapture probability (Fig. 3c). To date, chytridiomycosis-related mortality of spiny common toads has only been reported in metamorphs (e.g., Bosch et al. 2001;Bosch and Martínez-Solano 2006;Garner et al. 2009) and this alone might have driven this negative population trend due to lack of recruitment. However, our finding, that adult survival could be also compromised, reveals the complexity of the disease dynamics and the number of factors that must be taken into account in evaluating its impacts.

Evolutionary implications
As the pathogen becomes endemic, the long-term survival of many amphibian species in the wild will depend on evolutionary adaptation. Species can employ two non-exclusive strategies to cope with a disease: resistance, which limits pathogen burden, and/or tolerance, which reduces its fitness costs. These two responses are expected to have different evolutionary consequences. Alleles conferring disease tolerance tend to become fixed by selection since, in a tolerant population, disease incidence remains high and the alleles will remain advantageous. In contrast, alleles conferring resistance tend to be polymorphic, since the incidence of disease will decrease in a resistant population and the allele may no longer be advantageous (Roy and Kirchner 2000). Our toad population in Laguna Grande suffered mass mortalities early on, which might have caused strong selection for resistant and tolerant individuals. While tolerant variants might have spread through the population, resistant variants might be polymorphic, and susceptible individuals could still be suffering from the impact of Bd. In amphibians, less susceptible species or individuals present traits such as resistant MHC alleles, reduced skin inflammatory responses, and increased expression of genes involved in skin integrity. In general, more robust innate and adaptive immune responses are expected to happen at the critical early subclinical stage of infection (reviewed in Zamudio et al. 2020). An investigation of some of these patterns in our population, for instance a comparison between past and current MHC allele pools, could shed some light on the evolution of defenses against Bd infection. The long-term persistence of this population is likely to be influenced by factors other than selection for resistant/tolerant individuals. Since migration from other ponds is unlikely (due to the distances involved and the steepness of the mountains), demographic compensation could be another target of selection. High recruitment might counterbalance the mortality caused by Bd, as shown in other amphibian populations (e.g., Muths et al. 2011;Valenzuela-Sánchez et al. 2022). The fact that infection status did not influence recapture probability during the breeding season (Fig. 3c) might suggest that Bd-positive males have increased their reproductive efforts, aiming for compensatory recruitment. However, high-elevation amphibian populations are very limited in their ability to exhibit compensatory recruitment (Hardy et al. 2022).

3
In fact, low recruitment in this pond has been demonstrated (Bosch et al. 2014;this study) and the number of toad egg strings in the area has decreased over the last two decades , which would seem to indicate that this hypothesis is also unlikely. Finally, density dependence during larval stage could also contribute to compensate the population decline (Wilbur 1977;Briggs et al. 2010), but, to date, there is no evidence of success.

High survival and high longevity
The estimated survival probabilities in our amphibian population, although different between Bd-positive and Bd-negative individuals, were quite high (Fig. 3a). Indeed, the estimates presented here are higher than those reported for other common toad populations (e.g., Gittins 1983;Schmidt and Anholt 1999;Frétey et al. 2004;Tomasevic et al. 2008) and are more similar to those of other toad species inhabiting high mountains (e.g., Scherer et al. 2005;Muths et al. 2010). The documented high longevity of common toads in high altitudes (Hemelaar 1988;Cvetković et al. 2009) might be related to this higher survival, since the rate of mortality can increase with age (Kirkwood and Austad 2000). In highlands, male and female toads can live 11 and 15 years, respectively (e.g., Schabetsberger et al. 2000). In our population, we have recaptured some males for as long as 9 years and some females for 13 years. Assuming that they start to reproduce at 4-6 years of age (Hemelaar 1988;Reading and Clarke 1995), longevity in our population may reach 13-14 years for males and 18-19 years for females; values that are higher than those reported in previous studies (Hemelaar 1988;Cvetković et al. 2009;Schabetsberger et al. 2000). Several predictions stemming from classical evolutionary theories of senescence could support the evolution of such a high longevity (Kirkwood and Austad 2000;Brys et al. 2007;Stark and Meiri 2018). Aging evolves to be delayed in safe environments (those with low mortality) and adaptations that reduce mortality are generally linked to increased longevity. Cold environments, such as high altitudes, allow for slower growth and metabolic rates, reducing the strength of intrinsic drivers of mortality related to deleterious mutations or oxidative damage (Brys et al. 2007;Stark and Meiri 2018;Cayuela et al. 2021). In addition, spiny common toads are nocturnal and utilize chemical protection, which both reduce extrinsic mortality pressures (Stark and Meiri 2018). Finally, the fact that our study site is in high mountains within a national park might also have helped to reduce threats to this population in the last few decades.
The high degree of longevity that this toad population has evolved might play a role in its persistence. A long lifespan can allow them to survive at lower densities for a time and possibly to overcome short-lived threats (Pimm et al. 1988). For instance, yearly climatic conditions in the mountains are quite variable and usually affect toad survival (e.g., Muths et al. 2020); we hypothesize that this is the reason why our survival estimates and transition probabilities were highly variable over the duration of the study (Fig. 3a). Although high longevity might help toads survive environmental stochasticity (at least if the stochastic events are short in duration), this trait alone might be not sufficient to compensate for the prolonged impact of chytridiomycosis (Chichorro et al. 2019). In fact, the negative consequences caused by Bd could make the population more vulnerable to stochastic or climatic events.

Conclusions
In conclusion, our study demonstrates how amphibian survival can be affected by chytridiomycosis long after the initial disease outbreak. This fits with recent experimental work demonstrating significant mortality due to Bd during overwintering in temperate amphibians (Wetsch et al. 2022). The spiny common toad population of Laguna Grande within the Sierra de Guadarrama National Park is in decline due, at least in part, to Bd. The impact of Bd on adult and metamorph survival should impose pressure toward increasing recruitment and developing resistance and/or tolerance. However, we did not find any evidence that these mechanisms are counteracting the population decline. The high survival probability and longevity of these animals, together with their demonstrated potential for adaptation to Bd, might still give them the time and tools necessary to evolve some kind of defense. Although available mitigation strategies in-situ are scarce and complex (Pabijan et al. 2020), our findings indicate the urgent need to initiate conservation measures to prevent the population collapse, increasing the opportunity for this adaptation. For example, Hudson et al. (2016) found that in-situ itraconazole treatments of adult amphibians increased their probability of survival, and modified strategies from Bosch et al. (2015) using chemical mitigation at the pond-scale are now being tested at the study National Park. At present, though, this population of spiny common toads has an uncertain future, even when the species has not been considered until now to be threaten by chytridiomycosis (sensu Scheele et al. 2019). Our study reveals how the impacts of chytridiomycosis can be underestimated: Bd is well-known for causing mass mortalities and sharp population declines in very susceptible species, but can also be responsible for slow and steady long-term declines in other species considered less vulnerable. This work also demonstrates that long-term studies on the effect of chytridiomycosis, while still scarce (but see e.g., Russell et al. 2019), are essential in identifying the species or regions that are threatened by this emerging disease.