Recovered frog populations coexist with endemic Batrachochytrium dendrobatidis despite load‐dependent mortality

Abstract Novel infectious diseases, particularly those caused by fungal pathogens, pose considerable risks to global biodiversity. The amphibian chytrid fungus (Batrachochytrium dendrobatidis, Bd) has demonstrated the scale of the threat, having caused the greatest recorded loss of vertebrate biodiversity attributable to a pathogen. Despite catastrophic declines on several continents, many affected species have experienced population recoveries after epidemics. However, the potential ongoing threat of endemic Bd in these recovered or recovering populations is still poorly understood. We investigated the threat of endemic Bd to frog populations that recovered after initial precipitous declines, focusing on the endangered rainforest frog Mixophyes fleayi. We conducted extensive field surveys over 4 years at three independent sites in eastern Australia. First, we compared Bd infection prevalence and infection intensities within frog communities to reveal species‐specific infection patterns. Then, we analyzed mark‐recapture data of M. fleayi to estimate the impact of Bd infection intensity on apparent mortality rates and Bd infection dynamics. We found that M. fleayi had lower infection intensities than sympatric frogs across the three sites, and cleared infections at higher rates than they gained infections throughout the study period. By incorporating time‐varying individual infection intensities, we show that healthy M. fleayi populations persist despite increased apparent mortality associated with infrequent high Bd loads. Infection dynamics were influenced by environmental conditions, with Bd prevalence, infection intensity, and rates of gaining infection associated with lower temperatures and increased rainfall. However, mortality remained constant year‐round despite these fluctuations in Bd infections, suggesting major mortality events did not occur over the study period. Together, our results demonstrate that while Bd is still a potential threat to recovered populations of M. fleayi, high rates of clearing infections and generally low average infection loads likely minimize mortality caused by Bd. Our results are consistent with pathogen resistance contributing to the coexistence of M. fleayi with endemic Bd. We emphasize the importance of incorporating infection intensity into disease models rather than infection status alone. Similar population and infection dynamics likely exist within other recovered amphibian‐Bd systems around the globe, promising longer‐term persistence in the face of endemic chytridiomycosis.


Laboratory protocol to detect and quantify Bd infections
DNA was extracted from swab tips with Prepman ® Ultra (Applied Biosystems) using a standard protocol (Boyle et al. 2004, Hyatt et al. 2007) without bead-beating step (Brannelly et al. 2020). Swab tips were cut off using scissors into 1.5 ml tubes containing 70 l Prepman ® Ultra. In order to avoid contamination, all tubes remained closed until respective swab tips were cut and scissors were cleaned with freshly made 2% bleach solution and rinsed with Milli-Q water between each sample. A negative control swab was included with every run of extractions to test for contamination, of which we observed no evidence. After cutting swab tips, tubes were briefly vortexed to ensure mixing with Prepman ® Ultra and incubated at 100ºC for 10 min in an oven. Samples were allowed to cool for 2 min and centrifuged for 3 min at 13,000 rpm. Supernatant was recovered and 5 l was diluted 10 -1 with UltraPure Distilled Water (Invitrogen) for diagnosis using quantitative Polymerase Chain Reaction (qPCR).
Swab samples were run in duplicate and considered positive when at least one well amplified more than 0 ITS copies. Our assay consistently recovered small Bd loads which likely represent low infection levels (Hyatt et al. 2007, Briggs et al. 2010. We report infection intensity as log 10 ITS copies per swab because the ITS copy number of local Bd strains is unknown (Longo et al. 2013). To generate ITS copies per swab, estimated numbers of ITS copies per well were averaged over the duplicate runs and multiplied by the dilution factors in the extraction process.

Mark-recapture modeling
We fit a robust design, multisite, multistate Arnason-Schwarz (Arnason 1972, 1973, Schwarz et al. 1993, Kéry and Schaub 2012 using a continuous-time formulation of the ecological process (Miller and Andersen 2008, Conn et al. 2012, Glennie et al. 2022). Our robust design consisted of 19-21 primary occasions per site with 2-3 secondary surveys per primary occasion per site. We modeled the ecological process with hazard rates and matrix exponentials to account for unequal time intervals during the first two years of the study, and because this method is more appropriate for modeling processes occurring instantaneously in continuous-time (Ergon et al. 2018). We specified six-weeks as the baseline interval between primary occasions-note that this implicitly assumes that infection state transitions were not expected to occur more than once over six-weekly time periods. We conducted our analysis with NIMBLE 0.12.2 (de Valpine et al. 2017) in R 4.1.0 (R Core Team 2021).

Ecological process
We modeled three latent ecological states ( ): (1) alive and Bd-, (2) alive and Bd+, and (3) dead. Starting after an individual's first capture, the latent ecological state of individual during primary occasion at site was modeled conditional on the ecological state during primary occasion −1: where the transition probability matrix (TPM) of the ecological process TPM z = TRM z , where is the vector of site-specific primary occasion intervals in units of six-weeks and TRM z is the transition rate matrix (TRM) of the ecological process (Miller andAndersen 2008, Conn et al. 2012): Above, 1 is the mortality hazard rate of uninfected individuals, 2 is the mortality hazard rate of infected individuals, 12 is the hazard rate of gaining Bd infection, and 21 is the hazard rate of clearing Bd infection. We calculated TRM z = V D V −1 , where V is the matrix of eigenvectors of TRM z and D is the diagonal matrix of eigenvalues of TRM z (Miller andAndersen 2008, Conn et al. 2012). Note that when using NIMBLE, one can use nimbleRcall to call the expm::expm() function within the model code (Goulet et al. 2021), but this approach is approximately 10 times slower than using NIMBLE's built-in matrix functions.
All hazard rates were modeled at the level of site-specific primary occasions. We modeled mortality hazard rates as log-linear functions of body condition (scaled mass index, Peig and Green 2009), average daily temperature over the primary interval (a proxy for season), Bd infection status ( 2 only), Bd infection intensity ( 2 only), and all pairwise interactions. After standardizing predictors, including Bd infection intensity, the effect of Bd infection status can be interpreted as the effect of being infected with Bd with an average infection intensity on the log mortality hazard rate. Hazard rates of infection transitions were modeled as log-linear functions of body condition, average daily temperature over the primary interval, average six-weekly rainfall over the primary interval, their interaction, and Bd infection intensity ( 21 only). All mortality and infection transition parameters additionally included random site and survey effects, which were drawn from two multivariate normal distributions, respectively. Note that the mortality hazard rates are 'apparent' mortality hazard rates, as mortality could not be disentangled from permanent emigration from the sites.

Observation process
The observation process was modeled conditional on the latent ecological state, with three possible observed states ( , data): (1) seen/recaptured and Bd-, (2) seen/recaptured and Bd+, and (3) not seen/recaptured. Note that the observed state was assigned based on the collected swab samples. The observed state of individual , primary occasion , secondary survey , at site was drawn from a categorical distribution conditional on the latent ecological state: where 1 and 2 are the recapture probabilities of individuals uninfected and infected with Bd, respectively. We modeled recapture probabilities at the level of secondary surveys as logit-linear functions of body condition, sex, daily temperature on the survey day, cumulative rainfall over the primary interval, their interaction, and Bd infection status and intensity ( 2 only). We additionally included random site effects, random secondary-survey level effects, and random individual effects. We included the individual effects to account for heterogeneity due to differences in individual home ranges, transience, and position on the transect.

Priors
For hazard rates, we used Exponential(1) for the log intercepts of log-linear functions because this corresponds to a uniform distribution on the baseline hazard rates transformed to probabilities. For recapture probabilities, we used Beta(1, 1) for the logit intercept because this is a uniform distribution on the baseline recapture probabilities. We used weakly informative (Half-)Cauchy(0, 2.5) on all predictor variables and standard deviations of random effects (Gelman et al. 2008).

Reversible Jump MCMC
For predictor variable selection, we used NIMBLE's built-in reversible jump MCMC (RJMCMC, Green 1995). RJMCMC samples across models with different dimensionalities (e.g., more or fewer predictors), and readily facilitates model selection directly within the MCMC algorithms. When RJMCMC includes a predictor in the model, the coefficients are estimated; when RJMCMC excludes a predictor, the coefficient is toggled to 0. After running the MCMC, the inclusion probability of predictors can be retrieved by calculating the proportion of MCMC samples for which a predictor was included in the model. Important predictors will be included in nearly all iterations, and predictors with little predictive potential will be excluded. We constrained interaction effects to only be included in the presence of main effects to maintain the principle of marginality.

Imputation of missing values
NIMBLE does not allow missing values in predictor variables. Therefore, for predictors with missing values, in particular matrices of time-varying individual covariates, missing values need to be provided. Bayesian analysis readily facilitates the imputation of these missing values using submodels for the covariates (Gelman et al. 2013). We imputed missing values in three predictor variables: body condition, infection intensity, and sex. For body condition, we imputed missing values from a normal distribution of the observed body condition values with random individual effects to account for repeat measures over the study. For infection intensity, we similarly imputed missing values from a normal distribution with random individual effects, but also included temperature, rainfall, and their interaction to account for seasonal differences in infection intensity. For sex, we imputed 15 missing values from a Bernoulli distribution roughly centered on the observed sex ratio.

Posterior predictive checks
We assessed goodness-of-fit through posterior predictive checks (PPCs, Gelman et al. 1996) inspired by Rankin et al. (2016) and Kéry and Royle (2020). PPCs are conducted by simulating replicate datasets ( rep ) from the joint posterior distribution for each MCMC iteration, calculating some summary statistics for the observed and replicate datasets, and then comparing some discrepancy statistics between observed and replicated datasets. By comparing the discrepancy statistics from the data ( ) with the discrepancy statistics from the replicate datasets ( rep ), we assess to what extent our model predicts capture histories that are consistent with our observed data. Bayesian pvalues summarize the similarities between observed and replicated data, calculated as Pr( rep > ), where values around 0.5 imply good fit. Note that our PPCs exclude the primary occasion of first capture, because this occasion is not modeled and thus not simulated in the replicate datasets.
For the first PPC, we counted the number of observed infection state transitions (going from uninfected to infected and vice versa) for each individual in the study. Note that infection state transitions can only be observed if individuals are observed during consecutive primaries in different states. We then calculated Freeman-Tukey statistics on these observed state transitions to calculate a Bayesian p-value. This PPC was conducted to assess the fit of the ecological process.
For the second PPC, we counted the number of times each individual was observed as either uninfected or infected with Bd during a secondary survey. Our summary statistics were thus the total number of recaptures of each individual for each alive state. We computed Freeman-Tukey statistics and calculated Bayesian p-values. This PPC was conducted to assess the fit of the observation process.