1) Although research has focused on density-dependent responses of hosts to single-parasite infections, hosts are exposed to numerous parasites simultaneously under natural conditions and if these exposures lead to infections, they can threaten host populations and ecosystem stability. Moreover, spatiotemporal variation in abundance of co-occurring parasites might influence host infection intensity. If interactions are consistent between different co-infecting parasites, then these patterns could give managers another tool to control disease spread and even predict problematic disease emergences.

2) We investigated how parasite density and identity alter within-host co-infection dynamics. To test this, we simultaneously exposed Cuban treefrogs (Osteopilus septentrionalis) as a model amphibian species to all pairwise combinations of three problematic parasites that commonly co-infect amphibians: the fungus Batrachochytrium dendrobatidis (Bd), the nematode Aplectana hamatospicula, and Ranavirus. Hosts were exposed to one parasite at a fixed dose and another parasite at a range of five doses.

3) Higher doses of Bd decreased Ranaviral and A. hamatospicula loads, but Bd load was not influenced by the dose of either parasite. Ranaviral load was negatively associated with A. hamatospicula dose, but A. hamatospicula load was not affected by Ranaviral dose. We found that all the pairwise co-infections were dependent on parasite density and that pairwise interactions were highly asymmetric – strong in one direction and weak in the other – consistent with interactions dominating food webs.

4) Synthesis and applications: We also revealed that the exposure dose of A. hamatospicula was positively associated with host tolerance to Bd infection and negatively associated with Ranaviral load in hosts. Ranavirus and Bd cause mass die-offs in amphibians, but A. hamatospicula does not. Therefore, in systems where these parasites coexist, maintaining or increasing densities of A. hamatospicula could reduce the negative effects of Bd and Ranavirus infections. Additionally, if these asymmetric and density dependent patterns from community ecology are applicable to other amphibian co-infections or co-infections in other systems, this should allow conservation organizations and resource managers to predict outbreaks and manage host declines associated with deadly parasites by modifying the abundance of co-infecting parasites that might be easier to manage.

Cuban treefrog tadpoles (Osteopilus septentrionalis) were collected from kiddy pools filled with water (140 L) in the University of South Florida Botanical Gardens (Tampa, FL, USA) in August 2016. All animal husbandry throughout the experiment was carried out according to IUCAC protocol #W IS00002203. A. hamatospicula were collected from the gastrointestinal tract of euthanized Cuban treefrogs from Flatwoods Wilderness Park (Tampa, FL, USA). Identical Petri dishes, but with gastrointestinal content from uninfected frogs were used as a sham treatment. Ranavirus (FV3) was cultured in fathead minnow (Pimephales promelas) cells and maintained at -80°C in minimal essential medium (MEM). MEM without Ranavirus was used as a sham. Bd (SRS-JEL 212 strain) was grown in a 1% tryptone solution. SRS-JEL212 was chosen as it was isolated from the Southwestern US, where the experimental frogs were collected. Additionally, this strain successfully infects Cuban treefrogs. Identical plates, but with a sterile 1% tryptone solution were used as a sham. See supplemental methods for more details.

Experimental Design

To examine how the dose of parasite exposures affects host-parasite dynamics, we exposed a total of 174 frogs to one of 28 total parasite treatments. Treatments included exposure to Ranavirus, Bd, or A. hamatospicula alone, simultaneous infections with all pairwise combinations (6 total pairwise combinations) of these parasites at a range of densities (all n=6), and controls (no exposure; n=12). In each of the pairwise co-infection treatments hosts were exposed to one parasite at an intermediate density and the second parasite at one of a range of densities (one of four potential densities) in all possible combinations (Table 1).

Parasite exposures

Two days before parasite exposure, individuals were moved to a 17°C environmental chamber because Cuban treefrogs can clear Bd at higher temperatures (McMahon et al. 2014, Cohen et al. 2017). In the co-infection treatments, the host was exposed to one parasite at an intermediate, fixed dose and simultaneously to the other parasite at a dose that ranged from low to high for each of the tested parasites (Garner et al. 2009, Echaubard et al. 2010, Hoverman et al. 2010, Gervasi et al. 2013, Knutie et al. 2017a, Knutie et al. 2017c). For the microparasites Bd and Ranavirus, doses ranged from 103 to 105 zoospores or plaque forming units (PFU), respectively, and the intermediate dose was 104. For the macroparasite A. hamatospicula, doses ranged from 15 to 60 J3 larvae and the intermediate dose was 30 (Table 1). Doses were applied to frogs held in Petri dishes (25 x 100 mm) for 24h. For Bd and A. hamatospicula exposures, 1-mL of DI water containing Bd zoospores or A. hamatospicula larvae was pipetted directly onto the backs of the hosts. This simulates exposure through water or soil contact (Kilpatrick et al. 2010, Roznik et al. 2021). For Ranavirus exposures, a 69 µl aliquot of MEM with Ranavirus was applied directly into the host’s mouth to mimic fecal-oral or cannibalistic transmission that is common in natural settings (Hoverman et al. 2010). All hosts also received sham exposures for parasites to which they were not exposed. For example, hosts co-exposed to Bd and A. hamatospicula also received the Ranavirus sham treatment. Control individuals received sham treatments for all three parasites.

Assessing parasite load and host health

To assess loads of Ranavirus, hosts were swabbed five times around the mouth and cloaca on day 4 after exposure. To assess Bd load, hosts were swabbed five times from hip to toe on both rear legs on day 16 after exposure. These time points were chosen to measure load differences as they represent times where each parasite would have had time to establish, but not create high levels of mortality (Gray et al. 2009, Voyles et al. 2009, Knutie et al. 2017c). Swabs were stored at -80°C for later processing. DNA from each swab was extracted using a Qiagen DNEasy Blood & Tissue Kit and analyzed using quantitative polymerase chain reaction (qPCR; Boyle et al. 2004, Picco et al. 2007).  A. hamatospicula loads were assessed by counting adult worms in amphibian gastrointestinal tracts when they experienced mortality or at the end of the experiment.

To measure growth, frogs were weighed weekly for four weeks. Individuals were also checked twice daily for mortality. If mortality occurred, frogs were weighed and swabbed and/or dissected, depending on their treatment group. All surviving frogs were euthanized and dissected 28 days after initial parasite exposure.          

Statistical Analyses

All analyses were run with R version 3.6.1 (R Core Team 2019). Plots were created using the visreg package and visreg function (Breheny and Burchett 2019). The survival plot was created using the survminer package and ggsurvplot function (Kassambara et al. 2019).

To test how the interaction between the identity of a co-infecting parasite and its dose altered the load of Ranavirus or Bd, we conducted a generalized linear model with a negative binomial error distribution. For the A. hamatospicula load model, we used a binomial error distribution, and the dependent variable was defined as the proportion of larvae that successfully reached maturity (i.e., using the cbind function on “successes” and “failures”). The independent variables were dose, identity of the changing-dose parasite, and their interaction. Dose was expressed as a log 10-transformed proportion of the highest possible dose hosts were exposed to for each parasite. This allowed us to make comparisons across doses even though the exposure dose varied widely for macro- and microparasites.

To test how the interaction between the identity of a co-infecting parasite and its dose affect host weight, a generalized linear model was run with growth rate as the dependent variable. Growth rates were calculated as the final mass minus the initial mass divided by weeks spent alive (g/wk). To address changes in host survival, we conducted a survival analysis using the survival package and the coxph function (Therneau and Lumley 2019), with host survival as the dependent variable. The identity of the fixed dose parasite, the identity of the changing dose parasite, and the dose of the changing dose parasite were used as interacting independent variables in both host health analyses. To test how these same factors affected tolerance (measured as host growth rate or survival given a parasite burden), the above-described models were rerun, but the interacting predictor variables were the identity of the changing dose parasite, exposure dose, and parasite load of the fixed-dose parasite. For all analyses Tukey Post-hoc tests were run to compare among parasite identities when the variables were significant (multcomp package and glht function (Hothorn 2010). Table S1 outlines all above-described analyses and error distributions.